{ "cells": [ { "cell_type": "markdown", "id": "35a53b91-f1eb-4c06-983f-3aaa32ac5a87", "metadata": { "tags": [] }, "source": [ "# Steady-state and step response\n", "\n", "[![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/cghiaus/dm4bem_book/HEAD?labpath=%2Ftutorials%2F02_2_2Toy.ipynb)\n", "\n", "[Steady-state](https://en.m.wikipedia.org/wiki/Steady_state) analysis focuses on the behavior of a system after transients have settled, specifically when the system reaches a stable equilibrium at a constant input.\n", "\n", "[Step response](https://en.m.wikipedia.org/wiki/Step_response) analyses system's dynamic response to a sudden change in input (i.e., to a step input). Step response is concerned with transient behavior, such as how the system responds immediately after the input change and how it eventually settles into a new steady-state.\n", "\n", "**Objectives:**\n", "- Find the steady-state values for thermal circuit and state-space.\n", "- Estimate the time step and the settling time from eigenvalue analysis.\n", "- Simulate the step response.\n", "- Compare the steady-state results with step response after settling time.\n", "\n", "![thermal_circuit](../figures/03_therm_circ.svg)\n", "> Figure 1. Thermal circuit for the cubic building." ] }, { "cell_type": "code", "execution_count": 1, "id": "764b77e2-3f93-4484-9842-6b328305d517", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "\n", "import dm4bem" ] }, { "cell_type": "markdown", "id": "55996710-77de-4f3d-a33a-1006e6911424", "metadata": {}, "source": [ "The following assumptions are done:\n", "- The indoor air temperature is controlled or not (i.e., the building is in free running).\n", "- The heat capacity of air and of glass is neglected or not.\n", "- The time step is calculated from the eigenvalues or it is imposed at a value designated as $\\Delta t$." ] }, { "cell_type": "code", "execution_count": 2, "id": "6f842ded-7a14-4a9b-85a2-d4e815262048", "metadata": {}, "outputs": [], "source": [ "controller = False\n", "neglect_air_glass_capacity = False\n", "imposed_time_step = False\n", "Δt = 498 # s, imposed time step" ] }, { "cell_type": "markdown", "id": "adbccd8a-d608-44ff-a796-e26d4432110f", "metadata": {}, "source": [ "The thermal circuit was described in the section on modelling (Figure 1). It is given by a _thermal circuit_ imported from the file [./toy_model/TC.csv](./toy_model/TC.csv). \n", "\n", "A thermal ciruit file contains the incidence matrix $A$, the diagonals of matrices $G$ and $C$, and the vectors $b$, $f$, and $y$ (see example below). If there is no value in a cell, it is considered `0`." ] }, { "cell_type": "code", "execution_count": 3, "id": "68d9fdda-aa62-4e46-9bbb-10cf18766f4d", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Matrices and vectors for thermal circuit from Figure 1\n" ] }, { "data": { "text/html": [ "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
 Aθ0θ1θ2θ3θ4θ5θ6θ7Gb
0q010.0000000000001125.000000To
1q1-11.000000000000630.0000000
2q20-1.000000100000630.0000000
3q300.000000-11000030.3750000
4q400.0000000-1100030.3750000
5q500.00000000-110044.7868240
6q600.00000000-1010360.0000000
7q700.000000000-11072.0000000
8q800.000000000001165.789474To
9q900.00000000010-1630.0000000
10q1000.0000000000109.000000To
11q1100.0000000000100.000000Ti_sp
12C018216000.000000023958000324001.089E+06nannan
13fΦo0.00000000Φi0QaΦanannan
14y00.000000000010nannan
\n" ], "text/plain": [ "" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "print('Matrices and vectors for thermal circuit from Figure 1') \n", "df = pd.read_csv('./toy_model/TC_0.csv')\n", "df.style.apply(lambda x: ['background-color: yellow'\n", " if x.name in df.index[-3:] or c in df.columns[-2:]\n", " else '' for c in df.columns], axis=1)" ] }, { "cell_type": "markdown", "id": "0c51ac0b-bc04-4914-b090-31dec69a74c3", "metadata": {}, "source": [ "In the thermal circuit [./toy_model/TC.csv](./toy_model/TC.csv), the controller is _off_, i.e., the gain of the proportional controller is 0, $K_p = G_{11} = 0$; the building is in _free-floating_. The controller may be turned _on_ by setting $K_p = G_{11} \\rightarrow \\infty$." ] }, { "cell_type": "code", "execution_count": 4, "id": "afe67ba5-5f98-4793-949a-63e9f4020726", "metadata": { "tags": [] }, "outputs": [], "source": [ "# MODEL\n", "# =====\n", "# Thermal circuit\n", "TC = dm4bem.file2TC('./toy_model/TC.csv', name='', auto_number=False)\n", "\n", "# by default TC['G']['q11'] = 0, i.e. Kp -> 0, no controller (free-floating)\n", "if controller:\n", " TC['G']['q11'] = 1e3 # Kp -> ∞, almost perfect controller\n", "\n", "if neglect_air_glass_capacity:\n", " TC['C']['θ6'] = TC['C']['θ7'] = 0\n", " # or\n", " TC['C'].update({'θ6': 0, 'θ7': 0})" ] }, { "cell_type": "markdown", "id": "e4524cc6-07cf-415c-b778-20ede67cdd9b", "metadata": {}, "source": [ "Thermal circuit is transformed in state-space representation." ] }, { "cell_type": "code", "execution_count": 5, "id": "b9754e8c-721e-430d-b5b6-374598cc78a9", "metadata": {}, "outputs": [], "source": [ "# State-space\n", "[As, Bs, Cs, Ds, us] = dm4bem.tc2ss(TC)" ] }, { "cell_type": "markdown", "id": "39c88179-2905-469d-b7f6-847006bdf671", "metadata": { "tags": [] }, "source": [ "## Steady-state analysis\n", "[Steady-state](https://en.m.wikipedia.org/wiki/Steady_state) means that the temperatures do not vary in time. In the system of differential-algebraic equations (DAE):\n", "\n", "$$\\left\\{\\begin{array}{ll}\n", "C \\dot{\\theta} = -(A^T G A) \\theta + A^T G b + f\\\\ \n", "q = G (-A \\theta + b)\n", "\\end{array}\\right.$$\n", "\n", "the term $\\dot \\theta = 0$, or in the state-space representation,\n", "\n", "$$\\left\\{\\begin{array}{rr}\n", "\\dot{\\theta}_s=A_s \\theta_s + B_s u\\\\ \n", "y = C_s \\theta_s + D_s u\n", "\\end{array}\\right.$$\n", "\n", "the term $\\dot \\theta_s = 0$.\n", "\n", "In [steady-state](https://en.m.wikipedia.org/wiki/Steady_state), the model can be checked if it is incorrect. Let's consider that:\n", "- the controller is not active, $K_p \\rightarrow 0$,\n", "- the outdoor temperature is $T_o = 10 \\, \\mathrm{^\\circ C}$,\n", "- the indoor temperature setpoint is $T_{i,sp} = 20 \\, \\mathrm{^\\circ C}$,\n", "- all flow rate sources are zero." ] }, { "cell_type": "code", "execution_count": 6, "id": "85631d1b-66cd-4df6-ac99-0b6cddc3cac2", "metadata": {}, "outputs": [], "source": [ "bss = np.zeros(12) # temperature sources b for steady state\n", "bss[[0, 8, 10]] = 10 # outdoor temperature\n", "bss[[11]] = 20 # indoor set-point temperature\n", "\n", "fss = np.zeros(8) # flow-rate sources f for steady state" ] }, { "cell_type": "markdown", "id": "aaa04025-f87d-43d6-9ef6-2cddb1984b8a", "metadata": { "tags": [] }, "source": [ "*Note*: Steady-state analysis is a test of [falsification (refutability)](https://en.m.wikipedia.org/wiki/Falsifiability) of the model, not a [verification and validation](https://en.m.wikipedia.org/wiki/Verification_and_validation). If the model does not pass the steady-state test, it means that it is wrong. If the model passes the steady-state test, it does not mean that it is correct. For this example, the values of the capacities in matrix $C$ or of the conductances in matrix $G$ can be wrong, even when the steady-state test is passed. " ] }, { "cell_type": "markdown", "id": "d342c6f7-9b07-4e51-99bd-b0372e226b38", "metadata": { "tags": [] }, "source": [ "### Steady-state from differential algebraic equations (DAE)\n", "\n", "The temperature values in [steady-state](https://en.m.wikipedia.org/wiki/Steady_state) are obtained from the system of DAE by considering that $C \\dot{\\theta} = 0$:\n", "\n", "$$\\theta_{ss} = (A^T G A)^{-1}(A^T G b + f)$$\n", "\n", "For the conditions mentioned above, in steady-state, all temperatures $\\theta_0 ... \\theta_7,$ including the indoor air temperature $\\theta_6,$ are equal to $T_o = 10 \\, \\mathrm{^\\circ C}$ (Figure 1)." ] }, { "cell_type": "code", "execution_count": 7, "id": "8928e41b-0afc-467f-942b-869939ff828b", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "θss = [10. 10. 10. 10. 10. 10. 10. 10.] °C\n" ] } ], "source": [ "A = TC['A']\n", "G = TC['G']\n", "diag_G = pd.DataFrame(np.diag(G), index=G.index, columns=G.index)\n", "\n", "θss = np.linalg.inv(A.T @ diag_G @ A) @ (A.T @ diag_G @ bss + fss)\n", "print(f'θss = {np.around(θss, 2)} °C')" ] }, { "cell_type": "markdown", "id": "ace90f7f-df8d-4217-85c2-7f3d93f0fa3e", "metadata": {}, "source": [ "If the sum auxiliary heat sources in the room, $\\dot Q_a = f_6$, is not zero, the temperature in the room $\\theta_6$ has the highest value and the temperature of the outdoor surface of the wall $\\theta_0$ is almost equal to the outdoor temperature, $T_o = b_0$." ] }, { "cell_type": "code", "execution_count": 8, "id": "8b9d04cc-07b3-44de-9d57-6d18c488ae0c", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "θssQ = [ 0.14 0.39 0.65 5.88 11.12 5.57 12.26 4.41] °C\n" ] } ], "source": [ "bss = np.zeros(12) # temperature sources b for steady state\n", "\n", "fss = np.zeros(8) # flow-rate sources f for steady state\n", "fss[[6]] = 1000\n", "\n", "θssQ = np.linalg.inv(A.T @ diag_G @ A) @ (A.T @ diag_G @ bss + fss)\n", "print(f'θssQ = {np.around(θssQ, 2)} °C')" ] }, { "cell_type": "markdown", "id": "eef1fe7c-a6c8-4e5c-91cb-a908aac9f122", "metadata": { "tags": [] }, "source": [ "### Steady-state from state-space representation\n", "The input vector $u$ of the state-space representation is obtained by stacking the vectors $b_T$ and $f_Q$:\n", "\n", "$$u = \\begin{bmatrix} b_T \\\\ f_Q\\end{bmatrix}$$\n", "\n", "where:\n", "- $b_T$ is a vector of the nonzero elements of vector $b$ of temperature sources. For the circuit presented in Figure 1, $b_T = [T_o, T_o, T_o, T_{i,sp}]^T$ corresponding to branches 0, 8, 10 and 11, where:\n", " - $T_o$ - outdoor temperature, °C;\n", " - $T_{i,sp}$ - set-point temperaure for the indoor air, °C.\n", "- $f_Q$ - vector of nonzero elements of vector $f$ of flow sources. For the circuit presented in Figure 1, $f_Q = [\\Phi_o, \\Phi_i, \\dot{Q}_a, \\Phi_a]^T$, corresponding to nodes 0, 4, 6, and 7, where:\n", " - $\\Phi_o$ - solar radiation absorbed by the outdoor surface of the wall, W;\n", " - $\\Phi_i$ - solar radiation absorbed by the indoor surface of the wall, W;\n", " - $\\dot{Q}_a$ - auxiliary heat gains (i.e., occupants, electrical devices, etc.), W;\n", " - $\\Phi_a$ - solar radiation absorbed by the glass, W.\n", "\n", "*Note*: Zero in vectors $b$ and $f$ indicates that there is no source on the branch or in the node, respectively. However, a source can have the value zero." ] }, { "cell_type": "code", "execution_count": 9, "id": "a5752ab5-b289-41f5-9268-4478190bdc3b", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "uss = [10 10 10 20 0 0 0 0]\n" ] } ], "source": [ "bT = np.array([10, 10, 10, 20]) # [To, To, To, Tisp]\n", "fQ = np.array([0, 0, 0, 0]) # [Φo, Φi, Qa, Φa]\n", "uss = np.hstack([bT, fQ]) # input vector for state space\n", "print(f'uss = {uss}')" ] }, { "cell_type": "markdown", "id": "88e43402-4429-4c5d-a0cb-4a9753372394", "metadata": { "tags": [] }, "source": [ "The steady-state value of the output of the state-space representation is obtained when $\\dot \\theta_{C} = 0$:\n", "\n", "$$y_{ss} = (-C_s A_s^{-1} B_s + D_s) u$$" ] }, { "cell_type": "code", "execution_count": 10, "id": "06543d9d-44f5-41ba-8a42-21b9f2160977", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "yss = 10.00 °C\n" ] } ], "source": [ "inv_As = pd.DataFrame(np.linalg.inv(As),\n", " columns=As.index, index=As.index)\n", "yss = (-Cs @ inv_As @ Bs + Ds) @ uss\n", "\n", "yss = float(yss.values[0])\n", "print(f'yss = {yss:.2f} °C')" ] }, { "cell_type": "markdown", "id": "3df80dbd-51ed-4ad8-b876-c3f55e6448b3", "metadata": {}, "source": [ "The error between the steady-state values obtained from the system of DAE, $\\theta_6$, and the output of the state-space representation, $y_{ss}$, \n", "\n", "$$\\varepsilon = \\left | \\theta_6 - y_{ss} \\right |$$\n", "\n", "is practically zero; the slight difference is due to [numerical errors](https://en.m.wikipedia.org/wiki/Numerical_error)." ] }, { "cell_type": "code", "execution_count": 11, "id": "774b0ca8-08cf-4f61-8735-e5d40718d6ca", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Error between DAE and state-space: 8.88e-15 °C\n" ] } ], "source": [ "print(f'Error between DAE and state-space: {abs(θss[6] - yss):.2e} °C')" ] }, { "cell_type": "markdown", "id": "92614b66-cd57-44b5-94d6-1641f0b2a9c3", "metadata": {}, "source": [ "Likewise, the steady-state of the output of the state-space representation for input $\\dot Q_a$ is obtained from\n", "\n", "$$y_{ss} = (-C_s A_s^{-1} B_s + D_s) u$$" ] }, { "cell_type": "code", "execution_count": 12, "id": "729cc0e7-818d-4033-818d-96715f8796ef", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "yssQ = 12.26 °C\n" ] } ], "source": [ "bT = np.array([0, 0, 0, 0]) # [To, To, To, Tisp]\n", "fQ = np.array([0, 0, 1000, 0]) # [Φo, Φi, Qa, Φa]\n", "uss = np.hstack([bT, fQ])\n", "\n", "inv_As = pd.DataFrame(np.linalg.inv(As),\n", " columns=As.index, index=As.index)\n", "yssQ = (-Cs @ inv_As @ Bs + Ds) @ uss\n", "\n", "yssQ = float(yssQ.values[0])\n", "print(f'yssQ = {yssQ:.2f} °C')" ] }, { "cell_type": "markdown", "id": "163c2431-a64b-4a90-9182-96d395a7fefb", "metadata": {}, "source": [ "The steady-state values obtained from DAE and state-space representation are pratically equal. " ] }, { "cell_type": "code", "execution_count": 13, "id": "49b56eb6-b09c-4c49-8fdd-861b6bd27662", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Error between DAE and state-space: 7.11e-15 °C\n" ] } ], "source": [ "print(f'Error between DAE and state-space: {abs(θssQ[6] - yssQ):.2e} °C')" ] }, { "cell_type": "markdown", "id": "747a5518-38c0-429a-9fbe-22154079afd1", "metadata": { "tags": [] }, "source": [ "## Eigenvalues analysis\n", "\n", "Let's consider the eigenvalues $\\lambda$ of the state matrix $A_s$." ] }, { "cell_type": "code", "execution_count": 14, "id": "4b3d65d2-56ad-4efb-887b-c0795efe12b5", "metadata": {}, "outputs": [], "source": [ "# Eigenvalues analysis\n", "λ = np.linalg.eig(As)[0] # eigenvalues of matrix As" ] }, { "cell_type": "markdown", "id": "9c0d162e-389e-4de3-b59f-76580adb72ea", "metadata": { "tags": [] }, "source": [ "### Time step\n", "\n", "The condition for [numerical stability](https://en.m.wikipedia.org/wiki/Euler_method#Numerical_stability) of [Euler explicit integration](https://en.m.wikipedia.org/wiki/Numerical_methods_for_ordinary_differential_equations#Euler_method) method is\n", "\n", "$$\\left | \\lambda_i \\Delta t + 1 \\right | < 1, \\forall \\lambda_i, $$\n", "\n", "i.e., in the [complex plane](https://en.m.wikipedia.org/wiki/Complex_plane), $\\lambda_i \\Delta t$ is inside a circle of radius 1 centered in $(-1 + 0j) \\in \\mathbb{C}$, where:\n", "- $\\lambda_i$ are the eigenvalues of matrix $A_s$,\n", "- $\\Delta t$ - time step,\n", "- $j = \\sqrt {-1}$ - [imaginary unit](https://en.m.wikipedia.org/wiki/Imaginary_unit).\n", "\n", "For real negative eigenvalues $\\left \\{ \\lambda_i \\in \\mathbb{R} \\;| \\; \\lambda_i <0 \\right \\}$, which is the case of thermal networks, the above condition becomes:\n", "\n", "$$- \\lambda_i \\Delta t - 1 < 1, \\forall \\lambda_i, $$\n", "\n", "or\n", "\n", "$$ 0 < \\Delta t < -\\frac{2}{\\min \\lambda_i} = \n", "2 \\min \\left ( -\\frac{1}{\\lambda_i} \\right ) = \n", "2 \\min T_i$$\n", "\n", "where $T_i$ are the [time constants](https://en.m.wikipedia.org/wiki/Time_constant), $T_i = - 1 / \\lambda_i$. Let's chose a time step $dt$ smaller than $\\Delta t_{max} = 2 \\min (-1 / \\lambda_i) $, floor rounded to:\n", "- 1 s, 10 s,\n", "- 60 s (1 min), 300 s (5 min), 600 s (10 min), 1800 s (30 min),\n", "- 36000 s (1 h), 7200 s (2 h), 14400 s (4 h), 21600 s (6 h)." ] }, { "cell_type": "code", "execution_count": 15, "id": "df463621-1d41-4c2b-bc4d-68199da29c47", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Δtmax = 498 s = 8.3 min\n", "dt = 300 s = 5.0 min\n" ] } ], "source": [ "# time step\n", "Δtmax = 2 * min(-1. / λ) # max time step for stability of Euler explicit\n", "dm4bem.print_rounded_time('Δtmax', Δtmax)\n", "\n", "if imposed_time_step:\n", " dt = Δt\n", "else:\n", " dt = dm4bem.round_time(Δtmax)\n", "dm4bem.print_rounded_time('dt', dt)" ] }, { "cell_type": "code", "execution_count": 16, "id": "d8b2e723-c306-4b10-85e2-09bc19e18d41", "metadata": { "tags": [] }, "outputs": [], "source": [ "if dt < 10:\n", " raise ValueError(\"Time step is too small. Stopping the script.\")" ] }, { "cell_type": "markdown", "id": "0d83dfe8-fe63-42c9-b913-47b1019d1a00", "metadata": {}, "source": [ "### Settling time and duration\n", "The [settling time](https://en.m.wikipedia.org/wiki/Step_response) is roughly 4 times the largest time constant." ] }, { "cell_type": "code", "execution_count": 17, "id": "89a7aa43-6fb4-46d8-bb9b-5fbd15a7102e", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "t_settle = 176132 s = 48.9 h\n" ] } ], "source": [ "# settling time\n", "t_settle = 4 * max(-1 / λ)\n", "dm4bem.print_rounded_time('t_settle', t_settle)" ] }, { "cell_type": "markdown", "id": "1a6e5158-dda8-4cee-b81f-9b48258bbae7", "metadata": {}, "source": [ "The duration of the simulation needs to be larger than the estimated [settling time](https://en.m.wikipedia.org/wiki/Settling_time). This requires a corresponding number of time steps in the time vector." ] }, { "cell_type": "code", "execution_count": 18, "id": "6062f401-b335-460a-bbba-61ddeb700662", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "duration = 176400 s = 49.0 h\n" ] } ], "source": [ "# duration: next multiple of 3600 s that is larger than t_settle\n", "duration = np.ceil(t_settle / 3600) * 3600\n", "dm4bem.print_rounded_time('duration', duration)" ] }, { "cell_type": "markdown", "id": "1dac581d-272f-4cdd-8843-eeb44dd6a58d", "metadata": { "tags": [] }, "source": [ "## Step response to outdoor temperature" ] }, { "cell_type": "markdown", "id": "e1cf346d-6283-4095-aa29-6dc378f54eb6", "metadata": { "tags": [] }, "source": [ "### Input vector\n", "In dynamic simulation, the inputs are [time series](https://en.m.wikipedia.org/wiki/Time_series), e.g., the oudoor temperature will have $n$ values $T_o = [T_{o(0)}, T_{o(1)}, ..., T_{o(n-1)}]$ at [discrete time](https://en.m.wikipedia.org/wiki/Discrete_time_and_continuous_time#Discrete_time) $t = [t_0, t_1, ... , t_{n-1}]$.\n", "\n", "The input vector $u$ of the state-space representation is obtained by stacking the vectors $b_T$ and $f_Q$ of the system of Differential Algebraic Equations:\n", "\n", "$$u = \\begin{bmatrix} b_T \\\\ f_Q\\end{bmatrix}$$\n", "\n", "where:\n", "- vector $b_T$ consists of the nonzero elements of vector $b$ of temperature sources; for the circuit presented in Figure 1, \n", "\n", "$$b = [\\begin{matrix}\n", "T_o &0 &0 &0 &0 &0 &0 &0 &T_o &0 &T_o &T_{i,sp} \n", "\\end{matrix}]^T$$\n", "\n", "and \n", "\n", "$$b_T = [T_o, T_o, T_o, T_{i,sp}]^T$$\n", "\n", "corresponding to branches 0, 8, 10 and 11; \n", "- vector $f_Q$ is the nonzero elements of vector $f$ of flow sources; for the circuit presented in Figure 1,\n", "\n", "$$f = [\\begin{matrix}\n", "\\Phi_o &0 &0 &0 &\\Phi_i &0 &\\dot{Q_a} &\\Phi_a \n", "\\end{matrix}]^T$$\n", "\n", "and\n", "\n", "$$f_Q = [\\Phi_o, \\Phi_i, \\dot{Q}_a, \\Phi_a]^T$$\n", "\n", "corresponding to nodes 0, 4, 6, and 7.\n", "\n", "For the thermal circuit shown in Figure 1, the [time series](https://en.m.wikipedia.org/wiki/Time_series) of the input vector, $u = [u_0, u_1, ... , u_{n-1}]^T$, is:\n", "\n", "$$u = \n", "\\begin{bmatrix}\n", "T_o\\\\ \n", "T_o\\\\ \n", "T_o\\\\ \n", "T_{i,sp}\\\\ \n", "\\Phi_o\\\\ \n", "\\Phi_i\\\\ \n", "\\dot{Q}_a\\\\ \n", "\\Phi_a\n", "\\end{bmatrix}\n", "= \\begin{bmatrix}\n", "T_{o(0)} & T_{o(1)}& ... & T_{o(n-1)}\\\\ \n", "T_{o(0)} & T_{o(1)}& ... & T_{o(n-1)}\\ \\\\ \n", "T_{o(0)} & T_{o(1)}& ... & T_{o(n-1)}\\ \\\\ \n", "T_{i,sp(0)} & T_{i,sp(1)}& ... & T_{i,sp(n-1)}\\ \\\\ \n", "\\Phi_{o(0)} & \\Phi_{o(1)} & ... & \\Phi_{o(n-1)}\\\\\n", "\\Phi_{i(0)} & \\Phi_{i(1)} & ... & \\Phi_{i(n-1)}\\\\ \n", "\\dot{Q}_{a(0)} & \\dot{Q}_{a(1)} & ... & \\dot{Q}_{a(n-1)}\\\\ \n", "\\Phi_{a(0)} & \\Phi_{a(1)} & ... & \\Phi_{a(n-1)}\n", "\\end{bmatrix}$$\n", "\n", "where:\n", "- $T_o = [T_{o(0)}, T_{o(1)}, ..., T_{o(n-1)}]$ is the [time series](https://en.m.wikipedia.org/wiki/Time_series) of the oudoor temperature at [discrete time](https://en.m.wikipedia.org/wiki/Discrete_time_and_continuous_time#Discrete_time) $t = [t_0, t_1, ... , t_{n-1}]$.\n", "\n", "- $T_{i, sp} = [T_{{i, sp}(0)}, T_{{i, sp}(1)}, ..., T_{{i, sp}(n-1)}]$ is the [time series](https://en.m.wikipedia.org/wiki/Time_series) of the setpoint indoor temperature at [discrete time](https://en.m.wikipedia.org/wiki/Discrete_time_and_continuous_time#Discrete_time) $t = [t_0, t_1, ... , t_{n-1}]$.\n", "\n", "- $\\Phi_o = [\\Phi_{o(0)}, \\Phi_{o(1)}, ..., \\Phi_{o(n-1)}]$ is the [time series](https://en.m.wikipedia.org/wiki/Time_series) of the solar radiation absorbed by the outdoor surface of the wall at [discrete time](https://en.m.wikipedia.org/wiki/Discrete_time_and_continuous_time#Discrete_time) $t = [t_0, t_1, ... , t_{n-1}]$.\n", "\n", "- $\\Phi_i = [\\Phi_{i(0)}, \\Phi_{i(1)}, ..., \\Phi_{i(n-1)}]$ is the [time series](https://en.m.wikipedia.org/wiki/Time_series) of the solar radiation absorbed by the indoor surface of the wall at [discrete time](https://en.m.wikipedia.org/wiki/Discrete_time_and_continuous_time#Discrete_time) $t = [t_0, t_1, ... , t_{n-1}]$.\n", "\n", "- $\\dot{Q}_a = [\\dot{Q}_{a(0)}, \\dot{Q}_{a(1)}, ..., \\dot{Q}_{a(n-1)}]$ is the [time series](https://en.m.wikipedia.org/wiki/Time_series) of the auxiliary heat gains (i.e., occupants, electrical devices, etc.) at [discrete time](https://en.m.wikipedia.org/wiki/Discrete_time_and_continuous_time#Discrete_time) $t = [t_0, t_1, ... , t_{n-1}]$.\n", "\n", "- $\\Phi_a = [\\Phi_{a(0)}, \\Phi_{a(1)}, ..., \\Phi_{a(n-1)}]$ is the [time series](https://en.m.wikipedia.org/wiki/Time_series) of the solar radiation absorbed by the glass at [discrete time](https://en.m.wikipedia.org/wiki/Discrete_time_and_continuous_time#Discrete_time) $t = [t_0, t_1, ... , t_{n-1}]$.\n", "\n", "_Note_: In Pandas [time series](https://pandas.pydata.org/pandas-docs/version/0.9.1/timeseries.html), time is an index. Therefore, Pandas representation of vector of inputs $u$ in time is the transpose of the matrix presented above. \n", "\n", "Let's consider a [step response](https://en.m.wikipedia.org/wiki/Step_response) in the conditions used for steady-state analysis, i.e., $T_o = 10 \\, \\mathrm{^\\circ C}$, $T_{i,sp} = 20 \\, \\mathrm{^\\circ C}$, and all the flow sources zero (including the HVAC system)." ] }, { "cell_type": "code", "execution_count": 19, "id": "3da47c23-3af1-4515-b919-2b394a655a2f", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Create input_data_set\n", "# ---------------------\n", "# time vector\n", "n = int(np.floor(duration / dt)) # number of time steps\n", "\n", "# DateTimeIndex starting at \"00:00:00\" with a time step of dt\n", "time = pd.date_range(start=\"2000-01-01 00:00:00\",\n", " periods=n, freq=f\"{int(dt)}s\")\n", "\n", "To = 10 * np.ones(n) # outdoor temperature\n", "Ti_sp = 20 * np.ones(n) # indoor temperature set point\n", "Φa = 0 * np.ones(n) # solar radiation absorbed by the glass\n", "Qa = Φo = Φi = Φa # auxiliary heat sources and solar radiation\n", "\n", "data = {'To': To, 'Ti_sp': Ti_sp, 'Φo': Φo, 'Φi': Φi, 'Qa': Qa, 'Φa': Φa}\n", "input_data_set = pd.DataFrame(data, index=time)\n", "\n", "# inputs in time from input_data_set\n", "u = dm4bem.inputs_in_time(us, input_data_set)" ] }, { "cell_type": "markdown", "id": "c26d5657-c137-46c9-93e2-6340e13f9d7d", "metadata": { "tags": [] }, "source": [ "### Time integration" ] }, { "cell_type": "markdown", "id": "c6b0cfa0-301c-4bdb-8f36-8fa6d2bebee3", "metadata": { "tags": [] }, "source": [ "The state-space model\n", "\n", "$$\\left\\{\\begin{array}{rr}\n", "\\dot{\\theta}_C=A_s \\theta_C + B_s u\\\\ \n", "y = C_s \\theta_C + D_s u\n", "\\end{array}\\right.$$\n", "\n", "is integrated in time by using [Euler forward (or explicit) method](https://en.m.wikipedia.org/wiki/Numerical_methods_for_ordinary_differential_equations#Euler_method) for numerical integration:\n", "\n", "$$ \\theta_{s,k+1} = (I + \\Delta t A) \\theta_{s,k} + \\Delta t B u_k $$\n", "\n", "and [Euler backward (or implicit) method](https://en.m.wikipedia.org/wiki/Numerical_methods_for_ordinary_differential_equations#Backward_Euler_method) for numerical integration:\n", "\n", "$$\\theta_{s,k+1} = (I - \\Delta t A)^{-1} ( \\theta_{s,k} + \\Delta t B u_k )$$\n", "\n", "where $k = 0, ... , n - 1$." ] }, { "cell_type": "code", "execution_count": 20, "id": "c6e318a6-782e-4976-8b9e-3f5bd0d01797", "metadata": {}, "outputs": [], "source": [ "# Initial conditions\n", "θ_exp = pd.DataFrame(index=u.index) # empty df with index for explicit Euler\n", "θ_imp = pd.DataFrame(index=u.index) # empty df with index for implicit Euler\n", "\n", "θ0 = 0.0 # initial temperatures\n", "θ_exp[As.columns] = θ0 # fill θ for Euler explicit with initial values θ0\n", "θ_imp[As.columns] = θ0 # fill θ for Euler implicit with initial values θ0\n", "\n", "I = np.eye(As.shape[0]) # identity matrix\n", "for k in range(u.shape[0] - 1):\n", " θ_exp.iloc[k + 1] = (I + dt * As)\\\n", " @ θ_exp.iloc[k] + dt * Bs @ u.iloc[k]\n", " θ_imp.iloc[k + 1] = np.linalg.inv(I - dt * As)\\\n", " @ (θ_imp.iloc[k] + dt * Bs @ u.iloc[k])" ] }, { "cell_type": "markdown", "id": "dc2b5250-4741-4dc7-9e48-f1961e6070b0", "metadata": {}, "source": [ "Then, we obtain the outputs\n", "\n", "$$ y = C_s \\theta_s + D_s u$$\n", "\n", "for explicit and for implicit Euler methods, respectively." ] }, { "cell_type": "code", "execution_count": 21, "id": "696aa2ef-3242-4a7c-a50d-bbf694a15b28", "metadata": {}, "outputs": [], "source": [ "# outputs\n", "y_exp = (Cs @ θ_exp.T + Ds @ u.T).T\n", "y_imp = (Cs @ θ_imp.T + Ds @ u.T).T" ] }, { "cell_type": "code", "execution_count": 22, "id": "049a8d85-8f52-40dc-8e90-41878ee10b4e", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkQAAAHnCAYAAABZgEKHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABpSElEQVR4nO3dd3gU1f4/8PdszaZteoMAoXfQgBSlWSiKiH6vCkpTVFBpIhZEpfwQ0KuAygWlCFwVxSvoRa8NkSpVQkmkSSckIYVkN3Xr+f2xycqSBJLNJpvdfb+eZx8yZ87OfIbxmA9nzpwjCSEEiIiIiHyYzN0BEBEREbkbEyIiIiLyeUyIiIiIyOcxISIiIiKfx4SIiIiIfB4TIiIiIvJ5TIiIiIjI5zEhIiIiIp/HhIiIiIh8HhMiIiIi8nlMiIiIiMjnMSGiOidJUpU+27Ztw5o1ayBJEs6fP+/usKts9+7dmDVrFvLy8twdSpVZrVYEBgbixRdftJcJITBnzhxs3769TmM5fPgw7rvvPjRq1AgajQZhYWHo0aMHPvvss3J1CwoKMGXKFMTFxcHPzw+dO3fGl19+WeFxq1PXXerTfXCXlStXQpIkBAYGVrh///79GDBgAIKCghAYGIh+/frh999/d7oeURkmRFTn9uzZ4/C59957odFoypXfeuutuO+++7Bnzx7Exsa6O+wq2717N2bPnu1RCdGff/6JwsJCdO3a1V526tQpzJw5E+np6XUaS15eHuLj4zFv3jz88MMP+Pe//40mTZpg5MiRmDt3rkPdhx56CGvXrsXMmTPx448/omvXrhg+fDjWrVtX7rjVqesu9ek+uMPly5cxbdo0xMXFVbj/wIED6N27N4qLi/Hpp5/i008/RUlJCe666y7s2bOn2vWIHAgiNxs9erQICAhwdxgu889//lMAEOfOnXN3KFW2cuVKAUCcPn3aXvb5558LAOLUqVNujOxv3bp1E/Hx8fbt//3vfwKAWLdunUO9e+65R8TFxQmz2exUXXfyhPtQmwYPHizuv//+Sv+fMGDAABEdHS0KCwvtZXq9XkRERIiePXtWux7RtdhDRPXa9Y/MZs2aBUmScPToUTz88MPQarUICwvD1KlTYTabcfLkSQwcOBBBQUFo0qQJ3nnnnQqP+9dff+Gxxx5DVFQU1Go12rRpg3/96183jScrKwvPPPMM4uPjoVarERkZidtvvx2//vqrPb6XXnoJAJCQkODw+K865y27zkOHDuGhhx5CcHAwtFotRowYgaysLCf+Jv+2YsUKdOjQAX5+fmjfvj1+/vln7N+/H6GhoWjWrBkAIDExEY8//jgAoGXLlpAkCUFBQRBC1OjcNREREQGFQmHf/uabbxAYGIiHH37Yod4TTzyBtLQ07Nu3z6m6FbnZfXeGp96H2vLZZ59h+/btWLp0aaV1fv/9d/Tt2xf+/v72sqCgIPTu3Ru7d++296JVtV5lauN+U/2nuHkVovrnkUcewYgRIzBu3Dhs3rwZ77zzDkwmE3799Vc899xzmDZtGtatW4dXXnkFzZs3x0MPPWT/7rFjx9CzZ080atQI7733HmJiYvDzzz9j0qRJyM7OxsyZMys978iRI5GUlIS33noLLVu2RF5eHpKSkpCTkwMAeOqpp3D16lV8+OGH2Lhxo/1RX9u2bZ0674MPPohHHnkE48ePx59//ok33ngDx44dw759+6BUKgHYxmT16dPHnnTdyJQpU/Dxxx9j2rRpuPPOO3HixAmMHj0aKpUKXbp0sddbvnw5xo8fD6vVak/YNBoNJEkqd0whBCwWy03PDcAhobkZq9UKq9WK3Nxc/Oc//8HPP/+MJUuW2PenpKSgTZs25Y7ZsWNH+/6ePXtWu25Fbnbfy1T1XtTGfahLrr7nmZmZmDJlChYsWICGDRtWWs9oNEKtVpcrLytLTk5GbGxsletVpqr3m7yMezuoiG78yGz16tUOj59mzpwpAIj33nvPoV7nzp0FALFx40Z7mclkEpGRkeKhhx5yqDtgwADRsGFDodPpHMonTJgg/Pz8xNWrVyuNNTAwUEyZMuWG11PZI7PqnLfsOl944QWHumWPTz777DN7mVwuF3feeecNYxJCiK+//loAEF9++aVD+bx58wQA8dprrzmUR0VFiUmTJt30uFu3bhUAqvSpzmPEcePG2b+nUqnE0qVLHfa3aNFCDBgwoNz30tLSBAAxb948p+pWpCr3XYiq3Yvaug91ydX3/P/+7/9Ez549hdVqFUJU/v+Ezp07i5YtWwqLxWIvM5lMomnTpg6PRKtarzJVvd/kXdhDRB5p8ODBDttt2rTBkSNHMGjQIHuZQqFA8+bNceHCBXtZSUkJtmzZgmeffRb+/v4wm832fffeey+WLFmCvXv3OhznWrfddhvWrFmD8PBw3H333UhMTLT31NyIs+cte1xS5pFHHsHo0aOxdetW+75rj3Uj/+///T907doVjz76qEN527ZtAcChZ+LSpUvIzMxEYmLiTY+bmJiIAwcOVCmGygbLVuS1117DU089hczMTHz33XeYMGECCgsLMW3aNHudG/WUXL+vOnWvV9X7XpV7UVv3oS658p5v2LAB3333HQ4dOnTT+zBx4kSMHTsWEyZMwIwZM2C1WjF79mx7G5fJZNWqVxln2zl5OHdnZETO9BBlZWVV6Rh9+vQR7dq1s2+npqbe9F+0//73vyuNNSsrS0yePFk0btxYABCBgYFi5MiRIj093V6noh6i6p637DpTU1PLxRAdHS2GDh1aaYwVSU9PFwDEokWLyu1bsmSJACAuXbpkL/vmm28EAJGSknLTY1utVmEymar0qYnx48cLhUIhMjMzhRBCdO/eXXTt2rVcvZSUFAFAfPzxx/ay6tStSFXue1W48j4sXrxY/OMf/xDDhw8XQUFB4rbbbhPp6eli4sSJIjQ0VLRr106cP39eCCHExYsXxcCBA0VERITQarXi6aeftveePPfcc+LJJ58UQghhsVjEkCFDxMSJE294Ha665/n5+SI6Olq8+OKLIjc31/4ZPny4CAgIELm5uaKgoMDhOwsWLBCBgYH2dtOjRw/xyiuvCABi586d1a5XEVfdb/IsHFRNPiU0NBRyuRxjxozBgQMHKvzce++9lX4/IiICixcvxvnz53HhwgXMnz8fGzduxJgxY2rlvBkZGQ7bZrMZOTk5CA8Pr9Z1p6amAkCF4ybWrVuHmJgYh7EbBw8ehL+/P1q3bn3TY2/fvh1KpbJKn5rMJ3XbbbfBbDbj7NmzAIAOHTrg+PHj5XplkpOTAQDt27e3l1WnbkWcve/Xc+V9OHr0KPbv348XX3wRmZmZMJlMuOuuuzBkyBBkZmYiISEBa9asAQDk5+dj+vTpSEtLw5EjR/Djjz9iy5YtAIDp06fjq6++wsWLF/Hqq6/CYrFg0aJFN7wOV93z7OxsXLlyBe+99x5CQ0Ptny+++AKFhYUIDQ0t10v6yiuvIDs7G8nJyTh//jx2796N3NxcBAQEOPSkVbVeRVx1v8mz8JEZ+RR/f3/069cPhw4dQseOHaFSqZw+VqNGjTBhwgRs2bLFYcK3soGbxcXFNT7v559/7vA/76+++gpmsxl9+/atVqyRkZEAbIOHr31U8/XXX2P37t3lHkEePXoUrVu3hlwuv+mxa+uR2fW2bt0KmUyGpk2bArANOF+xYgU2bNjgcE1r165FXFwcunXrZi+rTt2bqey+V4Ur78PRo0cxe/Zs+38fzZo1Q7t27XD33XcDAFq3bm0f+Fz2OA4AGjdujO7duyM3NxcA0LBhQ4waNQpDhgwBAOzateum991V9zwmJgZbt24tV75gwQJs374dP/74IyIiIsrtV6vV9iT24sWLWL9+PZ5++mloNBqn6t1ITe43eRYmRORz3n//fdxxxx3o1asXnn32WTRp0gT5+fk4ffo0vvvuO/z2228Vfk+n06Ffv3547LHH0Lp1awQFBeHAgQP46aefHN5i69Chg/08o0ePhlKpRKtWrZw678aNG6FQKHDPPffY3zLr1KkTHnnkEXsdhUKBPn362P/FX5FGjRqha9euWLRoESIjI9GxY0fs2LED77//PgA4TAQIACEhIdi+fTv++9//Ijo6GrGxsWjcuHGFxw4KCnIY91JTzzzzDIKDg3HbbbchOjoa2dnZ+M9//oP169fjpZdesicVgwYNwj333INnn30Wer0ezZs3xxdffIGffvoJn332mcMv9erUvV5V7ztw87fMXHUfrFYrjh07hoEDB9rrHjt2zGGG6+PHj+Oxxx4DYOt9ev/993HmzBmYzWYUFBTgtddes9ft3Lkzli5dit27d1c6Q/S1XHXP/fz8Kkzu16xZA7lcXm5fSkoKNmzYgC5dukCtVuPIkSNYsGABWrRogf/3//5ftetVpDr3m7yMu5/ZEdXlGKIy586dE08++aRo0KCBUCqVIjIyUvTs2VPMnTu30jhLSkrE+PHjRceOHUVwcLDQaDSiVatWYubMmQ4TwAkhxPTp00VcXJyQyWQCgNi6dWu1zlt2nQcPHhT333+/CAwMFEFBQWL48OHiypUrDnUBiD59+lQa97XXPHDgQBEYGChCQkLE/fffL1atWiUAiP/9738Odc+cOSP69OkjAgICKnyrrzZ98sknolevXiIiIkIoFAoREhIi+vTpIz799NNydfPz88WkSZNETEyMUKlUomPHjuKLL76o8LjVqXutqt73/Px8AUAMGzbshsdzxX04efKkiIqKcohRpVI5jLeJj48XycnJ4ueffxatW7cWR44cEWazWWRmZoqAgABhMBiEEELs27dPNGjQQDz88MNi3LhxN/37qAuVteeTJ0+K3r17i7CwMKFSqUTz5s3F66+/Xm6cUVXrVaQ67Zy8iySEF87wReThZs2ahdmzZyMrK6vCRwZU//zwww8YPHgwjhw5Yu8lrC1ff/01li9fjl9++QUAkJSUhGHDhuHUqVMAbMufxMTEoKCgAIsWLcJvv/2Gr7/+GlevXsVTTz2F9PR0HD16FBcvXsQdd9yBf//732jVqhVatWqF5OTkSnsDibwZB1UTEbnA1q1bMWzYsFpPhgDbYPBOnTrZt48cOeKwnZycbJ+I8vHHH0dOTg6ioqIwatQotG3bFp06dUJ+fj4GDx6MmTNnom/fvoiNjcWIESPKrRdH5CvYQ0RUD7GHiIiobjEhIiIiIp/HR2ZERETk85gQERERkc9jQkREREQ+z+cmZrRarUhLS0NQUNBNFxIkIiKi+kEIgfz8fMTFxd10gV5n+FxClJaWhvj4eHeHQURERE64dOmSw5p/ruJzCVFQUBAA219ocHCwm6MhIiKiqtDr9YiPj7f/Hnc1n0uIyh6TBQcHMyEiIiLyMLU13IWDqomIiMjnMSEiIiIin8eEiIiIiHweEyIiIiLyeUyIiIiIyOcxISIiIiKfx4SIiIiIfB4TIiIiIvJ59Soh2rFjB+6//37ExcVBkiR8++23DvuFEJg1axbi4uKg0WjQt29f/Pnnn+4JloiIiLxGvUqICgsL0alTJyxZsqTC/e+88w4WLlyIJUuW4MCBA4iJicE999yD/Pz8Oo6UiIiIvEm9Wrpj0KBBGDRoUIX7hBBYvHgxZsyYgYceeggAsHbtWkRHR2PdunUYN25cXYZKREREXqRe9RDdyLlz55CRkYH+/fvby9RqNfr06YPdu3dX+j2DwQC9Xu/wISIiIrpWveohupGMjAwAQHR0tEN5dHQ0Lly4UOn35s+fj9mzZ9dqbERERJ7KahWwCAGLxQqLxQKzxQyr2QyzxQxhtcBiNsFqscBiMcNiMcFqscJiMUJYrLBYzBBWM6xmC6wWM6xWM6xWC2A1Q1issFrNgNUCYbUdy/azpfTn0jLxdzmsFkjCDGG1AlYLIKyAsEBYrSgsKq7VvwePSYjKXL/KrRDihivfTp8+HVOnTrVv6/V6xMfH11p8RETkHSxWAZPFCrNVwGw2w2Q0wmo2wmwsgcVshMVkgMVkhMVsgMVkgNVkhNVigtVssv1pMV/zswnCYi79mCAsptKkwQxYTaV/miFZbQmGZDXZkgOrCZLVAggzZFYzZMICmTBDEhbIS/+UCTPkwgKZsEAOs/1PubBABgsUpX/KYYUM1tLy0p+v+dNPsrr7r/yG9AZRq8f3mIQoJiYGgK2nKDY21l6emZlZrtfoWmq1Gmq1utbjIyKi6rNaBYwWKwwmK0pMJhiLi2AyFsNkLIbFaIDZWAKrsRgWUwksphJYjSUQ5hJYTSWAyQCr2QCYDRBmI2A1AhYTYLH9KVmMpQmFCZLVaEsoSv+UCxPkVpPtT2HbVsAMhTBDCTMUsP2phAVqmBEoWdz9V+U6lfch3JBFSLBK16ZWto8Fclilv7eFZEuxLJIcAjJYJVnpn3LbPkleWs+2LVD6pyS31xGSHHD4WYaCEjOAr1z6V3Etj0mIEhISEBMTg82bN+OWW24BABiNRmzfvh1vv/22m6MjIvJ8JosVRUYLSowmlBQVoqS4EMaSQphKCmExFsFSUgSrsRAWYwmEsbA0KSmCZCoGzMWQzMWQzAZIFgNkFgMkixFya9nHAIXVCIUwQSmMUMIIpTBBBRPUMEEDE7T1IemoYrJgFRKMki1tMkEBs6QoTQIUsJQmA1ZJYfvlX/anTFn6S19h+2UvU9g+kgJCrgBkCkBSAnI5IFPatuVKSDI5IFcAMiUkmcL2s1wJqbSOJFdCkitKP0rI5NeWKSErLbN95JDJlZDL5fZtuUwBuVwBmUIOuVwFmUIORek+SWaLFTI55JIEOQBlrd6Ayun1euAVH0mICgoKcPr0afv2uXPncPjwYYSFhaFRo0aYMmUK5s2bhxYtWqBFixaYN28e/P398dhjj7kxaiKiumGyWFFoMKOgxIjiogIYCvUwFeXDVKKHuaQIlpJ8WEsKYDUWQBgKIJkKIZmKIJmKIDOXQGYpgdxSDLnFAIW1BCqrAUqrASphgBoG+MEIDQzQSqa6uaAbJB9WSLClTSqYJBXMkhImyfazRaaERaaCWaaGVaayfeQqCJkSkCshZCpArgLkSkCugkyhKv1TCUlh25Yp1JApbX/KlSrIFSooVGrIlWoolLZ9irJ9pWUKpQqSQm07lkwOv7r5W6I6Uq8Soj/++AP9+vWzb5eN/Rk9ejTWrFmDl19+GcXFxXjuueeQm5uLbt264ZdffkFQUJC7QiYiqpAQAgazLYEpLDagKD8XxkIdDIV5MBfrYSrSwVqsgzAUQBgLAEOBLXExFUJuLoLCUgSFpRgqazH8rMXwE8XQoAT+MKChZHBtsDdITExQwAA1DJIaJpkKJskPJpkaZpkfLHI1LHINLAo/WOV+EArbR1L4AUq/0qTDD5LSDwqVH2RKPyhUGijUGihUflCpNVCWftR+GihU/oBCBSj8IJMp4CdJTDqozkhCiNodpVTP6PV6aLVa6HQ6BAcHuzscIqqnSkwW5BebUJCfh2J9Dorzr8JYqIOpKA/mIj0sJXqgRA8Y8iEz5kNhKoDSUgCVuRB+1kL4iyIEohiBKIa/qxOYUlZIKIYfSiQNjDI/GGUaGOX+MMs1MCv8YVX4Qyj9YVUGQFL6QVL6Q6bSQK72h1wdAIXaHwq/AKjU/lBpAqHWBEKtCYDSzx+SUgMoNLbHM0T1QG3//uZ/6UTktUpMFuiKDNDnXUWhLgfF+mwYCnJhKrgKS1EuRHEeZIY8KAx6KE16+Fny4W/JR4AoQDAKEYIiRDo7rqWCXpcSqFAk+aNEFoASeQBM8gCYFAGwKAMglAGAMgCSOgAydSBk6kAoNEFQ+gVA5R8MtX8w1AFB0ARoofQLBFQBkCk1CJAkBNTsr4mIwISIiDyAxSqgKzbhqj4f+pxMFOZdQYk+G5aCbIiiHEhFOVCU5EFlyoXGlAd/ix5BVj2CUYgIFCFaqmZH+HXJjAkKFEiBKJYFwCgPgFERALMyEBZlEIQ6CJI6CDKNFgpNMJT+WqgCQuAXaPtoAkMg12gBdRD85Eo+AiKqp5gQEVGdM5qtyC004GpeHgpz0lCUmw6T/gqs+ZlAUQ5kxVehNOZCbdIhwKxDsNAhBAVoLlVjYrbrkhoDVCiUBaJYHgSDIggmpRYWtRbw00LyD4VcEwJFQBiUgaHwCw6HJjgc/kHhkPmHQqnUIFSSEOravwYiqkeYEBGRS5gtVuQUGHA1OxP6nDQUXU2HUZcBS34mpMIsKEtyoDHmINCSizBrHiIkHaIlY9UOfk1yY4EMBVIQihRaFCtDYFKFwOIXCqsmHLKAcCgCI6AKioAmJBIB2gj4B4dDpgmBWukHzkhGRJVhQkREN1RisiAjrwhZmekoyL4Ew9XLsOjTIS9Ih7IoEwGGLASbsxEq8hAOHaKrMubmmlUUS6CGTh6KQmUYDKowmP3CAf8wyAIjoA6KgJ82Ev6hUQgKjYEyMBxyvxBoZTJoa++SicgHMSEi8lFCCOiLzcjIyYEu/SwKs1NhyEuD0KVBUXQFfiWZCDZlI1xcRRxy0eRmic41vTgFUgDy5aEoVoXB6BcBERAJWWAUVNoY+IXGIDAsFgFhcZAFRcFPFcBxNUTkdkyIiLxUsdGCy1fzkZV+Afor52HIvgBJnwpVYToCSzIQas5ELLLRSiq48YGuSXTyZCEoUEag2C8KZv9oiKAYKEIaQBMah6DIhggKj4M8MBKBSj8E1u7lERG5FBMiIg+lKzYhLTsPV9POoCjjNEw55yHTp8KvOB1aQwYiRTaa4CqaV7Zg4zWJTiH8kaeMRJE6Eib/aEhBsVCExCEgoiGCIxshIKIBpMAYhChUCKmTqyMiqltMiIjqKbPFivS8Yly+fAF5l/+CIesMpNzz8Cu8jDDjZcSJK2iFXMgqe6W8NOExQ45cRSQK/WJgCIiDpG0IRWg8AqKaICQmAerwRgjw03IuGyLyaUyIiNyoxGTBxex8XLl0GgVpJ2DJ+gtK3XkEFaci0pSOhlIW4it7E6s04SmGH3JUcSj0bwBLUEPIQ+PhH9kYobHNEBDZCIqgGETK5Iisu8siIvI4TIiIapnVKpCuK8alS+dx9eJxGDNPQX71LIKLziPWfBmNpStoKZnLf7H0TSwLZMiVRyLfvwFMQY2hCG8CTXQzhDZoAb/IZtAERKChVMUluomIqEJMiIhcpMhoxtkrOlw5fxxFl1OArJMI0p9BhPESGiMd3SuaVLA06TFCiWxVA+QHNIE1NAGqyGYIiWuB0AYtIQ9piAiFChF1ezlERD6FCRFRNRnMFpy9okPa2T+hv5QCZJ5AUP5pNDBdQAspHe2v7+0p7byxQkKOIgb5AY1hCW0KdXRLhMS3RVCDVlBp4xEnk9f9xRAREQAmRESVslgFLuYU4vz5M9CdOwRcSUGw7gQaGM+jmZSGNtfPy1Pa21MiqZGtboKikBaQIlshqGFbhDduC2V4U0Qq/TiWh4ioHmJCRATb4OaTablI/esoCi4egjLrT0QXnUIrXEA/Se9Y+drExy8BhtCWUMa0QWiTjgiKbw8/bSM0lMnKn4SIiOotJkTkc/QlJhxPzUH6qSSYLv2BwOyjaGg4jdZSKjpJpr8rlj7qskCGbHU88kPaQBbbASFNOiO0MRMfIiJvwoSIvJrRbMWJ9DycPXkUhWf3QZN1BE0MJ9FZuoBu1yY/pXlNsaRBdkALWKLaIbDxrQhtmgh5dBtEq/wR7Z5LICKiOsCEiLyGEAKpucVIPnsJuSd+hyptHxoUpKC9dAYdr33DqzT5KZQFIie4LSwxtyCkWReENE2EJjQB8ez1ISLyOUyIyGNZrAInM/Lx54ljyD+1E0GZf6Cd+TgGShf/nr257LV2SYWswNawxt6C0BbdEdi0GwJCExDA5IeIiMCEiDyI0WxF8uU8HDv+J0x/bUVUzj7cIo7jYSnn70ql+c1VVRzyo7rAv1lPRLTuCVVUWzSQK90TOBER1XtMiKjeEkLgREY+/vjzFPTHf0Nk1h50QwpGyjL/riSVDnoObA1r/G0Ib9MHqoSeCAuKQZj7QiciIg/DhIjqlQxdCX4/lY605G3QXtqKrpZDGCm7aNt5zVIW2dr2kBL6Irz9nZDHd0W0OtBdIRMRkRdgQkRuZbUKHEnNw+9HT6D42M9ok78H98iOIlgqslUoTYJyAltAatoHoe3ugbxxT0T7BbsvaCIi8jpMiKjOFRnN2PVXNg4dOQzN6e/Ry7wHz0lnbAOhS1evKFJoURjfDyGd7oOyeT+EB3J+ZyIiqj1MiKhOFBnN2HI8E/v/2I+Q8z+iv7QX/WXnbTtLe4HytK2hbjMImnb3wr9BIvy5thcREdURJkRUa0pMFmw7mYWdSUehPf1fDMZO3C+7YO8FskIGfXQ3BN76f1C0uQ8hwXHuDZiIiHwWEyJyKSEEki7m4dv9f8GY8h3utW7DHFky5DLbvEBWyFHU8HYEdP4/yNoMRkhAhJsjJiIiYkJELpKVb8DGpFTs37cLffTf4WX5LgRJxfbeoMLoLvDvOgKytg8g0J8vxBMRUf3ChIicJoTAHxdy8emuvyA7sQnDZb9inOyE/b+qkoCGUCc+BqnTMASEN3NvsERERDfAhIiqzWC24Psj6fh61xF0zdyINxS/IFKhBwBYJTksLe+FsttT8GvSG+DSGERE5AGYEFGV6UtM+Pfu8/jp9wP4P8O3WCXfBn+lAQBg8o+G8rYnIbt1FGQcHE1ERB6GCRHdVF6REZ/sOodfdu/HWPNX+Fa+CwqFFQBgjuoARa8pULZ9AOBaYURE5KGYEFGl9CUmfLz9DH74PQlPWjdgk3wrVAoLAMCa0AeyO6ZA0bQfIElujpSIiKhmmBBROUazFev2XcCqLUfxuPEr/Cj/GX4KEwBANO0H6c7XIWvYxc1REhERuQ4TIrITQuDnPzPw9g/H0EX3EzYq1iNSobPti+8O6a43IDW5w81REhERuR4TIgIAXMgpxJv//RO6v/ZgkXINOivPAgBEWHNIA+dBatGfj8aIiMhrMSHycUazFct3nMHK31IwGV9gtOoXyCQBoQqE1OcVSN3GAwqVu8MkIiKqVUyIfNjxdD2mfnUEIVf2YJNiORrJsmw7Og2HdPdsICjavQESERHVESZEPshiFVix8yyW/JKCF6XP8YTqZwCA0DaEdP8HQPO73BwhERFR3WJC5GMy9SWY8MUhZJ9PwXrlh2gnu2Db0WUspHtmA+og9wZIRETkBkyIfMjeszmYsO4Q7ijagjWqVfCXDBD+4ZCGfgS07O/u8IiIiNyGCZEPEML2iOzdn47hZdnneEr1o21Hk16QHloBBMe6N0AiIiI3Y0Lk5YxmK6ZvTMbmpBNYqfwQveXJth29Xwb6vgrI5O4NkIiIqB5gQuTFdMUmPPvZQVw6ewzfqt9BUykdQukPaegyoN1Qd4dHRERUbzAh8lKX84rxxOr9kGUew0bVAkRKeYA2HtKwdUBsR3eHR0REVK8wIfJC57ML8fjKfYjRHcZq9bsIRiEQ3R4YsZFzCxEREVWACZGXOZ1ZgMdX7kWrgv1Yrl4EPxiB+O7AY+sBTYi7wyMiIqqXmBB5kTNZBRi2fC9aFh3EStVCqGACmt8DPPJvQOXv7vCIiIjqLSZEXiItrxgjV+5DQuERfKJ+z5YMtboXeHgt1yIjIiK6CZm7A6Cayy4wYMSqfQjVH8ca9T/hBwPQ/G7g4TVMhoiIiKqAPUQerthowdg1B2DMPo+v1P9EAIqBJr2ARz8DFGp3h0dEROQRmBB5MKtVYNp/juB86mV84/cOIpBne5ts2OeAUuPu8IiIiDwGEyIPtvjXU/g5+RI+Uy1GU1wGghsAj/8H8NO6OzQiIiKPwjFEHurH5HR88NtpTFd8ge6yY4AqyJYMBce5OzQiIiKPw4TIA53PLsRLXx/FUNkujFWULtT64EdAdDv3BkZEROShmBB5mBKTBc9+noSGxrN4R7XSVtj7ZaDNYPcGRkRE5ME4hsjDzP7uGM6lZ+EHvyVQwWibeLHvdHeHRURE5NE8qofIbDbj9ddfR0JCAjQaDZo2bYo5c+bAarW6O7Q68duJK/hi/0W8rvwMTZEKBMYAD34MyDzqNhIREdU7HtVD9Pbbb+Ojjz7C2rVr0a5dO/zxxx944oknoNVqMXnyZHeHV6tyC414ZUMy+ssOYIR8CwDJNm4oINzdoREREXk8j0qI9uzZgwceeAD33XcfAKBJkyb44osv8Mcff7g5str3xn9TYMnPwjuaVYAAcPskoFk/d4dFRETkFTzqWcsdd9yBLVu24NSpUwCAI0eOYNeuXbj33nsr/Y7BYIBer3f4eJqfUtLx/dF0zFR+ihChB6LaAf1ed3dYREREXsOjeoheeeUV6HQ6tG7dGnK5HBaLBW+99RaGDx9e6Xfmz5+P2bNn12GUrlVgMGPWpmO4U5aEB+S/A5IMeGAJ1ygjIiJyIY/qIVq/fj0+++wzrFu3DklJSVi7di3effddrF27ttLvTJ8+HTqdzv65dOlSHUZcc+//egqF+qtYoF5tK+gxAWhwq3uDIiIi8jIe1UP00ksv4dVXX8WwYcMAAB06dMCFCxcwf/58jB49usLvqNVqqNWeucjp8XQ9Pvn9PF5TbECUyAHCmvIVeyIiolrgUT1ERUVFkF33irlcLvfK1+6FEHjzvylIEJcwRvGzrfDedwGVv3sDIyIi8kIe1UN0//3346233kKjRo3Qrl07HDp0CAsXLsSTTz7p7tBc7uc/M3Dg/FWsU38KOaxA68FA87vcHRYREZFX8qiE6MMPP8Qbb7yB5557DpmZmYiLi8O4cePw5ptvujs0lzJZrHj7p5PoL/sDPaVkQK4G+s91d1hERERey6MSoqCgICxevBiLFy92dyi16sv9F5GarcO//b6wFfScCIQluDcoIiIiL+ZRY4h8QYHBjMW//oVH5VsRjwwgIAq44wV3h0VEROTVmBDVM5/sOofCwnxMUX1rK+jzMqAOdGtMRERE3o4JUT2SX2LCql3nMEr+CyJELhDSCLi14ukEiIiIyHWqnRD99ttvaNu2bYVLYOh0OrRr1w47d+50SXC+5tO9F2At1uF55fe2gr7TOSM1ERFRHah2QrR48WI8/fTTCA4OLrdPq9Vi3LhxWLhwoUuC8yVFRjNW7jyHJ+Q/QYt8IKIV0PFRd4dFRETkE6qdEB05cgQDBw6sdH///v1x8ODBGgXli9btu4jiQj2eVJZOwtjnZUAmd29QREREPqLaCdGVK1egVCor3a9QKJCVlVWjoHyNwWzB8h1n8ah8G0KQD4QmAG2HujssIiIin1HthKhBgwZITk6udP/Ro0cRGxtbo6B8zf+OpiM3vxDjlT/YCm6fBMg9aoooIiIij1bthOjee+/Fm2++iZKSknL7iouLMXPmTAwePNglwfkCIQRW/34eQ2S7EYNs27xDnR5zd1hEREQ+pdrdEK+//jo2btyIli1bYsKECWjVqhUkScLx48fxr3/9CxaLBTNmzKiNWL3SwQu5SLmci0Xq72wFPZ4DlH7uDYqIiMjHVDshio6Oxu7du/Hss89i+vTpEEIAACRJwoABA7B06VJER0e7PFBvtfr387hd9ieaS5cBVSDQxfsWqiUiIqrvnBqo0rhxY/zwww/Izc3F6dOnIYRAixYtEBoa6ur4vNrlvGL89GcGlsl/sRV0fgzw07o3KCIiIh9Uo5G7oaGh6Nq1q6ti8Tnr9l1AjDUTd6sO2Qq6PuXegIiIiHxUtROi4uJiLF68GDqdDpMnT+YbZU4yW6z4zx+pGKPYAhmsQEIfILKVu8MiIiLySdV+y2zs2LE4ffo0wsPDcffdd9dGTD5h+6ks6PLz8Zhiq63gtqfdGxAREZEPq3YP0bZt27B582a0a9cOM2bMQGZmJqKiomojNq+2/sAl3Cfba5uIMbgh0HKQu0MiIiLyWdVOiPr06YP3338fLVu2RKNGjZgMOSEr34DfTmTiM8V2W0HiGE7ESERE5EbVfmS2cuVKNG7cGFeuXMGWLVtqIyavtzEpFbEiA91lxwFIQOfh7g6JiIjIp1W7WyIgIIATL9aAEAJf/XEJ/yffaSto2hfQNnRrTERERL6u2j1EVDPJl3U4m5WPf5QlRJ0fd29ARERExISorn1/NB3dZCfQUMoC1MFA6/vcHRIREZHPY0JUh6xWge+PpOEf8h22gnYPAip/9wZFRERETIjqUtLFXOTq8nCvbJ+tgI/LiIiI6oVqJ0SvvfYa9u/fXxuxeL3vjqThTtlh+EsGIKQxEH+bu0MiIiIiOJEQpaenY/DgwYiNjcUzzzyD//3vfzAYDLURm1exWAX+l5yBe+V7bQXtHgQkyb1BEREREQAnEqLVq1fjypUr+OqrrxASEoIXX3wREREReOihh7BmzRpkZ2fXRpweb9/ZHBQW6HCn/LCtoN1Qd4ZDRERE13BqDJEkSejVqxfeeecdnDhxAvv370f37t2xYsUKNGjQAL1798a7776Ly5cvuzpej/V9cjrulB2GBkYgtAkQ29ndIREREVEplwyqbtOmDV5++WX8/vvvSE1NxejRo7Fz50588cUXrji8x7NaBbYcv4L7+LiMiIioXnL5AlqRkZEYO3Ysxo4d6+pDe6yUNB3y9Tr0Ux+2FbQd6s5wiIiI6Dp87b4O/Ho8E3fKDkEjGYHQBCC2k7tDIiIiomswIaoDvx67ggHyA7aNdkP5uIyIiKieYUJUyy7nFeNU+lX0kR21FbTiUh1ERET1DROiWrbl+BV0kZ1CsFQE+EcADW51d0hERER0nRonRO+//z4A4OTJk7BarTUOyNv8ejwTd8mSbBst+gMyuXsDIiIionJq/JZZ+/btAQAvvPACTp8+jcDAQLRr1w7t27dH+/btcd99vvuIqMBgxt4zOZgpP2QraDnAvQERERFRhWqcEN11110AgB9++AEAoNfrkZKSgpSUFGzevNmnE6K9Z3LQwHoZzZTpEDIlpGZ3ujskIiIiqkC1EqL8/HzMnj0b33//PbKzs6HVatGqVSvcfvvt+Mc//oFWrVohODgYPXv2RM+ePWsrZo+x63Q27pTZeoekJrcDfsFujoiIiIgqUq2EaNSoUTh06BDGjRuHyMhIFBUV4eWXX8aFCxfw5ptv4r777sOyZcvQoEGD2orXo+w6nY3ZsrLHZQPdGwwRERFVqloJ0S+//ILff/8dnTt3tpfNmDED3333HRQKBd566y3cdttt2LVrFxISElwdq0dJ1xUjIzMTt6lP2Ao4foiIiKjeqtZbZtHR0SgsLKxwX6NGjfDxxx/j+eefx+TJk10SnCfb9Vc2usuOQylZgLBmQFhTd4dERERElahWQjR58mQ8+eSTOHLkSKV1Hn/8cfz22281DszT/X46G7fLUmwbTfu6NRYiIiK6sWo9Mps8eTKuXLmCxMRE3H333Rg6dCisViuka5ai+OKLLxAREeHyQD2JEAK7TufgOSZEREREHqHar93PmzcPDz74IN599128+OKLKC4uRvv27REVFQW9Xo+SkhKsWbOmFkL1HCcy8iEryEBLv8sQkCAl9HJ3SERERHQDTs1D1LVrV6xfvx5GoxFJSUk4deoU9Ho9IiIicOeddyIqKsrVcXqUXX/9/bhMiusMaELdGxARERHdUI0mZlSpVOjevTu6d+/uqni8wu4z2bhP/qdtg4/LiIiI6j0u7upiFqvAHxeu/j2gOqGPewMiIiKim2JC5GKnruQjynARsdJVCLkaaMTeMyIiovqOCZGLHTh/FT1ltsdlUqPugFLj5oiIiIjoZpgQudj+c38nRGjKx2VERESegAmRCwkhcOBcDrrITtoKGt/h3oCIiIioSmolIZLJZLjzzjtx8ODB2jh8vZWaWwxNwQVESnrb+KG4zu4OiYiIiKqgVhKiTz75BH369MGkSZNq4/D11v5zV9G1tHdIapAIKNRujoiIiIiqokbzEFVmzJgxAICZM2fWxuHrrT8uXEUX6ZRtg2+XEREReQyOIXKh/eeu/j1+qFEP9wZDREREVVajhGjnzp0YMWIEevTogcuXLwMAPv30U+zatcslwXmSnAIDdFlpaCZLh4AExHd1d0hERERURU4nRBs2bMCAAQOg0Whw6NAhGAwGAEB+fj7mzZvnsgA9xeFLeUgsGz8U1ZbrlxEREXkQpxOiuXPn4qOPPsKKFSugVCrt5T179kRSUpJLgqvI5cuXMWLECISHh8Pf3x+dO3euF2+zHUnV2QdUc/wQERGRZ3F6UPXJkyfRu3fvcuXBwcHIy8urSUyVys3Nxe23345+/frhxx9/RFRUFM6cOYOQkJBaOV91HLmUhxc4foiIiMgjOZ0QxcbG4vTp02jSpIlD+a5du9C0adOaxlWht99+G/Hx8Vi9erW97Przu4MQAqcuZaC9dN5WwB4iIiIij+L0I7Nx48Zh8uTJ2LdvHyRJQlpaGj7//HNMmzYNzz33nCtjtNu0aRO6dOmChx9+GFFRUbjllluwYsWKG37HYDBAr9c7fFzt0tViNDachEKyQgTFASHxLj8HERER1R6ne4hefvll6HQ69OvXDyUlJejduzfUajWmTZuGCRMmuDJGu7Nnz2LZsmWYOnUqXnvtNezfvx+TJk2CWq3GqFGjKvzO/PnzMXv27FqJp8yR1Dx0lM4AAKSGibV6LiIiInI9SQghqvslk8mE/v374+OPP0bDhg1x7NgxWK1WtG3bFoGBgbURJwBApVKhS5cu2L17t71s0qRJOHDgAPbs2VPhdwwGg/0NOADQ6/WIj4+HTqdDcHCwS+Ka+/0xdN43BYPl+4C7ZgK9prrkuERERGSj1+uh1Wpd+vv7Wk71ECmVSqSkpECSJPj7+6NLly6ujqtCsbGxaNu2rUNZmzZtsGHDhkq/o1aroVbX7hIaR1N1GC2dtW00YA8RERGRp3F6DNGoUaOwatUqV8ZyU7fffjtOnjzpUHbq1Ck0bty4TuO4ltliRerlS4iXZdkKuKArERGRx3F6DJHRaMTKlSuxefNmdOnSBQEBAQ77Fy5cWOPgrvfCCy+gZ8+emDdvHh555BHs378fy5cvx/Lly11+rqo6nVWAlpZTgBwQ4S0g+WndFgsRERE5x+mEKCUlBbfeeisAWy/NtSRJqllUlejatSu++eYbTJ8+HXPmzEFCQgIWL16Mxx9/vFbOVxVHL+nQqfRxmdTgVrfFQURERM5zOiHaunWrK+OossGDB2Pw4MFuOXdFDqfm4S5Z6fihOCZEREREnoir3ddQSmoeOspsr9xzQDUREZFncrqHaM6cOTfc/+abbzp7aI9htlihv3IekQo9hEwBKaaDu0MiIiIiJzidEH3zzTcO2yaTCefOnYNCoUCzZs18IiE6n1OI1tbTto2otoDSz70BERERkVOcTogOHTpUrkyv12PMmDF48MEHaxSUpziWno9OpY/LOKCaiIjIc7l0DFFwcDDmzJmDN954w5WHrbeOp+vRQeKAaiIiIk/n8kHVeXl50Ol0rj5svXQ8TYe2sgu2jdiO7g2GiIiInOb0I7MPPvjAYVsIgfT0dHz66acYOHBgjQPzBNnp5xEmFUBIckiRbdwdDhERETnJ6YRo0aJFDtsymQyRkZEYPXo0pk+fXuPA6rucAgMiCv8CVKUzVHNANRERkcdyOiE6d+6cK+PwOMfT89FWuggAkMXydXsiIiJP5vQYoosXL0IIUek+b3c8XY82ZeOHotu7NxgiIiKqEacTooSEBGRlZZUrz8nJQUJCQo2C8gTH0/VoU9pDhBgmRERERJ7M6YRICFHhIq4FBQXw8/P+8TRn0jKRIKXbNqL5yIyIiMiTVXsM0dSpUwHYVrR/44034O/vb99nsViwb98+dO7c2WUB1kdGsxWKnBOQKwQsmgjIg6LdHRIRERHVQLUTorIZqoUQSE5Ohkqlsu9TqVTo1KkTpk2b5roI66HTmQVoIWzjhzigmoiIyPNVOyHaunUrAOCJJ57A+++/j+DgYJcHVd+dupKPNpItIZI4foiIiMjjOf3a/erVqwEAx44dw8WLF2E0Gh32DxkypGaR1WOnMwvQW1Y6oJrjh4iIiDxejeYhGjp0KJKTkyFJkv0V/LKB1haLxTUR1kOnr+TjGb5hRkRE5DWcfsts0qRJSEhIwJUrV+Dv748///wTO3bsQJcuXbBt2zYXhlj/5F85g2CpGFaZEoho6e5wiIiIqIac7iHas2cPfvvtN0RGRkImk0Emk+GOO+7A/PnzMWnSJPvga29jNFsRmHcSUAKW8FaQyZXuDomIiIhqyOkeIovFgsDAQABAREQE0tLSAACNGzfGyZMnXRNdPXQhpxBNcRkAoIhp6+ZoiIiIyBWc7iFq3749jh49iqZNm6Jbt2545513oFKpsHz5cjRt2tSVMdYrpzML0FyWCgCQIvm4jIiIyBs4nRC9/vrrKCwsBADMnTsXgwcPRq9evRAeHo7169e7LMD65q/MAvSRbL1hiGzt3mCIiIjIJZxOiAYMGGD/uWnTpjh27BiuXr2K0NDQCpf08BZnrugxVrI9MmNCRERE5B2cGkNkMpnQr18/nDp1yqE8LCzMq5MhANBlnEeAZLC9YRbq/YvYEhER+QKnEiKlUomUlBSvT36uZ7EKKHJtSaA5pCkgd7qDjYiIiOoRp98yGzVqFFatWuXKWOq9y7nFaGy9BABQxLRxczRERETkKk53cRiNRqxcuRKbN29Gly5dEBAQ4LB/4cKFNQ6uvvkrMx/NSwdUyyJbuTkaIiIichWnE6KUlBTceuutAFBuLJG3Pko7nVmARFnZgGomRERERN7C6YSobNV7X3L6Sj6GSbY5iBDBhIiIiMhbOD2GCAB27tyJESNGoGfPnrh82dZz8umnn2LXrl0uCa6+yb5yCVqpCAIyILy5u8MhIiIiF3E6IdqwYQMGDBgAjUaDpKQkGAwGAEB+fj7mzZvnsgDrE8VV26NBo7YxoPRzczRERETkKk4nRHPnzsVHH32EFStWQKn8e4HTnj17IikpySXB1Sd5RUbEGi8AAORRfFxGRETkTZxOiE6ePInevXuXKw8ODkZeXl5NYqqXzucU2d8wU0TzlXsiIiJv4nRCFBsbi9OnT5cr37Vrl1cu7nohpxAtuGQHERGRV3I6IRo3bhwmT56Mffv2QZIkpKWl4fPPP8e0adPw3HPPuTLGeuF8dhGaykoXdQ1v4d5giIiIyKWcfu3+5Zdfhk6nQ79+/VBSUoLevXtDrVZj2rRpmDBhgitjrBfSsrIRLeXZNsK9rweMiIjIl9VoMa633noLM2bMwLFjx2C1WtG2bVsEBga6KrZ6xZhpezxoVIVApQl1czRERETkSjVendTf3x+JiYkAvHeGagBQ5J0FYFvUVeXmWIiIiMi1ajQx46pVq9C+fXv4+fnBz88P7du3x8qVK10VW72hKzIh0mgbUK2KbObmaIiIiMjVnO4heuONN7Bo0SJMnDgRPXr0AADs2bMHL7zwAs6fP4+5c+e6LEh3u3C1EI2lKwAARSQHVBMREXkbpxOiZcuWYcWKFRg+fLi9bMiQIejYsSMmTpzoVQnRuexCNJFl2DbC2UNERETkbZx+ZGaxWNClS5dy5YmJiTCbzTUKqr65kFOEBKk0IQrjG2ZERETexumEaMSIEVi2bFm58uXLl+Pxxx+vUVD1TfqVLESVvXLPhIiIiMjr1Ogts1WrVuGXX35B9+7dAQB79+7FpUuXMGrUKEydOtVeb+HChTWL0s0MWbZX7g2qUKg1Ie4NhoiIiFzO6YQoJSUFt956KwDgzJkzAIDIyEhERkYiJSXFXs8bXsVX6s4BACwh7B0iIiLyRk4nRFu3bnVlHPWWvsSEcEMqoASUUc3dHQ4RERHVghrNQ+QLLmT/PaBayVfuiYiIvFKNxhCVlJTg6NGjyMzMhNVqddg3ZMiQGgVWX5zPueaV+7AE9wZDREREtcLphOinn37CqFGjkJ2dXW6fJEmwWCw1Cqy+SM0tRvfSSRk5BxEREZF3cvqR2YQJE/Dwww8jPT0dVqvV4eMtyRAAZGVnIVLS2TbCmBARERF5I6cToszMTEydOhXR0dGujKfeMWfb3qArUYUBfsFujoaIiIhqg9MJ0T/+8Q9s27bNhaHUT4o82yv3Rm0T9wZCREREtcbpMURLlizBww8/jJ07d6JDhw5QKpUO+ydNmlTj4NxNCAH/wlRADsg4oJqIiMhrOZ0QrVu3Dj///DM0Gg22bdvmMAGjJElekRBlFxgRIzIBAJpITspIRETkrZxOiF5//XXMmTMHr776KmQy75zOKDW3CPFSFgBAHtbYzdEQERFRbXE6kzEajXj00Ue9NhkCbK/cNyxNiBDKhIiIiMhbOZ3NjB49GuvXr3dlLPVO6tVCNJRK51kKaeTeYIiIiKjWOP3IzGKx4J133sHPP/+Mjh07lhtUXRcr3M+fPx+vvfYaJk+ejMWLF7v8+LrsVKglE6yQQRbcwOXHJyIiovrB6YQoOTkZt9xyCwA4rG4P1M0K9wcOHMDy5cvRsWPHWjuHOfs8AKBYE4MAufLGlYmIiMhjeeRq9wUFBXj88cexYsUKzJ07t9bOI9NfAgCYg+Jr7RxERETkfjUaEb1z506MGDECPXv2xOXLlwEAn376KXbt2uWS4Crz/PPP47777sPdd99907oGgwF6vd7hUxVCCGgKUwEAMr5hRkRE5NWcTog2bNiAAQMGQKPRICkpCQaDAQCQn5+PefPmuSzA63355ZdISkrC/Pnzq1R//vz50Gq19k98fNV6e3IKjYixls5BFMVJGYmIiLyZ0wnR3Llz8dFHH2HFihUOA6p79uyJpKQklwR3vUuXLmHy5Mn47LPP4OfnV6XvTJ8+HTqdzv65dOlSlb6XmluMeMmWECk4SzUREZFXc3oM0cmTJ9G7d+9y5cHBwcjLy6tJTJU6ePAgMjMzkZiYaC+zWCzYsWMHlixZAoPBALlc7vAdtVoNtVpd7XOl5hahPV+5JyIi8glOJ0SxsbE4ffo0mjRp4lC+a9cuNG1aO8tc3HXXXUhOTnYoe+KJJ9C6dWu88sor5ZKhmrh8tQD9yxIiTspIRETk1ZxOiMaNG4fJkyfjk08+gSRJSEtLw549ezBt2jS8+eabrozRLigoCO3bt3coCwgIQHh4eLnymtJnXoRKssAiKSAPinXpsYmIiKh+cTohevnll6HT6dCvXz+UlJSgd+/eUKvVmDZtGiZMmODKGN3ClHMeAFCkiUWQzHU9T0RERFT/OJ0QAcBbb72FGTNm4NixY7BarWjbti0CAwNdFVuVbNu2rVaOK9fZBl9bgjkHERERkbdzOiG6ePEi4uPj4e/vjy5dupTb16iR5w5EFkIgoCgVkAGysCbuDoeIiIhqmdOv3SckJCArK6tceU5ODhISPPs19XyDGTHC9sq9X2QT9wZDREREtc7phEgIUeGaZQUFBVWeI6i+ytCVoKFkS/ZU4bXzxhwRERHVH9V+ZDZ16lQAtgVc33jjDfj7+9v3WSwW7Nu3D507d3ZZgO6QritBs9KECCEcQ0REROTtqp0QHTp0CICthyg5ORkqlcq+T6VSoVOnTpg2bZrrInSDK3kF6Ilc24a2oXuDISIiolpX7YSobJX7J554Au+//z6Cg4NdHpS76bLSoJQssEIGWWCMu8MhIiKiWub0W2arV692ZRz1iuGq7ZX7QlUkguQ1mpmAiIiIPIDTg6q9mdBdBgAY/Nk7RERE5AuYEFVAWZAGABDBcW6OhIiIiOoCE6IKaEoyAACKUL5hRkRE5AucSohMJhP69euHU6dOuToetys0mBFusa1yrwlnQkREROQLnEqIlEolUlJSKpyY0dNl6EsQI10FAPiFe+7yI0RERFR1Tj8yGzVqFFatWuXKWOqFDF0JYqUc2wbnICIiIvIJTr9TbjQasXLlSmzevBldunRBQECAw/6FCxfWODh3SM8rRLeySRk5qJqIiMgnOJ0QpaSk4NZbbwWAcmOJPPlRWn5WKhSSFRbIIQ+Mdnc4REREVAecTojKZqz2NoaciwCAQnUkgmVyN0dDREREdaFG0zDn5eVh1apVOH78OCRJQtu2bfHkk09Cq9W6Kr46x0kZiYiIfI/Tg6r/+OMPNGvWDIsWLcLVq1eRnZ2NhQsXolmzZkhKSnJljHVKXpAOALAGNXBzJERERFRXnO4heuGFFzBkyBCsWLECCoXtMGazGU899RSmTJmCHTt2uCzIuuRvn5SRb5gRERH5CqcToj/++MMhGQIAhUKBl19+GV26dHFJcHWtxGRBqDkTkAP+nIOIiIjIZzj9yCw4OBgXL14sV37p0iUEBQXVKCh3uaIvQWzZpIwRnKWaiIjIVzidED366KMYO3Ys1q9fj0uXLiE1NRVffvklnnrqKQwfPtyVMdaZdN3fCZHESRmJiIh8htOPzN59911IkoRRo0bBbDYDsC3p8eyzz2LBggUuC7AuXcktQFf7pIwcVE1EROQrnE6IVCoV3n//fcyfPx9nzpyBEALNmzeHv7+/K+OrU/nZqZBLAmYooAiIcnc4REREVEdqNA8RAPj7+6NDhw6uiMXtTLmXAAAFqkiEyJx+mkhEREQexqUTM7Zp0wZjx4713IkZdWkAgBINJ2UkIiLyJS6dmHHRokUePTGjvNCWEJkDY90cCREREdUlTsx4DXVxFgBAFsyEiIiIyJdwYsZrBJiyAQCqkDg3R0JERER1iRMzlioxWRBqsb1y7x/BhIiIiMiXcGLGUln5BkRLtoRIE8qEiIiIyJdwYsZSmfklaCHlAQCkII4hIiIi8iWcmLFUTp4OiVKRbSMw2r3BEBERUZ3ixIylCrIvAwCMkgoqPw+dR4mIiIicUq2EaOrUqVWuu3DhwmoH404lV21zEBUoIxAmSW6OhoiIiOpStRKiQ4cOOWwfPHgQFosFrVq1AgCcOnUKcrkciYmJrouwjlj16QCAEr9IN0dCREREda1aCdHWrVvtPy9cuBBBQUFYu3YtQkNDAQC5ubl44okn0KtXL9dGWQekggwAgNmfi7oSERH5Gqdfu3/vvfcwf/58ezIEAKGhoZg7dy7ee+89lwRXl1TFmQAAKYjrmBEREfkapxMivV6PK1eulCvPzMxEfn5+jYJyB43BNku1UstX7omIiHyN0wnRgw8+iCeeeAJff/01UlNTkZqaiq+//hpjx47FQw895MoYa53FKhBsyQEAaMIbuDkaIiIiqmtOv3b/0UcfYdq0aRgxYgRMJhOEEFAqlRg7diz++c9/ujLGWne10IhI5AEAApkQERER+RynEyJ/f38sXboU//znPx0mZgwICHBlfHUiM78EsaXLdsi1XLaDiIjI19RoYsYtW7Zgy5YtyMzMhNVqddj3ySef1CiwupSdl492UoFtI5CDqomIiHyN0wnR7NmzMWfOHHTp0gWxsbGQPHgyw/wc26SMZiig8A9zczRERERU12o0hmjNmjUYOXKkK+Nxi5KrtmU78pVhCPXgxI6IiIic4/RbZkajET179nRlLG5j1tlmqS5WcZZqIiIiX+R0QvTUU09h3bp1rozFffJts1Sb/JkQERER+SKnH5mVlJRg+fLl+PXXX9GxY0colUqH/Z60uKuiyDZLtQiMdnMkRERE5A5OJ0RHjx5F586dAQApKSkO+zxtgLWfIQsAoOAr90RERD7J6YTo2oVePV2QKQeQAHUoEyIiIiJf5PQYIm9RZDQjXNgmZQzgLNVEREQ+qdo9RFVdp2zjxo3VDsYdcgqMiC6dpVoTxh4iIiIiX1TthEir1dZGHG6Tk1+MDtADACQOqiYiIvJJ1U6IVq9eXRtxuI3+aibkkrBtBPC1eyIiIl/EMUR5tjmICmRBgFx5k9pERETkjXw+ITLkXQEAFCpC3RwJERERuYtHJUTz589H165dERQUhKioKAwdOhQnT56s0TEt+bZJGQ3qcFeESERERB7IoxKi7du34/nnn8fevXuxefNmmM1m9O/fH4WFhc4ftNCWEJk0ES6KkoiIiDyN0xMzusNPP/3ksL169WpERUXh4MGD6N27t1PHlBdn237wZ0JERETkqzwqIbqeTqcDAISFhVVax2AwwGAw2Lf1er3DflXJVQCALIhvmBEREfkqj3pkdi0hBKZOnYo77rgD7du3r7Te/PnzodVq7Z/4+HiH/f4mW0KkCuYcRERERL7KYxOiCRMm4OjRo/jiiy9uWG/69OnQ6XT2z6VLl+z7hBAIsuQBADRhsbUZLhEREdVjHvnIbOLEidi0aRN27NiBhg0b3rCuWq2GWq2ucJ++2Iyw0lmqA8NiXB4nEREReQaPSoiEEJg4cSK++eYbbNu2DQkJCTU6XnahAdGSbRySKpgJERERka/yqITo+eefx7p16/Df//4XQUFByMiwzTKt1Wqh0WiqfbzcPB2aSSW2jQC+ZUZEROSrPGoM0bJly6DT6dC3b1/ExsbaP+vXr3fqePk56QAAI5SAOtiVoRIREZEH8ageIiGES49XrLMt21GgCEWYJLn02EREROQ5PKqHyNWMpQlRkZLrmBEREfkyn06IREEWAMDEdcyIiIh8mk8nRFKhLSEyc9kOIiIin+bTCZGqJMf2QwCX7SAiIvJlPp0Q+Rlty3bIA5kQERER+TLfTojMZZMyMiEiIiLyZT6bEAkhEGi1Lduh0TIhIiIi8mU+mxAVGS3QinwAQEAIEyIiIiJf5rMJUW6hEaFSAQBAHcTX7omIiHyZzyZEusISBEtFAADJnwkRERGRL/PZhKhAl/33hl+I2+IgIiIi9/PZhKioNCEqlAIBuUct6UZEREQu5rMJkbHAlhAVKbjKPRERka/z2YTIlG+blLFEGeLeQIiIiMjtfDYhshbZEiKzOsS9gRAREZHb+WxChOI8AIDFL9S9cRAREZHb+WxCJDfkAQAkTZh7AyEiIiK389mESGm0rWMmD2BCRERE5Ot8NiHyM9nWMVMGRbg5EiIiInI3n02I/IWth8gvmAkRERGRr/PZhChIFAIAAkK5sCsREZGv89mEKESyJURqPjIjIiLyeb6bEMG20j0XdiUiIiKfXcRLLZkBSEAlr91bLBaYTKa6DYocKJVKyOVyd4dBREQ+wGcTIgAwQQmlKsChTAiBjIwM5OXluScochASEoKYmBhIkuTuUIiIyIv5dEJUpAiG9rpftGXJUFRUFPz9/fmL2E2EECgqKkJmZiYAIDY21s0RERGRN/PphMio1DpsWywWezIUHs6xRe6m0WgAAJmZmYiKiuLjMyIiqjU+O6gaACzXLexaNmbI39/fDdFQRcruBcdzERFRbfLphEim0VZYzsdk9QfvBRER1QWfTojk/iHuDoGIiIjqAZ9OiNSBIe4OwaM0adIEixcvtm9LkoRvv/22St+dNWsWOnfuXCtxERER1ZRPJ0SawFB3h+AyY8aMgSRJ5T4DBw6stXOmp6dj0KBBVao7bdo0bNmyxb49ZswYDB06tJYiIyIiqh6ffstMXskYIk81cOBArF692qFMrVbX2vliYmKqXDcwMBCBgYG1FgsREVFN+HQPEfyC3R2BS6nVasTExDh8QkNDsW3bNqhUKuzcudNe97333kNERATS09MBAH379sWECRMwYcIEhISEIDw8HK+//jqEEJWe7/pHZqmpqRg2bBjCwsIQEBCALl26YN++fQAcH5nNmjULa9euxX//+197T9a2bdtc/vdBRERUVT7dQwS/m/cQCSFQbLLUQTCONEq5y96w6tu3L6ZMmYKRI0fiyJEjOH/+PGbMmIEvvvjCYcLDtWvXYuzYsdi3bx/++OMPPPPMM2jcuDGefvrpm56joKAAffr0QYMGDbBp0ybExMQgKSkJVqu1XN1p06bh+PHj0Ov19h6tsLCKl1AhIiKqC76dEKlvnhAVmyxo++bPdRCMo2NzBsBfVb3b8/3335d7LPXKK6/gjTfewNy5c/Hrr7/imWeewZ9//omRI0fiwQcfdKgbHx+PRYsWQZIktGrVCsnJyVi0aFGVEqJ169YhKysLBw4csCc3zZs3r7BuYGAgNBoNDAZDtR67ERER1RbfToi87JFZv379sGzZMoeysuREpVLhs88+Q8eOHdG4cWOHt8XKdO/e3aFXqkePHnjvvfdgsVhuOkv04cOHccstt7Cnh4iIPJJvJ0TqmydEGqUcx+YMqINgyp+3ugICAirtlQGA3bt3AwCuXr2Kq1evIiAgoNK61VW2zAYREZEn8u2EqAo9RJIkVfvRVX105swZvPDCC1ixYgW++uorjBo1Clu2bIFM9ve4+r179zp8Z+/evWjRokWV1hDr2LEjVq5ciatXr1apl0ilUsFiqfuxWURERBXx7bfMqtBD5EkMBgMyMjIcPtnZ2bBYLBg5ciT69++PJ554AqtXr0ZKSgree+89h+9funQJU6dOxcmTJ/HFF1/gww8/xOTJk6t07uHDhyMmJgZDhw7F77//jrNnz2LDhg3Ys2dPhfWbNGmCo0eP4uTJk8jOzuZaZURE5Fae3/XhJAtkgMp1j4zqg59++snhrTEAaNWqFR577DGcP38e3333HQDb/EErV67EI488gnvuucf+OvyoUaNQXFyM2267DXK5HBMnTsQzzzxTpXOrVCr88ssvePHFF3HvvffCbDajbdu2+Ne//lVh/aeffhrbtm1Dly5dUFBQgK1bt6Jv375OXzsREVFNSOJGE814Ib1eD61Wi7TpsYidl+awr6SkBOfOnUNCQgL8/PzcFKF79O3bF507d65wsLU7+fI9ISKiv5X9/tbpdAgOdv0THp99ZGaUcxAwERER2fhsQmSWMSEiIiIiG58dQ2SS8/HLtbh0BhER+TKf7SGyKNhDRERERDY+mxCZ5f7uDoGIiIjqCZ9NiKzsISIiIqJSPpsQCQV7iIiIiMjGZxMiq5I9RERERGTjswmRUHrXLNVERETkPJ9NiKDkI7PqWrNmDUJCQuzbs2bNsi/7URWSJOHbb791eVxEREQ15bsJkcq7EqIxY8Zg6NChdXrOadOmYcuWLVWun56ejkGDBgEAzp8/D0mScPjw4VqKjoiIqOp8dmJGycsSIncIDAxEYGBglevHxMTUYjRERETO89keIpmq6r/IPU3fvn0xceJETJkyBaGhoYiOjsby5ctRWFiIJ554AkFBQWjWrBl+/PFH+3e2bdsGSZLwv//9D506dYKfnx+6deuG5OTkSs9T0SOzTz75BO3atYNarUZsbCwmTJhg33ftI7OEhAQAwC233AJJkrjSPRERuZXPJkSSuoo9REIAxsK6/whRo+tbu3YtIiIisH//fkycOBHPPvssHn74YfTs2RNJSUkYMGAARo4ciaKiIofvvfTSS3j33Xdx4MABREVFYciQITCZTFU657Jly/D888/jmWeeQXJyMjZt2oTmzZtXWHf//v0AgF9//RXp6enYuHFjja6XiIioJnz2kZlcVcW3zExFwLy42g2mIq+lAVWNsQKdOnXC66+/DgCYPn06FixYgIiICDz99NMAgDfffBPLli3D0aNH0b17d/v3Zs6ciXvuuQeALalq2LAhvvnmGzzyyCM3PefcuXPx4osvYvLkyfayrl27Vlg3MjISABAeHs5HaURE5HYe2UO0dOlSJCQkwM/PD4mJidi5c2e1j6Hw8+7X7jt27Gj/WS6XIzw8HB06dLCXRUdHAwAyMzMdvtejRw/7z2FhYWjVqhWOHz9+0/NlZmYiLS0Nd911V01DJyIiqnMe10O0fv16TJkyBUuXLsXtt9+Ojz/+GIMGDcKxY8fQqFGjKh9H4VfFMURKf1tvTV2r4bQASqXSYVuSJIcySZIAAFar9abHKqt7IxoNJ7okIiLP5XE9RAsXLsTYsWPx1FNPoU2bNli8eDHi4+OxbNmyah2nyj1EkmR7dFXXnyokIbVh79699p9zc3Nx6tQptG7d+qbfCwoKQpMmTar8Gr5KpQIAWCwW5wIlIiJyIY/qITIajTh48CBeffVVh/L+/ftj9+7dFX7HYDDAYDDYt/V6PQBAqfHet8xqYs6cOQgPD0d0dDRmzJiBiIiIKs9vNGvWLIwfPx5RUVEYNGgQ8vPz8fvvv2PixInl6kZFRUGj0eCnn35Cw4YN4efnB61W6+KrISIiqhqP6iHKzs6GxWKxj38pEx0djYyMjAq/M3/+fGi1WvsnPj4eAKD2D671eD3RggULMHnyZCQmJiI9PR2bNm2y9+bczOjRo7F48WIsXboU7dq1w+DBg/HXX39VWFehUOCDDz7Axx9/jLi4ODzwwAOuvAwiIqJqkYSo4fvddSgtLQ0NGjTA7t27HQb/vvXWW/j0009x4sSJct+pqIcoPj4eOp0OwcGOSVFJSQnOnTtnH7DtS7Zt24Z+/fohNzfXYXkOd/Ple0JERH/T6/XQarUV/v52BY96ZBYREQG5XF6uNygzM7Ncr1EZtVoNtVpdF+ERERGRh/KoR2YqlQqJiYnYvHmzQ/nmzZvRs2dPN0VFREREns6jeogAYOrUqRg5ciS6dOmCHj16YPny5bh48SLGjx/v7tA8Wt++feFBT0+JiIhcyuMSokcffRQ5OTmYM2cO0tPT0b59e/zwww9o3Lixu0MjIiIiD+VxCREAPPfcc3juuefcHQYRERF5CY8aQ1RXqjJ7M9UN3gsiIqoLHtlDVFtUKhVkMhnS0tIQGRkJlUpVpWUryPWEEDAajcjKyoJMJqvyXEhERETOYEJ0DZlMhoSEBKSnpyMtzQ3rl1E5/v7+aNSoEWQydmYSEVHtYUJ0HZVKhUaNGsFsNnOdLTeTy+VQKBTspSMiolrHhKgCZSvDX79iPBEREXknPocgIiIin8eEiIiIiHweEyIiIiLyeT43hqhseQq9Xu/mSIiIiKiqyn5v19YyUz6XEOXk5AAA4uPj3RwJERERVVdOTg60Wq3Lj+tzCVFYWBgA4OLFi7XyF1pfde3aFQcOHHB3GHWK1+zd9Ho94uPjcenSJQQHB7s7nDrjS/e4DK/Z+1WlPet0OjRq1Mj+e9zVfC4hKpvgT6vV+tT/ROVyuU9dL8Br9hXBwcE+dc2+eI95zb6jKu25tibq5aBqH/H888+7O4Q6x2smb+SL95jXTHVBErU1Oqme0uv10Gq10Ol0Ppl9E3kLtmUi71GV9lzbbd7neojUajVmzpwJtVrt7lCIqAbYlom8R1Xac223eZ/rISIiIiK6ns/1EBERERFdjwkRERER+TwmRB5g6dKlSEhIgJ+fHxITE7Fz506H/cePH8eQIUOg1WoRFBSE7t274+LFizc8ZnJyMvr06QONRoMGDRpgzpw55Wb/3L59OxITE+Hn54emTZvio48+cvm1XW/Hjh24//77ERcXB0mS8O2339r3mUwmvPLKK+jQoQMCAgIQFxeHUaNGIS0t7abHra/XC9z4mgGgoKAAEyZMQMOGDaHRaNCmTRssW7bspsetz9fsq3ypLQO+157Zlj2coHrtyy+/FEqlUqxYsUIcO3ZMTJ48WQQEBIgLFy4IIYQ4ffq0CAsLEy+99JJISkoSZ86cEd9//724cuVKpcfU6XQiOjpaDBs2TCQnJ4sNGzaIoKAg8e6779rrnD17Vvj7+4vJkyeLY8eOiRUrVgilUim+/vrrWr3eH374QcyYMUNs2LBBABDffPONfV9eXp64++67xfr168WJEyfEnj17RLdu3URiYuINj1mfr1eIG1+zEEI89dRTolmzZmLr1q3i3Llz4uOPPxZyuVx8++23lR6zvl+zL/K1tiyE77VntmXPxoSonrvtttvE+PHjHcpat24tXn31VSGEEI8++qgYMWJEtY65dOlSodVqRUlJib1s/vz5Ii4uTlitViGEEC+//LJo3bq1w/fGjRsnunfv7sxlOKWi/6Fcb//+/QKA/ZdKRTzleoWo+JrbtWsn5syZ41B26623itdff73S43jSNfsKX27LQvhee2Zb9jwe98jsRl3OQgjMmjULcXFx0Gg06Nu3L/7888+bHrO+dkcajUYcPHgQ/fv3dyjv378/du/eDavViv/9739o2bIlBgwYgKioKHTr1q1cN+2YMWPQt29f+/aePXvQp08fh1cXBwwYgLS0NJw/f95e5/rzDhgwAH/88QdMJpNLr7MmdDodJElCSEiIvczbrveOO+7Apk2bcPnyZQghsHXrVpw6dQoDBgyw1/HEa2ZbZlu+nre3Z29ty4B3tGePSojWr1+PKVOmYMaMGTh06BB69eqFQYMG2Z+xv/POO1i4cCGWLFmCAwcOICYmBvfccw/y8/MrPaZer8c999yDuLg4HDhwAB9++CHeffddLFy40F7n3LlzuPfee9GrVy8cOnQIr732GiZNmoQNGzbU6vVmZ2fDYrEgOjraoTw6OhoZGRnIzMxEQUEBFixYgIEDB+KXX37Bgw8+iIceegjbt2+314+NjUWjRo3s2xkZGRUes2zfjeqYzWZkZ2e79DqdVVJSgldffRWPPfaYwyRd3na9H3zwAdq2bYuGDRtCpVJh4MCBWLp0Ke644w57HU+7ZrZlG7blv/lCe/bGtgx4UXt2X+dU9d2oy9lqtYqYmBixYMEC+76SkhKh1WrFRx99VOkx63N35OXLlwUAsXv3bofyuXPnilatWtn3Dx8+3GH//fffL4YNG1bpce+55x7xzDPPOJSlpqYKAGLPnj1CCCFatGgh5s2b51Bn165dAoBIT0+vyWVVGW7QxW40GsUDDzwgbrnlFqHT6W54HE+5XiEqvuZ//vOfomXLlmLTpk3iyJEj4sMPPxSBgYFi8+bNlR6nvl8z27KNr7RlIXyvPftKWxbCe9qzx/QQ3azL+dy5c8jIyHDYr1ar0adPH+zevdte5kndkREREZDL5fZ/BZTJzMxEdHQ0IiIioFAo0LZtW4f9bdq0ueGbKTExMRUeE/j7Xx6V1VEoFAgPD3f6mlzBZDLhkUcewblz57B58+abTuHuyddbXFyM1157DQsXLsT999+Pjh07YsKECXj00Ufx7rvvVvq9+nzNbMt/8/W2DPhOe/bGtgx4V3v2mIToZl3OZf8xVLa/jCd1R6pUKiQmJmLz5s0O5Zs3b0bPnj2hUqnQtWtXnDx50mH/qVOn0Lhx40qP26NHD+zYsQNGo9Fe9ssvvyAuLg5NmjSx17n+vL/88gu6dOkCpVJZwytzXtn/PP/66y/8+uuvVWrsnn69JpOp3OrOcrkcVqu10u/V52tmW/6bL7dlwLfasze2ZcDL2rNT/UpucLMu599//10AEGlpaQ77n3rqKTFgwIBKj1vfuyPLXtVdtWqVOHbsmJgyZYoICAgQ58+fF0IIsXHjRqFUKsXy5cvFX3/9JT788EMhl8vFzp077cd49dVXxciRI+3beXl5Ijo6WgwfPlwkJyeLjRs3iuDg4Apf43zhhRfEsWPHxKpVq+rkNc78/Hxx6NAhcejQIQFALFy4UBw6dEhcuHBBmEwmMWTIENGwYUNx+PBhkZ6ebv8YDAaPvN6bXbMQQvTp00e0a9dObN26VZw9e1asXr1a+Pn5iaVLl3rkNbMt+0ZbFsL32rOvtWUhvKs9e0xCZDAYhFwuFxs3bnQonzRpkujdu7c4c+aMACCSkpIc9g8ZMkSMGjWq0uOOHDlSDBkyxKEsKSlJABBnz54VQgjRq1cvMWnSJIc6GzduFAqFQhiNxppcVpX861//Eo0bNxYqlUrceuutYvv27Q77V61aJZo3by78/PxEp06dys1pMXr0aNGnTx+HsqNHj4pevXoJtVotYmJixKxZs+zPZcts27ZN3HLLLUKlUokmTZqIZcuW1cr1XWvr1q0CQLnP6NGjxblz5yrcB0Bs3brVI69XiBtfsxBCpKenizFjxoi4uDjh5+cnWrVqJd577z2H+D3pmtmWfaMtC+F77dnX2rIQ3tWePSYhEsI2cOvZZ591KGvTpo3DwK23337bvs9gMFRp4FZISIjDv0gWLFhQbuBWmzZtHL43fvx4zvFA5CS2ZSLv4S3t2aMSopt1OS9YsEBotVqxceNGkZycLIYPHy5iY2OFXq+3H8PTuiOJvBHbMpH38Jb27FEJkRA37nK2Wq1i5syZIiYmRqjVatG7d2+RnJzs8H1P644k8lZsy0TewxvasyTEddM+EhEREfkYj3ntnoiIiKi2MCEiIiIin8eEiIiIiHweEyIiIiLyeUyIiIiIyOd5TEK0dOlSJCQkwM/PD4mJidi5c6d938aNGzFgwABERERAkiQcPny4Ssfctm0bJElCXl5e7QRNROVU1pZNJhNeeeUVdOjQAQEBAYiLi8OoUaOQlpZ202OyLRPVvRv9Xp41axZat26NgIAAhIaG4u6778a+fftuekx3tmWPSIjWr1+PKVOmYMaMGTh06BB69eqFQYMG2VeBLiwsxO23344FCxa4OVIiupEbteWioiIkJSXhjTfeQFJSEjZu3IhTp05hyJAh7g6biK5zs9/LLVu2xJIlS5CcnIxdu3ahSZMm6N+/P7Kystwc+Q04PYNRHbrtttvE+PHjHcpat24tXn31VYeysrVxDh06VKXjlq07k5ubK4QQIjs7WwwbNkw0aNBAaDQa0b59e7Fu3TqH7/Tp00dMnDhRvPTSSyI0NFRER0eLmTNnOntpRD6lqm25zP79+wUA++KYlWFbJqpb1W3LOp1OABC//vrrDY/rzrZc73uIjEYjDh48iP79+zuU9+/fH7t373bpuUpKSpCYmIjvv/8eKSkpeOaZZzBy5Mhy3Xxr165FQEAA9u3bh3feeQdz5szB5s2bXRoLkbdxpi3rdDpIkoSQkJBqnYttmaj2VLctG41GLF++HFqtFp06darWueqyLdf7hCg7OxsWiwXR0dEO5dHR0cjIyHDpuRo0aIBp06ahc+fOaNq0KSZOnIgBAwbgP//5j0O9jh07YubMmWjRogVGjRqFLl26YMuWLS6NhcjbVLctl5SU4NVXX8Vjjz2G4ODgap2LbZmo9lS1LX///fcIDAyEn58fFi1ahM2bNyMiIqJa56rLtlzvE6IykiQ5bAshypVVZvz48QgMDLR/KmOxWPDWW2+hY8eOCA8PR2BgIH755Rf7M9EyHTt2dNiOjY1FZmZmFa+EyLdVpS2bTCYMGzYMVqsVS5cutZezLRPVHzdry/369cPhw4exe/duDBw4EI888oi9fdXHtqyoVm03iIiIgFwuL/cvyMzMzHLZaWXmzJmDadOm3bTee++9h0WLFmHx4sX2N12mTJkCo9HoUE+pVDpsS5IEq9VapViIfFVV27LJZMIjjzyCc+fO4bfffnPoHWJbJnK/qrblgIAANG/eHM2bN0f37t3RokULrFq1CtOnT6+XbbneJ0QqlQqJiYnYvHkzHnzwQXv55s2b8cADD1TpGFFRUYiKirppvZ07d+KBBx7AiBEjAABWqxV//fUX2rRp41zwRGRXlbZclgz99ddf2Lp1K8LDwx2OwbZM5H7O/l4WQsBgMACon2253idEADB16lSMHDkSXbp0QY8ePbB8+XJcvHgR48ePBwBcvXoVFy9etM9XcvLkSQBATEwMYmJiqnye5s2bY8OGDdi9ezdCQ0OxcOFCZGRk8H+iRC5yo7ZsNpvxj3/8A0lJSfj+++9hsVjs/wINCwuDSqWq8nnYlolq143acmFhId566y0MGTIEsbGxyMnJwdKlS5GamoqHH364Wuepy7bsEQnRo48+ipycHMyZMwfp6elo3749fvjhBzRu3BgAsGnTJjzxxBP2+sOGDQMAzJw5E7Nmzar0uGXdaQqF7a/hjTfewLlz5zBgwAD4+/vjmWeewdChQ6HT6Wrpyoh8y43a8vnz57Fp0yYAQOfOnR2+t3XrVvTt27fS47ItE9WtG7XlkpISnDhxAmvXrkV2djbCw8PRtWtX7Ny5E+3atbvhcd3ZliUhhHD5UT3El19+iaeeegoFBQXuDoWIaoBtmcg7uLMte0QPkasZDAacOXMGS5Yswd133+3ucIjISWzLRN6hPrRlj3nt3pV+/PFHdOvWDQEBAfjggw/cHQ4ROYltmcg71Ie27NOPzIiIiIgAH+0hIiIiIroWEyIiIiLyeR6XEM2fPx9du3ZFUFAQoqKiMHToUPu8Q2WEEJg1axbi4uKg0WjQt29f/Pnnnw51DAYDJk6ciIiICAQEBGDIkCFITU11qJObm4uRI0dCq9VCq9Vi5MiRyMvLq+1LJCIiojrmcQnR9u3b8fzzz2Pv3r3YvHkzzGYz+vfvj8LCQnudd955BwsXLsSSJUtw4MABxMTE4J577kF+fr69zpQpU/DNN9/gyy+/xK5du1BQUIDBgwfDYrHY6zz22GM4fPgwfvrpJ/z00084fPgwRo4cWafXS0RERLXP4wdVZ2VlISoqCtu3b0fv3r0hhEBcXBymTJmCV155BYCtNyg6Ohpvv/02xo0bB51Oh8jISHz66ad49NFHAQBpaWmIj4/HDz/8gAEDBuD48eNo27Yt9u7di27dugEA9u7dix49euDEiRNo1aqV266ZiIiIXMvjeoiuVzZbZVhYGADg3LlzyMjIQP/+/e111Go1+vTpg927dwMADh48CJPJ5FAnLi4O7du3t9fZs2cPtFqtPRkCgO7du0Or1drrEBERkXfw6IRICIGpU6fijjvuQPv27QHAvvbRtSvulm2X7cvIyIBKpUJoaOgN61S08FxUVFS5FX6JiIjIs3n0TNUTJkzA0aNHsWvXrnL7JEly2BZClCu73vV1KqpfleMQERGRZ/HYHqKJEydi06ZN2Lp1Kxo2bGgvL1vd/vpenMzMTHuvUUxMDIxGI3Jzc29Y58qVK+XOm5WVVa73iYiIiDybxyVEQghMmDABGzduxG+//YaEhASH/QkJCYiJicHmzZvtZUajEdu3b0fPnj0BAImJiVAqlQ510tPTkZKSYq/To0cP6HQ67N+/315n37590Ol09jpERETkHTzuLbPnnnsO69atw3//+1+HN720Wi00Gg0A4O2338b8+fOxevVqtGjRAvPmzcO2bdtw8uRJBAUFAQCeffZZfP/991izZg3CwsIwbdo05OTk4ODBg5DL5QCAQYMGIS0tDR9//DEA4JlnnkHjxo3x3Xff1fFVExERUW3yuISosvE7q1evxpgxYwDYepFmz56Njz/+GLm5uejWrRv+9a9/2QdeA0BJSQleeuklrFu3DsXFxbjrrruwdOlSxMfH2+tcvXoVkyZNwqZNmwAAQ4YMwZIlSxASElJr10dERER1z+MSIiIiIiJX87gxRERERESuxoSIiIiIfB4TIiIiIvJ5TIiIiIjI5zEhIiIiIp/HhIiIiIh8HhMiIiIi8nlMiIiIiMjnMSEiIo8ya9YsdO7c2d1hEJGX4UzVRFRvVLY0T5nRo0djyZIlMBgMCA8Pr6OoiMgXMCEionojIyPD/vP69evx5ptv4uTJk/YyjUYDrVbrjtCIyMvxkRkR1RsxMTH2j1arhSRJ5cquf2Q2ZswYDB06FPPmzUN0dDRCQkIwe/ZsmM1mvPTSSwgLC0PDhg3xySefOJzr8uXLePTRRxEaGorw8HA88MADOH/+fN1eMBHVG0yIiMjj/fbbb0hLS8OOHTuwcOFCzJo1C4MHD0ZoaCj27duH8ePHY/z48bh06RIAoKioCP369UNgYCB27NiBXbt2ITAwEAMHDoTRaHTz1RCROzAhIiKPFxYWhg8++ACtWrXCk08+iVatWqGoqAivvfYaWrRogenTp0OlUuH3338HAHz55ZeQyWRYuXIlOnTogDZt2mD16tW4ePEitm3b5t6LISK3ULg7ACKimmrXrh1ksr//fRcdHY327dvbt+VyOcLDw5GZmQkAOHjwIE6fPo2goCCH45SUlODMmTN1EzQR1StMiIjI4ymVSodtSZIqLLNarQAAq9WKxMREfP755+WOFRkZWXuBElG9xYSIiHzOrbfeivXr1yMqKgrBwcHuDoeI6gGOISIin/P4448jIiICDzzwAHbu3Ilz585h+/btmDx5MlJTU90dHhG5ARMiIvI5/v7+2LFjBxo1aoSHHnoIbdq0wZNPPoni4mL2GBH5KE7MSERERD6PPURERETk85gQERERkc9jQkREREQ+jwkRERER+TwmREREROTzmBARERGRz2NCRERERD6PCRERERH5PCZERERE5POYEBEREZHPY0JEREREPu//A+2n27Igan0XAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# plot results\n", "y = pd.concat([y_exp, y_imp], axis=1, keys=['Explicit', 'Implicit'])\n", "# Flatten the two-level column labels into a single level\n", "y.columns = y.columns.get_level_values(0)\n", "\n", "ax = y.plot()\n", "ax.set_xlabel('Time')\n", "ax.set_ylabel('Indoor temperature, $\\\\theta_i$ / °C')\n", "ax.set_title(f'Time step: $dt$ = {dt:.0f} s; $dt_{{max}}$ = {Δtmax:.0f} s')\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "3dde9c4b-5343-4e4e-baf3-f8a20355f61e", "metadata": { "tags": [] }, "source": [ "> Figure 2. Step response to outdoor temperature $T_o$ by using Euler\n", "[implicit](https://en.m.wikipedia.org/wiki/Numerical_methods_for_ordinary_differential_equations#Backward_Euler_method)\n", "and\n", "[explicit](https://en.m.wikipedia.org/wiki/Numerical_methods_for_ordinary_differential_equations#Euler_method) integration.\n", "\n", "The results of explicit and implicit Euler integration are practically identical if $dt < dt_{max}$. If $dt > dt_{max}$, Euler explicit method of integration is numerically unstable.\n", "\n", "The value the indoor temperature obtained after the [settling time](https://en.m.wikipedia.org/wiki/Settling_time) is almost equal to the value obtained in steady-state." ] }, { "cell_type": "code", "execution_count": 23, "id": "d6801343-7156-4c91-b7cd-f93f93229ad2", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Steady-state indoor temperature obtained with:\n", "- DAE model: 10.0000 °C\n", "- state-space model: 10.0000 °C\n", "- steady-state response to step input: 9.9645 °C\n" ] } ], "source": [ "print('Steady-state indoor temperature obtained with:')\n", "print(f'- DAE model: {float(θss[6]):.4f} °C')\n", "print(f'- state-space model: {float(yss):.4f} °C')\n", "print(f'- steady-state response to step input: \\\n", "{y_exp[\"θ6\"].tail(1).values[0]:.4f} °C')" ] }, { "cell_type": "markdown", "id": "8c9bed2e-d665-49cc-87dc-60caf7961acf", "metadata": {}, "source": [ "## Step response to internal load\n", "\n", "Let's consider that the temperature sources are zero and that the flow-rate sources are zero, with the exception of the auxiliary heat gains, $\\dot Q_a$." ] }, { "cell_type": "code", "execution_count": 24, "id": "7fbeb653-de29-4d9f-9eaa-a254662c8c34", "metadata": {}, "outputs": [], "source": [ "# Create input_data_set\n", "# ---------------------\n", "# time vector\n", "n = int(np.floor(duration / dt)) # number of time steps\n", "\n", "# Create a DateTimeIndex starting at \"00:00:00\" with a time step of dt\n", "time = pd.date_range(start=\"2000-01-01 00:00:00\",\n", " periods=n, freq=f\"{int(dt)}s\")\n", "# Create input_data_set\n", "To = 0 * np.ones(n) # outdoor temperature\n", "Ti_sp = 20 * np.ones(n) # indoor temperature set point\n", "Φa = 0 * np.ones(n) # solar radiation absorbed by the glass\n", "Φo = Φi = Φa # solar radiation\n", "Qa = 1000 * np.ones(n) # auxiliary heat sources\n", "data = {'To': To, 'Ti_sp': Ti_sp, 'Φo': Φo, 'Φi': Φi, 'Qa': Qa, 'Φa': Φa}\n", "input_data_set = pd.DataFrame(data, index=time)\n", "\n", "# Get inputs in time from input_data_set\n", "u = dm4bem.inputs_in_time(us, input_data_set)" ] }, { "cell_type": "code", "execution_count": 25, "id": "aed79f67-572f-400d-b13c-e2a131f19d4d", "metadata": {}, "outputs": [], "source": [ "# Initial conditions\n", "θ_exp[As.columns] = θ0 # fill θ for Euler explicit with initial values θ0\n", "θ_imp[As.columns] = θ0 # fill θ for Euler implicit with initial values θ0\n", "\n", "I = np.eye(As.shape[0]) # identity matrix\n", "for k in range(u.shape[0] - 1):\n", " θ_exp.iloc[k + 1] = (I + dt * As)\\\n", " @ θ_exp.iloc[k] + dt * Bs @ u.iloc[k]\n", " θ_imp.iloc[k + 1] = np.linalg.inv(I - dt * As)\\\n", " @ (θ_imp.iloc[k] + dt * Bs @ u.iloc[k])" ] }, { "cell_type": "code", "execution_count": 26, "id": "8e6f5688-2f16-498f-8063-613f7759b828", "metadata": {}, "outputs": [], "source": [ "# outputs\n", "y_exp = (Cs @ θ_exp.T + Ds @ u.T).T\n", "y_imp = (Cs @ θ_imp.T + Ds @ u.T).T" ] }, { "cell_type": "code", "execution_count": 27, "id": "f49b1146-2cd4-4c8c-a095-976147fac7a2", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkQAAAHnCAYAAABZgEKHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABd7UlEQVR4nO3deXxMV/8H8M+dSWayjyyyEcS+U4kttbYILap+raWtrRRVWzWKopZHUa2gVWoruihtaauqrZTYGnssiS1o1JJExJIF2WbO748kw0hCMmsm83m/nnnMnHvuud+b6ZGvc869VxJCCBARERHZMJmlAyAiIiKyNCZEREREZPOYEBEREZHNY0JERERENo8JEREREdk8JkRERERk85gQERERkc1jQkREREQ2jwkRERER2TwmRERERGTzmBARERGRzWNCRGYnSVKJXrt378a6desgSRIuX75s6bBLLCoqCjNnzsTdu3ctHUqJaTQauLi44L333tOWCSEwe/Zs7Nmzx6yxnDhxAi+++CKqVKkCR0dHeHh4oHXr1vj2228L1c3IyMD48ePh7+8PBwcHNG3aFBs3biyy3dLUtZSy9D1YyurVqyFJElxcXIrcfvjwYYSGhsLV1RUuLi7o2LEj/vnnH73rERVgQkRmd+DAAZ3XCy+8AEdHx0LlzZo1w4svvogDBw7Az8/P0mGXWFRUFGbNmmVVCdHp06dx7949NG/eXFsWFxeHGTNmIDEx0ayx3L17FwEBAZg7dy62b9+Or7/+GtWqVcOAAQMwZ84cnbq9e/fG+vXrMWPGDPzxxx9o3rw5+vfvjw0bNhRqtzR1LaUsfQ+WcP36dYSFhcHf37/I7UeOHEG7du3w4MEDfPPNN/jmm2+QmZmJ559/HgcOHCh1PSIdgsjCBg0aJJydnS0dhtF88sknAoCIj4+3dCgltnr1agFAXLx4UVv23XffCQAiLi7OgpE91LJlSxEQEKD9/PvvvwsAYsOGDTr1OnfuLPz9/UVubq5edS3JGr4HU+revbvo0aNHsX8nhIaGCh8fH3Hv3j1tWVpamvDy8hIhISGlrkf0KI4QUZn2+JTZzJkzIUkSTp06hVdffRUqlQoeHh6YMGECcnNzcf78eXTt2hWurq6oVq0aFixYUGS7Fy5cwGuvvQZvb28olUrUq1cPX3zxxVPjuXnzJoYPH46AgAAolUpUrFgRzz77LP7++29tfBMnTgQABAYG6kz/lea4Bed5/Phx9O7dG25ublCpVHjjjTdw8+ZNPX6SD61atQqNGjWCg4MDGjZsiL/++guHDx+Gu7s7atSoAQAICgrC66+/DgCoXbs2JEmCq6srhBAGHdsQXl5esLOz037++eef4eLigldffVWn3pAhQ5CQkIBDhw7pVbcoT/ve9WGt34OpfPvtt9izZw+WLVtWbJ1//vkHHTp0gJOTk7bM1dUV7dq1Q1RUlHYUraT1imOK75vKPrunVyEqe/r06YM33ngDI0aMQEREBBYsWICcnBz8/fffGDVqFMLCwrBhwwZMmjQJNWvWRO/evbX7njlzBiEhIahSpQoWLlwIX19f/PXXXxg7dixSUlIwY8aMYo87YMAAREdH46OPPkLt2rVx9+5dREdH49atWwCAYcOG4fbt2/j888+xZcsW7VRf/fr19Truyy+/jD59+mDkyJE4ffo0pk+fjjNnzuDQoUOwt7cHkLcmq3379tqk60nGjx+PFStWICwsDM899xzOnTuHQYMGQaFQIDg4WFtv5cqVGDlyJDQajTZhc3R0hCRJhdoUQkCtVj/12AB0Epqn0Wg00Gg0uHPnDn788Uf89ddfWLp0qXZ7bGws6tWrV6jNxo0ba7eHhISUum5Rnva9Fyjpd2GK78GcjP2dJycnY/z48Zg/fz4qV65cbL3s7GwolcpC5QVlMTEx8PPzK3G94pT0+6ZyxrIDVERPnjJbu3atzvTTjBkzBACxcOFCnXpNmzYVAMSWLVu0ZTk5OaJixYqid+/eOnVDQ0NF5cqVRWpqqk756NGjhYODg7h9+3axsbq4uIjx48c/8XyKmzIrzXELzvPdd9/VqVswffLtt99qy+RyuXjuueeeGJMQQvz0008CgNi4caNO+dy5cwUA8cEHH+iUe3t7i7Fjxz613cjISAGgRK/STCOOGDFCu59CoRDLli3T2V6rVi0RGhpaaL+EhAQBQMydO1evukUpyfcuRMm+C1N9D+Zk7O/8//7v/0RISIjQaDRCiOL/TmjatKmoXbu2UKvV2rKcnBxRvXp1nSnRktYrTkm/bypfOEJEVql79+46n+vVq4eTJ0+iW7du2jI7OzvUrFkT//33n7YsMzMTO3fuxNtvvw0nJyfk5uZqt73wwgtYunQpDh48qNPOo1q0aIF169bB09MTnTp1QlBQkHak5kn0PW7BdEmBPn36YNCgQYiMjNRue7StJ/nf//6H5s2bo2/fvjrl9evXBwCdkYmrV68iOTkZQUFBT203KCgIR44cKVEMxS2WLcoHH3yAYcOGITk5Gb/99htGjx6Ne/fuISwsTFvnSSMlj28rTd3HlfR7L8l3YarvwZyM+Z1v3rwZv/32G44fP/7U72HMmDEYOnQoRo8ejalTp0Kj0WDWrFnaPi6TyUpVrzj69nOycpbOyIj0GSG6efNmidpo3769aNCggfbztWvXnvov2q+//rrYWG/evCnGjRsnqlatKgAIFxcXMWDAAJGYmKitU9QIUWmPW3Ce165dKxSDj4+P6NWrV7ExFiUxMVEAEIsWLSq0benSpQKAuHr1qrbs559/FgBEbGzsU9vWaDQiJyenRC9DjBw5UtjZ2Ynk5GQhhBCtWrUSzZs3L1QvNjZWABArVqzQlpWmblFK8r2XhDG/h8WLF4tXXnlF9O/fX7i6uooWLVqIxMREMWbMGOHu7i4aNGggLl++LIQQ4sqVK6Jr167Cy8tLqFQq8dZbb2lHT0aNGiXefPNNIYQQarVa9OzZU4wZM+aJ52Gs7zw9PV34+PiI9957T9y5c0f76t+/v3B2dhZ37twRGRkZOvvMnz9fuLi4aPtN69atxaRJkwQAsW/fvlLXK4qxvm+yLlxUTTbF3d0dcrkcgwcPxpEjR4p8vfDCC8Xu7+XlhcWLF+Py5cv477//MG/ePGzZsgWDBw82yXGTkpJ0Pufm5uLWrVvw9PQs1Xlfu3YNAIpcN7Fhwwb4+vrqrN04duwYnJycULdu3ae2vWfPHtjb25foZcj9pFq0aIHc3Fz8+++/AIBGjRrh7NmzhUZlYmJiAAANGzbUlpWmblH0/d4fZ8zv4dSpUzh8+DDee+89JCcnIycnB88//zx69uyJ5ORkBAYGYt26dQCA9PR0TJkyBQkJCTh58iT++OMP7Ny5EwAwZcoU/PDDD7hy5QomT54MtVqNRYsWPfE8jPWdp6Sk4MaNG1i4cCHc3d21r++//x737t2Du7t7oVHSSZMmISUlBTExMbh8+TKioqJw584dODs764yklbReUYz1fZN14ZQZ2RQnJyd07NgRx48fR+PGjaFQKPRuq0qVKhg9ejR27typc8O3goWbDx48MPi43333nc5f3j/88ANyc3PRoUOHUsVasWJFAHmLhx+dqvnpp58QFRVVaAry1KlTqFu3LuRy+VPbNtWU2eMiIyMhk8lQvXp1AHkLzletWoXNmzfrnNP69evh7++Pli1bastKU/dpivveS8KY38OpU6cwa9Ys7X8fNWrUQIMGDdCpUycAQN26dbULnwum4wCgatWqaNWqFe7cuQMAqFy5MgYOHIiePXsCAPbv3//U791Y37mvry8iIyMLlc+fPx979uzBH3/8AS8vr0LblUqlNom9cuUKNm3ahLfeeguOjo561XsSQ75vsi5MiMjmLFmyBG3atEHbtm3x9ttvo1q1akhPT8fFixfx22+/YdeuXUXul5qaio4dO+K1115D3bp14erqiiNHjuDPP//UuYqtUaNG2uMMGjQI9vb2qFOnjl7H3bJlC+zs7NC5c2ftVWZNmjRBnz59tHXs7OzQvn177b/4i1KlShU0b94cixYtQsWKFdG4cWPs3bsXS5YsAQCdGwECQIUKFbBnzx78+uuv8PHxgZ+fH6pWrVpk266urjrrXgw1fPhwuLm5oUWLFvDx8UFKSgp+/PFHbNq0CRMnTtQmFd26dUPnzp3x9ttvIy0tDTVr1sT333+PP//8E99++63OL/XS1H1cSb934OlXmRnre9BoNDhz5gy6du2qrXvmzBmdO1yfPXsWr732GoC80aclS5bg0qVLyM3NRUZGBj744ANt3aZNm2LZsmWIiooq9g7RjzLWd+7g4FBkcr9u3TrI5fJC22JjY7F582YEBwdDqVTi5MmTmD9/PmrVqoX//e9/pa5XlNJ831TOWHrOjsica4gKxMfHizfffFNUqlRJ2Nvbi4oVK4qQkBAxZ86cYuPMzMwUI0eOFI0bNxZubm7C0dFR1KlTR8yYMUPnBnBCCDFlyhTh7+8vZDKZACAiIyNLddyC8zx27Jjo0aOHcHFxEa6urqJ///7ixo0bOnUBiPbt2xcb96Pn3LVrV+Hi4iIqVKggevToIdasWSMAiN9//12n7qVLl0T79u2Fs7NzkVf1mdJXX30l2rZtK7y8vISdnZ2oUKGCaN++vfjmm28K1U1PTxdjx44Vvr6+QqFQiMaNG4vvv/++yHZLU/dRJf3e09PTBQDRr1+/J7ZnjO/h/PnzwtvbWydGhUKhs94mICBAxMTEiL/++kvUrVtXnDx5UuTm5ork5GTh7OwssrKyhBBCHDp0SFSqVEm8+uqrYsSIEU/9eZhDcf35/Pnzol27dsLDw0MoFApRs2ZNMW3atELrjEparyil6edUvkhClMM7fBFZuZkzZ2LWrFm4efNmkVMGVPZs374d3bt3x8mTJ7WjhKby008/YeXKldixYwcAIDo6Gv369UNcXByAvMef+Pr6IiMjA4sWLcKuXbvw008/4fbt2xg2bBgSExNx6tQpXLlyBW3atMHXX3+NOnXqoE6dOoiJiSl2NJCoPOOiaiIiI4iMjES/fv1MngwBeYvBmzRpov188uRJnc8xMTHaG1G+/vrruHXrFry9vTFw4EDUr18fTZo0QXp6Orp3744ZM2agQ4cO8PPzwxtvvFHoeXFEtoIjRERlEEeIiIjMiwkRERER2TxOmREREZHNY0JERERENo8JEREREdk8m7sxo0ajQUJCAlxdXZ/6IEEiIiIqG4QQSE9Ph7+//1Mf0KsPm0uIEhISEBAQYOkwiIiISA9Xr17VeeafsdhcQuTq6gog7wfq5uZm4WiIiIioJNLS0hAQEKD9PW5sNpcQFUyTubm5MSEiIiKyMqZa7sJF1URERGTzmBARERGRzWNCRERERDaPCRERERHZPCZEREREZPOYEBEREZHNY0JERERENo8JEREREdk8JkRERERk85gQERERkc1jQkREREQ2jwkRERER2Tybe7grkS0SGg2EEPmvwu+R/1mj0UBA5NUHIAr+77F9NBoBAU3+NpH/R34dCEDneHnbAd1jS9BAaAQ0AkD+e5F3QEBo8uLOb1/KiyavDvI+F8Sm/azR5O0Lkfe//LYkkXcueXULfiDiYft57/J/Tvl/5reZ/z/dfR45xsM2CyqKgj8exlIQ6yPlD9t85Dwf309bVPDmsVi18TxsR9uwtq7mkY+Pton8n6GATuHDBiBEfh2d9gr/zKRHtxWUPaz8yM8B+fUfHl8n6EfaKVyme6xH///xXYr3SAXx1MrQ/fIfb6fwhuKbLGKDKKq8uHpP2P5oeYnOqeh6ut/HY+WP1S/49Gj9R2s8Xv/xPR/fqu3LxdR/1L0HWcW0bRxMiKhMEUIgVyOQnatBVq4G2dk5yMm8j+zsB8jJzkROViZyc7KgzslCTnYmRG4WRG4OoM7M+zM3C1Bn5/+ZA2hyITS5gEYNoVEDmlxAqAF1/p/526BRQxK5kIQGklBDlv9eJtSQ8l9yoYYEDSA0kKCBlP9LveAXtiQ0QH65DHm/nCWhyd+WV1fCw7qy/L9YpUf2kSAgKygDIEGj/SyDgEw83FZUfQAPjwFAJj38BWWa50MTEZlHWlYJkz49MSEio8jMUePu/RzcfZCNu+n3cS/tFjJTbyHr3h2IB3chZaZCykqFPCsV9jlpsMvJgJ36AeTqB7BXZ0IhMuFQ8JKy4YQsOCMLHlKOpU+tbCnDWY1GSNoUT2j//fjwvYAEIeWdgNBJ24rfBzrbkV9fKvxZenRb4XqPtlNAaN8W3vaw7OF77X6Q8ouK2q+UdR5tP/8cpKfVf+RnWKhu/jFF/tsiz0kq6j+ix+KUHovt0eM99nMrdC5S0fEX+bOQiv756H5Pxf1MpEJlBZ8KjUJIhdstKrbCm4uu/9R2Hjt28fWLa6bk7T8az+O7FXvcItov7r+7J8ZRbJxG2u+xevceZAP4ppTHLDkmRPREmTlqXLl9Hwl3H+Dm3Xu4l3IF6jtXIKVehd39G1Bm3oRLzh24i7vwklLhI6WirpSh38GeMoyhhgw5sEMu7JEr2SFXste+1JIdciUF1JId1DIFNJI91DI7CMkekMmgkewAmQyQ7CBkckCyA2Ty/Jdd3kuSATI7COlhufa9VFAvrw1JJuWXySBJEmSSDJBkkGSP/imHTJIAmSxvu0wOSZIgyR6WS5IckiyvvkySAzIJkiTl7yvLP15+m5Is77h42L4kk0GSAEkmByBBkkmQkF8GCZDl/VBlkpRX95E6Iv9YBWUySco7hiQBUv62/M9Fvc/7ziTtX1pckEhEppSWlgaEMSEiE9JoBK7cvo+ziWmIv5mGtISLkG7FwSntEnyyrqCq7AZqSilog9uwkzSFGyjmN+F9yQmZcldk2bkix94NOQo3qBVu0ChVkBzcIFO6wE7pBDsHF9g5OsNO6QI7pTPsHVxg5+gCe0cXyBXOgL0jYKeEXCaH3LQ/CiIislFlKiHau3cvPvnkExw7dgyJiYn4+eef0atXLwBATk4Opk2bhu3bt+Pff/+FSqVCp06dMH/+fPj7+1s2cCui0QhcupmBo//dwdmrKbh/LQYut2JQR3MRjWX/4jnpOpRS7sMdHvsvJFeyR7rSF5nO/tC4+MNO5QuFygcOFXzhUMEPkos34OINOFSAk9wOTuY9PSIiIr2UqYTo3r17aNKkCYYMGYL/+7//09l2//59REdHY/r06WjSpAnu3LmD8ePHo2fPnjh69KiFIi77hBA4l5SOPXE3ceJSItRXDqFJ7km0lp1Bb+kylFJO3gjPI6M8uTIl7rlUg8azNhz868PBpxYk96qAKgB2Lj5wl3FyhIiIyhdJiJJeq2dekiTpjBAV5ciRI2jRogX+++8/VKlSpUTtpqWlQaVSITU1FW5ubkaKtmzJzFFjT9xN7D6fjHNnT+OZ+/+gk+wYgmRxuqM/ALLt3ZDr0wQOVYMhq9QM8G0EVKiSt26GiIiojDD17+8yNUJUWqmpqZAkCRUqVCi2TlZWFrKyHt67IC0tzQyRmV+uWoP9F1Ow9UQCzpw+gc7qfXhdfhQNZZcB+4f1sp18YVejPWTV2wFVWkPhUR2KUl8pQEREVL5YbUKUmZmJyZMn47XXXntipjhv3jzMmjXLjJGZV3J6Jr4/dBWbD8Uh+N5evCrfi9byM9opMCHJIAJaQ1a/B1CrCxQe1fW4VJKIiKh8s8qEKCcnB/369YNGo8GyZcueWHfKlCmYMGGC9nNaWhoCAgJMHaLJnU9KxxeRF3Ek9ixek/7CVnkEKijuAci/l0SNjpAa9IZUpxskZy8LR0tERFS2WV1ClJOTgz59+iA+Ph67du166jyiUqmEUqk0U3Smdz4pHZ/tvIATMafwjt0v+MRun3ZdkFBVgdRsIKSm/QFVZQtHSkREZD2sKiEqSIYuXLiAyMhIeHp6Wjoks0lOy8T8P84h8sQ5jJL/inDljocLpCs3B0LGQqr7IhdDExER6aFMJUQZGRm4ePGi9nN8fDxOnDgBDw8P+Pv745VXXkF0dDS2bdsGtVqNpKQkAICHhwcUCoWlwjapXLUGX/0Tj6V/n0Mf9e/Yo9gCN+lB3sZqbYHnpgFVWlk2SCIiIitXpi673717Nzp27FiofNCgQZg5cyYCAwOL3C8yMhIdOnQo0TGs6bL7uBvpCPvxJKTrxzDXfg0ayP7L2+DTCOg8E6jxPBdIExGRTbCpy+47dOiAJ+VnZSh3MykhBNbsj0f4n6cxUvoJ7yh/hRwCwqECpC7/A5q+kfeMKyIiIjKKMpUQEZCWmYOwH07i3NlT2GC/FE1ll/I2NOoDKXQu4FLRsgESERGVQ0yIypCLyRkYtv4Iqt2JwjbFUrhJ9yEcVJB6LAEavGzp8IiIiMotJkRlxOH423hr/RH0z9mC9xWbIIMAAlpCeuUrXkJPRERkYkyIyoDtMYmYsDEa06U1eN1+Z15hs0HAC58AduXnHkpERERlFRMiC/vtZAImbjyMhXZf4EX5YQhIkF74BGjxlqVDIyIishlMiCxo26kEvL/pCFbYLUR7+SkIuQJS71VAg16WDo2IiMimMCGykJ1nb+C9jUexVL4kLxmyd4LUbwNQo/B9mIiIiMi0mBBZwKlrdzF2wzF8Kv8CneXHIORKSP03AtXbWzo0IiIim8S7+5nZ1dv38ea6o5ggvkYP+UEImT2kvt8yGSIiIrIgJkRmdD87F8PWH0XnB9sx1O4PAIDUewVQu4uFIyMiIrJtTIjMRAiBab/EwvPmAfzPfl1eYcepQMP/s2hcRERExDVEZvPD0auIij6FP5SfwQ5qoNGrQLuJlg6LiIiIwITILM4lpWHWr6ewVvEF3KUMwK8J0HMpn1RPRERURnDKzMRy1Bq898NJjMSPaCk7B6FwBV5ZC9g7WDo0IiIiyseEyMS+3H0JrkkHMdruVwCA1GMx4FnDskERERGRDk6ZmdC5pDSs3BWLbXYr8x7W+swAoNErlg6LiIiIHsMRIhNRawTe/+kUxkg/oqosGcKtEhA619JhERERURGYEJnIj0evQrp+DEPttgMApO6LAAc3C0dFREREReGUmQmkZeYg/M8zWG+/CnIIoFEfoHaopcMiIiKiYnCEyAQ+33kBXbL+RD3ZVQhHd6DrfEuHRERERE/AESIj+/dmBn765zR22v8IAJA6TgWcPS0cFRERET0JR4iMbOGOOIySbYGHlAFUrAsEDbF0SERERPQUTIiM6GxiGs7ERmOw/K+8gtCPADkH4YiIiMo6JkRGtOTvC5hg9yPsJTVQszNQs5OlQyIiIqISYEJkJKcTUnHpzFG8KDuUV9BphmUDIiIiohJjQmQki/++gLF2WyCTBFCvB+DbyNIhERERUQkxITKCuBvpuHz22MPRofaTLBsQERERlQoTIiP4an88R4eIiIisGBMiA93KyMLR4xwdIiIismZMiAy04dAVvI4/IJMERK0uHB0iIiKyQrxJjgGyctXYcuAMtsl3AwCkVqMsGg8RERHphyNEBth2MhHPP/gLzlIWRMV6QPUOlg6JiIiI9MCEyAAbD8VjsF3eXaml1qMASbJwRERERKQPJkR6upicDq9rEagspUDj6Ak0etXSIREREZGemBDp6Yej1/C6/G8AgKz5m4C9o4UjIiIiIn0xIdJDjlqDg8eOoY38NAQkoNlAS4dEREREBmBCpIdd55LxfFbe6JCo3gGoUMWyAREREZFBmBDp4cfDl/GKfC8AQNZsgIWjISIiIkMxISqlm+lZyL0YiUrSLaiVFYA6L1o6JCIiIjIQE6JS+jM2Ea/IIgEA8iZ9AXsHC0dEREREhmJCVEqRJ+LQRXY07wOny4iIiMoFJkSlkJyWCe9rO6CQ1Mj2asDnlhEREZUTTIhK4Y/YJLwgOwgAUDT5PwtHQ0RERMZSphKivXv3okePHvD394ckSfjll190tgshMHPmTPj7+8PR0REdOnTA6dOnzRffibMIkeUfr8HLZjsuERERmVaZSoju3buHJk2aYOnSpUVuX7BgAcLDw7F06VIcOXIEvr6+6Ny5M9LT000eW1JqJnyv74CdpEG2d2PAo7rJj0lERETmYWfpAB7VrVs3dOvWrchtQggsXrwYU6dORe/evQEA69evh4+PDzZs2IARI0aYNLY/YhPRXTtd9opJj0VERETmVaZGiJ4kPj4eSUlJ6NKli7ZMqVSiffv2iIqKKna/rKwspKWl6bz0cejUabSUnc37wOkyIiKicsVqEqKkpCQAgI+Pj065j4+PdltR5s2bB5VKpX0FBASU+thpmTnwux4BmSSQ5dOMj+ogIiIqZ6wmISogSZLOZyFEobJHTZkyBampqdrX1atXS33MfXEp6CwdBgAoeXUZERFRuWM1CZGvry8AFBoNSk5OLjRq9CilUgk3NzedV2lFnb6E5rLzeR/qvlDq/YmIiKhss5qEKDAwEL6+voiIiNCWZWdnY8+ePQgJCTHZcTUaAXXcTthLatx3q8Gry4iIiMqhMnWVWUZGBi5evKj9HB8fjxMnTsDDwwNVqlTB+PHjMXfuXNSqVQu1atXC3Llz4eTkhNdee81kMZ28dhfNcw8DckBZv+gr4IiIiMi6lamE6OjRo+jYsaP284QJEwAAgwYNwrp16/D+++/jwYMHGDVqFO7cuYOWLVtix44dcHV1NVlMu88mYpDsBABAzukyIiKickkSQghLB2FOaWlpUKlUSE1NLdF6orDwlfg0bSKy7d2gmBwPyMtUDklERGQTSvv7u7SsZg2RJSSnZSLw9j4AgKZGJyZDRERE5RQToieIunQLz8mOAwAc6nO6jIiIqLxiQvQEZ86eRj3ZVWggB2o+b+lwiIiIyESYED2BdHk3ACDdqwng5GHZYIiIiMhkmBAV4+rt+6j74AQAwLHOc5YNhoiIiEyKCVExDlxKwbOy0wAARc0Olg2GiIiITIoJUTH+PRsNb+kuciQlENDC0uEQERGRCTEhKoIQAnb/5V1uf88nGLBTWjgiIiIiMqVSJ0S7du1C/fr1kZaWVmhbamoqGjRogH379hklOEuJT7mHRtknAADO9Xh1GRERUXlX6oRo8eLFeOutt4q8S6RKpcKIESMQHh5ulOAs5eClZLSSnQEA2NfoYNlgiIiIyORKnRCdPHkSXbt2LXZ7ly5dcOzYMYOCsrRrZw5BJd1HltwF8Gti6XCIiIjIxEqdEN24cQP29vbFbrezs8PNmzcNCsrSnK/vBwDc82vFx3UQERHZgFInRJUqVUJMTEyx20+dOgU/Pz+DgjKHa7fvF1memPoADbNPAQBc6vL+Q0RERLag1AnRCy+8gA8//BCZmZmFtj148AAzZsxA9+7djRKcKXVdUvTC7+jLt/CM7AIAQFGjjTlDIiIiIgsp9XzQtGnTsGXLFtSuXRujR49GnTp1IEkSzp49iy+++AJqtRpTp041RaxmcfVcNNykB8iSOULp3cDS4RAREZEZlDoh8vHxQVRUFN5++21MmTIFQggAgCRJCA0NxbJly+Dj42P0QM3m6kEAQKpnU3hz/RAREZFN0Os3ftWqVbF9+3bcuXMHFy9ehBACtWrVgru7u7HjM6vMHDX80k4AMkAR2NrS4RAREZGZGDQE4u7ujubNmxsrFos7nZCKZogDAKhqt7VwNERERGQupV5U/eDBA8ybNw+TJ09GYmKiKWIym4LpvgLn4uIQILsJDWSQAspPokdERERPVuqEaOjQobh48SI8PT3RqVMnU8RkNjlq3YTo/qV/AAC3XGoBSldLhEREREQWUOops927dyMiIgINGjTA1KlTkZycDG9vb1PEZnK5Gg0U+TmhEAIOyScBABr/IEuGRURERGZW6oSoffv2WLJkCWrXro0qVapYbTIE6I4QJaRmonrORUAOuNdsYcGoiIiIyNxKPWW2evVqVK1aFTdu3MDOnTtNEZPZqDUPE6KTV+6goSweAKAIaGapkIiIiMgCSj1C5OzsbNU3XnxUrlqjfX/137NQSfeRK9nDrmI9C0ZFRERE5lbqEaLyJOeREaLsa9EAgFTXWoCdwlIhERERkQXYdEJUMEIkhIDrrdi8Qr+mlguIiIiILMKmE6KCRdWJqZmonnsJAOBWPdiSIREREZEF2HRCVLCoOvbaXTTKX1BtX/kZS4ZEREREFlDqhOiDDz7A4cOHTRGL2eXkT5n9d/kC3KUMqCEHfPiEeyIiIltT6oQoMTER3bt3h5+fH4YPH47ff/8dWVlZpojN5HLzR4gyrsQAANKcqwJ2SkuGRERERBZQ6oRo7dq1uHHjBn744QdUqFAB7733Hry8vNC7d2+sW7cOKSkppojTJLSX3d/Ke6ArKtaxXDBERERkMXqtIZIkCW3btsWCBQtw7tw5HD58GK1atcKqVatQqVIltGvXDp9++imuX79u7HiNKkctkJyeCd+s/wAALpU4XUZERGSLjLKoul69enj//ffxzz//4Nq1axg0aBD27duH77//3hjNm4xaI3D6ehpqyfISN3tf3pCRiIjIFpX6TtVPU7FiRQwdOhRDhw41dtNGl6PRIPZ6KgZI+SNZnDIjIiKySTZ92X2uWuDS1euoIN3LK3APtGxAREREZBE2nhBpkJGUd0PGHKUHoHSxcERERERkCTadEKVn5cIu/UreB/dqFo2FiIiILMemE6ILN9JRGTcBAHae1SwbDBEREVmMwQnRkiVLAADnz5+HRqMxOCBzOpeUjgApGQAguVe1cDRERERkKQZfZdawYUMAwLvvvouLFy/CxcUFDRo0QMOGDdGwYUO8+OKLBgdpKmcT0zBEyhshQgUmRERERLbK4ITo+eefBwBs374dAJCWlobY2FjExsYiIiKiTCdEKRnZ8FPcyvtQIcCywRAREZHFlCohSk9Px6xZs7Bt2zakpKRApVKhTp06ePbZZ/HKK6+gTp06cHNzQ0hICEJCQkwVs1F5Sal5b1x8LBsIERERWUypEqKBAwfi+PHjGDFiBCpWrIj79+/j/fffx3///YcPP/wQL774IpYvX45KlSqZKl6jkkEDDykj74Ozt2WDISIiIospVUK0Y8cO/PPPP2jatKm2bOrUqfjtt99gZ2eHjz76CC1atMD+/fsRGFj2b3LogXTIoQEgAU6elg6HiIiILKRUV5n5+Pjg3r17RW6rUqUKVqxYgXfeeQfjxo0zSnCPy83NxbRp0xAYGAhHR0dUr14ds2fP1vvqtorS3bw3Tp6A3OhPMSEiIiIrUaqEaNy4cXjzzTdx8uTJYuu8/vrr2LVrl8GBFeXjjz/Gl19+iaVLl+Ls2bNYsGABPvnkE3z++ed6tcf1Q0RERASUcsps3LhxuHHjBoKCgtCpUyf06tULGo0GkiRp63z//ffw8vIyeqAAcODAAbz00kvaK9eqVauG77//HkePHtWrPS8UJEQVjRUiERERWaFS35hx7ty5OHDgAFQqFd577z08ePAADRs2RGBgIDw9PfG///0Pn3zyiSliRZs2bbBz507ExcUBAE6ePIn9+/fjhRdeKHafrKwspKWl6bwKaEeIuKCaiIjIpum1cKZ58+bYtGkTsrOzER0djbi4OKSlpcHLywvPPfccvL1Nk2BMmjQJqampqFu3LuRyOdRqNT766CP079+/2H3mzZuHWbNmFbnt4ZQZEyIiIiJbZtBKYoVCgVatWqFVq1bGiueJNm3ahG+//RYbNmxAgwYNcOLECYwfPx7+/v4YNGhQkftMmTIFEyZM0H5OS0tDQEDeTRi9pPzRImfTTPERERGRdbCqS6smTpyIyZMno1+/fgCARo0a4b///sO8efOKTYiUSiWUSmWR21xxP++NQwVThEtERERWwqqedn///n3IZLohy+VyvS+7d8WDvDdKV0NDIyIiIitmVSNEPXr0wEcffYQqVaqgQYMGOH78OMLDw/Hmm2/q1Z6LVDBCpDJilERERGRtrCoh+vzzzzF9+nSMGjUKycnJ8Pf3x4gRI/Dhhx/q1Z4LR4iIiIgIVpYQubq6YvHixVi8eLFR2nORmBARERGRidYQyWQyPPfcczh27JgpmjcariEiIiIiwEQJ0VdffYX27dtj7NixpmjeKOyRCwcpJ+8DEyIiIiKbZpIps8GDBwMAZsyYYYrmjcKpYHQIABRMiIiIiGyZVV12b0zOUmbeG3snPumeiIjIxhmUEO3btw9vvPEGWrdujevXrwMAvvnmG+zfv98owZkS1w8RERFRAb0Tos2bNyM0NBSOjo44fvw4srKyAADp6emYO3eu0QI0lWa+9nlvlG6WDYSIiIgsTu+EaM6cOfjyyy+xatUq2Nvba8tDQkIQHR1tlOBMyUHk35SRI0REREQ2T++E6Pz582jXrl2hcjc3N9y9e9eQmMzCiQkRERER5dM7IfLz88PFixcLle/fvx/Vq1c3KChzcBRcQ0RERER59E6IRowYgXHjxuHQoUOQJAkJCQn47rvvEBYWhlGjRhkzRpNw1HCEiIiIiPLofb35+++/j9TUVHTs2BGZmZlo164dlEolwsLCMHr0aGPGaBJKkZ33RuFs2UCIiIjI4vRKiHJyctClSxesWLECU6dOxZkzZ6DRaFC/fn24uLgYO0aTsEf+XarlSssGQkRERBanV0Jkb2+P2NhYSJIEJycnBAcHGzsuk1Mgf4TIjgkRERGRrdN7DdHAgQOxZs0aY8ZiVnYiN/+Ng2UDISIiIovTew1RdnY2Vq9ejYiICAQHB8PZWXctTnh4uMHBmZJ2ysxOYdlAiIiIyOL0TohiY2PRrFkzAEBcXJzONkmSDIvKDOwLFlVzhIiIiMjm6Z0QRUZGGjMOs7MTXENEREREeWz2aff2BWuIeJUZERGRzdN7hGj27NlP3P7hhx/q27RZ2ItsQAJHiIiIiEj/hOjnn3/W+ZyTk4P4+HjY2dmhRo0aZT4henjZPdcQERER2Tq9E6Ljx48XKktLS8PgwYPx8ssvGxSUOdiJXI4QEREREQAjryFyc3PD7NmzMX36dGM2axJ2ouCyeyZEREREts7oi6rv3r2L1NRUYzdrdLzsnoiIiAroPWX22Wef6XwWQiAxMRHffPMNunbtanBgpqYdIZLzxoxERES2Tu+EaNGiRTqfZTIZKlasiEGDBmHKlCkGB2Zq9topM44QERER2Tq9E6L4+HhjxmF29rwxIxEREeXTew3RlStXIIQodltZZwcuqiYiIqI8eidEgYGBuHnzZqHyW7duITAw0KCgzMEO6vw3nDIjIiKydXonREKIIh/impGRAQcHK0oyOEJERERk80q9hmjChAkA8p5oP336dDg5OWm3qdVqHDp0CE2bNjVagCbHZ5kRERHZvFInRAV3qBZCICYmBgrFw8vWFQoFmjRpgrCwMONFaEqSHJDrva6ciIiIyolSZwORkZEAgCFDhmDJkiVwc3MzelBmw/VDREREBAMuu1+7di0A4MyZM7hy5Qqys7N1tvfs2dOwyMzBjjdlJCIiIgPvQ9SrVy/ExMRAkiTtJfgFC63VarVxIjQljhARERERDLjKbOzYsQgMDMSNGzfg5OSE06dPY+/evQgODsbu3buNGKIJ8QozIiIiggEjRAcOHMCuXbtQsWJFyGQyyGQytGnTBvPmzcPYsWO1i6/LNF5hRkRERDBghEitVsPFxQUA4OXlhYSEBABA1apVcf78eeNEZ2ocISIiIiIYMELUsGFDnDp1CtWrV0fLli2xYMECKBQKrFy5EtWrVzdmjKbDNUREREQEAxKiadOm4d69ewCAOXPmoHv37mjbti08PT2xadMmowVoUnJ7S0dAREREZYDeCVFoaKj2ffXq1XHmzBncvn0b7u7uRT7So0yS8aaMREREpOcaopycHHTs2BFxcXE65R4eHtaTDAFMiIiIiAiAngmRvb09YmNjrSv5KYpMbukIiIiIqAzQ+yqzgQMHYs2aNcaMxfw4QkREREQwYA1RdnY2Vq9ejYiICAQHB8PZ2Vlne3h4uMHBmRxHiIiIiAgGJESxsbFo1qwZABRaS2TKqbTr169j0qRJ+OOPP/DgwQPUrl0ba9asQVBQUOkb4wgRERERwYCEqOCp9+Z0584dPPvss+jYsSP++OMPeHt749KlS6hQoYJ+DTIhIiIiIhiQEAHAvn37sGLFCvz777/48ccfUalSJXzzzTcIDAxEmzZtjBWj1scff4yAgACsXbtWW1atWjX9G2RCRERERDBgUfXmzZsRGhoKR0dHREdHIysrCwCQnp6OuXPnGi3AR23duhXBwcF49dVX4e3tjWeeeQarVq164j5ZWVlIS0vTeWlJXENEREREBiREc+bMwZdffolVq1bB3v7hHZ9DQkIQHR1tlOAe9++//2L58uWoVasW/vrrL4wcORJjx47F119/Xew+8+bNg0ql0r4CAgIebuSiaiIiIoIBCdH58+fRrl27QuVubm64e/euITEVS6PRoFmzZpg7dy6eeeYZjBgxAm+99RaWL19e7D5TpkxBamqq9nX16tWHGzllRkRERDAgIfLz88PFixcLle/fv99kD3f18/ND/fr1dcrq1auHK1euFLuPUqmEm5ubzkuLCRERERHBgIRoxIgRGDduHA4dOgRJkpCQkIDvvvsOYWFhGDVqlDFj1Hr22Wdx/vx5nbK4uDhUrVpVvwaZEBEREREMuMrs/fffR2pqKjp27IjMzEy0a9cOSqUSYWFhGD16tDFj1Hr33XcREhKCuXPnok+fPjh8+DBWrlyJlStX6tcg1xARERERAEkIIQxp4P79+zhz5gw0Gg3q168PFxcXY8VWpG3btmHKlCm4cOECAgMDMWHCBLz11lsl3j8tLQ0qlQqpk13h9vwEoNNM0wVLRERERqH9/Z2aqrv8xUgMnjNycnLS3iXaHA977d69O7p3726cxjhlRkRERDBgDREArFmzBg0bNoSDgwMcHBzQsGFDrF692lixmR4TIiIiIoIBI0TTp0/HokWLMGbMGLRu3RoAcODAAbz77ru4fPky5syZY7QgTYZriIiIiAgGJETLly/HqlWr0L9/f21Zz5490bhxY4wZM8ZKEiKOEBEREZEBU2ZqtRrBwcGFyoOCgpCbm2tQUGbDhIiIiIhgQEL0xhtvFHmH6JUrV+L11183KCiz4bPMiIiICAZeZbZmzRrs2LEDrVq1AgAcPHgQV69excCBAzFhwgRtvfDwcMOiNBWOEBEREREMSIhiY2PRrFkzAMClS5cAABUrVkTFihURGxurrWeOS/H1xkXVREREBAMSosjISGPGYRkcISIiIiIYeB8iq8eEiIiIiGDgGqLMzEycOnUKycnJ0Gg0Ott69uxpUGBmwYSIiIiIYEBC9Oeff2LgwIFISUkptE2SJKjVaoMCMwuuISIiIiIYMGU2evRovPrqq0hMTIRGo9F5WUUyBDAhIiIiIgAGJETJycmYMGECfHx8jBmPeXHKjIiIiGBAQvTKK69g9+7dRgzFApgQEREREQxYQ7R06VK8+uqr2LdvHxo1agR7e3ud7WPHjjU4OJNjQkREREQwICHasGED/vrrLzg6OmL37t06N2CUJMlKEiKuISIiIiIDEqJp06Zh9uzZmDx5MmQyK72dEZ9lRkRERDBgDVF2djb69u1rvckQwCkzIiIiAmBAQjRo0CBs2rTJmLGYHxMiIiIiggFTZmq1GgsWLMBff/2Fxo0bF1pUXWafcP8oJkREREQEAxKimJgYPPPMMwCg83R7oIw/4f5RXFRNRERE4NPuLR0BERERlQEGrYjet28f3njjDYSEhOD69esAgG+++Qb79+83SnAmxxEiIiIiggEJ0ebNmxEaGgpHR0dER0cjKysLAJCeno65c+caLUCT4ggRERERwYCEaM6cOfjyyy+xatUqnQXVISEhiI6ONkpwJscRIiIiIoIBCdH58+fRrl27QuVubm64e/euITGZD0eIiIiICAYkRH5+frh48WKh8v3796N69eoGBWU2TIiIiIgIBiREI0aMwLhx43Do0CFIkoSEhAR89913CAsLw6hRo4wZo+kwISIiIiIYcNn9+++/j9TUVHTs2BGZmZlo164dlEolwsLCMHr0aGPGaDp8lhkREREBkIQQwpAG7t+/jzNnzkCj0aB+/fpwcXExVmwmkZaWBpVKhdTJrnD78ArgWMHSIREREdFTaH9/p6bCzc3N6O3rPUJ05coVBAQEwMnJCcHBwYW2ValSxeDgTI5TZkRERAQD1hAFBgbi5s2bhcpv3bqFwMBAg4IyGyZEREREBAMSIiFEkc8sy8jIgIODg0FBmQ0TIiIiIoIeU2YTJkwAkPcA1+nTp8PJyUm7Ta1W49ChQ2jatKnRAjQp3piRiIiIoEdCdPz4cQB5I0QxMTFQKBTabQqFAk2aNEFYWJjxIjQRDWRAESNcREREZHtKnRAVPOV+yJAhWLJkiUlWepuDMOy5tkRERFSO6L2IZu3atcaMw+w0vAcRERER5bPZYRImRERERFTAZhMiIfEKMyIiIspjwwmRzZ46ERERPUavrCAnJwcdO3ZEXFycseMxGyZEREREVECvrMDe3h6xsbFF3pjRWghYb+xERERkXHoPkwwcOBBr1qwxZixmxREiIiIiKqD3yuLs7GysXr0aERERCA4OhrOzs8728PBwg4MzLY4QERERUR69h0liY2PRrFkzuLm5IS4uDsePH9e+Tpw4YcQQizdv3jxIkoTx48eXel9hxdN9REREZFx6jxAV3LHaUo4cOYKVK1eicePGerbAhIiIiIjyGHQznrt372LNmjU4e/YsJElC/fr18eabb0KlUhkrviJlZGTg9ddfx6pVqzBnzhy92uAaIiIiIiqgd1Zw9OhR1KhRA4sWLcLt27eRkpKC8PBw1KhRA9HR0caMsZB33nkHL774Ijp16vTUullZWUhLS9N5AbzKjIiIiB7Se4To3XffRc+ePbFq1SrY2eU1k5ubi2HDhmH8+PHYu3ev0YJ81MaNGxEdHY0jR46UqP68efMwa9aswhs4QkRERET5DBohmjRpkjYZAgA7Ozu8//77OHr0qFGCe9zVq1cxbtw4fPvtt3BwcCjRPlOmTEFqaqr2dfXqVQCcMiMiIqKH9B4hcnNzw5UrV1C3bl2d8qtXr8LV1dXgwIpy7NgxJCcnIygoSFumVquxd+9eLF26FFlZWZDLdR/aqlQqoVQqC7XFKTMiIiIqoHdC1LdvXwwdOhSffvopQkJCIEkS9u/fj4kTJ6J///7GjFHr+eefR0xMjE7ZkCFDULduXUyaNKlQMvRkTIiIiIgoj94J0aeffgpJkjBw4EDk5uYCyHukx9tvv4358+cbLcBHubq6omHDhjplzs7O8PT0LFT+VJwyIyIionx6J0QKhQJLlizBvHnzcOnSJQghULNmTTg5ORkzPpPhGiIiIiIqYNB9iADAyckJjRo1MkYsetm9e7eee3LKjIiIiPIY9caM9erVw9ChQ01+Y0aj4KM7iIiIKJ9Rb8y4aNEis9yY0Rg4ZUZEREQFrO7GjERERETGpndCdPToUZ1kCHh4Y8bg4GCjBGdSHCEiIiKifHpnBQU3ZnycKW/MaEycMiMiIqICemcFBTdm3LRpE65evYpr165h48aNGDZsmMluzGhcTIiIiIgoj1XdmNGYBK8yIyIiony2e2NGjhARERFRPqu/MaP+OEJEREREeUqVEE2YMKHEdcPDw0sdjFlxUTURERHlK1VCdPz4cZ3Px44dg1qtRp06dQAAcXFxkMvlCAoKMl6EJiI4QkRERET5SpUQRUZGat+Hh4fD1dUV69evh7u7OwDgzp07GDJkCNq2bWvcKE2Bi6qJiIgon97zRgsXLsS8efO0yRAAuLu7Y86cOVi4cKFRgjMpTpkRERFRPr2zgrS0NNy4caNQeXJyMtLT0w0Kyhx4lRkREREV0DsrePnllzFkyBD89NNPuHbtGq5du4affvoJQ4cORe/evY0Zo2lwxoyIiIjy6X3Z/ZdffomwsDC88cYbyMnJgRAC9vb2GDp0KD755BNjxmganDIjIiKifHonRE5OTli2bBk++eQTnRszOjs7GzM+k+GUGRERERUw6MaMO3fuxM6dO5GcnAyNRqOz7auvvjIoMJPjCBERERHl0zshmjVrFmbPno3g4GD4+flBsrLL2PksMyIiIipg0BqidevWYcCAAcaMx3w4QkRERET59M4KsrOzERISYsxYzIwjRERERJRH74Ro2LBh2LBhgzFjMSvBESIiIiLKp/eUWWZmJlauXIm///4bjRs3hr29vc52PtyViIiIrIXeCdGpU6fQtGlTAEBsbKzONqtYYG0NMRIREZFZ6J0QPfqgV6vEhIiIiIjy2fC8kQ2fOhEREeko9QhRSZ9TtmXLllIHY068DxEREREVKHVCpFKpTBGH+XFRNREREeUrdUK0du1aU8RhAUyIiIiIKI/NZgWcMiMiIqICNpsQccqMiIiICthuVsCEiIiIiPLZcFbAKTMiIiLKY8MJkQ2fOhEREemw3ayAA0RERESUz4YTIts9dSIiItJlu1kBEyIiIiLKZ7NZgWBCRERERPlsOCvgIiIiIiLKY7sJkcx2T52IiIh02XBWYMOnTkRERDpsNyvgs8yIiIgonw0nRLZ76kRERKTLdrMCJkRERESUz6qygnnz5qF58+ZwdXWFt7c3evXqhfPnz+vXGKfMiIiIKJ9VJUR79uzBO++8g4MHDyIiIgK5ubno0qUL7t27V/rGmBARERFRPjtLB1Aaf/75p87ntWvXwtvbG8eOHUO7du1K2ZpV5YJERERkQlaVED0uNTUVAODh4VFsnaysLGRlZWk/p6WlAQCEjCNERERElMdqh0mEEJgwYQLatGmDhg0bFltv3rx5UKlU2ldAQED+Frl5AiUiIqIyz2oTotGjR+PUqVP4/vvvn1hvypQpSE1N1b6uXr0KAJC4hoiIiIjyWeWU2ZgxY7B161bs3bsXlStXfmJdpVIJpVJZeAMTIiIiIspnVQmREAJjxozBzz//jN27dyMwMFD/xiROmREREVEeq0qI3nnnHWzYsAG//vorXF1dkZSUBABQqVRwdHQsXWNcVE1ERET5rGoN0fLly5GamooOHTrAz89P+9q0aZMerVnVqRMREZEJWdUIkRDCiK0xISIiIqI8tpsVcMaMiIiI8tlwQmS7p05ERES6bDcrYEJERERE+Ww3K2BCRERERPlsNiuQmBARERFRPpvNCgQTIiIiIspns1mBxMvMiIiIKJ/NJkS8UzUREREVsN2EiM8yIyIionw2nBBxhIiIiIjyWNWjO4zqKQmRWq1GTk6OmYKhotjb20Mu50geERGZng0nREUPjgkhkJSUhLt375o3HipShQoV4OvrC4kjekREZEI2nBAV/Qu2IBny9vaGk5MTfxFbiBAC9+/fR3JyMgDAz8/PwhEREVF5ZrMJUVE3ZlSr1dpkyNPT0wJR0aMcHR0BAMnJyfD29ub0GRERmYwNL6oufOoFa4acnJzMHQ0Vo+C74HouIiIyJdtNiGTFnzqnycoOfhdERGQONpsQSbZ76kRERPQY280K+CyzUqtWrRoWL16s/SxJEn755ZcS7Ttz5kw0bdrUJHEREREZyoazgvI1FTN48GBIklTo1bVrV5MdMzExEd26dStR3bCwMOzcuVP7efDgwejVq5eJIiMiIiodXmVWjnTt2hVr167VKVMqlSY7nq+vb4nruri4wMXFxWSxEBERGaL8ZQUl9YRF1dZKqVTC19dX5+Xu7o7du3dDoVBg37592roLFy6El5cXEhMTAQAdOnTA6NGjMXr0aFSoUAGenp6YNm0ahBDFHu/xKbNr166hX79+8PDwgLOzM4KDg3Ho0CEAulNmM2fOxPr16/Hrr79qR7J2795t9J8HERFRSdnsCFFJn2UmhMCDHLWJgynM0V5utCusOnTogPHjx2PAgAE4efIkLl++jKlTp+L777/XueHh+vXrMXToUBw6dAhHjx7F8OHDUbVqVbz11ltPPUZGRgbat2+PSpUqYevWrfD19UV0dDQ0Gk2humFhYTh79izS0tK0I1oeHh5GOVciIiJ92HBCVLIRogc5atT/8C8TB1PYmdmhcFKU7uvZtm1boWmpSZMmYfr06ZgzZw7+/vtvDB8+HKdPn8aAAQPw8ssv69QNCAjAokWLIEkS6tSpg5iYGCxatKhECdGGDRtw8+ZNHDlyRJvc1KxZs8i6Li4ucHR0RFZWVqmm3YiIiEzFZhOi8riGqGPHjli+fLlOWUFyolAo8O2336Jx48aoWrWqztViBVq1aqUzKtW6dWssXLgQarX6qXeJPnHiBJ555hmO9BARkVWy2YSopFNmjvZynJkdauJgij5uaTk7Oxc7KgMAUVFRAIDbt2/j9u3bcHZ21ju+xxU8ZoOIiMga2XBCVLIRIkmSSj11VRZdunQJ7777LlatWoUffvgBAwcOxM6dOyF7ZHH5wYMHdfY5ePAgatWqVaJniDVu3BirV6/G7du3SzRKpFAooFabf20WERFRUcrfvFEJlccps6ysLCQlJem8UlJSoFarMWDAAHTp0gVDhgzB2rVrERsbi4ULF+rsf/XqVUyYMAHnz5/H999/j88//xzjxo0r0bH79+8PX19f9OrVC//88w/+/fdfbN68GQcOHCiyfrVq1XDq1CmcP38eKSkpfFYZERFZlPUPfeirHCZEf/75p85VYwBQp04dvPbaa7h8+TJ+++03AHn3D1q9ejX69OmDzp07ay+HHzhwIB48eIAWLVpALpdjzJgxGD58eImOrVAosGPHDrz33nt44YUXkJubi/r16+OLL74osv5bb72F3bt3Izg4GBkZGYiMjESHDh30PnciIiJDSOJJN5oph9LS0qBSqRB3aAdqteissy0zMxPx8fEIDAyEg4ODhSK0jA4dOqBp06ZFLra2JFv+ToiI6KGC39+pqalwc3Mzevvlb5ikpMrhjRmJiIhIPzacFZSvZ5kRERGR/mx2DZEkY0L0KD46g4iIbJnNjhCJcriomoiIiPRjs1lBebzsnoiIiPRjs1kBEyIiIiIqYLtZgZGeJE9ERETWz2YTIkkq/bPCiIiIqHyy4YTIZk+diIiIHmO7WQGnzEpt3bp1qFChgvbzzJkztY/9KAlJkvDLL78YPS4iIiJD2XBCVL5OffDgwejVq5dZjxkWFoadO3eWuH5iYiK6desGALh8+TIkScKJEydMFB0REVHJ2eyNGcEbMxrMxcUFLi4uJa7v6+trwmiIiIj0V76GSUpBKsen3qFDB4wZMwbjx4+Hu7s7fHx8sHLlSty7dw9DhgyBq6sratSogT/++EO7z+7duyFJEn7//Xc0adIEDg4OaNmyJWJiYoo9TlFTZl999RUaNGgApVIJPz8/jB49Wrvt0SmzwMBAAMAzzzwDSZL4pHsiIrKo8psVPE1JH+4qBJB9z/wvIQw6vfXr18PLywuHDx/GmDFj8Pbbb+PVV19FSEgIoqOjERoaigEDBuD+/fs6+02cOBGffvopjhw5Am9vb/Ts2RM5OTklOuby5cvxzjvvYPjw4YiJicHWrVtRs2bNIusePnwYAPD3338jMTERW7ZsMeh8iYiIDGG7U2YlXUOUcx+Y62/aWIryQQKgcNZ79yZNmmDatGkAgClTpmD+/Pnw8vLCW2+9BQD48MMPsXz5cpw6dQqtWrXS7jdjxgx07twZQF5SVblyZfz888/o06fPU485Z84cvPfeexg3bpy2rHnz5kXWrVixIgDA09OTU2lERGRxVjlCtGzZMgQGBsLBwQFBQUHYt29fqduQSjpCZKUaN26sfS+Xy+Hp6YlGjRppy3x8fAAAycnJOvu1bt1a+97DwwN16tTB2bNnn3q85ORkJCQk4Pnnnzc0dCIiIrOzuhGiTZs2Yfz48Vi2bBmeffZZrFixAt26dcOZM2dQpUqVkjdU0hEie6e80Rpzs3cybHd7e53PkiTplEn5tx3QaDRPbUsqwS0KHB0dSxkhERFR2WF1wyTh4eEYOnQohg0bhnr16mHx4sUICAjA8uXLS9VOiW/MKEl5U1fmflnoPkkHDx7Uvr9z5w7i4uJQt27dp+7n6uqKatWqlfgyfIVCAQBQq9X6BUpERGREVjVClJ2djWPHjmHy5Mk65V26dEFUVFSR+2RlZSErK0v7OS0tDUDJRj1s0ezZs+Hp6QkfHx9MnToVXl5eJb6/0cyZMzFy5Eh4e3ujW7duSE9Pxz///IMxY8YUquvt7Q1HR0f8+eefqFy5MhwcHKBSqYx8NkRERCVjVSNEKSkpUKvV2vUvBXx8fJCUlFTkPvPmzYNKpdK+AgICAAAKJad4ijJ//nyMGzcOQUFBSExMxNatW7WjOU8zaNAgLF68GMuWLUODBg3QvXt3XLhwoci6dnZ2+Oyzz7BixQr4+/vjpZdeMuZpEBERlYokhIHXd5tRQkICKlWqhKioKJ3Fvx999BG++eYbnDt3rtA+RY0QBQQEIDU1FW5ubjp1MzMzER8fr12wbUt2796Njh074s6dOzqP57A0W/5OiIjoobS0NKhUqiJ/fxuDVU2ZeXl5QS6XFxoNSk5OLjRqVECpVEKpVJojPCIiIrJSVjVlplAoEBQUhIiICJ3yiIgIhISEWCgqIiIisnZWNUIEABMmTMCAAQMQHByM1q1bY+XKlbhy5QpGjhxp6dCsWocOHWBFs6dERERGZXUJUd++fXHr1i3Mnj0biYmJaNiwIbZv346qVataOjQiIiKyUlaXEAHAqFGjMGrUKEuHQUREROWEVa0hMpeS3L2ZzIPfBRERmYNVjhCZikKhgEwmQ0JCAipWrAiFQsEbOFqIEALZ2dm4efMmZDJZie+FREREpA8mRI+QyWQIDAxEYmIiEhIs8PwyKsTJyQlVqlSBrJw/jJeIiCyLCdFjFAoFqlSpgtzcXD5ny8Lkcjns7Ow4SkdERCbHhKgIBU+Gf/yJ8URERFQ+cR6CiIiIbB4TIiIiIrJ5TIiIiIjI5tncGqKCx1OkpaVZOBIiIiIqqYLf26Z6zJTNJUS3bt0CAAQEBFg4EiIiIiqtW7duQaVSGb1dm0uIPDw8AABXrlwxyQ+0rGrevDmOHDli6TDMiudcvqWlpSEgIABXr16Fm5ubpcMxG1v6jgvwnMu/kvTn1NRUVKlSRft73NhsLiEquMGfSqWyqb9E5XK5TZ0vwHO2FW5ubjZ1zrb4HfOcbUdJ+rOpbtTLRdU24p133rF0CGbHc6byyBa/Y54zmYMkTLU6qYxKS0uDSqVCamqqTWbfROUF+zJR+VGS/mzqPm9zI0RKpRIzZsyAUqm0dChEZAD2ZaLyoyT92dR93uZGiIiIiIgeZ3MjRERERESPY0JERERENo8JkRVYtmwZAgMD4eDggKCgIOzbt09n+9mzZ9GzZ0+oVCq4urqiVatWuHLlyhPbjImJQfv27eHo6IhKlSph9uzZhe7+uWfPHgQFBcHBwQHVq1fHl19+afRze9zevXvRo0cP+Pv7Q5Ik/PLLL9ptOTk5mDRpEho1agRnZ2f4+/tj4MCBSEhIeGq7ZfV8gSefMwBkZGRg9OjRqFy5MhwdHVGvXj0sX778qe2W5XO2VbbUlwHb68/sy1ZOUJm2ceNGYW9vL1atWiXOnDkjxo0bJ5ydncV///0nhBDi4sWLwsPDQ0ycOFFER0eLS5cuiW3btokbN24U22Zqaqrw8fER/fr1EzExMWLz5s3C1dVVfPrpp9o6//77r3BychLjxo0TZ86cEatWrRL29vbip59+Mun5bt++XUydOlVs3rxZABA///yzdtvdu3dFp06dxKZNm8S5c+fEgQMHRMuWLUVQUNAT2yzL5yvEk89ZCCGGDRsmatSoISIjI0V8fLxYsWKFkMvl4pdffim2zbJ+zrbI1vqyELbXn9mXrRsTojKuRYsWYuTIkTpldevWFZMnTxZCCNG3b1/xxhtvlKrNZcuWCZVKJTIzM7Vl8+bNE/7+/kKj0QghhHj//fdF3bp1dfYbMWKEaNWqlT6noZei/kJ53OHDhwUA7S+VoljL+QpR9Dk3aNBAzJ49W6esWbNmYtq0acW2Y03nbCtsuS8LYXv9mX3Z+ljdlNmThpyFEJg5cyb8/f3h6OiIDh064PTp009ts6wOR2ZnZ+PYsWPo0qWLTnmXLl0QFRUFjUaD33//HbVr10ZoaCi8vb3RsmXLQsO0gwcPRocOHbSfDxw4gPbt2+tcuhgaGoqEhARcvnxZW+fx44aGhuLo0aPIyckx6nkaIjU1FZIkoUKFCtqy8na+bdq0wdatW3H9+nUIIRAZGYm4uDiEhoZq61jjObMvsy8/rrz35/Lal4Hy0Z+tKiHatGkTxo8fj6lTp+L48eNo27YtunXrpp1jX7BgAcLDw7F06VIcOXIEvr6+6Ny5M9LT04ttMy0tDZ07d4a/vz+OHDmCzz//HJ9++inCw8O1deLj4/HCCy+gbdu2OH78OD744AOMHTsWmzdvNun5pqSkQK1Ww8fHR6fcx8cHSUlJSE5ORkZGBubPn4+uXbtix44dePnll9G7d2/s2bNHW9/Pzw9VqlTRfk5KSiqyzYJtT6qTm5uLlJQUo56nvjIzMzF58mS89tprOjfpKm/n+9lnn6F+/fqoXLkyFAoFunbtimXLlqFNmzbaOtZ2zuzLediXH7KF/lwe+zJQjvqz5QanSu9JQ84ajUb4+vqK+fPna7dlZmYKlUolvvzyy2LbLMvDkdevXxcARFRUlE75nDlzRJ06dbTb+/fvr7O9R48eol+/fsW227lzZzF8+HCdsmvXrgkA4sCBA0IIIWrVqiXmzp2rU2f//v0CgEhMTDTktEoMTxhiz87OFi+99JJ45plnRGpq6hPbsZbzFaLoc/7kk09E7dq1xdatW8XJkyfF559/LlxcXERERESx7ZT1c2ZfzmMrfVkI2+vPttKXhSg//dlqRoieNuQcHx+PpKQkne1KpRLt27dHVFSUtsyahiO9vLwgl8u1/wookJycDB8fH3h5ecHOzg7169fX2V6vXr0nXpni6+tbZJvAw395FFfHzs4Onp6eep+TMeTk5KBPnz6Ij49HRETEU2/hbs3n++DBA3zwwQcIDw9Hjx490LhxY4wePRp9+/bFp59+Wux+Zfmc2ZcfsvW+DNhOfy6PfRkoX/3ZahKipw05F/zHUNz2AtY0HKlQKBAUFISIiAid8oiICISEhEChUKB58+Y4f/68zva4uDhUrVq12HZbt26NvXv3Ijs7W1u2Y8cO+Pv7o1q1ato6jx93x44dCA4Ohr29vYFnpr+CvzwvXLiAv//+u0Sd3drPNycnp9DTneVyOTQaTbH7leVzZl9+yJb7MmBb/bk89mWgnPVnvcaVLOBpQ87//POPACASEhJ0tg8bNkyEhoYW225ZH44suFR3zZo14syZM2L8+PHC2dlZXL58WQghxJYtW4S9vb1YuXKluHDhgvj888+FXC4X+/bt07YxefJkMWDAAO3nu3fvCh8fH9G/f38RExMjtmzZItzc3Iq8jPPdd98VZ86cEWvWrDHLZZzp6eni+PHj4vjx4wKACA8PF8ePHxf//fefyMnJET179hSVK1cWJ06cEImJidpXVlaWVZ7v085ZCCHat28vGjRoICIjI8W///4r1q5dKxwcHMSyZcus8pzZl22jLwthe/3Z1vqyEOWrP1tNQpSVlSXkcrnYsmWLTvnYsWNFu3btxKVLlwQAER0drbO9Z8+eYuDAgcW2O2DAANGzZ0+dsujoaAFA/Pvvv0IIIdq2bSvGjh2rU2fLli3Czs5OZGdnG3JaJfLFF1+IqlWrCoVCIZo1ayb27Nmjs33NmjWiZs2awsHBQTRp0qTQPS0GDRok2rdvr1N26tQp0bZtW6FUKoWvr6+YOXOmdl62wO7du8UzzzwjFAqFqFatmli+fLlJzu9RkZGRAkCh16BBg0R8fHyR2wCIyMhIqzxfIZ58zkIIkZiYKAYPHiz8/f2Fg4ODqFOnjli4cKFO/NZ0zuzLttGXhbC9/mxrfVmI8tWfrSYhEiJv4dbbb7+tU1avXj2dhVsff/yxdltWVlaJFm5VqFBB518k8+fPL7Rwq169ejr7jRw5kvd4INIT+zJR+VFe+rNVJURPG3KeP3++UKlUYsuWLSImJkb0799f+Pn5ibS0NG0b1jYcSVQesS8TlR/lpT9bVUIkxJOHnDUajZgxY4bw9fUVSqVStGvXTsTExOjsb23DkUTlFfsyUflRHvqzJMRjt30kIiIisjFWc9k9ERERkakwISIiIiKbx4SIiIiIbB4TIiIiIrJ5TIiIiIjI5llNQrRs2TIEBgbCwcEBQUFB2Ldvn3bbli1bEBoaCi8vL0iShBMnTpSozd27d0OSJNy9e9c0QRNRIcX15ZycHEyaNAmNGjWCs7Mz/P39MXDgQCQkJDy1TfZlIvN70u/lmTNnom7dunB2doa7uzs6deqEQ4cOPbVNS/Zlq0iINm3ahPHjx2Pq1Kk4fvw42rZti27dummfAn3v3j08++yzmD9/voUjJaIneVJfvn//PqKjozF9+nRER0djy5YtiIuLQ8+ePS0dNhE95mm/l2vXro2lS5ciJiYG+/fvR7Vq1dClSxfcvHnTwpE/gd53MDKjFi1aiJEjR+qU1a1bV0yePFmnrODZOMePHy9RuwXPnblz544QQoiUlBTRr18/UalSJeHo6CgaNmwoNmzYoLNP+/btxZgxY8TEiROFu7u78PHxETNmzND31IhsSkn7coHDhw8LANqHYxaHfZnIvErbl1NTUwUA8ffffz+xXUv25TI/QpSdnY1jx46hS5cuOuVdunRBVFSUUY+VmZmJoKAgbNu2DbGxsRg+fDgGDBhQaJhv/fr1cHZ2xqFDh7BgwQLMnj0bERERRo2FqLzRpy+npqZCkiRUqFChVMdiXyYyndL25ezsbKxcuRIqlQpNmjQp1bHM2ZfLfEKUkpICtVoNHx8fnXIfHx8kJSUZ9ViVKlVCWFgYmjZtiurVq2PMmDEIDQ3Fjz/+qFOvcePGmDFjBmrVqoWBAwciODgYO3fuNGosROVNaftyZmYmJk+ejNdeew1ubm6lOhb7MpHplLQvb9u2DS4uLnBwcMCiRYsQEREBLy+vUh3LnH25zCdEBSRJ0vkshChUVpyRI0fCxcVF+yqOWq3GRx99hMaNG8PT0xMuLi7YsWOHdk60QOPGjXU++/n5ITk5uYRnQmTbStKXc3Jy0K9fP2g0Gixbtkxbzr5MVHY8rS937NgRJ06cQFRUFLp27Yo+ffpo+1dZ7Mt2paptAV5eXpDL5YX+BZmcnFwoOy3O7NmzERYW9tR6CxcuxKJFi7B48WLtlS7jx49Hdna2Tj17e3udz5IkQaPRlCgWIltV0r6ck5ODPn36ID4+Hrt27dIZHWJfJrK8kvZlZ2dn1KxZEzVr1kSrVq1Qq1YtrFmzBlOmTCmTfbnMJ0QKhQJBQUGIiIjAyy+/rC2PiIjASy+9VKI2vL294e3t/dR6+/btw0svvYQ33ngDAKDRaHDhwgXUq1dPv+CJSKskfbkgGbpw4QIiIyPh6emp0wb7MpHl6ft7WQiBrKwsAGWzL5f5hAgAJkyYgAEDBiA4OBitW7fGypUrceXKFYwcORIAcPv2bVy5ckV7v5Lz588DAHx9feHr61vi49SsWRObN29GVFQU3N3dER4ejqSkJP4lSmQkT+rLubm5eOWVVxAdHY1t27ZBrVZr/wXq4eEBhUJR4uOwLxOZ1pP68r179/DRRx+hZ8+e8PPzw61bt7Bs2TJcu3YNr776aqmOY86+bBUJUd++fXHr1i3Mnj0biYmJaNiwIbZv346qVasCALZu3YohQ4Zo6/fr1w8AMGPGDMycObPYdguG0+zs8n4M06dPR3x8PEJDQ+Hk5IThw4ejV69eSE1NNdGZEdmWJ/Xly5cvY+vWrQCApk2b6uwXGRmJDh06FNsu+zKReT2pL2dmZuLcuXNYv349UlJS4OnpiebNm2Pfvn1o0KDBE9u1ZF+WhBDC6K1aiY0bN2LYsGHIyMiwdChEZAD2ZaLywZJ92SpGiIwtKysLly5dwtKlS9GpUydLh0NEemJfJiofykJftprL7o3pjz/+QMuWLeHs7IzPPvvM0uEQkZ7Yl4nKh7LQl216yoyIiIgIsNERIiIiIqJHMSEiIiIim2d1CdG8efPQvHlzuLq6wtvbG7169dLed6iAEAIzZ86Ev78/HB0d0aFDB5w+fVqnTlZWFsaMGQMvLy84OzujZ8+euHbtmk6dO3fuYMCAAVCpVFCpVBgwYADu3r1r6lMkIiIiM7O6hGjPnj145513cPDgQURERCA3NxddunTBvXv3tHUWLFiA8PBwLF26FEeOHIGvry86d+6M9PR0bZ3x48fj559/xsaNG7F//35kZGSge/fuUKvV2jqvvfYaTpw4gT///BN//vknTpw4gQEDBpj1fImIiMj0rH5R9c2bN+Ht7Y09e/agXbt2EELA398f48ePx6RJkwDkjQb5+Pjg448/xogRI5CamoqKFSvim2++Qd++fQEACQkJCAgIwPbt2xEaGoqzZ8+ifv36OHjwIFq2bAkAOHjwIFq3bo1z586hTp06FjtnIiIiMi6rGyF6XMHdKj08PAAA8fHxSEpKQpcuXbR1lEol2rdvj6ioKADAsWPHkJOTo1PH398fDRs21NY5cOAAVCqVNhkCgFatWkGlUmnrEBERUflg1QmREAITJkxAmzZt0LBhQwDQPvvo0SfuFnwu2JaUlASFQgF3d/cn1inqwXPe3t6FnvBLRERE1s2q71Q9evRonDp1Cvv37y+0TZIknc9CiEJlj3u8TlH1S9IOERERWRerHSEaM2YMtm7disjISFSuXFlbXvB0+8dHcZKTk7WjRr6+vsjOzsadO3eeWOfGjRuFjnvz5s1Co09ERERk3awuIRJCYPTo0diyZQt27dqFwMBAne2BgYHw9fVFRESEtiw7Oxt79uxBSEgIACAoKAj29vY6dRITExEbG6ut07p1a6SmpuLw4cPaOocOHUJqaqq2DhEREZUPVneV2ahRo7Bhwwb8+uuvOld6qVQqODo6AgA+/vhjzJs3D2vXrkWtWrUwd+5c7N69G+fPn4erqysA4O2338a2bduwbt06eHh4ICwsDLdu3cKxY8cgl8sBAN26dUNCQgJWrFgBABg+fDiqVq2K3377zcxnTURERKZkdQlRcet31q5di8GDBwPIG0WaNWsWVqxYgTt37qBly5b44osvtAuvASAzMxMTJ07Ehg0b8ODBAzz//PNYtmwZAgICtHVu376NsWPHYuvWrQCAnj17YunSpahQoYLJzo+IiIjMz+oSIiIiIiJjs7o1RERERETGxoSIiIiIbB4TIiIiIrJ5TIiIiIjI5jEhIiIiIpvHhIiIiIhsHhMiIiIisnlMiIiIiMjmMSEiIqsyc+ZMNG3a1NJhEFE5wztVE1GZUdyjeQoMGjQIS5cuRVZWFjw9Pc0UFRHZAiZERFRmJCUlad9v2rQJH374Ic6fP68tc3R0hEqlskRoRFTOccqMiMoMX19f7UulUkGSpEJlj0+ZDR48GL169cLcuXPh4+ODChUqYNasWcjNzcXEiRPh4eGBypUr46uvvtI51vXr19G3b1+4u7vD09MTL730Ei5fvmzeEyaiMoMJERFZvV27diEhIQF79+5FeHg4Zs6cie7du8Pd3R2HDh3CyJEjMXLkSFy9ehUAcP/+fXTs2BEuLi7Yu3cv9u/fDxcXF3Tt2hXZ2dkWPhsisgQmRERk9Tw8PPDZZ5+hTp06ePPNN1GnTh3cv38fH3zwAWrVqoUpU6ZAoVDgn3/+AQBs3LgRMpkMq1evRqNGjVCvXj2sXbsWV65cwe7duy17MkRkEXaWDoCIyFANGjSATPbw33c+Pj5o2LCh9rNcLoenpyeSk5MBAMeOHcPFixfh6uqq005mZiYuXbpknqCJqExhQkREVs/e3l7nsyRJRZZpNBoAgEajQVBQEL777rtCbVWsWNF0gRJRmcWEiIhsTrNmzbBp0yZ4e3vDzc3N0uEQURnANUREZHNef/11eHl54aWXXsK+ffsQHx+PPXv2YNy4cbh27ZqlwyMiC2BCREQ2x8nJCXv37kWVKlXQu3dv1KtXD2+++SYePHjAESMiG8UbMxIREZHN4wgRERER2TwmRERERGTzmBARERGRzWNCRERERDaPCRERERHZPCZEREREZPOYEBEREZHNY0JERERENo8JEREREdk8JkRERERk85gQERERkc37f7OSZyCYAZwRAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# plot results\n", "y = pd.concat([y_exp, y_imp], axis=1, keys=['Explicit', 'Implicit'])\n", "# Flatten the two-level column labels into a single level\n", "y.columns = y.columns.get_level_values(0)\n", "ax = y.plot()\n", "ax.set_xlabel('Time')\n", "ax.set_ylabel('Indoor temperature, $\\\\theta_i$ / °C')\n", "ax.set_title(f'Time step: $dt$ = {dt:.0f} s; $dt_{{max}}$ = {Δtmax:.0f} s')\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "89360f32-edb1-4931-9408-2a92cbed8fc3", "metadata": {}, "source": [ "> Figure 3. Step response to auxiliary heat gains $\\dot Q_a$ by using Euler\n", "[implicit](https://en.m.wikipedia.org/wiki/Numerical_methods_for_ordinary_differential_equations#Backward_Euler_method)\n", "and\n", "[explicit](https://en.m.wikipedia.org/wiki/Numerical_methods_for_ordinary_differential_equations#Euler_method) integration." ] }, { "cell_type": "code", "execution_count": 28, "id": "e69e61a0-e3d8-46a4-a786-50e89d6144f3", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Steady-state indoor temperature obtained with:\n", "- DAE model: 12.2566 °C\n", "- state-space model: 12.2566 °C\n", "- steady-state response to step input: 12.2549 °C\n" ] } ], "source": [ "print('Steady-state indoor temperature obtained with:')\n", "print(f'- DAE model: {float(θssQ[6]):.4f} °C')\n", "print(f'- state-space model: {float(yssQ):.4f} °C')\n", "print(f'- steady-state response to step input: \\\n", "{y_exp[\"θ6\"].tail(1).values[0]:.4f} °C')" ] }, { "cell_type": "markdown", "id": "0ebf5b58-0cf9-4c35-b751-09382cfd4510", "metadata": { "jp-MarkdownHeadingCollapsed": true, "tags": [] }, "source": [ "## Discussion\n", "\n", "### Time step\n", "\n", "The time step depends on:\n", "- Capacities considered into the model:\n", " - if the capacities of the air $C_a =$ `C['Air']` and of the glass $C_g =$ `C['Glass']` are considered, then the time step is small;\n", " - if the capacities of the air and of the glass are zero, then the time step is large (and the order of the state-space model is reduced).\n", " \n", "- P-controller gain `Kp`:\n", " - if $K_p \\rightarrow \\infty$, then the controller is perfect and the time step needs to be small;\n", " - if $K_p \\rightarrow 0$, then, the controller is ineffective and the building is in free-running.\n", "\n", "#### Air and glass capacities\n", "\n", "In cell `[2]`, consider:\n", "\n", "```\n", "controller = False\n", "neglect_air_glass_capacity = False\n", "imposed_time_step = False\n", "```\n", "\n", "Note that the maximum time step for numerical stability of explicit Euler integration method is: `dtmax = 499 s = 8.3 min`.\n", "\n", "Now, in cell `[2]`, consider:\n", "\n", "```\n", "controller = False\n", "neglect_air_glass_capacity = False\n", "imposed_time_step = True\n", "Δt = 498 # s, imposed time step\n", "```\n", "\n", "and explain the simulation results.\n", "\n", "Determine the number of state variables, $\\theta_s$.\n", "\n", "In cell `[2]`, consider:\n", "\n", "```\n", "controller = False\n", "neglect_air_glass_capacity = True\n", "imposed_time_step = False\n", "```\n", "\n", "Note that the maximum time step for numerical stability of explicit Euler integration method is: `dtmax = 9587 s = 2.7 h`. Observe the simulation results, comment and correct them.\n", "\n", "Determine the number of state variables, $\\theta_s$ and compare with the case in which the capacities of the air and of the glass were considered.\n", "\n", "#### Controller gain\n", "\n", "The controller models an HVAC system able to maintain setpoint temperature by heating (when $q_{HVAC} > 0$) and by cooling (when $q_{HVAC} < 0$).\n", "\n", "In cell `[2]`, consider:\n", "\n", "```\n", "controller = True\n", "neglect_air_glass_capacity = True\n", "imposed_time_step = False\n", "```\n", "\n", "Maximum time step for numerical stability of explicit Euler integration method is: `dtmax = 8441 s = 2.3 h`. \n", "\n", "Observe, comment and improve the simulation results.\n", "\n", "In cell `[2]`, consider:\n", "\n", "```\n", "controller = True\n", "neglect_air_glass_capacity = False\n", "imposed_time_step = False\n", "```\n", "\n", "Comment on time step and simulation results.\n", "\n", "In cell `[4]`:\n", "```\n", "if controller:\n", " TC['G']['q11'] = 1e3 # Kp -> ∞, almost perfect controller\n", "\n", "```\n", "\n", "Change the gain of the controller to `1e5` in cell `[4]`:\n", "```\n", "if controller:\n", " TC['G']['q11'] = 1e5 # Kp -> ∞, almost perfect controller\n", "```\n", "Comment the error raised in cell `[16]`." ] }, { "cell_type": "markdown", "id": "08d18237-b075-4f6c-82df-fa70d47a5216", "metadata": { "tags": [] }, "source": [ "## References\n", "\n", "1. [C. Ghiaus (2013)](https://doi.org/10.1016/j.energy.2012.10.024). Causality issue in the heat balance method for calculating the design heating and cooling loads, *Energy* 50: 292-301, , open access preprint: [HAL-03605823](https://hal.archives-ouvertes.fr/hal-03605823/document)\n", "\n", "2. [C. Ghiaus (2021)](https://doi.org/10.1007/978-3-030-76477-7_5). Dynamic Models for Energy Control of Smart Homes, in *S. Ploix M. Amayri, N. Bouguila (eds.) Towards Energy Smart Homes*, Online ISBN: 978-3-030-76477-7, Print ISBN: 978-3-030-76476-0, Springer, pp. 163-198, open access preprint: [HAL 03578578](https://hal.archives-ouvertes.fr/hal-03578578/document)\n", "\n", "3. [J.A. Duffie, W. A. Beckman, N. Blair (2020)](https://www.eng.uc.edu/~beaucag/Classes/SolarPowerForAfrica/Solar%20Engineering%20of%20Thermal%20Processes,%20Photovoltaics%20and%20Wind.pdf). Solar Engineering of Thermal Processes, 5th ed. John Wiley & Sons, Inc. ISBN 9781119540281\n", "\n", "4. [Réglementation Thermique 2005. Méthode de calcul Th-CE.](https://pdfslide.fr/documents/rt2005-methode-de-calcul-th-ce.html). Annexe à l’arrêté du 19 juillet 2006\n", "\n", "5. H. Recknagel, E. Sprenger, E.-R. Schramek (2013) Génie climatique, 5e edition, Dunod, Paris. ISBN 978-2-10-070451-4\n", "\n", "6. [J.R. Howell et al. (2021)](http://www.thermalradiation.net/indexCat.html). Thermal Radiation Heat Transfer 7th edition, ISBN 978-0-367-34707-0, A Catalogue of Configuration Factors\n", "\n", "7. [J. Widén, J. Munkhammar (2019)](http://www.diva-portal.org/smash/get/diva2:1305017/FULLTEXT01.pdf). Solar Radiation Theory, Uppsala University" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.11.7" }, "toc-autonumbering": true, "toc-showcode": false, "toc-showmarkdowntxt": false, "toc-showtags": false }, "nbformat": 4, "nbformat_minor": 5 }