|
|
@@ -0,0 +1,709 @@ |
|
|
{ |
|
|
"cells": [ |
|
|
{ |
|
|
"cell_type": "markdown", |
|
|
"metadata": {}, |
|
|
"source": [ |
|
|
"# Heat loss calculation\n", |
|
|
"\n", |
|
|
"Example of a small residential building with a heated floor area of 100 m<sup>2</sup>.\n", |
|
|
"\n", |
|
|
"The example is taken from the SIMIEN wiki documentation pages." |
|
|
] |
|
|
}, |
|
|
{ |
|
|
"cell_type": "code", |
|
|
"execution_count": 692, |
|
|
"metadata": {}, |
|
|
"outputs": [ |
|
|
{ |
|
|
"data": { |
|
|
"text/markdown": [ |
|
|
"### Resulting HLC" |
|
|
], |
|
|
"text/plain": [ |
|
|
"<IPython.core.display.Markdown object>" |
|
|
] |
|
|
}, |
|
|
"metadata": {}, |
|
|
"output_type": "display_data" |
|
|
}, |
|
|
{ |
|
|
"data": { |
|
|
"text/markdown": [ |
|
|
"Transission losses: $H_{tr} = \\sum{U_i \\cdot A_i} = 66.8$ $W/K$ , where $i$ is each envelope surface part" |
|
|
], |
|
|
"text/plain": [ |
|
|
"<IPython.core.display.Markdown object>" |
|
|
] |
|
|
}, |
|
|
"metadata": {}, |
|
|
"output_type": "display_data" |
|
|
}, |
|
|
{ |
|
|
"data": { |
|
|
"text/markdown": [ |
|
|
"Infiltration losses: $H_{inf} = 0.33 \\cdot n_{inf} \\cdot H_b \\cdot A_{fl} = 14.4$ $W/K$ , where $n_{inf} = 0.07 n_{50}$ for balanced ventilation ref. NS 3031:2014" |
|
|
], |
|
|
"text/plain": [ |
|
|
"<IPython.core.display.Markdown object>" |
|
|
] |
|
|
}, |
|
|
"metadata": {}, |
|
|
"output_type": "display_data" |
|
|
}, |
|
|
{ |
|
|
"data": { |
|
|
"text/markdown": [ |
|
|
"Ventilation losses: $H_{ve} = 0.33 \\cdot \\dot V \\cdot (1-\\eta_{T}) \\cdot A_{fl} = 11.9$ $W/K$" |
|
|
], |
|
|
"text/plain": [ |
|
|
"<IPython.core.display.Markdown object>" |
|
|
] |
|
|
}, |
|
|
"metadata": {}, |
|
|
"output_type": "display_data" |
|
|
}, |
|
|
{ |
|
|
"data": { |
|
|
"text/markdown": [ |
|
|
"Total Heat Loss Coefficient (HLC) per heated floor area: $H'' = \\frac{H_{tr} + H_{inf} +H_{ve} }{A_{fl}} = 0.93$ $W/(K m^2)$" |
|
|
], |
|
|
"text/plain": [ |
|
|
"<IPython.core.display.Markdown object>" |
|
|
] |
|
|
}, |
|
|
"metadata": {}, |
|
|
"output_type": "display_data" |
|
|
} |
|
|
], |
|
|
"source": [ |
|
|
"# Heated floor area \"BRA\"\n", |
|
|
"A_fl = 100 # m^2\n", |
|
|
"# Facade area\n", |
|
|
"A_ew = 102.5 # m^2\n", |
|
|
"# Inddor height\n", |
|
|
"H_b = 2.5 # m^2\n", |
|
|
"\n", |
|
|
"# Windows\n", |
|
|
"A_wi = 20 # m^2\n", |
|
|
"#A_gl = A_wi * 0.8 # is the transparent part if frame factor is 20% of window\n", |
|
|
"\n", |
|
|
"# U-values\n", |
|
|
"U_ew = 0.18 # W/(m^2 K\n", |
|
|
"U_rf = 0.13 # W/(m^2 K)\n", |
|
|
"U_fl = 0.15 # W/(m^2 K)\n", |
|
|
"U_wi = 1.2 # W/(m^2 K)\n", |
|
|
"\n", |
|
|
"# Ventilation rate\n", |
|
|
"V_ve = 1.2 # m^3/(h m^2) ex. 1.25 m^3/(h m^2) equals a fresh air exhange rate of 0.5 per hour if H_b = 2.5 meter\n", |
|
|
"\n", |
|
|
"# Heat recovery efficiency \n", |
|
|
"n_T = 0.7 # 70 % efficiency\n", |
|
|
"\n", |
|
|
"# Infiltration rate at 50 Pa pressure difference\n", |
|
|
"n_50 = 2.5 # air change per hour [/h] at 50 Pa\n", |
|
|
"n_inf = n_50 * 0.07 # ach [/h], where 0.07 is a multiplier for balanced ventilation ref. NS 3031:2014\n", |
|
|
"\n", |
|
|
"# Transmission losses\n", |
|
|
"H_tr = U_ew * (A_ew - A_wi) + A_wi * U_wi + A_fl * U_rf + A_fl * U_fl\n", |
|
|
"\n", |
|
|
"#Infiltration\n", |
|
|
"H_inf = n_inf * A_fl * H_b *0.33 \n", |
|
|
"\n", |
|
|
"# Ventilation\n", |
|
|
"H_ve = V_ve * A_fl * 0.33 * (1-n_T)\n", |
|
|
"\n", |
|
|
"# Heat loss coiefficent\n", |
|
|
"HLC = (H_tr + H_inf + H_ve)/A_fl\n", |
|
|
"\n", |
|
|
"# Display formulas in notebook\n", |
|
|
"display(md(\"### Resulting HLC\"))\n", |
|
|
"display(md(\"Transission losses: $H_{tr} = \\\\sum{U_i \\cdot A_i} = %.1f$ $W/K$ , where $i$ is each envelope surface part\"%(H_tr)))\n", |
|
|
"display(md(\"Infiltration losses: $H_{inf} = 0.33 \\cdot n_{inf} \\cdot H_b \\cdot A_{fl} = %.1f$ $W/K$ , where $n_{inf} = 0.07 n_{50}$ for balanced ventilation ref. NS 3031:2014\"%(H_inf)))\n", |
|
|
"display(md(\"Ventilation losses: $H_{ve} = 0.33 \\cdot \\dot V \\cdot (1-\\\\eta_{T}) \\cdot A_{fl} = %.1f$ $W/K$\"%(H_ve)))\n", |
|
|
"display(md(\"Total Heat Loss Coefficient (HLC) per heated floor area: $H'' = \\\\frac{H_{tr} + H_{inf} +H_{ve} }{A_{fl}} = %.2f$ $W/(K m^2)$\"%(HLC)))" |
|
|
] |
|
|
}, |
|
|
{ |
|
|
"cell_type": "markdown", |
|
|
"metadata": {}, |
|
|
"source": [ |
|
|
"# Calculate the indoor heating need at design conditions\n", |
|
|
"\n", |
|
|
"Under stationary conditions and by neglecting heat gains from internal loads and solar irradiance." |
|
|
] |
|
|
}, |
|
|
{ |
|
|
"cell_type": "code", |
|
|
"execution_count": 701, |
|
|
"metadata": {}, |
|
|
"outputs": [ |
|
|
{ |
|
|
"data": { |
|
|
"text/markdown": [ |
|
|
"### Resulting heating need at design dimensions (-20, 20 ℃)" |
|
|
], |
|
|
"text/plain": [ |
|
|
"<IPython.core.display.Markdown object>" |
|
|
] |
|
|
}, |
|
|
"metadata": {}, |
|
|
"output_type": "display_data" |
|
|
}, |
|
|
{ |
|
|
"data": { |
|
|
"text/markdown": [ |
|
|
"Specific heating need: $q_{tot} = H'' (T_i - T_e)= 37.3$ $W/m^2$ , or in total: $3.7$ $kW$" |
|
|
], |
|
|
"text/plain": [ |
|
|
"<IPython.core.display.Markdown object>" |
|
|
] |
|
|
}, |
|
|
"metadata": {}, |
|
|
"output_type": "display_data" |
|
|
}, |
|
|
{ |
|
|
"data": { |
|
|
"text/markdown": [ |
|
|
"Next, we can separate the heating need between ventilation and room heating:" |
|
|
], |
|
|
"text/plain": [ |
|
|
"<IPython.core.display.Markdown object>" |
|
|
] |
|
|
}, |
|
|
"metadata": {}, |
|
|
"output_type": "display_data" |
|
|
}, |
|
|
{ |
|
|
"data": { |
|
|
"text/markdown": [ |
|
|
"Temperature after heat recovery: $\\eta_T = \\frac{T_{HR}-T_e}{T_i - T_e} \\rightarrow T_{HR} = \\eta_T(T_i - T_e) + T_e = 8.0$ $℃$" |
|
|
], |
|
|
"text/plain": [ |
|
|
"<IPython.core.display.Markdown object>" |
|
|
] |
|
|
}, |
|
|
"metadata": {}, |
|
|
"output_type": "display_data" |
|
|
}, |
|
|
{ |
|
|
"data": { |
|
|
"text/markdown": [ |
|
|
"Specific ventilation heating need: $q_{ve} = 0.33(20℃-T_{HR}) \\dot V_{ve} = 4.8$ $W/m^2$ , to reach supply temperature of 20℃" |
|
|
], |
|
|
"text/plain": [ |
|
|
"<IPython.core.display.Markdown object>" |
|
|
] |
|
|
}, |
|
|
"metadata": {}, |
|
|
"output_type": "display_data" |
|
|
}, |
|
|
{ |
|
|
"data": { |
|
|
"text/markdown": [ |
|
|
"The rest will need to be covered by room heating: $q_{sh} = 32.5$ $W/m^2$" |
|
|
], |
|
|
"text/plain": [ |
|
|
"<IPython.core.display.Markdown object>" |
|
|
] |
|
|
}, |
|
|
"metadata": {}, |
|
|
"output_type": "display_data" |
|
|
} |
|
|
], |
|
|
"source": [ |
|
|
"# Calculate the indoor heating need \n", |
|
|
"# neglecting indoor heat gains and sun\n", |
|
|
"\n", |
|
|
"# Outdoor temperature\n", |
|
|
"T_e = -20\n", |
|
|
"# Indoor temperature\n", |
|
|
"T_i = 20\n", |
|
|
"\n", |
|
|
"# Specific heating need at design dimensions \"DUT\"\n", |
|
|
"q_dim = HLC * (T_i-T_e)\n", |
|
|
"\n", |
|
|
"# First, we calculate the temperature after heat recovery \n", |
|
|
"T_HR = n_T * (T_i - T_e) + T_e\n", |
|
|
"\n", |
|
|
"# If we want to supply 20 degrees air to the room, the ventilation heating need is\n", |
|
|
"q_ve = 0.33 * V_ve *(20 - T_HR)\n", |
|
|
"\n", |
|
|
"# And the resulting specific heating need to be supplied by room heating is\n", |
|
|
"q_sh = q_dim - q_ve\n", |
|
|
"\n", |
|
|
"# Display formulas in notebook\n", |
|
|
"display(md(\"### Resulting heating need at design dimensions (%.0f, %.0f ℃)\"%(T_e,T_i)))\n", |
|
|
"display(md(\"Specific heating need: $q_{tot} = H'' (T_i - T_e)= %.1f$ $W/m^2$ , or in total: $%.1f$ $kW$\"%(q_dim, q_dim * A_fl / 1000)))\n", |
|
|
"display(md(\"Next, we can separate the heating need between ventilation and room heating:\"))\n", |
|
|
"display(md(\"Temperature after heat recovery: $\\\\eta_T = \\\\frac{T_{HR}-T_e}{T_i - T_e} \\\\rightarrow T_{HR} = \\\\eta_T(T_i - T_e) + T_e = %.1f$ $℃$\"%(T_HR)))\n", |
|
|
"display(md(\"Specific ventilation heating need: $q_{ve} = 0.33(20℃-T_{HR}) \\dot V_{ve} = %.1f$ $W/m^2$ , to reach supply temperature of 20℃\"%(q_ve)))\n", |
|
|
"display(md(\"The rest will need to be covered by room heating: $q_{sh} = %.1f$ $W/m^2$\"%(q_sh)))" |
|
|
] |
|
|
}, |
|
|
{ |
|
|
"cell_type": "markdown", |
|
|
"metadata": {}, |
|
|
"source": [ |
|
|
"## Dynamic heat balance\n", |
|
|
"The simplified dynamic heat balance can be expressed as: \n", |
|
|
"\n", |
|
|
"$C'' \\frac{dT_i}{dt} = q''_{int} + q''_{sol} - H'' (T_i - T_e)$\n", |
|
|
"\n", |
|
|
"Where: \n", |
|
|
"* Total heat capacity of the building envelope, $C''$ (to the inside) is calculated\n", |
|
|
"* Internal heat gains $q''_{int}$ from occupancy, lighting and plug-loads are specified in operating hours 8:00 to 16:00\n", |
|
|
"* For the calculation of solar gains $q''_{sol}$ Window glazing g-value, frame factor and simplified horizon factor is introduced.\n", |
|
|
"\n", |
|
|
"New asumptions representing summer conditions:\n", |
|
|
"* The average exterior solar gains hitting the windows is $300 W/m^2$ in operating hours 8-16\n", |
|
|
"* The building is ventilated with 4 air changes per hour ($10 m^3/hm^2$)\n", |
|
|
"* Outdoor temperature $T_e = 20 ℃$. " |
|
|
] |
|
|
}, |
|
|
{ |
|
|
"cell_type": "code", |
|
|
"execution_count": 669, |
|
|
"metadata": {}, |
|
|
"outputs": [ |
|
|
{ |
|
|
"data": { |
|
|
"text/markdown": [ |
|
|
"### Heat gains and losses" |
|
|
], |
|
|
"text/plain": [ |
|
|
"<IPython.core.display.Markdown object>" |
|
|
] |
|
|
}, |
|
|
"metadata": {}, |
|
|
"output_type": "display_data" |
|
|
}, |
|
|
{ |
|
|
"data": { |
|
|
"text/markdown": [ |
|
|
"Internal heat gains in operating hours: $q_{int} = 9.0$ $W/m^2$" |
|
|
], |
|
|
"text/plain": [ |
|
|
"<IPython.core.display.Markdown object>" |
|
|
] |
|
|
}, |
|
|
"metadata": {}, |
|
|
"output_type": "display_data" |
|
|
}, |
|
|
{ |
|
|
"data": { |
|
|
"text/markdown": [ |
|
|
"Solar heat gains to the indoors in operating hours: $q_{sol} = 9.6$ $W/m^2$" |
|
|
], |
|
|
"text/plain": [ |
|
|
"<IPython.core.display.Markdown object>" |
|
|
] |
|
|
}, |
|
|
"metadata": {}, |
|
|
"output_type": "display_data" |
|
|
}, |
|
|
{ |
|
|
"data": { |
|
|
"text/markdown": [ |
|
|
"Combined specific heat losses to the outside: $H'' = 4.1$ $W/(K m^2)$ , given the newly introduced air exchange rate of 4.0 per hour." |
|
|
], |
|
|
"text/plain": [ |
|
|
"<IPython.core.display.Markdown object>" |
|
|
] |
|
|
}, |
|
|
"metadata": {}, |
|
|
"output_type": "display_data" |
|
|
} |
|
|
], |
|
|
"source": [ |
|
|
"# We consider the same 100 m^2 case building in summer conditions, in the example below.\n", |
|
|
"\n", |
|
|
"# Heat capacity of envelope parts\n", |
|
|
"C_ew = 3 # Wh/(m^2 K) (82.5 m^2)\n", |
|
|
"C_fl = 41 # Wh/(m^2 K) (100 m^2)\n", |
|
|
"C_rf = 5 # Wh/(m^2 K) (100 m^2)\n", |
|
|
"\n", |
|
|
"# Temperature at 08:00 in the morning\n", |
|
|
"T_init = 20 # Celsius\n", |
|
|
"\n", |
|
|
"# Air exchange (summer)\n", |
|
|
"V_ve = 4 # ACH [/h]\n", |
|
|
"\n", |
|
|
"# Average solar irradiance flux from 08:00 to 16:00\n", |
|
|
"q_ver = 300 # W/m2\n", |
|
|
"\n", |
|
|
"# Window glazing\n", |
|
|
"f_gl = 0.25 # Solar heat glazing factor and solar shading \n", |
|
|
"f_wi = 0.20 # Frame factor e.g. 0.2 is 20 % frame\n", |
|
|
"f_horiz = 0.80 # Horizon factor\n", |
|
|
"\n", |
|
|
"# Outdoor temp\n", |
|
|
"T_e = 20 # Celsius\n", |
|
|
"\n", |
|
|
"# Internal heat gains 08:00 to 16:00\n", |
|
|
"q_light = 4 # W/m2\n", |
|
|
"q_equip = 3 # W/m2\n", |
|
|
"q_occup = 2 # W/m2\n", |
|
|
"\n", |
|
|
"# Internal heat gains\n", |
|
|
"q_int = q_light + q_equip + q_occup\n", |
|
|
"\n", |
|
|
"# Solar heat gains to the indoors\n", |
|
|
"q_sol = q_ver * A_wi * (1-f_wi) * f_gl * f_horiz / A_fl # W/(m^2)\n", |
|
|
"\n", |
|
|
"# We neglect infiltration losses by considering them part of the 4 ACH summer air exchange\n", |
|
|
"HLC = (H_tr + 0.34 * V_ve * A_fl * H_b) / A_fl\n", |
|
|
"\n", |
|
|
"# Display results\n", |
|
|
"display(md(\"### Heat gains and losses\"))\n", |
|
|
"display(md(\"Internal heat gains in operating hours: $q_{int} = %.1f$ $W/m^2$\"%(q_int)))\n", |
|
|
"display(md(\"Solar heat gains to the indoors in operating hours: $q_{sol} = %.1f$ $W/m^2$\"%(q_sol)))\n", |
|
|
"display(md(\"Combined specific heat losses to the outside: $H'' = %.1f$ $W/(K m^2)$ , given the newly introduced air exchange rate of %.1f per hour.\"%(HLC, V_ve)))" |
|
|
] |
|
|
}, |
|
|
{ |
|
|
"cell_type": "code", |
|
|
"execution_count": 689, |
|
|
"metadata": {}, |
|
|
"outputs": [ |
|
|
{ |
|
|
"data": { |
|
|
"text/markdown": [ |
|
|
"### Building envelope heat capacity" |
|
|
], |
|
|
"text/plain": [ |
|
|
"<IPython.core.display.Markdown object>" |
|
|
] |
|
|
}, |
|
|
"metadata": {}, |
|
|
"output_type": "display_data" |
|
|
}, |
|
|
{ |
|
|
"data": { |
|
|
"text/markdown": [ |
|
|
"Normalized internal heat capacity of building envelope $C'' = \\sum{C_i \\cdot A_i} \\cdot \\frac{1}{A_{fl}} = 48.5$ $Wh/(K m^2)$" |
|
|
], |
|
|
"text/plain": [ |
|
|
"<IPython.core.display.Markdown object>" |
|
|
] |
|
|
}, |
|
|
"metadata": {}, |
|
|
"output_type": "display_data" |
|
|
}, |
|
|
{ |
|
|
"data": { |
|
|
"text/markdown": [ |
|
|
"Time constant: $\\tau = \\frac{C''}{H''} = 11.9$ $h$" |
|
|
], |
|
|
"text/plain": [ |
|
|
"<IPython.core.display.Markdown object>" |
|
|
] |
|
|
}, |
|
|
"metadata": {}, |
|
|
"output_type": "display_data" |
|
|
} |
|
|
], |
|
|
"source": [ |
|
|
"# Normalized envelope heat capacity\n", |
|
|
"C_env = (C_ew * (A_ew-A_wi) + C_fl * A_fl + C_rf * A_fl)/A_fl\n", |
|
|
"\n", |
|
|
"# Building time constant aka NS3031:2014\n", |
|
|
"tau = C_env / HLC # hour\n", |
|
|
"\n", |
|
|
"# Display results\n", |
|
|
"display(md(\"### Building envelope heat capacity\"))\n", |
|
|
"display(md(\"Normalized internal heat capacity of building envelope $C'' = \\\\sum{C_i \\cdot A_i} \\cdot \\\\frac{1}{A_{fl}} = %.1f$ $Wh/(K m^2)$\"%(C_env)))\n", |
|
|
"display(md(\"Time constant: $\\\\tau = \\\\frac{C''}{H''} = %.1f$ $h$\"%(tau)))" |
|
|
] |
|
|
}, |
|
|
{ |
|
|
"cell_type": "code", |
|
|
"execution_count": 670, |
|
|
"metadata": {}, |
|
|
"outputs": [ |
|
|
{ |
|
|
"data": { |
|
|
"text/markdown": [ |
|
|
"### Solving the dynamic heat balance manually to aquire indoor temperature" |
|
|
], |
|
|
"text/plain": [ |
|
|
"<IPython.core.display.Markdown object>" |
|
|
] |
|
|
}, |
|
|
"metadata": {}, |
|
|
"output_type": "display_data" |
|
|
}, |
|
|
{ |
|
|
"data": { |
|
|
"text/markdown": [ |
|
|
"Infinite indoor temperature: $T_\\infty = 24.6$ $°C$ , if solar and internal gains are operating indefinitely" |
|
|
], |
|
|
"text/plain": [ |
|
|
"<IPython.core.display.Markdown object>" |
|
|
] |
|
|
}, |
|
|
"metadata": {}, |
|
|
"output_type": "display_data" |
|
|
} |
|
|
], |
|
|
"source": [ |
|
|
"# Infinite indoor temperature if solar and internal gains are operating indefinitely\n", |
|
|
"T_infinite = (q_int + q_sol + HLC * T_e)/HLC\n", |
|
|
"\n", |
|
|
"# Display results\n", |
|
|
"display(md(\"### Solving the dynamic heat balance manually to aquire indoor temperature\"))\n", |
|
|
"display(md(\"Infinite indoor temperature: $T_\\infty = %.1f$ $°C$ , if solar and internal gains are operating indefinitely\"%(T_infinite)))" |
|
|
] |
|
|
}, |
|
|
{ |
|
|
"cell_type": "markdown", |
|
|
"metadata": {}, |
|
|
"source": [ |
|
|
"$T_i(t) = T_\\infty + (T_i(0) - T_\\infty) \\cdot \\exp{(\\frac{-t}{\\tau})}$\n", |
|
|
"\n", |
|
|
"The building can be considered to be in thermal balance until 08:00 (as $T_i = T_e$), thus we can plot the temperature development from $T_i(0)$ at 8:00 in the morning, applying the heat gains indefinetely:" |
|
|
] |
|
|
}, |
|
|
{ |
|
|
"cell_type": "code", |
|
|
"execution_count": 682, |
|
|
"metadata": {}, |
|
|
"outputs": [ |
|
|
{ |
|
|
"data": { |
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXhV1b3/8feXJCSBQAIBAiRAmAWiouAATuBQtXpr7dXa1lpbp9ra/pzaq/XWtrb1ub221d721qotVqtYtHW2dboVcGhVEqYwySRDIEyBzHPy/f1xDhowJIeQk5OT/Xk9T56cs3P22d+9JJ9s11l7LXN3REQkOHrFugAREelaCn4RkYBR8IuIBIyCX0QkYBT8IiIBkxjrAiIxaNAgz83N7dC+VVVV9O3bt3MLijNqA7VB0M8fgtkGBQUFe9x98MHb4yL4c3Nzyc/P79C+CxYsYNasWZ1bUJxRG6gNgn7+EMw2MLPNrW1XV4+ISMAo+EVEAkbBLyISMAp+EZGAUfCLiASMgl9EJGAU/CIiARMX4/g77OXbmbrmLfgwI9aVxNTU0lK1QcDbIOjnD3HcBkOPhvN/1qlvqSt+EZGA6dlX/Of/jKWpwbtb72BLA3jH4sGC3gZBP39QG7SkK34RkYDp2Vf8IiLdhLtT29BMRV0DlbWNVNQ2Uln38ffK2gaq6puorGukqq6RyvD2G88ez5Th6Z1ai4JfRCRCTc1OWU0DpdX1lNY0UFbTQHn4e1l1+HltAxW1jZTXNlBe0/jx85oGGpvbX+M8KcHom5xI396JpCUnUlXX1OnnoeAXkUBqbnZKaxooqaxjb1V96Ku6nn1V9eytamBvVR37qj8O+X1V9ZTXNrb5nqlJCfRPTaR/ShL9U5MYlNab0YP60j81kX4pSfRLSaRf8seP08KP05ITSUtJpG9yAsmJCVE/dwW/iPQYjU3N7K2qZ1dFHXsq69hdUcfu8PdVG2p5aN277K2qZ09lPXur6jjUBXhaciID+iYxoE9vMvr0JndQXzJSk8jo05uMPkmhr9Te9E9NIr3FV+/E+PjYVMEvIt2eu1Ne00hxeQ3FpbXsKK9lZ3ktO8vr2FVey86K0OM9lXV4K2GelpxI34RmcpKbGTmwD8eNHMCgtN4M7Bv6GpSWzIA+occD+iZ1yVV3LCn4RSTmahuaKC6rZdu+GraVVrOtNPS4uKyGHWW1FJfVUtPwyb7uQWm9GdIvhaz+yeQNT2dIv2QG909hcFoyg/slM6RfMoPSkkntnRBeiGVmDM6u+1Hwi0jUNTY1U1xWy9a91WzdV82WvdVs3VvD1n2h73sq6w54vRlk9UthWEYKk4b1Z/ZRQxiWnsKw9FSGpqcwND0U7vHStdLdKPhFpFM0NDVTtK+GTSVVbNpTxeaSajaVhL5v3Vt9wIiWhF7G8IwURg7sw1lHDSF7QCrZGakMz0glZ0Ao3JMSFOrRouAXkcNSVt3A+t2VbNhdycbdVeHvlWwuOTDc05ITGZXZh8nD+3N+3lBGZfZhxMA+jBjQh2HpKSQq2GNGwS8irSqraWDdzgrW7qxk7c4K1u2qYN3OSnZVfNwtk5Rg5Gb2ZdyQND41ZShjBvVl9KC+5A7qS2bf3phZDM9ADkXBLxJwTc3Oh3uqWF1czpod5awurmBNcTnby2o/ek1qUgLjs9I4bfxgJmSlMW5IGmMHp5EzIFVX7nFIwS8SIPWNzazdWcGKbWUUbitjxbYy1uyooK6xGQj1vY8bnMYJowdy1ND+TByaxvgh/cjOSKVXL1299xQKfpEeqqnZWb+rkqVb9/Hyyjp+Wfg2H+yooL4pFPL9khOZkt2fy08axaRh/Zg0rD/js9J6/Bh2UfCL9Bi7K+pYsmUfS7aWsnRLKYXbyqisC00xkJoIx41K5Gun5JKXnU5edjqjBvbRVXxAKfhF4pB7qF8+f9M+Fm3ay6JNe9lUUg1AYi9j0rD+XHxcNlNHZDB1ZAabVyzizNknx7hq6S4U/CJxoLnZWburgn9tKOHdjSXkb9pHSVU9AAP6JDFt1EC+dNJIjh85gLzsdFKSDuyu2arRNdKCgl+kG3J3Nuyu4l8bS/jXhj28u3Eve8NBnzMglTMmDuaE3IGckDuAMYPS1GUjh0XBL9JN7Kuq5+31e3hz7W7eWreHHeWh4ZTD01OYNXEwM8ZkMmNsJjkD+sS4Uol3Cn6RGGlqdpZu3cfCD3azcN0elheV4g7pqUmcOm4Qp44fxIwxmYzK7KMboaRTRS34zWwE8CdgKNAMPOTu/9Pi598Bfg4Mdvc90apDpDupqG3grXV7+L9VO5n/wS72VTfQy2DqiAxuPGs8p08YzLE5GSSo60aiKJpX/I3Are6+2Mz6AQVm9rq7rwr/UTgH2BLF44t0C9tKa3h95Q7+sWYX724soaHJyeiTxOyJQzhr0hBOGzeY9D5JsS5TAiRqwe/uxUBx+HGFma0GsoFVwH3AfwDPR+v4IrG0paSal1cU8/cVO1i2tRSAsYP7ctUpozlrUhbHj8zQVAcSM+atLVfT2QcxywXeBPKAWcBZ7n6jmW0CprfW1WNm1wHXAWRlZU2bN29eh45dWVlJWlpah/btKdQGXdMGO6qaWbSjkfydTWwuD90dO7p/L6YPTWBaViJD+8Yu6PVvIJhtMHv27AJ3n37w9qgHv5mlAQuBu4FXgPnAp9y9rK3gb2n69Omen5/foeOHVt2Z1aF9ewq1QfTaYFd5LS8s287zS7dTuK0MgONHZvDpo4dx7pShjBjYPUbg6N9AMNvAzFoN/qiO6jGzJOBpYK67P2NmRwOjgWXhUQo5wGIzO9Hdd0SzFpHOUlnXyCsrdvD80m28s34PzQ5HZ6fz/QsmccExwxiWnhrrEkXaFM1RPQbMAVa7+70A7l4IDGnxmk1EcMUvEmvNzc67G0t4Mn8rr67cQW1DMyMGpnLD7HFcNDWbcUOC1YUg8S2aV/ynAFcAhWa2NLztDnf/exSPKdKpistq+Gt+EX8pKGLL3mr6pSRyybQcLj4um+NHDtD4eolL0RzV8zbQ5m+Fu+dG6/giHdXY1Mz/rd7Fk4u2sHDtbpodZozJ5JZzJnBe3tBPzIMjEm90565I2J7KOua9v4W5722huKyWof1T+OascVw6PYdRmX1jXZ5Ip1HwS6C5O0u3lvKnf23mb8uLqW9q5rTxg7jrM1M4a1KW7qCVHknBL4HU0NTMS8u388d3NrG8qIy05ES+dNJIrpgxirGD9UGt9GwKfgmUyrpG5r2/hYff/pDtZbWMG5LGTz6bx8XHZZOWrF8HCQb9S5dAKK1r5p5X1vD4u5spr23kpNED+enFecyaMERz2UvgKPilR9tSUs39C9bz1/wamtnAeXlDue70sUwdkRHr0kRiRsEvPdLmkir+9431PLNkGwm9jNNzEvnBZaeSO0ijc0QU/NKjtAz8xF7GV2aM4vozxrJ68bsKfZEwBb/0CFtKqvn1G+t4Nhz4V87I5fozxjCkfwoAq2Ncn0h3ouCXuFZSWcdv3ljP3Pc208s+Gfgi8kkKfolL1fWNPPz2hzywcCPV9Y1cdsIIbjp7AlkKfJF2KfglrjQ2NfOXgiLue30tuyrqOGdyFredN5FxQ/rFujSRuKHgl7jx1rrd3PXiKtbvqmTaqAHcf/nxTM8dGOuyROKOgl+6va17q/np31bx6sqdjMrswwNfnsa5U7I0JbJIByn4pduqbWjidws28MDCDfQy47vnTuTqU0drWmSRI6Tgl27H3Xl15U5+8tIqtpXWcOExw7jj05MYnqElDUU6g4JfupVtpTV8/9lC5n+wm4lZ/fjztSczY2xmrMsS6VEU/NItNDU7j/1rE/e8+gEAd144mStnjCIxoVdsCxPpgRT8EnNrd1Zw29PLWbKllNMnDObuz+YxYmCfWJcl0mMp+CVm6hqbuH/+Bu5fsJ605ER+ddlULpo6XKN1RKJMwS8xsWJbGTc/uZR1uyr57NTh3HnhZDLTkmNdlkggKPilSzU1Ow8s3MB9r68lM603f/zqCcw+akisyxIJFAW/dJktJdXc8tRS8jfv44Kjh3H3xXlk9Okd67JEAkfBL1Hn7vwlv4i7XlxJLzPuu+xYPjs1W335IjGi4Jeo2ltVz+1PL+e1VTs5ecxAfvn5qWTrRiyRmFLwS9QUbN7Lt55YQkllPXd8+iiuOXWMFjYX6QYU/NLp3J05b3/Iz15ew7CMFJ755kzystNjXZaIhCn4pVOV1TTw3b8s47VVO/nU5Cx+fumxpKcmxbosEWmhzeA3sxnAl4HTgGFADbAC+BvwuLuXRb1CiRuFRWV884kCiktr+f4Fk7j61NH6AFekGzpk8JvZy8B24HngbmAXkAJMAGYDz5vZve7+QlcUKt2Xu/PE+1u464VVZKb15smvz2DaqAGxLktEDqGtK/4r3H3PQdsqgcXhr1+a2aCoVSZxoaGpmR++sJIn3tvCGRMGc99lUxnYV2PzRbqztoI/w8wmuvs7LTea2WnAdnff0MofBgmQvVX1fOPxAt77cC/XnzGW7547kQSN2hHp9tqa8/ZXQEUr22vCP5MA+2BHBRf99m2WbC3lvsuO5fbzj1Loi8SJtq74c919+cEb3T3fzHKjVpF0e6+v2slN85bQJzmRJ687meNGqj9fJJ60dcWf0sbP2r310sxGmNl8M1ttZivN7Mbw9p+Y2XIzW2pmr5nZ8MMtWmLD3bl/wXqueyyfMYPTeOFbpyj0ReJQW8G/yMyuPXijmV0NFETw3o3Are4+CTgZuMHMJgM/d/dj3H0q8BLwgw7ULV2ssamZ255ezj2vfMCFxwznqa/PYFi6pl4QiUdtdfXcBDxrZpfzcdBPB3oDF7f3xu5eDBSHH1eY2Wog291XtXhZX8A7Urh0ner6Rr71xBLeWLOL/3fmOG4+Z4LG54vEMXNvO3fNbDaQF3660t3fOOyDhD4TeBPIc/dyM7sb+ApQBsx2992t7HMdcB1AVlbWtHnz5h3uYQGorKwkLS2tQ/v2FEfSBhX1zn0FtXxY1swVk3tz5sj4vAs36P8Ogn7+EMw2mD17doG7Tz94e7vBf6TMLA1YCNzt7s8c9LPvASnu/sO23mP69Omen5/foeMvWLCAWbNmdWjfnqKjbbB1bzVXPvw+RaU1/PoLx3Fe3tDOL66LBP3fQdDPH4LZBmbWavAfso/fzC41s+fM7Fkzu6yDB00CngbmHhz6YU8A/96R95boWrW9nM/97p/sqaxj7jUnxXXoi8iB2urjvw04Mfx4EfDk4byxhTqB5wCr3f3eFtvHu/u68NPPAGsO530l+v65YQ9f/1MBaSmJzP3GTCZk9Yt1SSLSidoK/seBP4Uf/6UD730KcAVQaGZLw9vuAK42s4lAM7AZuL4D7y1RMv+DXXz9sQJyM/vw6FUnauSOSA90yOB391+ZWV9CnwNUHu4bu/vbQGtDP/5+uO8lXeP1VTu5Ye5iJg7tx2NXn6j1cEV6qLZm5zR3r2pr5/BrNByzB3i5sJhv/3kJednpPHrViZpDX6QHa+sGrvlm9m0zG9lyo5n1NrMzzexR4Mrolidd4YVl2/nWn5dw7IgMHrtaoS/S07XVx38ecBXwZzMbDZQSmsYhAXgNuM/dl7axv8SBpwuK+O5flzE9dyB//OoJ9E3WomwiPV1bffy1wP3A/eFhmYOAGncv7ariJLqeXLSF258pZObYTH7/len06a3QFwmCiH7T3b2B8PQL0jM8tWgrtz1dyBkTBvPgFdNISUqIdUki0kV0iRdAf1tezO3PLOf0CYN56CvTSE5U6IsESVsf7koPNH/NLm56cgnTRg3gwS8r9EWCKKLgN7NRZnZ2+HGqmelWzjj07sYSrn+8gIlD+zHnqyeQ2luhLxJE7QZ/eE7+vwIPhjflAM9FsyjpfMuLSrnm0XxyBqTy6NdOpH+KhmyKBFUkV/w3EJp+oRwgPM/OkGgWJZ1rW0UzX3n4fTL6JDH3mpPJTEuOdUkiEkORBH+du9fvf2JmiWjxlLixuaSKn+fX0juhF3OvOYmh6W2tqCkiQRBJ8C80szuAVDM7h9CEbS9GtyzpDHsq67hizvs0NjuPX3MSozL7xrokEekGIgn+24DdQCHwdUKTrH0/mkXJkattaOKaR/PZWV7LzdNSNLWyiHykzXH8ZtYLWO7uecDvu6YkOVLNzc7NTy5lWVEpv7v8eFL2fBDrkkSkG2nzit/dm4FlB0/UJt3bz15Zw8srdvCfn57EeXnDYl2OiHQzkdy5OwxYaWbvAx9N0+zun4laVdJhj/1rEw+9uZErZ4zi6lNHx7ocEemGIgn+u6JehXSKN9bs5IcvrOSso4bwg3+bQmj1SxGRA7Ub/O6+sCsKkSOzYlsZ33piCZOH9+fXXzyOhF4KfRFpXbvBb2YVfDxuvzeQBFS5e/9oFiaRKy6r4apHFpGRmsTDV2pOfRFpWyRX/AeMAzSzzwInRq0iOSy1DU1c/1gB1fVNPP2NmQzprxu0RKRthz07p7s/B5wZhVrkMLk7P3h+BcuKyvjl549l4lCN1ReR9kXS1fO5Fk97AdPRlA3dwtz3tvBUfhHfPnMc504ZGutyRCRORNIZ/G8tHjcCm4CLolKNRCx/017uenElsycO5qazJ8S6HBGJI5EE/x/c/Z2WG8zsFGBXdEqS9uwsr+UbcxeTnZHKr76gETwicngi6eP/TYTbpAvUNTbxjccLqKpr5MErppOeqnn1ReTwHPKK38xmADOBwWZ2S4sf9Qe0dFOM3PXiKhZvKeX+y4/Xh7ki0iFtdfX0BtLCr2mZMOXAJdEsSlo37/0tPPHeFr4xayyfPlpz8IhIxxwy+MN37C40s0fcfXMX1iStWF1czg9eWMlp4wfxnU9NjHU5IhLHIvlwt9rMfg5MAT66O8jdNZa/i1TXN/LtPy8hPTWJ+y6bqg9zReSIRPLh7lxgDTCa0IRtm4BFUaxJDvLjF1exYXclv7psKoO0Xq6IHKFIgj/T3ecADe6+0N2vAk6Ocl0S9uKy7cxbtJVvnDGWU8YNinU5ItIDRNLV0xD+XmxmFwDbgZzolST7bd1bzR3PFHL8yAxuPkc3aYlI54gk+H9qZunArYTG7/cHbo5qVUJDUzPf+vMSMPifLxxHUsJhT6skItKq9tbcTQDGu/tLQBkwu0uqEn7x2gcs2xoarz9iYJ9YlyMiPUh7a+42AR1aYtHMRpjZfDNbbWYrzezG8Pafm9kaM1tuZs+aWUZH3r8ne3Ptbh5cuJEvnTRS4/VFpNNF0n/wTzP7XzM7zcyO3/8VwX6NwK3uPonQh8E3mNlk4HUgz92PAdYC3+tw9T3Q7oo6bnlqKROy0vjBhZNjXY6I9ECR9PHPDH//cYttTjtz8rt7MVAcflxhZquBbHd/rcXL3kV3AX/E3fn+c4WU1zYy95qTSUnSzBgi0vnMPfpT65tZLvAmoSv98hbbXwSedPfHW9nnOuA6gKysrGnz5s3r0LErKytJS0vr0L5d7Z/bG3loeR2XTezN+aM7b/K1eGqDaAl6GwT9/CGYbTB79uwCd5/+iR+4e5tfQBYwB3g5/HwycHV7+7XYPw0oAD530Pb/BJ4l/Menra9p06Z5R82fP7/D+3alHWU1fsyPXvWLf/u2NzY1d+p7x0sbRFPQ2yDo5+8ezDYA8r2VTI2kj/8R4FVgePj5WuCmSP7amFkS8DQw192fabH9SuBC4PJwcYHm7tzxTCF1jU384tJjNSWDiERVJME/yN2fApoB3L0RaGpvJzMzQv+nsNrd722x/TzgNuAz7l7doap7mL8WFPGPNbv4j3OPYszgYP2vqIh0vUg+3K0ys0zC6+ya2cmExvS35xTgCqDQzJaGt90B/BpIBl4P/W3gXXe//nAL7ymKy2r48YurOHH0QL46MzfW5YhIAEQS/LcALwBjzewdYDARjMRx97eB1vos/n5YFfZg7s5tTxfS5M4vLjmWXuriEZEu0G7wu/tiMzsDmEgoyD9w94Z2dpMIzFu0lTfX7uYnF01hZKbuzhWRrtFu8JtZCvBN4FRC3T1vmdkD7l4b7eJ6sqJ91fz0pVXMHJvJ5SeNinU5IhIgkXT1/Amo4OMF1r8IPAZcGq2iejp353vPFAJwzyXHqItHRLpUJME/0d2PbfF8vpkti1ZBQfDCsu28tW4PP75oCjkD1MUjIl0rkuGcS8IjeQAws5OAd6JXUs9WVtPAT15azbE56eriEZGYiOSK/yTgK2a2Jfx8JLDazAoB99BkaxKhe15Zw96qOh752gm6UUtEYiKS4D8v6lUExOIt+3ji/S1cdcpo8rLTY12OiARUJMM5N5vZAGBEy9e7++JoFtbTNDY185/PrmBo/xQtoygiMRXJcM6fAF8FNhC+e5cIpmWWAz3yz02sLi7ngS9PIy05kv/REhGJjkgS6PPAWHevj3YxPdW20hrufX0tZx01hHOnZMW6HBEJuEhG9awAtDziEbjrhZW4w10XTSE8P5GISMxEcsX/X4SGdK4A6vZvdPcOrcUbNK+v2slrq3Zy+/lHacy+iHQLkQT/o8B/A4WEp2aWyFTXN/LD51cwMasfV586OtbliIgAkQX/Hnf/ddQr6YEeWLCB7WW1/OWLx5GUEEmvmohI9EUS/AVm9l+EpmZu2dWj4Zxt2FZaw4NvbuQzxw7nhNyBsS5HROQjkQT/ceHvJ7fYpuGc7bjnlTUA3Hb+UTGuRETkQJHcwDW7KwrpSQo27+P5pdv59pnjyM5IjXU5IiIHaLfj2cyyzGyOmb0cfj7ZzK6OfmnxqbnZ+clLqxjSL5nrzxgb63JERD4hkk8cHwFeBYaHn68FbopWQfHuxeXbWbq1lO+eO5G+ukNXRLqhQwa/me1PrUHu/hThoZzu3gg0dUFtcaemvomfvbyGvOz+/PvxObEuR0SkVW1d8b8f/l5lZpmE5+kJz81fFu3C4tFDb26kuKyWH1w4RatqiUi31VZfxP7kuoXQUM6xZvYOMBi4JNqFxZsdZbU8sHADnz56KCeO1vBNEem+2gr+wWZ2S/jxs8DfCf0xqAPOBpZHuba4cs+ra2hqdr53/qRYlyIi0qa2gj8BSOPjK//9NOHMQZZtLeWZxdu4/oyxjBio5hGR7q2t4C929x93WSVxyt25+++rGZTWmxtma/imiHR/bX24q08nI/Dmuj28/+Fe/t9Z4+mXkhTrckRE2tVW8J/VZVXEKXfn56+uIWdAKl84YWSsyxERicghg9/d93ZlIfHolRU7WLGtnJvOnkDvRM2+KSLxQWnVQU3Nzi9e+4BxQ9K4+LjsWJcjIhIxBX8HPbtkGxt2V3HrORNI0M1aIhJHFPwdUNfYxH2vr+Xo7HTOyxsa63JERA6Lgr8Dnly0lW2lNXzn3IlaPF1E4o6C/zBV1zfy63+s58TRAzl9/KBYlyMictgU/Ifp0X9uZk9lHd/V1b6IxCkF/2Eoq2nggYUbmD1xsNbRFZG4FbXgN7MRZjbfzFab2UozuzG8/dLw82Yzmx6t40fDH97aSFlNA7d+amKsSxER6bBoLhHVCNzq7ovNrB9QYGavAyuAzwEPRvHYna6kso45b3/IBccMIy87PdbliIh0WNSC392LgeLw4wozWw1ku/vrQNz1jz/8zofUNDRx89njY12KiMgRMXeP/kHMcoE3gTx3Lw9vWwB8x93zD7HPdcB1AFlZWdPmzZvXoWNXVlaSlpbWoX33q2pwvrOwmrxBCdwwNeWI3isWOqMN4l3Q2yDo5w/BbIPZs2cXuPsnutSjvhq4maUBTwM37Q/9SLj7Q8BDANOnT/dZs2Z16PgLFiygo/vu95t/rKOmcS0/+vwMpgyPv26ezmiDeBf0Ngj6+YPaoKWojuoxsyRCoT/X3Z+J5rGipaqukYff+ZAzjxoSl6EvInKwaI7qMWAOsNrd743WcaLtife2sK+6gRtmj4t1KSIinSKaXT2nAFcAhWa2NLztDiAZ+A2hRdv/ZmZL3f3cKNbRYbUNTTz01kZmjs1k2qgBsS5HRKRTRHNUz9scehWvZ6N13M70l4IidlfU8T+XTY11KSIinUZ37h5CQ1MzDyzYwHEjM5gxNjPW5YiIdBoF/yE8t2Qb20pr+NbscXF3z4GISFsU/K1oanZ+t2ADk4b158yjhsS6HBGRThX1cfzx6OUVxWzcU8Vvv3S8rvZFAqChoYGioiJqa2tjXUqHpKSkkJOTQ1JSUkSvV/AfxN353zfWM2ZwX62uJRIQRUVF9OvXj9zc3Li72HN3SkpKKCoqYvTo0RHto66eg/xj9S7W7Kjgm7PGaS1dkYCora0lMzMz7kIfQvOeZWZmHtb/rSj4D/LgmxvIzkjloqnDY12KiHSheAz9/Q63dgV/C8u2lrJo0z6uOnU0SQlqGhHpmZRuLcx5+0PSkhP5/PScWJciIhI1Cv6w4rIa/l5YzGUnjKBfSmSfjIuIxCMFf9ij/9xMsztfnZkb61JEJIB+//vfM3XqVKZOnUqvXr0+enzLLbd0+rE0nJPQ1MtPvLeZ8/KGMmJgn1iXIyIBdO2113Lttdeybds2Zs6cydKlS9vfqYN0xQ88vbiI8tpGrj41sjGwIiLRsmLFCo4++uioHiPwV/zNzc4f39nEsSMyOH6kpl4WCbq7XlzJqu0RLxYYkcnD+/PDf5sS0WsLCwvJy8vr1OMfLPBX/G+s2cWHe6q4+tTRcT2OV0R6htau+O+8885OPUbgr/j/8PZGhqencL6mZxARiPjKPFoKCwu5+eabP3q+Y8cOGhsbO/UYgb7iX7m9jHc37uXKmbm6YUtEYq65uZl169Zx1FFHfbRtyZIlTJ3auYtBBTrt5rz9IX16J/CFE0fGuhQREdavX09OTg7JyckfbVu6dKmCv7PsKq/lxWXbuXRaDumpumFLRGJvwoQJrFq16oBt69evZ/z48Z16nMD28f/pX5tpbHa+doqGcIpI9zVnzpxOf89AXvHXNjQx973NnD0pi9xBfWNdjohIlwpk8L+0vJh91Q187ZTcWJciItLlAhn8c9/bzJjBfZkxJjPWpYiIdLnABf/K7WUs2VLK5SeN0g1bIhJIgQv+J97bQnJiL/79+OxYlyIiEhOBCv7Kuv3eYCAAAAepSURBVEaeW7KNC48ZTkaf3rEuR0QkJgIV/M8t2UZVfROXn6wbtkQkuAIT/O7O3Pe2MHlYf44bkRHrckREYiYwwb9kaymri8u5/OSR+lBXRAItMME/990t9O2dwEVT9aGuiHQ/Wnqxk5VW1/PS8u1cMi2HtORAnLKIdNTLt8OOws59z6FHw/k/a/MlWnqxkz29eBt1jc1cftKoWJciItImLb3YCUIf6m7muJEZTB7eP9bliEh3186VebRp6cVOsGZvMxt3V/FlXe2LSBzoiiv+qAW/mY0ws/lmttrMVprZjeHtA83sdTNbF/4e1RXO529tID01iQuOGRbNw4iIdIqDr/jdnTvvvJOXXnqJRx99lN/+9rdHfIxoXvE3Are6+yTgZOAGM5sM3A78w93HA/8IP4+K3RV1FOxs4pJpOaQkJUTrMCIinaK1pRcXL17M6aefTn5+PvX19SQmJlJaWnpEx4la8Lt7sbsvDj+uAFYD2cBFwKPhlz0KfDZaNTyVv5Umhy+dpDt1RaT7a23pxTFjxvDggw9SV1fHW2+9xcKFC+nf/8g+rzR3P9Ja2z+IWS7wJpAHbHH3jBY/2+fun+juMbPrgOsAsrKyps2bN++wj/tWUQMrd9Vx/fFpHay8Z6isrCQtTW0Q5DYI+vlD222Qnp7OuHHjuriizrV+/XrKysoO2DZ79uwCd5/+iRe7e1S/gDSgAPhc+HnpQT/f1957TJs2zTtq/vz5Hd63p1AbqA2Cfv7ubbfBqlWruq6QKGntHIB8byVTozqqx8ySgKeBue7+THjzTjMbFv75MGBXNGsQEZEDRXNUjwFzgNXufm+LH70AXBl+fCXwfLRqEBGJlHdBt3e0HG7t0bziPwW4AjjTzJaGvz4N/Aw4x8zWAeeEn4uIxExKSgolJSVxGf7uTklJCSkpKRHvE7U7d939beBQ02CeFa3jiogcrpycHIqKiti9e3esS+mQlJQUcnJyIn59j5+yQUSkPUlJSYwePTrWZXSZHj9lg4iIHEjBLyISMAp+EZGA6ZI7d4+Ume0GNndw90HAnk4sJx6pDdQGQT9/CGYbjHL3wQdvjIvgPxJmlu+t3bIcIGoDtUHQzx/UBi2pq0dEJGAU/CIiAROE4H8o1gV0A2oDtUHQzx/UBh/p8X38IiJyoCBc8YuISAsKfhGRgOnRwW9m55nZB2a23syitrZvd2JmD5vZLjNb0WJbly5wH0tmNsLM5pvZajNbaWY3hrcHqQ1SzOx9M1sWboO7wtsD0wYAZpZgZkvM7KXw80Cdf1t6bPCbWQLwW+B8YDLwxfBi7z3dI8B5B23rsgXuu4FG4FZ3nwScDNwQ/u8epDaoA85092OBqcB5ZnYywWoDgBsJrfW9X9DO/5B6bPADJwLr3X2ju9cD8wgt9N6jufubwN6DNnfZAvex5u7F7r44/LiC0C9+NsFqA3f3yvDTpPCXE6A2MLMc4ALgDy02B+b829OTgz8b2NrieVF4WxBluXsxhIIRGBLjerqEmeUCxwHvEbA2CHdzLCW0tOnr7h60NvgV8B9Ac4ttQTr/NvXk4G9tERiNXQ0IM0sjtN7zTe5eHut6upq7N7n7VCAHONHM8mJdU1cxswuBXe5eEOtauqueHPxFwIgWz3OA7TGqJdYCtcC9mSURCv257v5MeHOg2mA/dy8FFhD63CcobXAK8Bkz20Soi/dMM3uc4Jx/u3py8C8CxpvZaDPrDXyB0ELvQRSYBe7NzIA5wGp3v7fFj4LUBoPNLCP8OBU4G1hDQNrA3b/n7jnunkvo9/4Nd/8yATn/SPToO3fDi7v/CkgAHnb3u2NcUtSZ2Z+BWYSmoN0J/BB4DngKGAlsAS5194M/AO4RzOxU4C2gkI/7d+8g1M8flDY4htCHlwmELu6ecvcfm1kmAWmD/cxsFvAdd78wiOd/KD06+EVE5JN6clePiIi0QsEvIhIwCn4RkYBR8IuIBIyCX0QkYBT8EhhmlmlmS8NfO8xsW/hxpZndH6Vj3mRmXwk/XmBmn1js28yONrNHonF8kdYkxroAka7i7iWEZqvEzH4EVLr7L6J1PDNLBK4Cjm+nrkIzyzGzke6+JVr1iOynK34JPDOb1WLO9h+Z2aNm9pqZbTKzz5nZPWZWaGavhKeDwMymmdlCMysws1f3TwVwkDOBxe7e2GLbpeG58tea2Wkttr9I6C5TkahT8It80lhCU/peBDwOzHf3o4Ea4IJw+P8GuMTdpwEPA63dFX4KcPBEYYnufiJwE6G7qvfLB05DpAuoq0fkk1529wYzKyQ07cEr4e2FQC4wEcgDXg9NDUQCUNzK+wzjwIVAAPZPGlcQfq/9dgHDO6F2kXYp+EU+qQ7A3ZvNrME/ntekmdDvjAEr3X1GO+9TA6S09t5AEwf+/qWEXy8SderqETl8HwCDzWwGhKaBNrMprbxuNTAuwvecAKxo91UinUDBL3KYwkt5XgL8t5ktA5YCM1t56cvA6RG+7Wzgb51ToUjbNDunSBSZ2bPAf7j7ujZekwwsBE49aASQSFQo+EWiyMwmElrr9c02XjMeyHb3BV1WmASagl9EJGDUxy8iEjAKfhGRgFHwi4gEjIJfRCRgFPwiIgHz/wHpEL1cG6rUfgAAAABJRU5ErkJggg==\n", |
|
|
"text/plain": [ |
|
|
"<Figure size 432x288 with 1 Axes>" |
|
|
] |
|
|
}, |
|
|
"metadata": { |
|
|
"needs_background": "light" |
|
|
}, |
|
|
"output_type": "display_data" |
|
|
} |
|
|
], |
|
|
"source": [ |
|
|
"import numpy as np\n", |
|
|
"import matplotlib.pyplot as plt\n", |
|
|
"import math\n", |
|
|
"\n", |
|
|
"x = range(0,48) # for 48 hours\n", |
|
|
"y = [T_infinite+(T_init-T_infinite)*math.exp(-x/tau) for x in range(0,48)]\n", |
|
|
"y2 = np.full((1, 48), T_infinite)[0]\n", |
|
|
"\n", |
|
|
"%matplotlib inline\n", |
|
|
"fig = plt.figure()\n", |
|
|
"plt.plot(x,y,label='$T_i$')\n", |
|
|
"plt.plot(x, y2,label='$T_\\infty$')\n", |
|
|
"plt.legend(loc=\"lower right\")\n", |
|
|
"plt.grid()\n", |
|
|
"plt.xlabel('Time (h)')\n", |
|
|
"plt.ylabel('Temperature (°C)')\n", |
|
|
"plt.show()" |
|
|
] |
|
|
}, |
|
|
{ |
|
|
"cell_type": "markdown", |
|
|
"metadata": {}, |
|
|
"source": [ |
|
|
"## Calculating indoor temperature over multiple days\n", |
|
|
"\n", |
|
|
"By defining schedules for the internal and solar heat gains, we can simulate the indoor temperature duration over several days.\n", |
|
|
"\n", |
|
|
"A discretized differential equation is defined in the function Ti_calc with a time-step delta_t set to 1 hour." |
|
|
] |
|
|
}, |
|
|
{ |
|
|
"cell_type": "code", |
|
|
"execution_count": 679, |
|
|
"metadata": {}, |
|
|
"outputs": [ |
|
|
{ |
|
|
"data": { |
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO2deZgcV3Xof6dn1eyaVZpFmpE02iVrsyxbthnZYBvbYLM4wOP5QR7BJDEBB+exJo9AQh6BhCQQCJg4gQQHQ7ANNjZeMBob29iyte/SSLMvmn3fp8/7o7rH43HPqHvU1V3VfX/f1990V3dVnXunqs6955x7jqgqBoPBYDDMxhNtAQwGg8HgTIyCMBgMBkNAjIIwGAwGQ0CMgjAYDAZDQIyCMBgMBkNAEqMtQDjJz8/X8vLyBe07NDREenp6eAVyEab9pv2m/fHZ/v3793eqakGg72JKQZSXl/Paa68taN/q6mqqqqrCK5CLMO037Tftr4q2GFFBROrn+s6YmAwGg8EQEKMgDAaDwRAQoyAMBoPBEBCjIAwGg8EQEKMgDAaDwRAQoyAMBoPBEBCjIAwGg8EQEKMgDAaDwQYmp7z8eF8DY5NT0RZlwdimIESkTET2ishJETkuIp/0bf+6iJwSkSMi8oiI5Myxf52IHBWRQyKysNVvBoMhLPQOj3Pt1/by0rnOaIviGn5zqp3PPXyUJ4+1RVuUBWPnDGISuFdV1wG7gLtFZD3wDLBRVTcDZ4DPzXOMPaq6RVV32CinwWC4CC+f76ahe5hnT7ZHWxTXsK+2G4AjTX1RlmTh2KYgVLVVVQ/43g8AJ4ESVX1aVSd9P3sZKLVLBoPBEB78D7vDjb1RlsQ97Ktzf59JJEqOikg58DzWzKF/xvbHgJ+o6o8C7FML9AAKfE9V75vj2HcBdwEUFRVtf/DBBxck4+DgIBkZGQvaNxaIt/Yf7pikIiuBrBQB4q/9s7lY+7/40gj1/V6SPfAvb00jwSMRlM5+wv3/H5lU/vjXw3gEEsTZfbZnz579c1lpbE/WJyIZwEPAPbOUwxewzFAPzLHrblVtEZFC4BkROaWqz8/+kU9x3AewY8cOXWjCrXhO1gXx1f7OwTE+/Ne/5g+uruDPb1wPxFf7AzFf+wdGJ2h86mlW5KdzvnOIpWu3s744K7IC2ky4///Pn+lA2ce7tpXys/1Nru0zW6OYRCQJSzk8oKoPz9j+IeBW4IM6xxRGVVt8f9uBR4CddspqiB9eq+sB4KCLp/6RZH99D16Fj1xTAcDhJtNvF2NfbTcJHuHDV5UDcMSlfWZnFJMA9wMnVfUbM7bfBHwGeKeqDs+xb7qIZPrfAzcAx+yS1RBfvOazDR9r7mNiyhtlaZzPq3XdJHqE27eUkJOW5GqbeqTYV9vNxpJsNhRnkZWa6FqlaucMYjdwJ3CdL1T1kIjcDPwzkIllNjokIt8FEJFiEXnCt28R8IKIHAb2AY+r6pM2ymqII16t7yHBI4xNejndNhBtcRzPvtpuNpRkk56SyObSHA67OConEoxOTHGosZed5YsRES4ry+Fwozv7zDYfhKq+AATyyjwRYJvfpHSz7/154DK7ZDPEL8Pjkxxv7uOWTUt59HALBxt72ViSHW2xHMvoxBSHG/v48O5yALaUZvPt6nMMj0+SlhxT9cbCxpGmPsanvOysyAPgstIc/uW5c4yMT7EoOSHK0oWGWUltiCsONfYy6VXetbWEvPRkDjW4c+ofKQ439jI+5eXy8lwALivLYcqrHG/pv8ie8cu+2i4AdixfDMDm0mymvMqJVvfNIoyCMMQVr9X1IALbli9mS1mOa23DkeJVn7/m8nL/w85KfGD8EHPzSm03a4oyWZyeDMCWMqvPDrnQzGQUhCGueLXOunmzFyVxWVkO5zoG6R+diLZYjuWV2m7WLskkJ8162BVkplCSs8j4IeZgcsrLgfoedlbkTm8rzEplaXaqKyOZjIIwxA3+m9dvLtlSloMqHHHhyC4SzO4vP5eVZZsZxBycaO1naHyKyyve2GebS93ZZ0ZBGOKGU20Db7h5L/ObS1w4sosE/ofdzjc97HJo6B6me2g8SpI5F39Kkp1vUqo51HUN0zvsrj4zCsIQN8y2p2enJbEiP52DxlEdEP+CwjfNIIxinZMDDT2U5CxiSXbqG7b7+8xtifuMgjDEDa/VWTfv0uxF09u2lOVwqLGXSOQkcxsHGnoozk5908NuU2k2IsZRHYiDDb1sXfbmCgb+Pjvksj4zCsIQF6gqr9Z1T88e/GxZlkPn4Bhdo0ZBzOZgQy9bly9+0/aMlERWF2aamdcsWvtGaO0bZduyN/dZVmoSlYUZHGjoiYJkC8coCENc0Ng9QvvAGDtmmUv8IYjn+0zKjZm094/S3DvC1rKA9bzYtnwxBxp68HqNYvVzoN5SmNsCKFWA7csXc6DeXX1mFIQhLvCP3LbPunnXLskiOdHD+V73loW0A39/zfWw27F8MQOjk9R0DEZSLEdzoKGHlEQP65cGztq6bdli+kcnOeeiPjMKwhAXHGzoIS05gdVFmW/YnpzoYWNxFjW9ZgYxkwMNvSQneNgwR4pqv6L1O7INloLYVJJNcmLgx6p/9rq/3j19ZhSEIS441NjL5tLsgEVbdpTnUtfnZXTCzCL8HKjvYUNJFimJgXMHLc9LIy892VUPOzsZm5zieHP/nDMugPK8NHLTk3nNRX1mFIQh5hmdmOJEaz9bAzgPwRoNT6qV/tsA45Nejjb3BXS2+hGRaT+EAY419zM+5WVbgAgmPyLCtmWWH8ItGAVhiHmOt/QzMaXTDunZTJtLXHTj2snJ1n7GJr3zKgiw+q22c4jOwbEISeZcDvp9NkH02fnOIdcsMjQKwhDz+G/euSJy8jNSWJIm04WE4p3XHdRzj4bh9WylbhoR28XBhl5KchZRmJU67++2u6zPjIIwxDwHGy9+81YuTmB/fY9ZMIfloF6SlfqGBYWB2FiSTVKCsN+YmTjQ0DOv/8HP5tJsEj3u6TOjIAwxz6GGXrbMYxsGWLXYQ8/wBOc6hiIklXM52NBz0dkDQGpSAhtLsl0zGraL1xfIBddnG0qyXePcNwrCENO0D8y/4MvP6hwrWmd/fXybmdoHRmnqGbmoLd3P9mWLOdzUx9hk/EaATS+QC6XPGntdUQ/dKAhDTOOvGBcoP85MlqQLuenJvBrncf3+h93F+svPjvLFjE9647rCnH+B3Lo5FsjNZvvyxYy5pM+MgogBPnDfy/xL9bloi+FIDjb2kugRNhTPX3faH4Lolqm/XRxs7CEp4eL95cc/ao5nM9NhX13zuRbIzcZvvnPDtWYUhMvpGRrnd+e7eOxwS7RFcSSHGnpZX5xFatLFi8XvKDdhm0eb+li7JLj+AqtaWlnuorhdUT05Zc0ENpcGp1ABlmYvoiRnkSvMmUZBuJwjvsVdJ9v66RsxpTNnMuVVjjT1zrn+YTb+TK9uGNnZgderHG3qC+lhB7BjeS6v1XfHZQRYTccgIxNT0/UeguXy8sXsq3V+nxkF4XKO+oq2qGLi+Gdxtt2qIBesPd1vJojXfqzrGmJgbDJkBXFFRS6dg+OuSkIXLvwFgDaF2mcr8nx95uyoOaMgXM6Rpj5KchaRnOCZLndosPA7qLeUBRddkpKYwOaS7LhdUe1/2G0OcTR8xYo8AF4+H3/X35GmXjJTEqnISw9pvyt8ZVxfqe2yQ6ywYRSEyznS1MeO8sVsKcvhZaMg3sCR5j6yUhMpz0sLep/LK3I52tTH8PikjZI5kyNNfaQkeqgszAhpv/K8NAozU3glDq+/o019bCzJxhMgCeR8VOSnU5CZwisOV6pGQbiY9v5R2vpH2Vyaw86KXI419zE4Fn8Ptrnw37wiwd+8u1bkMenVuPRDHG3uZUNxFokJoT0WRIRdK/J45XyX423q4WR80svJ1gE2l4VmXgKrz66oyOWVWmf3mVEQLuZos98kkM0VK3KZ8mpchxvOZHzSy+m2ATaVhOpwXUyiR3j5vLOn/uFmyqsca+4P2bzk54oVubQPjFHXNRxmyZzL6bYBxqe8bC5ZaJ/lcaF/jHoH95lREC7mcFMfHoENxVls9z3YnG7TjBRnLlg3b6jOw/SURDaXZvO7c/HVj61DysjEVMgOaj9XVFh+iFfiSLEe9gWILLTPdrnAD2EUhIs52tRLZWEmacmJpCUnsrEk2/E2zUjhn12FOoMAy8x0pKmPoTgy19X2WakyFjqDWFmQTn5GfPkhjjb1sTgtidLF8yc1nItVhRnkpSc7+p41CsKlqCpHm/veMEK+YkUuh5t6TWU0LIdrVmoiy3KDd1D7iUc/RG2/l/TkBFbkhxaN42faph5HfojDTb1sLs0Jycc1ExFhZ0Wuo5WqbQpCRMpEZK+InBSR4yLySd/2r4vIKRE5IiKPiEjAIYuI3CQip0WkRkQ+a5ecbqWlb5TOwfE3TG+vqMhlYkpNlS+s6nChOqj97CiPPz9EXZ93QdE4M7liRS4tfaM0do+EUTJnMjI+xdn2wQWbl/xcUZFLc+8ITT3O9EPYOYOYBO5V1XXALuBuEVkPPANsVNXNwBngc7N3FJEE4NvA24H1wAd8+xp8HJ22f76uX3eU5yJC3K+HmHZQL/DmTUtO5LKyHH4XJwpifNJLQ7+Xy4JccT4Xfj/Eyw62qYeLE619THl1wSY5P/41JE41M9mmIFS1VVUP+N4PACeBElV9WlX9xt2XgdIAu+8EalT1vKqOAw8Ct9klqxs50tRHokdYuyRzeltWahLrl2bF1cg3ENMO6gX4H/zsWpEbN36IMxcGmNSF+WtmUlmYweK0JMc+7MLJ4cbXIwgvhTVFmeSkJTnWUZ0YiZOISDmwFXhl1lf/G/hJgF1KgMYZn5uAK+Y49l3AXQBFRUVUV1cvSMbBwcEF7xsNnjs6QkmG8PKLv33D9rLkMX5dO8lTz+4lJSF4c4Hb2j8f1Y1WTqrhplNUd58Jap/Z7V80MMWUV/m3R6vZVBCR2yRq+PtrpPkU1T3B9ddcrMj08tzJZqqr3WXmDPX6f+bIKDkpwskDL3PyEs9dkeFl7/FmqvOd12e2X/kikgE8BNyjqv0ztn8Bywz1QKDdAmwL6PlS1fuA+wB27NihVVVVC5Kzurqahe4baVSVT1Q/zS2bi6mq2vTGL5e28+S/v0raso1cU1kQ9DHd1P6L8dTDR8lKbeGOt+8J2gcxu/07xyf5p4NPM5xZSlXVWpskdQZPPnSE9KTGkPprLuqSavnLx06wcvNOyhYQIBAtQr3+/2p/NTtWpFNVdfkln7s2qZYvObTPbI1iEpEkLOXwgKo+PGP7h4BbgQ9q4JCHJqBsxudSwOSz9tHYPUL/6GRAk8DOilySEoQXajqjIJkzOOaL7rqUh11aciKXlebExXqIYy19lGd5Llk5AFxdmQ8Q09ff8Pgk5zuHgq6ZcTGu8fXZiw7sMzujmAS4Hzipqt+Ysf0m4DPAO1V1Ltf9q0CliFSISDLwfuBRu2R1G8dbLPvnhuI3V7BKS05k67LFvFQT+w+2QPgd1Bsv0Z4OVrjr0RhPXzIx5eVM2yDLsoKr/3AxVhZksCQrlRfOOu9hFy5OtQ2gGvj+WwgrCzIoykpxpFK1cwaxG7gTuE5EDvleNwP/DGQCz/i2fRdARIpF5AkAnxP748BTWM7tn6rqcRtldRUnWvtJ8AhrZjioZ3L1qnyOtfTRMzQeYcmiTzgc1H6uWpXHlFd5OYZnETXtg4xPeVmWGZ5HgYiwe1U+L57rxOuNzfUQ/lKh68OkIPx99tK5Lsf1mZ1RTC+oqqjqZlXd4ns9oaqrVLVsxrY/9P2+RVVvnrH/E6q6WlVXqupX7JLTjZxo6WdlQfqcVb92r8pHlbgJ05zJdH7+MCiI7csXsygpwZEju3BxwvewW54VvkfBNZX59A5PcKLV+TWXF8KJln6yFyVRkrOwFdSBuHpVPt1D447rM7OS2oWcaO1n/TwF0i8rzSYjJTGmH2xzcaxl4SuoZ5OSmMAVK3J5/mxHGCRzJida+0lN8rAk/dL9D36uWmXF9v82Rs1MJ1r6WL80Kyw+Gz9Xr3KmH8IoCJfRPTROa9/ovA6yxAQPu1bkOe5iiwQnWvpZXxy+m/eaygLOdww5dqXrpXKipZ81S7LwhPFhV5iZypqizJi8/ianvJxqGwib/8FPYVYqq4syHDeoMwrCZZwI0v65e1Ue9V3DNHbH5oMtEFNe5VRbP+vmmV2Fij/CJBadrqrKcd9oONxcXZnPvrrumMsLdr5ziLFJb9j8DzPZvSqffbXO6jOjIFzGiVbLxn6xh6B/yvrSudh7sM1FXdcQoxPesD7wKgutCJPfOmxkFw6ae61w6XCPhsG6/sYnvbxW57zFX5eCf4AWrhDXmVy9Kp+xSa+jaroYBeEyjrf0szQ7ldz05Hl/t6owg8LMFF6Io3DXkz4HXzhnECLCNZUFvFjTyZTDIkwulWBnowshVtfjHG/pIznRw4qChWW9nY8rVuSR6HFWnxkF4TJOtPQHNeITEa5elc+LNbEbbjibk639JHqEyqLQaipfDH9UzjFfjYlY4URrPyK8IZ9XuEhPsdbjvFATWw7+E639rF2SSVKIZVmDISMlka3LchzluzEKwkWMTkxxrmMwaBPKNaut0LljLbH1YJuLEy39rCrMICUxPIu+/OxeFZurg4+39LMiP520ZHsy7ly9Kp/jLf10x8h6HMtnM38E4aWye1U+R5qds4bJKAgXcbptAK8GbxK4trIAEag+HVujuLk42ToQVvOSn/yMFDYUZ/H8mdjqRyviK/y2dD/XVFrrcX4bI2HCrX2j9A5P2OKz8fOW1QVWnzlkMGIUhIs4HqKDLC8jhc2lOew93W6nWI6ge2ictv5R1i0Nv7kErKicAw09MZN2o294gubeEVtHw5tLc8hNT46ZAUq4V1AHYrrPTjnjnjUKwkWcaO0jMyUxpBq4VasLONTYGzPT/Lmww0E9k2srC5iYip20G/4Vu3Y+7BI8wltWF/DcmY6YcPCfaPH7bOzvs+ozHY7wHRoF4SJOtPSzLsRFYHvWFsbUNH8u7FYQO8oXk5acwG9iZDbmT/ho5wwCoGpNAd1D4xz2VUB0M8db+qjISyc9xd4qCf4+O+KAoIh5FYSIXCki3/bVj+4QkQYReUJE7hYR+4yXhjdhLQIbCPmG3lySTV4MTfPn4kRrP4WZKeRnpNhy/JTEBK6pzGfvqXYCZ6h3F/7+Ksi0p7/8vGV1AR7BMSaTS+FEqzVAs5trK60+2+uAPptTQYjIr4A/wMqoehOwFKs+9J8DqcAvROSdkRDSYC0CGx6fCtlB5vEI18bQNH8u7HJQz+S6tYW09o1ysnXA1vNEAn9KErvJSUtm27LF7HX5AKVvZIKmHnt9Nn4Wpyezddliqh0wW51vBnGnqn5EVR/1ZVqdVNVBVT2gqn+vqlXASxGSM+65FBPK9JQ1Bqb5gRif9FLTPmD7A2/PmkIA1zv9rf4atF2h+tmztpCjzX20D4xG5Hx2cOaCNSiwKwhiNnvWFHC4qY+OgbGInG8u5lMQOSKye/ZGEblGRFYCqKozYrHigNNtAyR4hFWFoS8C809ZY9XMVNM+yMSU2v7AK8xKZVNJNs+evGDreezmfOcgk161ZYFcIKrWWKVvn3Px9Xe6zVIQa2x0UM+kyjcYiXZo9XwK4h+BQHPpEd93hghyum2A8ry0OWtAzMfi9GS2lOU4YspqB/7Z1foIjO6uW1vIQZdHhb3+sIuMgli/NIvCzBRXD1BOtw2QmZJIcXZqRM63odjqs2jPVudTEOWqemT2RlV9DSi3TSJDQE5fGLikG7pqTSFHmvvoHIzulNUOTrb2k5LooTwv/PlxZnOdLyrMzcr2dNsAiR5hRX54U5LMhYiwZ00hz5/pYGLKG5FzhpvTbQOsXpIZ1hoQ8yEiVK0p4PkzHUxGsc/mUxDzqcrwlVIyXJTh8UkauodZU7Tw6a3/wfYbB0RGhBt/fpxEG/LjzGZTSTb5GSmu7sczFwZYUZBOcmLkotz3rC1gYGzSldldVa008pGacfnZs6aQ/tFJDjREz3c43xXyqoh8dPZGEfkIsN8+kQyzOXthEFVYs2ThI74NxVkUZ6fyzAl3288Dcbrt0mZXoeDxCHvWWFFhbh0Nn2obYHVRZB92V1cWkJzgcaX/5kL/GP2jkxHz2fi5ZrXVZ8+caIvoeWcyn4K4B/h9EakWkb/3vZ7DCn39ZGTEM4BlXoJLc5CJCDdsWMJvz3YwMu6cgiSXSufgGF1D4xF94F2/rpCB0Un2Oyhvf7AMjk3S1DMS8YddRkoiV63K4+kTF1y3juRUm+XjirRSdUKfzakgVPWCql4FfAmo872+pKpXqmr0VFoccrptgNQkzyXXWb5hfRGjE96YqrF85kJkHa5gjYaTEoRfu3A2diYMg42FcsP6JTR0D08PeNyC36kfaaUKVp/Vdw1z5sJgxM8NQaTaUNW9qvot3+s3kRDK8EbOXBigsjCTBM+lOcgur8glKzWRp4+778E2F2d9N86aCI7uMlIS2b0qn6dOtLluNDwdwRTh0TDAW9cXIgJPHXPX9Xe6bYCirBRy0uYv0mUH/j57+nh0xuTzraS+Q0R+LiKPiMj7IimU4Y2Ey2aclODh+nVF/ObUhahGRoST0xcGyF6UZHvKiNnctGEJjd0jrltVfbptgLTkhJASPoaLwsxUti1bzNNRtKkvBCuCMPIzLrD6bGtZDk9HabY63wziM8C7gfcAn46MOIbZdA+N0zEwFrbp7Q3ri+gZnuA1F9rPA3GmbYA1RZELP/Tz1vVFeASejNLIbqGcbhugsigTzyXORhfKDeuLON7ST2P3cFTOHyqTU17Otg+yJsxVCkPhxg1LONrcR3PvSMTPPZ+C+BHwH77Xf0dGHMNs/CaB1WFSENeuLiA50RMTZiZVtcxvUbh58zNS2FGey1PH3KMgVJXTFwZYGwXzkp8bNiwBcE00XV3XMOOT3qjNIGBGn0VhMDKfk/ofgY8Bf6iqX42cSIaZ+J2K4ZpBpKckcvWqfJ52of18Nv7ww0jHp/u5acMSTl8YoLZzKCrnD5XOwXG6h8aj1l8AFfnprC7KcI2ZKdz330KoyE+nsjAjKmam+XwQoqpDqjqn+1wiPa+PQ061WTb2wjDa2G9YX0RTzwin2txlP5+N/+atLIzOzXvjRmtk95RLzEyRTrExFzduWMK+2m5XpCs51TaAR1hQDrRwcsOGIl6p7aZ3OLJ9Np+Jaa+I/ImILJu5UUSSReQ6Efkh8CF7xTOc8aXYCKcuvn6dZT//lYvMI4HwK4jVUbIPl+QsYnNpNk+6pB9PRyEkOBA3rF+CV3HFornTbf2U56UvKAdaOLlh/RKmvMqzJyO7gn8+BXETMAX8WERaROSEiJwHzgIfAP5BVX8QARnjFlWddsKGk4LMFK6oyOPxIy2uNjOduTBAfkYyeTYVCQqGGzcs4VBjL619kXcghsrptn7y0pNtK6oULBtLsijJWeQKxRrJVfrzsakkm+LsVJ442hrR887ngxhV1e+o6m5gOXA9sE1Vl6vqR1X1UMSkjFNa+kYZGLPHxn7L5qWc6xhytZnp9IXBiK9unc1NPjOTG5z+TnnYiQg3b1rC82c76BueiLY4czIyPkV997Aj+szjEW7etDTifRZUti5VnVDVVlUNOmuUiJSJyF4ROSkix0Xkk77td/g+e0Vkxzz714nIURE5JCKvBXveWOK0b4m/HRfo2zcuwSPw+JHIjkjChder1FyIfE6h2awsyKCyMCPiI7tQ8XqVMw5QqH5u3VzMxJQ62n9ztn3AyoHmlD67zNdnEXTw25nOcRK4V1XXAbuAu0VkPXAMa33F80EcY4+qblHVORVJLHO6zYoPsOOmzstI4aqV+fzSpWam5t4RhsanHPHAu2XzUvbVddPW59yKaU09I4xMTDliNAywuTSbZblpPHakJdqizIl/lX6lA64xgMtKsynLXcQvIzios01B+GYcB3zvB4CTQImqnlTV03adN5Y4224t8c9elGTL8W/dvJS6rmGOt/Tbcnw7eT2nUHSjSwDeeVkxqvBLJz/s2qPr0J+NiHDr5qW8dK6LLofWKKnpGCTRIyzPu7QcaOFCRLhlUzEv1nTSE6EIsMRgfiQiy4FKVf21iCwCEn0P/aAQkXJgK/BKCLIp8LSIKPA9Vb1vjmPfBdwFUFRURHV1dQineJ3BwcEF72sXB8+NkJeEbXJljCsJAt/55SvcXDLuuPbPx5PnrRuk7cxhqmsvPcLrUv//y7M8PPDCaVZNNVyyLHbwVO38/RWN63/J+BRTXuWbDz/PnmX2DIKCJVD7Xz4xSmEavPjbYIwdkWHphNVn//Twc1SV2d9nF1UQvpoQdwG5wEqgFPgultP6oohIBvAQcI+qhjJU3a2qLSJSCDwjIqdU9U3/KZ/iuA9gx44dWlVVFcIpXqe6upqF7msHqkrH3qd5z7YSqqo22nae/27ax9HOQe5YneSo9l+MX1w4xNLsLm55256wHO9S//8f9Jzjb544RcWmy1kegcp2ofJ4x2HyMzrm7K9oXP+qyg/OPseZ0VS+VLUroueeTaD2f+m1ai4rz6Sqant0hAqAqvKDM89xdiyVv4xAnwVjYrob2A30A6jqWaAwmIOLSBKWcnhAVR8ORTBVbfH9bQceAXaGsr/baesfZXBs0vYFOrdsXkpj9wh1/e5K3mel2HCGbRjgls3FADx22JlmppqOQVYVOktxWWamYl6u7aK931n+m7HJKeq7hqK+QG42ftPc78510TFgv2kuGAUxpqrTBi8RScQy/8yLb5X1/cBJVf1GKEKJSLqIZPrfAzdgObfjhpp2y0G20uYL9Mb1S0hKEF5pnbT1POFkyqvUtA+y2kE3b0nOIi4vX8yjDlQQqlZ/Oe1hB/COzUtRxXFRYHWdw3g1+iuoA3Hr5mK8Ck8es7/PglEQz4nI54FFIvI2rMR9jwWx327gTuA6X6jqIRG5WUTeJSJNwJXA4yLyFICIFIvIE759i4AXROQwsA94XFWfDLFtrsYfQWH3BZqdlsRbVhfycqtl23QDTT3DjE16HXs2oVAAACAASURBVBHBNJN3XlbMmQuD0xXInELHwBgDo5OsKnDew66yKJM1RZk85rBwa/8AzYkKYs2STCoLM/jFIfsHI8EoiM8AHcBRrOR9TwB/frGdVPUFVRVV3ewLVd2iqk+o6iOqWqqqKapapKo3+n7foqo3+96fV9XLfK8NqvqVhTfRndR0DFp1DiKw6vU920roHVNeOtdp+7nCwbkO/+zKWSaTt29aSoJHHGdmev1h5yyF6ue2rcXsr++hvss5SQ/Ptg8gYq1zcSLv2lbCaxHos3kVhIh4gKOq+n1VvUNV3+t7746hpovxmwQikQ/xunWFpCXCwweabT9XODjXbt0UTrt58zNSuGplHo8ebsHroNlYTYc/nt9Z/eXnXVtLEHHW9VfTPkjp4kVRz8E0F7dvsfrskYP29tm8CkJVvcDh2Qn7DPZzrn0wYiaBlMQEdi5N5MljbQyNOd8XUdM+SH5GclRKQF6Md28robF7hFfruqMtyjRnLwySmZIY1ozA4WRp9iJ2r8zn4YNNjlGsNRG8/xZCcc4irlqZx8MHmm1d6BqMiWkpcFxEnhWRR/0v2yQy0DM0TtfQeETtn7uLExmZmHJFArVzHYOscOjNe+OGJWSkJPKz/U3RFmWamvZBVkZoNrpQ3rPdUqxOqHQ45VXOdw45KkouEO/eWkpD97CtfRaMgvgScCvwZeDvZ7wMNuE3CURSQazK8bAsN42HDzrnwTYX5zoGHWde8pOWnMgtm5by+NFWx8zGrBBXZ/aXnxs3LCE9OYGHHKBYG7utKnJOnkGAlSgyLTmBhw/Y12cXVRCq+lygl20SGaISQSEivGtrCS+d63J06uruoXF6hidYWeAsB/VM3rujlOFxZ8zG+kYm6BgYc7yCSEtO5O0+xTo6MRVVWSIVYn6ppKckctPGJfzyiH19dlEFISIDItLve42KyJSIOCuOL8aoaR8kNclDSc6iiJ733dtKUIWfH3RWFM5M3HDz7li+mPK8NEeYmaYHGw4fDYN1/Q2OTUaltOZMojGDXyjv2VbKwOgkv7ap+FIwM4hMVc3yvVKB9wD/bIs0BsC6qVfkZ+DxRNZmvDwvnR3LF/Oz/Y2OzfDqD3F18gNPRHjv9lJ+d76Lxu7hqMpyzsHx/LPZVZFHSc6iqJuZatoHKcy0L0lmONm1Io+l2am29VnI2VxV9efAdTbIYvARzVWvv7ejjHMdQ45wFgbiXPsgKYmRn12Fyru2lSICD9loHw6Gmo5BkhM9lOU6IyPpfHg8lpnzt2c7omrmdOqq80AkeITbt5bQ2DPC+GT40+UEY2J694zXe0XkqwSRasOwMIbGJmnuHYnaBXrrZUvJSEnkx/ucmZXUH8EU6dlVqJTkWKGbDx2IbuimNRtNJ8Hh/eXnfZeX4VX46avRUaxOTksyF5+8vpJn/vRakhPDX70hmCO+Y8brRmAAuC3skhgAON9hLQKL1gWalpzI7VuLefxIqyPLQdZ0DDraQT2TO3aU0tg9wkvnuqImg9sedmW5aVxTmc9PXm2ISuqXC/1jDI5NUumiPktNSrAthDkYBfGvqvr7vtdHfWkvKm2RxkBNh1VmI5o39fsvX8bYpJefH3LOylaA0YkpmnqiN7sKlRs3LGFxWhL/ta8+KucfnZiisWfYNf3l53/sXEZL3yjPnWmP+LndEAQRSYJREN8KcpshDNS0D5LgEcqjWFNgY0k2m0uz+fG+Bkc5q2s7h1B1XoqNuUhNSuCOHWU8ffxCVNJZn++w+sttCuKt64vIz0jhv15pjPi5a9qjP0BzEnMqCBG5UkTuBQpE5FMzXn8JODNBSQxQ0z7I8tw0W+yJofD+y5dxqm2AQ429UZVjJtNJ+lyiIAA+sHMZk17lJ69G4WHnonDNmSQlePi9HaX85tSFiDurazuHyExJjEiSTDcw31MoGcjAqjqXOePVD7zXftHik9rOIUekkXjnlmLSkhMc5ayuaR9EBCry3eGDAEvWayrz+fG+yNvUazuGECGqs9GF8v7Ll0XFWX2+c4iKgnRHpyWJJHMqCN+K6S8Bu1T1SzNe3/BVlTOEmSmvUtc1zAoHOGEzUhK5bUsxjx1upW/EGc7qcx1DlOQsYlGyuyawH7zCsqlXn46sTb22c5DibOdmJJ2PZXnRcVbXdg6xwkUDELsJxo4xLCJfF5EnROQ3/pftksUhLb1WLLNTRsgfvGI5IxNT/PdrkTePBOKcyyJy/Fy/rojCzBQeeCWyszFrNuqMa2kh+J3VvzkVGcU6OjFFc+8IFfnuu8bsIhgF8QBwCqjAStxXB7xqo0xxS22nFeLqFAWxsSSbneW5/Mfv6qNebc7rVc53OjdJ33wkJXh4/+Vl7D3dHrGV1apWRlKnXEsL4W3riyjOTuXfX6yNyPnqu4ZRhQoXK9VwE4yCyFPV+4EJn9npfwO7bJYrLqnzVYdy0hT3w7vLaegeZm+ERnFz0dw7wuiE15UKAuD9O5fhEeE/X45MyGvX0DgDo5OuVhCJCR7uvLKcl851RaSMa22n5dR30v0XbYJREH4DdKuI3CIiW4FSG2WKW853DJGenECBgwq73LC+iKXZqfzgpbqoynHeN7tyq8mkOGcRN21cwo/3NUQkDXidw2ajC+X9l5eRmuThBy/W2X4u/zVW7vI+CyfBKIi/FpFs4F7gz4B/Bf7UVqnilFoHRlBYo7jlvFDTydkLA1GTw//Ac/Po7iNXVzAwOhkRn860QnW5PX1xejLv2lrCIweb6Rkat/VctR1DFGamkJGSaOt53MTFalInAJWq2qeqx1R1j6puV1VTUc4GajuHHOkge//ly0hJ9ER1FlHb6bzZVahsW7aYrcty+PeX6mz36dR2DpGUIJQsdnZSw2D48FUVjE16+fGr9jr5a13us7GDi9WkngLeGSFZ4pqxySmaeoapyHNe1s3c9GRu31LCwweao5afqa5riOV5zppdLYSPXF1Bfdcwz9qUv99PbccQy3LTXJOkbz7WLMlk96o8/vN39UxMhT9jqR+3R33ZQTAmppdE5J9F5BoR2eZ/2S5ZnNHYPYzXwREUH95dzsjEFD96JTp5hepiZHR304YllOQs4v4X7I3McepsdKH8/lUVtPaN8iubqvQNTShdQ+MxcY2Fk2AUxFXABt5Yk/rv7BQqHvFncXXqTb1uaRZVawr49xdrI14ScmLKS2PPCOX5zptdhUpigocPXbWcV2q7OdbcZ8s5vF6ltiu2RsPXrS1kRUE633vunC35wS4MWTMTp95/0SKYinJ7ArxMwaAwM70GwsFpEf7wLSvpHBznvyNc8aupZ4Qpr7oyZUQg3nf5MtKTE7jv+fO2HL+lz1kLLsOBxyP84bUrOd7Sz2/Pdob9+G3DltKJpT4LB8EUDCoSkftF5Fe+z+tF5CP2ixZf1HUNkZeeTHaac8scXlGRy5ayHL7//HkmbbQFz8Yfnx4rN2/2oiQ+uGs5vzzSMh2dFU6ctuAyXNy2tZglWal8p7om7MduG/LiEVjmgsp7kSQYE9MPgKeAYt/nM8A9dgkUr5zvcL6NXUT4o6qVNHQP84RNtuBA1HZaq4+d3j+h8AdXV5CY4OF7z58L+7FrYyAkOBApiQn8wTUVvHy+mwMN4S2J2zbkpcwBWZSdRjC9ka+qPwW8AKo6CUTWCB0HuCXE7m3rilhZkM53q+2xBQeirnOIzNREctOTI3K+SFCYlcrv7SjlZ/ubwp7S2okLLsPFB3YuI3tREt+tDq9ivTCsrrj/Ik0wCmJIRPLw1aEWkV2APd61OGVwbJL2gTHHRjDNxOMRPnbtSk609vO8DbbgQNR1WcrT7SGus/nYtSvxKnz/+fBGNDlxwWW4SE9J5ENXlfP0iQthW7ipqrQNxZbPJlwEoyA+BTwKrBSRF4H/AP7EVqniDLetEr59awnF2an806/PRGQWUds5FDMO6pmU5aZx25Zi/mtfPV2DY2E7bqyFuM7mw1eVsygpge+EaRZxoX+MsSn33H+RJJgopgPAW7DCXT8GbFDVIxfbT0TKRGSviJwUkeMi8knf9jt8n70ismOe/W8SkdMiUiMinw2+Se7DbTlgkhM93H3dKg409No+ixibnKKld8Q1fRMqf1y1krFJL/8Wpoyl0wsuY7S/wFq4eeeVy/nFoebpGtKXwvnpIIjYVaoLJZgoplTgE8BfYaX7vtu37WJMAveq6jqs7K93i8h64BjwbuD5ec6ZAHwbeDuwHviAb9+YpNa3BsJNo+Q7tpdRkrOIf3jG3lnE9ALCGFgDEYhVhZncvHEpP3ixju4w5Bry91esj4Y/du0KUpMS+Oazl167bDrqywUm3kgTjInpP7AWyn0L+GesB/Z/XmwnVW31zT5Q1QHgJFCiqidV9fRFdt8J1KjqeVUdBx4EbgtCVldS12VVSnNT5a/kRA8fv24Vhxp7qT7dYdt5/BFMblKeoXLPWysZnpjiu89dusnk9QWXsdtfAHkZKXzoqnIeO9LCmUv0RdR2DJHkgaVZwYx744tg0hauUdXLZnzeKyKHQzmJiJQDW4FXgtylBJiZ8rIJuGKOY98F3AVQVFREdXV1KKJNMzg4uOB9L5XD50fISSRq54eFtb/AqxQsEr788H64MtUWp+iva63cT82nDtF33j6nazT//wBXLk3k3184zzpPK4tTFx5q+etaaxbSdOogPeeC769ot38hbPAoKR74/H+9yMe3Lvzh/trpUfJTleeffy6M0sUGwSiIgyKyS1VfBhCRK4AXgz2BiGQADwH3qGqwVT8CXdkB7Riqeh9wH8COHTu0qqoqWNHeQHV1NQvd91L5RPVTvHNdMVVVm6Jyflh4+7szG/n0Q0fwLlnP9euKwi7X0z1HWZzWyq037An7sWcSzf8/QMWmIa7/++c4OFbIl2/auODjPNV9hNz0C9zyttD6K9rtXyhnOM03f1ND4eptrC/OWtAxvnLgOZZkjLiy/XYTzFDlCqyEfXUiUgf8DniLiBwVkXmd1SKShKUcHlDVh0OQqwkom/G5FGgJYX/X0Ds8Tv/opGtNKO/aVsLyvDS+/tRpW1JY13UOxayDeibL89K5Y0cZP97XcEllSeu7huNqNfBHrl5BZmoi33jmzIL293qV+u5hCtNiLyQ4HASjIG7Cqkf9Ft+rArgZuBV4x1w7iWVvuB84qarfCFGuV4FKEakQkWTg/VihtjFHfZf1MHDrTZ2U4OH/3LiGU20DPHKwOezHr+sccnR+qnDyietXISKX5Hit7xpmuQNTxttFdloSH7t2Bb8+eYF9td0h739hYJTxSS+FaWYFdSCCCXOtB/qBbCDP/1LVet93c7EbuBO4TkQO+V43i8i7RKQJuBJ4XESeAhCRYhF5wnfOSeDjWCk+TgI/VdXjC2+mc6n3jRaXu/gheMumpVxWms3fP306rJleR8anaOkbjYsZBMDS7EXcuWs5Dx1o4nRb6I7X8UkvrX0jLHfpYGOhfOTqFSzJSuVvnjgZckSdf4BWZGYQAQkmzPWvgCPANwkh3beqvqCqoqqbVXWL7/WEqj6iqqWqmqKqRap6o+/3Lap684z9n1DV1aq6UlW/suAWOpyGLivqxK0zCLByNH3u5nW09o2GLZ4foL7bXetDwsHH96wiIyWRrzxxMuR9m3qsENdlLh5sLIRFyQl86obVHGrs5fGjrSHt2+BTEGYGEZhgeuX3gJWqWmXSfYef+q5hCjNTWJTsnhDXQOxakcdb1xXyL3vPhSWeH15fYR4vJiawajB/4vpKnj/TQfXp9pD2fX026t7BxkJ5z7ZS1i7J5G+fPMXYZPCz2PruIRI8Qm6qmUEEIhgFcQzIsVuQeKW+O3Zsxp+5aS1D45NhWbwEUOcb3S2P0UVyc/G/riynPC+Nrzx+MqS06v7RcLyZmAASPMLnb15HY/cI//m74Kse1ncNU5KziMQYKM1qB8EoiP+HFer6lIg86n/ZLVi80NA1zLLc2BghVxZl8r7Ll/Gjl+vDkkitoXuY3PRkslKdWyPDDpITPXz27es42z7Ig682XnwHH/VdwyxKis0srsFw7eoCrqnM51u/qaEnyFlsQwwN0OwgGAXxQ+Bvga/yug/i7+0UKl4YnZiirX80pi7QP7thNWnJCfzlY8cvOQVHY/cwZXE4Gga4cUMROyty+YdnztA/OhHUPg3dQyzLTYvJLK7B8ue3rGdwbJKvPXWxZA0W8RYWHCrBKIhOVf2mqu5V1ef8L9sliwMaY9BmnJeRwp/duIYXa7ouucB8PN+8IsL/vXU9PcPjfOPp4GL867uGWRZD19JCWLMkkw9fVc6DrzZwpKl33t/2DU/QNzIRU/dfuAlGQewXkf8nIleKyDb/y3bJ4gB/iF2sjZL/x85lrFuaxV//8gTD45MLOsbklJfm3hGW5S4Ks3TuYWNJNnfuWs5//K6OY83zl2DxetUyl8TYtbQQ7nlrJfkZKfzFz4/hnWfxpj9Kzs0h5nYTjILYipWN9W8IIczVcHGmo05i7KZOTPDwV7dtoKVvlH9ZYM7+1r5RprzK8hjxzyyUT92whtz0FL5wkYdd+8AYY5NeMxoGMlOT+PzNaznc1MdPXpvbhzMdBGH6bE6CWSi3J8DLhLmGgYauITJSYquUpp8d5bm8e2sJ33vu/IJy9jd0x+bsKlSyFyXxhVvWcrixd16Hdb1/PY0ZDQNw+5YSdpbn8rUnT83psI6FNUh2E8xCuSIRuV9EfuX7vF5EPmK/aLFPffdwTDsVP3fzOhYlJ/C5h4/MO/oNxHQKEjO64/YtJVxRkcvfPnmKzjkqz8XqbHShiAhfvn0DA6OT/NXjJwL+pr5rmILMFNKSg8lZGp8EY2L6AVbKi2Lf5zPAPXYJFE80xHjenILMFP78lnW8WtfDA/saQtq3oXuYpARhicnRj4jw17dvZHh8ki89Fvhh19A1TIJHKFkcvz6b2axdksUfVa3k4QPNARcd1hufzUWZU0GIiF+t5qvqTwEvTOdJCl/CnThlyqs09sR+1Ml7t5dyTWU+f/urU7T2jQS9X2P3MGWL00gwC5gAa43Jn1xXyWOHW3jq+Jujw+q7hynOSSUpwaSMmMnHr1vFyoJ0vvDIMQbH3hgw0WCivi7KfFfTPt/fIRHJw1ePQUR2AfOHVBguSmvfCBNTse+EFRH+5l2bmPIqf/7IsaDXRjTE8RqIufijqpWsX5rFFx45Ru/wG+3qDV1DMX8tLYSUxAS+9t7NtPSN8PUnT01v969Bcmua/Ugxn4LwD90+hZVqe6WIvIhVgvRP7BYs1mmIowiKstw07r1hNc+eaufnh4JLCV7fNWSch7NISvDw9Ts20zs8zpdnmZrqu81oeC62L8/lQ1eW88Pf1U+nBI/FNUh2MJ+CKBCRTwFVwCPA14BfAd8H3mq/aLGN36kYLw/B399dweXli/m/Pz9Oc+/8pqa+4Qn6Ryfjpm9CYUNxNn9ctZKHDzbz7MkLAPSNTNA7PGHs6fPwf25cw7LcNO7970MMjE64vg5LpJhPQSQAGUAmkI5VnjQBSPNtM1wC9V2WE7Y4Jz6cigke4Ru/twWvKvf+9NC8UU0mxHV+Pn5dJWuXZPKZh47QMTAWV7PRhZKeksg/vO8ymntG+PJjJ2KiDkskmC++q1VVvxwxSeKMhu4hSuPMCVuWm8YX37mBT//sCPe/UMtHr10R8HcNZvo/L8mJHr75ga2841sv8H9+dpj3bCsFiJmkj3axfXkud+9Zxbd+U8PKgnQyUxJZnBZfiSBDJRgfhMEG4jXP0B3bS7lxQxFff+o0J1v7A/7GzCAuzuqiTL5wyzqqT3fwj7+2cjUZH8TF+cT1lWwuzeZcxxDL8mJ3DVK4mE9BXB8xKeIMVY35NRBz4Y9qyk5L4u7/OvCm0EOwZld56clkpJgFTPNx567lXL+2kHMdQ+RnmP4KhqQED//4vi0sSkpgZUFGtMVxPHMqCFUNvQK4ISh6hicYGItfJ2xeRgrffP9W6jqH+MIjR98U+mpCXINDRPjb924mPyOFFeZhFzQrCjJ45O6r+PzN66ItiuMxQ44o0GhMKFy5Mo97b1jD1586zc6KXD54xfLp7xq6h9m2bHEUpXMP+RkpPPLHV0VbDNexdklWtEVwBWbZZRRo6rHCPMsWx6+CAPijt6zkLasL+NJjJ6bTWU9MeWnpHY3b2dVCKMtNi+vBhsE+jIKIAk091gyiNI5rHQB4PMI/vG8LeenJ/OGP9tM9NE5L7whTXjUPPIPBARgFEQWaekbIXpQUd7WWA5Gbnsx3/+d22gfG+KMf7ed8h0nBbDA4BeODiAJNPcOUmqyb01xWlsPX3rOZe35yaNr8Fo8RXgaD0zAziCjQ2DNiFMQsbt9awsfesoLm3hGSEzwUZZo03wZDtDEziAijqjT1DFO1uiDaojiOT9+4ltqOIXqGx/HE0Qpzg8GpGAURYbqGxhmd8JoZRAASPML37txOiMXnDAaDTRgFEWH8NvbSOA9xnQsRIcFMHgwGR2B8EBHGv0gu3kNcDQaD8zEKIsKYGYTBYHALtikIESkTkb0iclJEjovIJ33bc0XkGRE56/sbMKeCiNSJyFEROSQir9klZ6Rp6hlmcVqSSaxmMBgcj50ziEngXlVdB+wC7haR9cBngWdVtRJ41vd5Lvao6hZV3WGjnBGlqWfEzB4MBoMrsE1BqGqrqh7wvR8ATgIlwG3AD30/+yFwu10yOJFGs0jOYDC4BJmdatmWk4iUA88DG4EGVc2Z8V2Pqr7JzCQitUAPoMD3VPW+OY59F3AXQFFR0fYHH3xwQTIODg6SkWFvymRV5a5nhrl+WRLvX5ts67lCJRLtdzKm/ab98dr+PXv27J/TSqOqtr6w6lrvB97t+9w76/ueOfYr9v0tBA4D117sXNu3b9eFsnfv3gXvGywX+kd0+Wd+qT98qdb2c4VKJNrvZEz790ZbhKgSz+0HXtM5nqm2RjGJSBLwEPCAqj7s23xBRJb6vl8KtAfaV1VbfH/bgUeAnXbKGglej2AyJiaDweB8bAulEavY6/3ASVX9xoyvHgU+BHzV9/cXAfZNBzyqOuB7fwPwZbtkjRQmxNVgcDcTExM0NTUxOjoabVFCJjU1ldLSUpKSgs8ibWes5W7gTuCoiBzybfs8lmL4qYh8BGgA7gAQkWLgX1X1ZqAIeMRXUDwR+C9VfdJGWSPC9CI5M4MwGFxJU1MTmZmZlJeX43s+uQJVpauri6amJioqKoLezzYFoaovAHP14PUBft8C3Ox7fx64zC7ZokVTzwh56cmkJZs1EAaDGxkdHXWdcgArhU1eXh4dHR0h7WdWUkcQUwfCYHA/blMOfhYit1EQEaTZLJIzGAwuwiiICOH1qrWK2iTpMxgMLsEoiAjRMTjG+JTXzCAMBoNrMAoiQjT1mAgmg8Fw6Xz/+99ny5YtbNmyBY/HM/3+U5/6VNjPZcJpIsT0GogcoyAMBsPC+ehHP8pHP/pRmpubueqqqzh06NDFd1ogZgYRIVp6rYU1xUZBGAyGMHDs2DE2bdpk6znMDCJCtPSOkJOWRLqpA2EwxARfeuw4J1r6w3rM9cVZfPEdG4L67dGjR9m4cWNYzz8bM4OIEC29IyzNNrMHg8EQHgLNIP7iL/4irOcww9kI0dw7YhzUBkMMEexI3y6OHj3Kn/7pn05/bmtrY3JyMqznMDOICNHSO2L8DwaDISx4vV7Onj3L2rVrp7cdPHiQLVu2hPU8RkFEgMGxSfpHJ42CMBgMYaGmpobS0lJSUlKmtx06dMgoCDfS2muFuC7NTo2yJAaDIRZYvXo1J06ceMO2mpoaKisrw3oe44OIAM0+BVFiZhAGg8Em7r///rAf08wgIoBZA2EwGNyIURARoLVvhASPUJiZcvEfGwwGg0MwCiICNPeOUJSZQmKC6W6DweAezBMrApgQV4PB4EaMgogALb2jRkEYDAbXYRSEzXi9SlufURAGg8F9GAVhM51DVqGg4hyzBsJgMLgLoyBsZjrE1STqMxgMLsMoCJtp8S2SMyYmg8HgNoyCsJkWs4raYDCEEVNyNIZo6R0lLTmBrEWmqw0Gw6VjSo7GEP41ECISbVEMBkMMYUqOxgAtfWaRnMEQk/zqs9B2NLzHXLIJ3v7VoH5qSo7GAC29I5SYEFeDwRBmzAzC5YxOTNE5OG5qURsMsUiQI327mF1ydGJigi9+8YsMDw8zPj7Od77znUs+h5lB2Ehbn0nzbTAYwk+gkqP33XcfIyMj5OTkMDg4GJbz2DaDEJEy4D+AJYAXuE9V/0lEcoGfAOVAHfB7qtoTYP+bgH8CEoB/VdXoqusF8PoaCGNiMhgM4SNQydGDBw/y7W9/+w3bLhU7ZxCTwL2qug7YBdwtIuuBzwLPqmol8Kzv8xsQkQTg28DbgfXAB3z7ugpTSc5gMNhBoJKjt912Gx/+8If59Kc/zZNPPhmW89g2g1DVVqDV935ARE4CJcBtQJXvZz8EqoHPzNp9J1CjqucBRORB334nsIFffPV/UTZ6hsO/TQjrcVdMKQ8meyn7xbfB4WGuW3p7oTYn2mJEDdN+0/6g2r/x09DpQNdt0iLe8Y538I53vCOsh41IS0WkHNgKvAIU+ZQHqtoqIoUBdikBGmd8bgKumOPYdwF3ARQVFVFdXR2yfMk6TrIHBG/I+85HYgJkJwv9fX1hPa4dTE1N0dvbG20xooZpv2l/MO33er1MTk1FQKLQ8Oo4YwMDF/3d6OhoSM9I2xWEiGQADwH3qGp/kAvGAv1IA/1QVe8D7gPYsWOHVlVVhS5kVRXV1dUsaN8YwbTftN+0v+qiv2s9eZLEorUX/V00SA7iN6mpqWzdujXoY9oaxSQiSVjK4QFVfdi3+YKILPV9vxRoD7BrE1A243Mp0GKnrAaDwWB4I7YpCLGmCvcDJ1X1GzO+ehT4kO/9h4BfBNj9VaBSRCpEJBl4v28/g8FgiCqqAY0Z06MCuAAABfFJREFUjmchcts5g9gN3AlcJyKHfK+bga8CbxORs8DbfJ8RkWIReQJAVSeBjwNPASeBn6rqcRtlNRgMhouSmppKV1eX65SEqtLV1UVqamgh93ZGMb1AYF8CwPUBft8C3Dzj8xPAE/ZIZzAYDKFTWlpKU1MTHR0d0RYlZFJTUyktLQ1pHwfGaxkMBoMzSUpKoqKiItpiRAyTasNgMBgMATEKwmAwGAwBMQrCYDAYDAERt3nj50NEOoD6Be6eD3SGURy3Ydpv2m/aH58sV9WCQF/ElIK4FETkNVXdEW05ooVpv2m/aX/8tn8ujInJYDAYDAExCsJgMBgMATEK4nXui7YAUca0P74x7Te8CeODMBgMBkNAzAzCYDAYDAExCsJgMBgMAYl7BSEiN4nIaRGpEZE31ceONUSkTET2ishJETkuIp/0bc8VkWdE5Kzv7+Joy2onIpIgIgdF5Je+z3HTfhHJEZGficgp33VwZZy1/0991/4xEfmxiKTGU/tDIa4VhIgkAN8G3g6sBz4gIuujK5XtTAL3quo6YBdwt6/NnwWeVdVK4Fnf51jmk1ip5P3EU/v/CXhSVdcCl2H1Q1y0X0RKgE8AO1R1I5CAVW8mLtofKnGtIICdQI2qnlfVceBB4LYoy2Qrqtqqqgd87wewHg4lWO3+oe9nPwRuj46E9iMipcAtwL/O2BwX7ReRLOBarGJeqOq4qvYSJ+33kQgsEpFEIA2rWmU8tT9o4l1BlACNMz43+bbFBSJSDmwFXgGKVLUVLCUCFEZPMtv5R+DTgHfGtnhp/wqgA/h3n4ntX0UknThpv6o2A38HNACtQJ+qPk2ctD9U4l1BBCpoFBdxvyKSgVUv/B5V7Y+2PJFCRG4F2lV1f7RliRKJwDbgX1R1KzBEHJlTfL6F24AKoBhIF5H/GV2pnEu8K4gmoGzG51Ks6WZMIyJJWMrhAVV92Lf5gogs9X2/FGiPlnw2sxt4p4jUYZkUrxORHxE/7W8CmlT1Fd/nn2EpjHhp/1uBWlXtUNUJ4GHgKuKn/SER7wriVaBSRCpEJBnLWfVolGWyFRERLPvzSVX9xoyvHgU+5Hv/IeAXkZYtEqjq51S1VFXLsf7fv1HV/0n8tL8NaBSRNb5N1wMniJP2Y5mWdolImu9euB7LDxcv7Q+JuF9JLSI3Y9mkE4B/U9WvRFkkWxGRq4HfAkd53Qb/eSw/xE+BZVg30R2q2h0VISOEiFQBf6aqt4pIHnHSfhHZguWgTwbOA7+PNViMl/Z/CXgfVkTfQeAPgAzipP2hEPcKwmAwGAyBiXcTk8FgMBjmwCgIg8FgMATEKAiDwWAwBMQoCIPBYDAExCgIg8FgMATEKAiDYRYikicih3yvNhFp9r0fFJHv2HTOe0Tkf/neV4vIjgC/2SQiP7Dj/AZDIBKjLYDB4DRUtQvYAiAifwkMqurf2XU+X9K4/421onk+uY6KSKmILFPVBrvkMRj8mBmEwRAkIlI1o37EX4rID0XkaRGpE5F3i8jXROSoiDzpS2eCiGwXkedEZL+IPOVP5zCL64ADqjo5Y9sdIrJPRM6IyDUztj+GtQLcYLAdoyAMhoWzEitt+G3Aj4C9qroJGAFu8SmJbwHvVdXtwL8BgVbq7wZmJw9MVNWdwD3AF2dsfw24BoMhAhgTk8GwcH6lqhMichQrVcuTvu1HgXJgDbAReMZK+0MCVorp2SzljcWLwEoiB5biKJ+xvR0rC6nBYDtGQRgMC2cMQFW9IjKhr+et8WLdWwIcV9UrL3KcESA10LGBKd54n6b6fm8w2I4xMRkM9nEaKBCRK8FKsy4iGwL87iSwKshjrgaOhUk+g2FejIIwGGzCV8b2vcDfishh4BBW7YHZ/AqrDGgw7AEeD4+EBsP8mGyuBoMDEJFHgE+r6tl5fpMCPAdcPSviyWCwBaMgDAYH4CvgU6Sqz8/zm0qgRFWrIyaYIa4xCsJgMBgMATE+CIPBYDAExCgIg8FgMATEKAiDwWAwBMQoCIPBYDAExCgIg8FgMATk/wMikqfdlIooVgAAAABJRU5ErkJggg==\n", |
|
|
"text/plain": [ |
|
|
"<Figure size 432x288 with 1 Axes>" |
|
|
] |
|
|
}, |
|
|
"metadata": { |
|
|
"needs_background": "light" |
|
|
}, |
|
|
"output_type": "display_data" |
|
|
} |
|
|
], |
|
|
"source": [ |
|
|
"import numpy as np\n", |
|
|
"import matplotlib.pyplot as plt\n", |
|
|
"%matplotlib inline\n", |
|
|
"\n", |
|
|
"# Number of days, n, to be simulated\n", |
|
|
"n = 4 # days\n", |
|
|
"t = range(1,24*n) # time span in hours\n", |
|
|
"Ta = np.full((1, 24*n), 20.0)[0,:] # outdoor temperature\n", |
|
|
"Ti = np.full((1, 24*n), 20.0)[0,:] # indoor temperature\n", |
|
|
"q_int_t = q_int * np.array([0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0]*n)\n", |
|
|
"q_sol_t = q_sol * np.array([0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0]*n)\n", |
|
|
"P = (q_int_t + q_sol_t) # Combined heat gains\n", |
|
|
"\n", |
|
|
"delta_t = 1\n", |
|
|
"inputs = np.vstack((Ta, P))\n", |
|
|
"\n", |
|
|
"# Function that will calculate the indoor temperature for a given value of R and C\n", |
|
|
"def Ti_calc(inputs, R, C):\n", |
|
|
" y = np.zeros_like(Ti)\n", |
|
|
" y[0] = Ti[0]\n", |
|
|
" for t in range(1,len(Ti)):\n", |
|
|
" y[t] = 1/(C/delta_t+1/R) * (C/delta_t*y[t-1] + 1/R*Ta[t] + P[t])\n", |
|
|
" return y\n", |
|
|
"\n", |
|
|
"# Simulate indoor temperature\n", |
|
|
"Ti = Ti_calc(inputs, (1/HLC), C_env)\n", |
|
|
"\n", |
|
|
"# First plot of the data to make sure that nothing is missing\n", |
|
|
"plt.figure()\n", |
|
|
"plt.plot(Ti,label='$T_i$')\n", |
|
|
"plt.plot(Ta,label='$T_e$')\n", |
|
|
"plt.legend(loc=\"lower right\")\n", |
|
|
"plt.grid()\n", |
|
|
"plt.xlabel('Time (h)')\n", |
|
|
"plt.ylabel('Temperature (°C)')\n", |
|
|
"plt.show()" |
|
|
] |
|
|
}, |
|
|
{ |
|
|
"cell_type": "markdown", |
|
|
"metadata": {}, |
|
|
"source": [ |
|
|
"### Calculating building heat loss and thermal capacity from indoor temperature\n", |
|
|
"\n", |
|
|
"The simplified dynamic heat balance,\n", |
|
|
"\n", |
|
|
"$C'' \\frac{dT_i}{dt} = q''_{int} + q''_{sol} - H'' (T_i - T_e)$\n", |
|
|
"\n", |
|
|
"can also be used in an inverse modeling procedure, if the indoor temperature is known. \n", |
|
|
"\n", |
|
|
"In the example below we estimate the two model parameters $C''$ and $H''$ by optimizing the sum of squared residuals between a model prediction and indoor temperature measurements. In the following example $T_i$ is the simulated indoor temperature, and therefore showing a perfect fit when applied to the same model." |
|
|
] |
|
|
}, |
|
|
{ |
|
|
"cell_type": "code", |
|
|
"execution_count": 667, |
|
|
"metadata": {}, |
|
|
"outputs": [ |
|
|
{ |
|
|
"name": "stdout", |
|
|
"output_type": "stream", |
|
|
"text": [ |
|
|
"R: 0.25 K/(W m^2)\n", |
|
|
"H: 4.07 W/(K m^2)\n", |
|
|
"C: 48.48 Wh/(K m^2)\n" |
|
|
] |
|
|
}, |
|
|
{ |
|
|
"data": { |
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOy9d5gkV3no/Xu7J+ec487M5qzd1UoLyllCsnURIMIVGFu2EQYM/mxs7v2uccT2ta5zwMgEo2s+jCSDQSCBAkISknZXm2Z3dnLOOed+vz+qR5pd9cz2zHR1V0+f3/PUM93VdU69dXar3jrnTaKqGAwGg8FwOa5QC2AwGAwGZ2IUhMFgMBh8YhSEwWAwGHxiFITBYDAYfGIUhMFgMBh8EhVqAQJJVlaWlpWVhVoMg8FgCBtOnjw5oKrZvn7bVAqirKyMEydOhFoMg8FgCBtEpHWl38wSk8FgMBh8YhSEwWAwGHxiFITBYDAYfGIUhMFgMBh8YhSEwWAwGHxiFITBYDAYfGIUhMFgMBh8YhSEwWAw2MDC/BxvfOdRZmemQi3KurFNQYhIsYi8ICI1InJeRD7t3f8XInJRRM6KyFMikrZC+xYROScip0XERL8ZDCFkdLCXzi9upfqV/wq1KGHDuRe+zZHqL3LuJ98MtSjrxs4ZxALwOVXdARwFHhGRncCPgd2quheoA353lT5uVNX9qnrIRjkNBsMVaDzxLIXay8TZ74dalLBhtvFlABbaT4ZYkvVjm4JQ1W5VfdP7eRyoAQpV9VlVXfAe9hpQZJcMBoMhMMw1/gyA1KGzIZYkfMgctBY+wnnMgmKDEJEy4ADw+mU//RLwwxWaKfCsiJwUkYdX6fthETkhIif6+/sDIa4hAjjz/LcZ7O0ItRhhQ5b3YVc218DC/FyIpXE+E2PDbJlvYE7dYT1mtisIEUkCngA+o6pjy/Z/AWsZ6vEVmh5T1YPAnVjLU9f5OkhVv6yqh1T1UHa2z4SEBsMlDPZ2sO+lX6H+yT8KtShhwfjoEOULTbRLAfEyR+vFN0MtkuNpPvUCblHOpN8W1mNmq4IQkWgs5fC4qj65bP9DwD3Ah1RVfbVV1S7v3z7gKeCInbIaIofW0y8AkDZ0JsSShAfNbz6PW5SunR8HYLD21RBL5Hwm6l5iQV2k3/gbAAzW/TzEEq0PO72YBHgMqFHVR5ftvwP4HeBeVfXp/yUiiSKSvPQZuA2otktWQ2Qx1/QKAGVz9czPzYZYGucz2fAz5tXN7jt+mRGSoDN8ja7BIrXvOE3RlVTsuYYxEsN2zOycQRwDPgLc5HVVPS0idwF/ByQDP/bu+ycAESkQkae9bXOBl0XkDPAG8ANV/ZGNshoiiIzBkyyoiziZp7XmeKjFcTxpfcdpjq4gMTmNtrjtZI6dD7VIjmZmepLKuYsMZV6FuFy0xm4jczQ8x8y2gkGq+jIgPn562se+pSWlu7yfm4B9dslmiFymJkYpn2/kdMqNHBp/jsHaV6nc965Qi+VYZqYnqZir5c389wMwlbWPne1fZWpilISk1BBL50yazvyMnbJAXKVlNp3I2seOjq8zPTlOfGJyiKVbGyaS2hBRNJ15iWhZJPrABxgiBVeYTv2DRdPpl4iRBeIqLCUaV36EKPHQev61EEvmXEYvvghA2YGbAIgrPRS2Y2YUhCGiGK97GY8KZQduoi1+Jznj4Tn1DxajtT8FoPzgzQAU7Tpm7W8Iv4ddsEjseYNmVylpWXkAFO9+NwAjYThmRkEYIoqknjdocZeSmp7FdM5+ihc7GBsZDLVYjsV62JWRmpkLQFZeMT1kE91zOsSSOZOF+Tkqps/Tl3HVW/uyCkrpJZOoMBwzoyAMEcPC/BxbZi7Qn3EQgKQtV+MSpfXsyyGWzJm8/bA7eMn+rqQd5E2YmZcvmqtfI1FmcJdfe8n+rsQd5IXhbNUoCEPE0HLhjUtu3pI91tR/sunyAH8DvP2wiyo/dsn+uZz9FGovw/3dIZLMuQxeeBGAkv23XLJ/Jmc/RdrN6GBvCKRaP0ZBGCKGgQvWenrxPst4mJqRTbsUENt3KpRiOZbBmpcAKN5/8yX7kyuOAtBWbWZelxPTfZwesskpLL9kf/IWK863rfqVUIi1boyCMEQMMZ2v0002uUUVb+3rSdlN8dQF1OMJoWTOJLr7BD1kveNhV7rnWjwqTDW/ESLJnEvhxHk6k3e/Y3/JnnfhUWEizGarRkEYIgL1eCieOEtnyqXhNZ6Cq8hihJ72+hBJ5lwKxqvpTHrnwy4pJZ1WdwmJZuZ1Cb0djeQyyHz+O6sTpKRl0uYuJqE3vNyqjYIwRARdLbVkM8xi0dFL9mdstewRXedNfqHlDHS1kk8/8/lX+fy9L20fZTPn8SwuBlky59Jx1lqSy9juO/CyL20f5TMXwmrMjIIwRATd1Zb9IXvnpUmBS3ceYVajmW8Nr6m/3bSfexGAtG3HfP7uKjlKClO01ZlZxBLzra8xo9GU7Trq+4Diq0lhkva68HF3NQrCEBEstr3BlMZSuv3SN+KY2DiaoytNZtfLmG1+nTmNonz3NT5/z999PQC9XsVrgPTB0zTHVBETG+fz9/w9NwDQe+GlIEq1MYyCMEQE6SPnaI7dhjvqnenHRrIOsmWujpnpyRBI5kxSB0/RHF1JbFyCz98Lt+y0UpV0GEM1wOzMFOXzDYxmHljxmKItuxgmBWkPn9mqURCGTc/M9CRl842MZe33+XvslmuJkQVazoaXC6JdzM3OUD5Xz3Cm7/ECrCylCbvJGwvfcpqBpPncq8TIArHlV694jLhctCTsJm80fGarRkEYNj0t1T8nRhaJK/N985bsuwGAkdqfBVEq59Jy/nXiZJ7o0hXW0r3M5h2iWLtM6VZgpM56uSj2LiOtxIx3zMIlyNAoCMOmZ6TO8lBaSpp2OZm5RbRLAXE9pjYEwFCtFQBXtNdnld+3SNtmjWfbmRftFsnxxHSfpJtssgpKVz0udatl9G8NkzEzCsKw6YnuPknPFW7e7tT9lE5Vm4A5ILrrBH1kXBJQ6IuyvceYUzczzeFZTjOQFE5U05W854rHbdn7LubVzXRTeLhVGwVh2PQUTJynM2nX6gcVX00647TVmzX1/IlqOnwEyF1OXHwizdFVpA1EtqvrWwFyBe8MkLucuIQkmqMrSA2TMTMKwrCpGehp8wZ8HVz1uLwlF8TqF+0XysEM9LRRoH3M+YgG9sVw5gG2zNUxO+OzvHxE8FaA3Db/KhMOZRxgy+zFsKiHbhSEYVPT7r1506quXfW44sq9XhfE8CvqEkjazlhxDVcaryVit1xDrMzTXB25y0xvB8it7MG0nOiya4gLkzEzCmITcP5PruPn3/ifoRbDkcy0vM68uilbIeBriXB0QbSD2ZbXmFM3ZXv8UxBLXjsjtZGb2TV16CwtqwTIXc6S8X/oovO95oyCCHNGBnrYNXeGnNbvh1oUR5IycIaW6C3EJSRd8VjjtgnJQ9W0RpUTF5/o1/FZBaV0SS4xXZHpAbYwP0fpXCMj6Ve22SyRW1RBD9nEdDk/yNAoiDBnKb98+UIzo8MDIZbGWSwuLFA2W8tQ2pW9SwDStltvdpHqtulZXKR0ppahNP8fdgBdyXspnTwbkR5g7XWnSJBZoop8JzVciY6UfZROnHH8mBkFEeZMtlhvbi5Rmt/8SYilcRZttSetCnIlR/w6vnzvMWY1mtnGyFwu6WyqJlmmkYKV00X4wlN6jExGaQujJHSBor/OSpuRs331oMLL8ZR4x8zhXnNGQYQ5cX1n6SabOY1ipj58koAFg/6Llq957k7/vEti4xJoiqkifTA8XBADTW+NZTTN2ra6veZy8vdZ5TV7zj4XcJmcjnacZFzjKarwb5a6RP4+q0pfz1lnv9QZBRHmFE7V0Jmyj4aY7WQMROY68Epo1ynGSKRoyxViIJYxkn2ILfP1TE2M2iiZM1noeJMZjaZk29pmEEVbdtFPOlHt4RH8FUjSRy/QFrsVl9u9pnZFFXsYIA13m7PHzCiIMGagq5UchljI289ozmG2zDcwMTYcarEcQ8bIedpiqxCX///NE7feQLQs0vjm8zZK5kxSh6tpia4kKjpmTe3E5aI1+SAl46ccv6YeSOZmZyibb2Isc22zB1gaswOOHzOjIMKYDu+SQGrFEZK2XU+UeGg+9UKIpXIGc7MzlC60MJHu/+wBoOKqm5lXNxO1L9ojmEOxvHEa1uSNs5zFkmvJZpiOpvMBlsy5tNYcJ0YWiCn2L6jwchaKryWHITqbLgRYssBhFEQYM91ynEUVSncdpeLgTdaDrc4UcAFou3iCGFkgumRt3iWJyWk0Rm8lvTd8cvYHgvb6MyTILO6i1SPOVyJvr2WH6D7j7DX1QDJUZwVV5u/0L2bkcvL2WHaILgePmVEQYUzCwFna3CUkJKWSkJRKU3QlaX3O960OBkP11jjkrdG7BGA45wgV83VMjo8EWizH0l9rPexytq/vYVdStZcB0nC1Rk5NDek+xTDJ5JdUrat9ybYDDJGCtDl3zIyCCFPU46FoupaBlJ1v7RvKPkzFXC0zUxMhlMwZLBmoC8p2rLlt0jbLDtEUQXYI7XyTSY2juHLt6+lgram3Je2jaMzZa+qBJGv0PO1x29dk41qOuFy0JO6neNS5XnO2KQgRKRaRF0SkRkTOi8invfv/QkQuishZEXlKRNJWaH+HiNSKSIOIfN4uOcOV3o5GMhnFk/921a+EqncTI4s0GDsEGaMX1mygXiIS7RDpI9W0xFat2RtnOfPF15LHAF0ttQGUzJlMT45TstjGZNb6FOoSc0XXkEc/3a3OHDM7ZxALwOdUdQdwFHhERHYCPwZ2q+peoA743csbiogb+HvgTmAn8KC3rcFLd43lHpde9fYSSvnBW/GoMF4b2fEQSwbq8Yz1GVwTklJpjNlGRl9k2CGWvHHGM/ZuqJ8c75p655kfB0IsR9N6/jWixEN86eEN9bM0Zh2nnGmHsE1BqGq3qr7p/TwO1ACFqvqsqi54D3sNKPLR/AjQoKpNqjoHfAu4zy5Zw5GZ1pPMq5uSHW//B01Jy6QpagspPRGekdRroI4pXp/BFSLLDvGWQX8D4wVQuu0gwyQjLc5dUw8UIw3Wy0PhrvXZbJYo23GIEZLQFmdG7wfFBiEiZcAB4PJXsl8CfuijSSHQvux7h3efr74fFpETInKiv79/48KGCUmDZ2mNKntHUrWBrCNUzV5genI8RJKFno0YqJdI2nYDUeKh8eTmjw4e9KaLyNtxbEP9uNxumhP3Uzj2ZiDEcjRRPafpI4PsgrIN9eNyu2lO2Evh6MnACBZgbFcQIpIEPAF8RlXHlu3/AtYy1OO+mvnYp776V9Uvq+ohVT2UnZ0dCJEdj3o8lM7WMZj6Th//hO23ECMLNJx05pQ1GGzEQL1ExcGbmFM3k3UvBk4whyJdbzJKIgVl2zbc11zRtRRoH13NFwMgmXPJnqihK2Hj4wUwW/wuCrXXkWNmq4IQkWgs5fC4qj65bP9DwD3Ah1TV14O/Ayhe9r0I6LJT1nCiq6WWFCZhmYF6icrDtzKnbiZqIldBbMRAvUQk2SEyxmpoi926ofFaIv/AHQB0nHx6w305lamJUYoXO5nOXJ+N63LyD9wJQMebvhZTQoudXkwCPAbUqOqjy/bfAfwOcK+qrlSn8DhQJSLlIhIDfAD4nl2yhhu93iWB9Ip3RnAmJKXSELuTrL7ItENs1EC9nJGcq6mYr9/U6Uvm52YpWWhlMj0wPiAlW/fTRwZRLS8GpD8n0nbxBC5R4orf+YK2HpbGzN3ivCBXO2cQx4CPADeJyGnvdhfwd0Ay8GPvvn8CEJECEXkawGvE/iTwDJZx+9uqGjkx/FdgtuMMC+qiZLvvKOHR/GNULDQyMtATZMlCTyAM1Esk77yZKPHQ8MaPAiCZM2mvO02MLBBVuC8g/YnLRWvqEbZMnMSzuBiQPp3GaJNlL8jbtjEPpiWsMTvMlnHnjZmdXkwvq6qo6l5V3e/dnlbVSlUtXrbv17zHd6nqXcvaP62qW1W1QlX/2C45w5H4wQt0uItWrJKWvvtWXKI0HXfelNVu3jK4bsBAvUTVVTczpbHM1m7e5brBhhMAZFcF5mEHIBU3ksYETWFQc3ld9JxjlETyitcXQe2TLTeQzhhN1c6a+ZtI6jAkf7qegaStK/5euf86JjSe+YYIDJjrObNhA/USsXEJ1MfvJX9wkz7ogMXus0xrDEWVG4uBWE7ZEes9r//MMwHr00mkj12kI6YyIDabJcoP3w3AwFlnjZlREGHGcH83uQyykLNyBGdUdAz1iQcoHNr8BtbLSRu9SHtMRcBu3umS6ynxdDo20nWjJA/X0B5dhjsqKmB9ZuWV0OwqJanzZwHr0ykszM9RMt/MeIBsNktkFZTS4iohodNZ8RBGQYQZHTWWj39S6epFXWaL30WR9kRE2oMlFhcWrJs3bXvA+szzepi0n9h8Xjnq8VA818BwSuDGa4ne7GvYOlO96fKCdTScJU7mcRcEbsa1RE/WUbZOn2VmejLgfa8XoyDCjMk2KwipcPvqdZbz9y+5G0aOHaKzqZp4mcOVH7ibt3TbQcvDpPnFgPXpFHra6y136bzAP+zit91MrMxTf2JzBRoO1Fs2m6wA2myWiNt2M3EyT4ODgjONgggzovqq6SWT9Oz8VY8r2XaAftId6TpnF/0NlndJxpaNezAtIS4XrWlXUzFxgsWFhSs3CCN6aq0StakBHK8lKg/ftinjcRY6TzOr0QG12SxRefh25tXN+AXnjJlREGFG1kQt3QkrG6iXEJeLltTDbBk/4TjXObuY6zzLvLopXmNN5Sux5JXTeHZz5Ria6TiDR4WS7euriLYaiclpNMTuJLtvc41Z0kgNbdFlRMfEBr7vlHQaYraT2eccpwijIMKImakJihc7mMnwz0AmFTeRztime7CtRMLQBTrcxcTGJQS03/IjlofJkMM8TDZK3EA17e5CEpJSbel/NP8YWxaaGO7vtqX/YKMeD0WzDQwnBybFhi9G8o9ROV/vmBgmoyDCiLaLJ3GLElvsX1DTlqP34lFh4PQPbJbMGeRPNzCYFEDfdC+ZuUU0uCtI7tpcadTzpuoZSLzybHS9ZOy5HZcoja//l23nCCa9nU2kMYHaYLNZImPvnbhEaXDImBkFEUYMN1oGstytV/t1fEZOIQ3RVWR0vmijVM5guL+bHIZYyHlnAsNA0J9zDVWzFzZN2o3RoX7y6WcuOzD5hHxRuf86hkmBumdtO0cw6bloeRCmlgfeZrPE0pipQ8bMKIhwoucs4xpPfqn/b32D+ddTNV+7aab5K9HpvXkTSwKTH+dyknfeZlXr2yRpN5bcpe0aLwB3VBSNKVdTMfb6pjDwT7efxqNC8Y7AezAtsTRmlaOvOcJ2aBREGJE6Wkt77NoiODMP3GNN81/b3LkOJ9pOA1AQoPw4l1N1+FYr7UbN5nAbHm+xPL4KruAuvWG23kY6Y9SfDn9vuriBajpd+SQm+6ySHDi8Y9ZwJvSBhkZBhAlWEFgTY6lrC2qq3PduhkiB+s1dBtLdd55+0snM9VWgcOPExiVQm3SYssGXUY/HlnMEE3dfNf2kk5VXfOWDN0DVNfexqMLIJrCD5U7X02ejzWaJyqP3sqjC4Knv236uK2EURJjQ2VRNgsziKlhb1k2X203TJprmr0TGRB3dcZW2nmN+y63kMkjT+TdsPU8wyByvpcvm8QJIzcylLmYHmd3hPYMYHR6gQPuYy7bHxrWctKw86h0yZqsqCBEpEpHfEpHvishxEXlJRP5BRO4WEaNcgsiGgsCWpqynN5cXzhJzszMUL7QxmbHxBH2rseWaXwSg/2R4L9fNzc5QvNjOlM3jtcRo4Y1ULTYw0NMWlPPZQVeddf8lFNvnwbSc4cLr2bpQx0BP+5UPtpEVH/Ii8lXgX4E54M+AB4FPAD8B7gBeFpHrgiGkAea6qllQF0Vb125UXJqyDp3ZfPmEYKmmwSLRhfbevFkFpdS7K0nreN7W89hNZ8MZomWR6Hz7PJiWk33wHgCafv7doJzPDsZazwKQW+m7BkugyT7wHgCaQ2w7XG0W8Jeqepuq/o2qvqqqDaparapPqupvADdgyoAGjbihWjrdBcTFJ665bVpWHvXR2x0xZbWDoUbr7S67wv6bd6DgRrbOXwxrr7DB5jMAZG6xz4NpOVt2H6WfdKIaw9gO1neBcY0nt6giKKer2HMN/aTjCvGYraYgekXkHSG7IrJLRLJVdU5VG2yUzbCM7OlGBhO2rLv9cOH1VM7XM9jbEUCpnMFi9zlmNJrCCvvfiDMPvMcKZHr1KdvPZRfzXdXMq5vCysBUkbsS4nLRnH4tVeNvMD83G5RzBprk0To6Y8oDWgNiNcTlojntGqrG32Bhfi4o5/TFalf7t0C2j/1FwF/bI47BF1MToxR4epndwJpx9kHrwdb46pMBlMwZWPlxyomKjrH9XJX73sUAabganBHItB7iR+rodBcSExsXtHNGb7+dZJmm7rhzEtH5i3o8FM03M5oS+Cj91YjafjspTFJ3MnRLmqspiD2q+o41CVV9BgiOpcYAQEfdaatIeuH6PSgq9lxLD1lE128OP/7lFMw2MZJkv0cOeL3C0q6lavz1sH0bzpluZDAxOEslS1Rdcy9zGsX4WWekkFgL/d2tVlr0nMAWCboSW6+9jzmNYuzUfwb1vMtZTUFEr/M3Q4AZabXWjLMq1h/iLy4Xrdk3sH3yBNOT44ESLeQM9naQwRie7OB45ABE77iDFKaoC8NaBxNjw5a7ZmbgiwStRlJKOjXxByjuez7s4kh6vB5MSUHyYFpiacyK+l4I2ZitpiDqReSuy3eKyJ1Ak30iGS7H03OeaY3ZcJ3lxH33ES9zXHwlfL1JLqe74RQAicUrl2ANNNbbsJvxM+Hn7tpRZxWcirfZ48sXMxV3Uqi9tNQcD/q5N8JUh+XBVLQtOB5My5mpuJMi7aHl4smgnxtWVxC/CfyViHxNRH7Du30dy/7w6eCIZwBIHKmjI6pkw3WDtx25nTESWTgfftP8lZhsrwYgv9K+BGqXY73ZHaQkDN+Gx1qs2Wh2RWBrZvhDxbsfwKNCz+tPBP3cG8E9UEMfGaRm5gb93G+P2XeCfm5YRUGoah2wB/gpUObdfgrs9f5mCBL5AVpjj46JpS7lWipHXwmpZ0RA6a9hlEQybU4ZcTmzVXdToL1hF1Xt6b3AlMauKeFjoMjKK6EuZgfZneFlqE6faKQnbv0ehBshK6+EuujtZHWEZsxWC5QTVZ1V1a+q6ue827+q6szyY4IjZuQy3N9NFiMsBmiN3bXzHtIZpzYMvUl8kTJWT2d08NwPl6h41wMsqtD3xn8E9bwbJXm0jo7oUlxud0jOP1JyK5WLjXS11Ibk/GtlYX6O4oU2plKD68G0nJHS26habKCnrT7o517trnrBu6xUsnyniMSIyE3e5aaH7BXP0OVdM04oCswa+7Zjv8CsRjN+OnSeEYFCPR4K51sZTwmOB9NyMnOLqI3dTV5n+AR/qcdDwVwzI8mhe9gVXvMAAG2vhodi7Wy6QKzM4w5S1Lkvlsas5ZXgj9lqCuIOYBH4dxHpEpELItIM1GOl3fg/qvq1IMgY0Uy0Wway/K2BMZAlJqdxMeEgJf2h84wIFG+7HwbPg2k5Y2V3UO5ppb3hXEjOv1YG+zpJZwxPdnDdNZdTXLmHFlcJyS3hUb51sMlygkgvD07UuS+sMSsmKQRjtpoNYkZV/0FVjwGlwM3AAVUtVdVfUdXTQZMykum7wCiJZOWVXPlYP5mtvJMC7aP5Qnh5k1xOr9eDKSlAs6u1Uvau9wPQ+eq3Q3L+tdLtnY0mBdHjy6ccBbewffZcWKQrme2qZlGFoqrQKQiA7vyb2T5zltHB3qCe16+FW1WdV9VuVR2xWyDDpaSON9AZsyWga+xbjr2XRRV6Xw+PB9tKTHZYb+75Ibp580qqqI+qIr0tPN6GJzsCOxtdL1mH7sctSv3LofHMWQtxQxfpcuUTl5AUUjkyr7qfKPFQ97PgjplJ2e1g1OOhcK6Z8QCH+GflFXMxdi8FnT8K62UmV/9FBkklI6cwZDIMFN/GtoVaejsaQyaDv7j6LjBEim1Flfylcu8xesgmpi70BXGuRPZUI/0JwbdxXU7lvndZmRBqgxt7YxSEg+ntaCRZpm0J8Z+ofA+lno6wXmZKm2igO7Y8pDIUXvM+AFpedv5sLG2iga7Y0LhrLkdcLlpyb2Hn1HFGh/pDLc6KTE+OU+DpYTbIUee+cLndtOTeGvQx80tBiEipiNzi/RwvIsl+tCkWkRdEpEZEzovIp737H/B+94jIoVXat4jIORE5LSIn/L2gzUSvt0hQSmngo16rbvigtcz0838PeN/BwLO4SOF8GxNBTqB2OSVb91sGxCZnl9T0LC5SNN8a8vFaIuPqB4mRRepe/L+hFmVFOupO4RIltsD+KnL+kHH1B4iRRWp/Grx79ooKQkR+BfgO8M/eXUWAPz6SC8DnVHUHcBR4xJs+vBq4H/CnvNmNqrpfVVdUJJuZqY7zABTYsGackVPIhbj9FHeF5zJTT3sDiTKDhMiDaTndRXeyY7aavs7mUIuyIt2tdVbJ2tzQeTAtp2r/u+mUXGLrnJv2ZaTVsnFlhtCDaTlV+6+jS3KJqw3emPkzg3gEOAaMAahqPZBzpUZeo/ab3s/jQA1QqKo1qhoeUTIhJmqw1grxT8+ypf/prfdRpN00nvu5Lf3bSV+j5cGUasPsaq0UvfvDuERpevHfQi3KivR7iwSllITWg2kJcbloy7+DndOnGOrrDLU4Plnou8i8uikoD/1LCHgTbubdzs7pNxkZ6AnKOf1RELOq+lZeBhGJAnQtJxGRMuAA8PoaminwrIicFJGHV+n7YRE5ISIn+vudu565HlInm+mLDZx76+Vsvf5B5tVN/2vfsu0cdjHdac2u8quCn1Pocoqr9tHgriCj2bk5rma6LgCQX+mMt2GAnGs+SJR4qA/ikslaiBtpoNNdQHRMbKhFeYvsox+wvJmCtDTnj4L4qYj8HhAvIrcC/wH4fSeISBLwBNsAU4gAACAASURBVPAZVR1bg2zHVPUgcCfW8pTP+teq+mVVPaSqh7KzfdU3Ck/U46Fgvp3JZPvy9qdl5VETf5DSnmfCbpkpeqCGXjJtm12tlYHy97B1oY7OpvOhFsUnrsF6BkgjNcM598iWXUdocxWSVO/MrLhZ0y0MxYfWCeJyKvZcQ7sUkFAfnGUmfxTE7wD9wDngV4Gngf/hT+ciEo2lHB5X1TWVMlPVLu/fPuAp4Mha2oc7fV3NJMk0ZG+z9Twz2+6lQHtpOPOyrecJNGmTTfTGOefmLbvuwwC0veTMZabUiSZ6Y+ybja4HcbnoLLyTHbNnGehqDbU4lzA7M0WBp5vZ9NC7uC5HXC46Cu9gx8wZBnrabT/fqgpCRFzAOVX9F1V9QFXf6/18xSUmbyK/x4AaVX10LUKJSOKSp5SIJAK3YRm3I4a+JiuoKbHIXqPitusfZE7dDL7mXG+Sy1lcWKAoxAnULievpIqa6F3ktznPm0k9HvIX2phIDr2L6+UUXPtBq8b3Tx8PtSiX0NV0Hrco0bmhd3G9nLxrP4hblMaf2n/PrqogVNUDnLk8YZ+fHAM+AtzkdVU9LSJ3icgvikgHcA3wAxF5BkBECkTkaW/bXOBlETkDvAH8QFV/tA4ZwpZJrwdTXoW9heVTM7I5n3g1W3qfYXFhwdZzBYqetlriZB53rjOMh0uMVd5HmaeN5vNrMbXZz2BPOylMoVn2zkbXQ+mOq2h2lZLW5Cz7zVCL9T6aXuoMo/5yyncepsVVTEqD/Qk3/VliygfOi8hzIvK9pe1KjVT1ZVUVVd3rdVXdr6pPq+pTqlqkqrGqmquqt3uP71LVu7yfm1R1n3fbpap/vLHLDD9ksM6qc5Bjf9Sr7n0/OQxx4VVn3aQrMdBiuR8m2zy7WiuVN3yIBXXR86qzZmM9TZYHU+IGaprbSU/pe9g+f8FR9pu5ngt4VCisDL2XnC+6S+9jRxDGzB8F8UXgHuAPgL9cthlsJHm8ie7o0qDUOdh1w/sYI5HZE86a5q/EdNdFAPIrnHXzZuYWcSH+AKVdP8SzuBhqcd5istPyYMqrtHc2ul623PxLeFRoe+GroRblLWKG6+lx5YQ8B9NKlN/0MWvMXvyaree54tNHVX/qa7NVKgN5c62MJQXHCBsbl0BN5i3sHH2JyXHn52N0DdYxSGpISkBeibmdD1CgvVx849lQi/I2/bWMa3xAMwIHktyiCs7H7ae043uOUazpk830x5WGWowVySuu5ELcPkrav2urB6I/kdTjIjLm3WZEZFFE1uKualgjIwM9ZDCGJzN4ZSFTj3yYBJnlwnPOn0UkTzQ7ziNniZ03fZAJjWfi9W+EWpS3SBprpCu6JOhV99bC7K73W4r1eOgLMC0uLFC42Mm0g5wgfDG9430U2jxm/swgklU1xbvFAf8N+DvbJDLQ3WitGccXBG+NfdvhW+iUXOJrnJ90Lm++nfEgza7WSkJSKhcybmL38HOOmY3lzrUyluQ8D6bl7Lzpg0xqnCMUa3fLRauKXI7zjPrL2Xnzh5jSWMZft8+1es2vFKr6n8BNNshi8DLebhmessuDt8YuLhdtRfeyc+aMo1NXD/d3k844munct7vkow85ZjY2Ojxg1TTPcO54gVexpt/IzqHnmZmaCKksbzlBFDvTqL9EYnIa59NuYPvQc7aNmT9LTPcv294rIl9ijak2DGvD01/LtMaQVxLcm7rkxo9ZOYWe/1pQz7sW3p5dOcvFdTnbD99Kh+STcOH/C7UodDdYhR/jgjgbXS/xhz9EkkxT/UJoU284MS3JSsQf+hApTHH+RXv+r/kzg3jPsu12YBy4zxZpDAAkjDbSFVWEy+0O6nkLt+yiJnonBS1POjb1xkSHdfMGc3a1VsTlor30F9g1d4au5oshlWVpNppV5jx//svZefQuq5DQudAqVvdQPf2kOyaNy2rsuOZueskk6pw9+dT8URBfUdWPebdf8cYkOHu+GuZkz7QwnBCaNfaJHR+g1NPhCGOhLzz9dcxoNHnFzkqBcDnlN30cjwqtLzwWUjkW+2qZ1Wjyy5wXEXw5Lreb5sJ72DV9IqTLnCkTTfTGOteDaTnuqCiaCu4mbbabudmZgPfvj4L4Wz/3GQLA5PgI+fQzH6I14123fdTywnk1tA+2lYgfa6QzBLOrtZJXUmW5brZ/N6Sum/GjDXS6C3FHRYVMhrVQesuv4Ral6dl/vvLBNqAeD4XzbUw6MC3JShz4yJco+R9niYmNC3jfKyoIEblGRD4HZIvIZ5dtvw84++4MY7oaLQNZbF5o3vgSklI5n3U7e0aed2Q5yOyZVkYSykIthl/M7n6QAu3l/Cuhq72cNdPCcKIzPb58UVC+nbNxV1He9kRIUr/0d7d6k2Q6f8a1RFx8om0uzKv1GgMkAVFA8rJtDHivLdIYGG2zcsBkhHDNOPO6h4mTeS4++5WQyeCLmakJ8jx9zKWHxwrn7ls+zDDJzL8RmnGcmZogP4zGa4mF/f+dPAaofuk7QT93b2NwkmSGCysqCG/E9BeBo6r6xWXbo96qcgYbmO+9yIK6KNiyO2QyVO47Rn1UFTl133KUsbqr6TwuUWJyne2fvkRcfCK1efeyZ/yVkKSz7mystsYrRLPR9bLnpgcZIA3Pia8H/dyTncFJkhku+DMvmRKRvxCRp0Xk+aXNdskilLiRBrpdebasJ66FoW0PUu5poe7NF0Mqx3KGvbOrtJLQKc+1UnjLrxMti9Q/8w9BP/dwm7VcmVEaPuMFEB0TS33Bfeyd/HnQjdUy2MC4xgclSWY44I+CeBy4CJRjJe5rAY7bKFNEkz7dxkB8WajFYNftv8SUxjL68r+EWpS3mOu5iEeFgi3ODmBaTnHlHs7FHqS89TtBX1Of763Ho0J+efiM1xIlt/x6SIzV8RMt9EQXOTotSTDxZxQyVfUxYN677PRLwFGb5YpIFhcWyF/sZjYl9EbFpJR0qjNvY/fwc4wOD4RaHABihhvocWUTn5gcalHWxMLBj5LHAOd+Gtw19eiRRnoly7EZSVejcMuOkBirs2bbGU0IDxfXYOCPgpj3/u0WkbtF5ABg5l820NvRSKzM48pyho9/xvW/ToLMUvN08JdHfJE21UJ/mPinL2f3jR+gn3TkxL8G9bwpU20MxBYH9ZyBZMlYfe6F4OQHs5wg+plPs68OfLjhj4L4IxFJBT4H/BbwFeA3bZUqQhlstQxkiQXOMCpW7jvGhejdFDd8M+TV5jyLixQsdDCdGn43b3RMLI1F97Nn6g26WmqDck71eMhb6GAquSwo57ODvTd/kB6yiDoRnGWm7pYaXKJE54aX15edXKkmtRuoUtVRVa1W1RtV9SpVvWJFOcPameqpAyC3zDlrxjMHf5lC7eXci/8RUjl62huIlzkkK3gp0ANJ2e2fwIPQ9qO/Dsr5hvq7rDKjmc6Yja6HqOgYmrd8kN2zp4NSxnWk3Urjklbk3DxfweZKNakXgXuDJEvEowMNTGocmXnOWRbYe8uH6CUT9/HQRLYuMeCtEZxUGJ43b15xJWeSr2Nnz1NBSQPe32I97BLywsMleCV23PUI0xpD/0/+xvZzzXhf0PLC0KhvF/4sMb0qIn8nIu8WkYNLm+2SRSAJ4y10RxU6yoMiKjqGpvIH2TN7itaakyGTY9qBs6u1knjDp0hhinPft9+mM+Yty5pZGr7jBZCWlce5zNvZO/QMIwM9tp7LPdxEP+kkpaTbep5wwp8n0bXALi6tSf2/7RQqUsmcbXOkB8X2uz7JjEbTE4S3uJXQwUbHza7WyvZDN1MbtZ2iuq/bbtNZ7K9nTt1BTxlvB9m3fJo4mafmB/amgEuebKEvxvjfLMefinI3+thMwaAAMzszRZ6nj/nU0Lu4Xk56dj5nM25jz8APQ5afKX68hZ6oAkfNrtbDxMGHKdIezj5vT3rmJWJHm+l254dNkr7VKN95mOrY/Wxp/nfm52ZtO0/ufAeTSWW29R+O+FMwKFdEHhORH3q/7xSRj9svWmTR03IRtyhROc5848u6+dOWy+v3/k9Izp8528FofPjOHpbYd+tHrJoHJ/7J1vOkT7cxFOfMut3rYeHwr5HLIGd+bE9J0tGhfqsOfEb4ecnZiT+vY18DngEKvN/rgM/YJVCkMtRmGRVTHOpBsWX31ZyJO0xVyzeDXhJyfm6WPE8vsw6cXa2VqOgYWio/xK65czScecWWc3gWFylY7HJEwGWg2Hvj+2hzFZL25j/akh+st8VyMY8LkzxfwcIfBZGlqt8GPACqugCELsH9JmW21/lG2KjrPksmo5wJgpF1OT1tdUSJB7dDAgg3yo67f4NJjWPkJ39pS/+9HQ2OCrgMBC63m97dv0rlYiPVP/vPgPc/1lEDQHqJyeK6HH8UxKSIZOKtQy0iR4FRW6WKQFzDTQyRQmpGdqhFWZGdR++gNmobRTVfYWF+LmjnHWqz3u5SCjbH211qehbn8v8bB8aep6OhOuD9D7Ras1GnBFwGir13/Qp9ZOB6JfDLnAt99SyqhEXlvWDij4L4LPA9oEJEXgG+AfyGrVJFIIkTLfRGO9uDQlwupo58ikLt5fSzwUvFPN1jZZfP2UT+6ZX3/g4LRNH19JcC3vdUtxWt7eTZ6HqIjUugqeqj7Jo7y8UTzwW07+jRJrpduSHPouw0/PFiehO4Hsvd9VeBXap61m7BIo2cuQ7GE8tCLcYV2Xfzg7S6ikh/8x+CVitChpoYI4H0rPygnC8YZBWUcjrrbvYPPh3wlNZODLgMFHvu/TSjJDL9fGCX51Kn2hgM47xVduGPF1Mc8CngD7HSfT/i3WcIEBNjw2QzzGK68+vgutxu+nY/TMViE+deeioo50wYb6HXYQGEgaD4Pb+HC6X5e38W0H6dGHAZKBKT07hQ/CAHpl4JWOCmejxWnq9NZNQPFP78D/oGVqDc3wJ/B+wE/s1OoSKNnmavB0VeeOQZ2nf3r9JDFjEv/0VQZhGZcx2Mxm8el80lCsq2cSrtVvb1PsVQX2fA+nVqwGWg2H7vbzGlsfT/8E8D0l9/dysJMouEcd4qu/BHQWxT1Y+r6gve7WHgik8yESkWkRdEpEZEzovIp737H/B+94jIoVXa3yEitSLSICKf9/+Swo8RrwdFWlF4GMhiYuNo3fXrbF+osX0WMTszRa6nn/m0zfl2l3Pn54llntrv/nlA+lsKuJxLc/5sdL2kZ+dzNv+9HBj9Ca21pzfcX583z1dCfnjcf8HEHwVxyuu5BICIXA3448C9AHxOVXdgFRh6RER2AtXA/cBLKzX0ZpH9e+BOrBnLg962m5L5XssIG06Vvw7c+0m6ySb25T+zdRaxFEAYnb053+5Ktx/kdPJ17O34FsP93Rvu7+3xcmbAZaDYev8XmCWG/u//wYb7muyyjPrZZZv2EbNu/FEQV2Ml7GsRkRbg58D1InJORFY0Vqtqt9fAjaqOAzVAoarWqOqVkuIfARpUtUlV54BvAff5IWtYEj3SRA/ZYVX5KyY2jvbdn2DbQi1nX7SvUtpwu5V0LnmTuLj6IvPu/0U8s9Q+8Ycb7svpAZeBIiOnkDOF7+fg2PO01JzYUF860MCMRpNTuHlnXevFHwVxB1Y96uu9WzlwF3AP8B5/TiIiZcABwN+k7oVA+7LvHd59vvp+WEROiMiJ/v7Q5AnaKClTrfTHOtvF1RcH7n2ELskl4dU/t20WMdO7+VMwl+64ipNpt7K/+9v0d7VsqK+3Ai438XgtseP+LzBFHEM/+OKG+okbb6HbXYDL7Q6QZJsHf9xcW4ExIBXIXNpUtdX726qISBLwBPAZVR3zUy7xJcoK8n1ZVQ+p6qHsbOcGma1G3kInU0nhZ1SMjomlc88jVC3Uc8amspAy3MQwyaRm5trSv1MouO/3ceOh6cmNPexcw00Mk0JqelaAJHMuaVl5VBd/kIMTL9F47rX19zPTwUicz/fPiMcfN9c/BM4Cf8Ma032LSDSWcnhcVZ9cg1wdwHKn5CKgaw3tw4bRwV5SmEQzwnN6u/+eX6ND8kh55U9tSWGdONFKX9Tmv3kLt+zizcy7OdD/3Q2VJU2YbKcvavPEi1yJHff/LmMkMPbD9SlWz+Ii+Ys9zCaH3wtaMPBniel9QIWq3rCWdN8iIsBjQI2qPrpGuY4DVSJSLiIxwAework3Hb2t1hp7bHZ4ZpGMjoml9/Bvs8XTwpvfD3yG0qzZDsYSN5+Lqy/K7v99FBcd//n76+4jc7aT8YTICfhKzcjmfNlDHJh6lQuv/WjN7fu7W4iVeSRjc3rJbRR/FEQ1kLaOvo8BHwFuEpHT3u0uEflFEekArgF+ICLPAIhIgYg8DW8lBPwkVhbZGuDbqnp+HTI4nrFua804rSh8jbAH7/gYdVFbKTn9aEAzvU5PjpPHAAub2GVzOblFFZzK+29cNfxDmi8cX3P7udkZcnSAhZTIehve/8AX6CMD90/+3zXbwgbarNlaQt7m9vpaL/4oiD/FcnV9RkS+t7RdqZGqvqyqoqp7VXW/d3taVZ9S1SJVjVXVXFW93Xt8l6retaz906q6VVUrVPWP13+Jzma+30qzkFcavj7Y4nIxf9MXyWWQU/8RmOAlgB5vXeXonM3p4uqL7Q98kUlJYPx7aw/96W2rxS2KOysyFOoS8YnJtOz7TbYt1PLmj766prZT3jxfGcXh+4JmJ/4oiK8DfwZ8ibdtEPbkKY5A3COt9JNOfGJyqEXZELuuvYvTCdewu+mxgPjzA4x0WMtvqYXhqzzXSlpWHheqfo29Myc4+8La3IeHOqzZaFJ+5L0NX/WeT9DkKiP/+J8xOzPld7uFwSYW1EVOUeS8hKwFfxTEgKr+jTeK+qdLm+2SRQhJU+0MRBdc+cAwIP09f0wCM9R++38GpL/Z3gYAciMsgOnge3+bDskn5We/v6a06jPe8coujhyFuoQ7KorJ6/8XBdrLqSf88qEBIGaslV5XNtExsTZKF774oyBOisifisg1InJwabNdsggha66TiU1iVCzdcRUnMt/DVX1PBiSRmoy2MkwKKWmZAZAufIiJjaPv6Bco87Rz8qm/9rudDjUzpbGbMourP+y5/n7Oxl3Fjvp/YmSgx682KdMdDMVsfi+59eKPgjiAlSrjT1ijm6thdWamJshhiIW0slCLEjCqPvAlpiWO8f/87IaD5xIm2umPyguQZOHFgVs/xIWYPVRd+BvGRgb9ahM30UavO29TZnH1l+R7/4xEnab2//4/fh2fs9DFVFJkeMmtB38C5W70sV3RzdVwZXpbLQ+K6E1kVMzIKaRmx6fYPXuaU89srKhQxlwn4/HhF2EeCMTlIubuL5Gq49Q8/tt+tTEBX1C+8zAn8t7H4cH/ov7UiuneABgd6ieVSTS9LDjChSH+BMrlishjIvJD7/edIvJx+0Xb/Ax3LhkVN5eB7Kr7P0uju5yC1/+IqYn1VaddmJ8j19PPXErkvt1V7nsXJ7Lv51DfEzSceXnVYz2Li+SZgC8Adn3wTxmSVPQHn8OzuLjicX2tVhbl2JzIM+r7iz9z0a9hxSMsWVLrgM/YJVAkMdNnGRVzSjaXUTEqOob52/6cPAY4863fX1cffR1NRImHqMzNM7taD9s/9OcMSyqe//rNVR92Az1txJmALwCSUzNoOfh5ti7UcWIVG86oN4trWqFRECuxooIQkSjvxyxV/TbggbeC2Fb+n2rwGxlqZkLjN1UpzSW2X30bx1Nv46r2b6wrZ/9QhzeAKXdzza7WSmp61lsPu+NP/tWKxw20WS7B8RE+Xktcdc+vciF6N1XVj65osF6KQcoN4xgku1ltBvGG9++kiGTiTZbnrQ2xvnUDwyXETbTRG5W/aY2K5Q8+yrTEMvWdT6z69uuLyR5rdpVZHB5V9uzkqnt+lfMxe9h+/i8Z7O3weczEUsBXGEfkBxJxuYj/hUdJ0inq/+1TPo9xj7QyQBoJSalBli58WO3JtJRR9bNYeZAqROQVrBKkv2G3YJFAxmwno3Gb1wiblVdM7b7fZcf8eY4/sbbYSs9QC3PqJrvALJmIy0XSL/418TpDyzd933qLA1bAV26JWS5ZonzX1ZwofojDo8/4DDpMmmqnf5PEINnFagoiW0Q+C9wAPAX8OfBD4F+AW+wXbXOzuLBArqeX2U1uhD183yOciz3IrvOP0tvR6He72PFWel25uKOirnxwBFC64ypOlv0KV40/z6lnv/mO32PGWukzAV/v4OCH/5hWVxE5P/08E2PDl/y2mWKQ7GI1BeEGkoBkIBGI8u5L8O4zbIC+zkZiZBFXmKb59hdxucj8wD/iwkP345/wOzYiebqT4djIdtm8nEMf+gMa3VsofvULjA72XvJb8nQHgzHmbfhyYuMSmL7jr8jRAc5/43Nv7X87BsnMUFdjNQXRrap/oKpf9LUFTcJNymC7ZYRNzNv8RsWC8u2c3foI+6df4+T3/9mvNrmL3Uwnmbe75UTHxMJ9f0+qjlP3jU9e8psJ+FqZ7Udu5Y2c93L1wBNvpQTfjDFIduCPDcJgA1NeI2xGhOTNOfz+L1ATvYttJ79IT1v9qseODvVbRZTSjE//5VTsvZaTxR/l8OiznH7uWwCMDg+QxgS6iSLyA82e//6/6ZRc0p75FOOjQ2/FICVHYGLDtbCagrg5aFJEIIuDzcypm9yi8CwUtFbcUVGkfvBfcaEMfvOXVg9g8rpsxoRpESW7OfiRP6HJVUbxz36bgZ72ZQFfZrxWIjE5jfE7/55cTx81X33k7Rik0h0hlszZrKggVHUomIJEGrHjLRFnhC0o3875/V9g19xZ3vj3P1zxuLEu6+ZNKzQurr6IiY3D9cBjJOoUnV/7GGOd1nJJSoFxcV2N7Udu5Y3ij3Fk5GnyG77FuMaTtslrnW+UzemAHwakRKgR9vB9n+RU4rs4WP+3NFW/7vOYuYEmAHJKjIJYibIdhziz87fYN3OcrJP/B4DcUqMgrsSh//4l6qOqKPV00BtVsGljkAKFGZ0QoB4PuQtdEWmEFZeL0oe+zJgkEfXkx97hegjgHmlmiBSSUtJDIGH4cOSB3+Z0/FFKPR0MkmrGyw+iY2KJe/9jTGksIwlloRbH8RgFEQJGBntJlmk0PTJd7DJyCum99R8oXOzi4ld++R2urwmT7fRHbb70I4FGXC6KP/oYA6TRG2M8mPyluGofve/7L0o+YApjXgmjIEJAf7vlQRGbHbkudruO3c3x8k9waOwnvPHEo5f8ljHXzbgJYPKLzNwi5j76LGkffCzUooQV5buuJqcwMl/Q1oJRECFgvMdaY0/ZZGm+18qRj/whZ+MOs7/6SzSceQWA+blZcjz9zG/yCPNAUlC2jYIyY38wBB6jIELA/GAzANnFke2D7XK7Kf74vzEiKST850cZ7u+mr6OBKPHgzjRvdwZDqDEKIgTIaDujJEZcrWVfpGfnM3rvV8n0DNP1Lw8w0HIegESTttpgCDlGQYSAuMkO+t2RWWvZF1sPXs+5q/6IXXPnyH3p9wDIKjFLJgZDqDEKIgSkzXYzHme8dJZz6N5f4+f5HyaPfuY0iuz8slCLZDBEPEZBBBn1eMhZ7GU2AmMgrsSRj/81pxKO0RizDZfbHWpxDIaIJ3LyPDiEof4uMmUO0oyXzuW4o6LY/1vfx+NnSnCDwWAvZgYRZAY7rTxDcdnGS8cX4nJFVH4qg8HJGAURZMa7LQWRmh+5QXIGgyE8MAoiyMwNtgKQXWwS0RkMBmdjm4IQkWIReUFEakTkvIh82rs/Q0R+LCL13r8+M4yJSIuInBOR0yJywi45g41rtI1hkk1iNYPB4HjsnEEsAJ9T1R3AUeAREdkJfB54TlWrgOe831fiRlXdr6qHbJQzqMRPdjLgNjnoDQaD87FNQahqt6q+6f08DtQAhcB9wNe9h30d+AW7ZHAiaXNdjMeZ4vIGg8H5BMUGISJlwAHgdSBXVbvBUiJAzgrNFHhWRE6KyMOr9P2wiJwQkRP9/f2BFTzAWDEQfcwlmxgIg8HgfGxXECKSBDwBfEZVx9bQ9JiqHgTuxFqeus7XQar6ZVU9pKqHsrOzAyCxfQz2dRAn80i6iYEwGAzOx1YFISLRWMrhcVV90ru7V0Tyvb/nA32+2qpql/dvH/AUcMROWYPBQEc9AHHZZaEVxGAwGPzATi8mAR4DalR1eUWY7wEPeT8/BHzXR9tEEUle+gzcBlTbJWuwmOi16kCkRngdCIPBEB7YOYM4BnwEuMnrqnpaRO4CvgTcKiL1wK3e74hIgYg87W2bC7wsImeAN4AfqOqPbJQ1KMwPWHUgciK8DoTBYAgPbMtpoKovA7LCzzf7OL4LuMv7uQnYZ5dsocI12s4QKWQkpYZaFIPBYLgiJpI6iMRPdTIQZepAGAyG8MAoiCCSPtvNhKkDYTAYwgSjIIKEZ3GRXE8fsyYGwmAwhAlGQQSJwd52YmQBV3ppqEUxGAwGvzAKIkgMvhUDYepAGAyG8MAoiCCxFAORlmfqQBgMhvDAKIggMT/cAUBWUUWIJTEYDAb/MAoiSLjGOhghicTktFCLYjAYDH5hFESQiJ3sYtDt7GSCBoPBsByjIIJEymwP47EmSM5gMIQPRkEEiSxPP7MJJkjOYDCED0ZBBIGJsWFSmERTikItisFgMPiNURBBYLDTcnGNyjBR1AaDIXwwCiIIjPZYab6TTKEgg8EQRhgFEQSmB1oBSCswQXIGgyF8MAoiCHhGO1hQF1l5pha1wWAIH4yCCAJR450MSAZR0TGhFsVgMBj8xiiIIJAw3c1wdE6oxTAYDIY1YRREEEib72XSFAoyGAxhhlEQNuNZXCTbM8h8UmGoRTEYDIY1YRSEzQz1d1qFgtKMgjAYDOGFURA2M9RlBcnFZppKcgaDIbwwCsJmJnqtILnkXFNJzmAwhBdGQdjM3FAbAFmFplCQwWAIL4yCsJvRwgHcXwAAB55JREFUTqY0lpS0zFBLYjAYDGvCKAibiZnsot+dg7jMUBsMhvDCPLVsJnm2h9GY3FCLYTAYDGvGKAibyVjoY8YUCjIYDGGIURA2MjM9SSajLCabGAiDwRB+GAVhIwOdlotrVLopFGQwGMIP2xSEiBSLyAsiUiMi50Xk0979GSLyYxGp9/5NX6H9HSJSKyINIvJ5u+S0k5EeK0guPtsEyRkMhvDDzhnEAvA5Vd0BHAUeEZGdwOeB51S1CnjO+/0SRMQN/D1wJ7ATeNDbNqyY6rcKBaXnm0JBBoMh/Iiyq2NV7Qa6vZ/HRaQGKATuA27wHvZ14EXgdy5rfgRoUNUmABH5lrfdBTtkrf/Dq4jW2YD3u8UzDkBWgYmiNhgM4YdtCmI5IlIGHABeB3K9ygNV7RYRX4USCoH2Zd87gKtX6Pth4GGAkpL1VWwbTSzD5ZlbV9vVGALqM3ZwTVxCwPs2GAwGu7FdQYhIEvAE8BlVHRMRv5r52Ke+DlTVLwNfBjh06JDPY67Eoc8+sZ5mBoPBsKmx1YtJRKKxlMPjqvqkd3eviOR7f88H+nw07QCWu/4UAV12ymowGAyGS7HTi0mAx4AaVX102U/fAx7yfn4I+K6P5seBKhEpF5EY4APedgaDwWAIEnbOII4BHwFuEpHT3u0u4EvArSJSD9zq/Y6IFIjI0wCqugB8EngGqAG+rarnbZTVYDAYDJdhpxfTy/i2JQDc7OP4LuCuZd+fBp62RzqDwWAwXAkTSW0wGAwGnxgFYTAYDAafGAVhMBgMBp8YBWEwGAwGn4jqumLLHImI9AOt62yeBQwEUJxww1y/uX5z/ZFJqapm+/phUymIjSAiJ1T1UKjlCBXm+s31m+uP3OtfCbPEZDAYDAafGAVhMBgMBp8YBfE2Xw61ACHGXH9kY67f8A6MDcJgMBgMPjEzCIPBYDD4xCgIg8FgMPgk4hWEiNwhIrUi0iAi76iPvdkQkWIReUFEakTkvIh82rs/Q0R+LCL13r/poZbVTkTELSKnROT73u8Rc/0ikiYi3xGRi97/B9dE2PX/pvf/frWI/LuIxEXS9a+FiFYQIuIG/h64E9gJPCgiO0Mrle0sAJ9T1R3AUeAR7zV/HnhOVauA57zfNzOfxkolv0QkXf9fAz9S1e3APqxxiIjrF5FC4FPAIVXdDbix6s1ExPWvlYhWEMARoEFVm1R1DvgWcF+IZbIVVe1W1Te9n8exHg6FWNf9de9hXwd+ITQS2o+IFAF3A19Ztjsirl9EUoDrsIp5oapzqjpChFy/lyggXkSigASsapWRdP1+E+kKohBoX/a9w7svIhCRMuAA8DqQq6rdYCkRICd0ktnOXwG/DXiW7YuU698C9ANf9S6xfUVEEomQ61fVTuB/A21ANzCqqs8SIde/ViJdQfgqaBQRfr8ikoRVL/wzqjoWanmChYjcA/Sp6slQyxIiooCDwD+q6gFgkghaTvHaFu4DyoECIFFEPhxaqZxLpCuIDqB42fcirOnmpkZEorGUw+Oq+qR3d6+I5Ht/zwf6QiWfzRwD7hWRFqwlxZtE5JtEzvV3AB2q+rr3+3ewFEakXP8tQLOq9qvqPPAkcC2Rc/1rItIVxHGgSkTKRSQGy1j1vRDLZCsiIljrzzWq+uiyn74HPOT9/BDw3WDLFgxU9XdVtUhVy7D+vZ9X1Q8TOdffA7SLyDbvrpuBC0TI9WMtLR0VkQTvvXAzlh0uUq5/TUR8JLWI3IW1Ju0G/lVV/zjEItmKiLwL+BlwjrfX4H8Pyw7xbaAE6yZ6QFWHQiJkkBCRG4DfUtV7RCSTCLl+EdmPZaCPAZqAj2G9LEbK9X8ReD+WR98p4Jfh/2/vfl5sCuM4jr8/UWZhJ0kk5dcCpUhNKGwtlNha+AusbNmSnbIUZW8hjR8Lw1KIhoSNlJKyU5MGX4tz5Lod15jmMDXvV926PT2d89zF6dM9957Pw3IWyef/G4s+ICRJ3Rb7LSZJ0m8YEJKkTgaEJKmTASFJ6mRASJI6GRDSkCQrkjxpX++TvGvff0pysadznkxyvH0/mWRXx5ztSS73cX6py9L/vQBpoamqj8AOgCRngE9Vdb6v87WlcSdonmgeta6pJGuTrKuqt32tR/rBbxDSLCXZP7B/xJkkV5LcTvImyZEk55JMJbnZ1pmQZGeSe0keJbn1o85hyEHgcVV9GRg7luRBkldJ9g2MX6d5AlzqnQEhzd0Gmtrww8BV4G5VbQemgUNtSFwAjlbVTuAS0PWk/h5guDxwaVXtBk4CpwfGHwL7kP4BbzFJczdRVTNJpmiqWm6241PAemALsA2409T+sISmYnrYan7dvAiaEjlogmP9wPgHmhZSqXcGhDR3nwGq6luSmfrZW/ON5toK8Lyqxv9wnGlgrOvYwFd+vU7H2vlS77zFJPXnJbAyyTg0NetJtnbMewFsnOUxNwPP5ml90kgGhNSTdhvbo8DZJE+BJzR7DwyboNkGdDYOADfmZ4XSaLa5SgtAkmvAqap6PWLOMuAesHfoH09SLwwIaQFoN/BZVVX3R8zZBKypqsl/tjAtagaEJKmTv0FIkjoZEJKkTgaEJKmTASFJ6mRASJI6fQdmj6Jx3WxyTwAAAABJRU5ErkJggg==\n", |
|
|
"text/plain": [ |
|
|
"<Figure size 432x288 with 1 Axes>" |
|
|
] |
|
|
}, |
|
|
"metadata": { |
|
|
"needs_background": "light" |
|
|
}, |
|
|
"output_type": "display_data" |
|
|
} |
|
|
], |
|
|
"source": [ |
|
|
"# Code snippet by S. Rouchier available on:\n", |
|
|
"# https://buildingenergygeeks.org/ssm_rc.html\n", |
|
|
"\n", |
|
|
"# Optimisation\n", |
|
|
"import scipy.optimize as so\n", |
|
|
"\n", |
|
|
"# The curve_fit function will find the optimal value for the parameters R and C\n", |
|
|
"p_opt, p_cov = so.curve_fit(f=Ti_calc, xdata=inputs, ydata=Ti, p0=(1/HLC, C_env))\n", |
|
|
"R_opt = p_opt[0]\n", |
|
|
"H_opt = 1/p_opt[0]\n", |
|
|
"C_opt = p_opt[1]\n", |
|
|
"print('R: %.2f K/(W m^2)' % R_opt)\n", |
|
|
"print('H: %.2f W/(K m^2)' % H_opt)\n", |
|
|
"print('C: %.2f Wh/(K m^2)' % C_opt)\n", |
|
|
"\n", |
|
|
"#%% Plot the results to see if the calculated temperature fits well with measurements\n", |
|
|
"plt.figure()\n", |
|
|
"plt.plot(Ti)\n", |
|
|
"plt.plot(Ti_calc(inputs, R_opt, C_opt))\n", |
|
|
"plt.xlabel('Time (h)')\n", |
|
|
"plt.ylabel('Temperature (C)')\n", |
|
|
"plt.show()" |
|
|
] |
|
|
}, |
|
|
{ |
|
|
"cell_type": "code", |
|
|
"execution_count": 668, |
|
|
"metadata": {}, |
|
|
"outputs": [], |
|
|
"source": [ |
|
|
"# Try to rerun the above code after increasing or reducing the internal heat gains by multiplications or additions to see changes in the parameter estimation:\n", |
|
|
"P = (q_int_t + q_sol_t) * 1.2 + 1.1 # Combined heat gains" |
|
|
] |
|
|
}, |
|
|
{ |
|
|
"cell_type": "markdown", |
|
|
"metadata": {}, |
|
|
"source": [ |
|
|
"### Building energy calculation methododology and norms: \n", |
|
|
"\n", |
|
|
"To view or print the standards, you have to use the following link via VPN or the campus network:\n", |
|
|
"https://innsida.ntnu.no/wiki/-/wiki/English/standards\n", |
|
|
"\n", |
|
|
"* **NS 3031:2014** \"Energy performance of buildings — Calculation of energy needs and energy supply\" (Available in Norwegian only)\n", |
|
|
"\n", |
|
|
"* **SN-NSPEK 3031:2020** \"Energy performance of buildings — Calculation of energy needs and energy supply\" (Available in Norwegian only)\n", |
|
|
"\n", |
|
|
"* **NS-EN ISO 52000-1:2017** \"Energy performance of buildings - Overarching EPB assessment - Part 1: General framework and procedures (ISO 52000-1:2017)\"\n", |
|
|
"\n", |
|
|
"* **NS-EN ISO 52000-1:2017** \"Energy performance of buildings - Energy needs for heating and cooling, internal temperatures and sensible and latent heat loads - Part 1: Calculation procedures (ISO 52016-1:2017)\"\n", |
|
|
"\n", |
|
|
"* **NS-EN ISO 13790:2008** \"Energy performance of buildings — Calculation of energy use for space heating and cooling (ISO 13790:2008)\"\n", |
|
|
"\n", |
|
|
"## Learn more about inverse modeling\n", |
|
|
"\n", |
|
|
"More tutorials on building energy simulation and statistical learning by S. Rouchier:\n", |
|
|
"https://buildingenergygeeks.org/\n", |
|
|
"\n", |
|
|
"A great introduction to RC-networks and inverse modeling techniques in the report:\n", |
|
|
"\"Thermal performance characterisation using time series data – statistical guidelines\"\n", |
|
|
"https://www.teknologisk.dk/_/media/76411_EBC_Annex_58_Final_Report_ST3b.pdf\n" |
|
|
] |
|
|
}, |
|
|
{ |
|
|
"cell_type": "code", |
|
|
"execution_count": null, |
|
|
"metadata": {}, |
|
|
"outputs": [], |
|
|
"source": [] |
|
|
} |
|
|
], |
|
|
"metadata": { |
|
|
"kernelspec": { |
|
|
"display_name": "Python 3", |
|
|
"language": "python", |
|
|
"name": "python3" |
|
|
}, |
|
|
"language_info": { |
|
|
"codemirror_mode": { |
|
|
"name": "ipython", |
|
|
"version": 3 |
|
|
}, |
|
|
"file_extension": ".py", |
|
|
"mimetype": "text/x-python", |
|
|
"name": "python", |
|
|
"nbconvert_exporter": "python", |
|
|
"pygments_lexer": "ipython3", |
|
|
"version": "3.7.7" |
|
|
} |
|
|
}, |
|
|
"nbformat": 4, |
|
|
"nbformat_minor": 4 |
|
|
} |
kristss revised this gist