{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 0,
     "link": "[](https://jckantor.github.io/cbe30338-2021/04.04-State-Estimation-Copy1.html)",
     "section": ""
    }
   },
   "source": [
    "<!--NOTEBOOK_HEADER-->\n",
    "*This notebook contains material from [cbe30338-2021](https://jckantor.github.io/cbe30338-2021);\n",
    "content is available [on Github](https://github.com/jckantor/cbe30338-2021.git).*\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 0,
     "link": "[](https://jckantor.github.io/cbe30338-2021/04.04-State-Estimation-Copy1.html)",
     "section": ""
    }
   },
   "source": [
    "<!--NAVIGATION-->\n",
    "< [4.4 Anomaly Detection](https://jckantor.github.io/cbe30338-2021/04.04-Anomaly-Detection.html) | [Contents](toc.html) | [Tag Index](tag_index.html) | [4.6 Lab Assignment 6: Anomaly Detection](https://jckantor.github.io/cbe30338-2021/04.06-Lab-Assignment-Anomaly-Detection.html) ><p><a href=\"https://colab.research.google.com/github/jckantor/cbe30338-2021/blob/master/docs/04.04-State-Estimation-Copy1.ipynb\"> <img align=\"left\" src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open in Colab\" title=\"Open in Google Colaboratory\"></a><p><a href=\"https://jckantor.github.io/cbe30338-2021/04.04-State-Estimation-Copy1.ipynb\"> <img align=\"left\" src=\"https://img.shields.io/badge/Github-Download-blue.svg\" alt=\"Download\" title=\"Download Notebook\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 1,
     "link": "[4.4 Estimation and Anamoly Detection](https://jckantor.github.io/cbe30338-2021/04.04-State-Estimation-Copy1.html#4.4-Estimation-and-Anamoly-Detection)",
     "section": "4.4 Estimation and Anamoly Detection"
    }
   },
   "source": [
    "# 4.4 Estimation and Anamoly Detection\n",
    "\n",
    "This is an advanced topic ... \n",
    "\n",
    "This notebook provides what you will need to implement one meaans for controller diagnositics."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 1,
     "link": "[4.4 Estimation and Anamoly Detection](https://jckantor.github.io/cbe30338-2021/04.04-State-Estimation-Copy1.html#4.4-Estimation-and-Anamoly-Detection)",
     "section": "4.4 Estimation and Anamoly Detection"
    }
   },
   "source": [
    "![](figures/FeedbackControlDiagram2.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 1,
     "link": "[4.4 Estimation and Anamoly Detection](https://jckantor.github.io/cbe30338-2021/04.04-State-Estimation-Copy1.html#4.4-Estimation-and-Anamoly-Detection)",
     "section": "4.4 Estimation and Anamoly Detection"
    }
   },
   "source": [
    "In the context of control, the control variable (CV) is heater temperature, the measured process variable (PV) is the sensor temperature, and we haven't yet established the relationship between them."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[4.4.1 State-Space Model](https://jckantor.github.io/cbe30338-2021/04.04-State-Estimation-Copy1.html#4.4.1-State-Space-Model)",
     "section": "4.4.1 State-Space Model"
    }
   },
   "source": [
    "## 4.4.1 State-Space Model\n",
    "\n",
    "The same model can be written in a **state-space** matrix/vector format:\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{dx}{dt} & = A x + B_u u + B_d d \\\\\n",
    "y & = C x\n",
    "\\end{align}\n",
    "\n",
    "where\n",
    "\n",
    "$$x = \\begin{bmatrix} T_{H,1} \\\\ T_{S,1} \\end{bmatrix}\n",
    "\\qquad\n",
    "u = \\begin{bmatrix} u_1 \\end{bmatrix}\n",
    "\\qquad\n",
    "d = \\begin{bmatrix} T_{amb} \\end{bmatrix}\n",
    "\\qquad\n",
    "y = \\begin{bmatrix} T_{S,1} \\end{bmatrix}$$\n",
    "\n",
    "and\n",
    "\n",
    "$$A = \\begin{bmatrix} -\\frac{U_a+U_b}{C^H_p} & \\frac{U_b}{C^H_p} \\\\ \\frac{U_b}{C^S_p} & -\\frac{U_b}{C^S_p} \\end{bmatrix}\n",
    "\\qquad\n",
    "B_u = \\begin{bmatrix} \\frac{\\alpha P_1}{C^H_p} \\\\ 0 \\end{bmatrix}\n",
    "\\qquad\n",
    "B_d = \\begin{bmatrix} \\frac{U_a}{C_p^H} \\\\ 0 \\end{bmatrix}\n",
    "\\qquad\n",
    "C = \\begin{bmatrix} 0 & 1 \\end{bmatrix}$$\n",
    "\n",
    "By common convention, in state-space models the vector $x$ contains all variables representing the **state** of the system. The input vector $u$ is reserved for all manipulated variables (MV's), and the input vector $d$ for all disturbance variables (DV's). The vector $y$ holds all measured process variables (PV's)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[4.4.1 State-Space Model](https://jckantor.github.io/cbe30338-2021/04.04-State-Estimation-Copy1.html#4.4.1-State-Space-Model)",
     "section": "4.4.1 State-Space Model"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "# parameter estimates.\n",
    "alpha = 0.00016       # watts / (units P * percent U1)\n",
    "P1 = 200              # P units\n",
    "P2 = 100              # P units\n",
    "Ua = 0.050            # heat transfer coefficient from heater to environment\n",
    "CpH = 2.2             # heat capacity of the heater (J/deg C)\n",
    "CpS = 1.9             # heat capacity of the sensor (J/deg C)\n",
    "Ub = 0.021            # heat transfer coefficient from heater to sensor\n",
    "Tamb = 21             # ambient room temperature\n",
    "\n",
    "# state space model\n",
    "A = np.array([[-(Ua + Ub)/CpH, Ub/CpH], [Ub/CpS, -Ub/CpS]])\n",
    "Bu = np.array([[alpha*P1/CpH], [0]])     # single column\n",
    "Bd = np.array([[Ua/CpH], [0]])           # single column\n",
    "C = np.array([[0, 1]])                   # single row\n",
    "\n",
    "u_initial = np.array([0])                # initial manipulable input\n",
    "d_initial = [Tamb]                       # initial disturbance input\n",
    "x_initial = np.array([Tamb, Tamb])       # initial steady state"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 3,
     "link": "[4.4.1.1 Prediction and Correction](https://jckantor.github.io/cbe30338-2021/04.04-State-Estimation-Copy1.html#4.4.1.1-Prediction-and-Correction)",
     "section": "4.4.1.1 Prediction and Correction"
    }
   },
   "source": [
    "### 4.4.1.1 Prediction and Correction\n",
    "\n",
    "Making the substitution $\\hat{y} = C\\hat{x}$, the observer equation becomes\n",
    "\n",
    "\\begin{align*}\n",
    "\\frac{d\\hat{x}}{dt} & = \\underbrace{A \\hat{x} + B_u u + B_d \\hat{d}}_{\\text{prediction}} - \\underbrace{L(\\hat{y} - y)}_{\\text{correction}}\n",
    "\\end{align*}\n",
    "\n",
    "As before, we can try approximating the derivative with a first order difference formula over small intervals of time. One clever detail is to choose the start of the interval, time $t_{k-1}$, to evaluate the model prediction part of the obsever equation, and use the process measurement available at end of the interval, time $t_k$, to evaluate the correction term.\n",
    "\n",
    "$$\\frac{\\hat{x}_{k} - \\hat{x}_{k-1}}{t_k - t_{k-1}} = \\underbrace{A \\hat{x}_{k-1} + B_u u_{k-1} + B_d \\hat{d}_{k-1}}_{\\text{prediction at time $t_{k-1}$}} - \\underbrace{L(\\hat{y}_k - y_k)}_{\\text{correction at time $t_k$}}$$\n",
    "\n",
    "For this purpose we introduce some extra notation. We use $\\hat{x}_{k-1}$ to represent the state estimate at time $t_{k-1}$ using all information up to time $t_{k-1}$. Given knowledge of the manipulated inputs $u_{k-1}$ and an estimate of the disturbances $\\hat{d}_{k-1}$, we will use the model to create a predicted state $\\hat{x}_k^{pred}$ for time $t_k$. A measurement is then taken at time $t_k$ which is used to produce a new estimate $\\hat{x}_k$. \n",
    "\n",
    "At each time step $t_k$ there are two calculations to perform:\n",
    "\n",
    "* **Model Prediction:** Use the model to update the state to the next time step, i.e., $\\hat{x}_{k-1} \\rightarrow \\hat{x}_{k}^{pred}$ with the equation\n",
    "\n",
    "\\begin{align}\n",
    "\\hat{x}_k^{pred} & = \\hat{x}_{k-1} + (t_k - t_{k-1}) ( A \\hat{x}_{k-1} + B_u u_{k-1} + B_d \\hat{d}_{k-1}) \\\\\n",
    "\\hat{y}_k^{pred} & = C \\hat{x}_k^{pred}\n",
    "\\end{align}\n",
    "\n",
    "* **Measurement Correction:** Use measurement $y_k$ to update $\\hat{x}_{k}^{pred} \\rightarrow \\hat{x}_{k}$ with the equation\n",
    "\n",
    "$$\\hat{x}_{k} = \\hat{x}_{k}^{pred} - (t_k - t_{k-1})L (\\hat{y}_{k}^{pred} - y_k)$$ \n",
    "\n",
    "While the notation may be a little daunting, the underlying concepts are quite straightforward. At each time step we have two calculations to perform. First, use all available information to predict the current state using the process model. Second, when a measurement becomes available, compare the measure to the model prediction and apply an appropriate correction to the state estimate. \n",
    "\n",
    "Let's see how this works in practice. The next cell introduces `tclab_observer`, a Python generator that creates instances of observers using the two-state model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "nbpages": {
     "level": 3,
     "link": "[4.4.1.1 Prediction and Correction](https://jckantor.github.io/cbe30338-2021/04.04-State-Estimation-Copy1.html#4.4.1.1-Prediction-and-Correction)",
     "section": "4.4.1.1 Prediction and Correction"
    }
   },
   "outputs": [],
   "source": [
    "def tclab_observer(L, x_initial=[Tamb, Tamb]):\n",
    "    # initialize variables\n",
    "    t_now = 0\n",
    "    x_now = x_initial\n",
    "    d = d_initial                         # estimate of unmeasured disturbance\n",
    "    y_err = None\n",
    "    \n",
    "    while True:\n",
    "        # yield current state, get MV for next period\n",
    "        t_next, Q, T_measured = yield x_now, y_err\n",
    "        u = np.array([Q])\n",
    "        d = np.array([Tamb])\n",
    "        y = np.array([T_measured])\n",
    "        \n",
    "        # model prediction\n",
    "        x_predict = x_now + (t_next - t_now)*(np.dot(A, x_now) + np.dot(Bu, u) + np.dot(Bd, d))\n",
    "        y_predict = np.dot(C, x_predict)\n",
    "        \n",
    "        # measurement correction\n",
    "        y_err = y_predict - y\n",
    "        x_correct = x_predict - (t_next - t_now)*np.dot(L, y_err)\n",
    "        \n",
    "        # update time and state\n",
    "        t_now = t_next\n",
    "        x_now = x_correct"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[4.4.2 Testing](https://jckantor.github.io/cbe30338-2021/04.04-State-Estimation-Copy1.html#4.4.2-Testing)",
     "section": "4.4.2 Testing"
    }
   },
   "source": [
    "## 4.4.2 Testing\n",
    "\n",
    "The first experiment is to test if the observer can converge to the measured temperature. We will start with a very bad initial estimate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[4.4.2 Testing](https://jckantor.github.io/cbe30338-2021/04.04-State-Estimation-Copy1.html#4.4.2-Testing)",
     "section": "4.4.2 Testing"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm0AAAFgCAYAAADpfU79AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAxOAAAMTgF/d4wjAAA1x0lEQVR4nO3de7xcdX3v//dnTy47OyHsZGcHJAGDCijlaLikQrVqRElq668qHlCB2kgrFoqnJ8rBg54e66nKKZr2USRHqci9CpHUI17YHjWIYilBiFxaCQibkAAht51NLvs28/n9MWvCzJ77zJpZa+15PR+PPLJnrZm1PvPJmtnvfNfN3F0AAACIt66oCwAAAEB1hDYAAIAEILQBAAAkAKENAAAgAQhtAAAACUBoAwAASABCGwAAQAJMi7qAVpk+fbofccQRUZcRG6Ojo5o5c2bUZcQKPSlGTwrRj2L0pBD9KEZPim3btm3M3ZtuSmxCm5kNShoJ/kjSF939NjNbKOkmSa+WNCrpY+7+i2rL6+3t1datW1tVbuIMDAxoxYoVUZcRK/SkGD0pRD+K0ZNC9KMYPSlmZjvCWE5sQlvg/e7+6KRpV0q6z91XmtkySd82s1e7+0QE9QEAAEQibqGtlHMkHStJ7r7RzLZLerOkuyu9yD2tXTs317SC3t4lSk2b0WSZAAAArRO30HarmXVJ+jdJ/11SRlKXu+cPKw5KOqbagg6kX9LFd55T00qPmjZHX/rgTwluwBSSzqQ1NDpUNL13Zq9SXamaXvdS+iXtOrir7mUAQCtYXG4Yb2bHuPsWM5su6W8l/SdJF0ja4u6z8563TtKd7n7TpNevlrQ697i7d9qiP7zqtTWv/08WrtbMmUc1+S7ia2RkRN3d3VGXESv0pFgzPcl4Rvsz+0OuqDEZZXTL0C0a9/GiedNtus7vPV9dJU6en/y6jGfUZcXPm5eapwt6Lyg5b6rjc1OIfhSjJ8VWrly5zd0XN7uc2IS2fGb2Ckmb3f0wM9svaUlutM3M7pf039z97krLWLBgvj/+m/sqrmdo7zO64p7LJElr3327+hYcH0b5scSBocXi2JNyo0PtsmHDBi1fvrzu16U9rct+dplG0iPVn5wgw3uHNffwuSXnrT1zrfpm9bW5oujF8XMTJfpRjJ4UM7NQQlssdo+a2WxJ0919KJj0QUkPBT+vk3SJpM8GJyIcKanq2aNmqSkdwhCNVoaqOASf4b3DWveTdZGtP2zdqW5d9darlLJUXf3Nve6eu+8pCLFDo0O64hdXHPqZ3aQA2ikWoU3SEZLuMLOUJJP0lKQ/CeZdLulmM3tC0pikCzhzFPUoF7TKHa9UdjkxCFVxlx+S4mByqPrGim/UFLpzrzssdVjZ0bQrfnFFwfslwAFotViENnd/StLJZeZtl3RWeytCXDQ7slUpaMV1VCnK4NPo7tGcuAeXVFeqqV2avTN7ddTso/Tc/uckSSPpEV3600sliQAHoOViEdrQuSqFsriObLU6VEX5C7/SyBKyoe9/v+lKPf38b/S5X1+pkfTooXn7xsf1FwMfkyR1p2bqr9/wKaUspbkzDqu4rfT2HanUNL6KAVTHNwVaKi6hrFTQanRUiVGUZEhPTGho1wtNLWN0/5B2bX/5ziqZ9IR2XH+epmdG9D/keqlLypi0pm9Co/by6/ZLuvyFCyRJM11avWuauvLO+TosI6WUfcH2rm71r7pVXanGvo57e2YoZcHKe+ZLbJvAlEVoQ8Oq7boMK5SFMbJVKmgxqhQ/YQQtqTBcNeP4TFq7HyvcbqYHf6dk6s1kf/6bF6eVDXCjJn1xQeFhuP1p6RM7pykl0/TMiIauO7vhGofMdExfT1DcLL30R9fK6/ysMNoHJAOfUtRkckALc5SsWihjZCsZmg1cYQWtnOnVn9Kw8RKjY/3B31d7WsNjLynt6aJdqDm7Jb105l9p4o7LNW9i9NCoWyMy7hrcmbs+3n7pn95X9zLqGe3LH30k7AHtxacNBUqNnjUb0AhlyZCemCjaHVirsAJX2EGrVLiqx7333qs3velNRdOrhZWFwd83HX1Lwecp/5IhX9p8jXTS0ZrmqUPHv+VUOw5OktLuuuKORzUykVZKaV320v/WTC8OiLWoZ7Qvf/SxnrBHwAOaxycIh6QzaX3yZ588dGZcLWrZdUkoa72wRrmOn9hftDuwVmEFrmaDVr5mg8LM2b3qO6Lx62FOPlt18tmnMtOEZfTXj3yh4HW1non6tb9YrD0HxrIPMu+UHdxdV33Nhu16wl7Tx+4R+gBCW6epdBza0OhQxcBWKqARyEKWSUsHsr940+4ayv1CrviS8Ea50k0tIZzANZV/Oae6UvrSW7+kodGhiiPY5S4lki/32VswZ+bLE+f21F3T/Mt/Vlfgv/fee3XG6W+se5tr9ti9ZkOfNLW3LXQGtt4OkM6k9VL6Jb144MWad3N+4c1fUO/M3oJpBLTG1DoKZp7WYd/7qDR+UJK0ZdcBZWq8zVxYo1xjmqlXXHgboyEtlD/6Nvliv6WCXH6Ay5cLc33dfU19LlPTptU1mjhzdq/6j1pSc9gL7T8VTYY+qfbgV3BGbj7OzkXE+Had4nK7PH+z5zda99PaLiR71OyjtGTuEgJaCfXuhqz3F1bt92coFsYo1/0PPqzXH7WkiSpQj1IX+80FuWrHkubCXLWRuJbVXkfYq3c0L1+YJ6jUGvwKzsgtWEDx2bn1HAfKf2rQLLaeKWby7s9Suzw5MWCSvF2SOaV2TTb6y6PeUbBRm6mrDrtcaaXUPS2lL5x9Uun/9U8Sxi+Erq7Hmno9mldpJE4qHo2rNhIXh895vaN5kzUT+qT6P7uFZ+TmKz47t9RlYcphFy+axb/8FJALatX+Z57b5dlpoazc6Njo/iHtfuGZgl2SOeV2TTa6G7KeUTCfNV9fDv595vXMUKqr8ctBINnK3XbrGyu+oV0ju2oaiSsnTqGummZDn1Rb8Ms/I7dg/U2enZvTzl28tSAAJg//WglX6xmf81Lzkr/LM+SD9I/PpLXrsVTduyQb2Q3JlyPClOpKaWHPwppG4sppNNRlPNNw3VGqNfgVnJGbr8zZueUuC1Pw0gh28dYizACYj2v5tQ6dTJhGd3/e/7P7ExPYSo2Mtfog/fxdkjnldk3yJYS4qDQS1+zt48qFuoPDB3XqgVPrvkNJnEbuKkl1WeEZuflKnJ1b62Vh2r2LtxZhBsB8jV7Lrxqf9fKJIJ26F4LfPAkQxu7PLutqR6ll1XoAf6UvprAP0s//H3L+LsmcTv1SQPKVC3M5zYS6cR+vOEJXTrO3o0tK6CunXbt4a9GKAFhOmMGw0eN965E7c7jc3pyo/9NOaIu5Wnd/RnHGZxhBrJRKI2NhHqTf7IVTgaRqJNTlwtywhhtaZ7XdsdWEcQ/iUpIUBsMIfjlhBcBSGr2WXzUzfVSfGf7cocd7/ym0RR8yZKaj58/Ss7sPltybU8vIYSuDHaEt5mrd/dmKL55KoSzMIFZKuePGOEgfaL1Ku13/5a5/0fLly2teVlj3KW429JXTbBh8Kf2Sdh1sZj9AoXaFyDAD4GT1XsuvmmYOj6lXxl3P7DpQdn5NI4cf/b8t623Noc3M3ijp1fmvcfebWlFUpyp1t4L8x2Gc/ZnOeN5tb9Jlb3tTSygLK4iVEvUQNIBiqa6UDksdVnGUrpRKu2OrCSv0ldNsGBzeO6x1P6ntGpi1aNWIYqtU+n0UajD8yA8OnYi2pMYT0epR6szh/L057dylXElNvxXN7P9IWiFpk16+041LIrSFpJbdoL0ze6t+WVa6vMWLLzxb1w2mawllBDEA1VTbHVtNM6GvnFaHwUa1akSxVUqFzLBHHw9J5R2bPXdWuIuW9MU/Xaahg+OHpvXOmq6uYG9Ol6T+//Iv2rv7xarLOmxef6i15av1N+g7JJ3o7vHauqeIdCatweHB8oHNXQumz1d67z7tGj5Y+jmqfnmLvY+ldHmdtVULZQQxAK3WbOgrJ4wwuGHDhrp2F5cT1xBZTamQGfboY9Ks7V+rPpU5A7lJtf62fT5pgS3jrp37Kl8IcejAmHK7xfccGJOXeX7u2KmCXYshSXtaf3Pf5dq+b9uhaVe8/hM6fPrhkqRMZkK7vvkXmjfxpPZufF+5xRxSy+hYV/4tWkrcliUfoQzAVBVGGGxkd3E5rRhRbJWkhsykq/jb2MzeFfz4SzO7XdK3JB36F3L3H7SwtqbsG5dWXb+x4nN6fKvGurMXivz0dx7VARsq+bxZ01P6hw8s1V99a5MOjhdeKbvL0zrMX2q4zozt1b5Z96grSI/9aWn2o5fKlR2SNUkLDv1Um2qXtyi4GXLPfM1PyJlTADCVtWpEsVXKhcywRh+Tqndmb8uWXW0I5bJJj/PHQF1SbENbmA6Op3XRzb8qmt7laV2+70r1p6vv4y5nqMv1+e5sYPv4rmk6akJKlQhozR47xuUtAABhKhcywxx9RKFqCeBv3P3udhQStjnTpetXLav4nKHdvbrix9kDGz//npPUO/+4gvnpjOvSf36oYHRt1vSUrv7QyUp1mWz/Ds1dv1/S7JpqSsu1V4WnKWfMNWPafsm69MoPXa/e7tIbOrspAQDobNVSwBpJp7SjkLB1WYVbkQRsZIZyewnn9cxQX4nnf/OjpxdcImOe7VPKgt2htk/KXSfsvV+Tesr/zyLtaf23+z6n5w5sL65DkrqmqXfhYv53AgAASqoW2jr+yqUpZbRAw1JmQvr2R6TxMhfd6+mT5hSf5pu79trQ6LCeG90lldm9edTso1q6HxwAACRbtdC2wMwuLjfT3deGXE+8ZNLSug9LQ89Wfl7v0VLP/KLJ5a69lrtIbsEiEnQrFQAA0H7VQluPpHIHhrXmHhJxkElnr7x8YFdxYJveI73/G1JXXut65kuTAle5a69FcY9QAACQfNVC2xZ3X9WWSuKi3Oha7pi1EgFtslIjbGHcggoAAHQuTkec7MDu4sDWe7S04LiqYS1n8k3eGV0DAADNqhbavtqWKuKqjtE1Kf+kg6FD077w5i8Q2AAAQNMqhjZ3/1q7ComcT0j7dmSPY8spc0ZoKeVOOmB3KAAACAO7R3Pu/C/SeOP3FZ28S1TiMh4AACA8hLac8RFJXS8/LnMZj8nK7RLlpAMAABAmQttkTZ4lKmV3iXJnAwAAEKau6k8pZmY/NrPvmtnvhV1Q5HLHsdUwQsYuUQAA0C6NjrT9d0mvlHSOpF+GV05ysUsUAID6ZDIZuSf/Wv1mpq6uhsbB6tJQaHP3jZI2Svp2uOUkF7tEAQCozdjYmLZs2aLx8fGoSwnN9OnTdcwxx2jGjBktW0dNoc3MPidpjaS9kr4n6Y2SLnL3O1pWWeH6j5N0o6QFkoYk/am7/3s71g0AAMK1ZcsWHXbYYerr65OZRV1O09xdu3bt0pYtW/Sa17ymZeupdaTtj939r83snZImJL1J0jcltSW0SfqapGvd/QYze7+k6ySd0fRSZ82TUjOkdHCpjxrPGAUAAI3JZDIaHx9XX1+fpk2bOudD9vX1affu3cpkMi3bVVprtzLB32+VtM7dH29XMjazhZJOkXRWMOkOSV8xsyXuPtjUwru6pHnHSpkJ6U1flOa/uuZbVQEAgPrljmGbCiNs+XLvp5XH6FWMgmb25eDH/Wb2KUkfkPT/zKxLUut22hY6WtJz7j4hSZ7txhZJx4SydJOUmibN7iOwAQDQgZYuXaqlS5fqxBNP1LRp0w49Pvfcc7Vv3z6tWLFCCxYs0IIFCyKt0yolQjN70N1PCY4pu0TSPe6+3sxeI+k/u/sXW16g2amSbnL338mbtlHSJ9z9nrxpqyWtzj3u6elZtH79+orLfin9kq7bc50k6cJ5F+qw1GEVn5/xjPZn9h96vN/361tD36r59VEaGRlRd3d31GXECj0pRk8K0Y9i9KQQ/ShWS0/6+/t17LHHtuWMy3o888wzWr58uZ566qlD00ZHR3Xfffdp3rx5es973lMwL18mk9HTTz+tHTt2FM1buXLlNndf3Gx9NYW2ZlfSjGD36BOS+tx9wrLjj89LOr3S7tH+/n4v1bh8uw7u0sU/uViStPbMtRXP/ix3Id2caq+P2sDAgFasWBF1GbFCT4rRk0L0oxg9KUQ/ilXrSTqd1ubNm3X88ccrlUopnXHtOdD4bSRrMa9nhlJd1XfHDg4O6rTTTtPOnTvrmicVv698ZhZKaKt2TNsJZnZ/uZnu/rvNFlCNu79oZg9JOl/SDZLOljTY9PFsdSp1Id0cLqgLAEBj9hwY06rrN7Z0HdevWqYFc2a2dB3tUC20PSfpsnYUUsVFkm4wsyskDUv6cJTF5C6km8MFdQEAQKtVC20vufvP2lJJBe7+uMK4xEdIuJAuAADhmNczQ9evWtbydUwF1ULb1DofFwAAxEqqy6bErst2qHbaxrvaUkWMpTNp7Tq4S0OjQ1GXAgAAInDKKafojDPO0J49e7R48WJdcMEFkdRRcaTN3Z9vVyFxVO2MUQAAMHUsWbKk5NmhDz74YATVFIvXBVJiptQZo5wpCgAAojB1bvrVYrkzRjlTFAAARKGhkTYzO83Mjgy7mDjLnTFKYAMAAFFodPfolZIeNrNrwywGAAAApTW0e9Td3xHcNP4NIdcDAACAEho+ps3dM5IeCrEWAAAAlNHw2aNmdmOYhcQJ12YDAABx08zZo8tDqyJGuDYbAACdZenSpZKksbExbd68WSeddJIk6YQTTtBtt90WYWWFKoY2M3ux3CxJvaFXEwNcmw0AgM6yadMmSdLg4KBOO+20Q4/jppZ7j54paW+J6fe2pKIY4dpsAAC0WCYtHdjd2nX0zJfq/D2+Y8cOnXfeeXr++edlZjr11FN1/fXXt6jA2lQLbb+S1OfuD0+eYWYvtKak+Mhdmw0AALTIgd3Sre9v7TrO+7Y0p7+ul9xyyy1asmSJfvSjH0mSdu9ucbCsQbUTEc5WmRE1dz81/HIAAACid/rpp+uuu+7SJz7xCX33u9/V7Nmzoy6p6kjb1939g22pBAAAdJ6e+dmRsFavo05nnHGGNm3apB//+Me644479JnPfEYPPfSQUqnoDpeqFtpe25YqAABAZ+pK1b3rsh2efvppLVq0SOecc45WrlyphQsXat++fTr88MMjq6laaPO2VAEAABAjd999t9asWaNUKqV0Oq2rrroq0sAmVQ9t/6nMZT9Mkrv7whbUBAAA0HZLlizRzp07JUmrVq3SqlWrIq6oULXQtlnSu9pRCAAAAMqrFtpG3f2ZtlQCAACAsqpd8sPaUkUMpJ37jQIAgPiqONLm7ie3q5CoXfazyzSSHom6DAAAgJKqjbR1jMmBjfuNAgCAOKl2TFvH4X6jAAAgjhhpmyR3v1ECGwAAnWHp0qVaunSpTjzxRE2bNu3Q43PPPVd33323TjvttKhLlMRIGwAA6HCbNm2SJA0ODuq000479FjKXmQ3LghtAAAgMulMuuVXbmj2kKeJiQldfPHFuvfeezUxMaEbb7wxktE3QhsAAIjM0OiQLv7JxS1dx9oz16pvVl/Dr3/sscf09a9/XWvXrtVXv/pVffrTn9bAwECIFdaGY9oAAAAqOOGEEw6NrJ1xxhn67W9/G0kdjLQBAIDI9M7s1doz17Z8Hc3o7u4+9HMqldLExESTFTWG0AYAACKT6ko1teuyk7B7FAAAIAEIbQAAAJKWLFminTt3Fkx729vepgceeODQ45NOOkmDg4NtriyL0AYAAJAAhDYAAIAEILQBAAAkQOShzcxuMLOtZrYp+HNV3rwuM7vazH5rZk+aWWuvvgcAAFrKzKIuoaVa+f7icsmPK939KyWmny/pREnHSzpc0oNm9lN3/01bqwMAAKEwM5mZxsfHlUo1fmupuBkfHz/03lolLqGtnHMlfdXd05J2m9ntkj4g6bORVgUAABpiZurt7dX27du1aNGiKTHy5u7avn27ent7OyK0rTazj0raIukz7r4pmH6MpGfynjcoqf13aAUAAKFZuHChnnnmGT3xxBNRlxKa7u5uLVy4sKXrMHdv7QrMfi7pdWVmnywpI+l5d8+Y2XslrZV0nLvvM7NHJH3E3TcGy7pE0qnu/pES61ktaXXucU9Pz6L169dXrC3jGd08dLP2pPdIkual5umC3gvUZZEf6he6kZGRgttwgJ6UQk8K0Y9i9KQQ/ShGT4qtXLlym7svbnY5LQ9t9TKzxyV9yN1/ZWbfl3SDu68L5v2dpAPu/tlqy+nv7/cdO3ZUXV86k9bQ6JCk7L3JUl1TZ/96voGBAa1YsSLqMmKFnhSjJ4XoRzF6Uoh+FKMnxcwslNAW+ZCSmS3O+/l0SX2SngwmrZN0kZmlzGy+sse43Rbm+nP3POub1TdlAxsAAEi+OBzTdoOZHSEpLemgpP/s7nuDeTdLWiZpc/D4Knf/jwhqBAAAiFTsdo+GZfr06X7EEUdEXUZsjI6OaubMmVGXESv0pBg9KUQ/itGTQvSjGD0ptm3bNrl706eVxmGkrSV6e3u1devWqMuIDY4xKEZPitGTQvSjGD0pRD+K0ZNiZpYOYzmRH9MGAACA6qbsSFsU0hnXngNjUZdR0vCYa+e+0ajLiBV6UoyeFCrXj3k9M5TqSv4FQQEkC6EtJOmM65JbH9S2oYNRl1LS3uG0bhrcGHUZsUJPitGTQuX6sah3lq457xSCG4Ba7QtjIeweDcmeA2OxDWwAwrVt6GBsR9UBxFIooY2RthZYc84bNG/2jKjLKLBhwwYtX74s6jJihZ4UoyeFJvdjz/4xrb791xFWBKCTEdpaYN7sGVowJ16nO8+dYbGrKWr0pBg9KUQ/AMQJu0cBAAASgNAGAACQAIQ2AACABCC0AQAAJAChDQAAIAEIbQAAAAlAaAMAAEgAQhsAAEACENoAAAASgNAGAACQAIQ2AACABIhVaDOz/2lmbmYnBY8XmtldZvaEmT1qZm+OukYAAIAoxCa0mdkpkk6XtCVv8pWS7nP34yStknSrmXGTewAA0HFiEdrMbKakayRdLMnzZp0TTJe7b5S0XRKjbQAAoOPEIrRJ+pykW9z96dwEM+uT1OXuO/KeNyjpmDbXBgAAEDlz9+rPamUBZmdI+rykM93dzWxQ0h9Jel7SFnefnffcdZLudPebSixntaTVucc9PT2L1q9f3+ryDxkec139cFqSdOnrU5o7w9q27lqMjIyou7s76jJihZ4UoyeFJvcj7p/zdmAbKUQ/itGTYitXrtzm7oubXU4cjg97q6TXSnrazCRpsaQBSX8mSWbWnzfa9koVHvN2iLuvkbQm97i/v99XrFjRwrIL7dw3qpsGN0qSli9fpgVzZrZt3bUYGBhQO/uRBPSkGD0pNLkfcf+ctwPbSCH6UYyetE7ku0fd/Up3P8rdl7j7EklbJa1w9x9KWifpEkkys2WSjpT0i8iKBQAAiEgcRtoquVzSzWb2hKQxSRe4+0TENQEAALRd7EJbMNqW+3m7pLOiqwYAACAeIt89CgAAgOoIbQAAAAlAaAMAAEgAQhsAAEACENoAAAASgNAGAACQAIQ2AACABCC0AQAAJAChDQAAIAEIbQAAAAlAaAMAAEgAQhsAAEACENoAAAASgNAGAACQAIQ2AACABCC0AQAAJAChDQAAIAEIbQAAAAnQ8tBmZje2eh0AAABTXTtG2pa3YR0AAABT2rQwFmJmL5abJak3jHUAAAB0slBCm7Lh7ExJe0tMvzekdQAAAHSssELbryT1ufvDk2eY2QshrQMAAKBjhRXazpY0Jklm1i/poLvvkyR3PzWkdQAAAHSsUE5EcPf9kv7czJ6TtF3SXjN72MzeIUlm1hvGegAAADpVKKHNzP5c0l9KulDSfEl9kj4l6ctmdpakn4SxHgAAgE4V1u7Rj0ta6e5b8qb9wMz+XdJmSWtCWg8AAEBHCus6bV2TApskyd0HJQ26+6dCWg8AAEBHCiu0zTCz7skTzWxWLesws24z+46ZbTazTWZ2l5ktCeYtDB4/YWaPmtmbQ6oZAAAgMcIKbesl3Zx/woGZzZN0k6Q7alzGtZJOcPelkr4XPJakKyXd5+7HSVol6VYzC2u3LgAAQCKEFdo+I2lc0lYze8jMHpT0rKSJYF5F7j7i7j9wdw8m3SfpVcHP50i6JnjeRmXPTmW0DQAAdJRQRqzcfVzSh8zs1ZJOCSY/5O5PNrjIj0u608z6lD1ebkfevEFJxzRcLAAAQALZy4Nb8WBmV0h6t7K3xZolaYu7z86bv07Sne5+06TXrZa0Ove4p6dn0fr169tTtKThMdfVD6clSZe+PqW5M6xt667FyMiIuruLDjvsaPSkGD0pNLkfcf+ctwPbSCH6UYyeFFu5cuU2d1/c7HJidWyYmX1S0vskvcPdD0g6YGYys/680bZXSip1puoa5V1apL+/31esWNGOsiVJO/eN6qbBjZKk5cuXacGcmW1bdy0GBgbUzn4kAT0pRk8KTe5H3D/n7cA2Uoh+FKMnrRPWMW1NC0bKPijpne4+lDdrnaRLgucsk3SkpF+0vUAAAIAIxWKkzcwWS/qypKckbTAzSRp19zdKulzZM1OfUPb+phe4+0RkxQIAAEQgFqHN3bdKKnlwiLtvl3RWeysCAACIl9jsHgUAAEB5hDYAAIAEILQBAAAkAKENAAAgAQhtAAAACUBoAwAASABCGwAAQAIQ2gAAABKA0AYAAJAAhDYAAIAEILQBAAAkAKENAAAgAQhtAAAACUBoAwAASABCGwAAQAIQ2gAAABKA0AYAAJAAhDYAAIAEILQBAAAkAKENAAAgAQhtAAAACUBoAwAASABCGwAAQAIQ2gAAABKA0AYAAJAAhDYAAIAEILQBAAAkAKENAAAgARIR2szsODP7pZltNrP7zezEqGsCAABop0SENklfk3Stux8v6e8kXRdxPQAAAG01LeoCqjGzhZJOkXRWMOkOSV8xsyXuPljudRl37dw32oYKs/bsH2vbugBEr1M/88Nj7f1ujTv6UazTezKvZ4ZSXdaSZcc+tEk6WtJz7j4hSe7uZrZF0jGSBsu9aN+4tOr6je2pEEDHWX37r6MuIRJ7h9O6aZDv1hz6UazTe3L9qmVaMGdmS5adhNAmST7pcVGENbPVklbnHk+bM097h/e2uq4ifd2mjb/YoC5rTcpu1MjIiAYGBqIuI1boSTF6UmhyPzLumjaW0a6RyV9JncMzHsl3a1zRj2Kd3pMNGzZo7ozWZABzj/eXT7B79AlJfe4+YWYm6XlJp1faPdq3YIE/PritTVW+rJXDos0YGBjQihUroi4jVuhJMXpSqFQ/0hnXngOduWtUyv5CWr58edRlxAb9KNbpPSmVA8xsm7svbnbZsR9pc/cXzewhSedLukHS2ZIGKwU2Seoya9nwJIDOlerq7O+WuTM6+/1PRj+K0ZPWiX1oC1wk6QYzu0LSsKQPR1wPAABAWyUitLn745LOiLoOAACAqCQitDViaGhIixc3vft4yhgdHdXMmQxX56MnxehJIfpRjJ4Uoh/F6ElJi8JYyJQNbb29vdq6dWvUZcQGB5gXoyfF6Ekh+lGMnhSiH8XoSTEzS4exnCkb2mqVfyZYXM/8BAAAmLKhrZY7IqQzrkv/+SEdHM8G4EW9s3TNeacQ3AAAQOxM2dDWyB0Rtg0d1J4DY5yqDAAAwrQvjIUk5YbxAAAASRVKaJuyI21zpmfv/1WLdMb1Zzc+0OKKAAAAGjdlQ1s9d0SoduwbAABA1Ng9CgAAkACENgAAgAQgtAEAACQAoQ0AACABCG0AAAAJQGgDAABIAEIbAABAAhDaAAAAEoDQBgAAkACENgAAgASIRWgzs24z+46ZbTazTWZ2l5ktCeYtDB4/YWaPmtmbIy4XAACg7WIR2gLXSjrB3ZdK+l7wWJKulHSfux8naZWkW81syt4zFQAAoJRYhDZ3H3H3H7i7B5Puk/Sq4OdzJF0TPG+jpO2SGG0DAAAdJRahrYSPS7rTzPokdbn7jrx5g5KOiaQqAACAiNjLg1vxYGZXSHq3pDMlzZK0xd1n581fJ+lOd79p0utWS1qde9zT07No/fr1Na1zeMx19cNpSdKlr09p7gxr9m3EzsjIiLq7u6MuI1boSTF6Uoh+FKMnhehHMXpSbOXKldvcfXGzy4nVsWFm9klJ75P0Dnc/IOmAmcnM+vNG214pacvk17r7Gklrco/7+/t9xYoVNa13575R3TS4UZK0fPkyLZgzs7k3EkMDAwOqtR+dgp4UoyeF6EcxelKIfhSjJ60Tm92jwUjZByW9092H8matk3RJ8Jxlko6U9Iu2FwgAABChWIy0mdliSV+W9JSkDWYmSaPu/kZJl0u62cyekDQm6QJ3n4isWAAAgAjEIrS5+1ZJJQ8kc/ftks5qb0UAAADxEpvdowAAACiP0AYAAJAAhDYAAIAEILQBAAAkAKENAAAgAQhtAAAACUBoAwAASABCGwAAQAIQ2gAAABKA0AYAAJAAhDYAAIAEqBjazCxlZo+0qxgAAACUVjG0uXta0lYzm9WmegAAAFDCtBqes1nSz83sdkn7chPdfW3LqgIAAECBWkLbXEmPSHpd3jRvTTkAAAAopWJoM7OUpBfd/fI21QMAAIASajmm7XfbVAsAAADKqOWSH3ea2eVmttDMenJ/Wl4ZAAAADqnlmLYvBX9/Udlj2Sz4O9WqogAAAFCoamhzdy7ACwAAELGaApmZLTWzDwU/95rZK1pbFgAAAPJVDW1m9jFJN0r6X8GkPkm3trIoAAAAFKplpO0iSadLGpYkd/+tpIWtLAoAAACFagltY+5+cNK0iVYUAwAAgNJqCW07zOx4BXdBMLMLJD3b0qoAAABQoJZLfvyVpH+WdIKZDUo6IOndLawJAAAAk9RyyY8nzex0SScoe422x4M7JUiSzOxN7n5vC2sEAADoeLWMtMndM5L+o8zsqyWdElpFAAAAKBLGhXMthGUAAACggjBCm4ewjIrM7Dgz+6WZbTaz+83sxFavEwAAIE5q2j0aA1+TdK2732Bm75d0naQzIq4JEUpnXHsOjEmS5vXMUKqr8QHf/GXla2S51ZZVbn6p9U1+bq311FtDtfWW00zf6+lDLctotDfNbjvNCuvfYvLrm9kOqz0344X/T6+lp/Wsr5R6+9Tq7SH/dZP70Ywwv9fKLTdfrd8JcRBmjUl4v+VUDW1m1uvuQ5WeEl45Jde/UNlj5s4KJt0h6StmtsTdB1u57nar5RdZvV/iOcNjrp37RpstMRbSGdel//yQDo5nz4eZNT2lqz90ct0fuuEx1/bhkYJl5at3uZPrmrysf/jAUv3VtzaVnD95faWWVUs9jdSQv9yhUdcHr72vbI3l6q1HpRprXW4j20AjPW3l56ZaPdX6VEq17ayW7bDacw/uS+uU3xupeTut9X2U+7dopE+t2h5KvS6/H80I63ut2nLz1fKd0IiwPzeN/lu1elnltDIEmlf5X4KZ7ZD0HUlfcfdfl5h/obtf15Lqsss/VdLN7n5i3rT7JX3S3e8p97r+/n7fsWNHTevYuW9Uq67fKEm6ftUyLZgzs7miVf//kOd2T9fHv/mQtg1Nvo5xVj1ftqXsHd6rw+ceXvfrpjJ6UoyeFKIfxehJIfpRrNN7UipHmNk2d1/c7LJr2T36GkkXSvq2mb2g7Nmid+Qu+9HKwJZncrIsirBmtlrS6tzjnp4eDQwM1LTw4THX3uFsENqwYYPmzmgsIWfctW9cyrj0T4+lNZZpaDEl7ZX0gWt+2vDrPePaO7w3vIJiYHrwzzTe4N6J/J7M6JL+/HdS6rLm//2qLSt/vlR5fTO6pAtPTOm6f6+vnlpqKLXcXE8m15gvrO27nj6U08g2UE9P2/G5qVZPpX+LnGrbWS3bYbVl5Wos1ZNaelrP+sq9vlqfGv2sNPK63LYX9jbS7PdaOY1+JzSiVZ+bMGsMc1mTNZMjqqk60lbwZLM/lPR/JKUkrZX0D+6+vyWVvbzOhZKekNTn7hNmZpKel3R6pd2j7R5pS2dcl9z6YNmRsnqsOecNmjd7xqHlhjGUu2HDBi1fvrzp2uJkXk+2R43sLpYKe9LssUST66r3eJswjkFqpIbJ03M9qTa830x/ytXYyHIb2Qbq6Wk7Pje1HodWTZjHfVXaXtZ/f6CoJ40cs1lpfZVe36rPSqPbXql+NKPZ77VKy23kO6ERrfrchFljmMsqt+x87Rxpk5kdJulPJV0s6TFJ/yTpTEl3Sfr9ZouoxN1fNLOHJJ0v6QZJZ0sajNvxbHsOjBUFtkaOj1jUO0uv6p9T8JpvfvT0pg+anDvDQtntG0eNvq9KPUl1hdevWpZV7TnN1lPu9ZOn17qdhNmfMJbbyGtqWVc7Pzet+jduZPmVtpdWfG7CqK2dNeRr1TbS6u2u1u+ERrT6c9Pu7+c4qeVEhK9K+mNlTwB4j7s/Hsxab2blLrgbtosk3WBmV0galvThNq23IbmRslrDVX4oK/WapG1UAAAgfLWMtD0p6bXuXmoH9dtDrqekICgm5hIf82bPqCtkEcoAAEA1tdx79EsV5j0fbjkAAAAoJSkX142t3EGMe/aHeyAjAABAPkJbE8I8YxQAAKASQlsdJp8avGd/8Rmji3pnHTplGwAAICyEthpVG1Wr94xRAACAehDaalTqOmw5pa6tBgAAECZCWwPy71ggtfbmsAAAABKhrSH1XocNAACgWYS2KrikBwAAiANCWwVc0gMAAMRFV9QFxFmpkw+4pAcAAIgCI2014pIeAAAgSoS2GnHyAQAAiBK7RwEAABKA0AYAAJAAhDYAAIAE4Ji2Erg2GwAAiBtC2yRcmw0AAMQRu0cn4dpsAAAgjhhpq4BrswEAgLggtFXAtdkAAEBcsHsUAAAgAQhtAAAACUBoAwAASABCGwAAQAIQ2gAAABKA0AYAAJAAhDYAAIAEILQBAAAkQOShzcy+YGb/YWa/NrP7zeztefO6zOxqM/utmT1pZhe3up69B8ZbvQoAAIC6xeGOCD+X9L/c/aCZvUHS3Wb2CncfkXS+pBMlHS/pcEkPmtlP3f03rSrmb7//H61aNAAAQMMiH2lz9x+6e+4O7Y9ISklaEDw+V9JX3T3t7rsl3S7pA2HXMK9nhhb1ziqYxk3iAQBAnMRhpC3fKkm/dfetweNjJD2TN39Q0mlhrzTVZbrmvFO058DYoWncJB4AAMSJuXtrV2D2c0mvKzP7ZHd/NnjemZKul/ROd388mPaIpI+4+8bg8SWSTnX3j5RYz2pJq3OPe3p6Fq1fvz7U95JkIyMj6u7ujrqMWKEnxehJIfpRjJ4Uoh/F6EmxlStXbnP3xc0up+WhraYizN4q6WZJ73b3X+dN/76kG9x9XfD47yQdcPfPVltmf3+/79ixo0UVJ8/AwIBWrFgRdRmxQk+K0ZNC9KMYPSlEP4rRk2JmFkpoi/yYNjN7i7KB7Y/zA1tgnaSLzCxlZvOVPcbttnbXCAAAELXIR9rM7AlJcyU9nzf5And/xMxSkv5R0spg+t+7+1dqXO6EpBdCLTbZ5kjaF3URMUNPitGTQvSjGD0pRD+K0ZNiR7p70+cRRB7aWsXMtoYxFDlV0I9i9KQYPSlEP4rRk0L0oxg9KRZWTyLfPQoAAIDqCG0AAAAJMJVD25qoC4gZ+lGMnhSjJ4XoRzF6Uoh+FKMnxULpyZQ9pg0AAGAqmcojbQAAAFMGoQ0AACABplxoM7PjzOyXZrbZzO43sxOjrqndzGzQzH5jZpuCP+cG0xea2V1m9oSZPWpmb4661lYws38MeuBmdlLe9LLv38x6zOybZvZksO28L5rqW6NCT+42s6fytpX/mjdvyvbEzLrN7DvB+9oUbBdLgnkduZ1U6UlHbieSZGY/MrOHg/f9czNbGkzv1O2kXD86dhuRJDP7n/nfry3bPtx9Sv2R9FNJfxr8/H5J/xp1TRH0YFDSSSWmf0PSZ4Ofl0l6RtK0qOttwft/i6TFk/tQ6f1L+mtlb5kmSccqe2HmeVG/lzb05G5Jf1TmNVO2J5K6Jb1LLx/X+5eSftTJ20mVnnTkdhK8p968n98j6cEO307K9aOTt5FTJP0w2AZOauX2MaVG2sxsobLNuyWYdIekY3P/W4TOkXSNJLn7RknbJU250TZ3v8fdt5aYVen9n5s372lJ90j649ZX2x4VelLJlO2Ju4+4+w88+NaUdJ+kVwU/d+R2UqUnlUzZnkiSuw/lPTxcUib4uVO3k6G8h/n9qGTK9sPMZir73i6WlH9mZ0u2jykV2iQdLek5d5+QpODLZ4ukYyKtKhq3mtkjZvZ1M+s3sz5JXe6+I+85g+qQ3tTw/o9R9n9CpeZNdVcF28ptZpb/S7qTevJxSXeynRT4uKQ78x537HZiZjeZ2bOS/lbShzt9O5ncj7xZnbiNfE7SLUH4ktTa3zdTLbRJhUlXkiySKqL1Fnd/g7Kjjrsk3RhM7/TeVHv/XmHeVHWBu79O0usl/VzS9ybNn/I9MbMrJB0n6dPBpI7fTkr0pKO3E3f/E3c/WtJnJF2VmzzpaR2znZTpR8dtI2Z2hrK7PteWmN2S7WOqhbZnJS02s2mSZGam7OjblkirajN33xL8PS7pHyT9vrvvkiQz68976ivVIb2p4f1vkbSkzLwpy92fDf52d/+KpFcF/0uUOqAnZvZJSe+T9AfufoDtpLgnEttJjrvfKGl57nEnbyfSy/0ws74O3UbeKum1kp42s0FljxsekPS7Umu2jykV2tz9RUkPSTo/mHS2pEF3H4ysqDYzs9lm1ps36YPK9kSS1km6JHjeMklHSvpFWwuMVqX3nz/vWGU/jN+NoMa2MbNpZnZE3uOzJW3PBRdN8Z6Y2WplPx/vnHScTsduJ6V60snbiZnNNbOj8h6/V9m9F7vVgdtJhX4Md+I24u5XuvtR7r7E3ZdI2ipphbv/UK3aPqI+6yLsP5JOkPSvkjZLekDS70RdU5vf/6uUDWkPS3pE0v+VtCSYd4SkH0l6QtJjkt4adb0t6sE1wYdnQtmzcp6s9v4lzZZ0m6Qng23n/VG/j1b3JHjPDwTbya8l/UTSGzqhJ8r+j9gl/VbSpuDPv3XydlKuJx2+nRwt6f689/5jSUs7dTsp149O3kYm9WdQL5892pLtg9tYAQAAJMCU2j0KAAAwVRHaAAAAEoDQBgAAkACENgAAgAQgtAEAACTAtKgLANB5zGxT8OMMScdLejR4/Hjw5zF3v63FNXxf0ufc/d8mTf+YXr6P4ExJv3L381pZSzXB/ZMfcPcFUdYBIFqENgBt5+5LpYIwsrSd6zezOZJep+w1p/Knnybpk5J+1913B3dVObmdtQFAOeweBRArZnaDmf1l8PNnzeybZvY9M3vSzG43s5PN7Kdm9pSZrcl73ZHB/PvN7GEz+1yF1fyBpLu8+EKVR0vaK2lYOnRLngfz1rEsWPcDZvZgcOX33Lw/NLONZvZrM9tkZm8Mpq8Mnvuwmf3MzE4Mpr8teN7a4DWPBaExt7xLgvf8c0l/lje938x+FNyY+2Ezu77+LgNIIkbaAMTdacGffZIelHSlsqFrmrL3/Puqu2+WdKOkz7v7PcH9h79nZu91938pscz3SrqhxPQBSZ+Q9KyZ/UzZ287c6u57gtvDfU3SH7r782a2QNKvzOxeSXMlXSfpLe6+2cymS+oxs4WSbpG03N0fMbPzJN0u6aRgfb8j6c/c/eJgt+znJa0ws9cre7P2k919u5nl35D6fGVvz3eWJJnZ/PraCSCpGGkDEHcD7r7X3dPK3p7t/7n7qLvvV/b4t1eZ2WxJb5f0j8Hxcg9Ieo2yN3MuEASq35O0YfI8z94g/fclvUvSL5W9cfrDQTD6PWVvE/fDYB0/lmTK3jrvnZJ+EIRHufu4u++V9EZJm9z9kWD6rZIWm9krglU+7u4PBD//q6RXBz+/TdL33X178PjavDLvk7TSzL5sZv+fpP019BDAFMBIG4C4G8n7OV3i8TRl/wPqkpa5+3iV5b1d0r3lnhfsMn1I0kNmdrWkf1c2RI1Ketjd3zL5NWZ20uRpuVlBXUWrCf4u9V5yryvJ3f/VzJZKeoeksyX9rZmdHIRaAFMYI20AEs/dX5L0c0mfyk0zs6PMbHGJp79HUqldpjKz1wa7JnOOltQv6SllR96OM7O35z1/qZnNUHa36h+Y2fHB9Olmdriyo2dLzex1wfQPSNrq7i9UeUsbJL0r2L0qSRfmrfNYSfvc/XZJlyp79u2cKssDMAUw0gZgqjhP0hozeyR4vE/SxyRtzT0hOBt0haTLyiyjR9Lfm9mRkg4qO+L1KXffFLz+3ZKuMrO/lzRd0hZJ73H3J83sQknfDHa/piVd5O73m9kFkm41s5SkIUnnVHsj7v6wmX1B0i/N7AVJ38+b/TZJq80sLSkl6bJgVyyAKc6KT54CgKnJzE6X9Bl3/6OoawGAehHaAAAAEoBj2gAAABKA0AYAAJAAhDYAAIAEILQBAAAkAKENAAAgAQhtAAAACUBoAwAASABCGwAAQAL8/3k0KjVUiPjjAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 640x360 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TCLab Model disconnected successfully.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<tclab.historian.Historian at 0x7fb0e234d7f0>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm0AAAFgCAYAAADpfU79AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAxOAAAMTgF/d4wjAAA1x0lEQVR4nO3de7xcdX3v//dnTy47OyHsZGcHJAGDCijlaLikQrVqRElq668qHlCB2kgrFoqnJ8rBg54e66nKKZr2USRHqci9CpHUI17YHjWIYilBiFxaCQibkAAht51NLvs28/n9MWvCzJ77zJpZa+15PR+PPLJnrZm1PvPJmtnvfNfN3F0AAACIt66oCwAAAEB1hDYAAIAEILQBAAAkAKENAAAgAQhtAAAACUBoAwAASABCGwAAQAJMi7qAVpk+fbofccQRUZcRG6Ojo5o5c2bUZcQKPSlGTwrRj2L0pBD9KEZPim3btm3M3ZtuSmxCm5kNShoJ/kjSF939NjNbKOkmSa+WNCrpY+7+i2rL6+3t1datW1tVbuIMDAxoxYoVUZcRK/SkGD0pRD+K0ZNC9KMYPSlmZjvCWE5sQlvg/e7+6KRpV0q6z91XmtkySd82s1e7+0QE9QEAAEQibqGtlHMkHStJ7r7RzLZLerOkuyu9yD2tXTs317SC3t4lSk2b0WSZAAAArRO30HarmXVJ+jdJ/11SRlKXu+cPKw5KOqbagg6kX9LFd55T00qPmjZHX/rgTwluwBSSzqQ1NDpUNL13Zq9SXamaXvdS+iXtOrir7mUAQCtYXG4Yb2bHuPsWM5su6W8l/SdJF0ja4u6z8563TtKd7n7TpNevlrQ697i7d9qiP7zqtTWv/08WrtbMmUc1+S7ia2RkRN3d3VGXESv0pFgzPcl4Rvsz+0OuqDEZZXTL0C0a9/GiedNtus7vPV9dJU6en/y6jGfUZcXPm5eapwt6Lyg5b6rjc1OIfhSjJ8VWrly5zd0XN7uc2IS2fGb2Ckmb3f0wM9svaUlutM3M7pf039z97krLWLBgvj/+m/sqrmdo7zO64p7LJElr3327+hYcH0b5scSBocXi2JNyo0PtsmHDBi1fvrzu16U9rct+dplG0iPVn5wgw3uHNffwuSXnrT1zrfpm9bW5oujF8XMTJfpRjJ4UM7NQQlssdo+a2WxJ0919KJj0QUkPBT+vk3SJpM8GJyIcKanq2aNmqSkdwhCNVoaqOASf4b3DWveTdZGtP2zdqW5d9darlLJUXf3Nve6eu+8pCLFDo0O64hdXHPqZ3aQA2ikWoU3SEZLuMLOUJJP0lKQ/CeZdLulmM3tC0pikCzhzFPUoF7TKHa9UdjkxCFVxlx+S4mByqPrGim/UFLpzrzssdVjZ0bQrfnFFwfslwAFotViENnd/StLJZeZtl3RWeytCXDQ7slUpaMV1VCnK4NPo7tGcuAeXVFeqqV2avTN7ddTso/Tc/uckSSPpEV3600sliQAHoOViEdrQuSqFsriObLU6VEX5C7/SyBKyoe9/v+lKPf38b/S5X1+pkfTooXn7xsf1FwMfkyR1p2bqr9/wKaUspbkzDqu4rfT2HanUNL6KAVTHNwVaKi6hrFTQanRUiVGUZEhPTGho1wtNLWN0/5B2bX/5ziqZ9IR2XH+epmdG9D/keqlLypi0pm9Co/by6/ZLuvyFCyRJM11avWuauvLO+TosI6WUfcH2rm71r7pVXanGvo57e2YoZcHKe+ZLbJvAlEVoQ8Oq7boMK5SFMbJVKmgxqhQ/YQQtqTBcNeP4TFq7HyvcbqYHf6dk6s1kf/6bF6eVDXCjJn1xQeFhuP1p6RM7pykl0/TMiIauO7vhGofMdExfT1DcLL30R9fK6/ysMNoHJAOfUtRkckALc5SsWihjZCsZmg1cYQWtnOnVn9Kw8RKjY/3B31d7WsNjLynt6aJdqDm7Jb105l9p4o7LNW9i9NCoWyMy7hrcmbs+3n7pn95X9zLqGe3LH30k7AHtxacNBUqNnjUb0AhlyZCemCjaHVirsAJX2EGrVLiqx7333qs3velNRdOrhZWFwd83HX1Lwecp/5IhX9p8jXTS0ZrmqUPHv+VUOw5OktLuuuKORzUykVZKaV320v/WTC8OiLWoZ7Qvf/SxnrBHwAOaxycIh6QzaX3yZ588dGZcLWrZdUkoa72wRrmOn9hftDuwVmEFrmaDVr5mg8LM2b3qO6Lx62FOPlt18tmnMtOEZfTXj3yh4HW1non6tb9YrD0HxrIPMu+UHdxdV33Nhu16wl7Tx+4R+gBCW6epdBza0OhQxcBWKqARyEKWSUsHsr940+4ayv1CrviS8Ea50k0tIZzANZV/Oae6UvrSW7+kodGhiiPY5S4lki/32VswZ+bLE+f21F3T/Mt/Vlfgv/fee3XG6W+se5tr9ti9ZkOfNLW3LXQGtt4OkM6k9VL6Jb144MWad3N+4c1fUO/M3oJpBLTG1DoKZp7WYd/7qDR+UJK0ZdcBZWq8zVxYo1xjmqlXXHgboyEtlD/6Nvliv6WCXH6Ay5cLc33dfU19LlPTptU1mjhzdq/6j1pSc9gL7T8VTYY+qfbgV3BGbj7OzkXE+Had4nK7PH+z5zda99PaLiR71OyjtGTuEgJaCfXuhqz3F1bt92coFsYo1/0PPqzXH7WkiSpQj1IX+80FuWrHkubCXLWRuJbVXkfYq3c0L1+YJ6jUGvwKzsgtWEDx2bn1HAfKf2rQLLaeKWby7s9Suzw5MWCSvF2SOaV2TTb6y6PeUbBRm6mrDrtcaaXUPS2lL5x9Uun/9U8Sxi+Erq7Hmno9mldpJE4qHo2rNhIXh895vaN5kzUT+qT6P7uFZ+TmKz47t9RlYcphFy+axb/8FJALatX+Z57b5dlpoazc6Njo/iHtfuGZgl2SOeV2TTa6G7KeUTCfNV9fDv595vXMUKqr8ctBINnK3XbrGyu+oV0ju2oaiSsnTqGummZDn1Rb8Ms/I7dg/U2enZvTzl28tSAAJg//WglX6xmf81Lzkr/LM+SD9I/PpLXrsVTduyQb2Q3JlyPClOpKaWHPwppG4sppNNRlPNNw3VGqNfgVnJGbr8zZueUuC1Pw0gh28dYizACYj2v5tQ6dTJhGd3/e/7P7ExPYSo2Mtfog/fxdkjnldk3yJYS4qDQS1+zt48qFuoPDB3XqgVPrvkNJnEbuKkl1WeEZuflKnJ1b62Vh2r2LtxZhBsB8jV7Lrxqf9fKJIJ26F4LfPAkQxu7PLutqR6ll1XoAf6UvprAP0s//H3L+LsmcTv1SQPKVC3M5zYS6cR+vOEJXTrO3o0tK6CunXbt4a9GKAFhOmMGw0eN965E7c7jc3pyo/9NOaIu5Wnd/RnHGZxhBrJRKI2NhHqTf7IVTgaRqJNTlwtywhhtaZ7XdsdWEcQ/iUpIUBsMIfjlhBcBSGr2WXzUzfVSfGf7cocd7/ym0RR8yZKaj58/Ss7sPltybU8vIYSuDHaEt5mrd/dmKL55KoSzMIFZKuePGOEgfaL1Ku13/5a5/0fLly2teVlj3KW429JXTbBh8Kf2Sdh1sZj9AoXaFyDAD4GT1XsuvmmYOj6lXxl3P7DpQdn5NI4cf/b8t623Noc3M3ijp1fmvcfebWlFUpyp1t4L8x2Gc/ZnOeN5tb9Jlb3tTSygLK4iVEvUQNIBiqa6UDksdVnGUrpRKu2OrCSv0ldNsGBzeO6x1P6ntGpi1aNWIYqtU+n0UajD8yA8OnYi2pMYT0epR6szh/L057dylXElNvxXN7P9IWiFpk16+041LIrSFpJbdoL0ze6t+WVa6vMWLLzxb1w2mawllBDEA1VTbHVtNM6GvnFaHwUa1akSxVUqFzLBHHw9J5R2bPXdWuIuW9MU/Xaahg+OHpvXOmq6uYG9Ol6T+//Iv2rv7xarLOmxef6i15av1N+g7JJ3o7vHauqeIdCatweHB8oHNXQumz1d67z7tGj5Y+jmqfnmLvY+ldHmdtVULZQQxAK3WbOgrJ4wwuGHDhrp2F5cT1xBZTamQGfboY9Ks7V+rPpU5A7lJtf62fT5pgS3jrp37Kl8IcejAmHK7xfccGJOXeX7u2KmCXYshSXtaf3Pf5dq+b9uhaVe8/hM6fPrhkqRMZkK7vvkXmjfxpPZufF+5xRxSy+hYV/4tWkrcliUfoQzAVBVGGGxkd3E5rRhRbJWkhsykq/jb2MzeFfz4SzO7XdK3JB36F3L3H7SwtqbsG5dWXb+x4nN6fKvGurMXivz0dx7VARsq+bxZ01P6hw8s1V99a5MOjhdeKbvL0zrMX2q4zozt1b5Z96grSI/9aWn2o5fKlR2SNUkLDv1Um2qXtyi4GXLPfM1PyJlTADCVtWpEsVXKhcywRh+Tqndmb8uWXW0I5bJJj/PHQF1SbENbmA6Op3XRzb8qmt7laV2+70r1p6vv4y5nqMv1+e5sYPv4rmk6akJKlQhozR47xuUtAABhKhcywxx9RKFqCeBv3P3udhQStjnTpetXLav4nKHdvbrix9kDGz//npPUO/+4gvnpjOvSf36oYHRt1vSUrv7QyUp1mWz/Ds1dv1/S7JpqSsu1V4WnKWfMNWPafsm69MoPXa/e7tIbOrspAQDobNVSwBpJp7SjkLB1WYVbkQRsZIZyewnn9cxQX4nnf/OjpxdcImOe7VPKgt2htk/KXSfsvV+Tesr/zyLtaf23+z6n5w5sL65DkrqmqXfhYv53AgAASqoW2jr+yqUpZbRAw1JmQvr2R6TxMhfd6+mT5hSf5pu79trQ6LCeG90lldm9edTso1q6HxwAACRbtdC2wMwuLjfT3deGXE+8ZNLSug9LQ89Wfl7v0VLP/KLJ5a69lrtIbsEiEnQrFQAA0H7VQluPpHIHhrXmHhJxkElnr7x8YFdxYJveI73/G1JXXut65kuTAle5a69FcY9QAACQfNVC2xZ3X9WWSuKi3Oha7pi1EgFtslIjbGHcggoAAHQuTkec7MDu4sDWe7S04LiqYS1n8k3eGV0DAADNqhbavtqWKuKqjtE1Kf+kg6FD077w5i8Q2AAAQNMqhjZ3/1q7ComcT0j7dmSPY8spc0ZoKeVOOmB3KAAACAO7R3Pu/C/SeOP3FZ28S1TiMh4AACA8hLac8RFJXS8/LnMZj8nK7RLlpAMAABAmQttkTZ4lKmV3iXJnAwAAEKau6k8pZmY/NrPvmtnvhV1Q5HLHsdUwQsYuUQAA0C6NjrT9d0mvlHSOpF+GV05ysUsUAID6ZDIZuSf/Wv1mpq6uhsbB6tJQaHP3jZI2Svp2uOUkF7tEAQCozdjYmLZs2aLx8fGoSwnN9OnTdcwxx2jGjBktW0dNoc3MPidpjaS9kr4n6Y2SLnL3O1pWWeH6j5N0o6QFkoYk/am7/3s71g0AAMK1ZcsWHXbYYerr65OZRV1O09xdu3bt0pYtW/Sa17ymZeupdaTtj939r83snZImJL1J0jcltSW0SfqapGvd/QYze7+k6ySd0fRSZ82TUjOkdHCpjxrPGAUAAI3JZDIaHx9XX1+fpk2bOudD9vX1affu3cpkMi3bVVprtzLB32+VtM7dH29XMjazhZJOkXRWMOkOSV8xsyXuPtjUwru6pHnHSpkJ6U1flOa/uuZbVQEAgPrljmGbCiNs+XLvp5XH6FWMgmb25eDH/Wb2KUkfkPT/zKxLUut22hY6WtJz7j4hSZ7txhZJx4SydJOUmibN7iOwAQDQgZYuXaqlS5fqxBNP1LRp0w49Pvfcc7Vv3z6tWLFCCxYs0IIFCyKt0yolQjN70N1PCY4pu0TSPe6+3sxeI+k/u/sXW16g2amSbnL338mbtlHSJ9z9nrxpqyWtzj3u6elZtH79+orLfin9kq7bc50k6cJ5F+qw1GEVn5/xjPZn9h96vN/361tD36r59VEaGRlRd3d31GXECj0pRk8K0Y9i9KQQ/ShWS0/6+/t17LHHtuWMy3o888wzWr58uZ566qlD00ZHR3Xfffdp3rx5es973lMwL18mk9HTTz+tHTt2FM1buXLlNndf3Gx9NYW2ZlfSjGD36BOS+tx9wrLjj89LOr3S7tH+/n4v1bh8uw7u0sU/uViStPbMtRXP/ix3Id2caq+P2sDAgFasWBF1GbFCT4rRk0L0oxg9KUQ/ilXrSTqd1ubNm3X88ccrlUopnXHtOdD4bSRrMa9nhlJd1XfHDg4O6rTTTtPOnTvrmicVv698ZhZKaKt2TNsJZnZ/uZnu/rvNFlCNu79oZg9JOl/SDZLOljTY9PFsdSp1Id0cLqgLAEBj9hwY06rrN7Z0HdevWqYFc2a2dB3tUC20PSfpsnYUUsVFkm4wsyskDUv6cJTF5C6km8MFdQEAQKtVC20vufvP2lJJBe7+uMK4xEdIuJAuAADhmNczQ9evWtbydUwF1ULb1DofFwAAxEqqy6bErst2qHbaxrvaUkWMpTNp7Tq4S0OjQ1GXAgAAInDKKafojDPO0J49e7R48WJdcMEFkdRRcaTN3Z9vVyFxVO2MUQAAMHUsWbKk5NmhDz74YATVFIvXBVJiptQZo5wpCgAAojB1bvrVYrkzRjlTFAAARKGhkTYzO83Mjgy7mDjLnTFKYAMAAFFodPfolZIeNrNrwywGAAAApTW0e9Td3xHcNP4NIdcDAACAEho+ps3dM5IeCrEWAAAAlNHw2aNmdmOYhcQJ12YDAABx08zZo8tDqyJGuDYbAACdZenSpZKksbExbd68WSeddJIk6YQTTtBtt90WYWWFKoY2M3ux3CxJvaFXEwNcmw0AgM6yadMmSdLg4KBOO+20Q4/jppZ7j54paW+J6fe2pKIY4dpsAAC0WCYtHdjd2nX0zJfq/D2+Y8cOnXfeeXr++edlZjr11FN1/fXXt6jA2lQLbb+S1OfuD0+eYWYvtKak+Mhdmw0AALTIgd3Sre9v7TrO+7Y0p7+ul9xyyy1asmSJfvSjH0mSdu9ucbCsQbUTEc5WmRE1dz81/HIAAACid/rpp+uuu+7SJz7xCX33u9/V7Nmzoy6p6kjb1939g22pBAAAdJ6e+dmRsFavo05nnHGGNm3apB//+Me644479JnPfEYPPfSQUqnoDpeqFtpe25YqAABAZ+pK1b3rsh2efvppLVq0SOecc45WrlyphQsXat++fTr88MMjq6laaPO2VAEAABAjd999t9asWaNUKqV0Oq2rrroq0sAmVQ9t/6nMZT9Mkrv7whbUBAAA0HZLlizRzp07JUmrVq3SqlWrIq6oULXQtlnSu9pRCAAAAMqrFtpG3f2ZtlQCAACAsqpd8sPaUkUMpJ37jQIAgPiqONLm7ie3q5CoXfazyzSSHom6DAAAgJKqjbR1jMmBjfuNAgCAOKl2TFvH4X6jAAAgjhhpmyR3v1ECGwAAnWHp0qVaunSpTjzxRE2bNu3Q43PPPVd33323TjvttKhLlMRIGwAA6HCbNm2SJA0ODuq000479FjKXmQ3LghtAAAgMulMuuVXbmj2kKeJiQldfPHFuvfeezUxMaEbb7wxktE3QhsAAIjM0OiQLv7JxS1dx9oz16pvVl/Dr3/sscf09a9/XWvXrtVXv/pVffrTn9bAwECIFdaGY9oAAAAqOOGEEw6NrJ1xxhn67W9/G0kdjLQBAIDI9M7s1doz17Z8Hc3o7u4+9HMqldLExESTFTWG0AYAACKT6ko1teuyk7B7FAAAIAEIbQAAAJKWLFminTt3Fkx729vepgceeODQ45NOOkmDg4NtriyL0AYAAJAAhDYAAIAEILQBAAAkQOShzcxuMLOtZrYp+HNV3rwuM7vazH5rZk+aWWuvvgcAAFrKzKIuoaVa+f7icsmPK939KyWmny/pREnHSzpc0oNm9lN3/01bqwMAAKEwM5mZxsfHlUo1fmupuBkfHz/03lolLqGtnHMlfdXd05J2m9ntkj4g6bORVgUAABpiZurt7dX27du1aNGiKTHy5u7avn27ent7OyK0rTazj0raIukz7r4pmH6MpGfynjcoqf13aAUAAKFZuHChnnnmGT3xxBNRlxKa7u5uLVy4sKXrMHdv7QrMfi7pdWVmnywpI+l5d8+Y2XslrZV0nLvvM7NHJH3E3TcGy7pE0qnu/pES61ktaXXucU9Pz6L169dXrC3jGd08dLP2pPdIkual5umC3gvUZZEf6he6kZGRgttwgJ6UQk8K0Y9i9KQQ/ShGT4qtXLlym7svbnY5LQ9t9TKzxyV9yN1/ZWbfl3SDu68L5v2dpAPu/tlqy+nv7/cdO3ZUXV86k9bQ6JCk7L3JUl1TZ/96voGBAa1YsSLqMmKFnhSjJ4XoRzF6Uoh+FKMnxcwslNAW+ZCSmS3O+/l0SX2SngwmrZN0kZmlzGy+sse43Rbm+nP3POub1TdlAxsAAEi+OBzTdoOZHSEpLemgpP/s7nuDeTdLWiZpc/D4Knf/jwhqBAAAiFTsdo+GZfr06X7EEUdEXUZsjI6OaubMmVGXESv0pBg9KUQ/itGTQvSjGD0ptm3bNrl706eVxmGkrSV6e3u1devWqMuIDY4xKEZPitGTQvSjGD0pRD+K0ZNiZpYOYzmRH9MGAACA6qbsSFsU0hnXngNjUZdR0vCYa+e+0ajLiBV6UoyeFCrXj3k9M5TqSv4FQQEkC6EtJOmM65JbH9S2oYNRl1LS3uG0bhrcGHUZsUJPitGTQuX6sah3lq457xSCG4Ba7QtjIeweDcmeA2OxDWwAwrVt6GBsR9UBxFIooY2RthZYc84bNG/2jKjLKLBhwwYtX74s6jJihZ4UoyeFJvdjz/4xrb791xFWBKCTEdpaYN7sGVowJ16nO8+dYbGrKWr0pBg9KUQ/AMQJu0cBAAASgNAGAACQAIQ2AACABCC0AQAAJAChDQAAIAEIbQAAAAlAaAMAAEgAQhsAAEACENoAAAASgNAGAACQAIQ2AACABIhVaDOz/2lmbmYnBY8XmtldZvaEmT1qZm+OukYAAIAoxCa0mdkpkk6XtCVv8pWS7nP34yStknSrmXGTewAA0HFiEdrMbKakayRdLMnzZp0TTJe7b5S0XRKjbQAAoOPEIrRJ+pykW9z96dwEM+uT1OXuO/KeNyjpmDbXBgAAEDlz9+rPamUBZmdI+rykM93dzWxQ0h9Jel7SFnefnffcdZLudPebSixntaTVucc9PT2L1q9f3+ryDxkec139cFqSdOnrU5o7w9q27lqMjIyou7s76jJihZ4UoyeFJvcj7p/zdmAbKUQ/itGTYitXrtzm7oubXU4cjg97q6TXSnrazCRpsaQBSX8mSWbWnzfa9koVHvN2iLuvkbQm97i/v99XrFjRwrIL7dw3qpsGN0qSli9fpgVzZrZt3bUYGBhQO/uRBPSkGD0pNLkfcf+ctwPbSCH6UYyetE7ku0fd/Up3P8rdl7j7EklbJa1w9x9KWifpEkkys2WSjpT0i8iKBQAAiEgcRtoquVzSzWb2hKQxSRe4+0TENQEAALRd7EJbMNqW+3m7pLOiqwYAACAeIt89CgAAgOoIbQAAAAlAaAMAAEgAQhsAAEACENoAAAASgNAGAACQAIQ2AACABCC0AQAAJAChDQAAIAEIbQAAAAlAaAMAAEgAQhsAAEACENoAAAASgNAGAACQAIQ2AACABCC0AQAAJAChDQAAIAEIbQAAAAnQ8tBmZje2eh0AAABTXTtG2pa3YR0AAABT2rQwFmJmL5abJak3jHUAAAB0slBCm7Lh7ExJe0tMvzekdQAAAHSssELbryT1ufvDk2eY2QshrQMAAKBjhRXazpY0Jklm1i/poLvvkyR3PzWkdQAAAHSsUE5EcPf9kv7czJ6TtF3SXjN72MzeIUlm1hvGegAAADpVKKHNzP5c0l9KulDSfEl9kj4l6ctmdpakn4SxHgAAgE4V1u7Rj0ta6e5b8qb9wMz+XdJmSWtCWg8AAEBHCus6bV2TApskyd0HJQ26+6dCWg8AAEBHCiu0zTCz7skTzWxWLesws24z+46ZbTazTWZ2l5ktCeYtDB4/YWaPmtmbQ6oZAAAgMcIKbesl3Zx/woGZzZN0k6Q7alzGtZJOcPelkr4XPJakKyXd5+7HSVol6VYzC2u3LgAAQCKEFdo+I2lc0lYze8jMHpT0rKSJYF5F7j7i7j9wdw8m3SfpVcHP50i6JnjeRmXPTmW0DQAAdJRQRqzcfVzSh8zs1ZJOCSY/5O5PNrjIj0u608z6lD1ebkfevEFJxzRcLAAAQALZy4Nb8WBmV0h6t7K3xZolaYu7z86bv07Sne5+06TXrZa0Ove4p6dn0fr169tTtKThMdfVD6clSZe+PqW5M6xt667FyMiIuruLDjvsaPSkGD0pNLkfcf+ctwPbSCH6UYyeFFu5cuU2d1/c7HJidWyYmX1S0vskvcPdD0g6YGYys/680bZXSip1puoa5V1apL+/31esWNGOsiVJO/eN6qbBjZKk5cuXacGcmW1bdy0GBgbUzn4kAT0pRk8KTe5H3D/n7cA2Uoh+FKMnrRPWMW1NC0bKPijpne4+lDdrnaRLgucsk3SkpF+0vUAAAIAIxWKkzcwWS/qypKckbTAzSRp19zdKulzZM1OfUPb+phe4+0RkxQIAAEQgFqHN3bdKKnlwiLtvl3RWeysCAACIl9jsHgUAAEB5hDYAAIAEILQBAAAkAKENAAAgAQhtAAAACUBoAwAASABCGwAAQAIQ2gAAABKA0AYAAJAAhDYAAIAEILQBAAAkAKENAAAgAQhtAAAACUBoAwAASABCGwAAQAIQ2gAAABKA0AYAAJAAhDYAAIAEILQBAAAkAKENAAAgAQhtAAAACUBoAwAASABCGwAAQAIQ2gAAABKA0AYAAJAAhDYAAIAEILQBAAAkAKENAAAgARIR2szsODP7pZltNrP7zezEqGsCAABop0SENklfk3Stux8v6e8kXRdxPQAAAG01LeoCqjGzhZJOkXRWMOkOSV8xsyXuPljudRl37dw32oYKs/bsH2vbugBEr1M/88Nj7f1ujTv6UazTezKvZ4ZSXdaSZcc+tEk6WtJz7j4hSe7uZrZF0jGSBsu9aN+4tOr6je2pEEDHWX37r6MuIRJ7h9O6aZDv1hz6UazTe3L9qmVaMGdmS5adhNAmST7pcVGENbPVklbnHk+bM097h/e2uq4ifd2mjb/YoC5rTcpu1MjIiAYGBqIuI1boSTF6UmhyPzLumjaW0a6RyV9JncMzHsl3a1zRj2Kd3pMNGzZo7ozWZABzj/eXT7B79AlJfe4+YWYm6XlJp1faPdq3YIE/PritTVW+rJXDos0YGBjQihUroi4jVuhJMXpSqFQ/0hnXngOduWtUyv5CWr58edRlxAb9KNbpPSmVA8xsm7svbnbZsR9pc/cXzewhSedLukHS2ZIGKwU2Seoya9nwJIDOlerq7O+WuTM6+/1PRj+K0ZPWiX1oC1wk6QYzu0LSsKQPR1wPAABAWyUitLn745LOiLoOAACAqCQitDViaGhIixc3vft4yhgdHdXMmQxX56MnxehJIfpRjJ4Uoh/F6ElJi8JYyJQNbb29vdq6dWvUZcQGB5gXoyfF6Ekh+lGMnhSiH8XoSTEzS4exnCkb2mqVfyZYXM/8BAAAmLKhrZY7IqQzrkv/+SEdHM8G4EW9s3TNeacQ3AAAQOxM2dDWyB0Rtg0d1J4DY5yqDAAAwrQvjIUk5YbxAAAASRVKaJuyI21zpmfv/1WLdMb1Zzc+0OKKAAAAGjdlQ1s9d0SoduwbAABA1Ng9CgAAkACENgAAgAQgtAEAACQAoQ0AACABCG0AAAAJQGgDAABIAEIbAABAAhDaAAAAEoDQBgAAkACENgAAgASIRWgzs24z+46ZbTazTWZ2l5ktCeYtDB4/YWaPmtmbIy4XAACg7WIR2gLXSjrB3ZdK+l7wWJKulHSfux8naZWkW81syt4zFQAAoJRYhDZ3H3H3H7i7B5Puk/Sq4OdzJF0TPG+jpO2SGG0DAAAdJRahrYSPS7rTzPokdbn7jrx5g5KOiaQqAACAiNjLg1vxYGZXSHq3pDMlzZK0xd1n581fJ+lOd79p0utWS1qde9zT07No/fr1Na1zeMx19cNpSdKlr09p7gxr9m3EzsjIiLq7u6MuI1boSTF6Uoh+FKMnhehHMXpSbOXKldvcfXGzy4nVsWFm9klJ75P0Dnc/IOmAmcnM+vNG214pacvk17r7Gklrco/7+/t9xYoVNa13575R3TS4UZK0fPkyLZgzs7k3EkMDAwOqtR+dgp4UoyeF6EcxelKIfhSjJ60Tm92jwUjZByW9092H8matk3RJ8Jxlko6U9Iu2FwgAABChWIy0mdliSV+W9JSkDWYmSaPu/kZJl0u62cyekDQm6QJ3n4isWAAAgAjEIrS5+1ZJJQ8kc/ftks5qb0UAAADxEpvdowAAACiP0AYAAJAAhDYAAIAEILQBAAAkAKENAAAgAQhtAAAACUBoAwAASABCGwAAQAIQ2gAAABKA0AYAAJAAhDYAAIAEqBjazCxlZo+0qxgAAACUVjG0uXta0lYzm9WmegAAAFDCtBqes1nSz83sdkn7chPdfW3LqgIAAECBWkLbXEmPSHpd3jRvTTkAAAAopWJoM7OUpBfd/fI21QMAAIASajmm7XfbVAsAAADKqOWSH3ea2eVmttDMenJ/Wl4ZAAAADqnlmLYvBX9/Udlj2Sz4O9WqogAAAFCoamhzdy7ACwAAELGaApmZLTWzDwU/95rZK1pbFgAAAPJVDW1m9jFJN0r6X8GkPkm3trIoAAAAFKplpO0iSadLGpYkd/+tpIWtLAoAAACFagltY+5+cNK0iVYUAwAAgNJqCW07zOx4BXdBMLMLJD3b0qoAAABQoJZLfvyVpH+WdIKZDUo6IOndLawJAAAAk9RyyY8nzex0SScoe422x4M7JUiSzOxN7n5vC2sEAADoeLWMtMndM5L+o8zsqyWdElpFAAAAKBLGhXMthGUAAACggjBCm4ewjIrM7Dgz+6WZbTaz+83sxFavEwAAIE5q2j0aA1+TdK2732Bm75d0naQzIq4JEUpnXHsOjEmS5vXMUKqr8QHf/GXla2S51ZZVbn6p9U1+bq311FtDtfWW00zf6+lDLctotDfNbjvNCuvfYvLrm9kOqz0344X/T6+lp/Wsr5R6+9Tq7SH/dZP70Ywwv9fKLTdfrd8JcRBmjUl4v+VUDW1m1uvuQ5WeEl45Jde/UNlj5s4KJt0h6StmtsTdB1u57nar5RdZvV/iOcNjrp37RpstMRbSGdel//yQDo5nz4eZNT2lqz90ct0fuuEx1/bhkYJl5at3uZPrmrysf/jAUv3VtzaVnD95faWWVUs9jdSQv9yhUdcHr72vbI3l6q1HpRprXW4j20AjPW3l56ZaPdX6VEq17ayW7bDacw/uS+uU3xupeTut9X2U+7dopE+t2h5KvS6/H80I63ut2nLz1fKd0IiwPzeN/lu1elnltDIEmlf5X4KZ7ZD0HUlfcfdfl5h/obtf15Lqsss/VdLN7n5i3rT7JX3S3e8p97r+/n7fsWNHTevYuW9Uq67fKEm6ftUyLZgzs7miVf//kOd2T9fHv/mQtg1Nvo5xVj1ftqXsHd6rw+ceXvfrpjJ6UoyeFKIfxehJIfpRrNN7UipHmNk2d1/c7LJr2T36GkkXSvq2mb2g7Nmid+Qu+9HKwJZncrIsirBmtlrS6tzjnp4eDQwM1LTw4THX3uFsENqwYYPmzmgsIWfctW9cyrj0T4+lNZZpaDEl7ZX0gWt+2vDrPePaO7w3vIJiYHrwzzTe4N6J/J7M6JL+/HdS6rLm//2qLSt/vlR5fTO6pAtPTOm6f6+vnlpqKLXcXE8m15gvrO27nj6U08g2UE9P2/G5qVZPpX+LnGrbWS3bYbVl5Wos1ZNaelrP+sq9vlqfGv2sNPK63LYX9jbS7PdaOY1+JzSiVZ+bMGsMc1mTNZMjqqk60lbwZLM/lPR/JKUkrZX0D+6+vyWVvbzOhZKekNTn7hNmZpKel3R6pd2j7R5pS2dcl9z6YNmRsnqsOecNmjd7xqHlhjGUu2HDBi1fvrzp2uJkXk+2R43sLpYKe9LssUST66r3eJswjkFqpIbJ03M9qTa830x/ytXYyHIb2Qbq6Wk7Pje1HodWTZjHfVXaXtZ/f6CoJ40cs1lpfZVe36rPSqPbXql+NKPZ77VKy23kO6ERrfrchFljmMsqt+x87Rxpk5kdJulPJV0s6TFJ/yTpTEl3Sfr9ZouoxN1fNLOHJJ0v6QZJZ0sajNvxbHsOjBUFtkaOj1jUO0uv6p9T8JpvfvT0pg+anDvDQtntG0eNvq9KPUl1hdevWpZV7TnN1lPu9ZOn17qdhNmfMJbbyGtqWVc7Pzet+jduZPmVtpdWfG7CqK2dNeRr1TbS6u2u1u+ERrT6c9Pu7+c4qeVEhK9K+mNlTwB4j7s/Hsxab2blLrgbtosk3WBmV0galvThNq23IbmRslrDVX4oK/WapG1UAAAgfLWMtD0p6bXuXmoH9dtDrqekICgm5hIf82bPqCtkEcoAAEA1tdx79EsV5j0fbjkAAAAoJSkX142t3EGMe/aHeyAjAABAPkJbE8I8YxQAAKASQlsdJp8avGd/8Rmji3pnHTplGwAAICyEthpVG1Wr94xRAACAehDaalTqOmw5pa6tBgAAECZCWwPy71ggtfbmsAAAABKhrSH1XocNAACgWYS2KrikBwAAiANCWwVc0gMAAMRFV9QFxFmpkw+4pAcAAIgCI2014pIeAAAgSoS2GnHyAQAAiBK7RwEAABKA0AYAAJAAhDYAAIAE4Ji2Erg2GwAAiBtC2yRcmw0AAMQRu0cn4dpsAAAgjhhpq4BrswEAgLggtFXAtdkAAEBcsHsUAAAgAQhtAAAACUBoAwAASABCGwAAQAIQ2gAAABKA0AYAAJAAhDYAAIAEILQBAAAkQOShzcy+YGb/YWa/NrP7zeztefO6zOxqM/utmT1pZhe3up69B8ZbvQoAAIC6xeGOCD+X9L/c/aCZvUHS3Wb2CncfkXS+pBMlHS/pcEkPmtlP3f03rSrmb7//H61aNAAAQMMiH2lz9x+6e+4O7Y9ISklaEDw+V9JX3T3t7rsl3S7pA2HXMK9nhhb1ziqYxk3iAQBAnMRhpC3fKkm/dfetweNjJD2TN39Q0mlhrzTVZbrmvFO058DYoWncJB4AAMSJuXtrV2D2c0mvKzP7ZHd/NnjemZKul/ROd388mPaIpI+4+8bg8SWSTnX3j5RYz2pJq3OPe3p6Fq1fvz7U95JkIyMj6u7ujrqMWKEnxehJIfpRjJ4Uoh/F6EmxlStXbnP3xc0up+WhraYizN4q6WZJ73b3X+dN/76kG9x9XfD47yQdcPfPVltmf3+/79ixo0UVJ8/AwIBWrFgRdRmxQk+K0ZNC9KMYPSlEP4rRk2JmFkpoi/yYNjN7i7KB7Y/zA1tgnaSLzCxlZvOVPcbttnbXCAAAELXIR9rM7AlJcyU9nzf5And/xMxSkv5R0spg+t+7+1dqXO6EpBdCLTbZ5kjaF3URMUNPitGTQvSjGD0pRD+K0ZNiR7p70+cRRB7aWsXMtoYxFDlV0I9i9KQYPSlEP4rRk0L0oxg9KRZWTyLfPQoAAIDqCG0AAAAJMJVD25qoC4gZ+lGMnhSjJ4XoRzF6Uoh+FKMnxULpyZQ9pg0AAGAqmcojbQAAAFMGoQ0AACABplxoM7PjzOyXZrbZzO43sxOjrqndzGzQzH5jZpuCP+cG0xea2V1m9oSZPWpmb4661lYws38MeuBmdlLe9LLv38x6zOybZvZksO28L5rqW6NCT+42s6fytpX/mjdvyvbEzLrN7DvB+9oUbBdLgnkduZ1U6UlHbieSZGY/MrOHg/f9czNbGkzv1O2kXD86dhuRJDP7n/nfry3bPtx9Sv2R9FNJfxr8/H5J/xp1TRH0YFDSSSWmf0PSZ4Ofl0l6RtK0qOttwft/i6TFk/tQ6f1L+mtlb5kmSccqe2HmeVG/lzb05G5Jf1TmNVO2J5K6Jb1LLx/X+5eSftTJ20mVnnTkdhK8p968n98j6cEO307K9aOTt5FTJP0w2AZOauX2MaVG2sxsobLNuyWYdIekY3P/W4TOkXSNJLn7RknbJU250TZ3v8fdt5aYVen9n5s372lJ90j649ZX2x4VelLJlO2Ju4+4+w88+NaUdJ+kVwU/d+R2UqUnlUzZnkiSuw/lPTxcUib4uVO3k6G8h/n9qGTK9sPMZir73i6WlH9mZ0u2jykV2iQdLek5d5+QpODLZ4ukYyKtKhq3mtkjZvZ1M+s3sz5JXe6+I+85g+qQ3tTw/o9R9n9CpeZNdVcF28ptZpb/S7qTevJxSXeynRT4uKQ78x537HZiZjeZ2bOS/lbShzt9O5ncj7xZnbiNfE7SLUH4ktTa3zdTLbRJhUlXkiySKqL1Fnd/g7Kjjrsk3RhM7/TeVHv/XmHeVHWBu79O0usl/VzS9ybNn/I9MbMrJB0n6dPBpI7fTkr0pKO3E3f/E3c/WtJnJF2VmzzpaR2znZTpR8dtI2Z2hrK7PteWmN2S7WOqhbZnJS02s2mSZGam7OjblkirajN33xL8PS7pHyT9vrvvkiQz68976ivVIb2p4f1vkbSkzLwpy92fDf52d/+KpFcF/0uUOqAnZvZJSe+T9AfufoDtpLgnEttJjrvfKGl57nEnbyfSy/0ws74O3UbeKum1kp42s0FljxsekPS7Umu2jykV2tz9RUkPSTo/mHS2pEF3H4ysqDYzs9lm1ps36YPK9kSS1km6JHjeMklHSvpFWwuMVqX3nz/vWGU/jN+NoMa2MbNpZnZE3uOzJW3PBRdN8Z6Y2WplPx/vnHScTsduJ6V60snbiZnNNbOj8h6/V9m9F7vVgdtJhX4Md+I24u5XuvtR7r7E3ZdI2ipphbv/UK3aPqI+6yLsP5JOkPSvkjZLekDS70RdU5vf/6uUDWkPS3pE0v+VtCSYd4SkH0l6QtJjkt4adb0t6sE1wYdnQtmzcp6s9v4lzZZ0m6Qng23n/VG/j1b3JHjPDwTbya8l/UTSGzqhJ8r+j9gl/VbSpuDPv3XydlKuJx2+nRwt6f689/5jSUs7dTsp149O3kYm9WdQL5892pLtg9tYAQAAJMCU2j0KAAAwVRHaAAAAEoDQBgAAkACENgAAgAQgtAEAACTAtKgLANB5zGxT8OMMScdLejR4/Hjw5zF3v63FNXxf0ufc/d8mTf+YXr6P4ExJv3L381pZSzXB/ZMfcPcFUdYBIFqENgBt5+5LpYIwsrSd6zezOZJep+w1p/Knnybpk5J+1913B3dVObmdtQFAOeweBRArZnaDmf1l8PNnzeybZvY9M3vSzG43s5PN7Kdm9pSZrcl73ZHB/PvN7GEz+1yF1fyBpLu8+EKVR0vaK2lYOnRLngfz1rEsWPcDZvZgcOX33Lw/NLONZvZrM9tkZm8Mpq8Mnvuwmf3MzE4Mpr8teN7a4DWPBaExt7xLgvf8c0l/lje938x+FNyY+2Ezu77+LgNIIkbaAMTdacGffZIelHSlsqFrmrL3/Puqu2+WdKOkz7v7PcH9h79nZu91938pscz3SrqhxPQBSZ+Q9KyZ/UzZ287c6u57gtvDfU3SH7r782a2QNKvzOxeSXMlXSfpLe6+2cymS+oxs4WSbpG03N0fMbPzJN0u6aRgfb8j6c/c/eJgt+znJa0ws9cre7P2k919u5nl35D6fGVvz3eWJJnZ/PraCSCpGGkDEHcD7r7X3dPK3p7t/7n7qLvvV/b4t1eZ2WxJb5f0j8Hxcg9Ieo2yN3MuEASq35O0YfI8z94g/fclvUvSL5W9cfrDQTD6PWVvE/fDYB0/lmTK3jrvnZJ+EIRHufu4u++V9EZJm9z9kWD6rZIWm9krglU+7u4PBD//q6RXBz+/TdL33X178PjavDLvk7TSzL5sZv+fpP019BDAFMBIG4C4G8n7OV3i8TRl/wPqkpa5+3iV5b1d0r3lnhfsMn1I0kNmdrWkf1c2RI1Ketjd3zL5NWZ20uRpuVlBXUWrCf4u9V5yryvJ3f/VzJZKeoeksyX9rZmdHIRaAFMYI20AEs/dX5L0c0mfyk0zs6PMbHGJp79HUqldpjKz1wa7JnOOltQv6SllR96OM7O35z1/qZnNUHa36h+Y2fHB9Olmdriyo2dLzex1wfQPSNrq7i9UeUsbJL0r2L0qSRfmrfNYSfvc/XZJlyp79u2cKssDMAUw0gZgqjhP0hozeyR4vE/SxyRtzT0hOBt0haTLyiyjR9Lfm9mRkg4qO+L1KXffFLz+3ZKuMrO/lzRd0hZJ73H3J83sQknfDHa/piVd5O73m9kFkm41s5SkIUnnVHsj7v6wmX1B0i/N7AVJ38+b/TZJq80sLSkl6bJgVyyAKc6KT54CgKnJzE6X9Bl3/6OoawGAehHaAAAAEoBj2gAAABKA0AYAAJAAhDYAAIAEILQBAAAkAKENAAAgAQhtAAAACUBoAwAASABCGwAAQAL8/3k0KjVUiPjjAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 640x360 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "from tclab import setup, Historian, Plotter, clock\n",
    "\n",
    "# input\n",
    "def u1(t):\n",
    "    return 50 if 20 <= t <= 200 else 0\n",
    "\n",
    "def experiment(observer, u1, t_final=400, t_step=2):\n",
    "    # execute the event loop\n",
    "    TCLab = setup(connected=False, speedup=20)\n",
    "    with TCLab() as lab:\n",
    "        sources = [('T1', lambda: lab.T1), ('Q1', lab.Q1),\n",
    "                   ('Th', lambda: Th), ('Ts', lambda: Ts), \n",
    "                   ('y_err', lambda: y_err)]\n",
    "        h = Historian(sources)\n",
    "        p = Plotter(h, t_final, layout=[('T1','Ts', 'Th'), ['Q1'], ['y_err']])\n",
    "\n",
    "        # initial input\n",
    "        U1 = u1(0)\n",
    "        lab.Q1(U1)\n",
    "        next(observer)\n",
    "\n",
    "        for t in clock(t_final, t_step):\n",
    "            # get new measurement\n",
    "            T1 = lab.T1\n",
    "            x, y_err = observer.send([t, U1, lab.T1])\n",
    "\n",
    "            # unpack state for historian\n",
    "            Th, Ts = x\n",
    "            p.update(t)\n",
    "\n",
    "            # set input for next period\n",
    "            U1 = u1(t)\n",
    "            lab.Q1(U1)\n",
    "            \n",
    "    return h\n",
    "\n",
    "# create observer/estimator instance\n",
    "L = np.array([[0.4], [0.2]])\n",
    "experiment(tclab_observer(L, x_initial=[50, 50]), u1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[4.4.2 Testing](https://jckantor.github.io/cbe30338-2021/04.04-State-Estimation-Copy1.html#4.4.2-Testing)",
     "section": "4.4.2 Testing"
    }
   },
   "source": [
    "<hr>\n",
    "\n",
    "**Study Question:**\n",
    "\n",
    "<hr>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[4.4.3 Choosing L](https://jckantor.github.io/cbe30338-2021/04.04-State-Estimation-Copy1.html#4.4.3-Choosing-L)",
     "section": "4.4.3 Choosing L"
    }
   },
   "source": [
    "## 4.4.3 Choosing L\n",
    "\n",
    "As this example has made clear, we need a notation that distinguishes the actual state of the physical system from our best estimate of that state. We will use $\\hat{x}$ to denote our best estimate of $x$, and use $\\hat{y}$ to denote the measurement that would result from $\\hat{x}$. We can write two systems of equations corresponding to the dynamics of the actual process and to a model of the process.\n",
    "\n",
    "* **Process**\n",
    "\\begin{align*}\n",
    "\\frac{dx}{dt} & = A x + B_u u + B_d d \\\\\n",
    "y & = C x\n",
    "\\end{align*}\n",
    "\n",
    "* **Prediction Model**\n",
    "\\begin{align*}\n",
    "\\frac{d\\hat{x}}{dt} & = A \\hat{x} + B_u u + B_d \\hat{d}\\\\\n",
    "\\hat{y} & = C \\hat{x}\n",
    "\\end{align*}\n",
    "\n",
    "We have made a bold assumption that the equatons and coefficients in these systems are the same. Any differences are due to different values for the states $x$ and $\\hat{x}$, and differences between unmeasured process disturbances $d$ and our estimate for those disturbances, $\\hat{d}$.\n",
    "\n",
    "As we learned from the example, the problem with using this formulation for the feedback control of unmeasured states is that we are neglecting available measurements. This is a obvious oversight. How can this be addressed?\n",
    "\n",
    "A brilliant answer to this question was given by David Luenberger in the early 1960's while a graduate student at Stanford University. The answer was to modify the second system of equations to include a corrective action due to differences between actual process measurements and predictions of what those measurments should be.\n",
    "\n",
    "* **Observer**\n",
    "\\begin{align*}\n",
    "\\frac{d\\hat{x}}{dt} & = \\underbrace{A \\hat{x} + B_u u + B_d \\hat{d}}_{\\text{prediction}} - \\underbrace{L(\\hat{y} - y)}_{\\text{correction}}\\\\\n",
    "\\hat{y} & = C \\hat{x}\n",
    "\\end{align*}\n",
    "\n",
    "The idea is that if predicted measurements $\\hat{y}$ are different from the actual measurements $y$, then corrective action will be imposed to modify the state estimates. In a sense, this is imposing feedback control on our process model to cause the model states to track the actual process states. We will need find values for a matrrix $L$ (after Luenberger?) to make this work.\n",
    "\n",
    "Taking the difference between predicted and actual measurements, we find\n",
    "\n",
    "\\begin{align*}\n",
    "\\hat{y} - y & = C \\hat{x} - C x \\\\\n",
    "& = C\\underbrace{(\\hat{x} - x)}_e\n",
    "\\end{align*}\n",
    "\n",
    "where $e = \\hat{x} - x$ is error in our prediction of the state. Subtracting the new model equations from the process model, we get an expression for the dynamics of the model error $e = \\hat{x} - x$\n",
    "\n",
    "\\begin{align*}\n",
    "\\frac{de}{dt} & = \\frac{d\\hat{x}}{dt} - \\frac{dx}{dt} \\\\\n",
    "\\\\\n",
    "\\frac{de}{dt} & = (A\\hat{x} + B_u u + B_d \\hat{d}) - (Ax + B_u u B_d d + L(C\\hat{x} - C x)) \\\\\n",
    "\\\\\n",
    "\\frac{de}{dt} & = (A - LC)(\\hat{x} - x) + B_d(\\hat{d} - d) \\\\\n",
    "\\\\\n",
    "\\implies \\frac{de}{dt} & = (A - LC) e + B_d(\\hat{d} - d) \n",
    "\\end{align*}\n",
    "\n",
    "The choice of $L$ determines observer performance. What we seek is a matrix $L$ such that $A - LC$ results in error $e$ that quickly decays to zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[4.4.3 Choosing L](https://jckantor.github.io/cbe30338-2021/04.04-State-Estimation-Copy1.html#4.4.3-Choosing-L)",
     "section": "4.4.3 Choosing L"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnQAAAFgCAYAAAAyxuTfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAxOAAAMTgF/d4wjAAA4hUlEQVR4nO3de5ycdX33/9dnZrO72Zx22WzEJIZ44CBya0CCIipE0I22/qriD7RAbbQtLRTbxgMWaetNFekvmPu+i/DDUzkotRKJbfG0LRqUQ5EgxCBVAuoSkiAhIZslh93NznzuP+aaZGbnfLzm2nk/H488stdhvt/PfPba2c9+v9fB3B0RERERia5Y2AGIiIiISG1U0ImIiIhEnAo6ERERkYhTQSciIiIScSroRERERCJOBZ2IiIhIxKmgExEREYm4jrADCMOMGTP8RS96UdhhtIzx8XG6urrCDqOlKCe5lJNsykcu5SSXcpJN+ci1ffv2CXevOSmRKOjMbBgYC/4BfNbdv2FmC4BbgZcD48Cfuvu9pdrr7e1l27ZtjQo3coaGhhgcHAw7jJainORSTrIpH7mUk1zKSTblI5eZPVePdiJR0AXe6+4/n7LuGuABd19pZsuBb5rZy919MoT4REREREIRpYIun/OAlwK4+0YzexZ4I3B3sRe5J9i9a0vB7b29S4l3dNYxTBEREZHGiVJBd5uZxYCfAH8NJIGYu2cOVQ4DS0o1dCDxApfceV7B7Qs7ZnPt+3+ook6kThLJBCPjI2GHUVcvJF5g98HdAPR29RKPxUOOSETambl72DGUZGZL3H2rmc0APg38D+AiYKu7z8rYbx1wp7vfOuX1q4HV6eXu3o5Fv7PmhKJ9/sGC1XR1Lazju2hdY2NjdHd3hx1GSwkjJ0lPsj+5v6l9VqLak5mTJPnayNc45IcaEFV4kp4kZqkbBcywGVzYeyGxKTcOmBWbdXifdqDPklzKSTblI9fKlSu3u/viWtuJREGXycxeDGxx9zlmth9Ymh6lM7MHgY+7+93F2pg//yh//JcP5Kwf2fsUV/z4YwDc8M7b6Z9/XL3Db0k6STVbIpngW9//FitWrGhen57gYz/6GGOJsdI7h2R07yhz580NO4yWUU4+uuPdrDlzDXE7Mno3nUfz9FmSSznJpnzkMrO6FHQtP+VqZrOAGe4+Eqx6P/BI8PU64FLgU8FFEUcDJa9yNYu3TbE2nU2dxkv/oqxlei9dWO3cs5N1P1hXn0AFyF/cRNmGDRt481lvLlqIjyXGuOyHl2WtWzhrIdeeee20LepEJBwtX9ABLwLuMLM4YMCvgT8Itl0OfNXMngAmgIt0hWu0lVuM5RvR6o538w9v/gcu//HlLT3SVUwrFz0bNmyoadQySiNTiclJRnb/tug+nWMJ4i9MsObkqxideCH79Z7gqp9dw1hiPOd1W0eeYtNvHmDejHl5253bOafs739vTydxs9wNPUdBRHItIvXR8gWdu/8aOLnAtmeBtzU3IqlULUVaJcYSY/zFhr+o6rVTzbAZXPeW65peWLVy0TMnPof+mf1N7TORdPYcmEgtJBPYwecb3mcyMclzN13AjGTx4/C4ZILnH8v/vYoBf4PzQsbpcy/E4B/7U39v/v0P/qxguwMJ+MiuDuLkKdSmGDFjSX9P7oYZM3nhd7+I1+H47e0/mnhHy/+qEGl7+imVmhUr2Bp9blihUblaR7oe/NGDLOhZUK8wp71GFF4Jd6644+eMTSaIk+BjL/wDXZ474tUIM+rQRhyjN3lkeU7SGUjAcyUOyefisKMjtX+mOUlyirykO8O78l1Isx++9J4qI8/2bKybgVW3EYuX/nUxvn+E3c+mbtruM4uPEvb1dBKPlS5aRaQ8KuikpGYXbJUUY+kRrX8a/Ke859NVqx2uTCxnWjFT5i/rrHYaWHhdXpdWqnOoRCFz3333ccYZZ1TU5mc8kTM9m7b30F6u3vw5AK4fyN3eHe/ib1/zCeIWJ+HO1d/5JZOHurCMK2sbUfjOSI4x8pVzy9o3c9Ry3LpYM+dyEuT/OezuiHP1uSflnzIu4vA0s6aVRbKooJPD8hVu9SrYqinSKhGPxZs+HRiGSouwQsqdVsxUbIqxkYVXLHNasY5TiaWUmmrsmtVL/4sqvzCt0LhvIplgydZj2LF/R97tkyT520evPrLiJdAR6+LK067O/rny25mXSNR8ukA1x0imLh/nytGriu6z90uVt3t4mrlBx4KmmCWqdNS2uUQywQuJF9h5YGdNhVupgq2Vzw1rllqLsVp/wU5Vj2nFTI0qvLJO/O85iqOm6XEUj8W59sxry/+jymDSx/m7n3wkp6163S7lqMt/VNExe9999/HGN7yeOd/+Ezh0sOB+W3cfIFnlLbOOTDPXb1o5UyVTzOXIN7KtolEaQUfUNFfOdOnOPTtZ98Pit+ho14ItzBGxfOpdhJWaVsxUaoqxXQqvRio00jz1lIJSI+f5bpdSzSh5vKOjolHIrlm9HHX0MfDB78KBwudRLnVnJH3OZZnSU/uHJicaej5lJVPM5cg3sl3vorEUFZDtQd/haSqRTLB7bHfFo26FPvSjWLDVUoyN7x/huR3DLTkiVkkRVkolH/TVTjFK7fIVelOLPChe6OUr8gqp+RSJWBxm5zkRMBAH+ueUFUqWL/zZ4tTFN8m31v2K53qPgBdT76KxlGYWkCoew6OsR1w9znvLvEVHqxdu5RZptX44H5dMMPJYPNQRsUL0gSlQv9G8fKot/jKfb1tKdefKGvNnB4+fm5vndi01qnSKuRyZI9vNLBozNbOALFU8Frq4qhqHr6Qu48r6dvjcnN7vbppKF3HlflCX+ms7rFt0VDqCVumHYdRHxESqUe5oXj61Fn+je0fLfsJKK9xEe2pRWekUczmmjmw3omgsJIwCslTxWOziqkqNWxefm/0xPrJvTckp+GZPc5ej3r8PWuedSVkSyQQf/dFHC14Jlyn9gdnf3V/0L+Fab9FRzdRmtR80lRZp1RZjmX9VqwiTqKvkKvBGFn+ZKhkFbJRmFJV5Ry3nziz6mkKPMWzGhS3VCqN47PJxrnjh02Xt2+xp7rL8yb/V9Q8I/ZaKmJHxkZxirpbz3hJJZ3TC2bVvvKobwtbyQ1ztCFolRVq1xZjOF5N2VUvxV87j4Rp9s/FKNKOorGTUMq3uN0wvUUDWQwwY+Itvsff5nUX3e/jBBznttNNq6ss8wez//AjzDo0fudl2nivrE55gZGwPu7/+Z/RNjufcmDsx5Wku+W7eXc2+YVFB14KKXZk6Mj4CDiQnuXr5J+jtnEdv11ziSUsdvBlXjo1Q+LYBcOSqsZ179vDtX3+/6ivHapnarGYETSNmIq1havFX7uPhyh0FbJRWKirzKfQYw1YY1azV6Pgo/7Z5qPaGXtRPdyzOmtP+JlXgzuyD2JGKK+t7fNJL6PD44Rtzp7f//ZTnLWfevDtTJftWYk5f4QuHqtGWvxWTHoxITTFyYIL0rZH2HJjA8+yT3VD9ny2Z8AT/c9Nn+e3B/H/lGE5s5CnwJLP/7a+YS4wkkKS6eztdTqqAjI/WNuVQ7dSmijOR9tMKNwJvVlFZzqhlWr5Cs9BoXdszGPMEl/3kU2Xsa0zalBtzA8SgI3ZkSCLn5t3V7lumGwZuoJ+umtrI1Ja/SfcdglU3bcxZ3+PbmOhOPXzxk//6cw7YyOFtMU8wx488sqdRz5YciTnbByZL7jeQgD27DjBahyFfA5bOn5VaqPKGsCrMRCRKmlVUljtqmVboMYZhj2rWSyUFbiGVjLB2xjoBmEjmv+9hJQVzqxfX+g1chpgnuHzfNQwkip8bUI2p8/KZX394dwdzkrmvAejyLj479xM5z0ms5vmImx68j+NXvCW1oBvCioiEplCh2QqjmvVQaYFbSLkFbm9XL0DBfSspmOtdXKdjq5e2LOhmz4CbVi3PWT/yfC9X3JWqqD7zrpPoPepYAGz/c8xdvx+YldtYDY84KjS9mm5pycqr6CvwDfeZR/G5PIVXX08n8Vhlo3aTnb1FbwIqIiLSSiotcEvtW0l7rVpct2VBF7OMm1NmsLFO0gNbfT2d9B/epwvSRdK7vwA9Gd/IKka00hc9vDA+wq7JPXTMyL2sYOGshbx08Stb+ia/IiIi0hrKLujM7HXAyzNf4+63NiKolpFMpJ5HeCDjHkI9/VWNZpW6GfDVb7w6a/i11Z/YICIiIq2jrILOzP5/YBDYBCSC1Q5M34LOJ2HdB2Dk6ZqbKnUz4IWzFrJ07lIVcCIiIlKVckfozgFOdPfWu6yjUQ6O5BZzvS+BnqNKvnTqfeRK3QxYo3EiIiJSi3ILumfaqpibKn3eXE/wIOAiSo3GpadWVcSJiIhIvRQt6MzsHcGX95vZ7cC/AIcLO3f/bgNjax0VnDeXbzQuTVOrIiIi0gilRug+NmU585kjDrRHQVdAvkd0ZS7rQgcRERFphlIF3f9097ubEUjUlJpahVQB14r3qhEREZHpJVZi+9qmRFEDMzvWzO43sy1m9qCZndjI/hLJBLsP7mZ4dLhoMbdw1sK63wVaREREJJ9SI3S1Pyi08b4AfNHdbzaz9wJfAU5vREeFRuWmTq2CpldFRESkeUoVdPPN7JJCG939hjrHUxEzWwCcArwtWHUH8HkzW+ruw/XuL98FD7rQQURERMJWqqDrAXIfepridY6lGi8Bdrj7JIC7u5ltBZYAw/XsaOoFELr9iIiIiLQKcy9cl5nZI+5+chPjqYiZvRa41d1flbFuI/ARd/9xxrrVwOr0ck9Pz6L169fntDc+voNbd6ZOG/xQ7x/xht/cDMBDJ1zOV/b/O3sSew7v+6G+DzEnPqfebykUY2NjdHd3hx1GS1FOcikn2ZSPXMpJLuUkm/KRa+XKldvdfXGt7ZT9LNcW9TSw2Mw63H3SzIzUqN3WzJ3cfS0ZF3gMDAz44OBgTmO7d23htjv/DwCnnbacebtTRd+y008h8cC3mMtcIDXN+u4z3z1tRuaGhobIl492ppzkUk6yKR+5lJNcykk25aNxShV0NzYliiq5+04zewS4ELgZOBcYbsT5c2lXv/FqnTMnIiIiLaVoQefuX2hWIDW4GLjZzK4ARoEPNLIznTMnIiIirSbqU664++M06DYlIiIiIlFQ6sbCIiIiItLiVNCJiIiIRFxVBZ2Z3WVm/25mb6h3QCIiIiJSmWrPoftr4BjgPOD++oUjIiIiAu5++F+UmRmxWOMnRKsq6Nx9I7AR+GZ9w2ktCZwRnJGJvWGHIiIi0haSySQ7d+5kZGQk8sVc2owZM1iyZAmdnZ0N66Osgs7MriJ1Y969wLeB1wEXu/sdDYssZImDe/jojH3ssCRsvAbikb8gWEREpOU99dRTxGIxli5dyowZM8IOp2buzu7du9m6dSuveMUrGtZPuVXK77n735rZW4FJ4Azg68C0LehGfnQ1O2Yks9YtnLWQ3q7ecAISERGZ5pLJJGNjYxx77LF0dEyfgZT+/n6ef/55kslkw6Zfy81WurI5E1jn7o+nnrI1zczsg3gnJCaOrIt3cvWbPkvvzH7dVFhERKSB0lOs063GSL+fRk4hFy0TzexzwZf7zewTwPuA/zSzGNC4ieCwxGLQ91I46uXwO2tT//e9lN6Z/fTP7FcxJyIi0maWLVvGsmXLOPHEE+no6Di8fP7557Nv3z4GBweZP38+8+fPDzXOUiN0K4L/VwGXAh9392fN7BXAbQ2NLCxG6ny5nj6dNyciItLmNm3aBMDw8DCnnnrq4WWA8fFxPv7xj9Pf388555wTToCBsioWd38C+MuM5SeBzzYoJhERERESSWfPgYnSO9agr6eTeKy6Kd6uri7OPvtshoeH6xtUFUoVdMeb2YOFNrr7aXWOR0RERASAPQcmWHXTxob2cdOq5cyf3dXQPpqhVEG3A/hYMwIRERERkeqUKuhecPcfNSUSERERkQx9PZ3ctGp5w/uYDkoVdNPrumERERGJjHjMpsV0aDOUurvdO5oShYiIiEgEnXLKKZx++uns2bOHxYsXc9FFF4USR9EROnd/plmBiIiIiLSqpUuXsmvXrpz1Dz/8cAjR5GrM8ydEREREpGlU0ImIiIhEXFUFnZmdamZH1zsYEREREalctSN01wCbzeyL9QxGRERERCpX1cNK3f0cM4sBr6lzPCIiIiJSoaqfPu/uSeCROsYiIiIiIlWo+qIIM7ulnoEU6ONmM9tmZpuCf2sytsXM7Doz+5WZPWlmlzQ6HhEREZFWVPUIHbCiblEUd427fz7P+guBE4HjgHnAw2b2Q3f/ZZPiEhERkWlu2bJlAExMTLBlyxZOOukkAI4//ni+8Y1vhBhZtqIFnZntLLQJ6K17NJU5H7jR3RPA82Z2O/A+4FOhRiUiIiLTxqZNmwAYHh7m1FNPPbzcasp5luvZwN486+9rSES5VpvZnwBbgSvdfVOwfgnwVMZ+w8CpTYpJREREGi2ZgAPPN7aPnqMgFq/oJc899xwXXHABzzzzDGbGa1/7Wm666aYGBVieUgXdT4F+d988dYOZ/bbWzs3sHuCVBTafDHwSeMbdk2b2buB7Znasu+8L9vHM5or0sxpYnV7u6elhaGgoZ78XEi8wuncUgPvvu//w1xs2bGBOfE7Z7ytqxsbG8uajnSknuZSTbMpHLuUkl3KSrZx8DAwMMDo6SiwWww7soudbH2hoTAfefQveM7/kfqOjo7g7e/fu5Utf+hILFy5k3bp1AOzZs4e9e6eOfR2RTCY5ePAgd911V93inqpUQXcucCjfBnd/ba2du/ubKtj3W2Z2DXA8qUJzK7AU2BjsckywLt9r1wJr08sDAwM+ODiYs9/ug7tZ94PUN+cNZ7yB7933PQBWrFhB/8z+ckONnKGhIfLlo50pJ7mUk2zKRy7lJJdykq1UPhKJBFu2bGHu3LnE43GIT0C8stGzSs2dMxdmzyu939y5mBnz5s1jxYoV3HjjjVx11VWceeaZDA4O0tXVVfC1iUSCmTNncs4556TeVwOUKui+7O7vb0jPZTCzxe6+Lfj69UA/8GSweR1wsZmtJ3VRxPnAylACFRERkfrrOQou+Gbj+6jQ6aefzqZNm7jrrru44447uPLKK3nkkUcaVqyVo1RBd0JToijsZjN7EZAADgL/r7unxzS/CiwHtgTLa9z9FyHEKCIiIo0Qi8PsgbCjyPGb3/yGRYsWcd5557Fy5UoWLFjAvn37mDev9Ehfo5Qq6LzE9oZy93OKbEsAlzYxHBERERHuvvtu1q5dSzweJ5FIsGbNmlCLOShd0P2PArcuMcDdfUEDYhIRERFpKUuXLmXXrl0ArFq1ilWrVoUcUbZSBd0W4B3NCEREREREqlOqoBt396dK7CMiIiIiISr1LNeC93YTERERkdZQtKBz95ObFYiIiIiIVKfUCJ2IiIiItDgVdCIiIiIRp4JOREREpIBly5axbNkyTjzxRDo6Og4vn3/++dx9992ceuqpYYcIlL7KVURERKRtbdq0CYDh4WFOPfXUw8uQusFwq1BBJyIiIi0pkUwwMj7S0D56u3qJx6p/Buvk5CSXXHIJ9913H5OTk9xyyy2hjNqpoBMREZGWNDI+wiU/uKShfdxw9g30z+yv+vWPPfYYX/7yl7nhhhu48cYb+eQnP8nQ0FAdIyyPzqETERERqdLxxx9/eETu9NNP51e/+lUocWiETkRERFpSb1cvN5x9Q8P7qEV3d/fhr+PxOJOTkzVGVB0VdCIiItKS4rF4TdOh7URTriIiIiIRp4JOREREpISlS5eya9eurHVnnXUWDz300OHlk046ieHh4SZHlqKCTkRERCTiVNCJiIiIRJwKOhEREZGIU0EnIiIiLcHMAHD3kCNpjPT7awQVdCIiItISYrEY8XicsbGxsEOpq0OHDmFmDS3odB+6AkbHR8MOQUREpO0MDAywfft2Fi1aRHd3d0OLoGZwd5599ll6e3tV0IVhzUNrwg5BRESk7fT19QGwY8cOEolEyNHUR3d3NwsWLGhoH6EXdGb2QeCvgFcCf+nun8/YFgP+D/AOwIG17n5DxvYrgVXB4j+7+9/UEktvVy8LZy1kx/4dh9ctnLWw5seCiIiISPn6+vro6+sjmUxG/nw6MyMWa/wZbqEXdMBPgfOAv86z7ULgROA4YB7wsJn90N1/aWZvBt4PvBqYBO4zs3vdfajaQOKxONeeeS0j4yOH1/V29RKPxattUkRERKrUjEJougg9U+7+M3f/BZDMs/l84EZ3T7j788DtwPsytt3s7vvdfRz4J1IFXk3Sz41L/1MxJyIiIq0u9IKuhCXAUxnLw8G6UttERERE2kbDp1zN7B5S58flc7K7P12iiczJ86mXhxTblhnDamB1ermjo4PFixeX6LZ9jI+P09XVFXYYLUU5yaWcZFM+ciknuZSTbMpHXovq0UjDCzp3f1MNL98KLAU2BsvHBOsyt5Fn29QY1gJr08sDAwO+bdu2GsKaXoaGhhgcHAw7jJainORSTrIpH7mUk1zKSTblI5eZ1eVS3lafcl0HXGxmcTM7itR5c9/I2PYBM5tlZl3AB4F/CSlOERERkdCEfpWrmV0IXAP0Ab9nZp8A3unujwBfBZYDW4Ld1wQXUODud5vZ7cCjwbZ/cffvNzf68iWSzp4DE2GHkdfohLNr33jYYbQU5SSXcpItMx99PZ3EY9G++amIRFvoBZ27fw34WoFtCeDSIq+9CriqQaHVTSLpXHrbw2wfORh2KHntHU1w6/DG0ju2EeUkl3KSLTMfi3pncv0Fp6ioE5Fq7KtHI60+5Tot7Dkw0bLFnIjUbvvIwZYdgReRlleXgi70Ebp2s/a819A3qzPsMLJs2LCBFSuWhx1GS1FOcikn2TZs2MCy017D6tt/FnYoIiIq6Jqtb1Yn82e31iXbczut5WIKm3KSSznJNrfTWu6PMxFpX5pyFREREYk4FXQiIiIiEaeCTkRERCTiVNCJiIiIRJwKOhEREZGIU0EnIiIiEnEq6EREREQiTgWdiIiISMSpoBMRERGJOBV0IiIiIhGngk5EREQk4iJT0JnZ35mZm9lJwfICM/u+mT1hZj83szeGHaOIiIhIGCJR0JnZKcDrga0Zq68BHnD3Y4FVwG1m1hFGfCIiIiJhavmCzsy6gOuBSwDP2HResB533wg8C2iUTkRERNpOyxd0wFXA19z9N+kVZtYPxNz9uYz9hoElTY5NREREJHTm7qX3ComZnQ58Bjjb3d3MhoHfBZ4Btrr7rIx91wF3uvutedpZDaxOL/f09Cxav359o8M/bHTCuW5zAoDLXh1nbqc1re9yjI2N0d3dHXYYLUU5yaWcZBsbG2Mi1tXSP9vNpmMkl3KSTfnItXLlyu3uvrjWdlr9nLMzgROA35gZwGJgCPgjADMbyBilO4bsc+wOc/e1wNr08sDAgA8ODjYw7Gy79o1z6/BGAFasWM782V1N67scQ0NDNDMfUaCc5FJOsg0NDXHGGWe19M92s+kYyaWcZFM+Gqelp1zd/Rp3X+juS919KbANGHT37wHrgEsBzGw5cDRwb2jBioiIiISk1Ufoirkc+KqZPQFMABe5+2TIMYmIiIg0XaQKumCULv31s8DbwotGREREpDW09JSriIiIiJSmgk5EREQk4lTQiYiIiEScCjoRERGRiFNBJyIiIhJxKuhEREREIk4FnYiIiEjEqaATERERiTgVdCIiIiIRp4JOREREJOJU0ImIiIhEnAo6ERERkYhTQSciIiIScSroRERERCJOBZ2IiIhIxKmgExEREYk4FXQiIiIiEaeCTkRERCTiQi3ozOyWMPsXERERmQ7CHqFbEXL/IiIiIpHX0egOzGxnoU1Ab6P7FxEREZnuGl7QkSrczgb25ll/XxP6FxEREZnWmlHQ/RTod/fNUzeY2W+b0L+IiIjItNaMgu5cYALAzAaAg+6+D8DdX9uE/kVERESmtYZfFOHu+4E/NrMdwLPAXjPbbGbnAJhZb6NjEBEREZnOGl7QmdkfA38OfAg4CugHPgF8zszeBvyg0TGIiIiITGfNmHL9MLDS3bdmrPuumf03sAVY24QYRERERKatZtyHLjalmAPA3YeBYXf/RBNiEBEREZm2mlHQdZpZ99SVZjaznP7NrNvM/tXMtpjZJjP7vpktDbYtCJafMLOfm9kb6x++iIiISGtrRkG3Hvhq5sUPZtYH3ArcUWYbXwSOd/dlwLeDZYBrgAfc/VhgFXCbmTVjGllERESkZTSjoLsSOARsM7NHzOxh4GlgMthWlLuPuft33d2DVQ8ALwu+Pg+4PthvI6mraDVKJyIiIm2l4aNZ7n4I+H0zezlwSrD6EXd/ssomPwzcaWb9pM7Pey5j2zCwpOpgRURERCLIjgx8tT4zuwJ4J6lHic0Etrr7rIzt64A73f3WKa9bDaxOL/f09Cxav359c4IGRiec6zYnALjs1XHmdlrT+i7H2NgY3d05pzm2NeUkl3KSbWxsjIlYV0v/bDebjpFcykk25SPXypUrt7v74lrbicz5Zmb2UeA9wDnufgA4YGaY2UDGKN0xQL4rateScXuUgYEBHxwcbEbYAOzaN86twxsBWLFiOfNndzWt73IMDQ3RzHxEgXKSSznJNjQ0xBlnnNXSP9vNpmMkl3KSTflonGacQ1ezYITt/cBb3X0kY9M64NJgn+XA0cC9TQ9QREREJEQtP0JnZouBzwG/BjaYGcC4u78OuJzUFbRPkHpe7EXuPhlasCIiIiIhaPmCzt23AXlPTHH3Z4G3NTciERERkdYSiSlXERERESlMBZ2IiIhIxKmgExEREYk4FXQiIiIiEaeCTkRERCTiVNCJiIiIRJwKOhEREZGIU0EnIiIiEnEq6EREREQiTgWdiIiISMSpoBMRERGJOBV0IiIiIhGngk5EREQk4lTQiYiIiEScCjoRERGRiFNBJyIiIhJxKuhEREREIk4FnYiIiEjEqaATERERiTgVdCIiIiIRp4JOREREJOJU0ImIiIhEnAo6ERERkYhTQSciIiIScSroRERERCJOBZ2IiIhIxKmgExEREYk4FXQiIiIiERf5gs7MjjWz+81si5k9aGYnhh2TiIiISDNFvqADvgB80d2PA/4/4CshxyMiIiLSVB1hB1ALM1sAnAK8LVh1B/B5M1vq7sOFXpd0Z9e+8SZEmLJn/0TT+hKRcOjnHEYnmvvZGgXKSTbl44i+nk7iMatbe5Eu6ICXADvcfRLA3d3MtgJLgOFCL9p3CFbdtLE5EYpIW1h9+8/CDiF0e0cT3Dqsz9ZMykk25eOIm1YtZ/7srrq1F/WCDsCnLOeUu2a2GlidXu6Y3cfe0b2NjitHf7ex8d4NxKx+FXk9jI2NMTQ0FHYYLUU5yaWcZBsbG2PjvRvomEiye2zqx1B78qSH8tnaypSTbMrHERs2bGBuZ/3qAXOP7gdRMOX6BNDv7pNmZsAzwOuLTbn2z5/vjw9vb1KUR9R7eLVehoaGGBwcDDuMlqKc5FJOsqXzkUg6ew5ouhVSv6BWrFgRdhgtRTnJpnwcka4JzGy7uy+utb1Ij9C5+04zewS4ELgZOBcYLlbMAcTM6jrMKSLtKx7T50na3E7lYirlJJvy0TiRLugCFwM3m9kVwCjwgZDjEREREWmqyBd07v44cHrYcYiIiIiEJfIFXTVGRkZYvLjm6eppY3x8nK4uDYFnUk5yKSfZlI9cykku5SSb8pHXono00pYFXW9vL9u2bQs7jJahk91zKSe5lJNsykcu5SSXcpJN+chlZol6tNOWBV0lNxZu1StTRURERNLasqCr5MbCi3pncv0Fp6ioExERkZY1HZ7l2lDbRw7qHlMiIiLSKPvq0UhbjtDNnpF65EYxe/ZP6FE+IiIi0mgq6KqlGwuLiIjIdKIpVxEREZGIU0EnIiIiEnEq6EREREQiTgWdiIiISMSpoBMRERGJuJYv6Mys28z+1cy2mNkmM/u+mS0Nti0Ilp8ws5+b2RtDDldERESk6Vq+oAt8ETje3ZcB3w6WAa4BHnD3Y4FVwG1m1pa3YhEREZH21fIFnbuPuft33d2DVQ8ALwu+Pg+4PthvI/AsoFE6ERERaSstX9Dl8WHgTjPrB2Lu/lzGtmFgSShRiYiIiITEjgx8tT4zuwJ4J3A2MBPY6u6zMravA+5091unvG41sDq93NPTs2j9+vVF+xqdcK7bnADgslfHmdtp9XobLWdsbIzu7u6ww2gpykku5SSb8pFLOcmlnGRTPnKtXLlyu7svrrWdyJxvZmYfBd4DnOPuB4ADZoaZDWSM0h0DbJ36WndfC6xNLw8MDPjg4GDR/nbtG+fW4Y0ArFixfFo/KmxoaIhS+Wg3ykku5SSb8pFLOcmlnGRTPhonElOuwQjb+4G3uvtIxqZ1wKXBPsuBo4F7mx6giIiISIhafoTOzBYDnwN+DWwwM4Bxd38dcDnwVTN7ApgALnL3ydCCFREREQlByxd07r4NyHsCm7s/C7ytuRGJiIiItJZITLmKiIiISGEq6EREREQiTgWdiIiISMSpoBMRERGJOBV0IiIiIhGngk5EREQk4lTQiYiIiEScCjoRERGRiFNBJyIiIhJxKuhEREREIk4FnYiIiEjEqaATERERiTgVdCIiIiIRp4JOREREJOJU0ImIiIhEnAo6ERERkYhTQSciIiIScVUXdGYWN7NH6xlMlXEca2b3m9kWM3vQzE4MOyYRERGRZuqo9oXunjCzbWY2090P1jOoCn0B+KK732xm7wW+ApweYjyHJZLOngMTh5f7ejqJx6ymNqa2VWh7uf0lks7ohLNr33jVMVYSb6NUkuuw9q338VDoGCjVb6Z8bRRqN+le1vvN11Y576ccfT2dACX7LdV3vn0LtVvo9Zk/N1NjLPWzWaiPcnJaz/dW7edLoe35cpIv1nJirEc8leSp3L4KxVvofRbLSSXHSD2O2XJ+1qr5/Ci0b77t5f6+qeazoZbP8Wo/t8o5vsuJsR6qLugCW4B7zOx2YF96pbvfUGO7ZTGzBcApwNuCVXcAnzezpe4+3IwYCkkknUtve5jtI0dq3Zkz4lz3+yeX/Q1NJJ3L/vkRDh5K5GybOSPO/37fMv7yXzbl3V5Of+n2f7s7wa3DG6uKsZJ4q223mn4L9VfuvqMTzrOjY3Vrt5IYy32fhY6BUv1mytdGoXYP7ktwyhvGir6ffG0Veo+lYsunsyM1qTAxmSzab7G+8/VbqN1ir8/8uclUzs9mvj7KzWm93luxGKvdvnc0NydTYy03xnrEU0meyu2rULyF3me+46TSYyTfvtUcs+X8rFX6+VHuZ1wlv2+q+Wyo5TM/32druZ9b5Rzfhfqtd4FnPuWv7opebHZTntXu7h+sPqSK+n8t8FV3PzFj3YPAR939x4VeNzAw4M8991zRtnftG2fVTamD7qZVy5k/uytnn2J/QezZP8Hq239WztsI3d7RvcybOy/sMFqKcpJLOcmmfORSTnIpJ9mUjyPStYWZbXf3xbW2V/UInZnFgZ3ufnmtQdRoakWaU+6a2WpgdXq5p6eHoaGhoo2OTjh7R1MV9oYNG5jbmWo26c6+Q5B0+NJjCSby/0GZ5QMnxPn6lvL2zaczBn/8qjgxy99v5vZUjOXHBjAD5/eX7OMr/119jJXE2yidMfjQifGy3kepfT3p7B3dW/d2K2mrkBnB9/lQxpFfTr+ljpFS7WbmpFS/+doqZGps+UyNt1C/tX6vKnn9DJwLj9lXNKeF3luxn4tmv7dKPl9KbR8fn6Sra1/J91gsxlLHaTnxVJKnWj5by+kr8zip9Rip9Zgt9rNW6pgs93dMqXgq+X1TzWdDsbZKfq8q+NwqtW85vwcza4t6qHWEboO7r6hbNJX3vwB4Auh390kzM+AZ4PXFplyrHaHLN41ayqLemVx/wSlA4fNDSqnmXKxKzj/YeO8G3r5yZcXnLFQbb6PU47yUtA0bNrBixYq6t1tJW8X6gOrOicpU6blW678zdDgnpfotdU5UqdjyqeT8wWq/V5W8Pv1zUyjGUu+t2vPB6v3eqjkPqtD2oaEhBgcHS77HUjFmqjaeep7DWCiecvqaepzU47zcao/Zas9preY8sULbK/l9U81nQ6m2Kj03slh7xfYt5/dgep/QR+gCd5rZ5cBNZJ9Dd6DGdsvi7jvN7BHgQuBm4FxguFHnz+05MJFTzJU6Dyrzm5pv2rYa8ZiVbKucfdJiZhW/phKNarce/RXad25n7vp6tFtNW4Xke309YizUbr6cVNNWtfL1UajfWvNQzuvTPzf17Lec7Y1+b/U8duvxWdWMPFXymkpzOvU4qUf+qj1mq+2vmn0LbW/E75t6H5OVxFXPfmtVa0F3bfD/Z0lNfVrwf7zGditxMXCzmV0BjAIfaEana897DX2zOht+1YqIiIhIKTUVdO4e+o2J3f1xQrhNSd+szqaOOomIiIgUUnNBZmbLzOz3g697zezFtYclIiIiIuWqqaAzsz8FbgH+PljVD9xWa1AiIiIiUr5aR+guBl5P6tw13P1XwIJagxIRERGR8tVa0E3keezXZI1tioiIiEgFai3onjOz4whu7mtmFwFP1xxVC0oknT37G38vNREREZFK1Xrbkr8E/hk43syGgQPAO2tss+VUc0NhERERkWap9bYlT5rZ64HjSd2D7nF3P/xEWjM7w93vqzHG0E29ofCi3pmH7xItIiIiErZaR+hw9yTwiwKbrwNOqbWPVrL2vNfwsoHZupmwiIiItIxG3xh42lU9fbP0ZAgRERFpLY0u6LzB7YuIiIi0vdAf3SUiIiIitan1SRG9pXappX0RERERKa3WEbonzOxLZvaaAts/X2P7IiIiIlJCrQXdK0hd4fpNM7vHzM4zs3h6o7t/pcb2RURERKSEmgo6d9/r7mvd/VjgGuBaYKuZfdLMZtUlQhEREREpquaLIsxsjpldRqqYewy4DFgIfL/WtkVERESktFovirgR2ELqSRHvcve3u/t6d78UmF9rcGZ2tZn9wsx+ZmYPmtlbMrbFzOw6M/uVmT1pZpfU2p+IiIhIFNX6pIgngRPcfW+ebW/Js65S9wB/7+4Hgwsv7jazF7v7GHAhcCJwHDAPeNjMfujuv6xDvyIiIiKRUes5dNcWKOZw92dqaTto43vunn6I6qNAnCMjf+cDN7p7wt2fB24H3ldrnyIiIiJRE6UbC68CfuXu24LlJcBTGduHg3UiIiIibcXcw3s6l5ndA7yywOaT3f3pYL+zgZuAt7r748G6R4EPuvvGYPlS4LXu/sE8/awGVqeXe3p6Fq1fv75obKMTznWbEwD84Qlxbv5l6uvLXh1nbuf0ul/y2NgY3d3dYYfRUpSTXMpJNuUjl3KSSznJpnzkWrly5XZ3X1xrO7WeQ1cTd39TqX3M7ExSxdw708VcYCuwFNgYLB8TrMvXz1pgbXp5YGDABwcHi/a7a984tw6nmn7DGa/mWzs2A7BixXLmz+4qFXakDA0NUSof7UY5yaWcZFM+ciknuZSTbMpH47T0lKuZvRn4KvB77v6zKZvXARebWdzMjiJ1Tt03mh2jiIiISNhCHaErw1eALuAms8PTnBe5+6OkCr3lpG6bArDG3X/R/BBFREREwtXSBV3wBIpC2xLApU0MR0RERKQltfSUq4iIiIiUpoJOREREJOJU0ImIiIhEnAo6ERERkYhTQSciIiIScSroRERERCJOBZ2IiIhIxKmgExEREYk4FXQiIiIiEaeCTkRERCTiVNCJiIiIRJwKOhEREZGIU0EnIiIiEnEq6EREREQiTgWdiIiISMSpoBMRERGJOBV0IiIiIhEXiYLOzM4ys4SZ/XnGupiZXWdmvzKzJ83skjBjFBEREQlLR9gBlGJmc4B/AL43ZdOFwInAccA84GEz+6G7/7LJIYqIiIiEKgojdGuBNcCuKevPB25094S7Pw/cDryv2cGJiIiIhK2lCzozezvQ6+7fzLN5CfBUxvJwsE5ERESkrZi7h9e52T3AKwtsPhn4NvBWd99pZjcDD7n754PXPgp80N03BsuXAq919w/m6Wc1sDq93NPTs2j9+vVFYxudcK7bnADgD0+Ic/MvU19f9uo4czutkrfZ8sbGxuju7g47jJainORSTrIpH7mUk1zKSTblI9fKlSu3u/viWtsJ9Rw6d39ToW1m9kbgxcCDZgYwH3inmQ24+98BW4GlwMbgJccE6/L1s5bU1C0AAwMDPjg4WDS2XfvGuXU41fSrTn4l83b8AoAVK5Yzf3ZXGe8uOoaGhiiVj3ajnORSTrIpH7mUk1zKSTblo3Fa9qIId78XWJBenjpCB6wDLjaz9aQuijgfWNmIWD79nV80olkRERGRumjpc+hK+CrwOLCF1CjdGnevW+XV19PJot6ZWesW9c6kr6ezXl2IiIiI1EXLjtBN5e5/OGU5AVzaqP7iMeP6C05hz4GJw+v6ejqJx6bX+XMiIiISfZEp6MIQj9m0O19OREREpp8oT7mKiIiICCroRERERCIv1PvQhcXMJoHfhh1HC5kN7As7iBajnORSTrIpH7mUk1zKSTblI9fR7l7zKXDteg7db+txE7/pwsy2KR/ZlJNcykk25SOXcpJLOcmmfOQys231aEdTriIiIiIRp4JOREREJOLataBbW3qXtqJ85FJOcikn2ZSPXMpJLuUkm/KRqy45acuLIkRERESmk3YdoRMRERGZNlTQiYiIiERcWxV0Znasmd1vZlvM7EEzOzHsmJrNzIbN7Jdmtin4d36wfoGZfd/MnjCzn5vZG8OOtVHM7B+DPLiZnZSxvmAOzKzHzL5uZk8Gx897wom+/ork424z+3XGsfJXGdumbT4AzKzbzP41eG+bguNiabCt7Y6TEvlo5+PkP8xsc/C+7zGzZcH6tjtGoGg+2vYYSTOzv8v8jG3IMeLubfMP+CHwh8HX7wX+K+yYQsjBMHBSnvX/BHwq+Ho58BTQEXa8DcrBm4HFU3NRLAfA3wI3B1+/lNSNqfvCfi8NzsfdwO8WeM20zUfwnrqBd3DkPOM/B/6jXY+TEvlo5+OkN+PrdwEPt+sxUiIfbXuMBO/rFOB7wXFwUqOOkbYZoTOzBaSS+rVg1R3AS9N/ZQrnAdcDuPtG4FlgWo7SufuP3T3fjRyL5eD8jG2/AX4M/F7jo228IvkoZtrmA8Ddx9z9ux58ogIPAC8Lvm6746REPoqZlvlIc/eRjMV5QDL4uu2OESiaj2KmbT4AzKyL1Pu7BMi8CrXux0jbFHTAS4Ad7j4JEHwwbQWWhBpVOG4zs0fN7MtmNmBm/UDM3Z/L2GeYNspNGTlYQuovqHzbprM1wbHyDTPL/AXebvn4MHCnjpPDPgzcmbHctseJmd1qZk8DnwY+0O7HyNR8ZGxq12PkKuBrQWEGNO73TTsVdJBdHQNYKFGE683u/hpSo5W7gVuC9cpN6Rx4kW3T0UXu/krg1cA9wLenbG+LfJjZFcCxwCeDVW19nOTJR1sfJ+7+B+7+EuBKYE169ZTd2uYYKZCPtjxGzOx0UtOpN+TZXPdjpJ0KuqeBxWbWAWBmRmrUbmuoUTWZu28N/j8E/G/gTe6+G8DMBjJ2PYY2yk0ZOdgKLC2wbVpy96eD/93dPw+8LPjLEtokH2b2UeA9wNvd/UC7HydT8wE6TtLc/RZgRXq5XY+RtHQ+zKy/jY+RM4ETgN+Y2TCpc5WHgNOg/sdI2xR07r4TeAS4MFh1LjDs7sOhBdVkZjbLzHozVr2fVE4A1gGXBvstB44G7m1qgOErloPMbS8l9YP67yHE2BRm1mFmL8pYPhd4Nl3Q0Ab5MLPVpH5G3jrl3KC2PE7y5aOdjxMzm2tmCzOW301q1uN52vAYKZKP0XY9Rtz9Gndf6O5L3X0psA0YdPfv0YhjJOyrP5r5Dzge+C9gC/AQ8KqwY2ry+38ZqQJuM/Ao8G/A0mDbi4D/AJ4AHgPODDveBubh+uAHa5LU1UNPlsoBMAv4BvBkcPy8N+z30ch8BO/3oeA4+RnwA+A17ZCP4P0tJjXl8StgU/DvJ+16nBTKRzsfJ6RmeB7MeO93Acva+BjJm492Pkby5GiYI1e51v0Y0aO/RERERCKubaZcRURERKYrFXQiIiIiEaeCTkRERCTiVNCJiIiIRJwKOhEREZGI6wg7ABGRTGa2KfiyEzgO+Hmw/Hjw7zF3/0aDY/gOcJW7/2TK+j/lyDMZu4CfuvsFjYyllOB51A+5+/ww4xCRcKmgE5GW4u7LIKtQWdbM/s1sNvBKUvfUylx/KvBR4DR3fz542szJzYxNRKQQTbmKSGSY2c1m9ufB158ys6+b2bfN7Ekzu93MTjazH5rZr81sbcbrjg62P2hmm83sqiLdvB34vufepPMlwF5gFA4/xujhjD6WB30/ZGYPB3fET2/7HTPbaGY/M7NNZva6YP3KYN/NZvYjMzsxWH9WsN8NwWseCwrKdHuXBu/5HuCPMtYPmNl/BA9B32xmN1WeZRGJIo3QiUiUnRr82wc8DFxDqiDrIPX8xBvdfQtwC/AZd/9x8Dznb5vZu939W3nafDdwc571Q8BHgKfN7EekHtNzm7vvCR6p9wXgd9z9GTObD/zUzO4D5gJfAd7s7lvMbAbQY2YLgK8BK9z9UTO7ALgdOCno71XAH7n7JcFU72eAQTN7NfBJ4GR3f9bMMh/8fSGpRxq+DcDMjqosnSISVRqhE5EoG3L3ve6eIPVIu/9093F330/qfLuXmdks4C3APwbn5z0EvILUQ7OzBMXWG4ANU7d56mH0bwLeAdxP6iH1m4Oi6Q2kHq33vaCPuwAj9bjBtwLfDQpL3P2Qu+8FXgdscvdHg/W3AYvN7MVBl4+7+0PB1/8FvDz4+izgO+7+bLD8xYwwHwBWmtnnzOz/AfaXkUMRmQY0QiciUTaW8XUiz3IHqT9cHVju7odKtPcW4L5C+wXTsI8Aj5jZdcB/kyqwxoHN7v7mqa8xs5OmrktvCuLK6Sb4P997Sb8uL3f/LzNbBpwDnAt82sxODgpeEZnGNEInItOau78A3AN8Ir3OzBaa2eI8u78LyDcNi5mdEEx3pr0EGAB+TWrE7lgze0vG/svMrJPUVO3bzey4YP0MM5tHatRtmZm9Mlj/PmCbu/+2xFvaALwjmLIF+FBGny8F9rn77cBlpK4Snl2iPRGZBjRCJyLt4AJgrZk9GizvA/4U2JbeIbhqdRD4WIE2eoD/ZWZHAwdJjZR9wt03Ba9/J7DGzP4XMAPYCrzL3Z80sw8BXw+mdBPAxe7+oJldBNxmZnFgBDiv1Btx981mdjVwv5n9FvhOxuazgNVmlgDiwMeC6V0RmeYs90IuEZH2Y2avB650998NOxYRkUqpoBMRERGJOJ1DJyIiIhJxKuhEREREIk4FnYiIiEjEqaATERERiTgVdCIiIiIRp4JOREREJOJU0ImIiIhEnAo6ERERkYj7vwa7t0l60fbIAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 640x360 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TCLab Model disconnected successfully.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<tclab.historian.Historian at 0x7fb0e270e9d0>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnQAAAFgCAYAAAAyxuTfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAxOAAAMTgF/d4wjAAA4hUlEQVR4nO3de5ycdX33/9dnZrO72Zx22WzEJIZ44CBya0CCIipE0I22/qriD7RAbbQtLRTbxgMWaetNFekvmPu+i/DDUzkotRKJbfG0LRqUQ5EgxCBVAuoSkiAhIZslh93NznzuP+aaZGbnfLzm2nk/H488stdhvt/PfPba2c9+v9fB3B0RERERia5Y2AGIiIiISG1U0ImIiIhEnAo6ERERkYhTQSciIiIScSroRERERCJOBZ2IiIhIxKmgExEREYm4jrADCMOMGTP8RS96UdhhtIzx8XG6urrCDqOlKCe5lJNsykcu5SSXcpJN+ci1ffv2CXevOSmRKOjMbBgYC/4BfNbdv2FmC4BbgZcD48Cfuvu9pdrr7e1l27ZtjQo3coaGhhgcHAw7jJainORSTrIpH7mUk1zKSTblI5eZPVePdiJR0AXe6+4/n7LuGuABd19pZsuBb5rZy919MoT4REREREIRpYIun/OAlwK4+0YzexZ4I3B3sRe5J9i9a0vB7b29S4l3dNYxTBEREZHGiVJBd5uZxYCfAH8NJIGYu2cOVQ4DS0o1dCDxApfceV7B7Qs7ZnPt+3+ook6kThLJBCPjI2GHUVcvJF5g98HdAPR29RKPxUOOSETambl72DGUZGZL3H2rmc0APg38D+AiYKu7z8rYbx1wp7vfOuX1q4HV6eXu3o5Fv7PmhKJ9/sGC1XR1Lazju2hdY2NjdHd3hx1GSwkjJ0lPsj+5v6l9VqLak5mTJPnayNc45IcaEFV4kp4kZqkbBcywGVzYeyGxKTcOmBWbdXifdqDPklzKSTblI9fKlSu3u/viWtuJREGXycxeDGxx9zlmth9Ymh6lM7MHgY+7+93F2pg//yh//JcP5Kwf2fsUV/z4YwDc8M7b6Z9/XL3Db0k6STVbIpngW9//FitWrGhen57gYz/6GGOJsdI7h2R07yhz580NO4yWUU4+uuPdrDlzDXE7Mno3nUfz9FmSSznJpnzkMrO6FHQtP+VqZrOAGe4+Eqx6P/BI8PU64FLgU8FFEUcDJa9yNYu3TbE2nU2dxkv/oqxlei9dWO3cs5N1P1hXn0AFyF/cRNmGDRt481lvLlqIjyXGuOyHl2WtWzhrIdeeee20LepEJBwtX9ABLwLuMLM4YMCvgT8Itl0OfNXMngAmgIt0hWu0lVuM5RvR6o538w9v/gcu//HlLT3SVUwrFz0bNmyoadQySiNTiclJRnb/tug+nWMJ4i9MsObkqxideCH79Z7gqp9dw1hiPOd1W0eeYtNvHmDejHl5253bOafs739vTydxs9wNPUdBRHItIvXR8gWdu/8aOLnAtmeBtzU3IqlULUVaJcYSY/zFhr+o6rVTzbAZXPeW65peWLVy0TMnPof+mf1N7TORdPYcmEgtJBPYwecb3mcyMclzN13AjGTx4/C4ZILnH8v/vYoBf4PzQsbpcy/E4B/7U39v/v0P/qxguwMJ+MiuDuLkKdSmGDFjSX9P7oYZM3nhd7+I1+H47e0/mnhHy/+qEGl7+imVmhUr2Bp9blihUblaR7oe/NGDLOhZUK8wp71GFF4Jd6644+eMTSaIk+BjL/wDXZ474tUIM+rQRhyjN3lkeU7SGUjAcyUOyefisKMjtX+mOUlyirykO8O78l1Isx++9J4qI8/2bKybgVW3EYuX/nUxvn+E3c+mbtruM4uPEvb1dBKPlS5aRaQ8KuikpGYXbJUUY+kRrX8a/Ke859NVqx2uTCxnWjFT5i/rrHYaWHhdXpdWqnOoRCFz3333ccYZZ1TU5mc8kTM9m7b30F6u3vw5AK4fyN3eHe/ib1/zCeIWJ+HO1d/5JZOHurCMK2sbUfjOSI4x8pVzy9o3c9Ry3LpYM+dyEuT/OezuiHP1uSflnzIu4vA0s6aVRbKooJPD8hVu9SrYqinSKhGPxZs+HRiGSouwQsqdVsxUbIqxkYVXLHNasY5TiaWUmmrsmtVL/4sqvzCt0LhvIplgydZj2LF/R97tkyT520evPrLiJdAR6+LK067O/rny25mXSNR8ukA1x0imLh/nytGriu6z90uVt3t4mrlBx4KmmCWqdNS2uUQywQuJF9h5YGdNhVupgq2Vzw1rllqLsVp/wU5Vj2nFTI0qvLJO/O85iqOm6XEUj8W59sxry/+jymDSx/m7n3wkp6163S7lqMt/VNExe9999/HGN7yeOd/+Ezh0sOB+W3cfIFnlLbOOTDPXb1o5UyVTzOXIN7KtolEaQUfUNFfOdOnOPTtZ98Pit+ho14ItzBGxfOpdhJWaVsxUaoqxXQqvRio00jz1lIJSI+f5bpdSzSh5vKOjolHIrlm9HHX0MfDB78KBwudRLnVnJH3OZZnSU/uHJicaej5lJVPM5cg3sl3vorEUFZDtQd/haSqRTLB7bHfFo26FPvSjWLDVUoyN7x/huR3DLTkiVkkRVkolH/TVTjFK7fIVelOLPChe6OUr8gqp+RSJWBxm5zkRMBAH+ueUFUqWL/zZ4tTFN8m31v2K53qPgBdT76KxlGYWkCoew6OsR1w9znvLvEVHqxdu5RZptX44H5dMMPJYPNQRsUL0gSlQv9G8fKot/jKfb1tKdefKGvNnB4+fm5vndi01qnSKuRyZI9vNLBozNbOALFU8Frq4qhqHr6Qu48r6dvjcnN7vbppKF3HlflCX+ms7rFt0VDqCVumHYdRHxESqUe5oXj61Fn+je0fLfsJKK9xEe2pRWekUczmmjmw3omgsJIwCslTxWOziqkqNWxefm/0xPrJvTckp+GZPc5ej3r8PWuedSVkSyQQf/dFHC14Jlyn9gdnf3V/0L+Fab9FRzdRmtR80lRZp1RZjmX9VqwiTqKvkKvBGFn+ZKhkFbJRmFJV5Ry3nziz6mkKPMWzGhS3VCqN47PJxrnjh02Xt2+xp7rL8yb/V9Q8I/ZaKmJHxkZxirpbz3hJJZ3TC2bVvvKobwtbyQ1ztCFolRVq1xZjOF5N2VUvxV87j4Rp9s/FKNKOorGTUMq3uN0wvUUDWQwwY+Itvsff5nUX3e/jBBznttNNq6ss8wez//AjzDo0fudl2nivrE55gZGwPu7/+Z/RNjufcmDsx5Wku+W7eXc2+YVFB14KKXZk6Mj4CDiQnuXr5J+jtnEdv11ziSUsdvBlXjo1Q+LYBcOSqsZ179vDtX3+/6ivHapnarGYETSNmIq1havFX7uPhyh0FbJRWKirzKfQYw1YY1azV6Pgo/7Z5qPaGXtRPdyzOmtP+JlXgzuyD2JGKK+t7fNJL6PD44Rtzp7f//ZTnLWfevDtTJftWYk5f4QuHqtGWvxWTHoxITTFyYIL0rZH2HJjA8+yT3VD9ny2Z8AT/c9Nn+e3B/H/lGE5s5CnwJLP/7a+YS4wkkKS6eztdTqqAjI/WNuVQ7dSmijOR9tMKNwJvVlFZzqhlWr5Cs9BoXdszGPMEl/3kU2Xsa0zalBtzA8SgI3ZkSCLn5t3V7lumGwZuoJ+umtrI1Ja/SfcdglU3bcxZ3+PbmOhOPXzxk//6cw7YyOFtMU8wx488sqdRz5YciTnbByZL7jeQgD27DjBahyFfA5bOn5VaqPKGsCrMRCRKmlVUljtqmVboMYZhj2rWSyUFbiGVjLB2xjoBmEjmv+9hJQVzqxfX+g1chpgnuHzfNQwkip8bUI2p8/KZX394dwdzkrmvAejyLj479xM5z0ms5vmImx68j+NXvCW1oBvCioiEplCh2QqjmvVQaYFbSLkFbm9XL0DBfSspmOtdXKdjq5e2LOhmz4CbVi3PWT/yfC9X3JWqqD7zrpPoPepYAGz/c8xdvx+YldtYDY84KjS9mm5pycqr6CvwDfeZR/G5PIVXX08n8Vhlo3aTnb1FbwIqIiLSSiotcEvtW0l7rVpct2VBF7OMm1NmsLFO0gNbfT2d9B/epwvSRdK7vwA9Gd/IKka00hc9vDA+wq7JPXTMyL2sYOGshbx08Stb+ia/IiIi0hrKLujM7HXAyzNf4+63NiKolpFMpJ5HeCDjHkI9/VWNZpW6GfDVb7w6a/i11Z/YICIiIq2jrILOzP5/YBDYBCSC1Q5M34LOJ2HdB2Dk6ZqbKnUz4IWzFrJ07lIVcCIiIlKVckfozgFOdPfWu6yjUQ6O5BZzvS+BnqNKvnTqfeRK3QxYo3EiIiJSi3ILumfaqpibKn3eXE/wIOAiSo3GpadWVcSJiIhIvRQt6MzsHcGX95vZ7cC/AIcLO3f/bgNjax0VnDeXbzQuTVOrIiIi0gilRug+NmU585kjDrRHQVdAvkd0ZS7rQgcRERFphlIF3f9097ubEUjUlJpahVQB14r3qhEREZHpJVZi+9qmRFEDMzvWzO43sy1m9qCZndjI/hLJBLsP7mZ4dLhoMbdw1sK63wVaREREJJ9SI3S1Pyi08b4AfNHdbzaz9wJfAU5vREeFRuWmTq2CpldFRESkeUoVdPPN7JJCG939hjrHUxEzWwCcArwtWHUH8HkzW+ruw/XuL98FD7rQQURERMJWqqDrAXIfepridY6lGi8Bdrj7JIC7u5ltBZYAw/XsaOoFELr9iIiIiLQKcy9cl5nZI+5+chPjqYiZvRa41d1flbFuI/ARd/9xxrrVwOr0ck9Pz6L169fntDc+voNbd6ZOG/xQ7x/xht/cDMBDJ1zOV/b/O3sSew7v+6G+DzEnPqfebykUY2NjdHd3hx1GS1FOcikn2ZSPXMpJLuUkm/KRa+XKldvdfXGt7ZT9LNcW9TSw2Mw63H3SzIzUqN3WzJ3cfS0ZF3gMDAz44OBgTmO7d23htjv/DwCnnbacebtTRd+y008h8cC3mMtcIDXN+u4z3z1tRuaGhobIl492ppzkUk6yKR+5lJNcykk25aNxShV0NzYliiq5+04zewS4ELgZOBcYbsT5c2lXv/FqnTMnIiIiLaVoQefuX2hWIDW4GLjZzK4ARoEPNLIznTMnIiIirSbqU664++M06DYlIiIiIlFQ6sbCIiIiItLiVNCJiIiIRFxVBZ2Z3WVm/25mb6h3QCIiIiJSmWrPoftr4BjgPOD++oUjIiIiAu5++F+UmRmxWOMnRKsq6Nx9I7AR+GZ9w2ktCZwRnJGJvWGHIiIi0haSySQ7d+5kZGQk8sVc2owZM1iyZAmdnZ0N66Osgs7MriJ1Y969wLeB1wEXu/sdDYssZImDe/jojH3ssCRsvAbikb8gWEREpOU99dRTxGIxli5dyowZM8IOp2buzu7du9m6dSuveMUrGtZPuVXK77n735rZW4FJ4Azg68C0LehGfnQ1O2Yks9YtnLWQ3q7ecAISERGZ5pLJJGNjYxx77LF0dEyfgZT+/n6ef/55kslkw6Zfy81WurI5E1jn7o+nnrI1zczsg3gnJCaOrIt3cvWbPkvvzH7dVFhERKSB0lOs063GSL+fRk4hFy0TzexzwZf7zewTwPuA/zSzGNC4ieCwxGLQ91I46uXwO2tT//e9lN6Z/fTP7FcxJyIi0maWLVvGsmXLOPHEE+no6Di8fP7557Nv3z4GBweZP38+8+fPDzXOUiN0K4L/VwGXAh9392fN7BXAbQ2NLCxG6ny5nj6dNyciItLmNm3aBMDw8DCnnnrq4WWA8fFxPv7xj9Pf388555wTToCBsioWd38C+MuM5SeBzzYoJhERERESSWfPgYnSO9agr6eTeKy6Kd6uri7OPvtshoeH6xtUFUoVdMeb2YOFNrr7aXWOR0RERASAPQcmWHXTxob2cdOq5cyf3dXQPpqhVEG3A/hYMwIRERERkeqUKuhecPcfNSUSERERkQx9PZ3ctGp5w/uYDkoVdNPrumERERGJjHjMpsV0aDOUurvdO5oShYiIiEgEnXLKKZx++uns2bOHxYsXc9FFF4USR9EROnd/plmBiIiIiLSqpUuXsmvXrpz1Dz/8cAjR5GrM8ydEREREpGlU0ImIiIhEXFUFnZmdamZH1zsYEREREalctSN01wCbzeyL9QxGRERERCpX1cNK3f0cM4sBr6lzPCIiIiJSoaqfPu/uSeCROsYiIiIiIlWo+qIIM7ulnoEU6ONmM9tmZpuCf2sytsXM7Doz+5WZPWlmlzQ6HhEREZFWVPUIHbCiblEUd427fz7P+guBE4HjgHnAw2b2Q3f/ZZPiEhERkWlu2bJlAExMTLBlyxZOOukkAI4//ni+8Y1vhBhZtqIFnZntLLQJ6K17NJU5H7jR3RPA82Z2O/A+4FOhRiUiIiLTxqZNmwAYHh7m1FNPPbzcasp5luvZwN486+9rSES5VpvZnwBbgSvdfVOwfgnwVMZ+w8CpTYpJREREGi2ZgAPPN7aPnqMgFq/oJc899xwXXHABzzzzDGbGa1/7Wm666aYGBVieUgXdT4F+d988dYOZ/bbWzs3sHuCVBTafDHwSeMbdk2b2buB7Znasu+8L9vHM5or0sxpYnV7u6elhaGgoZ78XEi8wuncUgPvvu//w1xs2bGBOfE7Z7ytqxsbG8uajnSknuZSTbMpHLuUkl3KSrZx8DAwMMDo6SiwWww7soudbH2hoTAfefQveM7/kfqOjo7g7e/fu5Utf+hILFy5k3bp1AOzZs4e9e6eOfR2RTCY5ePAgd911V93inqpUQXcucCjfBnd/ba2du/ubKtj3W2Z2DXA8qUJzK7AU2BjsckywLt9r1wJr08sDAwM+ODiYs9/ug7tZ94PUN+cNZ7yB7933PQBWrFhB/8z+ckONnKGhIfLlo50pJ7mUk2zKRy7lJJdykq1UPhKJBFu2bGHu3LnE43GIT0C8stGzSs2dMxdmzyu939y5mBnz5s1jxYoV3HjjjVx11VWceeaZDA4O0tXVVfC1iUSCmTNncs4556TeVwOUKui+7O7vb0jPZTCzxe6+Lfj69UA/8GSweR1wsZmtJ3VRxPnAylACFRERkfrrOQou+Gbj+6jQ6aefzqZNm7jrrru44447uPLKK3nkkUcaVqyVo1RBd0JToijsZjN7EZAADgL/r7unxzS/CiwHtgTLa9z9FyHEKCIiIo0Qi8PsgbCjyPGb3/yGRYsWcd5557Fy5UoWLFjAvn37mDev9Ehfo5Qq6LzE9oZy93OKbEsAlzYxHBERERHuvvtu1q5dSzweJ5FIsGbNmlCLOShd0P2PArcuMcDdfUEDYhIRERFpKUuXLmXXrl0ArFq1ilWrVoUcUbZSBd0W4B3NCEREREREqlOqoBt396dK7CMiIiIiISr1LNeC93YTERERkdZQtKBz95ObFYiIiIiIVKfUCJ2IiIiItDgVdCIiIiIRp4JOREREpIBly5axbNkyTjzxRDo6Og4vn3/++dx9992ceuqpYYcIlL7KVURERKRtbdq0CYDh4WFOPfXUw8uQusFwq1BBJyIiIi0pkUwwMj7S0D56u3qJx6p/Buvk5CSXXHIJ9913H5OTk9xyyy2hjNqpoBMREZGWNDI+wiU/uKShfdxw9g30z+yv+vWPPfYYX/7yl7nhhhu48cYb+eQnP8nQ0FAdIyyPzqETERERqdLxxx9/eETu9NNP51e/+lUocWiETkRERFpSb1cvN5x9Q8P7qEV3d/fhr+PxOJOTkzVGVB0VdCIiItKS4rF4TdOh7URTriIiIiIRp4JOREREpISlS5eya9eurHVnnXUWDz300OHlk046ieHh4SZHlqKCTkRERCTiVNCJiIiIRJwKOhEREZGIU0EnIiIiLcHMAHD3kCNpjPT7awQVdCIiItISYrEY8XicsbGxsEOpq0OHDmFmDS3odB+6AkbHR8MOQUREpO0MDAywfft2Fi1aRHd3d0OLoGZwd5599ll6e3tV0IVhzUNrwg5BRESk7fT19QGwY8cOEolEyNHUR3d3NwsWLGhoH6EXdGb2QeCvgFcCf+nun8/YFgP+D/AOwIG17n5DxvYrgVXB4j+7+9/UEktvVy8LZy1kx/4dh9ctnLWw5seCiIiISPn6+vro6+sjmUxG/nw6MyMWa/wZbqEXdMBPgfOAv86z7ULgROA4YB7wsJn90N1/aWZvBt4PvBqYBO4zs3vdfajaQOKxONeeeS0j4yOH1/V29RKPxattUkRERKrUjEJougg9U+7+M3f/BZDMs/l84EZ3T7j788DtwPsytt3s7vvdfRz4J1IFXk3Sz41L/1MxJyIiIq0u9IKuhCXAUxnLw8G6UttERERE2kbDp1zN7B5S58flc7K7P12iiczJ86mXhxTblhnDamB1ermjo4PFixeX6LZ9jI+P09XVFXYYLUU5yaWcZFM+ciknuZSTbMpHXovq0UjDCzp3f1MNL98KLAU2BsvHBOsyt5Fn29QY1gJr08sDAwO+bdu2GsKaXoaGhhgcHAw7jJainORSTrIpH7mUk1zKSTblI5eZ1eVS3lafcl0HXGxmcTM7itR5c9/I2PYBM5tlZl3AB4F/CSlOERERkdCEfpWrmV0IXAP0Ab9nZp8A3unujwBfBZYDW4Ld1wQXUODud5vZ7cCjwbZ/cffvNzf68iWSzp4DE2GHkdfohLNr33jYYbQU5SSXcpItMx99PZ3EY9G++amIRFvoBZ27fw34WoFtCeDSIq+9CriqQaHVTSLpXHrbw2wfORh2KHntHU1w6/DG0ju2EeUkl3KSLTMfi3pncv0Fp6ioE5Fq7KtHI60+5Tot7Dkw0bLFnIjUbvvIwZYdgReRlleXgi70Ebp2s/a819A3qzPsMLJs2LCBFSuWhx1GS1FOcikn2TZs2MCy017D6tt/FnYoIiIq6Jqtb1Yn82e31iXbczut5WIKm3KSSznJNrfTWu6PMxFpX5pyFREREYk4FXQiIiIiEaeCTkRERCTiVNCJiIiIRJwKOhEREZGIU0EnIiIiEnEq6EREREQiTgWdiIiISMSpoBMRERGJOBV0IiIiIhGngk5EREQk4iJT0JnZ35mZm9lJwfICM/u+mT1hZj83szeGHaOIiIhIGCJR0JnZKcDrga0Zq68BHnD3Y4FVwG1m1hFGfCIiIiJhavmCzsy6gOuBSwDP2HResB533wg8C2iUTkRERNpOyxd0wFXA19z9N+kVZtYPxNz9uYz9hoElTY5NREREJHTm7qX3ComZnQ58Bjjb3d3MhoHfBZ4Btrr7rIx91wF3uvutedpZDaxOL/f09Cxav359o8M/bHTCuW5zAoDLXh1nbqc1re9yjI2N0d3dHXYYLUU5yaWcZBsbG2Mi1tXSP9vNpmMkl3KSTfnItXLlyu3uvrjWdlr9nLMzgROA35gZwGJgCPgjADMbyBilO4bsc+wOc/e1wNr08sDAgA8ODjYw7Gy79o1z6/BGAFasWM782V1N67scQ0NDNDMfUaCc5FJOsg0NDXHGGWe19M92s+kYyaWcZFM+Gqelp1zd/Rp3X+juS919KbANGHT37wHrgEsBzGw5cDRwb2jBioiIiISk1Ufoirkc+KqZPQFMABe5+2TIMYmIiIg0XaQKumCULv31s8DbwotGREREpDW09JSriIiIiJSmgk5EREQk4lTQiYiIiEScCjoRERGRiFNBJyIiIhJxKuhEREREIk4FnYiIiEjEqaATERERiTgVdCIiIiIRp4JOREREJOJU0ImIiIhEnAo6ERERkYhTQSciIiIScSroRERERCJOBZ2IiIhIxKmgExEREYk4FXQiIiIiEaeCTkRERCTiQi3ozOyWMPsXERERmQ7CHqFbEXL/IiIiIpHX0egOzGxnoU1Ab6P7FxEREZnuGl7QkSrczgb25ll/XxP6FxEREZnWmlHQ/RTod/fNUzeY2W+b0L+IiIjItNaMgu5cYALAzAaAg+6+D8DdX9uE/kVERESmtYZfFOHu+4E/NrMdwLPAXjPbbGbnAJhZb6NjEBEREZnOGl7QmdkfA38OfAg4CugHPgF8zszeBvyg0TGIiIiITGfNmHL9MLDS3bdmrPuumf03sAVY24QYRERERKatZtyHLjalmAPA3YeBYXf/RBNiEBEREZm2mlHQdZpZ99SVZjaznP7NrNvM/tXMtpjZJjP7vpktDbYtCJafMLOfm9kb6x++iIiISGtrRkG3Hvhq5sUPZtYH3ArcUWYbXwSOd/dlwLeDZYBrgAfc/VhgFXCbmTVjGllERESkZTSjoLsSOARsM7NHzOxh4GlgMthWlLuPuft33d2DVQ8ALwu+Pg+4PthvI6mraDVKJyIiIm2l4aNZ7n4I+H0zezlwSrD6EXd/ssomPwzcaWb9pM7Pey5j2zCwpOpgRURERCLIjgx8tT4zuwJ4J6lHic0Etrr7rIzt64A73f3WKa9bDaxOL/f09Cxav359c4IGRiec6zYnALjs1XHmdlrT+i7H2NgY3d05pzm2NeUkl3KSbWxsjIlYV0v/bDebjpFcykk25SPXypUrt7v74lrbicz5Zmb2UeA9wDnufgA4YGaY2UDGKN0xQL4rateScXuUgYEBHxwcbEbYAOzaN86twxsBWLFiOfNndzWt73IMDQ3RzHxEgXKSSznJNjQ0xBlnnNXSP9vNpmMkl3KSTflonGacQ1ezYITt/cBb3X0kY9M64NJgn+XA0cC9TQ9QREREJEQtP0JnZouBzwG/BjaYGcC4u78OuJzUFbRPkHpe7EXuPhlasCIiIiIhaPmCzt23AXlPTHH3Z4G3NTciERERkdYSiSlXERERESlMBZ2IiIhIxKmgExEREYk4FXQiIiIiEaeCTkRERCTiVNCJiIiIRJwKOhEREZGIU0EnIiIiEnEq6EREREQiTgWdiIiISMSpoBMRERGJOBV0IiIiIhGngk5EREQk4lTQiYiIiEScCjoRERGRiFNBJyIiIhJxKuhEREREIk4FnYiIiEjEqaATERERiTgVdCIiIiIRp4JOREREJOJU0ImIiIhEnAo6ERERkYhTQSciIiIScSroRERERCJOBZ2IiIhIxKmgExEREYk4FXQiIiIiERf5gs7MjjWz+81si5k9aGYnhh2TiIiISDNFvqADvgB80d2PA/4/4CshxyMiIiLSVB1hB1ALM1sAnAK8LVh1B/B5M1vq7sOFXpd0Z9e+8SZEmLJn/0TT+hKRcOjnHEYnmvvZGgXKSTbl44i+nk7iMatbe5Eu6ICXADvcfRLA3d3MtgJLgOFCL9p3CFbdtLE5EYpIW1h9+8/CDiF0e0cT3Dqsz9ZMykk25eOIm1YtZ/7srrq1F/WCDsCnLOeUu2a2GlidXu6Y3cfe0b2NjitHf7ex8d4NxKx+FXk9jI2NMTQ0FHYYLUU5yaWcZBsbG2PjvRvomEiye2zqx1B78qSH8tnaypSTbMrHERs2bGBuZ/3qAXOP7gdRMOX6BNDv7pNmZsAzwOuLTbn2z5/vjw9vb1KUR9R7eLVehoaGGBwcDDuMlqKc5FJOsqXzkUg6ew5ouhVSv6BWrFgRdhgtRTnJpnwcka4JzGy7uy+utb1Ij9C5+04zewS4ELgZOBcYLlbMAcTM6jrMKSLtKx7T50na3E7lYirlJJvy0TiRLugCFwM3m9kVwCjwgZDjEREREWmqyBd07v44cHrYcYiIiIiEJfIFXTVGRkZYvLjm6eppY3x8nK4uDYFnUk5yKSfZlI9cykku5SSb8pHXono00pYFXW9vL9u2bQs7jJahk91zKSe5lJNsykcu5SSXcpJN+chlZol6tNOWBV0lNxZu1StTRURERNLasqCr5MbCi3pncv0Fp6ioExERkZY1HZ7l2lDbRw7qHlMiIiLSKPvq0UhbjtDNnpF65EYxe/ZP6FE+IiIi0mgq6KqlGwuLiIjIdKIpVxEREZGIU0EnIiIiEnEq6EREREQiTgWdiIiISMSpoBMRERGJuJYv6Mys28z+1cy2mNkmM/u+mS0Nti0Ilp8ws5+b2RtDDldERESk6Vq+oAt8ETje3ZcB3w6WAa4BHnD3Y4FVwG1m1pa3YhEREZH21fIFnbuPuft33d2DVQ8ALwu+Pg+4PthvI/AsoFE6ERERaSstX9Dl8WHgTjPrB2Lu/lzGtmFgSShRiYiIiITEjgx8tT4zuwJ4J3A2MBPY6u6zMravA+5091unvG41sDq93NPTs2j9+vVF+xqdcK7bnADgslfHmdtp9XobLWdsbIzu7u6ww2gpykku5SSb8pFLOcmlnGRTPnKtXLlyu7svrrWdyJxvZmYfBd4DnOPuB4ADZoaZDWSM0h0DbJ36WndfC6xNLw8MDPjg4GDR/nbtG+fW4Y0ArFixfFo/KmxoaIhS+Wg3ykku5SSb8pFLOcmlnGRTPhonElOuwQjb+4G3uvtIxqZ1wKXBPsuBo4F7mx6giIiISIhafoTOzBYDnwN+DWwwM4Bxd38dcDnwVTN7ApgALnL3ydCCFREREQlByxd07r4NyHsCm7s/C7ytuRGJiIiItJZITLmKiIiISGEq6EREREQiTgWdiIiISMSpoBMRERGJOBV0IiIiIhGngk5EREQk4lTQiYiIiEScCjoRERGRiFNBJyIiIhJxKuhEREREIk4FnYiIiEjEqaATERERiTgVdCIiIiIRp4JOREREJOJU0ImIiIhEnAo6ERERkYhTQSciIiIScVUXdGYWN7NH6xlMlXEca2b3m9kWM3vQzE4MOyYRERGRZuqo9oXunjCzbWY2090P1jOoCn0B+KK732xm7wW+ApweYjyHJZLOngMTh5f7ejqJx6ymNqa2VWh7uf0lks7ohLNr33jVMVYSb6NUkuuw9q338VDoGCjVb6Z8bRRqN+le1vvN11Y576ccfT2dACX7LdV3vn0LtVvo9Zk/N1NjLPWzWaiPcnJaz/dW7edLoe35cpIv1nJirEc8leSp3L4KxVvofRbLSSXHSD2O2XJ+1qr5/Ci0b77t5f6+qeazoZbP8Wo/t8o5vsuJsR6qLugCW4B7zOx2YF96pbvfUGO7ZTGzBcApwNuCVXcAnzezpe4+3IwYCkkknUtve5jtI0dq3Zkz4lz3+yeX/Q1NJJ3L/vkRDh5K5GybOSPO/37fMv7yXzbl3V5Of+n2f7s7wa3DG6uKsZJ4q223mn4L9VfuvqMTzrOjY3Vrt5IYy32fhY6BUv1mytdGoXYP7ktwyhvGir6ffG0Veo+lYsunsyM1qTAxmSzab7G+8/VbqN1ir8/8uclUzs9mvj7KzWm93luxGKvdvnc0NydTYy03xnrEU0meyu2rULyF3me+46TSYyTfvtUcs+X8rFX6+VHuZ1wlv2+q+Wyo5TM/32druZ9b5Rzfhfqtd4FnPuWv7opebHZTntXu7h+sPqSK+n8t8FV3PzFj3YPAR939x4VeNzAw4M8991zRtnftG2fVTamD7qZVy5k/uytnn2J/QezZP8Hq239WztsI3d7RvcybOy/sMFqKcpJLOcmmfORSTnIpJ9mUjyPStYWZbXf3xbW2V/UInZnFgZ3ufnmtQdRoakWaU+6a2WpgdXq5p6eHoaGhoo2OTjh7R1MV9oYNG5jbmWo26c6+Q5B0+NJjCSby/0GZ5QMnxPn6lvL2zaczBn/8qjgxy99v5vZUjOXHBjAD5/eX7OMr/119jJXE2yidMfjQifGy3kepfT3p7B3dW/d2K2mrkBnB9/lQxpFfTr+ljpFS7WbmpFS/+doqZGps+UyNt1C/tX6vKnn9DJwLj9lXNKeF3luxn4tmv7dKPl9KbR8fn6Sra1/J91gsxlLHaTnxVJKnWj5by+kr8zip9Rip9Zgt9rNW6pgs93dMqXgq+X1TzWdDsbZKfq8q+NwqtW85vwcza4t6qHWEboO7r6hbNJX3vwB4Auh390kzM+AZ4PXFplyrHaHLN41ayqLemVx/wSlA4fNDSqnmXKxKzj/YeO8G3r5yZcXnLFQbb6PU47yUtA0bNrBixYq6t1tJW8X6gOrOicpU6blW678zdDgnpfotdU5UqdjyqeT8wWq/V5W8Pv1zUyjGUu+t2vPB6v3eqjkPqtD2oaEhBgcHS77HUjFmqjaeep7DWCiecvqaepzU47zcao/Zas9preY8sULbK/l9U81nQ6m2Kj03slh7xfYt5/dgep/QR+gCd5rZ5cBNZJ9Dd6DGdsvi7jvN7BHgQuBm4FxguFHnz+05MJFTzJU6Dyrzm5pv2rYa8ZiVbKucfdJiZhW/phKNarce/RXad25n7vp6tFtNW4Xke309YizUbr6cVNNWtfL1UajfWvNQzuvTPzf17Lec7Y1+b/U8duvxWdWMPFXymkpzOvU4qUf+qj1mq+2vmn0LbW/E75t6H5OVxFXPfmtVa0F3bfD/Z0lNfVrwf7zGditxMXCzmV0BjAIfaEana897DX2zOht+1YqIiIhIKTUVdO4e+o2J3f1xQrhNSd+szqaOOomIiIgUUnNBZmbLzOz3g697zezFtYclIiIiIuWqqaAzsz8FbgH+PljVD9xWa1AiIiIiUr5aR+guBl5P6tw13P1XwIJagxIRERGR8tVa0E3keezXZI1tioiIiEgFai3onjOz4whu7mtmFwFP1xxVC0oknT37G38vNREREZFK1Xrbkr8E/hk43syGgQPAO2tss+VUc0NhERERkWap9bYlT5rZ64HjSd2D7nF3P/xEWjM7w93vqzHG0E29ofCi3pmH7xItIiIiErZaR+hw9yTwiwKbrwNOqbWPVrL2vNfwsoHZupmwiIiItIxG3xh42lU9fbP0ZAgRERFpLY0u6LzB7YuIiIi0vdAf3SUiIiIitan1SRG9pXappX0RERERKa3WEbonzOxLZvaaAts/X2P7IiIiIlJCrQXdK0hd4fpNM7vHzM4zs3h6o7t/pcb2RURERKSEmgo6d9/r7mvd/VjgGuBaYKuZfdLMZtUlQhEREREpquaLIsxsjpldRqqYewy4DFgIfL/WtkVERESktFovirgR2ELqSRHvcve3u/t6d78UmF9rcGZ2tZn9wsx+ZmYPmtlbMrbFzOw6M/uVmT1pZpfU2p+IiIhIFNX6pIgngRPcfW+ebW/Js65S9wB/7+4Hgwsv7jazF7v7GHAhcCJwHDAPeNjMfujuv6xDvyIiIiKRUes5dNcWKOZw92dqaTto43vunn6I6qNAnCMjf+cDN7p7wt2fB24H3ldrnyIiIiJRE6UbC68CfuXu24LlJcBTGduHg3UiIiIibcXcw3s6l5ndA7yywOaT3f3pYL+zgZuAt7r748G6R4EPuvvGYPlS4LXu/sE8/awGVqeXe3p6Fq1fv75obKMTznWbEwD84Qlxbv5l6uvLXh1nbuf0ul/y2NgY3d3dYYfRUpSTXMpJNuUjl3KSSznJpnzkWrly5XZ3X1xrO7WeQ1cTd39TqX3M7ExSxdw708VcYCuwFNgYLB8TrMvXz1pgbXp5YGDABwcHi/a7a984tw6nmn7DGa/mWzs2A7BixXLmz+4qFXakDA0NUSof7UY5yaWcZFM+ciknuZSTbMpH47T0lKuZvRn4KvB77v6zKZvXARebWdzMjiJ1Tt03mh2jiIiISNhCHaErw1eALuAms8PTnBe5+6OkCr3lpG6bArDG3X/R/BBFREREwtXSBV3wBIpC2xLApU0MR0RERKQltfSUq4iIiIiUpoJOREREJOJU0ImIiIhEnAo6ERERkYhTQSciIiIScSroRERERCJOBZ2IiIhIxKmgExEREYk4FXQiIiIiEaeCTkRERCTiVNCJiIiIRJwKOhEREZGIU0EnIiIiEnEq6EREREQiTgWdiIiISMSpoBMRERGJOBV0IiIiIhEXiYLOzM4ys4SZ/XnGupiZXWdmvzKzJ83skjBjFBEREQlLR9gBlGJmc4B/AL43ZdOFwInAccA84GEz+6G7/7LJIYqIiIiEKgojdGuBNcCuKevPB25094S7Pw/cDryv2cGJiIiIhK2lCzozezvQ6+7fzLN5CfBUxvJwsE5ERESkrZi7h9e52T3AKwtsPhn4NvBWd99pZjcDD7n754PXPgp80N03BsuXAq919w/m6Wc1sDq93NPTs2j9+vVFYxudcK7bnADgD0+Ic/MvU19f9uo4czutkrfZ8sbGxuju7g47jJainORSTrIpH7mUk1zKSTblI9fKlSu3u/viWtsJ9Rw6d39ToW1m9kbgxcCDZgYwH3inmQ24+98BW4GlwMbgJccE6/L1s5bU1C0AAwMDPjg4WDS2XfvGuXU41fSrTn4l83b8AoAVK5Yzf3ZXGe8uOoaGhiiVj3ajnORSTrIpH7mUk1zKSTblo3Fa9qIId78XWJBenjpCB6wDLjaz9aQuijgfWNmIWD79nV80olkRERGRumjpc+hK+CrwOLCF1CjdGnevW+XV19PJot6ZWesW9c6kr6ezXl2IiIiI1EXLjtBN5e5/OGU5AVzaqP7iMeP6C05hz4GJw+v6ejqJx6bX+XMiIiISfZEp6MIQj9m0O19OREREpp8oT7mKiIiICCroRERERCIv1PvQhcXMJoHfhh1HC5kN7As7iBajnORSTrIpH7mUk1zKSTblI9fR7l7zKXDteg7db+txE7/pwsy2KR/ZlJNcykk25SOXcpJLOcmmfOQys231aEdTriIiIiIRp4JOREREJOLataBbW3qXtqJ85FJOcikn2ZSPXMpJLuUkm/KRqy45acuLIkRERESmk3YdoRMRERGZNlTQiYiIiERcWxV0Znasmd1vZlvM7EEzOzHsmJrNzIbN7Jdmtin4d36wfoGZfd/MnjCzn5vZG8OOtVHM7B+DPLiZnZSxvmAOzKzHzL5uZk8Gx897wom+/ork424z+3XGsfJXGdumbT4AzKzbzP41eG+bguNiabCt7Y6TEvlo5+PkP8xsc/C+7zGzZcH6tjtGoGg+2vYYSTOzv8v8jG3IMeLubfMP+CHwh8HX7wX+K+yYQsjBMHBSnvX/BHwq+Ho58BTQEXa8DcrBm4HFU3NRLAfA3wI3B1+/lNSNqfvCfi8NzsfdwO8WeM20zUfwnrqBd3DkPOM/B/6jXY+TEvlo5+OkN+PrdwEPt+sxUiIfbXuMBO/rFOB7wXFwUqOOkbYZoTOzBaSS+rVg1R3AS9N/ZQrnAdcDuPtG4FlgWo7SufuP3T3fjRyL5eD8jG2/AX4M/F7jo228IvkoZtrmA8Ddx9z9ux58ogIPAC8Lvm6746REPoqZlvlIc/eRjMV5QDL4uu2OESiaj2KmbT4AzKyL1Pu7BMi8CrXux0jbFHTAS4Ad7j4JEHwwbQWWhBpVOG4zs0fN7MtmNmBm/UDM3Z/L2GeYNspNGTlYQuovqHzbprM1wbHyDTPL/AXebvn4MHCnjpPDPgzcmbHctseJmd1qZk8DnwY+0O7HyNR8ZGxq12PkKuBrQWEGNO73TTsVdJBdHQNYKFGE683u/hpSo5W7gVuC9cpN6Rx4kW3T0UXu/krg1cA9wLenbG+LfJjZFcCxwCeDVW19nOTJR1sfJ+7+B+7+EuBKYE169ZTd2uYYKZCPtjxGzOx0UtOpN+TZXPdjpJ0KuqeBxWbWAWBmRmrUbmuoUTWZu28N/j8E/G/gTe6+G8DMBjJ2PYY2yk0ZOdgKLC2wbVpy96eD/93dPw+8LPjLEtokH2b2UeA9wNvd/UC7HydT8wE6TtLc/RZgRXq5XY+RtHQ+zKy/jY+RM4ETgN+Y2TCpc5WHgNOg/sdI2xR07r4TeAS4MFh1LjDs7sOhBdVkZjbLzHozVr2fVE4A1gGXBvstB44G7m1qgOErloPMbS8l9YP67yHE2BRm1mFmL8pYPhd4Nl3Q0Ab5MLPVpH5G3jrl3KC2PE7y5aOdjxMzm2tmCzOW301q1uN52vAYKZKP0XY9Rtz9Gndf6O5L3X0psA0YdPfv0YhjJOyrP5r5Dzge+C9gC/AQ8KqwY2ry+38ZqQJuM/Ao8G/A0mDbi4D/AJ4AHgPODDveBubh+uAHa5LU1UNPlsoBMAv4BvBkcPy8N+z30ch8BO/3oeA4+RnwA+A17ZCP4P0tJjXl8StgU/DvJ+16nBTKRzsfJ6RmeB7MeO93Acva+BjJm492Pkby5GiYI1e51v0Y0aO/RERERCKubaZcRURERKYrFXQiIiIiEaeCTkRERCTiVNCJiIiIRJwKOhEREZGI6wg7ABGRTGa2KfiyEzgO+Hmw/Hjw7zF3/0aDY/gOcJW7/2TK+j/lyDMZu4CfuvsFjYyllOB51A+5+/ww4xCRcKmgE5GW4u7LIKtQWdbM/s1sNvBKUvfUylx/KvBR4DR3fz542szJzYxNRKQQTbmKSGSY2c1m9ufB158ys6+b2bfN7Ekzu93MTjazH5rZr81sbcbrjg62P2hmm83sqiLdvB34vufepPMlwF5gFA4/xujhjD6WB30/ZGYPB3fET2/7HTPbaGY/M7NNZva6YP3KYN/NZvYjMzsxWH9WsN8NwWseCwrKdHuXBu/5HuCPMtYPmNl/BA9B32xmN1WeZRGJIo3QiUiUnRr82wc8DFxDqiDrIPX8xBvdfQtwC/AZd/9x8Dznb5vZu939W3nafDdwc571Q8BHgKfN7EekHtNzm7vvCR6p9wXgd9z9GTObD/zUzO4D5gJfAd7s7lvMbAbQY2YLgK8BK9z9UTO7ALgdOCno71XAH7n7JcFU72eAQTN7NfBJ4GR3f9bMMh/8fSGpRxq+DcDMjqosnSISVRqhE5EoG3L3ve6eIPVIu/9093F330/qfLuXmdks4C3APwbn5z0EvILUQ7OzBMXWG4ANU7d56mH0bwLeAdxP6iH1m4Oi6Q2kHq33vaCPuwAj9bjBtwLfDQpL3P2Qu+8FXgdscvdHg/W3AYvN7MVBl4+7+0PB1/8FvDz4+izgO+7+bLD8xYwwHwBWmtnnzOz/AfaXkUMRmQY0QiciUTaW8XUiz3IHqT9cHVju7odKtPcW4L5C+wXTsI8Aj5jZdcB/kyqwxoHN7v7mqa8xs5OmrktvCuLK6Sb4P997Sb8uL3f/LzNbBpwDnAt82sxODgpeEZnGNEInItOau78A3AN8Ir3OzBaa2eI8u78LyDcNi5mdEEx3pr0EGAB+TWrE7lgze0vG/svMrJPUVO3bzey4YP0MM5tHatRtmZm9Mlj/PmCbu/+2xFvaALwjmLIF+FBGny8F9rn77cBlpK4Snl2iPRGZBjRCJyLt4AJgrZk9GizvA/4U2JbeIbhqdRD4WIE2eoD/ZWZHAwdJjZR9wt03Ba9/J7DGzP4XMAPYCrzL3Z80sw8BXw+mdBPAxe7+oJldBNxmZnFgBDiv1Btx981mdjVwv5n9FvhOxuazgNVmlgDiwMeC6V0RmeYs90IuEZH2Y2avB650998NOxYRkUqpoBMRERGJOJ1DJyIiIhJxKuhEREREIk4FnYiIiEjEqaATERERiTgVdCIiIiIRp4JOREREJOJU0ImIiIhEnAo6ERERkYj7vwa7t0l60fbIAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 640x360 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# create observer/estimator instance\n",
    "L = np.array([[0.4], [0.2]])\n",
    "experiment(tclab_observer(2*L, x_initial=[50, 50]), u1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[4.4.4 Eigenvalues](https://jckantor.github.io/cbe30338-2021/04.04-State-Estimation-Copy1.html#4.4.4-Eigenvalues)",
     "section": "4.4.4 Eigenvalues"
    }
   },
   "source": [
    "## 4.4.4 Eigenvalues\n",
    "\n",
    "Given a system\n",
    "\n",
    "$$\\frac{dx}{dt} = A x$$\n",
    "\n",
    "1. An eigenvalue $\\lambda_i$ and eigenvector $u_i$ satisfy a relationship\n",
    "\n",
    "$$A u_j = \\lambda_j u_j$$\n",
    "\n",
    "2. Complex eigenvalues of real matrices come in complex conjugate pairs determine can be written\n",
    " \n",
    "$$\\lambda_{j,j+1} = \\sigma \\pm i \\omega$$\n",
    " \n",
    "3. Each eigenvalue corresponds to a 'mode'\n",
    " \n",
    "    * Real eigenvalues \n",
    "     \n",
    "    $$e^{\\sigma t}$$\n",
    "     \n",
    "    * Each pair of complex conjugate eigenvalues corresponds to two modes\n",
    "     \n",
    "    $$e^{\\sigma t}\\sin(\\omega t)$$\n",
    "     \n",
    "    $$e^{\\sigma t}\\cos{\\omega t}$$\n",
    "     \n",
    "3. The system is asympototically stable if and only if the eigenvalues of $A$ have negative real part.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[4.4.4 Eigenvalues](https://jckantor.github.io/cbe30338-2021/04.04-State-Estimation-Copy1.html#4.4.4-Eigenvalues)",
     "section": "4.4.4 Eigenvalues"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "eigenvalue = -0.03643\n",
      "time constant =  27.45 seconds\n",
      "eigenvector = [-0.91681931  0.39930233]\n",
      "     A u = [ 0.03339978 -0.01454661]\n",
      "lambda u = [ 0.03339978 -0.01454661]\n",
      "\n",
      "eigenvalue = -0.00690\n",
      "time constant =  145.03 seconds\n",
      "eigenvector = [-0.35205837 -0.93597805]\n",
      "     A u = [0.00242755 0.00645385]\n",
      "lambda u = [0.00242755 0.00645385]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "eigenvalues, eigenvectors = np.linalg.eig(A)\n",
    "print \n",
    "\n",
    "for val, vec in zip(eigenvalues, eigenvectors.T):\n",
    "    print(f\"eigenvalue = {val : 8.5f}\")\n",
    "    print(f\"time constant = {-1/val : 5.2f} seconds\")\n",
    "    print(f\"eigenvector = {vec}\")\n",
    "    print(f\"     A u = {np.dot(A, vec)}\")\n",
    "    print(f\"lambda u = {val*vec}\")\n",
    "    print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[4.4.4 Eigenvalues](https://jckantor.github.io/cbe30338-2021/04.04-State-Estimation-Copy1.html#4.4.4-Eigenvalues)",
     "section": "4.4.4 Eigenvalues"
    }
   },
   "source": [
    "Application to state space model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[4.4.4 Eigenvalues](https://jckantor.github.io/cbe30338-2021/04.04-State-Estimation-Copy1.html#4.4.4-Eigenvalues)",
     "section": "4.4.4 Eigenvalues"
    }
   },
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'tclab_phase_plane' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-6-00b63113dbdb>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     45\u001b[0m \u001b[0mx_ic\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mx_ss\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;36m25\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcos\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx_ss\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;36m30\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0ma\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinspace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     46\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 47\u001b[0;31m \u001b[0mtclab_phase_plane\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx_ic\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;36m50\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mTamb\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     48\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     49\u001b[0m \u001b[0;31m#def tclab_sim(x_ics, [50], [21]):\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'tclab_phase_plane' is not defined"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.gridspec as gridspec\n",
    "\n",
    "u_ss = np.array([50])                # constant manipulable input\n",
    "d_ss = [Tamb]                       # constanat disturbance input\n",
    "x_initial = np.array([Tamb, Tamb])  # initial steady state\n",
    "\n",
    "def tctlab_phase_plane(ic, u_ss, d_ss):\n",
    "    \n",
    "    def deriv(t, x):\n",
    "        return np.dot(A, x) + np.dot(Bu, u_ss) + np.dot(Bd, d_ss)\n",
    "\n",
    "    fig = plt.figure(figsize=(12, 6))\n",
    "    gs = gridspec.GridSpec(2, 3)\n",
    "    ax = fig.add_subplot(gs[:, 0])\n",
    "\n",
    "    bx = fig.add_subplot(gs[0, 1:])\n",
    "    cx = fig.add_subplot(gs[1, 1:])\n",
    "\n",
    "    for ic in x_ic:\n",
    "        soln = solve_ivp(deriv, [0, 500, 0], ic, t_eval=np.linspace(0, 500, 501))\n",
    "        t = soln.t\n",
    "        x = soln.y\n",
    "        ax.plot(x[0], x[1], 'b', label=\"trajectory\")\n",
    "        ax.plot(ic[0], ic[1], 'bo', label=\"initial condition\")\n",
    "        \n",
    "        bx.plot(t, x[0])\n",
    "        cx.plot(t, x[1])\n",
    "    ax.set_xlabel('Th: heater temp / deg C')\n",
    "    ax.set_ylabel('Ts: sensor temp / deg C')\n",
    "    ax.axis('equal')\n",
    "    ax.axis('square')\n",
    "    ax.grid(True)\n",
    "    ax.plot(x_ss[0], x_ss[1], 'r.', ms=20, label=\"steady state\")\n",
    "\n",
    "    eigenvaalues, eigenmatrix = np.linalg.eig(A)\n",
    "    for val, vec in zip(eigenvalues, eigenmatrix.T):\n",
    "        r = 20\n",
    "        ax.plot([x_ss[0], x_ss[0] + r*vec[0]], [x_ss[1], x_ss[1] + r*vec[1]], 'r--', lw=2)\n",
    "        ax.plot([x_ss[0], x_ss[0] - r*vec[0]], [x_ss[1], x_ss[1] - r*vec[1]], 'r--', lw=2)\n",
    "    ax.set_xlim(Tamb-2, 90)\n",
    "    ax.set_ylim(Tamb-2, 90)\n",
    "\n",
    "x_ss = np.linalg.solve(A, -(np.dot(Bd, d_ss) + np.dot(Bu, u_ss)))\n",
    "x_ic = [[x_ss[0] + 25*np.cos(a), x_ss[1] + 30*np.sin(a)] for a in np.linspace(0, 2*np.pi, 1)]\n",
    "\n",
    "tclab_phase_plane(x_ic, [50], [Tamb])\n",
    "\n",
    "#def tclab_sim(x_ics, [50], [21]):\n",
    "#    pass\n",
    "    #_, bx = plt.subplots(2, 1)\n",
    "        #bx[0].plot(t, x[0])\n",
    "        #bx[0].set_title('heater temp')\n",
    "        #bx[0].grid(True)\n",
    "        #bx[1].plot(t, x[1])\n",
    "        #bx[1].set_title('sensor temp')\n",
    "        #bx[1].grid(True)\n",
    "        #plt.tight_layout()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[4.4.4 Eigenvalues](https://jckantor.github.io/cbe30338-2021/04.04-State-Estimation-Copy1.html#4.4.4-Eigenvalues)",
     "section": "4.4.4 Eigenvalues"
    }
   },
   "outputs": [],
   "source": [
    "from scipy.signal import place_poles\n",
    "\n",
    "evals, evecs = np.linalg.eig(A)\n",
    "results = place_poles(A.T, C.T, 3*evals)\n",
    "print(3*evals, results.computed_poles)\n",
    "L = results.gain_matrix.T\n",
    "print(L)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[4.4.4 Eigenvalues](https://jckantor.github.io/cbe30338-2021/04.04-State-Estimation-Copy1.html#4.4.4-Eigenvalues)",
     "section": "4.4.4 Eigenvalues"
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(2, 1)\n",
    "ax[0].plot(h.logdict[\"Time\"], h.logdict[\"Ts\"], h.logdict[\"Time\"], h.logdict[\"T1\"])\n",
    "ax[1].plot(h.logdict[\"Time\"], np.array(h.logdict[\"Ts\"]) - np.array(h.logdict[\"T1\"]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[4.4.4 Eigenvalues](https://jckantor.github.io/cbe30338-2021/04.04-State-Estimation-Copy1.html#4.4.4-Eigenvalues)",
     "section": "4.4.4 Eigenvalues"
    }
   },
   "source": [
    "**Study Question:** Rerun the estimator simulation using:\n",
    "\n",
    "* Change parameter value\n",
    "* Change in d\n",
    "* Change sensor offset\n",
    "* "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpages": {
     "level": 2,
     "link": "[4.4.4 Eigenvalues](https://jckantor.github.io/cbe30338-2021/04.04-State-Estimation-Copy1.html#4.4.4-Eigenvalues)",
     "section": "4.4.4 Eigenvalues"
    }
   },
   "source": [
    "<!--NAVIGATION-->\n",
    "< [4.4 Anomaly Detection](https://jckantor.github.io/cbe30338-2021/04.04-Anomaly-Detection.html) | [Contents](toc.html) | [Tag Index](tag_index.html) | [4.6 Lab Assignment 6: Anomaly Detection](https://jckantor.github.io/cbe30338-2021/04.06-Lab-Assignment-Anomaly-Detection.html) ><p><a href=\"https://colab.research.google.com/github/jckantor/cbe30338-2021/blob/master/docs/04.04-State-Estimation-Copy1.ipynb\"> <img align=\"left\" src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open in Colab\" title=\"Open in Google Colaboratory\"></a><p><a href=\"https://jckantor.github.io/cbe30338-2021/04.04-State-Estimation-Copy1.ipynb\"> <img align=\"left\" src=\"https://img.shields.io/badge/Github-Download-blue.svg\" alt=\"Download\" title=\"Download Notebook\"></a>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}