Simulate_Dynamics.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulate Dynamic Systems in Python\n",
    "\n",
    "A step response is a common evaluation of the dynamics of a simulated system. A linear time invariant (LTI) system can be described equivalently as a transfer function, a state space model, or solved numerically with and ODE integrator. This tutorial shows how to simulate a first and second order system in Python."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### First Order System Simulation\n",
    "\n",
    "Consider a first order differential equation with constants $K_p=3$ and $\\tau_p=2$, input $u$, and output response $y$.\n",
    "\n",
    "$\\tau_p \\frac{dy}{dt} = -y + K_p u$\n",
    "\n",
    "Three methods to represent this differential equation is as a (1) transfer function, (2) state space model, and (3) semi-explicit differential equation. Source code is included below to simulate a step response in these three forms.\n",
    "\n",
    "1. Transfer Function\n",
    "\n",
    "$\\frac{Y(s)}{U(s)} = \\frac{K_p}{\\tau_p \\,s + 1}$\n",
    "\n",
    "2. State Space Model\n",
    "\n",
    "$\\dot x = A x + B u$\n",
    "\n",
    "$y = C x + D u$\n",
    "\n",
    "$A = -\\frac{1}{\\tau_p} \\quad B = \\frac{K_p}{\\tau_p} \\quad C=1 \\quad D=0$\n",
    "\n",
    "3. Differential Equation\n",
    "\n",
    "$\\frac{dy}{dt} = -\\frac{1}{\\tau_p} y + \\frac{K_p}{\\tau_p} u$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA6FUlEQVR4nO3dd3RU5dbH8e9OIQkJzQDSqwGkhAChKy2CNAVU7quigl5FEbBXVATFe/HqFS9iQ0VQsWMBBFR6lxqQDkoLICVACgmp+/1jhhgggQAZTpLZn7XOypz+S5s9zynPEVXFGGOM9/JxOoAxxhhnWSEwxhgvZ4XAGGO8nBUCY4zxclYIjDHGy/k5HeBClS1bVmvUqOF0DGOMKVRWr159RFXL5TSv0BWCGjVqsGrVKqdjGGNMoSIiu3ObZ4eGjDHGy1khMMYYL2eFwBhjvJwVAmOM8XJWCIwxxst5rBCISKCIrBCRdSKyUURG5rCMiMhYEdkhIutFpKmn8hhjjMmZJy8fTQE6qWqiiPgDi0Vkpqouz7ZMNyDMPbQE3nV/NcYYR6WmQmYmZGS4vp4azhzPzISyZSEw8PT1t28/fVnV07+eeq0KjRpBUNDf66anw9Klf89XhZIloVkzz3yvHisE6urfOtE96u8ezuzzuhfwiXvZ5SJSWkQqquoBT+UyxuQvVdcbl7//6dNjYmDfPkhJ+XtITf37a/bXzZpB+/anr//2hwn8tuYEkloSTS1OejqcyDhOgs9uSClJsaSapKdDWkY61941nw7t4bpa12Wt33ngXH7f9RcBMVFw4krS0+Fk6WhSrlyCz6HGyN5ryMiAdP+j9Hnlfdo1L82g5oNc31BGBtVveoNjSYfwXT4En+RQfMjEJ2wavtXn4bO9Kz67r3VNK7ONTk++Q6uwmjzY7IGsd/neQ54hg2R85ozCJ6MYguITPgmpuBqfdbfj81cTBEUqrqLpXZ/QpnoEdzW+EzIziY9L5vm3nkbSiyFzRiMoDa5Wmi1qDKGh+f47FE8+j0BEfIHVwFXA26r69BnzpwOjVXWxe3wO8LSqrjpjuYHAQIBq1ao127071/sijDE5yMyEEycgPt41JCS4viYmuoYrroDu3U9f5/PvjvHxlH2knyiNHq9CUhIkpsVzpOyPpCUHwsa+JCdDSorS6oF/0b7jHp6IHEpZnxA4eZKBo79k3V8/E7ClGwG7WxNACsWu2EhA5JsUO1qdYqv+STFSCSCF4E4vIJrMcy2fIEB9IC2NiXOncdJ3O/5/tsfvWHX8ScO/5B/4VVqB//FK+P3VCD/S8ZMU/Kosphg+tCjfxFWVMjLYuu9PRJLwi6+Ib7o/fqTj65/gGlID8E0PwJcMfEnHxzcFXwW/TPc37+NDKplkCmRmBJKBH5n4kOmXQqZ/ChlpwWSmh6AIGb7pZAYfoZh/ABVKVAQfHxBh29E/yfRR9GgdMtW1vpaMITPoOJlx1dGToShCZtAx9Io/uaL4FYSVrQsipGZk8tuB5aj6oHuvQRFKlBCazX4VWrS4qL8BEVmtqpE5zfPoncWqmgFEiEhp4HsRaaiqG7Jny2m1HLYzHhgPEBkZaU/SMV4pMxOOHYPY2NOHo0ehWEAGtw44TqBfIMHFggEY//Vunnl/DimHK3Py984U5wQlfQ9S8rohlEhPoeSc5wghkRASqdD+dRa/tZn76t1OTb+ycOIEFZavZGDKKoIPVaf4kRoEkUxxn+MEldpG8VRfgk4MJJCTBJBC2geQ/LESUnIKBAVDYCBPHj/GMb/DpCQeIIWFpBBASnIiKfv/Ii05gxRWkEoxUilGZnIGSf5Khp8vFC8J/v7ElC7PPr8DpAXWI51WpOFPWtAm0qr/RUZAM9L+uo80/MmQDMp3fo6q1QNp0W0M+PmBnx8PDR/PzqSdZKwaTMbRq0nHj4wKS8ms+xPpuzuSsb03GfiSEXic60eNpV3zcjx6zROuN3KgVu/3+ev4Mfx/vxe/1LL4+EBmlcVkll9LsYNtCTjaFBEg5C8GvDKDa5tUoHvY39W0a89vOJmaQcl9N+FLMXx8IPmKlWQU30dwXDMCUqoiAmmB+3n439E0DasEFSIASD2RwYO3zUfwpUzNDohA/frQ7OJqwHl5tEVw2o5EXgROqOrr2aa9D8xX1S/c41uBDuc6NBQZGanWxYQpKtLT4eBBOHIEGjc+fd6y1YkMenk1x48JJza141hsJiX1OGU73U+o/x6umPMEoenJlOEYFRpNICRoPd3Ltaa2hMLx4xzYuYe0hD2USvYjJC2TJIoTT0kSQveTEAAJ+zuSSAgJlCCl/mKOldrDDU1vpU71plC8OBOW/M7Mo9NJ3tuJE5vuJoniJAekktT+NZJPViR54b9JJogUAoi4ZxJ33RdHv0b9KBfs6s5m+Fub+XLWTkql1aNUZi0CAsA3MIn04D0E+YZwhV8VihWDYsWgRdtk+t7kj6/4IuL6fDh7tus4u7+/673d3//v19nHfX2hdm2oVu30n9/u3a5DT76+rsHPz/XVx+fv16eGgADX16LsXC0CjxUCESkHpKnqcREJAn4BXlXV6dmW6QEMAbrjOkk8VlXPWfOsEJjCIjMT9u93vSHt3g1797qGzYc3syfhTxI2t+bo7pKU1UNUqfcBbaI+pWfJSK4LbgQHD7L79y38setXyscFUO5YKa7gKImEEFsqkSMh6Rw92J6j6VWIJZT4Or9xos46ujW/nfYRvaB0aaZEH+Kxn74k7kgLEtY8TlCwL6VKgdb/muBiwVRJ7kaJEB9KlIAqVx3n+WG+BBcLxkdcn4j/+ANWr4bixSE42PU1KOjvr6eGgICsD9GmAHPq0FBFYJL7PIEP8LWqTheRBwBU9T1gBq4isANIAu72YB5j8t2JE7Bjh2vo1k1Rv6SsQzN/7Eqg/QP3Uz1jP9XmDqIae6jDXqKqfU6V9FgqHwylnMYRSygH9vhx4Od9FK/hD81rQa1aZNSoxyu/HODwyTAOL3ibI5QluJQ/QeE/UaKED9VpR/nSwYSGwhWhyr+HC5LtYGu3SGh24y2UKgUlSrg+Bbv8I4fvpPRZU2rXdg2m6Ltsh4byi7UIzOWm6jp8s3EjbN4Ma7ccI3rfJvZvqcTxLRWozR+EXfkN9SJG0y69It1868LOneiePRz0S2F3adizvzd7tBZ7qUpMwwXEVN7PvvXP8NeBnpQp50+5WvsZ9dZe2tavSfng8gCkpcHChVCunGsIDXUdRjHmYjh2stiYwq7d7ctZc2ANgcuvo8HJA9RnE+FVP+T//NdQd18pruQkO6nJ9uPl+SMulU21hG79HoaaNZEaNWh3x3foyTI0bB1FzaoBVK8Kbas8QpUqULkyVKhw6s29knv4m78/REU58V0bb2OFwHittDT4/XdYsQKWr0ijRueZlKq8hUdLXQ/R0bB+PS8t/oh6R+MISi3ORhqzkQZsSavHzw3i2Vr8LvZtfpZqtfyoHZbO/Q8l0Pv6MqftY9uUfs58c8ZcACsExmscPgxLlsDCpSf5dfNytu8+QtjvdWnBClrzG81+/pCrDyvptT/Cr0kzaNyYac1v4Y7EBI7ueJLwcpHUrw9XXw0d6kHdulCz5qkbqfyAMudJYEzBZIXAFGlLlsD4yYdZsjKRw6uuoA1LaRP8PW+W/4Dm+4QDhLGCFqykOTMjb6b+PyvxcKdhXBlyJQAP74Ghma5LE+3KGFNUWSEwRUZm5hlv1sePs/Kjl2iyaQwP7yhNGOmspDlLTrTm9avq81utNpTa+R4tIn1p3hxubfMQbdqcvs0zr003piiyQmAKtb17Yfp017Dp5Pf0aDGSuw9WpNmGWNi8mX71w3mjuB+DIyJZv/gnIloU45pr4IE2rzCxFVx5pdPfgTHOs0JgChVV1yWcX09J5bPF8zi+vBzd4jfSn2l09v2JnZuS2NE8gWavfght2hDqX4zgV1P5d+sAWrY8vYdHY4yLFQJTKGzbBl9+6RoSN++hT/3b+OjEUiJOBDCH7kynJw8Ve5Gbn9/CK/+MgsDSgOtOxuefDXA0uzEFnRUCU+C1uXURu469Td/fKvFR3DLC2M7UmCa81rYa80KfpW3oA9xwAzzfA2rVauB0XGMKHSsEpuA6cQK++443fnueun/t4cfQeoyMG8OyoCg6d/Gj303C5O5QqpTTQY0p3KwQmAJh/34YPx5WxU2jashwhm0pT9XZK6BtWxjwEFXmx9LmygHc26cOU3q6OkEzxuQPKwTGUcuWwdixMOubBG7L+JR/Bb9MUMBfLOzWiH6bN0OFCrQC9j9in/yN8RQrBOayy8iAqVPhxfdXk+T7LENW+DAuYyXz6MgjGW8R+tgRJjzRDwJKZK1jRcAYz7FCYC6btDT49FMYPRqu2L6cF0sPpl3qGj68qhIRR6KpcU1VhgyBPn2sl01jLicrBOaymPJ9BoP+9y3hOzfwwZ6FVGMPb8QP5e5ru9E9bDDTJlUkIsLplMZ4JysExvNUSVj8Cj/+8SKlk/wYJR8wq8QdDHzQj60PQcWKTgc0xrtZITAeExMfQ5X1u+DZZ+kfG8uAek34PuU+nux6O+8O8aNkSacTGmPACoHJZ7//Dg8/FU9oRHfu/WEFlZIq4PPSy8gdd/Dsdl/ergIhIU6nNMZkZx3rmnxx7Bg8+CBENT7CTbOe5d03V7Cwli9Lfp0A/fuDry/16lkRMKYgskJgLokqfDwxkxq9PoAPX2OD1kfxoV7QUlo/vpdr61zndERjzHnYoSFz0bZvh/vug6ST/Ziz50tOlKhAp6NzqdS5IfPfgIYNnU5ojMkLaxGYC5aeDq+/Dq0bJdJnwcP8uGI2/wsP5eYGY3hpSkN+/tmKgDGFibUIzAXZsQN6D17JlfvfY0XKPBZxLY1lC3eFl2LXcD87B2BMIWSFwFyQvXt3Myi2NTfsyeD+cv/hr8pPMusjaNrU6WTGmIvlsUNDIlJVROaJyGYR2SgiD+ewTAcRiRORaPcw3FN5TD5YvZqOg66nXlojIlo8QIuhQ1ixwoqAMYWdJ1sE6cDjqrpGREoAq0XkV1XddMZyi1S1pwdzmEuQmAgzd31Nq6+XUfW9yTB2LB363srirVC/vtPpjDH5wWOFQFUPAAfcrxNEZDNQGTizEJgCSBX++18Y99XHjI2/hyMpAVyxZA3BYfXxxYqAMUXJZblqSERqAE2A33KY3VpE1onITBHJ8TmDIjJQRFaJyKrDhw97MqrB9WCwW2+Fr55cyfxVI9gdUpHFn4ym+FVXOx3NGOMBHj9ZLCIhwBTgEVWNP2P2GqC6qiaKSHfgByDszG2o6nhgPEBkZKR6NrF327ULuvRfS+tFa5nB0wxkPPv9b2RaPV9EnE5njPEEj7YIRMQfVxGYrKrfnTlfVeNVNdH9egbgLyJlPZnJ5G7JEmhy+wQGn4zkuaDHac8CKg7qw6KFvpQv73Q6Y4yneKxFICICfARsVtU3clmmAnBQVVVEWuAqTLGeymRy99VXMOSueL70eY/MykqrNrcxuu/VDBzodDJjjKd58tBQW+BO4HcRiXZPGwZUA1DV94BbgEEikg4kA7eqqh36uYxU4dVX4e1n9zKHHizmGkboZ0wfWYc2bZxOZ4y5HDx51dBi4JxHlVV1HDDOUxnMuWVmwqDHjvHbju4s8/mTMZlPMaPuYyyfIdSq5XQ6Y8zlYncWezFVSNr7EL/MW86DHa/iSMZjLP1OKFPG6WTGmMvJOp3zYr6zfmLSop95pNMNnCw/l59nWREwxhtZi8AL7YvfR8XpC/B57DF8pk3jo/CWFCsGvr5OJzPGOMEKgRdJTob/Td7K3nmtGDUzhZLzluLbOIIgp4MZYxxlhcBLJCdDr15Qbs1nvJYSx2PPRPJO/bpWBIwxdo7AG6SlQd++UPbXz3kt9iOuC/qSf7SaR5C/lQFjjLUIiryMDOh7z358Fn/LG/yLKOZw65AGdItyOpkxpqCwQlCEqcK9Q48Sl9CGr9L20CN0Ir0GNmC4PfXBGJONFYIibORI2PD+Zn7y/4t/dK1Oo8o9eOUVp1MZYwoaKwRF1Mcfw6cj/2Axt3Bv6heElriWD98KtR5EjTFnsZPFRdAvvyhP/WciM+jOSF4ktUsfvpxQFh/7bRtjcmAtgiImMRFufXU0U08M44c6V7M86AEWfQP+/k4nM8YUVPYZsYgJCVaWpv5KTIg/r5V/kZ9+ghIlnE5ljCnIrEVQ1Lz6KvVSEombsIt5wZWoXNnpQMaYgs4KQRFxNPkoCd99SfW33oIVK2hZuZLTkYwxhYQdGioCfpiazsOv9aT4/YNZ+b+nsWaAMeZCWIugkFu3Du69PZYlvht5+foSPNG5t9ORjDGFjLUICrH4eOh7i/LWiUdYFH8L8//aQtWS1ZyOZYwpZKxFUEipwl0Dj9FpxxfUZxNRxZez+KMgu2HMGHPBrBAUUuPeSWPPyfZ8ELCFtrKIsR8EUa+e06mMMYWRFYJCaONGePGVNaxI3siQqNK0qV6b2293OpUxprCyQlDIpKTA7bfDawfGs8D/H6xNeIr1b5R1OpYxphCzQlDIPP881F7/He1ZQGuftcx/twSBgU6nMsYUZlYICpG5c+Hztf1Z4/cDvdJnMfy1EjRo4HQqY0xh57HLR0WkqojME5HNIrJRRB7OYRkRkbEiskNE1otIU0/lKQpe/2YuE3Z9wrutkgjsWZ0hQ5xOZIwpCjx5H0E68LiqXg20AgaLSP0zlukGhLmHgcC7HsxT6E1v8gd1U6vzfuC7fDG+kl0qaozJFx47NKSqB4AD7tcJIrIZqAxsyrZYL+ATVVVguYiUFpGK7nVNdjEx+Dw3jBpz57L36kb42UE9Y0w+uSx3FotIDaAJ8NsZsyoDe7ONx7innbn+QBFZJSKrDh8+7LGcBdWKmN9IvfduGDoUGlkRMMbkL4+/pYhICDAFeERV48+cncMqetYE1fHAeIDIyMiz5hdlX087xqIvrid4fSL+H4+hjtOBjDFFjkcLgYj44yoCk1X1uxwWiQGqZhuvAuz3ZKbCZPVqGHLHLjakJfHQ4IZ8XuHMUyzGGHPpzlkIRCQQ6AlcC1QCkoENwE+quvE86wrwEbBZVd/IZbGpwBAR+RJoCcTZ+QGXjAx44AH4T/z/+JQhxG4Zjo9YH4HGmPyXayEQkRHADcB8XMf2DwGBQB1gtLtIPK6q63PZRFvgTuB3EYl2TxsGVANQ1feAGUB3YAeQBNx9Sd9NEfL++xC0aiFRzKFJsU0sH2PPmzTGeMa5WgQrVXVELvPeEJHyuN/Uc6Kqi8n5HED2ZRQYfL6Q3ubgQXjmyzdYVuJlHkkZw9BhJbjqKqdTGWOKqlwLgar+BCAiDVV1Qw7zD+FqJZh89siTyQzyeYE9VyaxrFIok592OpExpijLy8ni90SkGDAR+FxVj3s0kZdbuhQWf3qEaPGnReQjTHzuButLyBjjUec9+6iq1wD9cF3ds0pEPheRzh5P5oUyM+HRR+HfPMu7OoSm1cfQpYvTqYwxRV2eLh9V1e0i8jywChgLNHFfFTQsl8tCzUWYPFmRHZ/SkXmEF9vKqv84ncgY4w3OWwhEJBzX1Tw9gF+BG1R1jYhUApYBVgjywYkT8PAHk5kR1J9hDbtw/w0h1KzpdCpjjDfIS4tgHPABrk//yacmqup+dyvB5JMnQmfht0P4scz/sfdZp9MYY7zFeQuBqrY7x7xP8zeO9wqWJIatXsj60d/wbeU+lLDbBowxl8m5biibhqt/n1mqmnbGvFrAAGCXqk7waEJv8dZb0Lw54Xfd7HQSY4yXOVeL4D7gMeBNETkKHMZ1Z3EN4A9gnKr+6PGEXuCD2f+h/3/+TbEly52OYozxQue6oewv4CngKXc30hVx9TW0TVWTLk+8ou+Vt3fj/9WzfFIjkw7l/bAbiI0xl1teLx/dBezyaBIvtHo1vP1IGr8TwG339+We0lYGjDGXn3Vn6aBhw+CF9P8yIX0wJQ5Mwsd+G8YYB9izrhyyaBHs+GUHffmGq2UrC0c5ncgY463y9BlURIJEpK6nw3iTJ15bz3P1WjAuNIqe/UO5+mqnExljvNV5C4GI3ABEA7Pc4xEiMtXDuYq0JUvgcOK/uXH3McY2KcPzdlueMcZBeWkRjABaAMcBVDUa1yWk5iKNHAnD5gfyTunuXFdtOLVrO53IGOPN8nKOIF1V41x9zJlLtWwZbPt1F32YSt2/trPsmSucjmSM8XJ5KQQbROR2wFdEwoCHgKWejVV0Pf+vIwyTl3lPH6D7HVcQFuZ0ImOMt8vLoaGhQAMgBfgCiAce8WCmImvVKtie8U9uDpjImKqteO45pxMZY0zeOp1LAp4DnhMRXyBYVU96PFkRFBuXwlMHlvBRE6XtVRHUteuwjDEFQF6uGvpcREqKSDCwEdgqIk96PlrRc31EIoN3Z1LmjlmMfraq03GMMQbI26Gh+qoaD/QGZgDVgDs9GarIeucdpE8f7nuwi903YIwpMPJysthfRPxxFYJxqpomIurZWEXP8u3zaDHuLXzmL3A6ijHGnCYvLYL3cXU4FwwsFJHquE4YmzzatP0EU566gVmhx9hTKdjpOMYYc5rzFgJVHauqlVW1u7rsBjqebz0RmSAih0RkQy7zO4hInIhEu4fhF5G/wEtMhLbXHWDQ3AwmdqhDpWA7N2CMKVjy8vD6AOBmXHcTZ1/+pfOsOhHX844/Occyi1S15/kyFGaTJkHUnmgO0JRV82ZiN+YZYwqavJwj+BGIA1bjupcgT1R1ofuBNl4rM9P1BMoP+B9v8BiPDymJr6/TqYwx5nR5KQRVVLWrh/bfWkTWAfuBJ1R1Y04LichAYCBAtWrVPBQl//3yCxRnONXlT+aF9GJSf6cTGWPM2fJysnipiDTywL7XANVVtTHwFvBDbguq6nhVjVTVyHLlynkgime88sEGhga9zNvtjnLHvemEhDidyBhjzpaXFsE1wAAR2Ynr0JAAqqrhl7Jj970Jp17PEJF3RKSsqh65lO0WFNu2wdb5cfRO8OeqdneyYnCg05GMMSZHeSkE3TyxYxGpABxUVRWRFrhaJ7Ge2JcTxo2DgUfn8S39aVP8fetq2hhTYOWlr6HdItIYuNY9aZGqrjvfeiLyBdABKCsiMcCLgL97m+8BtwCDRCQdSAZuVdUicaNafDx8OiGNDbxLN2byxkN2pZAxpuDKy+WjDwP3Ad+5J30mIuNV9a1zraeqt51n/jhcl5cWORMnZXB9rZvYvqUSGWHhREU5ncgYY3KXl0ND/wRaquoJABF5FViG6wSvOYMqjJn6CxMSpvP2deV5sIfavQPGmAItL1cNCZCRbTzDPc3kQAQm/tOPRoeLMz/tGe68035UxpiCLS8tgo+B30Tke1wFoBfwkUdTFXLtV/0MQ4by1yuP2g1kxpgCLy8ni98Qkfm4LiMFuFtV13o0VWF28iR88gksXWpFwBhTKOTl0NApAih2WChXGZkZfPBcN46EVSGzdi2n4xhjTJ7k5Qllw4FJQBmgLPCxiDzv6WCFjSo8PO5nwqbM58W6+5yOY4wxeZaXcwS3AU1OPadYREbj6h5ilCeDFTZLl8LsZ0ryfHoIG4NG4iMX0tgyxhjn5OXdaheQvX+EAOAPj6QpxCZMgH8m/8jEtMFclfKA03GMMSbP8tIiSAE2isivuM4RdAYWi8hYAFV9yIP5CoUTJ2DKV+ls5jM6MJ9J9zqdyBhj8i4vheB793DKfM9EKby++05p27Abu34PJrNaDVq2dDqRMcbkXV4uH5106rWIlAGqqup6j6YqZN6asoYnEmYzqW0wd3UQ7EZiY0xhkpe+huYDN7qXjQYOi8gCVX3Ms9EKhz17YMfMslyfHsT9gS/y+13FnI5kjDEXJC8ni0u5nx1wE/CxqjYDrvNsrMLj00/h1tSfmJl5I81Dn6RKFacTGWPMhclLIfATkYrAP4DpHs5TqKi6Hk7fn0lMoj8DBjidyBhjLlxeCsFLwM/AH6q6UkRqAds9G6twWLYMfMv3p2qxzSwLjaR3b6cTGWPMhcvLyeJvgG+yjf8J3OzJUIXFxK+P0F8/5bNm0KNRJsWLO53IGGMuXF66mKgjInNEZIN7PNy6mHD576gQHthallXVHmJw/yudjmOMMRclL4eGPgCeBdIA3JeO3urJUIVFiXWrKV2+LF9/MYY2bZxOY4wxFycvhaC4qq44Y1q6J8IUOl98Abffjt04YIwpzPJSCI6ISG1c3UsgIrcABzyaqhCYsGI8CZM/5s8uLZyOYowxlyQvXUwMBsYD9URkH7AT6OfRVAXcsmWw4KPXaFg8iZgSidiTB4wxhdl5WwSq+qeqXgeUA+oBHfj7aWVe6bXXoMMXdfmmdHvSN3d3Oo4xxlySXAuBiJQUkWdFZJyIdAaSgP7ADlw3l3mluDiY89NJeiUtZfKmz2ncIPD8KxljTAF2rkNDnwLHgGXAfcBTQDGgt6pGez5awfT999ApdSbraEyFJpWoW9fpRMYYc2nOVQhqqWojABH5EDgCVFPVhLxsWEQmAD2BQ6raMIf5AvwP6I6rtTFAVddcYP7L7sOpvzO08iN8efImbrvN6TTGGHPpznWOIO3UC1XNAHbmtQi4TQS6nmN+NyDMPQwE3r2AbTvi+HFYe3wy1x/Zw3dNj9K3r9OJjDHm0p2rRdBYROLdrwUIco8LoKpa8lwbVtWFIlLjHIv0Aj5RVQWWi0hpEamoqgX20tRp0+C6lTVYG1iT0MxB1KjhdCJjjLl0uRYCVfX18L4rA3uzjce4p51VCERkIK5WA9WqVfNwrNxNmQI3Jy5lCo/RP6qVYzmMMSY/5eWGMk/J6XZczWlBVR2vqpGqGlmuXDkPx8pZYiLMnZVKT6bzPX242brdM8YUEU4WghigarbxKsB+h7Kc18yZ0KZRX7YEX0Hp8FDq1HE6kTHG5A8nC8FU4C5xaQXEFeTzA5Om/snNmVOZ0jKGPn2cTmOMMfknL11MXBQR+QLXXchlRSQGeBHwB1DV94AZuC4d3YHr8tG7PZUlP4wZWYryjUvQ5/oHGXuL3URmjCk6PFYIVPWcV9m7rxYa7Kn957ew3euhThhzvx3tdBRjjMlXTh4aKlymTMHOEBtjiiIrBHkwdcuPxH/zGfuirMtpY0zRY4XgPHbuhO8/+Q8H0+P4xX+P03GMMSbfWSE4j3ffhaofVGNW+Xqc3NDT6TjGGJPvrBCcx7RpcMORHXy35R2qlCnvdBxjjMl3VgjOYccOiNuyn9r8werAa4iKcjqRMcbkPysE5/Dj1Ax6ht3PrKBI2nf2o3hxpxMZY0z+s0JwDpMXrOAGmc609mu58Qan0xhjjGdYIcjFsWOwbWlx2v/hz6zjD9CzZ0595BljTOHnsTuLC7tZs6DTkd2s4hrCMl6mYkWnExljjGdYiyAX06bBDUxjGjdw441OpzHGGM+xQpCDjAyYvn4hPXy/Y7pfFD3t9gFjTBFmhSAHq1ZBjaovkFTyKIevWUXjxk4nMsYYz7FCkIMrr4TBacWZXa0c3epcj9h5YmNMEWYni3NQowbcn5kCoyYwsHtlp+MYY4xHWYsgJwkJsHIldOiAj/2EjDFFnL3N5WDlZ/8hNbIphIQ4HcUYYzzODg2dYfuhGFZOGMWP5QIYkZmOn4/9iIx3SktLIyYmhpMnTzodxVyAwMBAqlSpgr+/f57XsXe5bDIzoW3Ho6zcG8hDzduRkuyHX7DTqYxxRkxMDCVKlKBGjRqIXTFRKKgqsbGxxMTEULNmzTyvZ4eGslm/HkpvCsQ3IZTFX80k0J5Rb7zYyZMnCQ0NtSJQiIgIoaGhF9yKs0KQzaxZ0JVZzKIr11/vg6+v04mMcZYVgcLnYn5nVgiy+XHufrr6fc8sutK1q9NpjDHm8rBC4HbiBESnvcO1PvOZ02Y1Xbo4ncgY7xYbG0tERAQRERFUqFCBypUrZ42npqbm6762bNlCREQETZo04Y8//riobQwYMICaNWtmZRw7dmy+ZvQkO1nstngxtIjbxZYrfCgZ3IoKFZxOZIx3Cw0NJTo6GoARI0YQEhLCE088kTU/PT0dP7/8eQv74Ycf6NWrFyNHjszzOhkZGfiecfz4tdde45ZbbsmXTJeTR1sEItJVRLaKyA4ReSaH+R1EJE5Eot3DcE/mOZfZsyFqbS3m/PU4N9bv7lQMYwqsESNAJG/DwIFnrz9w4OnLjBhx4RkGDBjAY489RseOHXn66adZsWIFbdq0oUmTJrRp04atW7cCMHHiRG666Sa6du1KWFgYTz31FOB68x4wYAANGzakUaNGjBkzhhkzZvDmm2/y4Ycf0rFjRwA+++wzWrRoQUREBPfffz8ZGRkAhISEMHz4cFq2bMmyZcvylHnWrFk0bdqUxo0bE+V+3u2IESO455576NChA7Vq1XK89eCxFoGI+AJvA52BGGCliExV1U1nLLpIVR3v33POHHiLObzISB6Oyvv1t8aYy2vbtm3Mnj0bX19f4uPjWbhwIX5+fsyePZthw4YxZcoUAKKjo1m7di0BAQHUrVuXoUOHcujQIfbt28eGDRsAOH78OKVLl+aBBx7IanFs3ryZr776iiVLluDv78+DDz7I5MmTueuuuzhx4gQNGzbkpZdeyjHbk08+yahRowD49NNPqVChAvfddx8LFy6kZs2aHD16NGvZLVu2MG/ePBISEqhbty6DBg26oGv/85MnDw21AHao6p8AIvIl0As4sxA4LjYWtm/ZRTjrWe7Tlu/aO53IGJObvn37Zh2SiYuLo3///mzfvh0RIS0tLWu5qKgoSpUqBUD9+vXZvXs3DRo04M8//2To0KH06NGDLjmcDJwzZw6rV6+mefPmACQnJ1O+fHkAfH19ufnmm3PNduahoWnTptGuXbusa/qvuOKKrHk9evQgICCAgIAAypcvz8GDB6lSpcrF/lguiScPDVUG9mYbj3FPO1NrEVknIjNFpEFOGxKRgSKySkRWHT58ON+DzpsH7dp3YmXVVGp33EnJkvm+C2MKvREjQDVvw/jxZ68/fvzpy1zMoSGA4OC/7/J84YUX6NixIxs2bGDatGmnXT8fEBCQ9drX15f09HTKlCnDunXr6NChA2+//Tb33nvvWdtXVfr37090dDTR0dFs3bqVEe6wgYGBZ50XOBdVzfVyzpzyOcWThSCn717PGF8DVFfVxsBbwA85bUhVx6tqpKpGlitXLn9TAkfikuh8+BBzamXSrVWtfN++McYz4uLiqFzZ9fly4sSJ513+yJEjZGZmcvPNN/Pyyy+zZs2as5aJiori22+/5dChQwAcPXqU3bt3X1S+1q1bs2DBAnbu3Jm1rYLIk4eGYoCq2carAPuzL6Cq8dlezxCRd0SkrKoe8WCuszzwz+Lo/2qx4s6RVOxstxMbU1g89dRT9O/fnzfeeINOnTqdd/l9+/Zx9913k5mZCcC///3vs5apX78+o0aNokuXLmRmZuLv78/bb79N9erVLzhfuXLlGD9+PDfddBOZmZmUL1+eX3/99YK342mieuaH9HzasIgfsA2IAvYBK4HbVXVjtmUqAAdVVUWkBfAtrhZCrqEiIyN11apV+Rv20CGoUweOHIF8uhzNmMJu8+bNXH311U7HMBchp9+diKxW1ciclvfYu56qpovIEOBnwBeYoKobReQB9/z3gFuAQSKSDiQDt56rCHhK5uzZ+LRvb0XAGOOVPPrOp6ozgBlnTHsv2+txwDhPZjifXcd3seDNu/FrHEE/J4MYY4xDvLqLCVW4Z8QSrt2eyvxq9hAaY4x38upCsGMH7PjftYTEhfLlhLdw3zxojDFexasLwYIF0I6FLNQOdGhY37qdNsZ4Ja8uBAsXQnsWsID2tGvndBpjjHGGVxeCmbu/oV3Jz1lY6yTtrVsJYwqcV155hQYNGhAeHk5ERAS//fYbAG+++SZJSUnnXT+vy2W3fPlyWrZsSUREBFdffXXWXcVFmddeL7lrF/gF/EC51CS2VlWaNHE6kTEmu2XLljF9+nTWrFlDQEAAR44cyXoOwZtvvskdd9xB8eLFz7mNvC6XXf/+/fn6669p3LgxGRkZWT2aFmVe2yJYuBCuXdiexf7NaF7yFhzq9M+YQkNGCjLy9J5jbvjiBmSkMG3rtKxp41ePR0YKA6f93Rf1/oT9yEih0n8r5Xl/Bw4coGzZsll98pQtW5ZKlSoxduxY9u/fT8eOHbO6jR40aBCRkZE0aNCAF198ESDH5X755Rdat25N06ZN6du3L4mJiWft99ChQ1SsWBFw9QFUv359wNV19J133kmnTp0ICwvjgw8+ACAxMZGoqCiaNm1Ko0aN+PHHH7O29cknnxAeHk7jxo258847ATh8+DA333wzzZs3p3nz5ixZsiTPPxOPUdVCNTRr1kzzwz33qL7FYH2M13XUqHzZpDFFyqZNm04bZwTKCE6b1vPznsoIdOqWqVnT3l/1vjICvW/qfVnT9sXvU0agFV+vmOf9JyQkaOPGjTUsLEwHDRqk8+fPz5pXvXp1PXz4cNZ4bGysqqqmp6dr+/btdd26dWctd/jwYb322ms1MTFRVVVHjx6tI0eOPGu/I0eO1NKlS2vv3r31vffe0+TkZFVVffHFFzU8PFyTkpL08OHDWqVKFd23b5+mpaVpXFxc1j5q166tmZmZumHDBq1Tp07W/k9lvO2223TRokWqqrp7926tV69enn8meXXm705VFVilubyveu2hoYUL4REWMIn+jLHzA8acl7549k3/026bdta0gc0GMrDZ6U+mqVSiUo7rn0tISAirV69m0aJFzJs3j//7v/9j9OjRDBgw4Kxlv/76a8aPH096ejoHDhxg06ZNhIeHn7bM8uXL2bRpE23btgUgNTWV1q1bn7Wt4cOH069fP3755Rc+//xzvvjiC+bPnw9Ar169CAoKIigoiI4dO7JixQp69OjBsGHDWLhwIT4+Puzbt4+DBw8yd+5cbrnlFsqWLQv83QX17Nmz2bTp79744+PjSUhIoESJEhf088lPXlkI9u+Ho6VeoZrvdjaGlsTd7bgxpoDx9fWlQ4cOdOjQgUaNGjFp0qSzCsHOnTt5/fXXWblyJWXKlGHAgAGndUd9iqrSuXNnvvjii/Put3bt2gwaNIj77ruPcuXKERsbC3BWl9IiwuTJkzl8+DCrV6/G39+fGjVqcPLkyVy7oM7MzGTZsmUEBQVdwE/Cs7zyHMGCBXBtufdZViOFuq0Ok61bcGNMAbF161a2b9+eNR4dHZ3VA2iJEiVISEgAXJ+og4ODKVWqFAcPHmTmzJlZ62RfrlWrVixZsoQdO3YAkJSUxLZt287a708//YS6uzzbvn07vr6+lC5dGoAff/yRkydPEhsby/z582nevDlxcXGUL18ef39/5s2bl9VldVRUFF9//XVWETnVBXWXLl0YN+7vnnVOPZfZSV7ZIujeHap/0pS1VWswsKc1B4wpiBITExk6dCjHjx/Hz8+Pq666ivHuJ94MHDiQbt26UbFiRebNm0eTJk1o0KABtWrVyjr0k9NyEydO5LbbbiMlJQWAUaNGUadOndP2++mnn/Loo49SvHhx/Pz8mDx5ctbDaFq0aEGPHj3Ys2cPL7zwApUqVaJfv37ccMMNREZGEhERQb169QBo0KABzz33HO3bt8fX15cmTZowceJExo4dy+DBgwkPDyc9PZ127drx3nvv4SSPdUPtKfnWDXWzZjB2LGT7ozHG/M26oT7diBEjsp5rXNBdaDfUXnloiPh42LoVOzlgjDFeemjoh4nPck296oT4ZGLPIzPG5EVRvsPY61oE23cls+6H9/gwaBOpGalOxzHGGMd5XSHo/8+TtF5fg41Bndm9raTTcYwxxnFeVQjS0mDdkpK0jD3CzF8/JzTU6UTGGOM8ryoE0dFQO2UjB6hIiRplqZT3bk+MMabI8qpCsGSJ0ubKd1jqF2FXjRpTCMTExNCrVy/CwsKoXbs2Dz/8cFYPpPPnz6dUqVI0adKEunXr0q5dO6ZPn5617ogRI6hcuTIRERFZw/Hjx0/b/q5du2jYsOF5c/zrX//K1+8rN9HR0cyYMeP8C+YzryoEv6zeTtsr32dpx5m0bl247p8wxtuoKjfddBO9e/dm+/btbNu2jcTERJ577rmsZa699lrWrl3L1q1bGTt2LEOGDGHOnDlZ8x999FGio6OzhlN3CF+o/CwE6enpuc67mEJwru3lldcUAlVYtS6JNn8GsSSjE23bnt0HiDGm4Jg7dy6BgYHcfffdgKvfoTFjxjBhwoQcHzYTERHB8OHDT+u+4UJMnDiRm266ia5duxIWFsZTTz0FwDPPPENycjIRERH069cPgM8++4wWLVoQERHB/fffT4b7gecfffQRderUoUOHDtx3330MGTIEgAEDBvDYY4/RsWNHnn76aVasWEGbNm1o0qQJbdq0YevWraSmpjJ8+HC++uorIiIi+Oqrrzh69Ci9e/cmPDycVq1asX79esDV2hk4cCBdunThrrvuuqjvNzuvuY9g717w+b0CZQhk329TyENr0BiTXQ4dqF2yc/RssHHjRpo1a3batJIlS1KtWrWs/oLO1LRpU1577bWs8TFjxvDZZ58BUKZMGebNm3fOONHR0axdu5aAgADq1q3L0KFDGT16NOPGjcvqE2jz5s189dVXLFmyBH9/fx588EEmT57Mddddx8svv8yaNWsoUaIEnTp1onHjxlnb3rZtG7Nnz8bX15f4+HgWLlyIn58fs2fPZtiwYUyZMoWXXnqJVatWZRWzoUOH0qRJE3744Qfmzp3LXXfdlZVj9erVLF68OF86r/OaQrB0KbRhKctoTcvWPvh5zXduTD65zN3R5NZ7Z27TT83L7tFHH72gLiGioqIoVaoUAPXr12f37t1UrVr1tGXmzJnD6tWrae7umSA5OZny5cuzYsUK2rdvn9XddN++fU/r1K5v375ZfRbFxcXRv39/tm/fjoiQlpaWY57FixczZcoUADp16kRsbCxxcXEA3HjjjfnWg6lHDw2JSFcR2SoiO0TkmRzmi4iMdc9fLyJNPZVl4dJk2hSbxVJa06aNp/ZijMkvDRo04Mx+xeLj49m7dy+1a9fOcZ21a9deUv9IAdm6Ivb19c3x+Luq0r9//6zzDlu3bmXEiBFnFaEzBQcHZ71+4YUX6NixIxs2bGDatGk5dpt9al9nOlUEs2/vUnmsEIiIL/A20A2oD9wmIvXPWKwbEOYeBgLveirPL9vn0ab8ByzpOMsKgTGFQFRUFElJSXzyyScAZGRk8PjjjzNgwIAcn0G8fv16Xn75ZQYPHpzvWfz9/bM+tUdFRfHtt99y6NAhwNW99O7du2nRogULFizg2LFjpKenZ32Sz0lcXByVK1cGXOcmTsnebTZAu3btmDx5MuC6Sqps2bKULJn/N8J6skXQAtihqn+qairwJdDrjGV6AZ+4n6S2HCgtIhXzO0hmJoSW3Ufjg7BSmtOqVX7vwRiT30SE77//nm+++YawsDDq1KlDYGDgaVfwLFq0KOvy0cGDBzN27FiioqKy5o8ZM+a0y0d37dp1UVkGDhxIeHg4/fr1o379+owaNYouXboQHh5O586dOXDgAJUrV2bYsGG0bNmS6667jvr162cdZjrTU089xbPPPkvbtm2zTjQDdOzYkU2bNmWdLB4xYgSrVq0iPDycZ555hkmTJl1U/vPxWDfUInIL0FVV73WP3wm0VNUh2ZaZDoxW1cXu8TnA06q66oxtDcTVYqBatWrNTj344YIkJpLx7bdsavUPGtU7+9OEMeZ01g31hUtMTCQkJIT09HT69OnDPffcQ58+fS57joLUDXVOZ3POrDp5WQZVHa+qkaoaWa5cuYtLExKC74ABVgSMMR4zYsQIIiIiaNiwITVr1qR3795OR8oTT147EwNkP91eBdh/EcsYY0yh8Prrrzsd4aJ4skWwEggTkZoiUgy4FZh6xjJTgbvcVw+1AuJU9YAHMxljLkBhe4KhubjfmcdaBKqaLiJDgJ8BX2CCqm4UkQfc898DZgDdgR1AEnC3p/IYYy5MYGAgsbGxhIaG5nrdvilYVJXY2FgCAy/skVve+8xiY8w5paWlERMTk+s17qZgCgwMpEqVKvj7+582/Vwni+3+WmNMjvz9/alZs6bTMcxl4DWdzhljjMmZFQJjjPFyVgiMMcbLFbqTxSJyGLiIW4sBKAscycc4nlaY8hamrFC48hamrFC48hamrHBpeaurao535Ba6QnApRGRVbmfNC6LClLcwZYXClbcwZYXClbcwZQXP5bVDQ8YY4+WsEBhjjJfztkIw3ukAF6gw5S1MWaFw5S1MWaFw5S1MWcFDeb3qHIExxpizeVuLwBhjzBmsEBhjjJfzmkIgIl1FZKuI7BCRZ5zOkxsRqSoi80Rks4hsFJGHnc6UFyLiKyJr3U+dK7BEpLSIfCsiW9w/49ZOZzoXEXnU/XewQUS+EJEL61bSw0RkgogcEpEN2aZdISK/ish299cyTmY8JZesr7n/FtaLyPciUtrBiKfJKW+2eU+IiIpI2fzYl1cUAhHxBd4GugH1gdtEpL6zqXKVDjyuqlcDrYDBBThrdg8Dm50OkQf/A2apaj2gMQU4s4hUBh4CIlW1Ia7u3G91NtVZJgJdz5j2DDBHVcOAOe7xgmAiZ2f9FWioquHANuDZyx3qHCZydl5EpCrQGdiTXzvyikIAtAB2qOqfqpoKfAn0cjhTjlT1gKqucb9OwPVGVdnZVOcmIlWAHsCHTmc5FxEpCbQDPgJQ1VRVPe5oqPPzA4JExA8oTgF7gp+qLgSOnjG5F3DqKeuTgN6XM1Nucsqqqr+oarp7dDmupyQWCLn8bAHGAE+Rw2N9L5a3FILKwN5s4zEU8DdXABGpATQBfnM4yvm8iesPM9PhHOdTCzgMfOw+jPWhiAQ7HSo3qroPeB3XJ78DuJ7g94uzqfLkylNPGnR/Le9wnry6B5jpdIhzEZEbgX2qui4/t+sthSCnxysV6OtmRSQEmAI8oqrxTufJjYj0BA6p6mqns+SBH9AUeFdVmwAnKDiHLc7iPrbeC6gJVAKCReQOZ1MVTSLyHK7DspOdzpIbESkOPAcMz+9te0shiAGqZhuvQgFrYmcnIv64isBkVf3O6Tzn0Ra4UUR24Trk1klEPnM2Uq5igBhVPdXC+hZXYSiorgN2quphVU0DvgPaOJwpLw6KSEUA99dDDuc5JxHpD/QE+mnBvrGqNq4PBevc/29VgDUiUuFSN+wthWAlECYiNUWkGK4TblMdzpQjcT0c9iNgs6q+4XSe81HVZ1W1iqrWwPVznauqBfJTq6r+BewVkbruSVHAJgcjnc8eoJWIFHf/XURRgE9uZzMV6O9+3R/40cEs5yQiXYGngRtVNcnpPOeiqr+ranlVreH+f4sBmrr/ri+JVxQC98mgIcDPuP6RvlbVjc6mylVb4E5cn6yj3UN3p0MVIUOBySKyHogA/uVsnNy5Wy7fAmuA33H9vxaoLhFE5AtgGVBXRGJE5J/AaKCziGzHdXXLaCcznpJL1nFACeBX9//ae46GzCaXvJ7ZV8FuCRljjPE0r2gRGGOMyZ0VAmOM8XJWCIwxxstZITDGGC9nhcAYY7ycFQJjciEiodku4f1LRPa5XyeKyDtO5zMmv9jlo8bkgYiMABJV9XWnsxiT36xFYMwFEpEOp567ICIjRGSSiPwiIrtE5CYR+Y+I/C4is9zdhSAizURkgYisFpGfT3XBYExBYIXAmEtXG1c33L2Az4B5qtoISAZ6uIvBW8AtqtoMmAC84lRYY87k53QAY4qAmaqaJiK/43p4zCz39N+BGkBdoCGubgxwL3PAgZzG5MgKgTGXLgVAVTNFJC1bD5aZuP7HBNioqgX6sZjGe9mhIWM8bytQ7tTzkUXEX0QaOJzJmCxWCIzxMPfjUW8BXhWRdUA0heO5AsZL2OWjxhjj5axFYIwxXs4KgTHGeDkrBMYY4+WsEBhjjJezQmCMMV7OCoExxng5KwTGGOPl/h9ChxSzKDVbtwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy import signal\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.integrate import odeint\n",
    "\n",
    "# Simulate taup * dy/dt = -y + K*u\n",
    "Kp = 3.0\n",
    "taup = 2.0\n",
    "\n",
    "# (1) Transfer Function\n",
    "num = [Kp]\n",
    "den = [taup,1]\n",
    "sys1 = signal.TransferFunction(num,den)\n",
    "t1,y1 = signal.step(sys1)\n",
    "\n",
    "# (2) State Space\n",
    "A = -1.0/taup\n",
    "B = Kp/taup\n",
    "C = 1.0\n",
    "D = 0.0\n",
    "sys2 = signal.StateSpace(A,B,C,D)\n",
    "t2,y2 = signal.step(sys2)\n",
    "\n",
    "# (3) ODE Integrator\n",
    "def model3(y,t):\n",
    "    u = 1\n",
    "    return (-y + Kp * u)/taup\n",
    "t3 = np.linspace(0,14,100)\n",
    "y3 = odeint(model3,0,t3)\n",
    "\n",
    "plt.figure(1)\n",
    "plt.plot(t1,y1,'b--',linewidth=3,label='Transfer Fcn')\n",
    "plt.plot(t2,y2,'g:',linewidth=2,label='State Space')\n",
    "plt.plot(t3,y3,'r-',linewidth=1,label='ODE Integrator')\n",
    "plt.xlabel('Time')\n",
    "plt.ylabel('Response (y)')\n",
    "plt.legend(loc='best')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Second Order System Simulation\n",
    "\n",
    "1. Laplace Domain, Transfer Function\n",
    "\n",
    "$\\frac{Y(s)}{U(s)} = \\frac{K_p}{\\tau_s^2 s^2 + 2 \\zeta \\tau_s s + 1}e^{-\\theta_p s}$\n",
    "\n",
    "2. State Space Form\n",
    "\n",
    "$\\begin{bmatrix}\\dot x_1\\\\\\dot x_2\\end{bmatrix} = \\begin{bmatrix}0&1\\\\-\\frac{1}{\\tau_s^2}&-\\frac{2 \\zeta}{\\tau_s}\\end{bmatrix} \\begin{bmatrix}x_1\\\\x_2\\end{bmatrix} + \\begin{bmatrix}0\\\\\\frac{K_p}{\\tau_{s}^2}\\end{bmatrix} u\\left(t-\\theta_p\\right)$\n",
    "\n",
    "$y = \\begin{bmatrix}1 & 0\\end{bmatrix} \\begin{bmatrix}x_1\\\\x_2\\end{bmatrix} + \\begin{bmatrix}0\\end{bmatrix} u$\n",
    "\n",
    "3. Second Order Differential Equation\n",
    "\n",
    "$\\tau_s^2 \\frac{d^2y}{dt^2} + 2 \\zeta \\tau_s \\frac{dy}{dt} + y = K_p \\, u\\left(t-\\theta_p \\right)$\n",
    "\n",
    "* $K_p    = 2.0$  # gain\n",
    "* $\\tau   = 1.0$  # time constant\n",
    "* $\\zeta  = 0.25$ # damping factor\n",
    "* $\\theta = 0.0$  # no time delay"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEICAYAAABS0fM3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABEkklEQVR4nO3deVxUVf/A8c9hVQF33EBFEfcFFfd0NNRc0zRLnxaXJy1/2VO2mG0uZbY+ZZZlPpXaapaVVmqZS5ppioq7KLIIrriBiMj2/f0xEyGCgjIMy/f9et0Xc+eeufd7GZjv3HPOPceICEoppUovJ0cHoJRSyrE0ESilVCmniUAppUo5TQRKKVXKaSJQSqlSThOBUkqVcnZLBMaYMsaYLcaYncaYvcaY6TmUMcaY2caYcGPMLmNMG3vFo5RSKmcudtz3ZeBWEUk0xrgCfxhjVojI5ixl+gIBtqUD8IHtZ66qVq0qfn5+dgpZKaVKpm3btp0WEe+cttktEYj1TrVE26qrbcl+99og4FNb2c3GmIrGmJoicjy3/fr5+RESEmKXmJVSqqQyxkTnts2ubQTGGGdjTChwClglIn9lK+IDxGRZj7U9p5RSqpDYNRGISLqIBAK+QHtjTPNsRUxOL8v+hDFmnDEmxBgTEhcXZ4dIlVKq9CqUXkMich5YB/TJtikWqJ1l3Rc4lsPr54lIkIgEeXvnWMWllFLqBtmtjcAY4w2kish5Y0xZoCfwWrZiy4AJxphFWBuJ46/VPqCUKjypqanExsaSnJzs6FBUPpQpUwZfX19cXV3z/Bp79hqqCSw0xjhjvfJYLCI/GWMeAhCRucByoB8QDiQBo+0Yj1IqH2JjY/Hy8sLPzw9jcqrFVUWNiHDmzBliY2OpV69enl9nz15Du4DWOTw/N8tjAR62VwxKqRuXnJysSaCYMcZQpUoV8tuWqncWK6VypUmg+LmR98yeVUOlwicLUzntFsKttzgRdKkSNGzo6JCUUipf9IrgJqxbBw/PfY2ILzpTpXUXaNsWHnsMUlMdHZpSxd6ZM2cIDAwkMDCQGjVq4OPjk7mekpJSoMc6cOAAgYGBtG7dmsOHD9/QPkaNGkW9evUyY5w9e3aBxmhPekVwgxIT4en7jhEe+wGxTf15699VeHfySrj3Xi5bejCi4WheeurfNGvm6EiVKp6qVKlCaGgoANOmTcPT05Mnn3wyc3taWhouLgXzEfbDDz8waNAgpk+/aki0XKWnp+Ps7HzFc2+88QZ33nlngcRUmPSK4AZNmgRjYqfzJf9ieFI4zzz6F1SqxF/TF/O67OetHx7g7od+JCnJ0ZEqVXKMGjWKxx9/nB49evD000+zZcsWOnfuTOvWrencuTNhYWEALFiwgCFDhtCnTx8CAgKYNGkSYP3wHjVqFM2bN6dFixa8/fbbLF++nFmzZvHRRx/Ro0cPAD7//HPat29PYGAgDz74IOnp6QB4enoyZcoUOnTowKZNm/IU88qVK2nTpg2tWrUiODgYsCa2MWPG0L17d+rXr+/4qwcRKVZL27ZtxdFWrRJp1GyqnKKKVOKMfPbZP9v27BFx6TVFFjSsJC+WHS///rfj4lTqZuzbt++K9alTRSBvy9ixV+9v7Ngry0ydmvdYpk6dKm+88YaMHDlS+vfvL2lpaSIiEh8fL6mpqSIismrVKhkyZIiIiMyfP1/q1asn58+fl0uXLkmdOnXkyJEjEhISIj179szc77lz567Y/9/nPWDAAElJSRERkfHjx8vChQtFRASQr7/+OscYR44cKX5+ftKqVStp1aqV7Nq1S06dOiW+vr4SEREhIiJnzpzJPF6nTp0kOTlZ4uLipHLlypnHKwjZ3ztb7CGSy+eqVg3lU0ICDH9pMR9mTOeNTt50rV6ee+75Z3uzZvDB3VOY8uC/2J7emeYfv8Di3jW56y7HxaxUSTJs2LDMKpn4+HhGjhzJoUOHMMaQmqV9Ljg4mAoVKgDQtGlToqOjadasGRERETzyyCP079+f3r17X7X/1atXs23bNtq1awfApUuXqFatGgDOzs4MHTo019iyVw39+OOPdOvWLbNPf+XKlTO39e/fH3d3d9zd3alWrRonT57E19f3Rn8tN0WrhvLpgw+gwc4M2h9x4xPn55k314XsvbX+PcaZrsMbMZ/RTOFFpr8oyFUjKCmlboSHh0fm4xdeeIEePXqwZ88efvzxxyvugnZ3d8987OzsTFpaGpUqVWLnzp10796dOXPm8MADD1y1fxFh5MiRhIaGEhoaSlhYGNOmTQOsd+1mbxe4FhHJtTtnTvE5iiaCfPr2G+H1+PeZmvQeL9/7CNWrX13GGHjvPZhVdRTDXOeT5tuf9esLP1alCtK0aXmtGIJ5865+/bx5V5axfbbelPj4eHx8rAMWL1iw4LrlT58+TUZGBkOHDuWll15i+/btV5UJDg7m22+/5dSpUwCcPXuW6OhcR3C+pk6dOvH7778TGRmZua+iSKuG8iE6GlK37aQ2MXzlPJq37s79xo2KFaHnsOq8eSCNmdErePPD41gsNQsvWKVKgUmTJjFy5Ejeeustbr311uuWP3r0KKNHjyYjIwOAV1555aoyTZs2ZcaMGfTu3ZuMjAxcXV2ZM2cOdevWzXd83t7ezJs3jyFDhpCRkUG1atVYtWpVvvdjb0aKWZ1FUFCQOGpimrfeTufcO72pfLoGv3T9nJUrrn0H35490G7Ex0TvnURXpy2sPeJPrVqFFKxSN2n//v00adLE0WGoG5DTe2eM2SYiQTmV16qhfFiwehO3p65h2W1rGHLH9cs3bw7tK/+br+Re/pX+aY6Xy0op5WiaCPLhuWEZNDzrzq6kCdxxR97G85gwARYwipEs5MN5yRTwDZFKKXXTNBHkw90pYVS4Ywinlz9HXufHGTwYjvtX40Llc1j6BJKeXryq4pRSJZ8mgvxYtgwGDryqu+i1uLrCupVV+LpzGv12HOTYpQj7xaeUUjdAE0Ee/bRjMWlr15DaKzjfr23coAy3TfmUe6PK4++WQ39TpZRyIE0EeRAVJfz4/kTWeyex4cKeG9pH13ZDcepmgW+/LeDolFLq5mgiyINHHk2j3YpqrPbzo1JCtxvf0ciRsGABJ09lFFxwSpVgL7/8Ms2aNaNly5YEBgby119/ATBr1iyS8jCiY17LZbV582Y6dOhAYGAgTZo0ybyruETLbRCioroU9qBz6ekiFbzS5Rg1pD7hcvDgje/r4Rm/yukyzlKr+wNy7FjBxaiUPeQ0cFlh+vPPP6Vjx46SnJwsIiJxcXFy9OhRERGpW7euxMXFXXcfeS2XVcOGDSU0NFRERNLS0mTv3r35jNzx8jvonF4RXEdYGPhe2Ecinlys7k+DBje+r01byrMiIJ3bZRk//FBgISpVIh0/fpyqVatmjslTtWpVatWqxezZszl27Bg9evTIHDZ6/PjxBAUF0axZM6ZOnQqQY7lff/2VTp060aZNG4YNG0ZiYuJVxz116hQ1a1pHAXB2dqZp06aAdejo++67j1tvvZWAgAD+97//AZCYmEhwcDBt2rShRYsWLF26NHNfn376KS1btqRVq1bcd999AMTFxTF06FDatWtHu3bt2Lhxoz1+ffmTW4YoqkthXxEsWCAypurj8qnznTJw4M3t653Z6XJnrUmywvSULCPhKlUkZf9WyTSEaVzx3IAvBwjTkGUHlmU+92HIh8I0ZOyyf8aiPppwVJiG1HyzZp6Pf+HCBWnVqpUEBATI+PHjZd26dZnbsn/T/3t457S0NLFYLLJz586rysXFxUnXrl0lMTFRREReffVVmT59+lXHnT59ulSsWFEGDx4sc+fOlUuXLomIdejoli1bSlJSksTFxYmvr68cPXpUUlNTJT4+PvMY/v7+kpGRIXv27JGGDRtmHv/vGEeMGCEbNmwQEZHo6Ghp3Lhxnn8neaVXBAVs3ZbTtK/7FltuXUpQ+5sbHXDoECdWHnuezvIXIWsSKKLjTylVJHh6erJt2zbmzZuHt7c3d999d64Dyy1evJg2bdrQunVr9u7dy759+64qs3nzZvbt20eXLl0IDAxk4cKFOQ4mN2XKFEJCQujduzdffvklffr0ydw2aNAgypYtS9WqVenRowdbtmxBRHj22Wdp2bIlPXv25OjRo5w8eZI1a9Zw5513UrVqVeCfIah/++03JkyYQGBgILfffjsJCQlcuHChAH5jN04HnbuOrXvO8p9ITz6p2oIZHW7u1+XjA807erFxcxd6Zaxg6dK7GD06HzclKOVAMvXqmyF/HPHjVc+NazuOcW3HXfFcLa9aOb7+epydnenevTvdu3enRYsWLFy4kFGjRl1RJjIykjfffJOtW7dSqVIlRo0adcVw1Jnxi9CrVy+++uqr6x7X39+f8ePHM3bsWLy9vTlz5gzAVUNKG2P44osviIuLY9u2bbi6uuLn50dycnKuQ1BnZGSwadMmypYtm4/fhH3pFcE1XLoERzb60uhsOjt/WYNtnoqbMmQILG3gxqBGD7Bwxd6b36FSJVRYWBiHDh3KXA8NDc0cAdTLyyvzW3RCQgIeHh5UqFCBkydPsmLFiszXZC3XsWNHNm7cSHh4OABJSUkcPHjwquP+/PPPiG0wzkOHDuHs7EzFihUBWLp0KcnJyZw5c4Z169bRrl074uPjqVatGq6urqxduzbzKiM4OJjFixdnJpG/h6Du3bs377333hXn5Wh6RXANoaHQIn0He2mGX6My2P4WbsqQIdD1t3LMXJ/Iv6v/wIULzfHyuvn9KlXSJCYm8sgjj3D+/HlcXFxo0KAB82wjN44bN46+fftSs2ZN1q5dS+vWrWnWrBn169enS5cumfvIXm7BggWMGDGCy5cvAzBjxgwaNmx4xXE/++wzJk6cSLly5XBxceGLL77InIymffv29O/fnyNHjvDCCy9Qq1Yt7rnnHgYOHEhQUBCBgYE0btwYgGbNmvHcc89hsVhwdnamdevWLFiwgNmzZ/Pwww/TsmVL0tLS6NatG3Pnzi2MX2mu7DYMtTGmNvApUAPIAOaJyDvZynQHlgKRtqe+E5EXr7XfwhyGetYsIfrpl2iQcpK/7pvDp58WzH4bdd3DZ1uG8UzKOzz4dW+dxlIVSToM9ZWmTZuGp6cnTz75pKNDua6iNAx1GvCEiDQBOgIPG2Oa5lBug4gE2pZrJoHCtnpHBO0bTmWLZQXt2xfcfkcEN2dpyr0M4ieWLy+4/Sql1I2wWyIQkeMist32+AKwH/Cx1/HsIblsJO1jDVs86tKhQ8Htt18/WMogBrGUlSuEDL3RWKkib9q0acXiauBGFEpjsTHGD2gN/JXD5k7GmJ3GmBXGmGaFEU9erZoRSP0MLz6avphWrQpuv0FBcKy+O+JxlvZ9h5PDPS1KKVVo7J4IjDGewBLgMRFJyLZ5O1BXRFoB7wI/5LKPccaYEGNMSFxcnF3jvcLWrZh27egS5I2bW8Ht1skJfvrRleVNEml5aillPVILbudKKZVPdk0ExhhXrEngCxH5Lvt2EUkQkUTb4+WAqzGmag7l5olIkIgEeed1RpibJCKwZQsF2jiQReemftQfPp7JiYG4OGnnLaWU49gtERjrnRQfA/tF5K1cytSwlcMY094Wzxl7xZQf249vZ82iV/nYZbfdjtFn3Gt47tiLuXjRbsdQSqnrsecVQRfgPuBWY0yobelnjHnIGPOQrcydwB5jzE5gNjBc7NWfNR9E4Jn3ttAyOpmNVdyxW0ReXtYrjrVr7XQApYq32NhYBg0aREBAAP7+/jz66KOk2Cb+XrduHRUqVKB169Y0atSIbt268dNPP2W+dtq0afj4+BAYGJi5nD9//or9R0VF0bx58+vGMXPmzAI9r9yEhoay3AFdCe3Za+gPETEi0jJL99DlIjJXRObayrwnIs1EpJWIdBSRP+0VT35ERsLBV24j5bI3374zw27HOX0aFnt6s+D5R/n3Y8fsdhyliiMRYciQIQwePJhDhw5x8OBBEhMTee655zLLdO3alR07dhAWFsbs2bOZMGECq1evztw+ceJEQkNDM5eKN3hXaEEmgrS03Mcsu5FEcK395ZUOMZGDnTuhOfvYndGaDv6N8zVHcX4cOgQvpR2j25FIloSusN+Vh1LF0Jo1ayhTpgyjR48GrOMOvf3223zyySc5TjYTGBjIlClTrhi+IT8WLFjAkCFD6NOnDwEBAUyaNAmAyZMnc+nSJQIDA7nnnnsA+Pzzz2nfvj2BgYE8+OCDpKenA/Dxxx/TsGFDunfvztixY5kwYQIAo0aN4vHHH6dHjx48/fTTbNmyhc6dO9O6dWs6d+5MWFgYKSkpTJkyha+//prAwEC+/vprzp49y+DBg2nZsiUdO3Zk165dgPVqZ9y4cfTu3Zv777//hs43K22lzMG+fdCUfeylGc3s2KG1fXuIOjYZt4u78Q5pxM6dEBhov+MpVZzs3buXtm3bXvFc+fLlqVOnTuZ4Qdm1adOGN954I3P97bff5vPPPwegUqVKrL1ONWxoaCg7duzA3d2dRo0a8cgjj/Dqq6/y3nvvZY4JtH//fr7++ms2btyIq6sr//d//8cXX3xBz549eemll9i+fTteXl7ceuuttMrS7/zgwYP89ttvODs7k5CQwPr163FxceG3337j2WefZcmSJbz44ouEhIRkJrNHHnmE1q1b88MPP7BmzRruv//+zDi2bdvGH3/8USCD12kiyMHWQ9EMrvs+fzh1pH1O90IXEGdnGNC4HytDh9I3dTsrVtyiiUAVXfa4NL7GZXBuo3fm9vzf27KaOHFivm4CCw4OpkKFCgA0bdqU6OhoateufUWZ1atXs23bNtrZRqG8dOkS1apVY8uWLVgslszhpocNG3bFoHbDhg3LHLMoPj6ekSNHcujQIYwxpKbm3IX8jz/+YMmSJQDceuutnDlzhvj4eABuv/32AhvBVKuGcrD75C6apUazt1EETe2YCAD69oWV9KEPK8kyaKJSRY9IwS/X0KxZM7KPK5aQkEBMTAz+/v45vmbHjh03NT7S37OhgbUqKqf6dxFh5MiRme0OYWFhTJs27aoklJ2Hh0fm4xdeeIEePXqwZ88efvzxxxyHzf77WNn9nQSz7u9maSLIJiMDjoe0pskJd/Yfmoi9x9y67Tb4zasxtzivYVvCRyRkv+VOqVIqODiYpKQkPrWN9piens4TTzzBqFGjKFeu3FXld+3axUsvvcTDDz9c4LG4urpmfmsPDg7m22+/5dSpU4B1eOno6Gjat2/P77//zrlz50hLS8v8Jp+T+Ph4fHysI+5knWwn67DZAN26deOLL74ArL2kqlatSvny5Qv69DQRZBcdDdVOp3E+oyruF0ZQpYp9j1e9OlTpfITdtS7TpdYbrFlj3+MpVVwYY/j+++/55ptvCAgIoGHDhpQpU+aKHjwbNmzI7D768MMPM3v2bIKDgzO3v/3221d0H42KirqhWMaNG0fLli255557aNq0KTNmzKB37960bNmSXr16cfz4cXx8fHj22Wfp0KEDPXv2pGnTppnVTNlNmjSJZ555hi5dumQ2NAP06NGDffv2ZTYWT5s2jZCQEFq2bMnkyZNZuHDhDcV/PXYbhtpe7D0M9c8/w5wBy3mMWbzS/ddC6eL/xOQkPJbeQoW4WoQP+4kPPrD/MZW6Hh2GOv8SExPx9PQkLS2NO+64gzFjxnDHHXcUehxFaRjqYmnvXqFp3Q/YV7EsjZsUTpIc2KccKw/M5bYzUfzyy3WrTpVSRdS0adMIDAykefPm1KtXj8GDBzs6pDzRXkPZhByMpW+ln/izkSetmhbOfMKdO0OYR1tqXDxBSmQs4eG+BAQUyqGVUgXozTffdHQIN0SvCLLp1csQdLI6Z9z60bp14RzTzQ263erEb2UD6dPoP+zapZcEqmgoblXH6sbeM70iyGbsXT7wwEW++3QuVCq84z71JHzjuo1eYecJuGU30LLwDq5UDsqUKcOZM2eoUqVKrv32VdEiIpw5c4YyZcrk63WaCLKLibEOBlepELMA0K2b4a8TIxnwwIccd3K//guUsjNfX19iY2Mp1DlA1E0rU6YMvr6++XqNJoJswv9YRp1GDXC9xt2L9vLUXbPgpdU0iDgPhTPtglK5cnV1pV69eo4OQxUCbSPIQkRY+MUkPri0geOJxx0TxG23wS+/OObYSqlSSRNBFpOev0jDaE8iqlTi1OGaDokh/hYLMV8sYuLMvQ45vlKq9NFEkMVfGzwJ2FufbcuXcfx44TeOXb4MtedsokLUfhasfYcTJwo9hEKRnAy//QYffQQLFsCXX0KWsbmUUoVM2wiy2L9PaMo+9tHU7oPN5cTdHRqUG8CG6nPoFZnOihVgG4q92LtwAebPh2/WHGL7+a/oGAF1YmrjTDrOrgn8ete7VO7blZl3z6OMS/56PCilbo4mApu4OHA/G00S5UjxqEy2kWcLzYhbOvPzslfoz1/8+HPxTwQisGgRPPVEBl2Pf83T5V+lW/Iu9pSvxCEGkoYL6STjtymKTkuicZ+2CQYOhPHjeXNtGepV8eWOOwxOeu2qlN1oIrDZtw8Ceg0hbO95ajfegpNTe4fE0b8/9JvUj2lM4z+/ZJCS4oSbm0NCuWn79sH9T+3Ebft6lpz4HIPw35SHuL/HJsqlDCO49kAAklNS6Tb8ETIan6HP5drw5Zekd+pE8/LxvNa4Ec+9uIk3X/ZiwAAHn5BSJZQmApv9+yEgPZrwapdpUqeaw+Jo0gSMnx9x0ZVpVP59/vhjArfe6rBwbsoT737HuBND6XPRmWec57Gmyiiee8GJ1weOp27drCVdgY7/rLZsycRyvUnYfDufbAhnZ+J9TBz4BvOHBPDOO5DPLtJKqevQC26bffugwerRhIc9TlBAHYfFYQz0GnCB5R0P0a/Bo3z/84Xrv6goOnOGH/fMolZ8NZp2HEflR4cRdsiJCRPIlgRy9tyDt1Kp5XHalt/BJjqxiU4ELXuMpv96jVmzU8nIsP8pKFVaaCKwOXQIGhBBeFJHmjR27K/ljn5e/Fy+Jf32evDz7ycdGkt+paSnkBa2Hzp1wqVzJ7qGHGXzW+/z7n+9yM98GtWrw9uvenFwZxPiRj9NM/YS4L2Uv/ZPZtE3Q+ndG2Jj7XceSpUmmghsDh+GBoQTTgMaNHBsLN27w7bfNxNwxpmL27w4fNix8eRVZHQaD87pS0KXdsgTT8Brr1GhogvNm9/4PqtWhU8+gcXrqvN4ow+Z0qY2P2zcRK/VT9OmeQqLFhVc/EqVVpoIgPR0iCj/Ef7O+zhcP4L69R0bT9my0L2nO7/Rkz6s5OefHRtPXqxdCx26ruGpl9bydkchdni/At2/xQKHVvbGv00ULWUfjTnAyviOvPbmVPZEHy3QYylV2mgiwDrOXLWqv5JYJh1T9wwFOCf0DevfH5bTj/6ui/nu16J9Z1l4ONw79BLfxbzIiuThnGj/G7UrFHz/W3d3ePUVJ5as9+axuj/wYdWurNr3Ip9PbMrFy4kFfjylSgvtNQTUrAmf3TaAxEOhPDuip6PDAayJ4LN79/LWt8vZMuZ94EVHh5Sj8+dhwMA03jw3hhhq81+Pz/lrtH2/X3TtCjt3GT5a9AyT47fz0seReAwdDh9/DNWrk5FBsb3vICUFDhywtlnFxEB0TBrHj7qQkAAJF4SWw5bSIziVwY0H4+rsCsDzbx1i49aLVHevh0+VCtSsaf2b9vGxNsz7+oKrq4NPTBVpdksExpjawKdADSADmCci72QrY4B3gH5AEjBKRLbbK6bcuLtDcIU06NqZp8c6uF7IpnZteOaF7hxe9zZNw8MdHU6ORGDoA4fpXr0LTQ6Xx2J28etSp0K5Ga98eXh8XA1E1iP/uQzTX4TAQFaOfZyHN7Zj8evdadvW/nEUhG++Eb74OYoDe8pyJLQyddIjaFTtG/ybzcTneE3aHuhMNU5RjiTKhf9JmZcEJ+/G4OIKrq4MOBVLO49TnD/VhZPnOnOc6uzwPsORtis4cqofx3dMp2ZtF+rXh2FjjjP+3ho6v4C6gj2vCNKAJ0RkuzHGC9hmjFklIvuylOkLBNiWDsAHtp+FLzwch7cSZ3Ob/22Y/3uJoMgYR4eSo3nz4HTc0ywKOUmX4NbMvb8MHQr53TPGYNzLwMyZnL+lOwH/6s/0Omn0GfE5D951D889Z21zKSqSk8Hd3TqDlAE4cIADHz1G/9O/Mi2sBo3SzxOLL2HxNQhPSuZY+RR20ZtTVCMRT6reNhO/Fmm8FvwKzsYFUlN5d/qbXEzdTKWzHalGFXw4SlD6BursDaVO3CGqyZscPeJDRGxtIhLWM2ORJ8//eyGmQQPw96f7xLWcOOZGiwqdaejnSd26UKeO9cuIr6816RZG3sjIsA5Fcv68dYmP/+dxQgI4O8OYsZcBcHexztmxYt1ZJr+9k7SL5XE62ZbLlyE1Fc7V+4h0ScVt94NkpDuRkQE1LT/x0LPhDGg4gAaVrf/r73x+kKnz1+F+oTFeZ7vh6gqu7qlc8lmFu1NZql3sgbu79ctivcAjPPJoGrW8amUOg7J+QwZr1hjc3QxubtbZBl1d/1lcXKyLszP4+XHVl5P9+6/s/SZy5ZKR8c/SsKH1PqOsVq+GqChrO2damvVn1iUjw/rzupOWiUihLMBSoFe25z4ERmRZDwNqXms/np6eMn/+fBERSUlJEYvFIp999pmIiFy8eFEsFossWrRIRETOnz8vFotFlixZIiIicXFxYrFYZNmyZSIicvz4cbFYLPLpkk9lf/cW8svUh8VisciqVatEROTw4cNisVhk3bp1IiJy4MABsVgssnHjRhER2b17t1gsFtmyZYuIiOzYsUMsFovs2LFDRES2bNkiFotFdu/eLSIiGzduFIvFIgcOHBARkXXr1onFYpHDhw+LiMiqVavEYrHIkSNHRERkxYoVYmnfXo5XrSqSlibLli0Ti8UicXFxIiKyZMkSsVgscv78eRERWbRokVgsFrl48aKIiHz22WdisVgkJSVFRETmz58vFotF/jZv3jwJDg7OXJ8zZ4706dMnc33WrFkycODAzPU33nhDhgwZIiIiYWEiHi7TpB+e8q+Gd8m4R8/Kiy++KPfcc09m+RdeeEFGjRqVuT558mQZO3Zs5voTTzwh//d//5e5/uijj8qjjz6auf5///d/8sQTT2Sujx07ViZPnpy5PmrUKHnhhRcy14N7jhB3v27ydnNvOUpN8aWzVKwwU779ViQjQ2TIkCHyxhtvZJYfOHCgzJo1K3O9T58+MmfOnH/2Fxws8+bNy1y3WCz5/tv75pslsnGjyIQJcVK2rEWadB8iE0dUlt29uojF1VVWVKsmmzp0lftuKSteHrXElZ/Fz0+kW48DUr1WZ3nggbXy6aciH3xwQFq1ssg331z9txcSIvLf/+6QRo0s8sgjO+Q//xHp2mullKvcUio2WCyuXBYfFkkl54YyqkdZWRjsLes6dxZLuXJy2N1djnk4y9veSGNTRZ7gERnNxxJY63YpW6mSePguEHe3DKladZl4lO8sc75fLlHnojL/9s6dOy8zZ4rce+8i8fe3yOuvX5TZs0XuueczqV/fIjNmpMjUqSJ9+86Xmj5dJepsjJy+eDrzb+/3LWekSvAn4lp3tMBt4kW81CZaajZoLpW9y8sAt8/kXj6VTgyWxuUryStdkN3DLPJKmzZyt6+vHGrVVtbWRcaVdZd+VJQDNJRI6sqjbsjdLshpKss5KshE3GSkMXLeHbnsVU6ecHeX/3N3l/Ne5SXWCxnt7CIjqSD7aCyhppkMLYcM9zSymh6ynD7Si7pyR6UysrAlcnz4QBnVuLG80KGDfNXNIlO6I808q0lnbpMJzJYHPaZKvequUjfAN8tH+hAp372RtP2wrRxLOJb5t9fpiTeEMZ0Fj9YCc6xlff8UapUXmgZneb1Fyg+qLc3mNJNLyZcy//bqPPyAMKGWUKaFwCJr2ebvC3WchDa32V4bJzh1FiBEcvlcLZQ2AmOMH9Aa+CvbJh8g69fdWNtzV0wGYIwZB4wDcHcv+Nm79hwNp/mB3Sxq6wZ4Fvj+b0rZsoi3N7/N/oIjboU7a1puUlPhvvugb9oyovEmyuVrQl4BR8/bXaO6Ew8168mSHav5ls00YwC1448z7f2NvL5wJuUK6d68y5dh1SrYfzCF+6auokFYbXqnL6EJ23hyQzIujdLYZvGAhDYwbRo+LXri9+ExGly8n6lT3Rg0CCIiXBkzxpV77zVYLBAWZh2zqVatq4/Xtq31G+eyZTBmDAQGwtatlXnqqUq8914TGjRw44cfajNzZk0enLmMwKC6/LXxL5g6lZQP/oflP9OpdH4z5+KcSceZbqwn4OIGvr94jkXnx1JDxvPlaU8+MBeoMrEfR2r6gGdTOHyY1OkvkrBzFudOuJB8uDYRk97HFTdO132NlPjjxL8xnnLxvpRjK+XcdrIquDZBFRrDGVeIiaHVPZ3ZGhfGkkTDbwjf4MM5KjHn2DH+kHQelI+Ix4clnCQhDS6WceJMHW8QbyhThoM9h/P6X2mEpwkXj/oyiP9yGXdO1xpEWkoSvx3fQLqU4RIvUsZnCzXfbsk9zf8F7y6CS5d4N2Awc9d+RUJILE7n/NnCJNxcT3O8yl244MqpxOdw5zIHmUWlipfZ2TyNHoGtIPIUeHmRmJGBcxq4ZmRQngQaEYZrxnG2X06F1ASis7xPl5yPs+14AmkZaZnPJThHQJ0/wSXLKJdlzkN6AnhdOSfKBY6yNy6GDPnnbspU19NQ7hg4ZflscEoHkwHOKXn+mzVi58mpjTGewO/AyyLyXbZtPwOviMgftvXVwCQR2Zbb/oKCgiQkJKTA4hMBD//tnDzSkU69X2Pz4ol4FqFcsPTHdHbMrEFFOc2m5kf5+qMcPgkK2bPTLrLr24d4Z+9G2rrsZt1WDwIDHR3VPzIyrENcT3k2jaEV7mTqiaV8X7s6r4Ztok63ekyeDH36FGx1x+nTsGIF/PQT/LI8nRaJf3JH82EMjj2JSazGD2n/Yhm3s6V6TT5ZnMSwroE4Gce3aIvAkSMQGWldIiKs6/vP7uJI8h7i93bE+Xg1qnGKaj7LaHv7B3St2JjhvrdBUhLn484yd+MruKa64LzlYZxJRzBkNPkBqRDD5f13czm+EZdxJ8l3J04dfqFT42CGdRwN5csTHp9O7/dmcz7Jn4Q/XqWspwsVK4JzwGo8PQVfOlLFy5MKFcCzfDqvveJ8xft27hyEhlqr/8qWJbMa5++qGWdna8cBJyfrYy+vK88/KQnOnrVWq6SmWpeUlH+Wy5f/Wby94ZZbrnz9unWwZm0GyZeF9FRnUlIg6XIKF+QE6WlOuCX7ZlbZNLbs5o47k2lVoxVuztYBxKa8fZgVG05QNrk+bik1MQbSXM6RVO4AbhmVKJ/SGGOs8Xe4fReDbneiqXfTzL+dV2afYsfuZDzwpoxLWZydAecUxDkZVyc33J3LWH8HzsLLM5y2iUhQTn8Hdk0ExhhX4CfgFxF5K4ftHwLrROQr23oY0F1Ecp0erKATwalT0KL6SfbSDP/ypzl/vnDqQ/NqyxYYMyWY5X/+TrPaqzmzw+LQQeiOHoWWYx5g56aPubeNhb591/H0046L51rOnYPJ08+yZO/TTNyWwIPnVrOcfrxZdjQ9343jpfv64+F2832Fk5OhRv9J9HD7igGrujEg/TdOUIPvfavxfYtkjhydwdB2FgYNgp49i1abRV4kJsLJk9YRetu04Yq/vwsXhBmvpHIx+TKu4pVZn53unIiLC5RxKUtZd+fMD+oHHuCK7tnp6RAdDRUqWNsitHeT/Rhjck0E9uw1ZICPgf05JQGbZcAEY8wirI3E8ddKAvaQ9Y5if/+ilQQA2rWD+MOLuXShM032leW336Bfwd6rlS8+PrBJTvBTvaoklHuHJ590XCzXU6kSfDirMhMP/I/XXoPGn53nwfQ5/JwxiLPPXWDlxg4MfeEr8PMjIiaFjz8yBLZwo3Jl6zdHT89/vjGePWv9xrxrt1Crxw/U9N/HE+7dMWvXUmbNGmL/WM+m2un8WN2Jl49twtSvzx13wPt3QMeO1m+jxZWnp3Xx9796m5eX4bWZbkD2byd5u6x2dsbhN3Aq+/Ya6gLcB+w2xoTannsWqAMgInOB5Vi7joZj7T5a6KPvh4dDg0orCU+sRX1/wdaXo8gwBoYPrsKSN4cyhO9YvLi9QxMBK1bQ8NBe6oZEMCTdq1h8wDVubJ0U58j0irz11nPcnVKN1h6zeDa6EnTqBC4uXKhVl9Tym/l+yUBi9z3OWSrj7HEU51texiOpLLU3jKQ2MdxKOO3WLqDh6TRSmjTFvVcfeOIJXr5tMj9ucGdEly4sHWAdVqOofalQKjd2byMoaAVdNTRlWhquq91IdxYudEziv68Wvev2rVvhofbbWMTdtPVZzemIug6pHjobd4TKQV2tN271LBo33t00EYiI4PXHn8f51CJ8ogLwPVGNSpwj3SmD9GoHSHJyJebYEGKoTQT1iem7mRbjK/FQ14nUrZiHoVSVKgIcUjVUXOyPPM+QY5X5ua4HFv+ilwQAgoLgeNtUnMMjadhmEKtWhdK/f+EdXwQ2Hd7HtlGtuKVREwKDg4vYddNNMAb8/Wk/8Su++2kev9S8SMrZGly4ABeSUkjzW0GVMtWo59qJGjWgdXMY03o87do5OnClCk6pTwRHD1WlQUQ9wiPeZczzjo4mZ8bAXbc2Z4GHE6MjYln0bTL9+xfevL6ffQZzZ8/n+wNpvDq3Fa1LYJ1H9+7QvbsXkLVbiRswyDEBKVWIrpkIjDFlgAFAV6AWcAnYA/wsInvtH579HQ4XAjjEIQJybAwrKu69y5Mhc/5kR1IfXjwiJL5LoXRzPX0anpyYzvKza3m68rM0O/eU/Q+qlCpUuXZkNsZMAzYCnbDeCPYhsBjr0BGvGmNWGWNaFkaQ9nLhAqTHnUIwJLpWLtJTILZtC+VqtyOEIIIvfM+XXxbOcSdNghFn3+MCXqz1mMH40RUL58BKqUJzrSuCrSIyLZdtbxljqmHrAVRcHT4M/sHDididSLXGa3B2DnZ0SLkyBsaPh08eG8NYpzlM/LwD48bZ9xJm3TpYtfcpdpiP6CKbmfO+KVI32ymlCkauVwQi8jOAMSbH+aVE5JSIFFz3HQc4cQL8iCKyShp1qldwdDjXNXIkbLE4E+i+gTYDH7LrsS5fhvum/sKc828yu8slGg+rxYABdj2kUspB8tJYPNcY4wYsAL4UkfN2jagQ9ekDPbeN5UR4FB9OauHocK6rYkXY9GNXFvd2pu/mGNIz0nF2sk9H/tdeg/abThLgVpF7az7N/re9rv8ipVSxdN1EICK3GGMCgDFAiDFmCzBfRFbZPbpC4HIkGt92gfg2KfjB7Oyhhld17pm1hgrD7gU73QKyZw/Mfime0LRnGZG6jJl3dsXHxz7HUko5Xp5GvRKRQ8DzwNOABZhtjDlgjBliz+AKRVSUdaDwYqRCh27WYSi/++76hfMpNRX+9VAsL6ZNZjn9SO3QlfHjC/wwSqki5LqJwBjT0hjzNrAfuBUYKCJNbI/ftnN8dvXTwZ84umsjP6bscXQo+ffCC1x45nkeeH4LaWnXL55XM169RLWqHRno8QnPV3yC+fOL9zg5Sqnry8sVwXvAdqCViDwstqkkReQY1quEYik1Fb5cs4UqcRfZZM44Opx8s/yQzoFLBzn317/44suM678gj3zq7OCTNcf5d48aPDbZ76oZkZRSJc91E4GIdBORz0TkUg7bPrNPWPYXEwPrn+nFOVORj1+4x9Hh5FunmsFMb1eVqZvP8+xzySQmFsx+x62dR4XB99M46FeefrJ4tJsopW7OtW4o+9EYM9A2p0D2bfWNMS8aY8bYNzz7iYoCv/OGqLTGBJQvfvfFPfeUByEhUaQk+tEhdiUzZ97c/kQEli6FDRuo8P67zJ7aSKuElColrnVFMBbr0BIHjDFbjTHLjTFrjDERWO8y3iYinxRKlHYQFQX1iCSSesWtrRiwjpf/2ssevMgUpjKdt97M4PDhG9tX4sV0Hpzdi6Qx9yELFhTO2BVKqSLjWjeUnRCRSSLiDwwDXgIeB5qLSC8RWVpYQdpDeGQq9Wp/SmTV5GKZCMA6b/DJdgNJcs1gXOuGPPhMeL73ERICjbt9zyMzV/NS1wzOtm16/RcppUqUvHYfjRKRTSISKiJJ9g6qMOyLjcGvym9EBa6lXj1HR3NjnJxg9ruG+ywBvLDnMHGnH2flyry//tQpGDI4g1nbF7ElrRupPZdRpVwV+wWslCqSSu0w1Mdi3KgX7cuXdGO4n6OjuXEdOkCbOnOYeCyNr9cdoPvWi3z/qwedOl37dRER0HfgZf59dCY+HOXh9HVs7KONw0qVRnm6IiiJjof5Uu+cM5GhLxbbqqG/zX65JqtOL+Mv6cSMxEfpfZuwaVPu5TdvhsBR/+Mht+rcaRYxlO9Y8JU7DRoUXsxKqaIjT4nAGFPWGNPI3sEUlpQUOBmbSi2OEUttatd2dEQ3p0YNWL0aplV5j878yX+rNeH7yBlXlTt3Dt57D27rkcjHsZNpcy6e7m3H8dZXNenb1wGBK6WKhLzcWTwQCAVW2tYDjTHL7ByXXcXEgE+l3zlhvPGu7eSQ+X8LWvPm8OM6L4Z3eRPPsmGMe/JlUnfuAODiRaHPyF3UavkXGx75mnXJXUk/0pUhNT/mh1lPMHy4g4NXSjlUXtoIpgHtgXUAIhJqjPGzX0j2FxUF9do/TuSBY1RquR7raBnFX/PmELq+HwtDP2b/Nz/TILgX1K9Phnd1Bp36iY/OlCGcDszgefbWv4O/PnPS6iClVJ4SQZqIxJsSNE+tszO0uOxKlGc5Arz9HB1OgXJygtFtxiCtR8OTZ+HgQb75ZBFhlzz4oM1AyqR+yd13GeaPhfLlHR2tUqooyEsi2GOM+RfgbBuO+j/An/YNy766d4fuq/uBy+3cP6W+o8OxC2MMVKkCnTrhk9CRxLBZPHa7KfYN40qpgpeXxuJHgGbAZeArIAF4zI4xFY7ISPDzowRd6OTqttsM//mPJgGlVM7yMjFNEvAc8JwxxhnwEJFku0dmRyKCiYyk2N5JppRSBSgvvYa+NMaUN8Z4AHuBMGPMU3l43SfGmFPGmBwH+zfGdDfGxBtjQm3LlPyHf2Pmh87nxJ7NvHXi+8I6pFJKFVl5qRpqKiIJwGBgOVAHuC8Pr1sA9LlOmQ0iEmhbXszDPm9aSgp8vjicShcz2BHvURiHVEqpIi0vicDVNhT1YGCpiKSSh9lyRWQ9cPbmwit4MTFw9NV7icmoy+rX/+PocJRSyuHykgg+BKIAD2C9MaYu1gbjgtDJGLPTGLPCGNMst0LGmHHGmBBjTEhcXNxNHTAqCvyIIUoa4F+j2k3tSymlSoK8zFA2W0R8RKSfWEUDPQrg2NuBuiLSCngX+OEaMcwTkSARCfL29r6pg0ZGQl2iicJP24qVUoo89BoyxrgDQwG/bOVvqk7f1u7w9+Plxpj3jTFVReT0zez3eg5HpeDX5E2iU6tS10+AUtB/VCmlriEvVUNLgUFAGnAxy3JTjDE1jO12ZWNMe1ssdp9Fft/RGOq6HCLK/yD1/DQJKKVUXu4s9hWR6/X+uYox5iugO1DVGBMLTAVcAURkLnAnMN4YkwZcAoaLyHUboW/WyQhv/MIbEZ14l95gpZRS5C0R/GmMaSEiu/OzYxEZcZ3t7wHv5WefBSH2cHnqXkokKvIBTQRKKUXeEsEtwChjTCTWYSYMICLS0q6R2cHlyxB3NAVv4jhhauHr6+iIlFLK8fKSCErMlCUxMeBbbz7HYivg7XcGN7fqjg5JKaUcLi/dR6OBisBA21LR9lyxExkJfo3eJbpWHFWb7HN0OEopVSTkZayhR4EvgGq25XNjzCP2DswefHygl1Mj4jzq0rNNQ0eHo5RSRUJeqob+DXQQkYsAxpjXgE1YbwIrVpo2haZtm0FQc4ZN9XF0OEopVSTk5T4CA6RnWU+nON+FFR2NdhdSSql/5OWKYD7wlzHme6wJYBDwsV2jspNzl85RJvwALvfdY72hQSmlVJ4ai98CRmMdSfQMMFpEZtk5Lrv4bv93nNy7hWfC5zo6FKWUKjLyUjX0N4N1+OliWS2UkgL/fT2Nmomw82gz7H8Ps1JKFQ956TU0BVgIVAKqAvONMc/bO7CCFhMDF7/uw6l0Xw4seKlUzFWslFJ5kZc2ghFA67/nKTbGvIp1COkZ9gysoFnnIYgiCj9tK1ZKqSzyUjUUBZTJsu4OHLZLNHYUHW2dhyCautSt6+holFKq6MjLFcFlYK8xZhXWNoJewB/GmNkAIlIs5nuMiErDr/PDRCWVp65fBvlrHlFKqZIrL4nge9vyt3X2CcW+9sXGMiApkU1VXejgp0lAKaX+dt1EICIL/35sjKkE1BaRXXaNyg5OH65D3d238JXXfdw1ydHRKKVU0ZGXXkPrjDHljTGVgZ1Yew29Zf/QClZ0lBN+6ceIPt9dG4uVUiqLvNSRVLDNLzwEmC8ibYGe9g2rYKWmwrGYdHyJ5Qh1qFPH0REppVTRkZdE4GKMqQncBfxk53js4uhRqN5qBmfdXfBqdhB3d0dHpJRSRUdeEsGLwC/AYRHZaoypDxyyb1gFKyoK6tZcSnTVJLzrnXR0OEopVaTkpbH4G+CbLOsRwFB7BlXQmjeHx2v3IuUoPHRHC0eHo5RSRcp1E4ExpiHwAVBdRJobY1oCt4tIsbmzuGpVGFrbC6r0oduYGo4ORymlipS8VA39D3gGSAWwdR0dbs+g7CIqSuchUEqpHOQlEZQTkS3ZnkuzRzD2En0+muO7N3G4fPr1CyulVCmTl0Rw2hjjj3V4CYwxdwLH7RpVAfs9ej1Jh/bx3umfHR2KUkoVOXlJBA8DHwKNjTFHgceAh+wZVEFKS4PxI6rjG+/E6k29SStW1zJKKWV/eZmhLEJEegLeQGOgO3DL9V5njPnEGHPKGLMnl+3GGDPbGBNujNlljGmTz9jz5NgxqLK/MXEZNTn5239wycvoSkopVYrkmghsw0o8Y4x5zxjTC0gCRgLhWG8uu54FQJ9rbO8LBNiWcVh7JhU4nYdAKaWu7Vrfjz8DzgGbgLHAJMANGCwiodfbsYisN8b4XaPIIOBTERFgszGmojGmpogUaPtDdDTUK7eVqEu1dR4CpZTKwbUSQX0RaQFgjPkIOA3UEZELBXRsHyAmy3qs7bmrEoExZhzWqwbq5HOgoIjIDPyCJhEpGfj4fQKUveGAlVKqJLpWG0Hq3w9EJB2ILMAkAJDTrME5TikvIvNEJEhEgry9vfN1kEOx5/CLK0dUmQoE+GkSUEqp7K6VCFoZYxJsywWg5d+PjTEJBXDsWKB2lnVf4FgB7PcKJyKqUG9/W6JWfaNtBEoplYNcq4ZExNnOx14GTDDGLAI6APEF3T4AEBlpbSyOxF8TgVJK5cBunSmNMV9h7Wpa1RgTC0wFXAFEZC6wHOiHtRdSEjC6oGNIS4OjUanU4hix+FKvXkEfQSmlij+7JQIRGXGd7YL1ZjW7iYmBml0f4MQ28GryO2XL9rLn4ZRSqlgq0bO4R0SAX9kdRFZJpVaNEn2qSil1w0r0fba33gpNt/+bxLRfefupdo4ORymliqQSnQiMgZoXzkDnIAK6lnd0OEopVSSV/PqSqCi0lVgppXJXohPBlqNbCN+2inUS5ehQlFKqyCrRieDbTVtwjz3BL3LI0aEopVSRVWITwblz8M74LlS74Mz/3huJ5Dh4hVJKqRLbWBwZCbVPe3GM2lRL7o3JaWQjpZRSJfeKICIC6hNBJPWoX9/R0SilVNFVYhPB4cOCv88XhHt6UK++1gsppVRuSmzV0L7oOFrW+ZRw9zL419d6IaWUyk2JvSKIPJJCgyhfwi92wt/f0dEopVTRVWKvCI4d8KXB8fKEH5+lbQRKKXUNJfKKIC0NjkRlUJ8IDuOvNxYrpdQ1lMhEEBsL1d13cI6KlK/hQblyjo5IKaWKrhKZCCIiwL/rGMLrnKBK23WODkcppYq0EpkILl+GJsmXCK9kaOBdx9HhKKVUkVYiG4v79oW+vw8htZwHQx/VBgKllLqWEnlFAEB4OK6NGlKhgt5DoJRS11KiEwENGjg6CqWUKvJKZCL4ZPvHJB3YzUfxax0dilJKFXklLhEkJsJP327lonMG+86nOTocpZQq8kpcIggLg1P/HUR4aguWvfwvR4ejlFJFXonrNbR/P/innCKcVrSoo11HlVLqekrcFcH+/dCAcA7jT5Mmjo5GKaWKPrsmAmNMH2NMmDEm3BgzOYft3Y0x8caYUNsy5WaPufXwYRr4fEF4nWM0bXqze1NKqZLPblVDxhhnYA7QC4gFthpjlonIvmxFN4jIgII67t6z22kgkYQ32q9XBEoplQf2vCJoD4SLSISIpACLgEF2PB6pqXByexsCTpYlPGIcjRvb82hKKVUy2DMR+AAxWdZjbc9l18kYs9MYs8IY0+xmDhgeDhXOVIR0NzxS7sXD42b2ppRSpYM9ew3lNLZD9smDtwN1RSTRGNMP+AEIuGpHxowDxgHUuUZPoCsaipvq0BJKKZUX9rwiiAVqZ1n3BY5lLSAiCSKSaHu8HHA1xlTNviMRmSciQSIS5O3tnesB9+5Lp5nfHPZ5VNH2AaWUyiN7JoKtQIAxpp4xxg0YDizLWsAYU8MYY2yP29viOXOjBwwJj6Z5jc/Y3WazJgKllMojuyUCEUkDJgC/APuBxSKy1xjzkDHmIVuxO4E9xpidwGxguIhkrz7Ks65dof3xWhx1707Lljd7BkopVTqYm/jcdYigoCAJCQnJvUDNmrBlC9SunXsZpZQqZYwx20QkKKdtJevO4tOnISkJfH0dHYlSShUbJSoRRG9cTnqzpmC0x5BSSuVViUkEIsJ7Hz/I/1I2E3cxztHhKKVUsVFiEsGzL16gcUxZDlXzIi76qh6oSimlclFiEsEfv5WncWgTQlYsJSZGq4aUUiqvSkwi2L9PaM4e9tJC7yFQSql8KBGJIC4Oyp6NIYlyJHtU1Z6jSimVDyUiEezfDy26D2G37wXqdNiunYaUUiofSkQi2LcPmqcfZI/vRRrVqeTocJRSqlgpEXMW79wJnf4YwNryDWk7oK6jw1FKqWKlRFwRbNgALWQ/e+L707lTiTglpZQqNMX+U/PMGTiwN41GhHHQuSkdOjg6IqWUKl6KfSL4809o0HsgxzxcqN85QWclU0qpfCr2bQRrN16kfdlf2F4Xurb3dHQ4SilV7BT7RDBzehliB/RjfUNv7hnm5ehwlFKq2Cn2iaCMmxMNDuykwerV0NDR0SilVPFT7NsIOHTI+jPgqjnvlVJK5UGxTgSp6aksfOt+wlr5UrzmWVNKqaKjWCeCj37egcf6v/i4chRGx5VQSqkbUmwTQVoaPPVAbXoc9GDFkcc4d87RESmlVPFUbBuLd+6EgLgTxOHL+cNPU7GioyNSSqniqdheEfzxBwSzmtUEc8stOk2xUkrdqGKbCFZujia4wsesKe9P166OjkYppYqvYpkIjhyB9ac+psulA6ztsV4TgVJK3YRimQgefBC6bfHjkKcX5RMm0KyZoyNSSqniq9g1Fp85A4e3nSOU6Uzgc756uSdOxTKdKaVU0WDXj1BjTB9jTJgxJtwYMzmH7cYYM9u2fZcxps319hlzNJX/MZalDML/0dvp1Mk+sSulVGlhtysCY4wzMAfoBcQCW40xy0RkX5ZifYEA29IB+MD2M1eVyu4moGwMz1WKZtsM+8SulFKliT2vCNoD4SISISIpwCJgULYyg4BPxWozUNEYU/NaO/VJFIZ3as+cj7zw1FGnlVLqptkzEfgAMVnWY23P5bcMxphxxpgQY0zIuYpVGd7zU3r2LPB4lVKqVLJnIsjpFq/sY8PlpQwiMk9EgkQkyLteXaY+Wa1AAlRKKWXfRBAL1M6y7gscu4EySiml7MieiWArEGCMqWeMcQOGA8uylVkG3G/rPdQRiBeR43aMSSmlVDZ26zUkImnGmAnAL4Az8ImI7DXGPGTbPhdYDvQDwoEkYLS94lFKKZUzu95QJiLLsX7YZ31ubpbHAjxszxiUUkpdm96Tq5RSpZwmAqWUKuU0ESilVCmniUAppUo5Y22vLT6MMReAMEfHUYiqAqcdHUQh0vMtuUrTuULRO9+6IuKd04ZiNww1ECYiQY4OorAYY0L0fEuu0nS+pelcoXidr1YNKaVUKaeJQCmlSrnimAjmOTqAQqbnW7KVpvMtTecKxeh8i11jsVJKqYJVHK8IlFJKFaBilQiuNwdySWOMiTLG7DbGhBpjQhwdT0EzxnxijDlljNmT5bnKxphVxphDtp+VHBljQcnlXKcZY47a3t9QY0w/R8ZYkIwxtY0xa40x+40xe40xj9qeL6nvb27nWyze42JTNWSbA/kgWeZABkZkmwO5RDHGRAFBIlKU+iIXGGNMNyAR63SlzW3PvQ6cFZFXbcm+kog87cg4C0Iu5zoNSBSRNx0Zmz3YppytKSLbjTFewDZgMDCKkvn+5na+d1EM3uPidEWQlzmQVTEiIuuBs9meHgQstD1eiPWfqdjL5VxLLBE5LiLbbY8vAPuxTkNbUt/f3M63WChOiSBP8xuXMAL8aozZZowZ5+hgCkn1vycnsv0s6fOSTjDG7LJVHZWIapLsjDF+QGvgL0rB+5vtfKEYvMfFKRHkaX7jEqaLiLQB+gIP26oXVMnxAeAPBALHgf86NBo7MMZ4AkuAx0QkwdHx2FsO51ss3uPilAhK3fzGInLM9vMU8D3W6rGS7qStvvXvetdTDo7HbkTkpIiki0gG8D9K2PtrjHHF+qH4hYh8Z3u6xL6/OZ1vcXmPi1MiyMscyCWGMcbD1uiEMcYD6A3sufarSoRlwEjb45HAUgfGYld/fyDa3EEJen+NMQb4GNgvIm9l2VQi39/czre4vMfFptcQgK3r1Sz+mQP5ZcdGZD/GmPpYrwLAOjjglyXtfI0xXwHdsY7SeBKYCvwALAbqAEeAYSJS7BtZcznX7lirDASIAh78u/68uDPG3AJsAHYDGbann8Vab14S39/czncExeA9LlaJQCmlVMErTlVDSiml7EATgVJKlXKaCJRSqpTTRKCUUqWcJgKllCrlNBEolQtjTJUso0aeyDKKZKIx5n1Hx6dUQdHuo0rlQUkeKVQpvSJQKp+MMd2NMT/ZHk8zxiw0xvxqmz9iiDHmdds8Eittww5gjGlrjPndNoDgL9nuOFXKoTQRKHXz/IH+WIdY/hxYKyItgEtAf1syeBe4U0TaAp8AJeoucVW8uTg6AKVKgBUikmqM2Y11+JOVtud3A35AI6A5sMo6JA3OWEeiVKpI0ESg1M27DCAiGcaYVPmn4S0D6/+YAfaKSCdHBajUtWjVkFL2FwZ4G2M6gXW4YmNMMwfHpFQmTQRK2ZltatU7gdeMMTuBUKCzQ4NSKgvtPqqUUqWcXhEopVQpp4lAKaVKOU0ESilVymkiUEqpUk4TgVJKlXKaCJRSqpTTRKCUUqWcJgKllCrl/h8eG0p6G8vhYQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy import signal\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.integrate import odeint\n",
    "\n",
    "# tau * dy2/dt2 + 2*zeta*tau*dy/dt + y = Kp*u\n",
    "Kp = 2.0    # gain\n",
    "tau = 1.0   # time constant\n",
    "zeta = 0.25 # damping factor\n",
    "theta = 0.0 # no time delay\n",
    "du = 1.0    # change in u\n",
    "\n",
    "# (1) Transfer Function\n",
    "num = [Kp]\n",
    "den = [tau**2,2*zeta*tau,1]\n",
    "sys1 = signal.TransferFunction(num,den)\n",
    "t1,y1 = signal.step(sys1)\n",
    "\n",
    "# (2) State Space\n",
    "A = [[0.0,1.0],[-1.0/tau**2,-2.0*zeta/tau]]\n",
    "B = [[0.0],[Kp/tau**2]]\n",
    "C = [1.0,0.0]\n",
    "D = 0.0\n",
    "sys2 = signal.StateSpace(A,B,C,D)\n",
    "t2,y2 = signal.step(sys2)\n",
    "\n",
    "# (3) ODE Integrator\n",
    "def model3(x,t):\n",
    "    y = x[0]\n",
    "    dydt = x[1]\n",
    "    dy2dt2 = (-2.0*zeta*tau*dydt - y + Kp*du)/tau**2\n",
    "    return [dydt,dy2dt2]\n",
    "t3 = np.linspace(0,25,100)\n",
    "x3 = odeint(model3,[0,0],t3)\n",
    "y3 = x3[:,0]\n",
    "\n",
    "plt.figure(1)\n",
    "plt.plot(t1,y1*du,'b--',linewidth=3,label='Transfer Fcn')\n",
    "plt.plot(t2,y2*du,'g:',linewidth=2,label='State Space')\n",
    "plt.plot(t3,y3,'r-',linewidth=1,label='ODE Integrator')\n",
    "y_ss = Kp * du\n",
    "plt.plot([0,max(t1)],[y_ss,y_ss],'k:')\n",
    "plt.xlim([0,max(t1)])\n",
    "plt.xlabel('Time')\n",
    "plt.ylabel('Response (y)')\n",
    "plt.legend(loc='best')\n",
    "plt.savefig('2nd_order.png')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Higher Order Simulation\n",
    "\n",
    "A simple higher order simulation is the combination of $n$ first order equations. The value of the time constant is $10/n$ in this example. The first equation is a first order differential expression.\n",
    "\n",
    "$\\tau \\frac{dy_1}{dt} = -y_1 + 1$\n",
    "\n",
    "Additional equations are also first order differential expressions for $i = 2\\ldots n$.\n",
    "\n",
    "$\\tau \\frac{dy_i}{dt} = -y_i + y_{i-1}$\n",
    "\n",
    "See [Dynamic Simulation in Python](https://apmonitor.com/pdc/index.php/Main/ModelSimulation) for solution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Distillation Exercise\n",
    "\n",
    "Over 40,000 distillation columns in the U.S. consume 40-60% of the total energy in the petrochemical industry and 6% of the total U.S. energy usage. Reducing distillation energy usage would have a major impact on energy independence, greenhouse gas emissions, and give a competitive advantage to any manufacturer that is better able to implement improvements.\n",
    "\n",
    "The objective of this exercise is to maintain an overhead mole fraction for the light component in a binary distillation column that separates cyclohexane and n-heptane. The nominal steady state for the column is a reflux ratio at 3.0 with an allowable range from 1.0 to 10.0. This reflux ratio produces a cyclohexane product with a purity of 0.935 (mole fraction), but the target composition is 0.970 with an allowable range of 0.960 to 0.980. The objective is to design a PID controller to meet this target value and stay within the specified product range.\n",
    "\n",
    "See [Distillation Control](https://apmonitor.com/pdc/index.php/Main/DistillationControl) for full problem statement and solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.integrate import odeint\n",
    "from IPython import display\n",
    "\n",
    "# define model\n",
    "def distill(x,t,rr,Feed,x_Feed):\n",
    "    # Inputs (3):\n",
    "    # Reflux ratio is the Manipulated variable\n",
    "    # Reflux Ratio (L/D)\n",
    "    #rr = p(1)\n",
    "\n",
    "    # Disturbance variables (DV)\n",
    "    # Feed Flowrate (mol/min)\n",
    "    #Feed = p(2)\n",
    "\n",
    "    # Mole Fraction of Feed\n",
    "    #x_Feed = p(3)\n",
    "\n",
    "    # States (32):\n",
    "    # x(0) - Reflux Drum Liquid Mole Fraction of Component A\n",
    "    # x(1) - Tray 1 - Liquid Mole Fraction of Component A\n",
    "    # .\n",
    "    # .\n",
    "    # .\n",
    "    # x(16) - Tray 16 - Liquid Mole Fraction of Component A (Feed)\n",
    "    # .\n",
    "    # .\n",
    "    # .\n",
    "    # x(30) - Tray 30 - Liquid Mole Fraction of Component A\n",
    "    # x(31) - Reboiler Liquid Mole Fraction of Component A\n",
    "\n",
    "    # Parameters\n",
    "    # Distillate Flowrate (mol/min)\n",
    "    D=0.5*Feed\n",
    "    # Flowrate of the Liquid in the Rectification Section (mol/min)\n",
    "    L=rr*D\n",
    "    # Vapor Flowrate in the Column (mol/min)\n",
    "    V=L+D\n",
    "    # Flowrate of the Liquid in the Stripping Section (mol/min)\n",
    "    FL=Feed+L\n",
    "    # Relative Volatility = (yA/xA)/(yB/xB) = KA/KB = alpha(A,B)\n",
    "    vol=1.6\n",
    "    # Total Molar Holdup in the Condenser\n",
    "    atray=0.25\n",
    "    # Total Molar Holdup on each Tray\n",
    "    acond=0.5\n",
    "    # Total Molar Holdup in the Reboiler\n",
    "    areb=1.0\n",
    "    # Vapor Mole Fractions of Component A\n",
    "    # From the equilibrium assumption and mole balances\n",
    "    # 1) vol = (yA/xA) / (yB/xB)\n",
    "    # 2) xA + xB = 1\n",
    "    # 3) yA + yB = 1\n",
    "    y = np.empty(len(x))\n",
    "    for i in range(32):\n",
    "        y[i] = x[i] * vol/(1.0+(vol-1.0)*x[i])\n",
    "\n",
    "    # Compute xdot\n",
    "    xdot = np.empty(len(x))\n",
    "    xdot[0] = 1/acond*V*(y[1]-x[0])\n",
    "    for i in range(1,16):\n",
    "        xdot[i] = 1.0/atray*(L*(x[i-1]-x[i])-V*(y[i]-y[i+1]))\n",
    "    xdot[16] = 1/atray*(Feed*x_Feed+L*x[15]-FL*x[16]-V*(y[16]-y[17]))\n",
    "    for i in range(17,31):\n",
    "        xdot[i] = 1.0/atray*(FL*(x[i-1]-x[i])-V*(y[i]-y[i+1]))\n",
    "    xdot[31] = 1/areb*(FL*x[30]-(Feed-D)*x[31]-V*y[31])\n",
    "    return xdot\n",
    "\n",
    "# Steady State Initial Conditions for the 32 states\n",
    "x_ss =np.array([0.935,0.900,0.862,0.821,0.779,0.738,\\\n",
    "0.698,0.661,0.628,0.599,0.574,0.553,0.535,0.521,    \\\n",
    "0.510,0.501,0.494,0.485,0.474,0.459,0.441,0.419,    \\\n",
    "0.392,0.360,0.324,0.284,0.243,0.201,0.161,0.125,    \\\n",
    "0.092,0.064])\n",
    "x0 = x_ss\n",
    "\n",
    "# Steady State Initial Condition\n",
    "rr_ss = 3.0\n",
    "\n",
    "# Time Interval (min)\n",
    "t = np.linspace(0,50,500)\n",
    "\n",
    "# Store results for plotting\n",
    "xd = np.ones(len(t)) * x_ss[0]\n",
    "rr = np.ones(len(t)) * rr_ss\n",
    "ff = np.ones(len(t))\n",
    "xf = np.ones(len(t)) * 0.5\n",
    "sp = np.ones(len(t)) * 0.935\n",
    "\n",
    "# Step in reflux ratio\n",
    "rr[10:] = 1.0\n",
    "rr[180:] = 5.0\n",
    "rr[360:] = 3.0\n",
    "\n",
    "plt.figure(figsize=(10,7))\n",
    "plt.ion()\n",
    "plt.show()\n",
    "\n",
    "# Simulate\n",
    "for i in range(len(t)-1):\n",
    "    ts = [t[i],t[i+1]]\n",
    "    y = odeint(distill,x0,ts,args=(rr[i],ff[i],xf[i]))\n",
    "    xd[i+1] = y[-1][0]\n",
    "    x0 = y[-1]\n",
    "\n",
    "    # Plot the results\n",
    "    display.clear_output(wait=True)\n",
    "    plt.clf()\n",
    "    #plt.figure(figsize=(10,7))\n",
    "    plt.subplot(2,1,1)\n",
    "    plt.plot(t[0:i+1],rr[0:i+1],'b--',linewidth=3)\n",
    "    plt.ylabel(r'$RR$')\n",
    "    plt.legend(['Reflux ratio'],loc='best')\n",
    "\n",
    "    plt.subplot(2,1,2)\n",
    "    plt.plot(t[0:i+1],sp[0:i+1],'k.-',linewidth=1)\n",
    "    plt.plot(t[0:i+1],xd[0:i+1],'r-',linewidth=3)\n",
    "    plt.ylabel(r'$x_d\\;(mol/L)$')\n",
    "    plt.legend(['Starting composition','Distillate composition'],loc='best')\n",
    "    plt.xlabel('Time (hr)')\n",
    "\n",
    "    plt.draw()\n",
    "    plt.show()\n",
    "\n",
    "# Construct results and save data file\n",
    "# Column 1 = time\n",
    "# Column 2 = reflux ratio\n",
    "# Column 3 = distillate composition\n",
    "data = np.vstack((t,rr,xd)) # vertical stack\n",
    "data = data.T             # transpose data\n",
    "np.savetxt('data.txt',data,delimiter=',')\n",
    "\n",
    "plt.savefig('distillation.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### TCLab Exercise\n",
    "\n",
    "**Objective**: Simulate a third-order transfer function, state space, and fifth-order time series model of $T_{C1}$ temperature response to heater $Q_1$.\n",
    "\n",
    "Higher-order models are better at predicting system dynamics for complex systems. This exercise is to create a third-order transfer function, state space, and a fifth-order time series model as linear representations of step response data. If heater 1 ($Q_1$) controls temperature 2 ($Q_1$), the generated heat first raises the heater 1 temperature ($T_{H1}$). The heat is transferred by convection to heater 2 ($T_{H2}$) and finally reaches temperature 2 sensor ($T_{C2}$) by conduction. This heat transfer pathway suggests a third-order dynamic model, although a higher-order model may improve the fit even further.\n",
    "\n",
    "See [TCLab Higher Order](https://apmonitor.com/pdc/index.php/Main/TCLabHigherOrder) for problem statement and solution details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "from gekko import GEKKO\n",
    "import tclab\n",
    "import time\n",
    "from scipy import signal\n",
    "\n",
    "# Import data\n",
    "try:\n",
    "    # try to read local data file first\n",
    "    filename = 'data.csv'\n",
    "    data = pd.read_csv(filename)\n",
    "except:\n",
    "    filename = 'http:https://apmonitor.com/pdc/uploads/Main/tclab_data3.txt'\n",
    "# read either local file or web file\n",
    "data = pd.read_csv(filename)\n",
    "\n",
    "U = 4.7052403301\n",
    "Us = 15.45761703\n",
    "alpha1 = 0.012321367852\n",
    "alpha2 = 0.005\n",
    "tau = 20.298826743\n",
    "Ta = 293.15\n",
    "mass = 0.004\n",
    "Cp = 500.0\n",
    "A = 0.001\n",
    "As = 0.0002\n",
    "eps = 0.9\n",
    "sigma = 5.67e-08\n",
    "\n",
    "Am = np.zeros((3,3))\n",
    "Bm = np.zeros((3,1))\n",
    "Cm = np.zeros((1,3))\n",
    "Dm = np.zeros((1,1))\n",
    "\n",
    "T0 = Ta\n",
    "c1 = U*A\n",
    "c2 = 4*eps*sigma*A*T0**3\n",
    "c3 = Us*As\n",
    "c4 = 4*eps*sigma*As*T0**3\n",
    "c5 = mass*Cp\n",
    "c6 = 1/tau\n",
    "\n",
    "Am[0,0] = -(c1+c2+c3+c4)/c5\n",
    "Am[0,1] = (c3+c4)/c5\n",
    "\n",
    "Am[1,0] = (c3+c4)/c5\n",
    "Am[1,1] = -(c1+c2+c3+c4)/c5\n",
    "\n",
    "Am[2,1] = c6\n",
    "Am[2,2] = -c6\n",
    "\n",
    "Bm[0,0] = alpha1/c5\n",
    "\n",
    "Cm[0,2] = 1\n",
    "\n",
    "# state space simulation\n",
    "ss = GEKKO(remote=False)\n",
    "x,y,u = ss.state_space(Am,Bm,Cm,D=None)\n",
    "u[0].value = data['Q1'].values\n",
    "ss.time = data['Time'].values\n",
    "ss.options.IMODE = 7\n",
    "ss.solve(disp=False)\n",
    "\n",
    "# state space simulation with scipy\n",
    "sys = signal.StateSpace(Am,Bm,Cm,Dm)\n",
    "tsys = data['Time'].values\n",
    "Qsys = data['Q1'].values.T\n",
    "tsys,ysys,xsys = signal.lsim(sys,Qsys,tsys)\n",
    "\n",
    "# print linearized models\n",
    "print('State Space Model')\n",
    "print(sys)\n",
    "\n",
    "print('Transfer Function Model')\n",
    "tf=sys.to_tf()\n",
    "print(tf)\n",
    "\n",
    "# Time series model\n",
    "ts = GEKKO(remote=False)\n",
    "tt = data['Time'].values\n",
    "tu = data['Q1'].values\n",
    "ty = data['T2'].values\n",
    "na = 5 # output coefficients\n",
    "nb = 5 # input coefficients\n",
    "print('Identify time series model')\n",
    "# diaglevel = 1 option to see solver output\n",
    "yp,p,K = ts.sysid(tt,tu,ty,na,nb,pred='model',objf=1000,diaglevel=0)\n",
    "\n",
    "# Create plot\n",
    "plt.figure(figsize=(10,7))\n",
    "ax=plt.subplot(2,1,1)\n",
    "ax.grid()\n",
    "plt.plot(data['Time'],data['T2'],'b.',label=r'$T_2$ measured')\n",
    "plt.plot(ss.time,np.array(y[0].value)+data['T2'][0],color='purple',linestyle=':',\\\n",
    "         linewidth=2,label=r'$T_2$ state space')\n",
    "plt.plot(tt,yp,'g--',\\\n",
    "         linewidth=2,label=r'$T_2$ time series')\n",
    "plt.ylabel(r'T ($^oC$)')\n",
    "plt.legend(loc=2)\n",
    "ax=plt.subplot(2,1,2)\n",
    "ax.grid()\n",
    "plt.plot(data['Time'],data['Q1'],'r-',\\\n",
    "         linewidth=3,label=r'$Q_1$')\n",
    "plt.ylabel('Heater (%)')\n",
    "plt.xlabel('Time (sec)')\n",
    "plt.legend(loc='best')\n",
    "plt.savefig('tclab_higher_order.png')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}