{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#  14章 離散時間システムの制御 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from control.matlab import *\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "#plt.rcParams['font.family'] ='sans-serif' #使用するフォント\n",
    "plt.rcParams['font.family'] = 'Times New Roman' # font familyの設定\n",
    "plt.rcParams['mathtext.fontset'] = 'cm' # math fontの設定\n",
    "plt.rcParams['xtick.direction'] = 'in' #x軸の目盛線が内向き('in')か外向き('out')か双方向か('inout')\n",
    "plt.rcParams['ytick.direction'] = 'in' #y軸の目盛線が内向き('in')か外向き('out')か双方向か('inout')\n",
    "plt.rcParams['xtick.major.width'] = 1.0 #x軸主目盛り線の線幅\n",
    "plt.rcParams['ytick.major.width'] = 1.0 #y軸主目盛り線の線幅\n",
    "plt.rcParams['font.size'] = 11 #フォントの大きさ\n",
    "plt.rcParams['axes.linewidth'] = 0.5 # 軸の線幅edge linewidth。囲みの太さ\n",
    "plt.rcParams['mathtext.default'] = 'it'#'regular'\n",
    "plt.rcParams['axes.xmargin'] = '0'\n",
    "plt.rcParams['axes.ymargin'] = '0.05'\n",
    "plt.rcParams['savefig.facecolor'] = 'None'\n",
    "plt.rcParams['savefig.edgecolor'] = 'None'\n",
    "\n",
    "plt.rcParams[\"legend.fancybox\"] = True     # 丸角\n",
    "# plt.rcParams[\"legend.framealpha\"] = 1    # 透明度の指定、0で塗りつぶしなし\n",
    "# plt.rcParams[\"legend.edgecolor\"] = 'gray' # edgeの色を変更\n",
    "plt.rcParams[\"legend.handlelength\"] = 1.8  # 凡例の線の長さを調節\n",
    "plt.rcParams[\"legend.labelspacing\"] = 0.4  # 垂直方向（縦）の距離の各凡例の距離\n",
    "plt.rcParams[\"legend.handletextpad\"] = 0.7 # 凡例の線と文字の距離の長さ\n",
    "plt.rcParams[\"legend.markerscale\"] = 1.0   # 点がある場合のmarker scale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def linestyle_generator():\n",
    "    linestyle = ['-', '--', '-.', ':']\n",
    "    lineID = 0\n",
    "    while True:\n",
    "        yield linestyle[lineID]\n",
    "        lineID = (lineID + 1) % len(linestyle)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_set(fig_ax, *args):\n",
    "    fig_ax.set_xlabel(args[0])\n",
    "    fig_ax.set_ylabel(args[1])\n",
    "    fig_ax.grid(ls=':', lw=0.5)\n",
    "    if len(args)==3:\n",
    "        fig_ax.legend(loc=args[2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bodeplot_set(fig_ax, *args):\n",
    "    fig_ax[0].grid(which=\"both\", ls=':', lw=0.5)\n",
    "    fig_ax[0].set_ylabel('Gain [dB]')\n",
    "\n",
    "    fig_ax[1].grid(which=\"both\", ls=':', lw=0.5)\n",
    "    fig_ax[1].set_xlabel('$\\omega$ [rad/s]')\n",
    "    fig_ax[1].set_ylabel('Phase [deg]')\n",
    "    \n",
    "    if len(args) > 0:\n",
    "        fig_ax[1].legend(loc=args[0])\n",
    "    if len(args) > 1:\n",
    "        fig_ax[0].legend(loc=args[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 離散化"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 例題14.3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\frac{0.3333 z + 0.3333}{z - 0.3333}\\quad dt = 1$$"
      ],
      "text/plain": [
       "TransferFunction(array([0.33333333, 0.33333333]), array([ 1.        , -0.33333333]), 1)"
      ]
     },
     "execution_count": 114,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P = tf(1, [1,1])\n",
    "Pd = c2d(P, 1, method='tustin')\n",
    "Pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\[\n",
       "\\left(\n",
       "\\begin{array}{rll|rll}\n",
       "-1\\phantom{.}&\\hspace{-1em}&\\hspace{-1em}\\phantom{\\cdot}&1\\phantom{.}&\\hspace{-1em}&\\hspace{-1em}\\phantom{\\cdot}\\\\\n",
       "\\hline\n",
       "1\\phantom{.}&\\hspace{-1em}&\\hspace{-1em}\\phantom{\\cdot}&0\\phantom{.}&\\hspace{-1em}&\\hspace{-1em}\\phantom{\\cdot}\\\\\n",
       "\\end{array}\\right)\n",
       "\\]"
      ],
      "text/plain": [
       "StateSpace(array([[-1.]]), array([[1.]]), array([[1.]]), array([[0.]]))"
      ]
     },
     "execution_count": 116,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Pss = ss(P)\n",
    "Pss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\[\n",
       "\\left(\n",
       "\\begin{array}{rll|rll}\n",
       "0.&\\hspace{-1em}333&\\hspace{-1em}\\phantom{\\cdot}&0.&\\hspace{-1em}667&\\hspace{-1em}\\phantom{\\cdot}\\\\\n",
       "\\hline\n",
       "0.&\\hspace{-1em}667&\\hspace{-1em}\\phantom{\\cdot}&0.&\\hspace{-1em}333&\\hspace{-1em}\\phantom{\\cdot}\\\\\n",
       "\\end{array}\\right)\n",
       "\\]"
      ],
      "text/plain": [
       "StateSpace(array([[0.33333333]]), array([[0.66666667]]), array([[0.66666667]]), array([[0.33333333]]), 1)"
      ]
     },
     "execution_count": 118,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Pssd = c2d(Pss, 1, method='tustin')\n",
    "Pssd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\frac{0.3333 z + 0.3333}{z - 0.3333}\\quad dt = 1$$"
      ],
      "text/plain": [
       "TransferFunction(array([0.33333333, 0.33333333]), array([ 1.        , -0.33333333]), 1)"
      ]
     },
     "execution_count": 123,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Pssd2 = ss(1/3, 1/3, 4/3, 1/3, 1)\n",
    "P2 = tf(Pssd2)\n",
    "P2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ゼロ次ホールドと双一次変換"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "連続時間システム \n",
      "  1\n",
      "-----\n",
      "s + 1\n",
      "\n"
     ]
    }
   ],
   "source": [
    "P = tf([0, 1], [1, 1])\n",
    "print('連続時間システム',P)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "離散時間システム（zoh） \n",
      "  0.1813\n",
      "----------\n",
      "z - 0.8187\n",
      "\n",
      "dt = 0.2\n",
      "\n",
      "離散時間システム（tustin） \n",
      "0.09091 z + 0.09091\n",
      "-------------------\n",
      "    z - 0.8182\n",
      "\n",
      "dt = 0.2\n",
      "\n"
     ]
    }
   ],
   "source": [
    "ts = 0.2\n",
    "\n",
    "Pd1 = c2d(P, ts, method='zoh')\n",
    "print('離散時間システム（zoh）', Pd1)\n",
    "\n",
    "Pd2 = c2d(P, ts, method='tustin')\n",
    "print('離散時間システム（tustin）',Pd2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaYAAACcCAYAAAA05w98AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAq1UlEQVR4nO2de3hTVdb/vztNk6YtUqRAgbapw8CDFxhwRJEZxEG8gWL1FRWKCrQojhdEKQVagSIVEFScB1+1hYJcH7yA6Kv4U7zUC6jvCOhAfRWB3hgHWqDU3tKmWb8/msQ0TdKT5pyec3b253nyNCc552StZn+zzl5n77UZEUEgEAgEAq1gUNsAgUAgEAg8EYFJIBAIBJpCBCaBQCAQaAoRmAQCgUCgKURgEggEAoGmEIFJIBAIBJrCqLYBwTJjxgxKTk5W2wxZcTgcMBj4ukbg0afc3Nz1RJShth1yIHSkD3j0SYqOdBeY4uPjsWTJErXNkJWqqirEx8erbYas8OhTbm5uhdo2yIXQkT7g0ScpOtJdKI6IiFDbBNkxm81qmyA7PPrEE0JH+oBHn6Sgu8DEW7cWAIxG3XVcO0TrPm3duhUpKSkwGAxISUnB1q1b1TapSxE60gda90kpHWnbax80NzerbYLs1NXVwWKxqG2GrGjZp61bt+KBBx5AfX09AKC0tBQPPPAAAGDKlCmoqqpCeXl5m0dubq6aJsuO0JE+UMunrVu3Ijs7G2VlZUhOTkZeXh7S0tLa7RNIR2fPnm2no3nz5kkzgIh09cjMzCTeqKysVNsE2VHDpy1btpDVaiXGGFmtVtqyZUu7fRwOByUmJhKAdg+j0UhRUVHtXo+MjKQTJ04QgCWkAQ3I8RA60gdq6Sg6OrqNBqKjo9voqSMdeR/vev3777+XpCPd9Zh4zI2bTCa1TZCdrvbJ19Vbeno6PvroI1x44YUYO3YsbrnlFhw9ehQVFb7vvdrtdsyePRtJSUltHr179+Yu9SV0pA/k9ElKLwgAsrOz3TpyUV9fj4ceegiRkZG46667UFlZGVBHjz76aDsd9enTR3q76yhyae2RnZ0t03WBdmhoaFDbBNmR06dAPaGmpiZ68803qXv37j6v3gCQxWKhZcuWERGRzWajHj16+NzParUGtAMc9ZiEjvSBXD511AtqaWmhnTt3Ul5enl8dAaB58+YRUWuPSUkdqS6QYB8iBaEP5PLJl6CMRiNNmjSJiIiam5vJZDL5FRJjjBwOR4fn9E5V+IKnwCR0pA868klK+pqIyGq1+tTHBRdc4N6nZ8+eBIAiIiJ87pucnNzus5XSkeoCCfYhBKUPQhFURUUFvfXWW7Rw4UKf93y8BXXkyBFKSkoK6upNqqA9EYFJ24SbjjoKDL/++iu9/fbblJOTE7AX5KK4uJhqamqCCjhK6Uh1gQT7mD9/foeO643z58+rbYLsBPLJV8OPjIykESNGtLmh6u/KzdUT6uicUq7egoGnwCR0pG2k/OD76wVFR0e3ec9gMJDRaFT8wk0qXAYmHnPjjY2NapsgG1IadN++fX2KJCIigu655x5as2YN7d+/nxoaGvyKz5eglBQTEV+BSehIu0i9yGKM+b1wmzRpEq1evZq+/PJLqqur65ILN6lwGZh4TEFUVVWpbYIs+Gr8FouFHn/8cXrsscfo6NGjRESSe0H+zqllQenlIXSkXfxdjPXq1YueeOIJOnToEBER9erVSzO9oGAQgUkn8JIb9ycoV4B6++23iYi65H6QEojApG20riOp7ThQT8hsNtOmTZuIiGjjxo1ksVg0cdEWDFwGpqysLHn/SxqgurpabRNCor6+ngoLCwP2hGw2m3t/LfWCgoGnwCR01LVIafM2m402b97sc3IqAOrTp0+74eNauWgLBi4DE4+5ca3Pv/DV+I8fP0779+8notbAFBMTQwaDQTP3g5SAp8AkdNS1+Msm9O/fn7744gsiIrLb7dSrVy/q1q1bu4E/erhwkwqXgck1wYsnzpw5o7YJfvF3pTdo0CAaOnSoe78TJ07otickFZ4Ck9BR1xIoPWe1Wt1z7UpKSqilpUWXF25SkaIj3dVZafWLLxwOh9om+CUrK8tneZKamhrs3LnT/VpKSgrS0tKQn58Pq9UKxhisVivy8/N9lj0RqIvQkXxIqbDdv39/n8f27dsXH3zwARhjAACr1QqDwYC0tDSUlJTg9OnTKCkpCTsN6a5WHm81ywDtlbavqanBG2+8ga1bt+LkyZM+9zl16hQGDBjQ7vW0tDSkpaWhuroacXFxClsq6CxCR/IQqMJ2amoqdu7ciS1btvisKxcdHY1Vq1Zh8ODBfs+vtd+GrkJ3rZPHLyoqKkqVz/W80rNarcjNzUV6ejr69u2LjIwMlJeXo3v37j6P7WhZbrV8EkhD6Ege/BU8nTVrFvr164f77rsPP//8MxYuXIgVK1YEnU0IVx3prnXyuI5MfX19lzdA7yu9srIyLFmyBCaTCffeey8yMjJw1VVXYdu2bW32A1qv9PLy8gKeXw2fBNIROpKHsrIyn6/X1ta6dTR69Gh3qi4rKyuo84erjnTXYxK5cXnwdaUHAH369MG6deswcuRIMMY6fd9Iy/fNBEJHHSF1ZVZ/mYPExERs2rQJ11xzjTsodYZw1ZHuAlMoX7JW6ap8PxHho48+QnNzs98rPV+5cNeNWIfDIflGLI/3MHhC6Mg/rmxCaWkpiMh938gVnIgIn376Kerr65GXl9duzaTo6GisWLFCFlvCVUe68zoyMlJtE2QnOjpa0fM3Nzdj48aNuPTSS3HDDTdg9+7dfq/0Orp3JBWlfRKEhtCRf/zdN1q4cCG2b9+O4cOHY+zYsdi2bRvS0tJQWFio2EjUcNWRooGJMfYwY2wmY2wxYyzO43XGGDvMGDvGGPuFMVbDGBvgfG8CY6zC+ZjifU673a6kyarQ2Ngo27m8BzTMmDEDAwcOxPTp02E2m7F582ZMnDgReXl57Rq9lHtHUpHTp3BH6EgacrU5f9mEsrIyTJkyBc3NzVi3bp07+HQmoyCVcNWRYoMfGGOXA7iKiO5jjP0VwGIAc5xv/wHARCI67tx3AxEdY635hVEAkshPEpzHnKtcPxK+BjRs2LABAwYMwHvvvYebb77ZncJxiUfKUsudgccfPjUQOpKOXG0uOTkZpaWl7V6PjIzEjh07cNttt3VZii1cdaTkfzcVwBHn80PObQAAER3zENNwAAedb/0JwE0AShljPn8hRW7cP/4GNDQ3N2P8+PHt/ndKXumFa25cAVIhdCQJudpcXl4ezGZzm9fMZjM2bNiA22+/vUvbdrjqSMnh4gkA/tf5vBFAPz/73QXgBQAgokMA/swYGwjgPcZYCRF95bnz559/3kZUmZmZyM3Nhd1uh81mAwDExMTA4XCgoaEBwO95WtePtsVigcFgQF1dHYDWRmc0Gt3bJpMJJpMJtbW17m2z2Yza2loQESIjI2E2m1FfXw+HwwGj0YioqCj3tsFgQHR0NBobG2G3293bNpsNzc3NYIwhNjYWNpsNTU1NaGpqQmNjo/u5ywepPjkcDpjNZr8piPLyctTU1OjKJy1+TyohdCRTm3vzzTfxzDPPoKKiAv3790d2djbuu+8+t08um+644w588cUXKCwshN1uR2JiIhYuXIgbb7xR6KirdNRRzaLOPgA8DWC283kPAOU+9mEACvwcfw+Ax71fF1WR23L27Fm6/PLL6fnnnw9qUT2l0XKl584CFWrlCR1JJ1Cb66iOY11dHY0aNYoWLVpEREQOh8Ndv05NwlVHSvYTdwIY5nw+BMA7jLE4xpjncKAxAIpcG4wxT3u6A9jrfVIec+OdmezoutqJi4vDFVdcgaSkJMUHNAQDjxM4VULoSCKB2py/kXYLFiwA0KqTESNGuMtsMcY0ke4MVx0pFpiI6CCAYsZYBoBxAHIALINHjtz5/G2P7acYY58wxuYC+JmIDitln5boSACeI+2Sk5Nxxx13ICkpCSdOnABjDK+++iruvPNOTRVR1YKoeUDoSDqB2lygNPfhw63/njVr1rhTe1ohXHWkaEkiIlrl9dIjXu8/7rWd29E5vSez8UBsbKzf97xH2pWXl6O8vBwjR45sd4MW+L2IqtoE8kkQHEJH0gjU5vyNtLNYLJpuq1q2TUl0N+SjpaVFbRNkx3UT0xf+Rtr9+uuv6NfP331w9Qnkk0B9wk1HeXl5sFgsbV4zm80oKChASkqKwpZ1nnDVkQhMGsDfSBVylkPxhb/UhFZQcRSbQAI86ciV6o6Li/NZ146IEBER0aaielJSEtavX6+J7EIgwlVHugtM4UJ5eTnGjRvn9325SgcJBHqmo7p2p06dwsSJEzF58mQMHjwY33//PYgIZWVlmg9K4YzuAlO45MYbGhpw5MgRTJ8+XTMj7YIhXHPjeoEXHfkbbZednQ2gdVTbgQMH8Pzzz2P//v0YOnSoGmZ2mnDVke4CE08pCBeu7vrZs2fx4osvgogwaNAglJSUoLCwUDMj7YIhXFMQeoEXHflLaZeWlsLhcCAxMRHHjh3DnDlzEBER0cXWhU646kgEJhXxzo0//vjjmDt3Lo4caa1A41ogTMnSQUoRroLSC7zoKFBK+5tvvgGg71Vgw1VHugtMvOArN/7WW29h2bJluOyyy9Q2TyDQBf4mlS9fvhxXX321SlYJQkV3gYmXdWT85cZffvlllSySl5iYGLVNEASAFx2lpaVhzZo17qHg/fv3R35+PubPn6+yZfIQrjpSdIKtEvBSSiXQmi88EK7l+vUCLzr67rvvsHz5cjQ3N2PFihWYN28eV9USwlVHuusx8ZAbf/XVV10FNtvByzDwcJ0YqBe0riPPMly+5ia59vnLX/4Cu92OL774Aunp6VwFJSB8daS7wKRnmpub8cgjj2DWrFkYOnRou5noehgGLhAoTUdzk1paWjB//nxMnToVI0eOxHfffYeRI0eqbLVATnQXmPSaGz979ixuuukmvPTSS5g7dy4OHDiAgoIC3Q0Dl0q45sb1gpZ1FGhuUk1NDVJTU7Fy5Uo8+OCD+PDDD9GrVy8AfLY5Hn2Sgu7uMflLgWmZhoYGXH311SgpKcHGjRtx//33A/i94GpdXR13DZCXexi8omUdBbr/OmbMGPzrX//C2rVr8fe//71N6o7HNsejT1LQXWDS481Ai8WCuXPn4rLLLvM5hLWhoYG7wMSjTzyhZR35qwSenJyMuXPnIiEhAdddd12793lsczz6JAXdpfK0judN2969e7uHrc6cOVPMqxAIJOBrbpLZbEZeXh7S0tJ8BiUBX+guMGk5N+5907ayshKrVq3yOaLIE28R8gCPPvGElnXkveCl2WxG3759MXny5IDH8djmePRJCroLTFrG101bh8PhLigpEAikkZaWhmPHjsHhcODkyZM4cOAADAbxcxUu6O6bbm5uVtsEv3R20qyvhQD1Do8+8YSWdQQAq1evxo033oiGhgb07NkTPXr06PAYHtscjz5JQXeBSasQEbp16+bzPV4mzQoEXcGyZcuQmZmJnj17tlncTxA+6C4wabGhOhwOPPbYY6ipqWlnn5RJs94TbXmAR594Qi0dBaroQETIycnBU089hXvvvRdbt24N6l4Yj22OR5+koGjrZIw9DKAJQD8ALxJRtdf7ewEMBlAB4C9oDZSLARwFkEBEK32cU0mTg4aI8PDDD+OVV17BE088gcsvvxzZ2dkoKytDcnKyeyRRIHjMnfPok1rwoiPX4CBXespV0QEApkyZggULFmDlypXIyMjAq6++GnQb4rHN8eiTJIhIkQeAywFscj7/K4AXvN5PBXCV12tzAKQ7n+cAuM37vJmZmaQVHA4HzZkzhwBQVlYWORyOTp2nsrJSZsvUh0efACwhhfTi78GTjqxWKwFo97BarZSbm0sAaNasWdTS0tKp8/PY5nj0SYqOlOwxpQI44nx+CMBmp2Bc/A3AfzHGvgXwEBGdAnA7gEyPY/4LwG7PkxYVFbW52svMzERubi7sdru74GFMTAwcDgcaGhoA/D7k0nWlZrFYYDAYUFdXB6B1joTRaHRvm0wmmEwm1NbWurfNZjNqa2tBRIiMjITZbMaiRYvwwgsv4IEHHsDixYtx7tw5OBwOGAwGREdHo7GxEXa73b1ts9nQ3NwMxhhiY2Nhs9nQ1NSE8+fPIzY2Fk1NTe6FwWJiYlTxqb6+Hg6HA0ajEVFRUe5t4ZNqC7alghMdBVptdvHixZg8eTKee+45VFdXizbHsU+S6ChytQY4JErZz+uYfAAznc+NAGw+9okA8BSAD53bPwMY6Hw+DsD/8z5m/vz5CsTw4Fm/fj0BoBkzZnT6Cs9FTU2NTFZpBx59Qog9pnDXkb8eEwCaNGkSNTc3h3R+Htscjz5J0ZHUBOaLjLG/M8biJO4PAKcAuGaHdQNw2nsHImohoqcBXCD1GK3kXFNTU7FkyRLk5+eHbJMWB3SECo8+yUBY68hXRQeLxYLU1FRs2bIl5DbDY5vj0ScpSG2dUwBsA3APY2wpY2yIhGN2AhjmfD4EwDuMsTjGWCQAMMYMzr9RAD73OGa45zHeJ1Vr/oXnaCKr1Yo9e/Zg8eLFiIiICPncri41T/DokwyEtY68KzokJSWhoKAAu3btgslkCvn8PLY5Hn2SgtRwHAvgQgCjAfwRwEnG2K0AdhHRj74OIKKDjLFixlgGgGS03oRdBqCIMfYJgC8YYx8B+AnA087D/hvAMsbYdLSOQnqzk37JivdoorKyMvdoIl6WqRB0CWGtI6BVL+PHj8dFF12EcePGCf0IfCI1MH0H4FcAy4joPcB9pXYIwFB/BxHRKq+XHvF4fomP/W34/aatT+TooQRLoPVh5BCWHFeLWoNHn2QgrHXkokePHti2bRuGDRsm63l5bHM8+iQFqam8/yaiq11i8jj2fxSwKSBqCKqzpYakwmPj49EnGQhrHZ07dw579+4FAIwfPx79+vWT9fw8tjkefZKCpMBERM/6eM1ORAvlNykwXT1st66uzu/sc7lKDbmGafIEjz6FSjjryGazITU1FbfddhtOn243FkMWeGxzPPokhfAc8iERIsL06dPR1NQEs9nsnjcASCs1JBAIfq+O8vnnn2Pbtm3o3bu32iYJNI42xl4HQVemIPLy8vDGG29g1apVWL9+vXs0kdVqRX5+vmw3bnnsrvPoE0/IqaNA9e8AYO3atVi/fj1ycnI6XFMpFHhsczz6JAXd9Zi6KjDt3r0bTz31FKZOnYonn3wSjDHFRhCZzWZFzqsmPPrEE3LpKFD9u7S0NHz88ceYM2cObrvtNuTm5srymf7gsc3x6JMUdNdj6orc+OHDhzF16lSMGDEC+fn5ihe85DGPzKNPPCGXjgKNWD1+/DjuuusuDB48GJs3b1Z8Ui+PbY5Hn6Sgu8CkNNXV1UhNTUW3bt2wa9euLik731qlgy949EnQnkAjVlNTUwG0Zh/8rVUmJzy2OR59koLuUnlKX3Wlp6ejtLQURUVF6N+/v6Kf5SKYNWf0Ao8+8YRcOkpOTkZpaWm716Ojo3H48GF88MEHGDBggCyf1RE8tjkefZKC7gKT0rWjZsyYgeuvvx6jRo1S9HM84TGPzKNPPCGXjvLy8trcYwJag9Ls2bMRHx+PG264QZbPkQKPbY5Hn6Sgu8CkVI2vhoYGWCwWTJgwQZHzB6K+vh5RUVFd/rlKwqNPPCGXjlwDglyLYyYmJmL58uWqlBrisc3x6JMUdHePSc6cq+cw127dumHmzJmynTsYHA6HKp+rJDz6xBNy6igtLQ0lJSWorq6GxWJBZWWlbOcOBh7bHI8+SUF3gUmu3LhrmGtpaSmICC0tLdiyZUu7ORhdAY+l7Xn0iSeUuFdrNBoxevRoXHHFFbKfW+rn8waPPklBd4FJri/K1zDXxsZGZGdny3L+YOCxq86jTzwh9w8eESE6Ohrr1q3DX//6V1nPLRUe2xyPPklBd4FJrty40oVZg8E7QPIAjz7xhJz3ak+cOIFRo0ahuLhYtnN2Bh7bHI8+SUF3gUmu3HhSUpLP1+UqzBoMPOaRefSJJ+TSkd1ux9SpU1FcXNxuddquhsc2x6NPUtBdYJKrCsOll17a7jW1CrNqZbl4OeHRJ56QS0fPPPMM9u3bh5dffhkpKSmynLOz8NjmePRJCrq7sybHhLNPPvkEe/bswbhx43D06FGUlZUhOTkZeXl5qgxzDfVKs7m5GRUVFWhsbJTJotAhIpw6dUptM4ImKioKiYmJ3E9slMO/b7/9FkuXLkVaWhqmTJkig1WhIXSkHULVke4Ck91uD+n42tpapKenY+DAgdi9e7fq6QegddBFKDc5Kyoq0K1bN6SkpChe108qdrtddyOKiAhnzpxBRUUFLrroIrXNUZRQdWSz2TB9+nQkJCRg7dq1MlkVGkJH2kAOHemunxhqznX+/PkoLS1FYWGhJoISEPqPRGNjI3r27KkZMQH6rPHFGEPPnj01dcWsFFJ0FGg5i6VLl6K4uBgFBQWIi4tT0FLpCB1pAzl0pGgoZow9DKAJQD8ALxJRtcd7zwK4HgABSCOiH52vTwDwqnO3eUS0zeucnbanqKgIL730EmbPnq3akFZfyJFH1pKYAO3ZIxUt2q2GjgItZzF48GCsXLkS06ZNw8033xy6gzIhdKQdQrVbsR4TY+xyAFcRUQGAjwEs9njvTwA2EdFwAAUAnna+zgCMApBERIneYgI6nxuvr6/HjBkzMGDAAM2tPKuVnpucdOZH4j//+U+X1lbTA2rpyN9yFgsXLsT06dPRp08fPP/886E5JzNCR63woCMlU3mpAI44nx9ybgMAiOh7Ijrs3NwH4N/O538CcBOAUsaYz1EIne2uWywWLFiwAIWFhYiJienUOZTCc8l2Xggm5Xr8+HEAQEJCAvbs2aOUSXolFSroyN98vvLycsyfPx/r169Hjx49JDnQVQgd8aMjJVN5CQD+1/m8Ea1pCF9cD+B5ACCiQwD+zBgbCOA9xlgJEX3lufOnn37appuYmZmJ3Nxc2O12d8OMiYmBw+FAQ0MDgNagxBhzrw9TV1cHg8GAuro6AK0VfI1Go3vbZDLBZDK5F+kymUwwm82ora0FESEyMhJmsxn19fVwOBwwGo2IiopybxsMBkRHR6OxsRF2u929bbPZ0NzcDMYYYmNjYbPZ0NTUhPPnz8NsNqOpqcm9gFtMTExAn1xXh/X19WhpaUFLSwsYY2hpaQHQeqUlZdvhcICIwBiDwWAIeZuI4HA40NLS0mYbaF011Xv7nXfewaFDh5Cdne2+OnRN/pTqgxI+VVVVtfueVEIVHfXv3x8VFRXtPiQxMdF9NS50JHSklI6UDEynALj61t0AnPbewZmK+CcRlXi+TkRHGWOLAIwA0EZQY8aMwTfffOPzA70XI3M1yhtuuAHp6enthoJ7LwLove09wse7BL33+8Fuu87HGENUVJTPEUW+fPIkOjoalZWV7qWyDQYDrr322nbn8eaWW27B3Llz3ftPmzYN06ZNw9mzZ3HnnXe22fezzz5zP/dOLRgMBrS0tGDNmjXo3bs3tmzZgoKCAqxfvx4DBgzAxx9/jKysLFgsFkyaNAkZGRn4+OOP0b17d6xduxbvv/8+Tp8+jX379uH48ePYtWsXduzYgdzcXJw7dw5xcXHYu3cv3nnnHfz444/IyMhASUkJNm/ejPXr1+Ozzz7DsWPHsGPHDiQkJGDfvn1Ys2YNvv32W/e+W7dude/72muvwWQyYfv27Xj22WcxePBgnz7Fx8e3+55UQhUdrVixot1yFgaDAddff32b/w0gdCR0JL+OlEzl7QQwzPl8CIB3GGNxjLFIAGCMXQagOxF9xhi7gDEWyxjztKc7gL3eJzWZTEEZUV1d7W6wWiU2NlZtE0KisLAQffv2xb333ouHH34YWVlZGDRoEO677z5Mnz4dM2bMwEUXXQSLxYIrr7wSb7zxBvbs2YOoqCiMGjUKw4YNw5gxYzB69GjU1NTAYrFg4MCBMJvNWLlyJcaOHYuioiKMGzfO/ZmjR492P8/IyEB6ejoyMjKQnJyMZ555xu++O3fuRFJSEtatW4cLLriga/5BoaGKjtLS0pCfnw+r1QrGGJKSkjB06NAuXacsWISO+NGRYj0mIjrIGCtmjGUASAaQA2AZgCLG2CkAmwBUOdMJBOBKAIsYY2MAvA/gO4/8uRtXd1Mq8fHx2Lu3nS41hc1mk/Wq3PPKLNj94+Pjgz7+4MGD7vTOxIkT8fLLL7uvUEeMGIEffvgBQGvK4cILLwTg++a753wNz31jYmICpgC+/fZbtziuvPLKgPNqcnNzMWvWLPTo0QNbtmwJxk1VUFNHaWlpbbIMWh+6LHTUCg86UnQeExGtIqJ1RLSIiM4R0SNE9AYRfU5EKUR0hfMxglrJJaKxRLSaiD71dU4pgnLNv2CMITExEdu2bdP0sEsV71/IwsUXX4yCggK0tLTg1KlTsNls+PLLLwEAv/32W8Ch+a5cttSbvK6b9lVVVe5jhg8f7v68mpoa9+f52reurg5ff/01rrrqKrz22mud8LbrUUtHLtauXYvjx4+DMSZ0pCBCR7+jrynFEvCef3Hy5En3/As1yg2FAzNnzkRRUREGDRqECRMmYNeuXZg1axZWr16NhoYGvPLKKygtLUVpaSn2798Pm82G8+fPY9++fRgyZAjWrFmDv/3tb/j1119RXl6O0tJSfP311zh79iwqKirw448/oqqqClOmTMGtt96Km266CRMnTkRdXR1++eUXbNiwAQsXLsSxY8dQVVWFefPmAYDPfXNycnDHHXeAiHDrrbeq/J/TPl9++SUeffRRVFRUYMWKFWqbwzVCRx4Qka4e2dnZFAir1UpoTWm0eVit1oDHqUlDQ0NIxxcXF8tkiXy0tLSobUKn8ff/BLCENKABOR4d6YiIyG6309ChQyk5OZlqa2s73F9thI60RSg60l2PqaMUhJbWWZJKU1OTpgdndAbS+P2IcEdKKm/dunX44Ycf8Prrr2tu7p8vhI74QXe18joSVP/+/X2+rsY6S1LRe27cF+G6joxe6EhH1dXVyMnJwejRo9sNe9YqQkf8oLvA1BEDBgxo95pa6ywJBHrl6aefxpkzZ/Diiy9qesCDgE90F5gC1fg6cOAAPv/8c4wfP949/8JqtSI/P1/TAx/0kCYJFtdERYE2CaSjn376Cf/4xz+Qnp6O4cOHd6FVoSF0xA+6u8fkr2tLRJg9ezbi4+Oxbds2dO/evYst6zyhluvXIuGaG9cLgVJETz75JCwWC5YtW9aFFoWO0BE/6K7H5C83/v333+Prr79GXl6eroISIIpPCroefzr65ZdfsHfvXixatAh9+vTpYqtCQ+iIH3QXmPwxbNgw/PTTT5gxY4bapgh8cOTIEVxxxRUoKSlR2xRBAP74xz/i559/xqOPPqq2KYIg4GGpC090F5h85cYrKysBAH/4wx90mZPt6tx4oJVJ5cL7e7j00ks7vZaWQH5WrVrV7rt36Sg5OVntwrWdggcduZaucNHR7xlPS114orvA5J1zraurw5AhQ5CTk6OSRaHTld11V2WM0tJSEJF7ZVK5g5Ov3Lgef+x4xvO7b2pqwsiRI3XdU9K7jo4dO4bs7Ow2rwW6x/Tuu++2+Tw9XpT7Q3eByfsGZ0REBObMmYNbbrlFJYtCx7U2TFfgb2VSb0EEw8aNG7F27Vps3LgRQ4cOda/ts2nTJtx///0oLi527/vRRx9h8uTJePDBBzv9eQL58PzuH3vsMfeaZXpE7zoqKirCwYMHsXLlSlx77bUgIuTl5WHatGkAWmsWbt++HaNGjUJjYyPeffddfPfddygqKsKGDRswceJENDQ0YP78+XjwwQeRlZWFP//5zzh58mQorqqC7gKTN1FRUcjKysLIkSPVNkUXKFEZY8KECXjkkUcQFxeH7t27w2QyYeDAgW3K9bu48sorsX37dq7SDnqnrKwMJpMJs2fPxnXXXae2ObpACR2NHTsWCQkJuPvuuwG0Fmb1XGbk9ddfx5AhQ1BQUICIiIiglrrQG7oLTJ73KVatWoUdO3aoaI08uFbR7Ar8VcAIpTJGr169UFZWhnnz5mH79u04cOCAe2SkZ7l+AO7XeUo76J24uDhs2LBBbTNCRu868oXn5OYlS5bg7rvvxvLly9uNquzsUhdaRXeByUV5eTlycnLw4Ycfqm2KrsjLy2sn4FArY7S0tCAtLQ0vvPACEhMTkZKSgq++al0w1V+5/nCdn6E1LBYL6uvr8d5776ltiq5QQkeupStMJpN7efozZ864752ZzWb88MMPMJlMeP/994Ne6kJP6C4wudaxX758OYgIixcvVtmi0PHOVSuJ98qkclTGWLp0KWw2G86dO4cNGzagX79+KCsrw+rVq1FQUIBXXnkFZWVlKC8vx1dffYWffvoJ58+fx/79+2X0TBAsVqsV48aNQ1NTE5YuXaq2OSGjdx316dMH586dw+7du5GQkIA777wTP//8M06fPo1///vfeOihh1BYWIg+ffrgmmuuwZAhQ/Duu++iqKgI+/bta7PURXFxsXupi3/+85/66zV1VH5ca4/MzEwqKyujyMhImjVrVnB12DVKZWVlSMdrsVx/U1OT2iZ0mnBY9iIzM5MqKyspJiaGJk+eHNL/SysIHWmLUHSkux6T0WjE8uXLAQALFixQ2Rp5sFgsapsgOwaD7ppWWGE0GvHcc8+hvr4eTz31lNrmyILQET/orlbeb7/9hnXr1iE9PV3TS1kEA4+NT1Sk1jYNDQ0oKCjAPffcg4svvlhtc2RB6IgfFA1MjLGHATQB6AfgRSKq9njvbwCucm6+T0Q/MMYiASwGcBRAAhGt9D7n2rVrAYAbMQGtk4R5u9praWnh8odCDZTQ0b59+7jqLQFCRzyhmMeMscsBXEVEBQA+RqtQXO9FAXgWwEoALwBYy1ovDR4BcIKIXgMQyRi7zd/5FyxYoEgpHb1CYpSbLGjt/6iUjg4cOMBVb0kutPb965VQ/49KhuJUAEeczw85t12MAlDlvBdmAxANIBnA7QAO+zmmDaHOstYSoZbqiYqKwpkzZzQlKj1e5RERzpw5o7XluVOhgI7sdjtXvSVA6EgryKEjJVN5CQD+1/m8Ea1pCM/3znpsu973fN37mHaUlpZi3rx5yM3Nhd1ud5e9j4mJgcPhcJcocc03cA0ntVgsMBgM7rkCZrMZRqPRvW0ymWAymVBbW+veNpvNqK2tBREhMjISZrMZ9fX1cDgcMBqNiIqKcm8bDAZER0ejsbERdrvdvW2z2dDc3AzGGGJjY2Gz2dDU1ASbzYbIyEg0NTW5h3XGxMRI9ik6OhrV1dU4ffq0e04DY8w9zyGYbdeomGC3gVYReY6siYiI8Pu+r21XPt1zuzM+hOJTVFQU4uLiUFVV1e57UglFdHTJJZfgkksucW+7ykgJHQkdaUFHSgamU2i9ggOAbgBO+3nP8/1Ax7TDarXi2Wef/f0k3bq1ed+72rD3hDjvfLT3tnfE974i834/2G3X+aqqqhAVFeXzCkOqT1pbO6eqqgrx8fFqmyELKhefVURHEyZMwJEjR7xfbj1A6EgzhKuOlOwn7gQwzPl8CIB3GGNxzhuzXwLoDwCMMTOAZiI65jxmuOcx/k4e6ixrLeEZXHmBR58AXKvCZyqiI1dlDp7gsc3x6BMk6IgpmU9ljGUCOIfWvPcLAJ4GUEREbzDGxqNVNBEA9hDRQae4lgEoBtAHwEryMpAxRlarFXl5eSHNstYSrm4wT3DsU5eP31VKR5x+P2qbISsc+xRQR4oGJiUQgtIHHPvExcQSoSN9wLFPAXWkuwm2AA4yxvym+HTKtYyxz9Q2QmZ49Gmi2gbIiNCRPuDRpw51pLsek0AgEAj4Rn+D5AUCgUDANSIwCQQCgUBTiMAkEAgEAk2hq8DEGHuYMTaTMbaYMRantj2hwhgbzBh7lzE2Rm1b5IAxdgFj7C3G2P8xxt5hjMV0fJS2YYwZGWOrGGNfMsb2OOvT6RqhI20jdKSjwMQCFLPUK0T0fwB+A8DFEGQA1wCYAeBiAHUApqtrjiwkAlhKRH8FEAngjyrbExJCR7og7HWkm8CEwMUs9YzO1jz2DxH9DxGdd06Q+RrAv9W2KVSIqISIfmOMRQP4gogOd3iQtkmF0JGmETrSV2AKqsCrQHVSEKCklJ5gjF0AIBfAA4yxP6ltT4gIHemLFIShjvQUmIIq8CpQD8bY7QBeICK72rbIARHVEFEmgCwAM9W2J0SEjnRCOOtIT4GpXTFL9UwR+IMxdgOAfxFRGWOsv9r2yMxPAH5U24gQETrSAeGuI92UJHIWpyxmjGWgtZhljto2hQpjzApgEIC/MMa+IaIGtW0KBcbY4wDmAKhkjBkAfAtglqpGhQhj7HoAcwHsAEAA8tW1KDSEjrSP0JEoSSQQCAQCjaGnVJ5AIBAIwgARmAQCgUCgKURgEggEAoGmEIFJIBAIBJpCBCaBQCAQaAoRmAQCgUCgKURgCiMYY39Q2waBQO8IHSmPCExhAmNsAIA8te0QCPSM0FHXIAJT+DAGwHDG2N1qGyIQ6Bihoy5ABKbw4RMA/yGiHWobIhDoGKGjLkAEJoFAIBBoChGYwgcCYHAWhRQIBJ1D6KgLEP/c8OEUgB4AHlDbEIFAxwgddQGiurhAIBAINIXoMQkEAoFAU4jAJBAIBAJNIQKTQCAQCDSFCEwCgUAg0BQiMAkEAoFAU4jAJBAIBAJNIQKTQCAQCDSFCEwCgUAg0BT/HyiU6oWNpM1MAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x165.6 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1,2, figsize=(6, 2.3))\n",
    "\n",
    "Tc = np.arange(0, 3, 0.01)\n",
    "y, t = step(P, Tc)\n",
    "ax[0].plot(t, y, ls='-.', label='continuous', color='k')\n",
    "ax[1].plot(t, y, ls='-.', label='continuous', color='k')\n",
    "\n",
    "\n",
    "Td = np.arange(0, 3, ts)\n",
    "y, t = step(Pd1, Td)\n",
    "ax[0].plot(t, y, ls='', marker='o', label='zoh', color='k')\n",
    "\n",
    "y, t = step(Pd2, Td)\n",
    "ax[1].plot(t, y, ls='', marker='o', label='tustin', color='k')\n",
    "\n",
    "for i in [0, 1]:\n",
    "    plot_set(ax[i], 't', 'y', 'best')\n",
    "    ax[i].set_xlim([0, 3])\n",
    "\n",
    "fig.tight_layout()\n",
    "# fig.savefig(\"c2d_step.pdf\", transparent=True, bbox_inches=\"tight\", pad_inches=0.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/minami/opt/anaconda3/lib/python3.8/site-packages/control/timeresp.py:293: UserWarning: return_x specified for a transfer function system. Internal conversion to state space used; results may meaningless.\n",
      "  warnings.warn(\n",
      "/Users/minami/opt/anaconda3/lib/python3.8/site-packages/control/timeresp.py:293: UserWarning: return_x specified for a transfer function system. Internal conversion to state space used; results may meaningless.\n",
      "  warnings.warn(\n",
      "/Users/minami/opt/anaconda3/lib/python3.8/site-packages/control/timeresp.py:293: UserWarning: return_x specified for a transfer function system. Internal conversion to state space used; results may meaningless.\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaYAAACcCAYAAAA05w98AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAwFklEQVR4nO2deXRTdfr/35+0TZeUfYfSFlAUDjsKSmWRGaRTQBHoIHTEKlAXGPGnlGJBKTDFA3qEYapIkU3aqSCLUhFGFK2iKDDAV1lEQCgUylIKdE+T3Of3B02mS5LeJDe5Sz+vc3JObu69n/s8yX3nuZ/nszEiAofD4XA4SkEntwEcDofD4VSHByYOh8PhKAoemDgcDoejKHhg4nA4HI6i4IGJw+FwOIqCByYOh8PhKAp/uQ1wh+eee47Cw8PlNkMyBEGATqetZwQt+rRw4cK1RDRNbjukgutI+WjRJzE6UmVgatmyJVJSUuQ2QzIKCgrQsmVLuc2QFC36tHDhwjy5bZASriPlo0WfxOhIlaHYz89PbhMkJTAwUG4TJEeLPmkNriPlo0WfxKDKwKS1qq2/vyorrk7Rok9ag+tI+WjRJzGo8s40mUxymyAppaWlcpsgOVr0SWtwHSkfLfokBlUGJg6Hw+FoF1UGJq3lxvV6vdwmSI4WfdIaXEfKR4s+iYEHJgWgxZtPiz5pDa4j5aNFn8SgysBUWVkptwmSUlJSIrcJkqNFn7QG15Hy0aJPYlBlYOJwOByOdlFlYOIpCOWjRZ+0BteR8tGiT2LggUkBaHEQnRZ90hpcR8pHiz6JQZWBiefGlY8WfdIaXEfKR4s+iUGVgUlrEJHcJkiOFn3iKBst3nNa9EkMqgxMWptKJSAgQG4TJEeLPmkNriPlo0WfxKDKO1Nr80dpMY+sRZ+0BteR8tGiT2KQJTAxxmYwxqYzxhYwxprW2reMMXaUMXaEMdbN3vlam+OrrKxMbhMkR4s+KQ2uo5po8Z7Tok9i8HlgYoz1AzCQiNYA+BrAgmr7egP4iIj6AlgDYLG9MrSWdxUEQW4TJEeLPikJrqO6aPGe06JPYpCjxjQWwImq98eqtgEARPR/RHS8avNHAFfsFaC13LjWUiqANn1SGGPBdVQDLd5zWvRJDHJ43RbAoar3FQDaOzhuBIB37e3IyckBY8y2nZiYiIULF8JsNsNoNAIADAYDBEFAeXk5ACAkJATA/6rGwcHB0Ol0tmnlAwMD4e/vb9vW6/XQ6/W27pp6vR6BgYEoKSkBESEgIACBgYEoKyuDIAjw9/dHUFCQbVun0yEkJAQVFRUwm822baPRCJPJBMYYQkNDYTQaUVZWBrPZjNDQUFRWVtq68RoMBu6TgnxSGFxHDeCe06JPYmC+rs4zxhYDKCCifzLGmgH4hYg61jqmN4BmRPStvTKSkpJo6dKl3jfWRxQWFqJ58+ZymyEpWvSJMbaQiFLktgPgOrKHFu85LfokRkdy1OW3A+hT9b4ngJ2MsaaMsQAAYIz1ANCEiL5ljDVmjIXWLoDnxpWPFn1SGFxHtdDiPadFn8Tg81QeER1ljJ1kjE0DEA5gPoB/AMhhjF0D8BGAgqoUAwEYULuM6ukHLaC1XD+gTZ+UBNdRXbR4z2nRJzHI0rJGRG/X+mhmtfeR9Z2vtUFn1nywltCiT0qD66gmWrzntOiTGFQZjs1ms2RlZWZmIjIyEjqdDpGRkcjMzJSsbLFUVFT4/JreRos+aQ0pdaQEtHjPadEnMaiyL6JUedfMzEwkJCTYeq7k5uYiISEBABAXFyfJNcSgtT8IQJs+aQ2ttV9o8Z7Tok9iUGWNSarc+Lx58+qMrC4rK8O8efMkKV8sWswjy+WTEmrAaoG3MSkfLfokBlXWmKTKjV+8eNGlz72FFvPIcviklBqwWuBtTMpHiz6JQZXhWKrqbXh4uEufewsFDt6sgTu1EF/6dO3aNcyePRvTp09XRA1YLWgtTaR0HbmDL30qLCxEcnIyxo4dizfffBO3b9/22bVro8rAJFVuPDU1tc4TSUhICFJTUyUpXyxKnkzTWgvJzc0FEdlqIfUFJ1/5dPDgQfTs2RMrV660jXqvja9rwGpBa21MStaRu/jKp+PHj6NXr15YtmwZzpw5g9TUVPTr1w+nT5/2yfVro8rAJBVxcXFIT09HREQEGGOIiIhAenq6z9M+Ss71u9sO5wufTCYTnnrqKYSGhuLo0aOIiIiwe5yva8ANESW07SlZR+7iC58EQUBcXBwEQcDPP/+MEydOYP/+/SgpKcHo0aNRWFjodRvqQESqe82bN4+0REVFhdwmOIQxRrg7QLPGizHm9Dxf+XT48GE6efIkERFlZGRQSEhIDTtDQkIoIyNDkmsBSCEF3P9SverTUUZGBkVERBBjjCIiIhx+j97+3sWiZB25i698OnHiBP33v/+1bWdkZFCbNm0IABkMBkl/SzE6kl0c7rzmzp0r2ZekBO7cuSO3CQ6JiIiwG5giIiKcnudtn0pKSux+bv0zBUD+/v60fv16ya6ptcDkTEeuBBt37xGpUbKO3EUOHdn77YODg336gKfKVJ7FYnHrvMTERERFRUlsjedYZwxWAo0aNcKYMWPw66+/AnC/Hc6bPhERoqKiMGvWrDr74uLicOHCBRQVFaG4uBjx8fFes0PtONORsxTuokWL0L9/f1vnCaX0blWSjtq3b4+RI0fi0KFD9R/sBG/7FBMTg2eeeabGZ/Z++/Lycp92IlJlYHKFX3/91daA2LdvX/To0cO2LyEhAV999ZVcpimC77//HuPHj7c1hE+dOhU//fQTBg4ciM8//1wx7XDVMZlMmDBhAgYPHuzwmEaNGiEoKAjl5eW4evWqD63TBs6CTZ8+fdCtWzfbWkGhoXXmhwXQcNr2MjMz0a5dO5s+MjMzER8fj19//RVRUVH4+OOP5TbRLhaLBWPGjMGjjz5a43NFPGjUV6VS4ktsG9OePXsoKCiIXn/99Tr7bt++Td26dSOdTkfLli0jQRCcliU23+4O5eXlkpXlKqtXr6bu3btTQUGB7bP8/Hx64IEHKDg4mA4fPuxWuXL6ZMVisdA999xD8fHxdfa583tCY6k8ZzoSm54rLy+3e6wcbUxy3HPOUp43b96kIUOGUEBAAOXk5LhVvhw+Ofrtg4ODJSlfjI5kF4c7LzFtTEeOHCGDwUB9+vSh69ev2z2muLiYYmNjCQDNnz/fYXDyduOu3Llxew2sV69epY4dO9L999/vlji85VNubi5t2bKFjEajqOPT09Np3759NT5z9/fUWmBytY3JUTtDeXk5DRkyxHZceHi4z4MSkTw6qi+AFxYWUteuXSk8PJyKiopcLt9bPl2/fp0yMjKorKyszj57v71er6fExERJrq3ZwFTfF7R27Vry9/cnANShQwenIrFYLDRt2jQCQIsXL7Z7jLcbd2/cuCFJOa4wd+5c2rZtm9Nj9uzZQwBowYIFLpfvLZ/mz59POp2OLl686HYZ7v6eWgtM9elow4YNFBAQIEpHgiDQq6++SgBo9uzZTsv1FnLoyN59hFq9Vn/88UdijNErr7zicvne8untt98mAHTixAm7+72ZIWqQgSkjI8MWlMQ+DVssFpoyZQoBsHucu12mxeJrQeXk5BAASkpKqvfY2NhYMhgMdPXqVZeu4Q2fTCYTtW/fnmJiYlw67+zZs7RkyRJbjdjd37OhBaY5c+YQANq5c6fT46wIgkAzZ84kAPTee++JOkdKfK2jI0eOOAxMtR9ypk+fTgEBAXT+/HmXruENnwRBoPvuu48GDRrk0nm5ubm0ePFiMpvNHl1fs4EpOTnZodPWvveuPg0bjUYaNmwY6fX6Ovlgb9eY7FWnvYXJZKLu3btTZGQkrVu3rt6notOnT5Ofnx/NmjXLpet4w6e9e/cSANq6datL56WnpxMAW3sZrzHVr6P//ve/xBijhIQEp99JbcxmM40ZM4Z0Oh3t2rXLpXM9xZc6EgSBBg4cSI0aNaLg4OB6H4QvXbpEQUFBdts7neENnw4ePEgA6MMPP3TpvM2bNxMA+uqrrzy6vmYDk6PcuMlkElW1dkRhYSHdf//91KxZMzp9+rTtc2+3MbmTe3aXjz76iADQrFmzRPu0ePFi2rJli0vX8YZPCQkJFBoa6rJYCwsLKSAggF577TUi4m1M1pcjHQmCQFFRUdSqVSu6deuW0+/EHsXFxdSvXz8yGAx09OjRGvu8mSLypY527txJAGjt2rWifXrnnXdo48aNLl3HGz4lJiaSv78/3bx506XzysrKKDQ0lKZOnerR9TUbmBylIFatWiW6au2Ic+fOUXR0NF2+fLnG594UlK9SECaTibp27Uq9e/em8PBwVbWbmUwmatmyJT311FNunT9mzBgKCwsji8VCRLxXHjnRUVZWFgGgNWvW1PudOOLKlSs0cuRI+uOPP2yfefsBz1c6EgSB+vXrR507dyaTyeTSua7ed1L7JAgCRUZGUnR0tFvn/+1vf6OmTZuK7nxkjwYXmCZOnEhdu3aV7OY3m80e/QBi8ZWgNm3aRABo27ZtLrezFBYW0sqVK0Xnl6X2yZrGq6/DhiMyMzMJAH333Xdu29BQAtPzzz9PvXv39rgtwYrFYnHYrVzJD0OOyM7OttWWXCEjI0NU2q86UvtkTeOtW7fOrfN37drlUrujPTQbmJzlxouKiiSp3Vhz5c8880y9Y5w8xRe5cbPZTPfddx/17NmTLBaLy38SW7duJQD07bffirqe1D49//zzZDAY3C63uLiYgoOD6aWXXnLbBq0Fpvp0JAWCINDkyZPpySef9HonIl/oSBAEeuCBB6hTp05UWVnp0rnuBGYpfcrIyKDGjRsTAAoLC3Prf9FoNFLz5s1p8uTJbtshRkeqnPnhrm//w2QyIT8/H8DdEf/WaWkEQcCFCxfcmqXAz88PDz74IAYOHOj1GX59sfzA5s2bcfr0abz55pvQ6XQuTzU0ZswYHD9+HEOHDhV1PSl9IiJkZ2cjOjoawcHBbpURGhqKv/zlL/jss880t9yDu9TWkSAIyMvLA3BXR1LAGMOAAQPw0EMPoWPHjnaPkWqGCF/8rl988QUOHz6MefPmubzQojszKkjlk3X5mqKiIgBAXl6eqOVraqPX6zF27Fjs2rXLu9Ml1Re5lPiqnYJYt24dBQYGOuyTLwVXrlwRfazceWR7jB49mnr06GFrYyFST7vZ9evX6cEHH/R4QtaNGzcSADp48KBb50NjNSZUPa1bf/dt27ZRQEAA/fzzz259P/XhTirLFXyho7i4OIqMjHS5tkTkXo1JKp+kTKN+9tlnBIC+/PJLt2wRoyPZxeHOq3ZgOnv2LKWkpHgt5Xb06FEKCQmhlStX1nusOw28vhCUyWSiCxcueFRGaWkpTZo0SVRu3VvjLzyhoKCA/Pz8yFkKyxlaDEzV78+LFy/SG2+8IVnbUm1Onz5NwcHB1KRJE1EPQ0p8wLNYLC6PRbIi53+DlGnUsrIyCgkJoRkzZrhli2YDk7t/LO5iNBrpiSeeIACUmJjoVLiuPJl4s8ZixWw2O1wiwlUEQaAePXrQww8/XO+xpaWlklyTiCTtgDJx4kRatGiRW+dqNTC5++TsKmazmf72t78RAHrxxRed1jrc+ROX8p6rjSAIkrS71V6WpT7NS+WT1B1Ppk6dSnPmzHHrXM0HJovFQjNmzKAjR4649QW5QmVlJb344osEgIYPH+6w9iH2ycRXi6t9/PHH1KpVqxrjsjzhH//4BwGod0ogqQRVUFBABoPB5fEf3kDLgQkAHThwQOJvrC4Wi4WSkpIIAD388MMO70t3/ki9GZj27NlDzZo1k+y/Ji0tjQDHUwJZkcqndevW1fku5Zhol0jDgcmayrP2FPv3v/8t3bdWDx9++CGFhoaSwWCg5OTkOlVtsYLydtdZK4cOHaLp06fXaFvyhNOnTxMAWr58udPjpEpBXLlyhV555ZU6AzU9QRCEGrOpi0XrgWnVqlUufyfukpWVRc2aNaOgoCB69dVXKT8/v8Z+d1JPUqbyamcz3nrrLXruueckq73n5+cTY4xSUlKcHieVTzdv3qRRo0ZRmzZtJMvQeFNHsovDnVdiYiIJgkB9+vShe++912s5cUecP3+e/vrXvxJjjAICAmj06NGUl5dHRHcnvhRTE/J211lv0qtXL4qKinJ6jBwTaoolJiaGHnnkEZfP02pg0ul01LRpU58vTX7lyhWaMmUK6XQ68vPzo5EjR9LZs2eJiNwaAC7VPeerbMaQIUOoR48eTo9Rso4mTZpUr/32EKMjVXYX9/f3x65du3Ds2DEkJyfDz8/Pp9ePjIzE5s2bcfz4ccyaNQunTp1CkyZNAAAnT54EESE8PByMMbRs2RLjxo1D9+7da3T9dNRFVqqusxaLBW+++SYuXbokSXnViY2NxQ8//IDLly87PMbdbt3VsVgs2L9/v22lVKmYMmUKpk6dKmmZaqVNmzYQBAEpKSkIDAz06bXbtWuHjRs34rfffkNSUhLOnTsHg8EAABg4cGCd4wMCAjB16tQ6K+9mZmYiMjISrVu3RmRkpMtdoGvjbPVeKZkwYQKOHz+O3377zeExUuiIiLB//37bgqlSMXHiRMycOdM73fTri1xKfCUnJ9OAAQPc7rbpTfbs2UOpqam27TFjxtieulq0aEEzZ86kEydOeP2pzDrLg6tz3Inh1KlTBMBpL0UpBgYeOHCAAFBWVpbHZUkBNFZjSk5Opscee4xat27t1fYZd/jmm29o7NixtnSawWCw6aRp06aUkJBAR48e9YqOfJXNuHz5MgGOl9shkkZHJ0+eJAD0wQcfeFyWFIjRkezicOc1btw4AkDp6emSfmHeQBAEysvLo02bNtFTTz1Fer2eGGM0efJk+uc//+mVXnmVlZXUuXNn6tOnj2RtS7Xp0aMHDR482OF+KVIQb7zxBul0OpcnmxTDpUuXKDs726VztBaYJk2aRADo7bffdul7kIv8/HzKysqip59+2jYeqva4KOvLk7ZaX7X/EhFFRUVRr169HO6XQkfvvPMOAaDc3FyPy6pNWloatWrVSvI5J2URBIAZAKYDWACgaa19jwKYW/Xq5eB88vPzU0RPLVe5ceMGzZ07l4KDg6lx48a0fv16hyvsusvq1asJAH3++eeSlludlJQUYozVabSWsgt8//793WoLEsOsWbMoMDCQiouLRZ+jtMDkqY46duxIbdq0UVxtSQy3bt2ilJQUuwHE09qNr9qYiIiWL19OAOj333+3u1+KwDR8+HC32oLqIyMjw7aQpCvfk2SBCUCYmONEltUPwEdV7x8BsLzaviAAhwAwAIEAvgPA7JQha3dHKTh37hwNHjyYANCePXvqPV7sH35RURG1b9+eHn74Ya/O8Xf69Gl66623avTKkVLQV65cIQC0ZMkSKc228c033xDg2qSwngYmJepoxYoVrn95CqJDhw5eqd1s2LDBtuCoN5eKv3TpEi1atKjOA54VT8dO3blzh/z9/d0ec+QMb65rJlYE2wC8VPupzJ0XgEUAkqrehwI4X23fcAC7q20fBhBhpwyvVq99hcVioa+//tqWR3bUK8qVP/zExEQCfDMmpTZSpkCs4y6OHTsmvaF0dyaM5s2b09NPPy36HAkCk6J0FBoaSuXl5S5/d0rCW7Ub63i9//znPxJZ6h6etjFt27aNAPGTL7uCN1eC9oc4JgMIBvAUY6w9gE+I6FeR59amLe4+zQFABYD2tfYVVtu27s91VFhubi7mzJmDhQsXwmw2w2g0AgAMBgMEQUB5eTkA2CYstfa2CQ4Ohk6nQ2lpKQAgMDAQ/v7+tm29Xg+9Xo+SkhLbdmBgIEpKSkBECAgIQGBgIMrKyiAIAvz9/REUFGTb1ul0CAkJQUVFBcxms23baDTCZDKBMYbQ0FA88MADuHLlCn7//XckJCQgKysL999/v80Hs9mMuXPn2u0l9Prrr+PJJ5+07cvNzcXy5csxefJk3HPPPSguLvaqT+Xl5fjuu+8QFRWF0NBQp5NUFhQU1PCpvt9px44daNeuHbp06YLy8nKv/E4jRoxAdnY2iouLYbFYnP5OVns9RFE6Cg8Pr9HrKzExUXU6io6Oxvvvv4833ngDeXl5aN26NYxGI1q3bm33ntu6dSuWLFmCvLw8dOjQAfPmzcOUKVNq+HTjxg2kpqZi9OjR6Nevn9d1VFlZif3796N3795o0aJFjXvuxo0baNKkCUJDQ1FZWWmbOFWsjrZv347GjRujd+/ekusoLCzMbq/fDh06oLCw0DMd1Re57gY4tABwL4BMAD8DeB5AMoBuYs6vVdZiALOq3jcDcKnavj8B2FFt+/8AdLFThiZqTFZu3LhBp06doujoaLsD1sQ8mRiNRurfvz8ZDAYKCwvz6jRHVqxru1jb+qSqMVVWVlLjxo1p+vTpXrD6f+zYsYMA0L59+0QdD89rTIrSkaP1mNTKjRs36Pz58zRy5Mg6C30Siatdmc1mGjp0KBkMBq90FrCHddC6vV6unrQxCYJA7dq1o9jYWE/Mc4g3V4IWK4ILAA4AGFXtMx2AX8ScX6usvgDWV70fAuA9AE0BBOBuPvxg1b5AAIcdlKH6Nqbq3Llzp8a20WikdevW2dqIxPzh79u3jxhjFBgY6JNGW6K7N/73339vG+AsVVrF2v7z6aefesNsG6WlpRQcHEwvv/yyqOMlCEyK0pGjpdXVSm0dmc1mWrt2re3+FKOjgwcPUmBgIG3YsMGXptP+/fvtroZb2ydXOHLkCAHwqi8ZGRm2wdAGg8G3vfIAzLHzmT+AJWLOt3NuIoBpuJsnbwYgDUBs1b4YAElVT5J9HZzv9dqAL6md5//www8JAMXFxVF5ebnoP3xvNQS7ghS98hYsWEABAQG2HnPenOz28ccfp/DwcFEdRSQITIrS0bx58zz67pRGbR1Zpyx74oknqKSkRHSbyPnz572+OKhYPGkDfPfddwkAXb16VUKL7PP0009Ts2bNRC01L1lgUtpLiymI6giCQKmpqQSABg4cSKdOnXL457xhwwbKzMwkInmmOSoqKqJXXnmFdu3a5dQnVxAEwTY1jbe77lo7WYiZnNPTwKS0l9Z1RES0cuVK0ul01KtXL2rXrp3DB7dt27bR6tWrZbD6LhUVFZSYmFhnQLxUOvI21k4WYtLiPDCpBEc339atW6lZs2ak1+vp73//O/3222+2JzlBEEgQBBoxYgTFxMSQIAg+HRhoxWw2U5s2berksZW4wJk9rl+/TiEhIfTRRx/VeywPTMrG0T33xRdfUOvWrYkxZusCXvshZ8KECfTII4/4fN5NK4Ig0D333EMjRoyo8bmS58qrTklJCTVu3Jjee++9eo/VbGDSem68OlevXqX4+HjbQLbmzZtTixYt6NChQ0R094/VOuOxLwcGVueFF14gg8FQo2uru7nxjIwMmjp1qi2F4YtaoNh0idYCU0PS0c2bN+mFF14gPz8/2z2k0+nI+h3cvHlTkul/PCEpKYn8/f1rzHTiro4+/fRTmjJlCt2+fVsq8+pFSh2pchJXX0/a6m2cTZ7Zpk0brF+/Hrm5uUhLS0NsbCzGjx9vO6dVq1bQ6/UAgLi4OKSnpyMiIgKMMURERCA9PR1xcXFetX/cuHEoLS3Fl19+KconZ+Tl5eHYsWO287092S0ABAUFAUCdyUG1TkPSUfPmzbFq1SpcvXoVq1atwosvvoipU6di3Lhxtv1STJjqCePHj4fZbEZ2drbtM3d1dPnyZRw6dAiNGjWSyrx6kVRH9UUuJb60loJwZ00TJVFZWUnNmjWrMVjVE5+qNzz7ohZoMplo0KBBNH/+fKfHQWM1Jq4jZSEIAnXs2JHGjBlj+0wqHfkCa9PCzJkznR4nRkeqrDFpjbu/lXoJCAjAE088gezsbNsAQHd8qqioAAAwxmyf+aIW6O/vj/79+6Nz586SlcnxPWrXEWMM48aNw5dffoni4mIA7uuIiGroyBcwxtC3b1/bBAGeoMrApNOp0myHBAQEyG2Cx4wbNw63b9/GN998A8A9n6ZNm4ahQ4fW+TwuLg4XLlyAIAi4cOGCV1KTK1euxLPPPit5uUqG60h5jB8/HkajEcnJyW6vMfXaa6+hX79+3lknqR6WLl2KGTNmeFyOKu9Mf3+xMympA18v0OYNRowYgdDQUGzfvh2A6z4ZjUZkZ2ejS5cu3jBPFGVlZfj+++9lu76v4TpSHoMGDUKjRo2watUq5ObmgoiQm5uLhIQEUcFJEATs2LEDnTt3lu3Bo6KiAjk5OR6VocrAJPVKjHJTex48NRIUFIRRo0bh008/hcVicdmnffv2oaioCOPHj/eShfUzZ84cjBw5Erdu3ZLNBl/CdaQ8/Pz8QER1OhCIXUH3p59+Qn5+vqw6Wrx4Mf785z/jypUrbpehysCk9lxybeSocnuDKVOmYNKkSSgvL3fZp08++QSNGjXCn//8Zy9ZVz/Tpk1DeXk5NmzYIJsNvoTrSJlYJ1atjaNJkq1Ly+t0OkRHR8PPzw+jRo3ypolOee6552CxWJCenu52GaoMTFrLjWslpRITE4MVK1YgNDTUJZ9KS0vxySefIDY2VtZ0TJ8+fTBo0CC8//77mvmTcwbXkTJxZYhEZmYmEhISbGk/a6eJzz//3Ks2OqNLly6Ijo5Genq627VyVd6ZWrkBrVj7/2sBQRCwb98+l9Iq27dvR0lJCeLj471nmEj+/ve/4+zZs9ixY4fcpngdriNlkpqaalu2wkpISAhSU1PrHDtv3rw6WrNYLKLSft7k5ZdfRn5+PjZt2uTW+aoMTDw3rlxOnjyJP/3pT/joo49En7NhwwZ07twZjzzyiBctE0dsbCzuu+8+LFiwQPO1Jq4jZRIXF4e33nrLtu1siISzNdDkZOTIkXjwwQexePFi2xASV1BlYOK5ceXSo0cPZGdni+7Sff78eezbtw/x8fE+H3dhDz8/PyxcuBAnTpzAmjVr5DbHq3AdKZeXX34Zu3btwvnz550OkfDFzCjuwBjDokWLcOHCBaxYscLl81UZmJTwByYlWsv1jx49GgaDQdSx69atg7+/v6LGEP31r3/FsGHDkJSUhLy8PLnN8RpcR8omJiYGjRs3dnqMK2k/XxMdHY2xY8ciJSUFZ86ccelcVSaZ7Q2kM5lMyMvLs80eoCaICNeuXZPbDBtBQUEICwvzaMDiZ599huzsbGzbts3pH+CCBQswatQohIWFuX0tqWGMYfXq1ejfvz9iY2Px448/au5PHOA68jZS6Gjv3r3YuHEjsrOz7c5taK1JJScn49KlSwgPD0dqaqrX58cUy7/+9S/06dMHTz75JI4dOya6XVOVgclsNtf5LC8vD40aNUJkZKTq/kTMZrNiGqKJCDdv3kReXh46derkdjl37tzBjh07sGfPHvzlL3+xe4zV74ceesjt63iLrl27YvPmzQgJCVHd/SQWriPvIZWOysrKsHv3bmzZsgWTJk2ye8zEiRMVE4hqExYWhq1bt6KsrMyl30aVdV97ueSKigq0aNFCdWIClJXrZ4yhRYsWHj8xx8bG4t5778XLL79st6wLFy6gU6dO2Lt3r0fX8SYxMTEYNmyY3GZ4Da4j7yGVjmJiYtC3b1/Mnj3b1hW8OteuXUOXLl3w6aefenQdbzJs2DDExMS4dI4qA5Mj0ahRTIDy7JbCnuDgYKSlpeHs2bNITEy0e0y3bt1w3333eXwtjntwHXkXKewJCAjA+++/j/z8fMyYMaNO8BUEAd27d5dk4lQlocrApIXJGqtTvdH26tWreOyxx2S0RhpCQkLw2GOP4dVXX0VaWhrefPNN21RFZWVliIyMxDPPPIMhQ4ZAp9O5PFElx3O0rCNAG1oKCQnBQw89hJSUFGzatAmvvvoqTCYTjEYjSkpK0K5dO+zevVtzgUkZCVkXsZcbVzPWWbM7d+6Mtm3bYvfu3XKb5DFGoxFBQUFYunQpCgsLsXjxYqxduxYVFRUYPHgwYmNjkZCQYBt7Yp2oEoBi8+VaQ4s60ul0+OOPPzSjJauO5s+fj4KCAqxYsQJZWVmwWCzo2bMn9u7dq7kFHwGV1pi0NF4BALKzs2vUFrRwo1kHb/r7+2PdunXYunUrhg0bhrFjx+K1116zO2Jd7ESVHGnQmo6ISHNasupIp9Nh5cqV+PzzzzFixAiMGjUKr7/+uur9c4Qqa0xiEdNwPXr0aMyePdt2fHx8POLj41FQUIAJEybYjvv222+dlmOxWLBixQq0bt0aGRkZWLNmDdavX49OnTrh66+/RlJSEoKDgxEbG4tp06bh66+/RpMmTZCWloZdu3ahoKAAOTk5+OOPP7Bjxw5s3rwZCxcuxK1bt9C0aVN89dVX2LlzJ06dOoVp06bhwoUL2LRpE9auXYtvv/0WZ86cwebNm9G2bVv8+OOPWLFiBQ4ePGj32I0bN0Kv1yMrKwvLli3zShqgen6dMYbx48fXmPFYqSPWOXXxpY4A97W0YsUKZGdn4/r165rRUu12qlGjRsk6QauvUGWNSa/Xy21CHdatW4d27drh6aefxowZM5CUlISuXbtiypQpePbZZ/Hcc8+hU6dOCA4OxoABA/DJJ59g9+7dCAoKwqBBg9CnTx8MHToUgwcPRlFREYKDg3HvvfciMDAQS5cuxfDhw5GTk1Nj9u3Bgwfb3k+bNg1Tp07FtGnTEB4ejiVLljg8dvv27ejYsSM+/PDDegfwuUtoaKjT/Uodsd6QUKKOAPe1FBISojkt1acjraLKGlPttUocIebpzNHxLVu2dOn8o0eP2hpaH3/8caxatcp2oz744IP45ZdfANxNLTRv3hyA/cbn6n39qx9rMBiczjl18OBB2/UGDBiAtLQ0h8cuXLgQL7zwApo1a4aMjAzRPrqC0Wh0OlN4ampqjTYmQDkj1hsKStQR4L6WavdY04KW6tORVlFljUmsoHxJt27dsGbNGlgsFly7dg1GoxH79+8HABQXFzudoJSIIAiC6Jy/tdG6oKDAdk7fvn1t1ysqKrJdz96xpaWl+OmnnzBw4EBs3LjRDW/rp76JG+Pi4pCeno6IiAgwxpxOVMnxDkrUEeC+lgRBAGNMU1pyZwJULaDKGpMSmT59OnJyctC1a1eMGjUKO3bswEsvvYR33nkH5eXl+OCDD5Cbm4vc3FwcOHAARqMRd+7cwY8//oiePXsiLS0Njz76KPLz83Hp0iXk5ubip59+QmFhIfLy8nDq1CkUFBRg8uTJGDNmDKKjo/H444+jtLQUZ8+exfr165GcnIxz586hoKAAc+bMAQC7x86fPx/jxo0DEWHMmDGyfWdxcXE8EHHq4K6WDhw4gJ49e2LFihUNTkuag4hU95o3bx7V5uTJk3U+UwsWi0VuE+rg6fdZXl4ukSXKAUAKKeD+l+rFdeR9uI7qIkZHPJWnAEhBU6lIRUNNQagJriPl01B1xAOTAtDaeBKg4QpKTXAdKZ+GqiNVBiYOh8PhaBefBibG2KOMsblVr1529k9mjB1gjJ1jjI1zVI7W5vjS4uhtsQsFclyH68g+XEfawWe98hhjQQCWARgAQA9gL2NsaFVjGBhjAQDMRPQwY6wngL0AttsrS2tVdi3mxrU2D5tS4DpyDNeRdvBljWkQgIKqjhlGACEAbMP8ichERFuq3v8K4HdHBfHcuPIxGo1ym6BVuI4cwHWkHbxSY2KMpQC4p9bHkwFkVduuANAeQK6d8x8F8E9H5efk5NSYQyoxMRHx8fGwWCy2m9PPz882cBX435T41bcZYzZxit0WBAFEBMYYdDqdS9unTp1CfHw8Pv74Y3Tq1Ak6nQ5EZLtGbZvr2/amTxaLBbdv30ZQUBDKyspsMzeHhISgoqICZrPZtm00GmEymcAYQ2hoqG1cCXB3SpXKykpbI67BYIDZbLYJzmAwQBAElJeXA7g7+wMA24wQwcHB0Ol0KC0tBQAEBgbC39/ftq3X66HX61FSUmLbDgwMRElJCYgIAQEBCAwMtPng7+/vtk++hutI/DYR4fLly5g+fTr+85//cB2pXUf19SeX6gXgTwB2VNv+PwBd7BzXHMDTzspKTk6u0zfenfECGRkZFBERQYwxioiIoIyMDJfLcJWHHnqIzp8/X+MzLY6/KCsrk8gS5QAFjGNqCDo6d+6cW8dbLBYym80eXVtquI7qIkZHvkzl7QfQAQAYY4EATER0jjEWwBhrWvV5cwCPEdGmqu329goiCXLJmZmZSEhIQG5uLojIth6QtxerszfvlRT+KA0tplUUgqZ1dO7cOZeWPqm+zAURaa4DREPVkc86PxCRkTGWwhhLAuAHYHrVricBDGGMpQLYA8DEGJsNoAWA4fbKkqJB0Nl6QO5Mk7NhwwaUlJQgNDQU7777LkaOHInQ0NAaU/V3794dALB3717s27cPjRs3xurVqyEIguYEVV5e3mB7FHkTresoJycHR48exdKlS7F792588803WLJkCc6cOYMNGzYgLS0NLVq0wL/+9S/s27evxjIXZ86cwc6dOx0uc9GhQweP/fU1DVVHPu0uTkRfENFSIlpCREerPttCRDOJKJ+IehPRA1WvTkR03lu2SL0e0KhRozBz5kw0bdoUTZo0gV6vrzNVv5UBAwYgKytL9atrcuRByzoaPnw42rZti4kTJwK4ux5RVFSUbf+WLVvQs2dPrFmzBn5+fjWWuYiKinK6zAVHPahygK0U4y+kXg+oVatWuHjxIubMmYOsrCwcOXLE7lT9ANCkSRMA/xt3YW181RLWxleOclGijuojJSUFEydOxFtvvVWnV2H19aVcWeZCyTRUHWnvH1EkqampdX50T9YDslgsiIuLw/LlyxEWFobIyMh6p+rXYtsSp2EhtY6svdv0er2tx1j1ZSYCAwPxyy+/QK/X44svvnB5mQuOOlBlYDKZTB6XIfV6QIsWLYLRaMStW7ewfv16tG/fHhcvXsQ777yDNWvW4IMPPsDFixdx6dIl/PDDDzh9+rRtqn4tiqp2uwNHeShRR23atMGtW7fw2WefoW3btpgwYQLOnDmD69ev48qVK3jxxRexbt06tGnTBkOGDEHPnj2RnZ2NnJwc/PDDDzWWuTh58qRtmYvDhw+rstbUYHVUX7c9Jb4SExPrdEFU83T9lZWVcptQB0+/zxs3bkhkiXKAArqLS/niOvI+XEd1EaMjVdaYqi+ZrAW02MYUHBwstwmceuA6Uj4NVUeq/CWrj1bXAlrzB9Dmn4TW0Np9pzV/gIarI1V6LUVuXElobc4yALaGa45y4TpSPg1VR6oMTI4g3stNEvj32LDhv7808O/RfVQZmOzNkhAUFISbN2+q8mZQUnWdiHDz5k0EBQV5VI69qZc4yoLryHtwHXmGKls/7d2AYWFhyMvLw40bN2SwyDOISFH58aCgIISFhXlUhtYa1rUI15F34TpyH1V6bS83HhAQgE6dOslgjecUFBSgZcuWcpshKaWlpQ22R5Fa4DpSPg1VR8qp+7rADz/8ILcJkrJs2TK5TZAcLfoEYJjcBkgJ15Hy0aJPEKEjpsZcMmOM1Gi3IxhjqszpO0PDPiknV+QhXEfKR8M+OdWRKmtMHA6Hw9EuqmxjAnCUMbZTbiMkZBhj7Fu5jZAYLfr0uNwGSAzXkfLRok/16kiVqTwOh8PhaBeeyuNwOByOouCBicPhcDiKggcmDofD4SgKVQUmxtgMxth0xtgCxlhTue2RAsbY/YyxbMbYULlt8RTGWGPG2DbG2G+MsZ2MMYPcNnkKY8yfMfY2Y2w/Y2w3Y8yzOWYUANeRsuE6UlFgYoz1AzCQiNYA+BrAAplNkgQi+g1AMQAtjI8ZAuA5AN0AlAJ4Vl5zJCEMwCIiegRAAIB7ZLbHI7iOVEGD15FqAhOAsQBOVL0/VrWtFdS35rMdiOhzIrpTNWrzJwBX5LbJU4joAhEVM8ZCAHxPRMfltslDxoLrSNFwHakrMLUFUFj1vgJAexlt4dRPJABNjJFhjDUGsBBAAmOst9z2eAjXkbqIRAPUkZoC0zUAIVXvGwG4LqMtHCcwxp4EsJyIzHLbIgVEVEREiQCSAEyX2x4P4TpSCQ1ZR2oKTNsB9Kl63xMaeYrQGoyxxwD8SkQXGWMd5LZHYk4DOCW3ER7CdaQCGrqOVDMlEREdZYydZIxNAxAOYL7cNkkBYywCQFcAUYyxn4moXG6b3IUx9gqA/wfgBmNMB+AggBdkNcpDGGMjAMwGsBkAAUiX1yLP4DpSPlxHfEoiDofD4SgMNaXyOBwOh9MA4IGJw+FwOIqCByYOh8PhKAoemDgcDoejKHhg4nA4HI6i4IGJw+FwOIqCB6YGBGOss9w2cDhqh+vI+/DA1EBgjHUBkCq3HRyOmuE68g08MDUchgLoyxibKLchHI6K4TryATwwNRz2AbhKRJvlNoTDUTFcRz6AByYOh8PhKAoemBoOBEBXNSkkh8NxD64jH8C/3IbDNQDNACTIbQiHo2K4jnwAn12cw+FwOIqC15g4HA6Hoyh4YOJwOByOouCBicPhcDiKggcmDofD4SgKHpg4HA6Hoyh4YOJwOByOouCBicPhcDiKggcmDofD4SiK/w/f5CaE2rHnDgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x165.6 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1,2, figsize=(6, 2.3))\n",
    "\n",
    "Tc = np.arange(0, 3, 0.01)\n",
    "Uc = 0.5 * np.sin(6*Tc) + 0.5 * np.cos(8*Tc)\n",
    "\n",
    "y, t, x0 = lsim(P, Uc, Tc)\n",
    "ax[0].plot(t, y, ls='-.', label='continuous', color='k')\n",
    "ax[1].plot(t, y, ls='-.', label='continuous', color='k')\n",
    "\n",
    "\n",
    "T = np.arange(0, 3, ts)\n",
    "U = 0.5 * np.sin(6*T) + 0.5 * np.cos(8*T)\n",
    "\n",
    "y, t, x0 = lsim(Pd1, U, T)\n",
    "ax[0].plot(t, y, ls='', marker='o', label='zoh', color='k')\n",
    "\n",
    "y, t, x0 = lsim(Pd2, U, T)\n",
    "ax[1].plot(t, y, ls='', marker='o', label='tustin', color='k')\n",
    "\n",
    "for i in [0, 1]:\n",
    "    plot_set(ax[i], 't', 'y', 'best')\n",
    "    ax[i].set_xlim([0, 3])\n",
    "    ax[i].set_ylim([-0.3, 0.3])\n",
    "\n",
    "fig.tight_layout()\n",
    "# fig.savefig(\"c2d_time.pdf\", transparent=True, bbox_inches=\"tight\", pad_inches=0.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 状態フィードバック"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 例14.4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A = [[0.9 1. ]\n",
      "     [0.  1.1]]\n",
      "\n",
      "B = [[0.]\n",
      "     [1.]]\n",
      "\n",
      "C = [[1. 0.]\n",
      "     [0. 1.]]\n",
      "\n",
      "D = [[0.]\n",
      "     [0.]]\n",
      "\n",
      "dt = 0.1\n",
      "\n"
     ]
    }
   ],
   "source": [
    "ts = 0.1\n",
    "A = [[0.9, 1], [0, 1.1]]\n",
    "B = [[0],[1]]\n",
    "C = np.eye(2)\n",
    "D = np.zeros([2,1])\n",
    "P = ss(A,B,C,D, ts)\n",
    "print(P)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[-0.02, -0.5 ]])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pole = [0.7, 0.8]\n",
    "F = - acker(P.A, P.B, pole)\n",
    "F"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.8 0.7]\n"
     ]
    }
   ],
   "source": [
    "G = ss(P.A+P.B*F, P.B, P.C, P.D, ts)\n",
    "print(G.pole())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaYAAACcCAYAAAA05w98AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAYXUlEQVR4nO2df4wjZ3nHP493195d+46964UmIQoEAk2LgtIIJSgQLlFToPyBDlWCFlVUpUkrIKqSVltSlSq3IkHtBYmiCokfQkpFRSXaJjSUk9orUUhIk6MtgZZEacOBgtLSXFb3K+u99a7PT/9Y29lz7F3PeO35zvr5SKvz6xmPP2PP9x7POzPvmLsTBEEQBCoUshYIgiAIgo1EYQqCIAikiMIUBEEQSBGFKQiCIJAiClMQBEEgRRSmIAiCQIrJrAXScvPNN/sll1yStUaww1lYWPiSu9+ctcewiBwFoyBpjnJbmPbu3cvBgwez1ngZi4uL7Nu3L2uNrqi6qXoBLCwsPJe1wzCJHCVH1U3VC5LnKLOuPDO7wsy+bmb7u0y70czuaP69KQu/IMgDkaNgJ5LZHpO7P21mLwK28XkzmwYOAdcAReCIme33jiEqJiYmRuaahGKxmLVCT1TdVL3yQORo9Ki6qXqlIeuTH1a7PHcdsOjr1IBZ4NLOmSJQyVF1U/XKEZGjEaLqpuqVBsVjTBcCJza0V4CLgWc3znT33Xdz9913d13A3NwczzzzDOVymUajwdmzZwGYnZ0FYHl5GYCZmRkKhQLVahWAUqnE5ORku10sFikWiywtLbXbpVKJpaUl3J2pqSlKpRLLy8s0Gg0mJydZXl5menqaRqNBoVBgdnaWlZUV6vV6u12r1VhbW8PMqFQq1Go1VlfX/2+pVCqsrq622+VymXq9Tq1Wa7fTrtPp06e54IILEq/T9PR0uz2MdarVaszNzY30e+p3nXJMXzk6cuTIeTman59nYWFh27Y5iBxFjpLnSLEwPc/6r7sWu4DjnTPNz89z6NChrgsws/MOApbL5fOmt76wFjMzM5u2p6enz2uXSqWe0+v1Onv37t309UmW12LXrl3ntdOu0+7du1M5JG0nWafFxUXK5fJIv6c07ZzRV47279/P0aNHuy5gu7a5FpGjyFG/ZN2V18bM5sxsCvg28KrmcyVgzd2Pdc4fXRDJUXVT9cojkaPho+qm6pWGzPaYzOzVwBuAt5rZUeAu4Fvu/jdmdtDMPgZMALd0e71qoDp/XSih6qbqlQciR6NH1U3VKw1ZnpX3LOsHaFvcumHaYeDwZq9v9bOqsbS0JLuBqLqpeuWByNHoUXVT9UqDTFfeTkH5xouqbqpeQXYobxOqbqpeachtYSoUNNWnpqayVuiJqpuq1zgQOUqOqpuqVxo0t8o+mJxUPKFQu59X1U3VaxyIHCVH1U3VKw25LUxra2tZK3Sldc2AIqpuql7jQOQoOapuql5pyG1hUu1PbTQaWSv0RNVN1WsciBwlR9VN1SsNuS1Mqn3jql0joOum6jUORI6So+qm6pUGza2yD1S/BOXRAlTdVL3GgchRclTdVL3SkNvCFH3jyVF1U/UaByJHyVF1U/VKQ24LU/SNJ0fVTdVrHIgcJUfVTdUrDbktTGa29UwZoNpnD7puql7jQOQoOapuql5pyO2aqF5M1jmSrxKqbqpe40DkKDmqbqpeachtYarX61krdGVlZSVrhZ6ouql6jQORo+Souql6pSG3hUm1P1U16KDrpuo1DkSOkqPqpuqVhtwWpugbT46qm6rXOBA5So6qm6pXGnK7JtE3nhxVN1WvcSBylBxVN1WvNOS2MKnutqa5v/2oUHVT9RoHIkfJUXVT9UpDbguTat+46gWLoOum6jUORI6So+qm6pWG3BYmVVT77EHXTdUryA7lbULVTdUrDZkNlGVmHwVWgYuBz7j7qY7p/wxcATwHvNXdz22cXiwWR2SajEqlkrVCT1TdVL3yQORo9Ki6qXqlIZM9JjO7GrjW3b8IfBO4s2P6AeCP3f0Sd39LZ5gAzp172VMSKPfzqrqpeqkTOcoGVTdVrzRktcd0AHiy+fh7wJeB2zdMvxH4VTP7DvBhd3++cwEPPvjgebuu8/PzLCwstA/mLi4uUi6XaTQanD17FnjprJXWYIczMzMUCgWq1SqwfgfIycnJdrtYLFIsFllaWmq3S6USS0tLuDtTU1OUSiWWl5dpNBpMTk6yvLxMvV6n0WhQKBSYnZ1lZWWFer3ebtdqNdbW1jAzKpUKtVqN1dVVYP1Xz+rqartdLpep1+vtjW6QdTp9+nR7PZKs0/T0dLs9jHWq1WpMTEyM9Hvqd53EOcAQc7Qd2xxEjiJHKXLk7iP/A74A3NJ8PAnUuswzAfwJ8E/dljE/P++9WF+tbHjhhRcye++tUHVT9XJ3Bw56Bhnp52/YOcoS5W1C1U3Vyz15jrI6+eF5oHXS/S7geOcM7n7O3T8B7O62gOgbT46qm6pXDogcZYCqm6pXGrIqTPcBVzUfXwk8YGZzZjYFYGaF5r/TwMPdFqDaN97azVZE1U3VKwdEjjJA1U3VKw2ZHGNy9yfM7Ckzuxm4FPg4cBfwLTN7EHjEzI4A/wV8otsyIlDJUXVT9VIncpQNqm6qXmnI7HRxd7+n46lbNzz+hUGWvWfPnk3P6d+zZw8nTpwY5C2CQIJh5igIsiKzwjQom43xtVXRGeaFaOVyeWjLHhRVN1WvcUB1rDzlbULVTdUrDbkd+UF1KBXVscdA103VaxyIHCVH1U3VKw25LUyqfePK176ouql6jQORo+Souql6pSG3hSkIgiDYmeS2MEXfeHJU3VS9xoHIUXJU3VS90pDbwrR+MbEeqn32oOum6jUORI6So+qm6pWG3BYm1QN9rXGqFFF1U/UaByJHyVF1U/VKQ24LUxAEQbAzyW1hUu0bb43oq4iqm6rXOBA5So6qm6pXGnJbmIIgCIKdSW4Lk+r97Vv3PlFE1U3VaxyIHCVH1U3VKw25LUxBEATBziS3hWlyUnOYv5mZmawVeqLqpuo1DkSOkqPqpuqVhtwWpmEOxDoIhYLuR6rqpuo1DkSOkqPqpuqVhtyuiWrfeLVazVqhJ6puql7jQOQoOapuql5pyG1hCoIgCHYmqQqTme3bbpGkTExMZK3QlVKplLVCT1TdVL22EzOrmNn7zOy6ZvsyM7shW6vIURpU3VS90rDlkU8zexcvL2DvAn5vkDc2s48Cq8DFwGfc/dSGaTcC1zabh939Pzpfr9qfqnowGXTdVL22mTuAo8AbzewtwJ8DnwceGmShkaPRo+qm6pWGftbkF4H/7XjuFYO8qZldDVzr7h80s7cBdwK3N6dNA4eAa4AicMTM9nvHaJOD9I0P89br1WpV9uwYVTdVr23m34EfuPvXzWwGeA8w0EpnnaNhorxNqLqpeqWhn8J0r7v/dOMTZnZkwPc9ADzZfPw94Ms0AwVcByw2A1Qzs1ngUuDZAd+zTZa3Xg/GlqeBXwc+6e5nga+Z2dyAyzxAhjnKK3v37uXkyZNZawSb0E9hut3Mvubu/9L8VXbU3Tv3oJJyIfCvzccrrHdDbJy2sXK0pp8XqIcffvi8AjI/P8/CwgL1er19J8dyuUyj0WiPutsaS6p1hfTMzAyFQqF9NkupVGrvDi8uLlIsFikWiywtLQFQLBYplUosLS3h7kxNTVEqlVheXqbRaLRfe+LECRqNBoVCgdnZWVZWVqjX6+12rVZjbW0NM6NSqVCr1VhdXQWgUqmwurrabpfL5W1Zp2q1SrVaTbVO09PT7fYw1qnRaFCtVlOtU2sdhrVO28iHgPsBzOx64HF3v3fAZUrnaJDvBzbP0eWXX566uMzNzbXfK22O9uzZI5mjcrksmaO5ublE31E/hekp1n/twXof+fuBv0r0Li/neaA14uAu4HiPad2mA3DTTTfx+OOPd134rl27zmt33kCrc7DDbru/+/a9dH7H9PT0edM6DzJunL6ysvKy+bdqb7a8FoOu08zMDJVKpb3spA7DXKfWZ5ZmnQZxSNoekCd5KUePMwY5GuT72WqvZs+ePdtyL6nIkWaO+jny+WrgL83sNtaPN80lfpeXcx9wVfPxlcADZjZnZlPAt4FXAZhZCVhz92OdC2j9alCj9WtDEVW3YXrt3bsXM0v9t41EjhJw8uRJ3L3nX9pjwNvBOOZo1PRTmH4I/DbwHOvdEQOfk+juTwBPmdnNwE3Ax4G7gAPuXgMOmtnHgD8Abhn0/YJ8M0hxATb9D26rv20kcrSBrb7TpF0/wc6in668s8A5d/9bM/s34N3b8cbufk/HU7dumHYYOLzZ61WvvygWi1kr9ETV7fWvfz2nTp3qOX27um0yJnK0gdYeUS/OnDkzQptkqOZI1SsNW+4xuft9wEXN5iuA8iazjwzVwqR8kVuWblv9QlbtttkuIkfJiBwlR9UrDX1dXefuP2j++/0uv9AyQbVvXLmfN8tjOdC7S+2ZZ54ZmpcS45ajzbaJPXv2bPracc3RIKh6pWHnXCosgnKX0yBuwzxLSvkzC9KzVXfdZihvE6puql5pyG1hGuZQKpuNDLHVqBBTU1PD0hqYQdwG+U9mK5Q/s52O6pBEytuEqpuqVxpyW5iGOS7UZoVnq1OIlft5N3PrZ49oWCh/Zjsd1fHVlLcJVTdVrzRo/lzqA9UxvlpXWSty0UUXpT6tepgnICh/ZjudQXK01XHFQX7MKG8Tqm6qXmnQ/LnUB6r9qY1GI2uFnpw6dUryc1P+zHY6g2wPw+zeVd4mVN1UvdKQ2z0m1b7xLLtG8nrRomp30jgQOUqOqpuqVxpyuyZZfQnDvGXGVgx6ZtzKysowtAZmm8ekCxKg+p+Z8jah6qbqlQbNrbIPsjrGlOUtMwbtOlleXpbceFW9xgHlY7Wq24Sqm6pXGjT34/tA8VgJrA+pv1l32t69ezd9/SAXJW6Fah+0qtc4oJoj5W1C1U3VKw25LUyqN/M7duzYloOCph0hYdAuQtXjCape48BmORrmWXdbobxNqLqpeqUht115qheTdd77pJMsx33byi0rVL3Ggc1yNMyz7rZCeZtQdVP1SkNuS2y9Xs9aoSuqJxiArpuq1zgQOUqOqpuqVxpyW5hU+1NVgw66bqpe40DkKDmqbqpeachtYVI9xqTcz6vqpuo1DkSOkqPqpuqVhtyuSV6PMWWJqpuq1zgQOUqOqpuqVxpGXpjM7EYzu6P596Ye89xlZs+Z2U/M7HXd5lHdba3Valkr9ETVTdVLmchRdqi6qXqlYaRn5ZnZNHAIuAYoAkfMbL9vOPXHzC4Cjrv7JZstS7VvXPWCRdB1U/VSJXKULapuql5pGPUe03XAoq9TA2aBSzvmeSfwETP7rpm9ecR+A6PaZw+6bqpewkSOMkTVTdUrDUPbYzKzg8DlHU9/APjrDe0V4GLg2dYT7n4vcK+ZvRu438ze6O5nOpf/yCOPnPdFzM/Ps7CwQL1eb+/SlstlGo0GZ8+eBV7qg20NDz8zM0OhUKBarQLr9zOZnJxst4vFIsVisX3L4mKxSKlUYmlpCXdnamqKUqnE8vIyjUaDyclJpqamOHHiBI1Gg0KhwOzsLCsrK9Tr9Xa7VquxtraGmVGpVKjVau1bXFcqFVZXV9vtcrm8beu0urrKmTNnEq/T9PR0uz2MdWr5jfJ76nedsmbYObrnnnu4557ud3mfm5tjcXExckTkaOQ52myUgu3+A34JuH9D+/vA6zaZ/3PAVd2m3XHHHa7I6dOns1boiaqbqpe7O3DQR5iRfv4iR9mi6qbq5Z48R6Puyvs28CoAMysBa+5+zMymzGyu+fxGpyrwZLcFnTt3bsiq6Wj9mlFE1U3VS5jIUYaouql6pWGkJz+4e83MDprZx4AJ4JbmpPcCbwduBb5hZieBo8Dn3H3nHNELgm0gchTsdEY+Vp67HwYOdzz3VeCrzce/0s9yisXi9sttA5VKJWuFnqi6qXopEznKDlU3Va805PYC2+iCSI6qm6rXOBA5So6qm6pXGqIwbTPKG4eqm6rXOBA5So6qm6pXGnJbmIIgCIKdSW4Lk+oYX+VyOWuFnqi6qXqNA5Gj5Ki6qXqlIbeFSXUoFdWxx0DXTdVrHIgcJUfVTdUrDbktTKp94wqjBfRC1U3VaxyIHCVH1U3VKw25LUxBEATBziS3hSn6xpOj6qbqNQ5EjpKj6qbqlYbcFiZ/aYR/KVT77EHXTdVrHIgcJUfVTdUrDbktTKoH+loj+yqi6qbqNQ5EjpKj6qbqlYbcFqYgCIJgZ5LbwqTaN966B4oiqm6qXuNA5Cg5qm6qXmnIbWEKgiAIdia5LUyq97dv3S1SEVU3Va9xIHKUHFU3Va805LYwBUEQBDuT3BamycmR30qqL2ZmZrJW6Imqm6rXOBA5So6qm6pXGnJbmMwsa4WuFAq6H6mqm6rXOBA5So6qm6pXGka+Jmb2SjP7lJn9UY/pl5jZJ83sd8zs/b2Wo9o3Xq1Ws1boiaqbqpcykaPsUHVT9UrDyAuTux8H/hso9Zjls8CX3P0LwAfM7KKRyQVBTogcBTuZrPb9ut5q0cxKwI3ufqz51DPAO7rNOzExMSS1wSiVev0/kT2qbqpeOSBylAGqbqpeaRjakU8zOwhc3vm8u//GJi/7GeDFDe0V4OJuMz700EPn9Y/Pz8+zsLBAvV5vD/9eLpdpNBrtoTpaF6C1TqucmZmhUCi0d4FLpRKTk5PtdrFYpFgssrS01G6XSiWWlpZwd6ampiiVSiwvL9NoNNoHkk+cOEGj0aBQKDA7O8vKygr1er3drtVqrK2tYWZUKhVqtVr7tsiVSoXV1dV2u1wub9s61Wo13D3xOk1PT7fbw1iniYkJqtXqSL+nftcpayJHkaNxzNHQCpO7H0zxskXO75rYBfy424zXX389jz32WNeF7Nq167x256i7nVdId57N0tmenp4+r935y2Tj9MXFRfbt27fp65Msr8V2rNPi4iK7d+9O5ZC0nWSdFhcXKZfLI/2e0rSzIHLU+/WRo3V2Yo4kTuMwsykzm3P3VeBRM3tNc9IbgH/s9ppHH310VHqJOHToUNYKPVF1U/VqckPWAv0SORoNqm6qXk1uSDKzjXrYezObA+4GLgQ+4u7Pm9n7gLe7+61mdhnwYeAY8KK7f6XHclxxyH4zk72VgKqbqhe03eTOqY4cZYeqm6oXJM/RyAvTdhGBSo6qm6oX6Bam7SJylBxVN1UvSJ4jzcu+++MJM3sga4ku3GBmD2Ut0QNVN1UvgPdkLTBkIkfJUXVT9YKEOcrtHlMQBEGwM5E4+SEIgiAIWkRhCoIgCKSIwhQEQRBIkbvCZGYfNbNbzOzO5imzMpjZz5rZj8zsOTP7CwGfK8zs62a2v9nua2DPUXs1n7ur+bn9xMxel5HXbjP7OzN72sweMLOyyme23USOEvlEjpJ5DZyjXJ2VZ2ZXA9e6+wfN7G3AncDtGWtt5BbgSneXGObX3Z82sxeB1mmanwV+392Pmdnfm9nD7v7TrL2aA4wed/dLRu3SwduBDwFngK8AvwX8MgKf2XYSOUpG5CgxA+cob3tMB4Anm4+/12xLYGYzrF/dfMzM/tTMVIr+KiQb2HNEbByA9J3AR8zsu2b25qyE3P0f3P1088Kex4EX0PrMtosDRI6SEjnqk+3IUd4K04XAiebjngNTZoG7n3X3m4CfA34e+MOMlTrpe2DPUePu97r7FcDHgfvNbHfWTsBrgMcQ/cwGJHKUnshRMl5DihzlrTA9D7RGJNwFHM/QpSvufhr4XeCarF066Dawp9Tn5+6HgW8Ar83Sw8zeC3wa+D/EP7OURI7SEznqk0FylLfCdB9wVfPxlYDUFetm7fsHvJL1DUOGJAN7jhoz27gdVnmpmykLl3cA/+nuPwEuQPQzG5DIUUoiR327DJQjlf7bvnD3J8zsKTO7GbiU9V1WCczsl4FPmdl9wA/d/YsCTq9mfSN4q5kdBW5jvQ/6GPBld39OxOs+MzsJHAU+5+6Z3O/bzG5j/SSAF5oh/w4in9l2EjlK7BQ5SuZ1GwPmKIYkCoIgCKTIW1deEARBsMOJwhQEQRBIEYUpCIIgkCIKUxAEQSBFFKYgCIJAiihMQRAEgRRRmIIgCAIpojAFQRAEUkRhGgPM7Boz+7GZZTp2VhDkmcjR6IjCNAa4+3eA/3H3H2XtEgR5JXI0OqIwjQFmdhlwbMsZgyDoSeRodMRYeWOAmf0m63e5/CHwPtZ/9f1ZtlZBkC8iR6MjCtMYYGZfAn4MfB44BZDVyMNBkFciR6MjuvLGg8tYvzfLp919LcIUBKmIHI2IKEw7HDO7CDju7vcDrzWzkpm9O2uvIMgTkaPREoVp53M18M3m428BvwY8lJlNEOSTyNEIiWNMQRAEgRSxxxQEQRBIEYUpCIIgkCIKUxAEQSBFFKYgCIJAiihMQRAEgRRRmIIgCAIpojAFQRAEUkRhCoIgCKT4fwu1kmaebYkyAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x165.6 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1,2, figsize=(6, 2.3))\n",
    "\n",
    "Td = np.arange(0, 2.1, ts)\n",
    "y, t = initial(G, Td, [1, -0.5])\n",
    "ax[0].step(t/ts, y[:,0], where='post', ls='-', lw=1, color='k')\n",
    "\n",
    "ax[1].step(t/ts, y[:,1], where='post', ls='-', lw=1, color='k')\n",
    "\n",
    "plot_set(ax[0], '$k$', '$x_1$')\n",
    "ax[0].set_xlim([0, 2/ts])\n",
    "ax[0].set_ylim([-1.1, 1.1])\n",
    "plot_set(ax[1], '$k$', '$x_2$')\n",
    "ax[1].set_xlim([0, 2/ts])\n",
    "ax[1].set_ylim([-1.1, 1.1])\n",
    "\n",
    "fig.tight_layout()\n",
    "#fig.savefig(\"disc_pole_response.pdf\", transparent=True, bbox_inches=\"tight\", pad_inches=0.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 最適レギュレータ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 例14.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[0.44896593, 1.56029484]])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Q = np.matrix([[1,0],[0,1]])\n",
    "R = 0.1\n",
    "X, E, F = dare(P.A, P.B, Q, R)\n",
    "F"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.33661977+0.j, 0.1030854 +0.j])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[2.23644233, 1.42321107],\n",
       "        [1.42321107, 2.75297806]])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[-0.44896593, -1.56029484]])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "FF = -np.linalg.inv(R+P.B.T*X*P.B)*(P.B.T*X*P.A)\n",
    "FF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.33661977 0.1030854 ]\n"
     ]
    }
   ],
   "source": [
    "G = ss(P.A-P.B*F, P.B, P.C, P.D, ts)\n",
    "print(G.pole())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaYAAACcCAYAAAA05w98AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAXb0lEQVR4nO2df4wjZ3nHP4/Xa++u9469a0LzSyGBiKZFQTRCCQokl6ihUP5AQZWgRRVVUdIKiKqk0pZUpcqtSFB7h0RRhcoPIaWiohK0CQ3lpPZKlB+kydGWQEuitOFAQWnTXFZ3XOLdW+/6/PSPtZ09x97zjHc939l5PtLq/HrG48/Y873H887MO+buBEEQBIEKpawFgiAIgmAjUZiCIAgCKaIwBUEQBFJEYQqCIAikiMIUBEEQSBGFKQiCIJCinLVAWm6++Wa/6KKLstYIdjgLCwtfdvebs/bYLiJHwThImqPcFqa9e/eyf//+rDVexeLiIuecc07WGn1RdVP1AlhYWHgua4ftJHKUHFU3VS9InqPMuvLM7HIz+6aZ7esz7QYzu6P99+Ys/IIgD0SOgp1IZntM7v60mb0M2MbnzWwKOABcBVSAw2a2z3uGqJiYmBibaxIqlUrWCgNRdVP1ygORo/Gj6qbqlYasT35Y7fPcNcCir9MAZoCLe2eKQCVH1U3VK0dEjsaIqpuqVxoUjzGdBxzf0F4BLgCe3TjT3Xffzd133913AXNzczzzzDPUajVarRanTp0CYGZmBoDl5WUApqenKZVKLC0tAVCtVimXy912pVKhUqlQr9e77Wq1Sr1ex92ZnJykWq2yvLxMq9WiXC6zvLzM1NQUrVaLUqnEzMwMKysrNJvNbrvRaLC2toaZMTs7S6PRYHV1/f+W2dlZVldXu+1arUaz2aTRaHTbadfp5MmTnHvuuYnXaWpqqtvejnVqNBrMzc2N9Xsadp1yzFA5Onz48Bk5mp+fZ2FhYcu2OYgcRY6S50ixML3A+q+7DruAY70zzc/Pc+DAgb4LMLMzDgLWarUzpne+sA7T09Obtqemps5oV6vVgdObzSZ79+7d9PVJltdh165dZ7TTrtPu3btTOSRtJ1mnxcVFarXaWL+nNO2cMVSO9u3bx5EjR/ouYKu2uQ6Ro8jRsGTdldfFzObMbBL4DnBh+7kqsObuR3vnjy6I5Ki6qXrlkcjR9qPqpuqVhsz2mMzsdcAbgbeb2RHgLuAhd/+6me03s48DE8At/V6vGqjeXxdKqLqpeuWByNH4UXVT9UpDlmflPcv6AdoOt26Ydgg4tNnrO/2satTrddkNRNVN1SsPRI7Gj6qbqlcaZLrydgrKN15UdVP1CrJDeZtQdVP1SkNuC1OppKk+OTmZtcJAVN1UvYpA5Cg5qm6qXmnQ3CqHoFxWPKFQu59X1U3VqwhEjpKj6qbqlYbcFqa1tbWsFfrSuWZAEVU3Va8iEDlKjqqbqlcacluYVPtTW61W1goDUXVT9SoCkaPkqLqpeqUht4VJtW9ctWsEdN1UvYpA5Cg5qm6qXmnQ3CqHQPVLUB4tQNVN1asIRI6So+qm6pWG3Bam6BtPjqqbqlcRiBwlR9VN1SsNuS1M0TeeHFU3Va8iEDlKjqqbqlcacluYzOzsM2WAap896LqpehWByFFyVN1UvdKQ2zVRvZisdyRfJVTdVL2KQOQoOapuql5pyG1hajabWSv0ZWVlJWuFgai6qXoVgchRclTdVL3SkNvCpNqfqhp00HVT9SoCkaPkqLqpeqUht4Up+saTo+qm6lUEIkfJUXVT9UpDbtck+saTo+qm6lUEIkfJUXVT9UpDbguT6m5rmvvbjwtVN1WvIhA5So6qm6pXGnJbmFT7xlUvWARdN1WvIhA5So6qm6pXGnJbmFRR7bMHXTdVryA7lLcJVTdVrzRkNlCWmX0MWAUuAD7r7j/rmf7PwOXAc8Db3f30xumVSmVMpsmYnZ3NWmEgqm6qXnkgcjR+VN1UvdKQyR6TmV0JXO3uXwK+DdzZM/0m4I/d/SJ3f1tvmABOn37VUxIo9/Oquql6qRM5ygZVN1WvNGS1x3QT8GT78feBrwC3b5h+A/DrZvZd4CPu/kLvAh544IEzdl3n5+dZWFjoHsxdXFykVqvRarU4deoU8MpZK53BDqenpymVSiwtLQHrd4Asl8vddqVSoVKpUK/Xu+1qtUq9XsfdmZycpFqtsry8TKvVolwus7y8TLPZpNVqUSqVmJmZYWVlhWaz2W03Gg3W1tYwM2ZnZ2k0GqyurgLrv3pWV1e77VqtRrPZ7G50o6zTyZMnu+uRZJ2mpqa67e1Yp0ajwcTExFi/p2HXSZyb2MYcbcU2B5GjyFGKHLn72P+ALwK3tB+XgUafeSaAPwH+qd8y5ufnfRDrq5UNL774YmbvfTZU3VS93N2B/Z5BRob52+4cZYnyNqHqpurlnjxHWZ388ALQOel+F3CsdwZ3P+3unwR291tA9I0nR9VN1SsHRI4yQNVN1SsNWRWme4G3tB9fAdxvZnNmNglgZqX2v1PAw/0WoNo33tnNVkTVTdUrB0SOMkDVTdUrDZkcY3L3J8zsKTO7GbgY+ARwF/CQmT0APGJmh4H/Aj7ZbxkRqOSouql6qRM5ygZVN1WvNGR2uri7H+x56tYNj39plGXv2bNn03P69+zZw/Hjx0d5iyCQYDtzFARZkVlhGpXNxvg6W9HZzgvRarXati17VFTdVL2KgOpYecrbhKqbqlcacjvyg+pQKqpjj4Gum6pXEYgcJUfVTdUrDbktTKp948rXvqi6qXoVgchRclTdVL3SkNvCFARBEOxMcluYom88Oapuql5FIHKUHFU3Va805LYwrV9MrIdqnz3ouql6FYHIUXJU3VS90pDbwqR6oK8zTpUiqm6qXkUgcpQcVTdVrzTktjAFQRAEO5PcFibVvvHOiL6KqLqpehWByFFyVN1UvdKQ28IUBEEQ7ExyW5hU72/fufeJIqpuql5FIHKUHFU3Va805LYwBUEQBDuT3Bamcjn9MH+dQV4H/e3duzf1sqenp1O/drtRdVP1KgKj5Gg7Ud4mVN1UvdKguVUOwSgDsW7nIK+lkm6tV3VT9SoC2zmg8SgobxOqbqpeacjtmqj2jS8tLWWtMBBVN1WvIhA5So6qm6pXGnJbmIIgCIKdSarCZGbnbLVIUiYmJrJW6Eu1Ws1aYSCqbqpeW4mZzZrZ+83smnb7UjO7PluryFEaVN1UvdJw1mNMZvZuXl3A3g38/ihvbGYfA1aBC4DPuvvPNky7Abi63Tzk7v/R+/rt7E/d7A64Z7v7rerBZNB1U/XaYu4AjgBvMrO3AX8OfAF4cJSFKudoFJS3CVU3Va80DLMmvwz8b89zrxnlTc3sSuBqd/+Qmb0DuBO4vT1tCjgAXAVUgMNmts97Rpvczr7xzQrP2Q4WLy0tyZ4do+qm6rXF/DvwQ3f/pplNA+8FRlpp9RyNgvI2oeqm6pWGYQrTPe7+/MYnzOzwiO97E/Bk+/H3ga/QDhRwDbDYDlDDzGaAi4FnR3zPLWGzvalhX3+2swKDHcnTwG8Cn3L3U8A3zGxuxGXeRE5zFASbMUxhut3MvuHu/9L+VXbE3Xv3oJJyHvCv7ccrrHdDbJy28X/uzvQzAvXwww+fUSDm5+dZWFig2Wx27+RYq9VotVrdUXc7Y0l1rpCenp6mVCp1z2apVquUy+Vuu1KpUKlUqNfr3fbzzz9PvV7H3ZmcnKRarbK8vEyr1aJcLtNsNimVSrRaLUqlEjMzM6ysrHSfv+yyy2RP0Q22lQ8D9wGY2bXA4+5+z4jLzG2OLrnkEk6cODHi6gc7lWEK01Os/9qD9T7yDwB/PeL7vgB0RhzcBRwbMK3fdABuvPFGHn/88b4L37Vr1xnt3hto9Q522Lv729uempo6o917kHHj9JWVlVfNv7Gd5d5SPzcFVL1gS6/zeZJXcvQ4BcjRhRdeOLD47NmzZ9N7QSlvE6puql6QPEfDHPl8HfBXZnYb68eb5hJbvZp7gbe0H18B3G9mc2Y2CXwHuBDAzKrAmrsf7V3A6urqFmhsPZ1fhYqouql6bTGFy9GJEydw975/Z/uBprxNqLqpeqVhmD2mHwF/CVzHenfEM6O+qbs/YWZPmdnNrPd7fwK4C3jI3b9uZvvN7OPABHDLqO8XBAJEjoJgSIYpTKeA0+7+t2b2b8B7tuKN3f1gz1O3bph2CDi02etVr7+oVCpZKwxE1U3Va4vZcTnau3fvpseJ9uzZk3rZytuEqpuqVxrO2pXn7vcC57ebrwFqm8w+NlQLk/JFbqpuql5byU7M0WZddcN0122G8jah6qbqlYahrshy9x+2//0B8INtNRoS5WNMqhuIqpuq11YTORoe5W1C1U3VKw2al33nmM3ONMoaVTdVryA7lLcJVTdVrzTktjCpDqUyOTmZtcJAVN1UvYpA5Cg5qm6qXmnQ3CqHQHVcKOVdaVU3Va8iEDlKjqqbqlcacluYVMf46lwNr4iqm6pXEYgcJUfVTdUrDbktTKr9qa1WK2uFgai6qXoVgc1ytHfvXsxs4N8op4OfDeVtQtVN1SsNmvvxQ6DaN67aNQK6bqpeRWCzHHVOB88C5W1C1U3VKw2a/7sPgeqXoDpWFei6qXoVgchRclTdVL3SkNvCFH3jyVF1U/UqApGj5Ki6qXqlIbeFKY4xJUfVTdWrCESOkqPqpuqVhtwWJtV7Gqke+wJdN1WvIhA5So6qm6pXGnK7JqoXk/Xeo0YJVTdVryIQOUqOqpuqVxpyW5iazWbWCn1ZWVnJWmEgqm6qXkUgcpQcVTdVrzTktjCp9qeqBh103VS9ikDkKDmqbqpeachtYYq+8eSouql6FYHIUXJU3VS90pDbNYm+8eSouql6FYHIUXJU3VS90jD2wmRmN5jZHe2/Nw+Y5y4ze87Mfmpmb+g3j+pua6PRyFphIKpuql7KRI6yQ9VN1SsNY73s28ymgAPAVUAFOGxm+3zDxRRmdj5wzN0v2mxZqn3jqhcsgq6bqpcqkaNsUXVT9UrDuPeYrgEWfZ0GMANc3DPPu4CPmtn3zOytY/YbGdU+e9B1U/USJnKUIapuql5p2LY9JjPbD1zW8/QHgb/Z0F4BLgCe7Tzh7vcA95jZe4D7zOxN7v5S7/IfeeSRM76I+fl5FhYWaDab3V3aWq1Gq9Xi1KlTwCt9sJ2hO6anpymVSiwtLQHr9zMpl8vddqVSoVKpUK/Xu+1qtUq9XsfdmZycpFqtsry8TKvVolwuMzk5yfHjx2m1WpRKJWZmZlhZWaHZbHbbjUaDtbU1zIzZ2VkajUb3Ftezs7Osrq5227VabcvWaXV1lZdeeinxOk1NTXXb27FOHb9xfk/DrlPWbHeODh48yMGDB/u+99zcHIuLi5EjIkdjz5G7j+0P+BXgvg3tHwBv2GT+zwNv6TftjjvucEVOnjyZtcJAVN1Uvdzdgf0+xowM8xc5yhZVN1Uv9+Q5GndX3neACwHMrAqsuftRM5s0s7n28xudloAn+y3o9OnT26yajs6vGUVU3VS9hIkcZYiqm6pXGsZ68oO7N8xsv5l9HJgAbmlPeh9wHXAr8C0zOwEcAT7v7jvniF4QbAGRo2CnM/absbj7IeBQz3NfA77WfvxrwyynUqlsvdwWMDs7m7XCQFTdVL2UiRxlh6qbqlcacnuBbXRBJEfVTdWrCESOkqPqpuqVhihMW4zyxqHqpupVBCJHyVF1U/VKQ24LUxAEQbAzyW1hUh3jq1arZa0wEFU3Va8iEDlKjqqbqlcacluYVIdSUR17DHTdVL2KQOQoOapuql5pyG1hUu0bVxgtYBCqbqpeRSBylBxVN1WvNOS2MAVBEAQ7k9wWpugbT46qm6pXEYgcJUfVTdUrDbktTP7KCP9SqPbZg66bqlcRiBwlR9VN1SsNuS1Mqgf6OiP7KqLqpupVBCJHyVF1U/VKQ24LUxAEQbAzyW1hUu0b79wDRRFVN1WvIhA5So6qm6pXGnJbmIIgCIKdSW4Lk+r97Tt3i1RE1U3VqwhEjpKj6qbqlYbcFqYgCIJgZ5LbwlQuj/1WUkMxPT2dtcJAVN1UvYpA5Cg5qm6qXmnIbWEys6wV+lIq6X6kqm6qXkUgcpQcVTdVrzSMfU3M7LVm9mkz+6MB0y8ys0+Z2e+a2QcGLUe1b3xpaSlrhYGouql6KRM5yg5VN1WvNIy9MLn7MeC/geqAWT4HfNndvwh80MzOH5tcEOSEyFGwk8lq36/vrRbNrArc4O5H2089A/xqv3knJia2SW00qtVB/09kj6qbqlcOiBxlgKqbqlcatu3Ip5ntBy7rfd7df2uTl/0c8PKG9gpwQb8ZH3zwwTP6x+fn51lYWKDZbHaHf6/VarRare5QHZ0L0DqnVU5PT1Mqlbq7wNVqlXK53G1XKhUqlQr1er3brlar1Ot13J3JyUmq1SrLy8u0Wq3ugeTjx4/TarUolUrMzMywsrJCs9nsthuNBmtra5gZs7OzNBqN7m2RZ2dnWV1d7bZrtdqWrVOj0cDdE6/T1NRUt70d6zQxMcHS0tJYv6dh1ylrIkeRoyLmaNsKk7vvT/GyRc7smtgF/KTfjNdeey2PPfZY34Xs2rXrjHbvqLu9V0j3ns3S256amjqj3fvLZOP0xcVFzjnnnE1fn2R5HbZinRYXF9m9e3cqh6TtJOu0uLhIrVYb6/eUpp0FkaPBr48crbMTcyRxGoeZTZrZnLuvAo+a2SXtSW8E/rHfax599NFx6SXiwIEDWSsMRNVN1avN9VkLDEvkaDyouql6tbk+ycw27mHvzWwOuBs4D/iou79gZu8HrnP3W83sUuAjwFHgZXf/6oDluOKQ/WYmeysBVTdVL+i6yZ1THTnKDlU3VS9InqOxF6atIgKVHFU3VS/QLUxbReQoOapuql6QPEeal30PxxNmdn/WEn243swezFpiAKpuql4A781aYJuJHCVH1U3VCxLmKLd7TEEQBMHOROLkhyAIgiDoEIUpCIIgkCIKUxAEQSBF7gqTmX3MzG4xszvbp8zKYGY/b2Y/NrPnzOwvBHwuN7Nvmtm+dnuogT3H7dV+7q725/ZTM3tDRl67zezvzOxpM7vfzGoqn9lWEzlK5BM5SuY1co5ydVaemV0JXO3uHzKzdwB3ArdnrLWRW4Ar3F1imF93f9rMXgY6p2l+DvgDdz9qZn9vZg+7+/NZe7UHGD3m7heN26WH64APAy8BXwV+B3gnAp/ZVhI5SkbkKDEj5yhve0w3AU+2H3+/3ZbAzKZZv7r5qJn9qZmpFP1VSDaw55jYOADpu4CPmtn3zOytWQm5+z+4+8n2hT2PAy+i9ZltFTcROUpK5GhItiJHeStM5wHH248HDkyZBe5+yt1vBH4B+EXgDzNW6mXogT3Hjbvf4+6XA58A7jOz3Vk7AZcAjyH6mY1I5Cg9kaNkXEKKHOWtML0AdEYk3AUcy9ClL+5+Evg94KqsXXroN7Cn1Ofn7oeAbwGvz9LDzN4HfAb4P8Q/s5REjtITORqSUXKUt8J0L/CW9uMrAKkr1s269w94LesbhgxJBvYcN2a2cTtc4pVupixcfhX4T3f/KXAuop/ZiESOUhI5GtplpByp9N8Ohbs/YWZPmdnNwMWs77JKYGbvBD5tZvcCP3L3Lwk4vY71jeDtZnYEuI31PuijwFfc/TkRr3vN7ARwBPi8u2dyv28zu431kwBebIf8u4h8ZltJ5CixU+QomddtjJijGJIoCIIgkCJvXXlBEATBDicKUxAEQSBFFKYgCIJAiihMQRAEgRRRmIIgCAIpojAFQRAEUkRhCoIgCKSIwhQEQRBIEYWpAJjZVWb2EzPLdOysIMgzkaPxEYWpALj7d4H/cfcfZ+0SBHklcjQ+ojAVADO7FDh61hmDIBhI5Gh8xFh5BcDMfpv1u1z+CHg/67/6/ixbqyDIF5Gj8RGFqQCY2ZeBnwBfAH4GkNXIw0GQVyJH4yO68orBpazfm+Uz7r4WYQqCVESOxkQUph2OmZ0PHHP3+4DXm1nVzN6TtVcQ5InI0XiJwrTzuRL4dvvxQ8BvAA9mZhME+SRyNEbiGFMQBEEgRewxBUEQBFJEYQqCIAikiMIUBEEQSBGFKQiCIJAiClMQBEEgRRSmIAiCQIooTEEQBIEUUZiCIAgCKf4f5kq9d+0Hw2EAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x165.6 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1,2, figsize=(6, 2.3))\n",
    "\n",
    "Td = np.arange(0, 2.1, ts)\n",
    "y, t = initial(G, Td, [1, -0.5])\n",
    "ax[0].step(t/ts, y[:,0], where='post', ls='-', lw=1, color='k')\n",
    "\n",
    "ax[1].step(t/ts, y[:,1], where='post', ls='-', lw=1, color='k')\n",
    "\n",
    "plot_set(ax[0], '$k$', '$x_1$')\n",
    "ax[0].set_xlim([0, 2/ts])\n",
    "ax[0].set_ylim([-1.1, 1.1])\n",
    "plot_set(ax[1], '$k$', '$x_2$')\n",
    "ax[1].set_xlim([0, 2/ts])\n",
    "ax[1].set_ylim([-1.1, 1.1])\n",
    "\n",
    "fig.tight_layout()\n",
    "#fig.savefig(\"disc_lqr_response.pdf\", transparent=True, bbox_inches=\"tight\", pad_inches=0.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 章末問題"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 章末問題7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A = [[0.5 1. ]\n",
      "     [0.  1. ]]\n",
      "\n",
      "B = [[0.]\n",
      "     [1.]]\n",
      "\n",
      "C = [[1. 0.]\n",
      "     [0. 1.]]\n",
      "\n",
      "D = [[0.]\n",
      "     [0.]]\n",
      "\n",
      "dt = 0.1\n",
      "\n"
     ]
    }
   ],
   "source": [
    "ts = 0.1\n",
    "A = [[0.5, 1], [0, 1]]\n",
    "B = [[0],[1]]\n",
    "C = np.eye(2)\n",
    "D = np.zeros([2,1])\n",
    "P = ss(A,B,C,D, ts)\n",
    "print(P)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.5, 1. ])"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P.pole()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[0., 1.],\n",
       "        [1., 1.]])"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Vc = ctrb(P.A, P.B)\n",
    "Vc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[-0.25, -1.5 ]])"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "regulator_pole = [0,0]\n",
    "F = -acker(P.A, P.B, regulator_pole)\n",
    "F"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 章末問題8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A = [[ 1.  1.]\n",
      "     [ 0. -1.]]\n",
      "\n",
      "B = [[0.]\n",
      "     [1.]]\n",
      "\n",
      "C = [[1. 0.]\n",
      "     [0. 1.]]\n",
      "\n",
      "D = [[0.]\n",
      "     [0.]]\n",
      "\n",
      "dt = 0.1\n",
      "\n"
     ]
    }
   ],
   "source": [
    "ts = 0.1\n",
    "A = [[1, 1], [0, -1]]\n",
    "B = [[0],[1]]\n",
    "C = np.eye(2)\n",
    "D = np.zeros([2,1])\n",
    "P = ss(A,B,C,D, ts)\n",
    "print(P)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "Q = np.matrix([[0.5,0],[0,0]])\n",
    "R = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "X, E, F = dare(P.A, P.B, Q, R)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[2., 1.],\n",
       "        [1., 1.]])"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.70710678+0.j, -0.70710678+0.j])"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.70710678 -0.70710678]\n"
     ]
    }
   ],
   "source": [
    "G = ss(P.A-P.B*F, P.B, P.C, P.D, ts)\n",
    "print(G.pole())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaYAAACcCAYAAAA05w98AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAWW0lEQVR4nO3dfWxkV3nH8e8zHs8de8bBdrdVNlmFBAhCpaA0qgINhCQqpRRVaFElKKhqpSpUAqKKSHWIGiTWIilityrlj6q8CDVtEQioEgolUrUQ5YU0CbSFtATBLumyEW26W3ffMvZ6Xnae/jEva49nd+fG47nnPPf5SFbmzlxfn198nj2+Z+49I6qKc845F4pC1g1wzjnnNvKByTnnXFB8YHLOORcUH5icc84FxQcm55xzQfGByTnnXFCKWTfgxbrtttt0z549WTfDGbe8vPw5Vb0t63bsFK8jNwlp6yjagWlxcZF9+/Zl3YyxWFlZYdeuXVk3YywsZQFYXl7+WdZt2EleR+GylCdtHWU2lScirxKRr4vIzUNeu1VE7up+vTaL9jkXA68jZ1FmZ0yq+iMReQGQjc+LSBnYD9wAlICDInKzDixRMTU1NbG27rRSqZR1E8bGUpYYeB2dZ63vWcuTRtYXPzSGPHcjsKIddWAWuGpwJy+oMFnKEhGvI+z1PWt50gjxPabLgRMbtteBK4CjG3e69957uffeey95sIWFBQ4dOgTA7OwsAGtrawDMzMxQKBRYXV0FIEkSisVif7tUKlEqlajVav3tJEmo1WqoKtPT0yRJwtraGu12m2KxSLlc7m8XCgVmZ2dZX1+n1Wr1t+v1Os1mExGhWq1y7NgxKpUKANVqlUajQaPR+bemUqnQarWo1+v97Xa7zdmzZ4PMdPr0afbs2UO9Xu9niDlTxEaqo4MHD26qo6WlJZaXl6P5/VitIxGhXq/3c8SeKa0QB6ZjdP6665kDjg/utLS0xP79+y95MBHZ8gZi7xfWMzMzc9Htcrm8aTtJkou+nna7UqlsauPg6wBzc3NbvmejkDIlSXLJ40FcmSI0Uh3dfPPNPPXUU0MPENvvx1odrayscNlll110f4gr06iynsrrE5F5EZkGvg1c2X0uAZqq+uzg/j4FESZLWWLkdWSHtTxpZHbGJCIvBV4JvEFEngLuAR5R1a+IyD4R+RAwBbx32PdbKqjBv0hiZilLDLyOzrPW96zlSSPLq/KO0nmDtuf2Da89CDx4se/vzbNaUKvVzHRCS1li4HV0nrW+Zy1PGsFM5eWZpQ9rtJTFxcVa37OWJ41oB6ZCIdqmbzE9PZ11E8bGUpY88DoKl7U8aUTbK4vFEC8ofHEsna5bypIHXkfhspYnjWgHpmazmXUTxqZ3n4EFlrLkgddRuKzlSSPagcnS/Gu73c66CWNjKUseeB2Fy1qeNKIdmCzNjVuaTrGUJQ+8jsJlLU8a0fZKS7+0yFcY2MRSljzwOgqXtTxpRDsw+dx4mCxlyQOvo3BZy5NGtAOTz42HyVKWPPA6Cpe1PGlEOzCJyKV3ioSleX5LWfLA6yhc1vKkEW1ySzefDa7+GzNLWfLA6yhc1vKkEe3A1Gq1sm7C2Kyvr2fdhLGxlCUPvI7CZS1PGtEOTJbmXy3942ApSx54HYXLWp40oh2YfG48TJay5IHXUbis5Ukj2uQ+Nx4mS1nywOsoXNbypBHtwGTpNLder2fdhLGxlCUPvI7CZS1PGtEOTJbmxi3d5GgpSx54HYXLWp40oh2YLLE0z28pi4uLtb5nLU8amS2UJSIfABrAFcAnVfXUwOvfBF4F/Ax4g6qe2/h6qVSaUEt3XrVazboJY2MpSwy8js6z1ves5UkjkzMmEbkeeJ2qfhb4FvCRgdf3Aner6h5Vff1gMQGcO7flqWhZmku2lCV0XkebWet71vKkkdUZ017gme7j7wN/B9yx4fVbgd8Wke8A71PVY4MHeOihhzad6i4tLbG8vEyr1er/QiuVCgsLCyOdEs/Pz/Pcc89RLBZZXV0FOn9NlkolarVafztJEmq1GqrK9PQ0SZKwtrZGu92mWCxSLpf724VCgdnZWdbX12m1Wv3ter1Os9lERKhWq5w8eZJGowF0/kpqNBr97UqlsiVTu93m7NmzwPkrd3oLPs7MzFAoFPoZkiSZaKbTp0+TJAn1et1EpsDtZUJ1FOrvx2odiUi/rRYypSVZLOIoIp8BvquqnxWRIrCqqsnAPlPAnwA3qepbBo9x55136v79+8fZpswWtFxZWWHXrl2Z/Oxxs5QFQESWVXVf1u0YJsQ6ypK1vmcpT9o6yurih2NA7yL9OeD44A6qek5VPwpcNuwAPjceJktZIuB1tIG1vmctTxpZDUz3A9d1H78G+JqIzIvINICIFLr/LQOPDjuApbnx3qm5BZayRMDraANrfc9anjQyeY9JVb8nIj8UkduAq4APA/cAj4jIQ8BjInIQ+DHw0WHH8IIKk6UsofM62sxa37OWJ43MLhdX1QMDT92+4fEvTrItzsXK68hZFO0NtpbW+KpUKlk3YWwsZckDr6NwWcuTRrQDk6WlVCytV2YpSx54HYXLWp40oh2YLM2NR3C/zMgsZckDr6NwWcuTRrQDk3POOZuiHZh8bjxMlrLkgddRuKzlSSPagSmrVRp2gqV5fktZ8sDrKFzW8qQR7cBk6Y3B3tpWFljKkgdeR+GylieNaAemcest9nqpr8XFxayb6pxzpmV2g+12jXtu/MSJEyPttxMf3tVbBdgCS1nywNJ7TNb6nrU8afgZk3POuaBEOzA1m82smzA2vc9LscBSljzwOgqXtTxpRDswOeecsynagalYjPbtsS1mZmaybsLYWMqSB15H4bKWJ41oB6aduAhhFKNevZfmCr5CIdpfwxaWsuRBVnW0E6z1PWt50og2eVZz4ydOnEBVR/o6efLkSMdcXV3d4VZPjqUseWDpPSZrfc9anjSiHZicc87Z9KIGJhHZNe6GpDU1NZV1Ey5p1Gm/a6+9Nuumjk2SJFk3IUgiUhWRd4rIjd3ta0TklmxbFUcdjcpa37OWJ41LvvMpIm9l6wD2VuCPtvODReQDQAO4Avikqp7a8NqtwOu6mw+q6r8Pfn8M869Z3rSbFUtvpo/ZXcBTwKtF5PXAXwCfBh7ezkHzUEejstb3rOVJY5Tkvwz898BzL9nODxWR64HXqervicgbgY8Ad3RfKwP7gRuAEnBQRG7WgdUmLc2Nz8/PjzQ4LSwsjDzYZWV1dTXXVxNdxL8CP1DVr4vIDPB2YFv/o7yONrPW96zlSWOUgek+VX1+4xMicnCbP3cv8Ez38feBv6NbUMCNwEq3gOoiMgtcBRzd5s8M1uHDh9m169Kzo4uLi5lejRj6oBi4HwHvBv5UVc8CXxWR+W0ecy9eR86gUQamO0Tkq6r6z92/yp5S1cEzqLQuB77bfbxOZxpi42sb/wXsvb6poB599NFN/0gvLS2xvLxMq9Xqf/JjpVKh3W73V+ntrT3Vu6N6ZmaGQqHQv/olSRKKxWJ/u1QqUSqVqNVq/e0kSajVaqgq09PTJEnC2toa7XabYrFIuVzubxcKBWZnZ1lfX6fVavW36/U6zWYTEaFardJoNFhZWQHobzcajX6GXqZDhw5llmn37t2mphwz8AfAAwAichPwpKret81jbruODhw4wIEDBzYddGFhgaNHj26qoyuvvHLLVabz8/McPnx4U5+79tprOXXq1Jb9nn/++U197uqrrx561er8/DxHjhzZVEfDjrmwsMCRI0c21dGw/VycRhmYfkjnrz3ozJG/C/j8Nn/uMaC3QuEccPwCrw17HYA3v/nNPPnkk0MPPjc3t2l78AO3BhdHHDxdHtwul8ubtgfflBx8Pe324uLipucGX4fsM4166fv6+vrQ9sdqjIPxM5yvoycJpI7uvvtu7rnnnk3PiQhzc3Ob+tzJkye3fHaTiGw605+ZmeHUqVND9yuXy5v6xbDj9fadn58HzvfBYcdcXFzs79ezsLBg6vOlLNVS2joa5Z3PlwJ/IyIfpPN+03zqVm11P3Bd9/FrgK+JyLyITAPfBq4EEJEEaKrqs4MH6J1RWND7S9ICS1nGzOtojIbdT3jo0KGsmzVWea6lUc6YfgL8FfAmOtMRh7f7Q1X1eyLyQxG5jc6894eBe4BHVPUrIrJPRD4ETAHv3e7Pcy4AXkcj6N1iMficy5dRBqazwDlV/XsR+RfgbeP4wap6YOCp2ze89iDw4MW+39L9F6VSKesmjI2lLGPmdTSC7VxgY63vWcuTxiWn8lT1fmB3d/MlQOUiu09MaAW1HZZupLOUZZy8jnaetb5nLU8aI93Bpao/6P73aeDpHW3RiGKdGx+mVquZ6YSWsoxbXutoUtNz1vqetTxp5PfW4oBYupLIUhY3HpO6/81a37OWJ41o1yOxtJTK9PR01k0YG0tZ8sDrKFzW8qQRba+0tI6UpdN1S1nywOsoXNbypBHtwGRpja/eqg0WWMqSB15H4bKWJ41oByZL86/tdjvrJoyNpSx54HUULmt50oj2PN7S3Lil6RRLWfJgO3UU2s2w1vqetTxpRJvc0i/NynpYYCtLHmynjkJbbd5a37OWJ41oTzt8bjxMlrLkgddRuKzlSSPagcnnxsNkKUseeB2Fy1qeNKIdmCx9NpCl98ssZckDr6NwWcuTRrTJLd18NvhZSjGzlCUPvI7CZS1PGtEOTK1WK+smjM36+nrWTRgbS1nywOsoXNbypBHtwGRp/tXSPw6WsuSB11G4rOVJI9qByefGw2QpSx54HYXLWp40ok3uc+NhspQlD7yOwmUtTxoTH5hE5FYRuav79doL7HOPiPxMRJ4TkZcP28fSaW69Xs+6CWNjKUvIvI62stb3rOVJY6LLJ4hIGdgP3ACUgIMicrNuuJlCRHYDx1V1z8WOZWlu3NJNjpayhMrraDhrfc9anjQmfcZ0I7CiHXVgFrhqYJ/fAN4vIv8mIr8y4fZlwtI8v6UsAfM6GsJa37OWJ40dO2MSkX3AKwaefg/wxQ3b68AVwNHeE6p6H3CfiLwNeEBEXq2qZwaP/9hjj236xS0tLbG8vEyr1eqfAlcqFdrtNmfPngXOz9n2lvqYmZmhUCiwuroKdD7/pFgs9rdLpRKlUolardbfTpKEWq2GqjI9PU2SJKytrdFutykWi5TL5f52oVBgdnaW9fV1Wq1Wf7ter9NsNhERqtUqqsrKygoA1WqVRqPR/8jrSqUSVaZms0m9Xqder/czxJwpa1nU0YUWZ+310ZB+P1brSESYnp7mzJkz+awjVZ3YF/BrwAMbtp8GXn6R/T8FXDfstbvuukutOH36dNZNGBtLWVRVgX06wRoZ5cvraDhrfc9SnrR1NOmpvG8DVwKISAI0VfVZEZkWkfnu8xvbtAo8M+xA586d2+GmTk7vLyALLGUJmNfRENb6nrU8aUz04gdVrYvIPhH5EDAFvLf70juANwG3A98QkZPAU8CnVDW/7wA6N4TXkbNu4h9qpKoPAg8OPPdl4Mvdx785ynFKpdL4G5eRarWadRPGxlKWkHkdbWWt71nLk0a0N9j6FESYLGXJA6+jcFnLk4YPTAGw1AEtZckDr6NwWcuTRrQDk3POOZuiHZgsrfFVqVSybsLYWMqSB15H4bKWJ41oByZLS6lYWq/MUpY88DoKl7U8aUQ7MFmaGw9hhYFxsZQlD7yOwmUtTxrRDkzOOedsinZg8rnxMFnKkgdeR+GylieNaAcmPb/Cf/QszfNbypIHXkfhspYnjWgHJktvDPZWA7bAUpY88DoKl7U8aUQ7MDnnnLMp2oHJ0tx473NTLLCUJQ+8jsJlLU8a0Q5MzjnnbIp2YGo27azi3/uESQssZckDr6NwWcuTRrQDk3POOZuiHZiKxYl/lNSOmZmZyboJY2MpSx54HYXLWp40oh2YRCTrJoxNoRDtr2ELS1nywOsoXNbypDHxP5dE5BeAO4H/U9WPDXl9D/B+4KfAaVX90rDjWJobX11dNfPXkaUsIfM62spa37OWJ42JD8mqehw4BCQX2OUvgc+p6meA94jI7ok1zrlIeB05y7I6Vxz60YwikgC3quqz3acOA28Ztu/U1NQONW3ykuRC/7bEx1KWCHgdbWCt71nLk8aOTeWJyD7gFYPPq+rvXuTbfg54YcP2OnDFsB0ffvjhTfPjS0tLLC8v02q1+svFVyoV2u12f2mP3g1rvcswZ2ZmKBQKrK6uAp2OUCwW+9ulUolSqUStVutvJ0lCrVZDVZmeniZJEtbW1mi32xSLRcrlcn+7UCgwOzvL+vo6rVarv12v12k2m4gI1WqVs2fP9ttcrVZpNBr9j1WuVCpRZWo0GpRKJer1ej9DzJmy5nWUzzoSEaampjhz5kwu62jHBiZV3fcivm2FzVMTc8CRYTvedNNNPPHEE0MPMjc3t2l7cJXewTuqB+dxB7fL5fKm7cG/ZAZfT7tdKBTYtWvXBV+HeDKtrKyQJMkljwfxZMqS19Ho25bqCDq1tDHPsP0hrkyjCuKyDxGZFpF5VW0Aj4vI1d2XXgn807DvefzxxyfVvB23f//+rJswNpaydN2SdQNG5XVkq+8Zy3NLmp1l0svei8g8cC9wOfB+VT0mIu8E3qSqt4vINcD7gGeBF1T1Cxc4jlpZsl9EzHz8gKUs0M8T3DXVXkdbGe17WTdjLNLW0cQHpnHxggqTpSwQ7sA0Ll5H4bKUJ20dxXzb9/dE5GtZN2JMbhGRh7NuxJhYygLw9qwbsMO8jsJlKU+qOor2jMk555xNQVz84JxzzvX4wOSccy4oPjA555wLSnQXP4jIB+gsxXIF8ElVPZVti9ITkVcBB4A/U9VHRl1wMzQichnw18Cr6azb9m5ggQizAIhIEfgY8Kt0Vk54B7CLSPNcjNdROLyOhlDVaL6A64G/7T5+I/CJrNu0jSxfAG7pPv4H4OUbHu/Oun0jZvgt4CWAAF8Ebo81S7e9VwNz3cffBH4p5jwXyel1FNCX19HWr9im8vYCz3Qff7+7HasGpFtwMzSq+o+qelo7ve1J4H+JNAuAqv5UVV8QkVngMTrtjzbPRezF6ygYXkdbxTYwXQ6c6D6+4MKUkRl5wc3AXQ08QeRZutMqy8AfAq8n8jwX4HUUrqvxOoruPaZjQG9FwjngeIZtGZeRF9wMlYi8A/gE8D9EnkVVzwBLIvI08C4iz3MBXkcB8jo6L7YzpvuB67qPXwNEf8e6plhwM0Qi8hbgP1T1OeDniTjLgB/Tme6ykmcjr6PAeB1tFt3KDyKyBJwErqLzpu3JjJuUmoi8lM6bnN8A/pzO1MolF9wMjYh8ELiDzpx4AfgO8HEizAIgIr8O/DHwJUCBzwN7iDTPxXgdhcPraMgxYhuYnHPO2RbbVJ5zzjnjfGByzjkXFB+YnHPOBcUHJuecc0Hxgck551xQfGByzjkXFB+YnHPOBcUHJuecc0HxgSkHROQGETkiIi/Lui3OxcrraHJ8YMoBVf0O8F+q+p9Zt8W5WHkdTY4PTDkgItfQWaPKOfcieR1Njq+VlwMi8vt0Ph3zJ8A76fzV9/FsW+VcXLyOJscHphwQkc/R+fyTTwOnAFS1mWWbnIuN19Hk+FRePlxD5zNRPqGqTS8m514Ur6MJ8YHJOBHZDRxX1QeAl4lIIiJvy7pdzsXE62iyfGCy73rgW93HjwC/AzycWWuci5PX0QT5e0zOOeeC4mdMzjnnguIDk3POuaD4wOSccy4oPjA555wLig9MzjnnguIDk3POuaD4wOSccy4oPjA555wLyv8DNKA1eZiHetMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x165.6 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1,2, figsize=(6, 2.3))\n",
    "\n",
    "Td = np.arange(0, 3.1, ts)\n",
    "y, t = initial(G, Td, [1, 0])\n",
    "ax[0].step(t/ts, y[:,0], where='post', ls='-', lw=1, color='k')\n",
    "\n",
    "ax[1].step(t/ts, y[:,1], where='post', ls='-', lw=1, color='k')\n",
    "\n",
    "plot_set(ax[0], '$k$', '$x_1$')\n",
    "ax[0].set_xlim([0, 3/ts])\n",
    "ax[0].set_ylim([-1.1, 1.1])\n",
    "plot_set(ax[1], '$k$', '$x_2$')\n",
    "ax[1].set_xlim([0, 3/ts])\n",
    "ax[1].set_ylim([-1.1, 1.1])\n",
    "\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 付録"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = [[1,1,0],[1,0,1],[2,1,2]]\n",
    "B = [[1],[0],[1]]\n",
    "C = np.eye(3)\n",
    "D = np.zeros([3,1])\n",
    "\n",
    "ts = 0.1\n",
    "\n",
    "P = ss(A,B,C,D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "Q = np.eye(3)\n",
    "R = 100\n",
    "X, E, F = dare(P.A, P.B, Q, R)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[269.91859465, 153.35336484, 177.00773162],\n",
       "        [153.35336484,  89.13772697, 100.33774379],\n",
       "        [177.00773162, 100.33774379, 119.07469571]])"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-0.52918501  0.64442964  0.34784429]\n"
     ]
    }
   ],
   "source": [
    "G = ss(P.A-P.B*F, P.B, P.C, P.D, ts)\n",
    "print(G.pole())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaYAAACcCAYAAAA05w98AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAvxElEQVR4nO2de3xU1bn3v09IZnKDBIKAyEVIUqwKchW5WMRKtdpW0CPvebW01ls/HrFWOUGUawoISovW8qrV0+OLHi8Fq5a+4AVaQLFcFEELSA0JcpFySQDNbSaZzHr/mIuTyUwylz17dibr+/nMJ3vtvdfavz2znjx7PWvttUQphUaj0Wg0ViEt2QI0Go1GowlEOyaNRqPRWArtmDQajUZjKbRj0mg0Go2l0I5Jo9FoNJZCOyaNRqPRWIr0ZAswkmnTpqnCwsJky2iGw+EgMzMz2TKaoTVFRmlp6R+UUnckW4cZaNuJDK0pMuK1nZRyTH369GHBggXJltGM6upqOnfunGwZzdCaIqO0tPRosjWYhbadyNCaIiNe20maYxKRe4AGoDfwW6XU2YBjE4HR3uQ6pdSnIpIBzAfKgF5KqUeDy0xLs15kMj3der5fa2rfiMgFwDLg10qpzUHHtO0kEK3JHJJSG0VkODBaKfUc8Fc8RuM7lgk8BjwKPA6sEBEBpgMHlVIrgQwRuT643MbGRjPkR0VtbW2yJbRAa2rfKKX2A9WABO7XtpN4tCZzSNZj0mRgr3d7tzftYyxQqTw4gWygHzAF2BMmDwBffvllIrRqNFakIcS+mG0nVnbv3k1NTY1RxWk0QPJCeb2AD73bDjzhvMBjpwPSvuOB+4PzAPDOO+/geUD0UFJSQmlpKS6XC6fTCUBOTg5ut5v6+noAsrOzAairqwMgKyuLtLQ0/1OI3W4nPT3dn7bZbNhsNr8x2mw27HY7NTU1KKXIyMjAbrdTV1eH2+3G6XTicDj86bS0NLKzs3E4HLhcLn/a6XTS2NiIiJCbm4vT6aShwfO/Jzc3l4aGBn86JycnrntyuVzU19fHfE/p6elkZmYaek9KKaqrq5P2O4W6p3ZIzLbz3nvvxWQ7P//5z7n99tu54YYbtO1EcE/adiIjWY7pBJ6nOYDOwMkwxwKPt5YH8HzB4SalDe4czMnJaZb2/Xg+srKyWk0Hj4Kx2+0hj/tGzASf31Y6XHmBxHpPPk2x3lOs6dbK82lK1u8ULt3OiNl2rrrqKrZt2xay0NZ+k/Lyco4fP0737t39+7TtaNuJl2Q9Fr4ODPVuDwbWiEi+t5N2C3AegIjYgUalVLk3z7DAPMGFnjlzJsGyo8eKYY5U19TQ0MCJEyf47LPP+OKLLwwr16oYYTu+p+9oOHPmDFVVVRw4cCBW6a2S6vXUKKyoKV6S0mJSSu0SkX0icgeeGPgcYBGwWSm1WkQWiMiDQCfgTm+2p4BFIvIzPPH114LLdblcnD17lvz8fFPuQ5N4nE4n+/fvp6qqitOnT3P69GnOnDnj3w6VdjqddO3alW7dujFt2jRmz56d7NswFBHpD3wLGCci2zHAdmKhrKyM9PR0ysrKjChOo/GTtHGGSqllQbumBxxbB6wLOt8JlLRWZteuXSkrK2PUqFGG6YwXm82WbAktsLImh8PB9u3b2bx5M5s2bWLHjh3079+fHj160K1bN7p16+Z3OoWFhS32devWjdzc3Gb9JamGUuoQnoEOPuK2nU6dOkWto6ysjHHjxrFnz562T44BK9dTK2FFTfGSUgPgCwoKOHDggKUcU3Bc1gpYSZPD4WDbtm1s2LCBDz74gA8//JALL7yQCRMmMGPGDMaPH09eXl6yZaY8sTqm8ePHs3PnTs6cOUPXrl0N1WSleupDazKHlHJMeXl5lgsr1NTUWK7iJFNTfX09W7du9beIdu7cyUUXXcTo0aOZOXMm48aNo0uXLknR1pGJpY+prKyM733vexQXFyfkgVDbTmRYUVO8pJRjys/Pt5xjsuLS9WZoqq6u5siRI/5PRUUF77//Prt27WLw4MFcccUVPPTQQ4wbN47OnTtTWVnZbGSXxvqUlZVxzz33UFxcnJAQeke1nWixoqZ4SSnHVFBQYDnHlJGRkWwJLYhXk9Pp5OjRo80cz+HDh5ulnU4nffv29X/69+/P3LlzGTt2LLm5uYZr0sRHtO+eKKUoKyujuLiYoqKihNidFeuE1mQOKeWYunfvzttvv51sGc2wYhM7Fk1ffPEFv/jFL9ixYwenT5+md+/efqfTr18/Bg8ezLXXXuvfV1BQENUABCt+Tx2JaOdbq6qqQilF9+7dKS4u5q9//avhmqxYJ7Qmc0gpx2Sz2XA6nQnpiI2Vuro6y720GY0mpRTPPvsss2fPpqSkhGeeeYaePXvG1FlulCaN8UQ7V56vtSQiFBcX88wzzxiuyYp1Qmsyh5RyTIA/3n3ppZcmWwoAbrc72RJaEKmmI0eOcPvtt3PmzBnee+89LrzwwqRr0iSGaPspfI4JvrE5o7FindCazKFdTggWjrS0NP8IIatgxSnp29KklOL5559n+PDhTJgwga1btybUKUWiSZNYou1jCnRM55xzDo2NjZw+fbqNXNFhxTqhNZlDSt1Renp6wp7eYsWKTezWNB07doy77rqLo0ePsmHDBi655JKka9Iknmj/uZWVlfGDH/wAwB/OO3DggKGRCivWCa3JHFKqxdTY2JiwEUKx4pvR10qE0qSU4qWXXmLYsGGMGDGCHTt2mOaUwmnSmEesfUw+EvFAaMU6oTWZQ0q1mJRSFBcX8/TTTydbih8rxn+DNZ04cYK7776bzz//nHXr1jFixIika9KYSzR9TIFDxX0k4oHQinVCazKHlGoxBYYUrIIV1/UJ1LR69WouueQSBg0axM6dO5PilII1acwnmqH9J0+eJCMjg27duvn3JcLurFgntCZzSKkWU0ZGRrOO2EDDSRbBa6BYgezsbCorK5k+fTq7d+/mzTff5LLLLku6Jk3yiOYlzeDWEngc01NPPWWoJivWCa3JHFLK1bpcLn+rySr9TA6HI9kSWvCnP/2JIUOGcN5557Fr166kOyWw5vfUkXC5XBGfG84xGW1zVqwTWpM5pFSLyRdr9YUVRo8enWRF0Rl8tDQ0NPDVV19x9uzZZp/gfYHpqqoq6uvrWbVqFePHj0+YtmhJ5PekaZto+ilCOabu3bvT1NREVVUVBQUFhmiyYp3QmswhpRyTL05upRZTIuK/b7zxBtOmTcPpdJKfn+//5OXltUgPGjSoxb7u3bvTu3dvw3XFQyrGydsT0fQxlZWVccMNN7TI73sgNMoxWbFOaE3mkFKOyRcnLyoq4t13302yGg9Gx3+dTicPPPAAf/7zn7nyyitjWhDPik3/VIyTtyei6WM6cOBAixYTfPNAaFSkwop1Qmsyh5Rytb4mrZVaTE6n09DynnrqKS6++GK++93vxrxKq9GajMCKmjoSkYaDlFJtOiajsGKd0JrMwXTHJCITRWSW9zMkxPGbRWSriJSLyA0B+3uKSIWIHBWR34UqO7iPyQpE++Jia5w5c4YlS5awdOnSuMoxUpNRWFGTVRGRe0TkThGZLyL5AftFRPZ4beeAiHwtIoXeY9d5beeoiNwcXGakfUzHjx8nKysr5KrCRr/LZMU6oTWZg6mhPBHJBB4DLgVswHoRmaC8b/eJSAbgUkqNEZHBwHrgdW/2O4HBSqnatq6TiI7YWIm1VROKpUuXcv3113PRRRfFVY6RmozCipqsiIgMB0YrpX4iIuOB+cD93sMDgR8ppSq85z6vlCoXz5c7Fuir4lxVLtTABx/FxcWsWLEinuKbYcU6oTWZg9ktprFApfLgBLKBfr6DSqlGpdQq7/Y/gM8BRCQLuAIoF5GlIhLSodpsNrznWyacF2pRvFg4fPgw//Vf/0VpaWncZRmlyUisqMmiTAb2erd3e9MAKKXKA5zSMGCX99AlwDXAIRG5JVShPttpi7YcU1lZmWErqlqxTmhN5pCwFpOILACKgnbfDLwSkHYAvYFDIfJPBH4LoJSqB64SkTzgBWAm8Ehwng0bNrB48WJ/euHChbz22mu4XC5/HDYnJwe32019fT3wTcehb76prKws0tLSqK31NMzsdjvp6en+tM1mw2azUVNT40/b7XZqampQSpGRkYHdbqeurg63243T6aSgoMCfTktLIzs7G4fDgcvl8qedTieNjY2ICLm5uTidThoaGgBPxZs1axa33norNpuN+vr6uO7J5XKRl5cX8z2lp6eTmZkZ9z01NDT400opMjMzk/Y7hboni9IL+NC77bOfUEwFHgdQSu0GRohIMbBWRL5QSn0QeHKw7ZSUlFBaWtqinu3bt4/evXtTWVnZ4jfJzMxEKcXnn39OQUGBZWwnsJ7l5ORo22kPtqOUMu0DfBd4IyD9CVAY4rxuwLQwZfQC3gx1rKSkRPmYO3eumjdvnko2p06diruMXbt2qZ49e6qvvvrKAEXGaDIaK2oCFigT7SOSD7AQuM+73RU4EuIcAZ4Lk//fgV8G7w+0nda44YYb1B//+Mewx0eOHKm2bt0aUVltYcU6oTVFRry2Y/Zj4RbgPAARsQONyhMDz/B14opIN+B7SqkXvene3r++QGoPYG1bF7LaLOPx8OCDDzJ37ly6dOmSbCma5PM6MNS7PRhYIyL53v5ZHxOAzb6EiATaeR6wIdaLtxbKA2uNiNW0X0x1TMrTr7RARB4EZuAZ0AAwBVgkIucCG4H/FJGPROQgYBeRScBuEZmPZwDEc6HKD4yTW8VA4o3/vvvuuxw8eJC77rrLIEXWjElbUZMVUUrtAvaJyB3AVcAcYBEBfU3e7TcD0nNF5G8i8p/A50qpPcHlRtLH5Ha7KS8vp6goOEL/DUbanRXrhNZkDqa/YKuUWgesC9q3CljlTYZaBOhgmP3NaGpq8m8HdsQmc9RKQ0NDzAt5ud1uZs6cyZIlS6J6ATKRmhKFFTVZFaXUsqBd04OO/zIo3eaImUDbCcexY8fo0qULnTt3DntOUVER69atC3s8GqxYJ7Qmc7BsD28sBBpXQUEBIkJVVVUSFeHvoIyFl156iaysrBbTv8RLPJoShRU1dSQicUxthfHA2HcIrVgntCZzSKkpiQIREX8/U/fu3ZMtJ2ocDgdz5szhpZdeSsn3FDTtj0gdkxUiFZr2TUq1mILDXVboZ8rJyYkp35NPPsnw4cMTMgN4rJoSiRU1dSQiCRVH4pgKCgpIS0ujsrIybk1WrBNakzmkVIspeFoVK0xNFMuU9FVVVSxbtowtW7YkQJE1p8m3oqaORCRTEpWVlUW0dpcvUnHOOefEpcmKdUJrMoeUajEFx8nNGDJ+8803s2bNmrDHY5lgcfHixdx0000MGjQoHmlhseKkj1bU1JEwqo8JjHsgtGKd0JrMIaVaTMEkOpRXVVXFmjVreP/99/nkk0+YPXt23G89Hzx4kJUrV7Jv3z6DVGo08eN2u6moqGh1qLgPK4TQNe2blGoxhepjOnDggGFzdwWzY8cORo8ezY4dO3jrrbe46aabqK6ubnZOtPHf2bNnc99999GzZ08jpTbDijFpK2rqSLTVx3TkyBEKCgoiWvvHqEiFFeuE1mQOKeWYgh2QkR2xodi6dSuXXXYZ5557Lhs3bqSgoIAxY8Y0C2NEs2T1Rx99xKZNm3jggQcSIddPNJrMwoqaOhJtPbxFGsYD40J5VqwTWpM5pJRjCtUJmMh+pm3btvk7g+12O7///e+55557GDdunH8FXd9Eim2hPPOVsWDBgoS/yR2pJjOxoqaORFsd6GVlZRGF8cC4WcatWCe0JnNIKccUikTFu91utz+U50NEuPvuu1m9ejW33norv/71ryM2zrfeeovjx49z2223Ga5Vo4mXaFpM3bp1Iz09nVOnTiVYlSZVSSnHFCpOnqgh4/v376egoIAePXq0OPad73yH7du388orrzB9+nT/dPPhaGpqYubMmTz66KOkpyd+PEok/QRmY0VNHYm2+piicUxgTKTCinVCazKHlHJMoUhUKC8wjBeKvn37smXLFkSE8ePHc+hQiyWn/KxcuZJu3brxwx/+0HCdGo0RROuYrPAOoab9klKOqbGxscW+RIXy2nJM4FmQ63e/+x233HILl112GZs3b25xTl1dHfPmzWPZsmWmTeHSVgsuGVhRU0cilO34cLlcfPHFFxQWFkZcnhF2Z8U6oTWZQ0o5plAkasj4tm3bGDNmTJvniQgzZszghRdeYOrUqaxYsaKZlieeeIKxY8c266vSaKzE4cOH6dGjB1lZWRHn0e8yaeIhpV6wDdU/E9gRG6o/KBaqq6upqKhgyJAhbZ7rM+ZJkybx97//ncmTJ7Nr1y6eeuopvv76a5YvX8727dsN0RUp0fyDMQsraupItNa3GW0YD4wJoVuxTmhN5pBSLaZwoTCj+5k+/PBDhg4dGtHiaoEzQRQWFrJ161bOnj3LhAkTuP/++7nllluiCpEYQbyzUyQCK2rqSLQWRo7FMRkRqbBindCazCGl7ihcnNzosEIk/Us+amtrm6Vzc3N57bXX+OEPf8jbb7/NnDlzDNMVKcGarIAVNXUkWutjisUxde3aFZvNxsmTJ2PWZMU6oTWZQ0yOSUTa1QJHRo8QisYxhUJEmD17NsePH497BmZN+0NEckVkqoiM9aYHiMgVyVUVnlgcE5gzibImNWnTMYnINSJybeAHmBfrBUVkoojM8n5CdtKIyCIROSoih0Wk0LvvHhG5U0Tmi0h+qHydOnUKeU0jW0xKqagck91uD3vMjHeWQtGapmRhRU0JZBZQD1wkIg8Ah4CHI83cli2IyAav/WwTkU4ikuG1qZ+KyIOhygxnOxC7Y4r3gdCKdUJrModI/jMOA44F7cuL5WIikgk8BlwK2ID1IjJBBQSiReRc4KRSqk/AvuHAaKXUT0RkPDAfuD+4/HCxViOf3A4ePIjNZqNPnz5tn0zynE9raE1JZyewRyn1FxHJAn4ERNSD3ZYtiMhkYLZSanvAvl8AB5VSK0Vkjohcr5T6c2C54WynsbGRw4cPM3DgwOjukPgfCK1YJ7Qmc4jkjv6vUupfgTtEZH2M1xsLVHodkVNEsoF+eJ4YfVwN/IeI3ArcpZT6CJgM7PUe3w28SAjHtGHDBh555BF/uqSkhNLSUnr16kVZWRmnTp0iNzcXt9vtn1/K99a0712ArKws0tLS/HFbu91Oenq6P71x40ZGjRrlnxjWZrNht9upqalBKUVGRgZ2u526ujrcbjc1NTX06tXLn05LSyM7OxuHw4HL5fKnnU4njY2NiAi5ubk4nU4aGhoAT79UQ0ODP52Tk4PL5fKvw5KTkxPVPdXV1XHOOef40zabDZvNRk1NTUT3lJ6eTmZmpqH35HQ66dKlS8z3FPw7GXFPCWQ/8L+BR5RS9cCb4aIAIZhM67YwEbhRRHYAdyulTgBTgJKAPDcCzRxTONvZv38/vXr1orq6GrfbHdVv0qdPH9auXdvMVqL5TbTtdGDbUUq1+sHTwhnr3R4PZLSVx3vuAuB/gj5u4KWAc7YAY8LkvxY4AnQBngXu9O5PB5yh8pSUlKhwFBQUqOPHj4c9Hin33nuvWrZsWcTnnzp1Ku5rGo3WFBnAAhVBXY/2AywLsKnLI7Up7/lt2gLQCZgLvOtNfw4Ue7evAt4JzhPOdtauXasmTZoU0/e3Y8cONXTo0JjyKmXNOqE1RUa8thOJa9uH5wkPYDvwvyJ0eAuUUj8O/ACTgMCJnToDIYftKKXWAWuBgcCJgHxh87QWJzcqnBftwIdIhpSbjdaUdPbyjU1tI0Kb8tKmLSilmpRSC/E81EWUJ5ztxNq/BPEPGbdindCazCESx9QfWCkiv8TT35Qfx/W2AOcBiIgdaFRKlXs7Z/O9+wM11eIx4teBod59g4GQa5m35piMGABRX1/P3r17GT58eMR5rFhptKakE49NtbAFEckXkQz4xn68/bnvBeQZFpgnuNBEOKb8/HwyMzM5ceJETPmtWCe0JnOIxDEdAG4HjgK3ATEPAVFKOYEF3pFBM4A7vYemAIu822tF5GURuQ94RinVqJTaBewTkTvwhCJCvvzji7mGwogh4x9//DHf/va3o5rN1xentRJaU9KJ2abC2MIiYLKIFAB7ROS33nIXerM9hWcE4M+ABuC14HLD2U48jgnieyC0Yp3QmswhksEP9UCTUuo1EfkIT99PzHhDdOuC9q0CVnm3vx8m37J4rltcXMwbb7wRTxFxv7+k0XiJy6ZC2ML0gO0LQ5zv5JvBD1ERr2PyhdAvv/zymMvQdDzabDEppV4HzvUm8wDLLjCf6D6mWByTFZvZWlNysaJNhbKdhoYGvvzySwYMGBBzufFEKqxYJ7Qmc4hoXJ9Sao/37yfxtlwSSVt9TPHO3RWLY7Liy29aU/Kxmk2Fsp2Kigr69u3b5iKCrRFPKM+KdUJrMoeUmiuvtT6meDtijx49isPhiHrCVSvGf7UmTTChbCfeMB7EF6mwYp3QmswhpRxTW8RjJNu3b+eyyy6LejG/eFpoiUJr0kTCgQMH4nZM8UQqrFgntCZzSCnH1NYbx/GEFWId+BBPGCRRaE2aYELZjhEtpry8PLKzszl+/HjUea1YJ7Qmc0gpx9TWnFHxdMTG6pisGP/VmjTBhLIdIxwTxP5AaMU6oTWZQ0o5ptbWlIHYDaSxsZFdu3YxatSoqPP65qeyElqTJphQtmOUY4o1hG7FOqE1mUNKOaa2Yq2xGsinn37KgAED6NKlS9snB+F2u6POk2i0Jk0wwbbjcDg4fvw4/fv3j7vsWCMVVqwTWpM5pJRjiqSPKZaO2HherLXilPRakyaYYNupqKigf//+hvwusUYqrFgntCZzSCnH1NYPFGtHbDyOKTMzM6Z8iURr0gQTbDtGhfEgdsdkxTqhNZlDSjmmtvqYIDYjiccxWTH+qzVpggm2HSMdU1FRUUyRCivWCa3JHFLKMUVS8aPtZ6qsrOTkyZN8+9vfjkmTFeO/WpMmmGDbMdIxdenShdzcXP71r3+1fXIAVqwTWpM5pJRjiuTl12g7Yrdv386ll14a86qMCV4JNSa0Jk0wwbZjpGOC2CIVVqwTWpM5pNQdRfKiWbQGEu+M4tEskWEWWpMmmGDbMdoxxTIi1op1Qmsyh5RyTC6Xq81zojWQrVu3xuWYHA5HzHkThdakCSbQdurq6qisrKRv376GlR/LkPFwdaK2tpbVq1cbIStqrFhPragpXlLKMUUSa41myHhTUxM7duxg9OjRMWuKxFmajdakCSbQdsrLyxkwYECrs/VHSyyhvHB14umnn2bq1Kl88sknRkiLCivWUytqipeUckyR9DFF0xH72Wef0bNnT7p37x6zJivGf7UmTTCBtmN0GA+M62NyOBwsX76c22+/nXnz5hklLy5NycaKmuLF9DsSkYkiMsv7GRLi+B9E5AsROSAiJ0TkNu/+niJSISJHReR3ocqOdDLDSI3EiBVrrRj/1ZraNyJyj4jcKSLzRSQ/6NhjIrJLRD4WkW8H7L/OaztHReTm4DIDbaesrIyioiJDNRcVFVFeXh7VCLJQdWLlypUMGzaMFStW8PHHH7Njxw4jZcakKdlYUVO8mOqYRCQTeAx4FHgcWCEBj2oikg78j1LqfKVUEfAy8Lr38J3AYKVUH6XUvaHKj7RJG2k/07Zt2xgzZkxEZYbD6XTGlT8RaE3tFxEZDoxWSj0H/BWYH3DsEuAFpdQw4DlgoXe/AGOBvl77eTm43EDbSUSLqXPnznTu3DmqIePBdcLlcvHoo4/y8MMPk5mZyZw5c5g7d66hOqPVZAWsqClezG4xjQUqlQcnkA308x1USrmUUhvB76SylVJnRSQLuAIoF5Gl3mMtiPRpLNKOWCNaTJG89Gs2WlO7ZjKw17u925sG/Kvh7vEm/w4c825fAlwDHBKRW0IVGmg7iXBMEH04L7hOrFq1ij59+jBu3DgAfvazn1FWVsb7779vqM5oNFkBK2qKl4RNsiQiC4DgeMDNwCsBaQfQGzgUooirgXcAlFL1wFUikge8AMwEHgnOsHnz5max8pKSEkpLS3G5XP6nipycHPr27cuqVauorKz0N4N9b09nZWWRlpbGv/71Lw4ePMj5559PfX09tbW1ANhsNmw2m3/VSJvNht1up6amBqUUGRkZ2O126urqcLvd1NbWkpub60+npaWRnZ2Nw+HA5XL5006nk8bGRkSE3NxcnE6nf1XR3NxcGhoa/OmcnJwW9+R2u6mvrwcIe0++e6ivr4/rntLT08nMzDT0nhoaGqiuro75nux2O+np6Ybek0XpBXzo3fbZTygmAcsBlFK7gREiUgysFZEvlFIfBJ4cbDurV69m7NixcdWz4N/k/PPPZ9++fVx88cVAdLbjcrlYtGgRS5cu5ezZs/569tBDDzFr1izefPNN0tLStO0Y8DtZwnaUUqZ9gO8CbwSkPwEKw5z7FGAPsb8X8GaoPLNnz1aRsHPnTjV48OBWz1m/fr26/PLLIyqvNRwOR9xlGI3WFBnAAmWifUTywROeu8+73RU4EuKcS4ArwuT/d+CXwft9tlNdXa2ysrJUU1OTkV+lUkqpxYsXq5KSkojPD6wTa9asUcOGDVNut7vZOY2Njepb3/qWWr9+vWE6I9VkFayoKV7bMfuxcAtwHoCI2IFGpVS5iGQEduKKSI73mDNgn+9xrgewNlThTU1NEYkoLi5usyPWiDAeWDP+qzW1a14Hhnq3BwNrRCRfRDIARORiIE8ptUlEuohIrogE2nkesCG4UJ/tHDhwgIEDByakxRjtu0y+OqGU4pFHHmHWrFktRt6mp6dTWlrKnDlzTFli3Ir11Iqa4sVUx+R1NAtE5EFgBp4BDQBTgEUBp04G3vQlRGQSsFtE5uMZAPFcqPIjdUyRdMQa5Zh8zW0roTW1X5RSu4B9InIHcBUwB4/tTBaR7wD/D1guIh/hGRxRC8wVkb+JyH8Cn6tv+qH8+GwnUf1LEH0fk69ObN68maqqKm688caQ502dOpXa2lrWrg35vGooVqynVtQUL6Yv5KGUWgesC9q3ClgVkH4p6Ph6POEJw/AZyXnnnRdKI9u2bePZZ5818pIajSEopZYF7ZoesH1+iCylkZadSMdUWFjoj1RE0yJbsmQJDz74YNgXftPS0li4cCFz587l2muvtXL/oCZCUuoXtNlsEZ/b2pDx8vJysrOz6d07XL9y5OTm5sZdhtFoTZpgfLaTSMfUuXNn8vLyOHbsWNsn46kTO3fuZN++fUybNq3Vc6+//no6derEG2+8YYTUVjVZDStqipeUckyRhvKg9Xi3UWE8sGYzW2vSBGNGKA+iC+c1NDSwZMkSZsyY0eZDp4iwcOFC5s2bF9X/gWixYj21oqZ46dCOKZyBaMdkPlbU1JGwomPau3cv7733HnfeeWfbJwPXXHMN+fn5vPrqq/FIbBUr1lMraoqXlHJM0dBaKM9Ix6TRtBe+/vprampqDAlhhyOa2f2ffPJJ7r33XnJyciI6X0RYtGgRCxYsSMmJTTsSKeWYIp0rD8LP3VVXV8dnn33GsGHDDNEUqVGZidakCSYjI8M/R14kkyHHSqRDxg8fPsy7777L9OnT2zw3kIkTJ9K3b19eeOGFWCW2ihXrqRU1xUtKOaZoJogM1xH78ccfc9FFF5GVlWWIJis+uWlNmmDcbnfCw3gQeSjvN7/5DdOmTaNr165RX2PhwoX86le/Ssj7PVasp1bUFC8p5Zii7fQMZSTxLgwYjBVfftOaNME0NTWZ4pgimWX81KlTvPjii9xxxx0xXWPcuHFceOGF/OEPf4hVZlisWE+tqCleUsoxRUuoeLfuX9J0VMxwTDk5OXTt2pWjR4+GPee3v/0tU6dOpVevXjFfZ+HChSxevNg/X5ymfZFSjimaPiZoGe9WShneYrJi/Fdr0gTj62NKtGOC1vuZvv76a5555hlmzpwZV50YMWIEo0eP5umnn465jFBYsZ5aUVO8pJRjinaurOBQ3tGjR2lqamLAgAGGaYqm38sstCZNMEopUx1TuH6mp59+mquvvpqBAwfGXSdKS0t57LHH/LNkx0tjY6Mll5hIRdtJKccUbSdgsIH4wnhGjkqyYihBa9IEU1NTg9PppGfPngm/Vrgh4/X19TzxxBPMmjXLn46HwYMHM3HiRH73u5ALXkeMy+Xi2WefpW/fvowaNYq33347rvKMJhVtJ6UcU7QEzt0Fun9J03E5c+YMxcXFCR0q7iNcKO/5559n1KhRDB482LBrLViwgOXLl3P27NmY8r/zzjsMHTqUl19+mbVr1zJv3jzuvfderrvuOv75z38aplPTnJRyTNH2MeXm5tK1a1e+/PJLIDGOybdIl5XQmjTBfP3116aE8SB0KM/lcrFs2TIeeugh/z4j6sSgQYP4wQ9+wOOPPx5Vvr179/L973+fe++9l8WLF7Nx40ZGjBjBlClT2Lt3L1deeSXjxo3jgQceiNnpGUUq2k5KOaZY8BlJQ0MDu3fvZuTIkcmWpNGYzunTp01zTIWFhVRUVDTrG3n11Vc5//zzGTNmjOHXmzdvHitWrKCqqqrNc0+ePMndd9/NxIkTueaaa9izZw/XX399s5akzWZjxowZ7Nu3j5qaGi644AKeeeaZlHyfKFmklGOKpWPSF+/+5JNPKCoqonPnzoZq8i1nbCW0Jk0wVVVVpjmmnJwcunXrxpEjRwBP5/2SJUuatZbAuDoxYMAAbrrpJpYtC14t5BscDgdLly7lwgsvJDMzk/3793Pfffe1mDw2UFOPHj149tlnefvtt3n11VcZPnw4f/vb3wzRHA2paDsp5ZhiwRfv1v1Lmo6Mr4/JLAL7mf7yl7+QmZnJpEmTEna9OXPm8Nxzz3H8+PFm+5VSvPrqq1xwwQXs2LGDrVu38vjjj9OtW7eIyx46dCgbN25k/vz53HHHHUyZMoXy8nKjb6FDYfpCgYkkPT362ykuLmbbtm0cO3aMq666ynBNRk1tZCRakyaYs2fPmu6YysrKuPLKK3nkkUd4+OGHWwy8MLJO9OnTh2nTprF06VKeeOIJwDPLywMPPEBDQwMrV65kwoQJbZYTTpOIcOONN3Ldddfx+OOPM3r0aG6//XZmz55Nly5dItaplOL06dNUVFRQUVFBeXk5J06coKioiCFDhjB48OAWTjMVbSelHFMsI4p8BuJwOJgzZ47hmqy4mqbWpAlGKUX37t1Nu57P7jZu3MhXX33FlClTWpxjdJ2YNWsWF110Ef/2b//GihUr+OCDD1i8eDE//vGPI75WW+dlZmby0EMPceutt/Lwww9zwQUXsHDhQm699Vb/CryNjY0cOnTI73wCnVBFRQXg6YcrLCxk4MCB9OvXj7179/LKK6+wZ88eunTpwuDBgxkyZAhDhgyhuLiYoUOHRrVQqtUx3TGJSA9gJlCllFoS4ngf4D+AL4CvlFJ/9O6/B2gAegO/VUqdDc4bSx9TYWEhZWVlZGdnM2jQoKjzt0Vtba3lnmi0pvZNa7YgIhOB0d7kOqXUpyKSAcwHyoBeSqlHg8vMz883Zai4j6KiIrZs2cKSJUuYNWtWyH/4RteJXr16ceedd3L11Vcza9Ys/vu//zvqEW2Rajr33HN5/vnn+eijj/jlL3/Jk08+Sffu3amoqODYsWP07t2bgQMH+p3PyJEjGThwIAMHDqRr165hfwu3282hQ4f4xz/+waeffsqaNWvYtWsXR44cadaq8v3t06cPSikaGhpoaGjA6XTidDqbbQenfduDBg0ybJWFaDHdMSmlTorI53iMKhT/B3hAKVUuIn8WkfeAc4HRSqmfiMh4PEZ2vxF6cnJyOOecc7j44ov1U7vG8ojIcMLYgohkAo8BlwI2YL2ITACmAweVUitFZI6IXK+U+nNgubHM4h0PxcXFbNq0ifz8fG6++WbTrvurX/2KmTNnRtWHFA8jR47k/fffZ/369YiIvwUU7astPtLS0hgwYAADBgzgRz/6EQCVlZXk5OTw2Wef+R3WE088waeffkplZSVNTU3YbDbsdrv/b2vbvvSUKVM6jmPyEnLJRRGxAxOVUr6ewzLge0AhsNe7bzfwIiEck6+pHC3FxcUJG/hgt9sTUm48aE3tmsmEt4WxQKXyzM3lFJFsoB8wBSgJyHMj0MwxVVdXN3tKLykpobS0FJfL5Z+9OicnB7fb7Z9pwNfa8I0Ky8rKIi0tjdraWsDzm6anp/vTNpsNm81GTU0NeXl5VFdXs2DBAqqrq1FKkZGRgd1up66uDrfbjdPpxOFw+NNpaWlkZ2fjcDhwuVz+tNPppLGxEREhNzfX/9QPnncVfa0F3z1kZGRQWVkZ0z25XC7q6+tD3pMvbbfbqampaXZPI0eOxO12k56eTlNTE9XV1Ybdk09Xv3796NevH1OnTvXfU0NDA3l5eYhITL9TZWVl2Hvy/S7p6elkZmY2+53iJWGOSUQWAEXB+5VSP24lWwFQHZB24GlZ9QI+DNrXgk2bNsVkXLfddhuFhYVUVlZG9aP50q39aE1NTWRkZBhuXPH8wwDiuqdQFTHee+rUqRPV1dUJ+ScY6z1ZlNZsoRdwOiAdaD+nw+QB4Kc//SmLFy8OecHgVyiCJw0NDocFh7qC05mZmQAsX76cn//85y3y+47X19eTmZnpTwcfD5cOfsgJPg6x31N9fT1ZWVlh7ylSDUbek09Ton6nSDSESsdDwhyTUmpBDNkqgcC77wwcBLK9H9++k6EyX3755WzdujVkwa39aD/5yU9anG/Uj1ZZWWk54/I5YCsZV2VlJV27dk0p40oQJwhvC4HHAo+3lgeIrX82Xu6/v/VovBX7HbUmc7DEY6GIZIhIvlKqAfhARM73HvoW8A7wOjDUu28wsCZUOR988EGClUbPY489lmwJLdCaIuaKZAsIQQtbEJF87wCHLcB54A+LN3rD4q8DwwLzBBeqbScytKaIuSKezBLtUhHxIiL5wGI84YX/UEqdEJGpwHeUUtNFZABwN1AOVCulXvbmKwHO4ImZP66UOhOibGX2/bSFiES9HEei0Zoiw6vJvKFqERJsC8BCYLNSarWIXIvH+XQC3lJK7fI6qUXAPqAn8GiwoWjbiQytKTLitR3THVMi0cYVGVpTZFjVMSUCbTuRoTVFRry2k1Iv2AK7RCRkmC+JXCEim5ItIgitKTJ+lGwBJqJtJzK0psiIy3ZSqsWk0Wg0mvaPJQY/aDQajUbjQzsmjUaj0VgK7Zg0Go1GYyna1eCHRExemWBNjwGTAAXcopT6zLv/OuD33tNm+obEm6HJe3wDcAFwFBiH5wElKd+TeKbq+AeQhed76gEM886VmOjv6QJgGfBrpdTmoGNJqU+JQttO/Jq8x7XtYILtKKXaxQcYDrzg3R6P510m37FMPNO0CJ6ZI97zbt8P3O49Zw5wvYmaLgEu9m7fDbzm3RY873GJ2d+Td99kPJOABu5L5vdUCAwMSD9vxvcUcL2XgSuC9iWlPiXwHrXtxKnJu0/bTnNtCbOd9hTKm0zzySsnBxzzT16plHLimX7FN3nlnjB5EqpJKfWJUsp37b8Dx7zblwDXAIdE5BaD9bSqyctE4E8i8rqI9PTuS+b3VK6UqgAQkWHALu+hRH9PPkJNKJys+pQoJqNtJy5NXrTtNCdhttOeHFNrE1HGPHllAjUFMglYDqCU2q2UGgF8F5gvIuPM1KSUug/oj6cSvxhJnkRrCmAq8KpXZ6K/p9ZIVn1KFNp2DNCkbSciDKlP7ckxJWTyygRqAkBELgE+Ukp9EbhfKVUGzANGma1JKdWklFoIdIk0T6I1eePl3ZVSzY4l8HtqjWTVp0ShbccgTdp22sSQ+tSeHFNCJq9MoCZE5GIgTym1SUS6iEiuiAR+53nABpM1pXn/ZuKJ//ryJO178jIB8HeimvA9tcAC9SlRaNsxRpO2nTAYXZ/a1cwPiZi8MlGa8DwhvIBnKQ/wjJq5FM8TzARgHbBTKbXRSD1taPob8D6wHvgn8KJSqjqZ35NSarX3+BPAHKVUjTc9n8R/T/2BV4C1eMJFy0hyfUoU2nbi1qRtp7mmhNpOu3JMGo1Go0l92lMoT6PRaDQdAO2YNBqNRmMptGPSaDQajaXQjkmj0Wg0lkI7Jo1Go9FYCu2YNBqNRmMptGPSaDQajaXQjqmDISKZIvJ6snVoNO0NbTvmoR1Tx2McsD/ZIjSadoi2HZPQjqkDISJXAbO822OSLEejaTdo2zEXPSVRB0NE/gr8SClVm2wtGk17QtuOeegWUwfCOytyJ21YGk10aNsxF+2YOhYjgQ+9U9SbuUaLRtPe0bZjItoxdSxOADbg+0qpD5MtRqNpR2jbMRHdx6TRaDQaS6FbTBqNRqOxFNoxaTQajcZSaMek0Wg0GkuhHZNGo9FoLIV2TBqNRqOxFNoxaTQajcZSaMek0Wg0GkuhHZNGo9FoLMX/BxEIYcNONNHjAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x165.6 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1,2, figsize=(6, 2.3))\n",
    "\n",
    "Td = np.arange(0, 1, ts)\n",
    "y, t = initial(G, Td, [0,0,1])\n",
    "# ax[0].step(t, y[:,0], where='post', ls='-', lw=1, color='k')\n",
    "ax[0].plot(t, y[:,0], ls='-', lw=1, color='k')\n",
    "\n",
    "# ax[1].step(t, y[:,1], where='post', ls='-', lw=1, color='k')\n",
    "ax[1].plot(t, y[:,1], ls='-', lw=1, color='k')\n",
    "\n",
    "plot_set(ax[0], '$t$', '$x_1$')\n",
    "ax[0].set_xlim([0, 1])\n",
    "plot_set(ax[1], '$t$', '$x_2$')\n",
    "ax[1].set_xlim([0, 1])\n",
    "\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}