{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 141,
   "id": "36f50fc8",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sp\n",
    "import numpy as np\n",
    "from scipy.integrate import odeint\n",
    "from matplotlib.pyplot import *\n",
    "import math\n",
    "import control as ctr"
   ]
  },
  {
   "attachments": {
    "2024-01-23_13-32-20.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAAg4AAAHuCAAAAADSfNv6AAAAIGNIUk0AAHomAACAhAAA+gAAAIDoAAB1MAAA6mAAADqYAAAXcJy6UTwAAAACYktHRAD/h4/MvwAAAAd0SU1FB+gBFwwuFamtYXoAAC8CSURBVHja7Z13mFXF+ce/M3Pu3UZnpbOAVClSpFiootKb9KKJsT4aY35qjDEmmsSoaNTExMQYY+x0kK4gSBMVRIr0LohIX8rusveeM+/vj92Fu7u3lz1zz53Pk+eJt+xhZs7nTnlnzgwjaDQlcLsToFEJrYPGB62Dxgetg8YHrYPGB62Dxgetg8YHrYPGB62Dxgetg8YHrYPGB62Dxgetg8YHrYPGB62Dxgetg8YHrYPGB62Dxgetg8YHrYPGB62DxgfD7gTEAQKY/xeaCGEOW1hP2oVYcEDt4ClkLKv4v4kVemWGy+4kJS1J33eQONy2UYOvIAGAWH73Oi33FL3QRE7S68Bl01anc98vMkBi+dd5TVvrFiNakl4HEO4EFl4wijpBU4CfwNI6REny68DRrwEOfAYJSHF0CWoOdkCm7CL5S45ZVUcAUwGAsOgk+te1kj9TduGIkhsPLDkmJBg+BCbanZpkxgE6CFzXFScXgyTf/jla9nFCnuzCCUVnYQIwHUxiViFGpJu6Ixk1TtCBYXgVrNzLhZwJMdYRWbILJ5Qdtxr1Q/5csNVb0L2DlJa0LIfF3isKJ+gAYBIwkzAVGA8YggvBdGAyGhwxhUWsoN0+trFdk0M1NzXIX7zsTM2bhjHpENMrlLJTWJSM08PMyhg1mZbmH8JNDbaM2YWsvNd6Tq+ttg9qBtLLFBmxpKwuGMYaWD4b+Ckm7Hr4aO6+UavutTtNIZOsIqXvvuTfnOmr9q/KP8R6rcrGyaZ7p40buBCAt82ejR3UzQixw+tGqlg/lC4xwsGdSMbqQWISTp3CrTiAYfDC4+qOfQpnhHBmg5LJK60DgzvN7hRFBcOQGsTYeDx87na4YGA/atidpqCIDLtT4Jey9Smp6GwY2bDqjIDs3pHcldMBL1+2suE1ag+i1RwIK11kEVEHGA8LBHhd39+F31XRqx4ixyE6kPDMQvZwcDB4XQf7HPzl3VLYnagkxBE6kOnB0p0YUNfigOXacMPeB15RdCSnOE7QgZiRdvJJ4G4AkGJWtx/+/A9IrUMUJP/CemLnnnB75n2PQT2kgOSvPoT3JklHeF7xJL8ODFXnHwZQ9a8ALPHWQ5UX9lA3AKU4Dig3ibYArlzczOJSbLvbPb+HB0QyOUfMNpP8tQMYXlt+ttnNGZYA8HuZ8Vi+APj5l4dYao0tVAxKl8UJOlCTOwFIAeJ5u9PlBg6AF+bana5y6bSY8nWxA3QAIwuMc4Ah4wtZXOaUDrUqh4uFVQGTq22EE3QAu5QLXsnutATAEnMmj7qlkwGYQuFGwxE6+HC5A6laobu37Xq90039WxhQuJdbTgdK1lms8hlRC8YvHp2/7K/X9+9bF4Bld3L8Y5R/w2kVhhIICICQ/92BeU16DbiuWmXVKq8iSt97Ap3xnFarC2Yv5JJxmRm1KucCAJPs/NYdU64aUFvN6qFsVSAXnytQU1xbIPfm6jmeOBQIuTcQCAQwKU9/efTMSLuz5pfSOjCIsXcrFr2xFclfbh+XxaOWmLkcDCBiqFz/2pGVVtudNb/ojkIIRHwCBayo7wC4s7vecktT7FZzNZS/kYXdaVIIAsWnQAhgELxq274DO3BI7rE7Z/4ppwNTbsBuM/EqDo+scWXPQd2yAJOrFxUpRjcWFQNDpQF3XF8LsBg3FF03C61DRcExcBAgiavdT9c6VBSMrDh1SxOI1qHCYElQ1sr7qqlItA4aH7QOGh+0DhoftA4aH7QOGh+0DhoftA4aH7QONqFmwauZKudj5tudAr9oHeyAod5gJee4kyCO7kAYXXGFko9s6trBFhipuR2Jrh3sgSlpg64dNL5oHUKg5q84UWgdQuBR83GpBKF1CApDs+xUqiB0VzIoHCNS6ieTQlnVhEbXDjZAxacKqdcKaR0qGJLERLEJZIGrFYDQOlQkZAkmAHjO5TNeNUMYACyVtjTTOlQckoSBsxs3b95x6EcJVGpw5dVXX9NcgCRXpYrQOlQUZBnInT1vZS6AKq2rwPph785FQLubx3YVUGVPDa1DBWEaxqHX/ncC2aO6tW1VK00AnrwD2zav2PDty71/PlIQKTHGK3tS3vwf9e4vCUAydvqlv+anTZzUpczOlwcXvbkRXX49EqYKv0wlnHQ8FmfvX/0sPbH3v30qSdOSVIQ0TTS+/5uP+64fNWK/ocK54VqHCsAUF0bdduTeXX9uYFnEDcFZEdwwSJro9+nH7T9qN00osO+O1iHxmMbBHrPar3q9oSlF2R2IGTdgUb/1T+WPe57Zf4681iHhmMbW7ptGrO1hSsN/aQtmup7+KOM3D4Pb7YPWIdGYxoZeRx6cnWkagcvaIHPYmjqv3Atuc3uhdUgwprGr3+nfvAoZdODADG+nlY3f+KXdW2lrHRKLNI4NPPXrZynkVqUub4tl9f72rM3dB61DQiGOn+z/2fPhlLPLe+X8rN8uFqad6dU6JBQLkz+59s3wDvV0eTu9jbuPG3bWD1qHRGIZO55wv8ms8ErZZY164MhvbO09aB1CENPNYXhKPtPGDDfqz/Bsg7fWChsX62odQhDL1LPFV85o+fPw1zNws8of8Qc774nWIQQXLkb/twx/w5MZZpFRREQAyDJN89L0hDRN07x8NpbATzsuWcPsqx60DkGRePerqHeQlnzjnOajSyoHxhgDWUwYhiGYJAAWuGEYBr90/5nJfoG/27iGUoVZVYUhFFpRdx8kfx/3p5VMXJ/1ID2Li5NrN11Ia9O9ASS3BDZ+fpjX7txdlCwqEBj91MwDTew7Q5x8sWjeG2SSpgST/ro82gKRdCHHOEgWEZGU1I/zn1De4zUAABn355OXvuhddA86rSz6FhF56Vd41b5boBuLEEQ/spBYe2hQo0s/9Dwpkd/z+dMAhFHwz2HSePO6FQDg5t/0mseL2wuO0ZhmXxOuG4uEQViNQfDpJuLUzRuAxo3ObzaFWPpC07uB7GZi+xm4Pfd1r1G0/QdDuys3Hqtt114gunZIGBxfosPllwSxcC16r9i+4utve1kWnrkLtd/cs3bNvt/D4z46pbjDyqz0rvlbbDv/ROuQKIjnbruiZakC5rj1s14ZwFUz61g871yDz++sxlD9D48C+LRkOEHoiq22RSa1DomCcPqHxlV8a31m1XsLXiJ4svuBC/ylaSERTNzHvdhdcicYWmGnbbdF65AoCIfQutTPXGBiVdPFGBPoBliNh8HNGASq1SOcLyj+KkMz7AwQeQgyJEDQAUPYidZdycRxAnXL1PrdSv6jFhg6pBePOjKqfQ+6HIms4joWQIdg/cv49D21DonjIjJLvTZRp+SuMQCNSt7nvlU0g8go8H89GezwTUFBup/cHWaStQ6J4yLSSr22Kmf76pBxqeoo/cs23F5LlB9pEvvh1YBVgMxYW6tpYYCPmbf+/4VZeWgdEkcaCku/IVzh/JnlyRR+6n6Gur8L9CckvJ81+FOuCPhxuE2J1iEEMbTJGSiI/I+ImQV1/X8kKgf+s1Ud8s5Viz23WocQmNHPNmfjx2hsOu+9Av6jkgFHCJaxdv35DX0Cr7PRtUNc4BhaNcrBOEMOtkeuA2EfWgXQIdDFyDj19baMz/rEvk2EjjsEhVHzWlHOHzDUqHPwAos0vkjYiVYRBqkl1u24kL/ux9gf2tE6BIdRtLNJTFZvfWxXxLMPDOvRNsJKhWHZaeR/93nsUx1ahxBEv5WXxLXYHOkfkShcl351ZLeF+KFNuaCTK2O/m1qHhMHQA4si7TwQtu3pUDeyGklizYGLDHkb98b8DJfWIWFw3FB/wfcRPpQtMRNjIqz0BVacBXjBjytingnVOiQMZlYeWTg3sltLRsE0DImsSpHY+u15EOHU52E/0REIrUPi4JiEf3qN4l+sgDAu3WYGw6foBUqWW1uYvX9Es8hWzkqsPOphIORt3RxrZ1LrkDi47DJo+2wUx7FyYZ28FNIqRAHyLn3xNHCKAIAMvIoHI6vyybi4JhcMTHhOLI+1tdBhqARC+OXCP9/qIgYGTOrozqyG4gWRLe/MON+7pE0Qdx3kNd1ggGV8sK5PH4qoyie2Yd95SJCFk+vPVo1tlaXeSDCRSD587t9+4S03c1V0z8rdOcnz2+3/rHdkN8A0/viPPAkvM8hV+f0bY7t7WodEYonNHbI2NbMEAEm4/LQmSYBdaqhLPjKNX/1l4vsRPnND7MfjGWb1l6+65TwzcyrHVjvoxiKRCKv900/f/ZmQvEwvjZX6yRV/5HUt/Ev25Ej/DUZ16gBo2Ko+AMS4Il93JRMKx1O9VjyIsGYuTNf22/BGfTPSW8JISi88F2FKGevzGVqHhMIsvFv/H39gYaxeNY0fBp/57Qgr8vqacc7BOTiPeUpT65BYhJmzsPLTz/CQ4WOvcajvgbufsfmGaB0SjGG2X1Tpd4+DB90CjLyuHb12TnojvF2kEofWIdEYZvfPak6+zTTMwBWEyVxLrz9473uw70n+IrQOCccwO69p+X6PDQY3/fcgLGlYk2/J/d3rttugdagADLPV2v5fdn70uMHMcqcUkEmCf9L5cfHWH0PvRZpwbE+A6kTwRFtADKvG4n/VeKnVcycMwaRpFT0mV3SeBTPYqoH9Nw3ZeodU4GbYnwK1IRbxakc/CIn7tv7i7BMt7l3t5YZgjBGBMW4YOPZOj16LO8ya18pS4Xw0HZUMiuSf1W0Vhzqck1X3b/e/+u4bbzTs27VNq2puBlgXv9/67Yq1XnR5cCKnyOatEoXWISiEdR1bxWOzBWZIavma9z+tD7/9NnidnCqwjh4+D6DJzeN7A1bYz0klFq1DCDLCfdo1FJxwdnHa4mort2zZfuhLAKJOx/bturRxARZXomqA1iEkMm4bsUix4vtRORgyBFbhyTzGsiu7AcCCUEUGrUOFMgUT4OXERWZO0RsWcaaOC9A6VBxSHJlbpy8MBhAIAINKp28XoXWIGQpv4bPkCy+OqmIJAMzGbYaDo+MOsWGZYCycp7xJYArG2J3cUOjaIRYkEwAxA6G37yW+bUXbG9RrHkqjdYgBybF00Y7zWVcN6M9CrTCVfDbGcdXXoWodoscSp366AACWvtrv3VrB7zQZ5jSMULbPUILuO0SNFGdvXAB0GXUN8EmfUyLoeieJL7b1aa38jvCqp09pHt6COkvXzfh6US1s/2XwrxKmYoJtW02HjdYhBIH37uMb3gKfeZMkGrAwHe+v5UHGF2ScnZU5UPm2QusQCk+gmyzxDnDbDV7G4Ol8D/BusKtILDs2tF6YByjaiPIJtBeGK2v6/02TUbAUmADOwAQmAEvyRbDpjSkYb3dmwkDrEBSOkR39lxHh4D7UbgcOgKN5IxzZFfjxaSkOz6t/YxIUtvopVBTCDi/qFm3bw6hGY3h2BtEBCzyjK1nKdx103CF6DgLNULTfiunOBg4H/CYJTMFou9MbDlqHqDkO1EXRs/pu5ADHAn6T+Ler21+neoAa0DpEDeEEcGrneQ4AVrVc4GTAcaTkMzEuZBhbBbQOUeOBmDq95IUUVsDTJsjwTsdw9YMO0DrEBLMuhxlFsHHF5ztvisd67MSjdYgaA+aQ+/KKGovKb85GwMMqCFMx0bbD7yLMkyZKagJNB5a8WA3UCHDHyTgzu1L/pGgrdNwhWhhqAz/CY1mWZRXiKFA7wDclPj0xrI76AWpA1w4x0ADYC1fxzoAHgfoBvsfwYVIEqAFdO0QNQ2uGoycYASCeewi8lf8vSv7d/Jw+SVLQyZFKBWFo0gRHdkACIOzbj3ot/RemxHxrdGYSBKgBrUPUMKtKD2AWJAATM4Ee1fzecjKSJUANaB1i4SfAmztdpmm6D/4LuN3/lwib1nbqlgwBakDrED2C+gxG/shdhmHsv/Usbunv/5F8iZkYj+iP26tY9Mgiaoi9tv7Y9s4DGx9edBY1/+V/w1cyCqdjaHIEHaB1iAFu5awc/e2F6QBw1fQr/c9QSbFmT78WSRGgBrQOsSCslt+8tXjHhaxW/e90B5yvnJIsAWpA6xCSYA/kCmncc0+hKdIB6d8GEqfmVOmXPD20pEmoTQTfKoyTKdOy0skM9DyNxJLTI2olSdABWodQsJPng/nADE5EzAh0vxmmYJzdeYgArUNQLHz4dYhnp1iQg1clP7CgSe8kKuTkSalNWLF0AyXm0pj05GkrtA6hiP7QZYAMTMVIu3MQCXpkEQoLJgGs5H8hIkqlYlHENnzVpUuyBKgBrUMIGFhVpJV+zyIGFmB7J2K+Qkg+E+ORDCuoS9A6BIWQ/3KXQrfhMlyZWZmVq2emZWSW3F2LWNluJLFDOZcX0JNxcTpLngA1oHUISVqvnielaZqe/CN5+fkmI7iqVa9WvW7takV33fKpKYjRhPpThVlcqFKs2j+oadIEqAGtQ0hE+9a+L2V+fu6PZ84cXnERaXXr5jSs4QIAKTljgBSLhmcOeKvepXNUp2BC8gSoAa1DSCgfJgdAAIEzXqlSrRYAIHOPHvxx6znKzGnUvE4GB6TkwNalZsvxL3Y1BQNInJhTPYkC1IDWITTc5zhcFGlBICZq1GgDIPfwnu3LTFfDdm1qcuD47oP7vu79y4fGkuSwjE/O/qxmMnUktQ4Rw0qGmgQiiGrV2gF5R7evneNu0rb70lNHeN7CHq/s/S2zBE+yADWgdYieog4kQRLPatZsqDyw+bMZe/MuckarOyw/9LrwuPcuatYzudqKJEutejAmDE5kmbzprb/7Wc0jADG+6ccTw4+6MRdj0pIoQA1oHeICY8IAefA19kKCJNu+vtr4r/BBcgWoAd1YhCT8k/Lc+UcNVslLAIzvl0ya3H5j105qnHAVPlqHoHCMzgqzBiW2+bOOQ2SBBTBm5O9r+yLGB1okpSxah6AwauB/hbSfr6L1I8agYneI4VRl77CkClADWodQsDAPLwEYVR3o+3rNj0ObJFWAGtA6hCT83zcjeenLZLn/lWwBakDrEE98TjmT7mPza96cfOO2pEtwciDx8fmRNZIs6ACtQ4Lg+BBj7U5EVOnWxB/Jdy9p0SMJCzf5UpwMSMzBWFfytRVah0RABk3DyKQLOkDrkBAkvtp4Q3vlD77yQxImWX0I05Ph4Cs/aB3iDxkXZroGJ2NboXVIABIrDg/OSboANaB1SAxJt4K6BK1D3JHi6Ee1bkrOkk3KRKuNxKL8kdWSMOgArUMCEJiCMXYnIkq0DvFGsh3LruqepAWbnKlWGYk5GGckZ1uhdQhFpAMEMuQ0jEjKoAO0DiGJ9L5KfLGlZ7tkDFADejVUCCRf1LBdRAElwjRMsH0FtSxdqYWdGq1DUAhbXe0iaS/IOD8rfZDtbUW0lZPWIQTpkZWQFMt/GN3A9se21xx1lUhMTLpvSQvz77QOIYh4XtL+g6+I4fHPfd9I3183zIdFtA5xRYoj8+r0tb+DngaXBQCMKl+AlZ4e7t9pHeKK5Asvjqpic1vBgL/nCgLArHpTfyPwQvVwe8Nah7iiSoD60n5Wm/4M64m7wx73ah3iieTbVrS9wf62orjHY7r39L2ACX8Of2ikdYgnks/GOG77uKJ4oGm5zw45je7voahyoDBiavaL7CDIMBUKUEtBw3bhyjnc4gCkyRizQo2TtA5xROKLbX1aKxKgJo5JK1FpXrYpAFjcgITgIY7sUyPpDoEwFRNVWUEt8fiHHLPamAYAKfY/1bVJx4f3iOCp0zrEDzLOzsocoEhbYYrXJhvyjVtMA4DF13T444+tc1+5+lMe1AetQwgiKCCJT48NrWcpUaRe46OfC/Pxuy0DAAnvXedfO7T4wL8u3pnPg40zlEi7ylw0I/jyFNsD1MWYrq/GMmvcc0X3V2LWrn73Q+K+YYcWBW3MtA5BYWhbN9zaX4rD8+vfqESJmsbBoR667j1YjIiIUKnDeJjcg3Y4GnRBj447BIVjYNg/GckXeEZXUiDoAGmcH3QczeYbJMAAGBg8mGDAwBfICWq31iEEYU4FAiQwBaPtTi4AEPeM2I6cj2sWlqjJGRjgdS1YlhN8fk3rEIKwxwnEt6xuf50KB19ZxuvLMgqq/PKH4m6COHP/I5aA17XpdrwQvPrSOsQLyWdhPFOhrQBOoQBbt15+vRsgy7WuX+5zY4Mv29M6xAkyvNOhxr6iDC2uqWrxS6cDi9z2ABlzR3uf/3WItk/rECek+HznTa2UeGxbYOLEMm+RwMuP4P2JkmkdKgTCFExU+LHtX79Q7aNeXhGia6x1iA9knJlTqb8SbQVQ7mEhKSa/0HxJY7hC/Z3WIT5IsfTEpDpqdCRRdjxkiW8eT/9n7cMuQFbNCvZ3KrR1ToCpfPAVw59wcUBmTt26Der/A94g39S1Q1yQ/Lv5OX0U/XERN2ve4PZyAOJUQx2VTDySz7NGZyrTVpSGwXjT52WwW66mz8kGGZiqRoA6RnTtEA+IbVrbqZsKAeqAKSz5Dx13SDySz8R4KNpWAAh76kU3FnGAjEJVAtQxonWIAxJr9vRvrkSAOkYckAUVSNp9RcugdYgdEqfmVO3niKJ0Qh4SShin7EosOT28VpLuFVcarUNwiLGQPigdoI4MrUNQiH1/OpQPku9f0KS3M0rSEZlIHBIzNoV6yk5iHo1Jd0RboXUIBQt1m8nAVIyyO5lxQusQK4QNX3XprHKAOgK0DrEiMRPjYcV+IRXQOsQIGRens6FOCFADWoeYkVi1f2BTJwSoAa1DHHBKgBrQOsSKFCfmVHdGgBrQOsQK4ZOzI2o6I+gArUOsMHzolAA1oHWIEcn3Lm7Wyzml6JiM2IPEXIxxO6at0DrEBBmY4pgANaB1CEnQ9Q6EdRuu7UjOCFADqq2kNkHGpYrXkhB228qCJ0HyGQocfBVH1NLBwOXNmEiICHZmShAcE9IDz06RkT9DDHZKgBpQSgdi598qlOMbFSlAbO42DL7a9uhvdpDPpFj53dAmticxjiikA0Pl/36LM5OLNkplF+4+gf52pwmgoD/+Dx0UoAbU6kqaeBhYUGgQAImPT2BwJ/t/eUGWv0hxbG72LUoVYayolBeBgdnY/jkkAIapwGhVdn/3j8Ti87dWd07QAWrpwKxag4GpIEDyQ8tQZ5Daa4w4pmCs3YmId5aUYgKwKNcgEObnYpDaU0OS717SsodqJRgbSmWGo08bHFkKCY6pwMTYr5hIJOZgrEtpYyNGKR2YZYwFpgEW2/glrlb7l0cGTcWtTgo6QDEdwDAqHZ9+LyRmmRhpmCqXtcRXm25or8jBV/FCrdxweVVvnF0A4Z2FtNGKJa4MhOmYoPbQJ3IUK3HCJGAG+PKd6HWV/UGHYCk1Lsx0OSpADSinA8fAelizFdOB8Wr/8iRWHB6So7SxUaBYdphVfTg8C7EYtQerHXQAPsR4RwWoAeV0ADAeWPPJUQzKVnoIJ8XRubVuVrD4YkO1/Ah0vwYb3gAm2J2S4Egsyh9ZVWljo0E1HWBhPI7NQ9ue6iXNF+G8ADWgoA4Mt2ZJYJRL7aAD27Hsqu7qlV6sKJchLpv0gzQUDzpIzME44bi2QsFClxgNqcw59/4hQzovQA2oqAPDRmCC2hsmSHzxba+2ShsbHYrlSJIUF6ej7hC1gw6EaaqHyaJDMR0443j/IIbXVKQjKTH7m/K3nYxzs9IHObCtUGnpLADaUostfATGA3ZqKolfutGEfVXKhx6l+OyH0Q1U3qA+WpTSgdhDq7POc9zXxsaiJl7KCLfflHyoepgsShTTIVueh+zzso2Vg+SvVB1YB7BYURr8zEpIcWRe3RtVa2fjglI6ALduuVB72OMuGycKGb7Z/16Xm7pWA0zuPxmSL7w4qooT2wq1dOCYMLYgk8PeaeNKh37YPKf1tTd2SIP/A2wEpmCMnSlMHOV0sGDaOGtLPEsWcmbjEE66L16UZ3KPbvyo9bXd27j9PMEt+bYV7a53ZFtRXoeqSLM3RcLuStglIay8vGMHNsxt1e7W9HLDSclnYyx3ZFtRRgeCNfV4vgPH05Gw2wUJcM+JE3t3n16y7n9lPibDnIoRTgw6oHztwFv2OevIajBcyPjSAhiTQGamPHu9p2xjIcXa7YpPqURPaR0YWNur7U6S3WRLgMidlV3v6t6dc94sO3lCmIaJTtriw5dyfYeL8Dozp2FCwiQYGTWy29zQqwWAc2VaBTLOzsoc4NC2orwODNyZ9WCYEDJrues269anvQGSZJQtDCk+PTaunjM7korFHVSA2Mnm7W/sXAUwORcWeLl6YIpDA9SA1qEcDI/VuBKwGC8qmgJPqY+lODy/fh9nBh2gdSgHo86wwIsbA46u9Ur1EyRf4BldyalthdahHEyyyzebUR+UGlM6OEANaB38UKohYFSmctiyusN1jm0rnJuxeFF6qzCJWRgH562gLkHrEAlkeKdjmFODDtA6RIbEmp03t3LaY9s+ODdniWGqw/YVLYPWIQJInJldub+D2wqtQyRILD05rI7l4DJzcNbiD8MUjLc7EQlF6xA+kh+cn9Pb0UXm5LzFG4n51phM5wYdoHWIABIOO/jKD1qHsCG28Ytruqn9LHGsaB3CpjhAbXcyEorWIVzIKHR2gBrQOoSPxJo9/Zs7OEANaB0iYQomOjlADWgdwobEqTlVnXXwlR8cnr34IbHk9PBajg46QOsQkpJHdhk+dHiAGtA6hIIYIwCQfP/CJr0cX1xOz1+MENtzjBWd6zmPxqQ7va3QOgRHYt63kAAZzg9QA1qHkBgCAAgb1nXt7OwANaB1CElRT1JiBsY7PEANaB3Cg4yCGWyIwwPUgNYhPCRW7x/Y1OEBakDrEC4fOj9ADWgdwkKKEx/VcHyAGtA6hAXhk7PD1T4fPk5oHcIgNQLUgNYhHCTfu7iZ2ufDx4tUyGOsED7CmLRUaCu0DmFgYCpGpUDQAVqHcODrNlzX0aH7ipbNqt0JSAZmOPPgKz9oHUIhcX6GSIUANaB1CAFBmlj53eDGKRCgBvRWYUGRMFCpKj7EuFQIUANah2BYAkd3LTt8bn6DQalSi2odAkEQ21+cng8ALblTd6gvS4pYHwXE/tH27ax7Bk1ogA2tVgvnL30BtA4BsfhLD4oXdv17/KD8xvd833ttavigdfCPJVY/mr7sV9VR8N/Tk/49WU46K1KhN+m0vkO87hnD0/hnTw832HI+Ao9tf+dvv4/XvuQqRzAcpgPFq6zZxuWd7iCXxDnc0Mk0nnjnncfddieSAkdGLSElrIC/BRZuI+AsHYihMC7Nn5mxHKNhGYStGAfIFjd9+s21F+NSPcg0Fp0PxAL/+wJZleGKPXGO0oHYkdsO8Hi0FyTOoDkYGWcXYDC45PUxKj7PYHHP9e+5ovGB2JnlAT+UVVYdPZEX4KrMqt43zH/QUTpI8e1nPFPGpcGwitp4V4uxOZJLALlxSmTejKdbRVc9WOcC/1l+b3E80O+AmWHfZUfpALhx+//OxaFOt6o898RBEKPMVQwcDMexqOf5eLRDxo1fRpc+Rtl3xLew/KUu4f9CBWOgSlyu0xezHuYoauYt8d3Shl1QOS4X5tFOlTMK9pfBK5xwDXScDgRvHPJErOu1a2eM9rgYMUjCi+aEDDMOtU7J/gBRwRIfKHecDmXPp4nyGpZ4asB9LdrDBJjg/32t9sN+jlCM7tI2F09QdFTSP8Lq/+TpXm97DcMQp359F96u5eSN6i/hwNohPnD6k/H0Hc8Pae7dOCu3yvv9HXtWYim0DgFghKdufm7BLgCVHni8QWrYoHUICCN5/fzd6783mnerjRSxQUkdiOLTH4wRJizeogUAWBXQp1cD9XSQTAEVihCQkhjjqSKDgjpIjhOHqWrdTCWc4KkwnPBBNR0k3/artfmisO9iI26T1ZqwUUwHyQ/0Pgl4wY1U6b0phWI6EJ456fZc1/aHIXanJDVRTAeBL+F59MWi/9RUOGp1lQj5F4Db4EmNB2TVQ63agYEIRmXEbVmiJjLUqh1My0NAIS6mxEMNCqKSDhJPtux9DNbwa1r9DqbdqUlJ1GosNu4DGO0CNqfKI9OKoZIODPd0xd/P8gfr5F6rBxa2oJQONHIk3jsrnqoOpMj2GqqhVKkzU56xQMdRmBJLjxREpdoBEMzFAAMuG22QEvC3TJUswPFzWmrpoAAU4IYTM4r/3+4UJhKtQ2mI7VyZJpmZPVSUef/MbA5WeG17R/ugdSiNFJ8+CADY2qZUb1aKGfcCwFPtHb0tkLObwmhwwyW4GwtLRz4YZsLN3U4Pn2sdykKQFln4qFTRSL53DbxSOj06pnXwB0ms/5b5TKtKLCwwyOkyaB38kw3DXOR78wVmwxJO7kQWoXUoD0fPZmap1kKy7V8g81pHDyqKs64pC0PTvuDrt1xuLSTme9ExBQ68UTCDKvwER8GwLrcWZGAOMKyG3alKPArqoEB/Lb9XHQ/mXiocwsZ1wHCv3clKPArqoACFRqnWQmIuoVOzC3YnK/GopQNDxuxVK+srkKrRMKySSBQZci5wawqcyK5AwZdGdO5xfbrdiQBHj5zLrYXE+k1wD1GurBKSccWQlgKL6hlq3AyxfnNRa0GYDVx7tXplFX+UyyIXaiRpFIQsai3I8CwARtmdoApBjbKPBWklZCjSvUlJayHxxXZkDbQ7nxVC8uvABbPib0RhpX4QX29iEiDMAro3TYFhpgN0uPj7WSeFYPHucTCMLG4tyMhfDIxJjYM0k1wHQsGcx25/anWe4DApjnUExw3NPZgLAYk1e1H9lmQvqfBI+tVQhth/dOOsq66/sZ2BOEY0mZnRf4/4emNHSZgJ9G6QGmu7K0wHmZDKlgwvoaDgxMGvPmjfo+eVcVysxDDy79xa2NFynf0EGFNRxWQzFaUDJWpFeg0Bxqy8vGP7VuS0HZYbt/qB49qrdmDukxyrDqF23+h3Fk8qKkgHydfNz0xEgfK8E8IiMObx5B7asTk/fodhmWkDd4gN33TCdOCmK1Jka6IK0oGhyUAjAfEBMnKnSYBBgmVWrtLYtSFu0+MMI17i1oJOx5alTltRYTrQFVck5spWFeKSkF65RuNOfW5e+1bcdODo1nYr5v5+9VHk9E6NcUXF9R1YPEeBlyFx3oJMy6hZt0Pv6+sC+XFMsWkM3iq+OfAx0L9KirQVFTeySMxWsgRk1ah8RZvuPVsA0nTH80fMMPx54O+fp8p8BeCAuAPxmzrf2N4ASR7nzYI5unTYJP5ponkKLJIsJsl1YKj6ypU1ADPeLgBgpjFkEys0MDg9VdqK5Need65hSRjxzAcvXr3LMQxSEEYWv88ufeJYkl4HSIrvCgkTF+EBADBc01FaVuvrikvJi4tw+LxmkjcWiL/QVzSqfiYHDACzxB2nq5++kxe3FfUaVT9d19n1Q/LrEF8ERo0AeJFkAg/eB7iKti0T+Pn9gHBAfRoErUNZSk+u+B5zLpzfodQ6lIN8HwQL/MKRaB3KwcJ74Ugc3RJqIkXroPFB66DxQeug8UHroPFB66DxQeug8UHroPFB66DxQeug8UHroPFB66DxQeug8UHroPEhJSe4iQgExuL27EfcL2gXqaeDlMaluyZlHHahl1JcuiBZyX12WqrpYDHOceHId+fyMq5o2DCNw2Kx3T9JgqPgx4On843s+o0zjZgvaCuppQNJgS1zPl93tuhlg+5Db6kJM4YyIMvAd7NXfnms6GV2t8EDGsHiSdtmpJQOlhBfPjsfqN3pyuzM/GO7tkydWuOeBxoQRft7toSx9/n3PKh6fbNalTzH929auNB1+y+uRvI+tUW+WDTvDTIpWTFpGe4kb6CPvXTmXqDF5G/yit84s+DODFT5N0WbZS/JP7lR6/E1ucVvXFj5SB3gSW/ANEhJ3bGbLLtLKhAppINJyxuj7n89RESm1+O1iIgOPQaMPBdYoWB4aff1cD9zhojI8nq8JhHRmVeqoOveQBfUOlQgQXUw6UPgZ6eJvJYsekdK0yT68mq02R+ND176qiZu2UPkNWXJFU0v0YFBqLYywAVV1yGJe8GRYYr3J+CN/1Y3ySjp6TEmhDS7fTl2W+/9RsSnE5jGV31O/fqTZiYZlwaaTBjkbbzgN7l9PzNMu3McDamig2ksuQ3T7rakUbrXzw0zY+q9hwaeExFuZGYZewfnP/M8rDIXZC4Lz75kDtoWuWAKkCI6SOPwOLw5xvQTJDIsvD5+110RPlNDwjvu5K9+K/0cwSzIevhPBRPzhQKnOEVKiujAcN+Zx+40DX/3XFh4s+2MD1hEtbvEbzcMfgF+Q06M05MjN/8xGXexTg0dTDZtUddnESAaIMzM1/Dk+Ug2OrTEphdrvA7Lf5XCJP5W5YXtkbY/CpASOpDhfRbPCTNQe2BYPW8/+J9IzrpimIyn63sDRZuEWf939KIKZ/5FSEroIPHplv43ymAR2Efx9wvhVw+S753a6M4gIV2Be+q9vSP5qoeU0AGYjp8Ga8qF1W7kwVXhN/YSCzAx0wzc+2RmlXsxP/mqh1TQgUT+8kq9QmR1BD4J/4ocizA46FiEYzAWJV/pJl2Co0Bi16EudWSwrDJ0xeqwC4P4yY112yL4BVs2/uY4S7bqwXEzmhLlBoyW2I6uCKFDo+Z7Djf0hCeElfbdyf6Vg1/Qyrp63u5anjKdTYKhdH/CYTowZKDcqawGvkPr4GEmZrmb7ikM93QUAyfQOkTHgNB23mmklf8gTeVNZBymg4mt710o+6uljMVoGPJPa+DNq/LDqx1k1mLkhPgOQ3XMPFsuLUwecCs8m+EwHTKxalWU+bQwOZJ/KfQCF4n33vP3flamLUUTFo7SQaDbu4fSyzXOlP7ButAhJoEHmxeEV49TxvKPQncBGEb3zCtX3fCC9jnqziM7SgfAuM3v2+fXfR/yT3PxUNOw/51mH4W6IOEcJg31/5m6nQeH6QCTyg/urPSG2BG850fCuy/TgCe8GyXTsrEtxE1l2IaqKCxXDzCK/9768cNpOvjLD0drfBU8qEDs0K72jcIdWUg0qr45L4uCCEGiYEtWS7jVrQj8omwrFtc8tqq//njQkw8J63FD2HNYXNbqeGQbgl9w974OdSjJbEgJHZhVqc+5NSHiBHPRP/wrSvTHoqAXlFiMAcm34iEVdAAwBu8E+1iKnVPr9wy/MDiG4gNPkBlQMvL+jaEK9xkDZiwF4Li51bw1InD4R+Jl/LyqGfbd47LlrXvfDtK4WHjr4Ph2yXdud9IlOBqYmf4bPIGAv2bT+PI/9e4NI7J0CcLjeOpEwPkHaRz/Ax5LvsohNXSAwO03rn460K9ZGt5f4E/Vw68cAGF1+fmPDyDAjCUBj5x6qEPyVQ7OeigvMCbtycIU8vj7zJJ0FwZHeEFJ+W3wFJnS32deegktc0lGeE0FSEKBo0GYzT7AhEUus3z1bnL26JuN/hfhKIBZGVMr/+Gvwk+NI8n49yPig6pW8rUVKdJYAIY17HUaNMXgZun6nbyGdddLtRZmR1qzC7PtXPF/zzBR7oKcP3sfPrrGVDj4GJjS9aZTGwsiMuk/wEN5Ps9oFj1S+e11aLwtumc0l2dhxPelLmh5iY6NRfqiJC3F1NGBTJpfC82mERFJr8frMSURnXjGQL/j0T7Bvak9qv41n4ik6fF6vJKI8t+sjTZborug/ZTuG0s+/4d7C502j1FSD1pp3z86De3vuvnK4vPv8r+a8+5Z/tyvWJRZNtPyf/8SGt49pGVG0RueLfPe/Q4PPFfZhjKMS9tUVoe5FyZWeEYqkm+emW+KZq2b1MzIP7530w+od8dDV8R0wcPPTD2HnLZX1srynjiwaT+qjPu/VnZnMnpK62CJj98ZcMG53UtJVT0rlpzxlLzObDmgjedCDLvHkVWZr593vKDkdVqjfp3ZORs2C5PpY7PiMGFWNpBS8HlhEo6PwsfilTKswtNnCjzurOrZBuUVxHjvLJaZRZ7c0/mFIqtqLTfy88mOIQW5e6bFfhUk3ZMAmkRSrsuTdHOyUUBFc9MsbpMKcb9gFMSlgdK1g8YH53YbNVGgddD4oHXQ+KB10PigddD4oHXQ+KB10PigddD4oHXQ+KB10PigddD4oHXQ+KB10PigddD4oHXQ+KB10Pjw/yFtjQP2hx5oAAAAJXRFWHRkYXRlOmNyZWF0ZQAyMDI0LTAxLTIzVDEyOjQ2OjAxKzAwOjAwnMgN0wAAACV0RVh0ZGF0ZTptb2RpZnkAMjAyNC0wMS0yM1QxMjo0NjowMSswMDowMO2VtW8AAAAodEVYdGRhdGU6dGltZXN0YW1wADIwMjQtMDEtMjNUMTI6NDY6MjErMDA6MDD4pZPNAAAALHRFWHRTb2Z0d2FyZQBZYW5kZXguRGlzayAoaHR0cDovL2Rpc2sueWFuZGV4LnJ1KVzqC+YAAAAASUVORK5CYII="
    }
   },
   "cell_type": "markdown",
   "id": "42f9b006",
   "metadata": {},
   "source": [
    "## EX. 1: Inverted Pendulum\n",
    "\n",
    "The system in this example consists of an inverted pendulum mounted to a motorized cart. The inverted pendulum system is an example commonly found in control system textbooks and research literature. Its popularity derives in part from the fact that it is unstable without control, that is, the pendulum will simply fall over if the cart isn't moved to balance it. Additionally, the dynamics of the system are nonlinear. The objective of the control system is to balance the inverted pendulum by applying a force to the cart that the pendulum is attached to. A real-world example that relates directly to this inverted pendulum system is the attitude control of a booster rocket at takeoff.\n",
    "\n",
    "![2024-01-23_13-32-20.png](attachment:2024-01-23_13-32-20.png)\n",
    "\n",
    "Let us consider the system with the following system parameters\n",
    "    \n",
    "    (M)       mass of the cart                         0.5 kg\n",
    "    \n",
    "    (m)       mass of the pendulum                     0.2 kg\n",
    "    \n",
    "    (l)       length to pendulum center of mass        0.3 m\n",
    "    \n",
    "    (b)       coefficient of friction for cart         0.1 N/m/sec\n",
    "    \n",
    "    (I)       mass moment of inertia of the pendulum   0.006 kg.m^2\n",
    "    \n",
    "    (F)       force applied to the cart\n",
    "    \n",
    "    (y)       cart position coordinate\n",
    "    \n",
    "    (theta)   angle between the pendulum and the vertical axis\n",
    "\n",
    "## TODO\n",
    "\n",
    "Inverted pendulum on the cart can be modeled as follows\n",
    "\n",
    "$$(M+m)\\ddot{y} + b\\dot{y} + ml\\ddot{\\theta}\\cos\\theta -ml\\dot\\theta^2\\sin(\\theta) = F$$\n",
    "\n",
    "$$ml\\cos(\\theta)\\ddot{y} + (I+ml^2)\\ddot{\\theta} - mgl\\sin\\theta = 0$$    \n",
    "\n",
    "Let $y_1 = \\dot{y}$ and $\\theta_1 = \\dot{\\theta}$\n",
    "\n",
    "1) Show that linearalised model have the following form\n",
    "\n",
    "$$\\dot x = Ax + Bu$$\n",
    "\n",
    "where state vector $x = (y,y_1,\\theta,\\theta_1)$, control vector $u=F$. \n",
    "\n",
    "$$\\left[\\begin{array}{c}\\dot{y} \\\\ \\dot{y1} \\\\ \\dot{\\theta} \\\\ \\dot{\\theta_1}\\end{array}\\right]=\n",
    "\\left[\\begin{array}{cccc}0 & 1 & 0 & 0 \\\\\n",
    "0 & \\frac{-\\left(I+m l^2\\right) b}{I(M+m)+M m l^2}& \\frac{-g m^2 l^2}{I(M+m)+M m l^2} & 0 \\\\\n",
    "0 & 0 & 0 & 1 \\\\\n",
    "0 & \\frac{m l b}{I(M+m)+M m l^2} & \\frac{m g l(M+m)}{I(M+m)+M m l^2} & 0\\end{array}\\right]\n",
    "\\left[\\begin{array}{c}y \\\\ y_1\\\\ \\theta \\\\ \\theta_1\\end{array}\\right]+\n",
    "\\left[\\begin{array}{c}0 \\\\ \\frac{I+m l^2}{I(M+m)+M m l^2} \\\\ 0\\\\ \\frac{-m l}{I(M+m)+M m l^2}\\end{array}\\right] u$$\n",
    "\n",
    "\n",
    "## Let us design a controller such $\\lim_{t->+\\infty}\\theta - 1 = 0$ and $\\lim_{t->+\\infty}y = 0.$\n",
    "\n",
    "2) Check for the code below, which design a PID controller to ensure that angular velocity tracks a refence input equal to 1 rad. Is it a good controller to achive desirable behavior of the system?\n",
    "\n",
    "3) Design state feedback controller for a given specification:\n",
    "    \n",
    "    3.1 Is the system controllable. Why is it important?\n",
    "    \n",
    "    3.2 Place arbitraty eignvalues in the system to stabilize it. \n",
    "\n",
    "    3.2 Use LQR regulator to stabilize the system around zero. Analyze how the different choice of weight matrices Q and R affects the closed-loop system behavior.\n",
    "    \n",
    "    The fucntions control.matlab.place and control.lqr from control library \n",
    "    https://python-control.readthedocs.io/en/0.9.4/\n",
    "    \n",
    "    might be useful for you.\n",
    "    \n",
    "    3.3 Add integral term to ensure that angular velocity tracks a refence input equal to 1 rad with zero steady-state error.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 147,
   "id": "5c8a9d97",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Answer to EX1 q1 using symbolic calculs\n",
    "# Define symbolic variables\n",
    "M, m, b, l, I, g, F = sp.symbols('M m b l I g F')\n",
    "y, y1, theta, theta1, doty1, dottheta1 = sp.symbols('y y1 theta theta1 doty1 dottheta1')\n",
    "\n",
    "# Define the differential equations of the system\n",
    "eq1 = (M+m)*doty1 + b*y1 + m*l*dottheta1*sp.cos(theta) - m*l*theta1**2*sp.sin(theta) - F\n",
    "eq2 = m*l*sp.cos(theta)*doty1 + (I+m*l**2)*dottheta1 - m*g*l*sp.sin(theta)\n",
    "\n",
    "# Solve for the first derivative of theta1 (angular velocity)\n",
    "dottheta1_sol = sp.solve(eq2, dottheta1)[0]\n",
    "\n",
    "# Solve for the first derivative of y1 (linear velocity)\n",
    "doty1_sol = sp.simplify(sp.solve(eq1.subs(dottheta1, dottheta1_sol), doty1)[0])\n",
    "dottheta1_sol = sp.simplify(dottheta1_sol.subs(doty1,doty1_sol))\n",
    "\n",
    "# Define the state-space representation of the system dynamics\n",
    "f1 = y1\n",
    "f2 = doty1_sol \n",
    "f3 = theta1\n",
    "f4 = dottheta1_sol\n",
    "f = sp.Matrix([f1, f2, f3, f4])\n",
    "\n",
    "# Define state and control variables\n",
    "variables_x = sp.Matrix([y,y1,theta,theta1])\n",
    "variables_u = sp.Matrix([F])\n",
    "\n",
    "# Compute the Jacobian matrices of the system\n",
    "jacobian_A = sp.simplify(f.jacobian(variables_x).subs([(theta,0), (theta1,0)]))\n",
    "jacobian_B = sp.simplify(f.jacobian(variables_u).subs([(theta,0), (theta1,0)]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "id": "30b8e0d8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0 & 1 & 0 & 0\\\\0 & - \\frac{b \\left(I + l^{2} m\\right)}{I M + I m + M l^{2} m} & - \\frac{g l^{2} m^{2}}{I M + I m + M l^{2} m} & 0\\\\0 & 0 & 0 & 1\\\\0 & \\frac{b l m}{I M + I m + M l^{2} m} & \\frac{g l m \\left(M + m\\right)}{I M + I m + M l^{2} m} & 0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[0,                                      1,                                    0, 0],\n",
       "[0, -b*(I + l**2*m)/(I*M + I*m + M*l**2*m),  -g*l**2*m**2/(I*M + I*m + M*l**2*m), 0],\n",
       "[0,                                      0,                                    0, 1],\n",
       "[0,           b*l*m/(I*M + I*m + M*l**2*m), g*l*m*(M + m)/(I*M + I*m + M*l**2*m), 0]])"
      ]
     },
     "execution_count": 102,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "jacobian_A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "id": "99695b56",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0\\\\\\frac{I + l^{2} m}{I M + I m + M l^{2} m}\\\\0\\\\- \\frac{l m}{I M + I m + M l^{2} m}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[                                  0],\n",
       "[(I + l**2*m)/(I*M + I*m + M*l**2*m)],\n",
       "[                                  0],\n",
       "[        -l*m/(I*M + I*m + M*l**2*m)]])"
      ]
     },
     "execution_count": 101,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "jacobian_B "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "id": "39e5f873",
   "metadata": {},
   "outputs": [],
   "source": [
    "def StateSpace(x, t, A, B, u):\n",
    "    n = A.shape[0]\n",
    "    B.reshape(n,1)\n",
    "    return np.dot(A,x) + np.dot(B,[u])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 150,
   "id": "ba84a7ee",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Position ')"
      ]
     },
     "execution_count": 150,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtwAAAE9CAYAAAA4WbXTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABaa0lEQVR4nO3deXxU1f3/8dcBNCwCGhFF3IlailRqAdF+29rWNrgg7hWqRmvrAlSxtVVroS183atiRaq1LmAFl7riQlx+Wr+tCyCigEsbFAuKRA0YBA0C5/fHmYRJMpPkJPfOPZN5Px+PeczMncmd90wuh0/unMVYaxERERERkXh0SDqAiIiIiEh7poJbRERERCRGKrhFRERERGKkgltEREREJEYquEVEREREYqSCW0REREQkRp2SDuCre/fudt999006RiNVVVUUFxcnHaMR5fKjXH5CzQXhZnvllVc+ttbukHSOXFGb7S/UbMrlR7n8hJorsjbbWptXl549e9oQHXnkkUlHyEi5/CiXn1BzWRtuNmC+DaAtzdVFbba/ULMplx/l8hNqrqjabGPzbOGbgQMH2kWLFiUdo5GKigpKSkqSjtGIcvlRLj+h5oJwsxljXrHWDk46R66ozfYXajbl8qNcfkLNFVWbrT7cIiIiIiIxyruCe/ny5UlHyOj8889POkJGyuVHufyEmgvCzlZI1Gb7CzWbcvlRLj+h5opK3hXcIiIiIiL5RAW3iIiIiEiM8q7gDnHKGIBRo0YlHSEj5fKjXH5CzQVhZyskarP9hZpNufwol59Qc0Ul72YpGTx4sJ0/f37SMUREWqXQZilRmy0i+axgZylZtmxZ0hEyKisrSzpCRsrlR7n8hJoLws5WSNRm+ws1m3L5US4/oeaKSt4V3Bs3bkw6QkZVVVVJR8hIufwol59Qc0HY2QqJ2mx/oWZTLj/K5SfIXBH2Asm7grvFNmyACy+EK66I9AMTERERkXZuwwY45JDIdtcpsj3lSFFRUcueePPNcNVV7vaAATBiRHyhgH79+sW6/9ZSLj/K5SfUXBB2tkLS4jY7x0I+PkLNplx+lMtPcLmuuAKefz6y3bXfQZPf+x48+6y7fcopMGNGvMFERFpAgyZFRAK3ZAl8/evw5ZcYKMxBk5WVlc0/yVpYsGDL/Xnz4guUMnXq1NhfozWUy49y+Qk1F4SdrZC0qM1OQMjHR6jZlMuPcvkJJtemTXDGGfDll3DWWZHtNu8K7urq6uaftHIlfPopdO0KnTrB22/DunWx5iovL491/62lXH6Uy0+ouSDsbIWkRW12AkI+PkLNplx+lMtPMLmuvx5efhn69oUrr4xst3lXcLfIm2+660GDoH9/d8Z70aJEI4mIiIhIwCoq4Le/dbdvvhl69oxs1+2z4H7rLXfdv78bMAnw738nl0dEREREwrV5M/zsZ/D55/DjH8MRR0S6+7wbNDlo0CC7cOHCpp908cVudOnkybB+PVx+Ofz+9/C738WWq6qqKsgljJXLj3L5CTUXhJut0AZNtqjNTkCoxweEm025/CiXn8Rz/eUvrs/2DjvAG29Ar15AAa80WVNT0/yTVq501336wJ57utsxr3ZWUVER6/5bS7n8KJefUHNB2NkKSYva7ASEfHyEmk25/CiXn0RzrVgBF1zgbt9wQ12xHaW8K7hX1hbTTfngA3e9885bCu53340vFDB58uRY999ayuVHufyEmgvCzlZIWtRmJyDk4yPUbMrlR7n8JJbLWjj7bFi7FkaOhBNPjOVl8m7hmxZJL7i7dXO3Yy64RURERCTPzJwJjz3mBkhOmwbGxPIy7bPgTu9S0qPHlm3WxvZBioiIiEgeWbUKzjvP3b7mGneiNiZ516Wkd+/eTT+hpgaqqtz82716QefOsO22bgLzqqrYco0dOza2fbeFcvlRLj+h5oKwsxWSZtvshIR8fISaTbn8KJefnOeyFs45Bz75BH7wA/jJT2J9ubybpaTZZYLffx922QV22mnLme7+/d1UgYsXb5kmUEQkAYU2S4mWdheRIM2aBaNHQ/furj7cbbeMTyvYWUqaHcVaexZ7++23bOvTx13HOHhnxIgRse27LZTLj3L5CTUXhJ2tkIQ6I0LIx0eo2ZTLj3L5yWmuDz+EcePc7WuuyVpsRynvCu5mffKJu06fy3Gnndz1hx/mPo+IiIiIhKG2K0lVFfzwh/DTn+bkZdtfwZ3QGW4RERERCdysWfDQQ64ryS235GwyjbwruLvVTvOXTUJnuIcMGRLbvttCufwol59Qc0HY2QpJs212QkI+PkLNplx+lMtPTnJ9+CH8/Ofu9rXX5qQrSa32N2jyyivhoovcikFXX+223XknnHoqjBrl5lsUEUmIBk2KiCTAWjjmGHj4YSgthSeeaNHZ7YIdNNnsqmW1Z7gzdSmJ8Qz3pEmTYtt3WyiXH+XyE2ouCDtbIQl1pcmQj49QsymXH+XyE3uumTNdsd2jR067ktTKu4J73bp1TT+htg93epeSHXd01zEW3PPmzYtt322hXH6Uy0+ouSDsbIWk2TY7ISEfH6FmUy4/yuUn1lwrV9bvSrLrrvG9VhZ5V3A3a/Vqd73ddlu27bCDu649+y0iIiIi7Z+1cPbZrj4cPjz2BW6yaX8Fd3W1u+7Zc8u22rPdn3wCmzfnPpOIiIiI5N5dd8EjjyTWlaRW+xs0OWwYvPwyvPACHHTQlu09e7pivKqq/tlvEZEc0qBJEZEcWbnSrTC+ejXcemurzm4X7KDJ6toz2Nmf4K579Ki/vVcvdx1Tt5I5c+bEst+2Ui4/yuUn1FwQdrZC0mybnZCQj49QsymXH+XyE3kua+Gss1yxfdhhcPrp0e7fU94V3JWVlU0/IVvBXTtryccfRx8KuPHGG2PZb1splx/l8hNqLgg7WyFpts1OSMjHR6jZlMuPcvmJPNf06TB7tuvh8Je/JNaVpFbeFdzNSugMt4iIiIgE4L334Nxz3e0bboBddkk2D+2t4N68Gdaudbe32ab+Y7UFd0xnuEVEREQkYZs3w2mnuXrw2GPh5JOTTgTkYcHdp3YRm0w++8xdd+sGHTvWfyzmLiUTJkyIZb9tpVx+lMtPqLkg7GyFpMk2O0EhHx+hZlMuP8rlJ7Jc118Pzz0HvXvDTTcl3pWkVt4V3EVFRdkfrD273bA7CcTepaSkpCSW/baVcvlRLj+h5oKws+WaMWZXY8yzxpg3jTFLjDHnpbYXG2OeMsb8J3W9XdrPXGyMqTDGvG2MKU3b/g1jzKLUY38ypun/zZpssxMU8vERajbl8qNcfiLJ9cYbcPHF7vYtt2xZhyUAeVdwL1u2LPuD2fpvQ+xnuMvKymLZb1splx/l8hNqLgg7WwI2Ar+01vYHhgFjjTFfBS4CnrHW7g08k7pP6rGTgAHAcGCaMab2a8M/A2cCe6cuw5t64Sbb7ASFfHyEmk25/CiXnzbn+vJLOOUUqKlx0/8ddVQ0wSKSdwV3k5oquDVoUkQKlLV2pbV2Qer2WuBNoC8wEpieetp04OjU7ZHA3dbaGmvtu0AFMNQY0wfoYa190bpFHGak/YyISHL+939hwQLYYw+47rqk0zTSPgvu7t0bP6ZBkyIiGGP2AL4OvAzsaK1dCa4oB3qnntYXWJ72YytS2/qmbjfcLiKSnLlz4dJLXX/tO+7IfOI1YZ2SDuCrR1MfYoJdSkpLS5t/UgKUy49y+Qk1F4SdLSnGmG2A+4Hx1trqJrpfZ3rANrG94euciet2Qrdu3RgxYkTdY9elzjydf/75ddtGjRrF6NGjKSsro6qqCoB+/foxZcoUpk6dSnl5ed1zp0+fTkVFBZMnT67bNnbsWIYPH17vdYYMGcLEiROZNGkS8+bNq9s+e/Zs5syZw7Jly+qeP2HCBEpKSup9pV1aWsq4ceMYP348S5cuBaC4uJjp06czc+ZMZs2aFdt7WrZsGXPmzPF+T+nzGMfxnrZP/T+ay99TS95TaWlpIr+n5t7T5s2b623P1e+pufeUfuzn8vfU3HuqzeX7nm669lqOnDiRXTZt4sG99uI7AwdSMXdupO8pCu1rafc77nArCZ16qpvwPN3KlbDzzrDjjvDhh7HnFBHJJKml3Y0xWwGPAuXW2mtT294GDrHWrkx1F3nOWruvMeZiAGvt5annlQO/B5YBz1prv5LaPir182dle10t7S4isTr3XDfX9le/Cq+8Ap07R7r74Jd2zzYqvsFzTGqUe4Ux5nVjzAHN7Xf58uXZH2zpGe4Y/sgYP3585PuMgnL5US4/oeaCsLPlWmomkVuBN2uL7ZRHgNrTUGXAw2nbTzLGFBlj9sQNjpyb6nay1hgzLLXPU9N+JqMm2+wEhXx8hJpNufwol59W5Xr6aVdsd+oEd94ZebEdpTi7lNSOil9gjOkOvGKMecpa+0bacw5jy0j3A3Gj3w9saqc1NTXZH2yq4N56a9e3e+1a+PRT2HZbj7fSvNqvSEKjXH6Uy0+ouSDsbAn4JnAKsMgYszC17TfAFcC9xpgzgP8CJwBYa5cYY+4F3sC15WOttZtSP3cOcAfQBXgidcmqyTY7QSEfH6FmUy4/yuXHO9eaNa5XA8DvfgcHNHvONlGxFdypMyG1g3HWGmNqR8WnF9wjgRmp0e4vGWO2Ncb0qR3E462pQZMAxcWu4K6qirzgFhEJlbX2n2Tufw3w/Sw/cylwaYbt84H9oksnItIK554LK1bA0KFw0UVJp2lWTmYpaTAqPl22kfBZderUxN8ItQvfZCu4a7uVpDrcR6m4uDjyfUZBufwol59Qc0HY2QpJk212gkI+PkLNplx+lMuPV6777nNdSLp0gRkzXJeSwMU+aDI1Kv4fwKXW2gcaPPYYcHnq7AvGmGeAX1trX2nwvLoR77169frGsGHD6h5LH806fuFCvr9iBS+eeSYH3Xxz49GsS5bA008zcehQXu3tZr8KaYRupvdUK4lR/HpPek96T9G/p0cffTSRQZNJ0aBJEYnU8uXwta+5LiU33ghjxsT6cpENdLfWxnYBtgLKgV9kefxmYFTa/beBPk3tc88997RZHX+8tWDtPfdkfvxHP3KPz5yZfR+tdNddd0W+zygolx/l8hNqLmvDzQbMtzG2u6FdmmyzExTq8WFtuNmUy49y+WlRrk2brP3ud10td8QR1m7eHHuuqNrsOGcpyTYqPt0jwKmp2UqGAZ/aZvpvVzXVHWTdOnfdtWvmx2u/rohhtcn0M2ohUS4/yuUn1FwQdrZC0mSbnaCQj49QsymXH+Xy06Jc11wDzz4LvXvDrbe6hW7yRJydXrKNit8NwFp7E/A4cDhu2eD1wOltesX16911t26ZH4+xD7eIiIiIxGTBArjkEnf79tvduip5JM5ZSpoaFV/7HAuMjexFawvuBM5wi4iIiEgM1q+HH/8YvvwSxo6Fww9POpG3vFtpcuDAgXbRokWZHxwwAN54A15/HQYObPz4jBlQVgYnn+xGt0aooqKCkpKSSPcZBeXyo1x+Qs0F4WZLaqXJpDTZZico1OMDws2mXH6Uy0+TucaOhWnToH9/t5pkly45yxX8SpOJaK5Lic5wi4iIiOSPRx91xfZWW8HMmTkttqOUdwV3k8sEN9elJMY+3OlTjoVEufwol59Qc0HY2QpJqEu7h3x8hJpNufwol5+MuVatgp/8xN2+/HIYNCinmaKUdwV3kxKcpUREREREImKtK7Y/+gi+/30I9A+Flmo/Bbe1iZ7hFhEREZGITJsGjz8O220H06dDh/wuWfMufdalP2tqXNG99dbZl/jcdlt3vXo1bNoUaa5Ro0ZFur+oKJcf5fITai4IO1shCXUZ6ZCPj1CzKZcf5fJTL9cbb8AFF7jbt9wCffsmEypCeTdLSdZlgj/5BHr1cn8JNXUGe7vt3HKgn3yypYuJiEiOFNosJVraXUS8fPEFDBsGr73mupTcemuicQp2lpJly5ZlfqC57iS1YurHXVZWFun+oqJcfpTLT6i5IOxshSRrm52wkI+PULMplx/l8lOX69e/dsV2v35w/fXJhopQ3hXcGzduzPxASwvumPpxh7p8sXL5US4/oeaCsLMVkqxtdsJCPj5CzaZcfpTLT1VVFTzyCNxwg5sC8O67YZttko4VmbwruLOqnaEk2xzctTRTiYiIiEhQtv/8czj9dHfn8sthcPvqeZd3BXdRUVHmBxI+w92vX79I9xcV5fKjXH5CzQVhZyskWdvshIV8fISaTbn8KJeHTZv47VtvudrssMPyfgrATNrPoMknn4TSUjj0UHjqqew7+PnPYepUmDIFzjsvtpwiIplo0KSISAOTJsHvfgc77eT6b/funXSiOgU7aLKysjLzAy3tUhLTGe6pU6dGur+oKJcf5fITai4IO1shydpmJyzk4yPUbMrlR7la6P/+D/7wB6wx8Le/BVVsRynvCu7q6urMDyQ8S0l5eXmk+4uKcvlRLj+h5oKwsxWSrG12wkI+PkLNplx+lKsFPvkERo+GzZu5r18/t6JkO5V3BXdWCffhFhEREZEWshbOOANWrIBhw5i5zz5JJ4pV+ym4NUuJiIiISH6YNg0efhh69oRZs9iU50u3NyfvBk0OGjTILly4sPEDl10Gl1wCF13kppPJZu5cOPBAN93MvHmR5aqqqgpyCWPl8qNcfkLNBeFmK7RBk1nb7ISFenxAuNmUy49yNeG111wtVlMD994LJ5wQRq4MCnbQZE1NTeYHEu7DXVFREen+oqJcfpTLT6i5IOxshSRrm52wkI+PULMplx/lymLdOjjpJFdsn3kmnHBCGLlilncF98qVKzM/kPAsJZMnT450f1FRLj/K5SfUXBB2tkKStc1OWMjHR6jZlMuPcmVx7rnw1lswYABcd13d5sRzxSzvCu6sWnqGu2dP6NABPv0UAl1yWERERKTdufNOuO026NzZLd3eXM3WjhRewd2hA2y3nbu9enW8mUREREQE3ngDzj7b3b7hBthvv2Tz5FjeFdy9s02I3tIuJbClH3eE3UrGjh0b2b6ipFx+lMtPqLkg7GyFJGubnbCQj49QsymXH+VKs26d66u9fj2ccoqbDjCEXDmUd7OUZF0mePhwKC+Hxx+Hww5reicHHQQvvQT/+hccfHA8QUVEMii0WUq0tLtIgbMWTjsNZsyA/v3dDHEtOTkaiIKdpSTrKNbaM9wt6Q8UwxnuESNGRLavKCmXH+XyE2ouCDtbIQl15oGQj49QsymXH+VKuf12V2x37Qr33Ze12A7184pK3hXcWX3+ubtuScFdO1OJFr8RERERiceiRVDbVWTaNDczSYFqPwX3F1+46y5dmn9uDGe4RURERCRl7VrXb/uLL+D006GsLOlEicq7grtbtn4/tQV3587N7ySGxW+GDBkS2b6ipFx+lMtPqLkg7GyFJGubnbCQj49QsymXn4LOZa2bkeTtt91sJFOnhpErQe1n0OQuu8D778Py5e52U268EcaNg3POcV9xiIjkiAZNiki7d/PNruDu1g3mz4evfCXpRK1WsIMms65alvAZ7kmTJkW2rygplx/l8hNqLgg7WyEJdaXJkI+PULMpl5+CzfXqq3Deee72zTe3uNgO9fOKSt4V3OtqZyNpyKfgjmF593nz5kW2rygplx/l8hNqLgg7WyHJ2mYnLOTjI9RsyuWnIHNVV8OJJ0JNDZx5Jvz4x2HkCkDeFdxZJXyGW0RERKRgWesWtKmogP33hylTkk4UlPZRcG/cCJs2QceO0KlT88+P4Qy3iIiISMG69lr4+9+he3c333ZLZo0rIO1j0ORnn7lfcLdu7nZzPv0Utt0WttnGTVsjIpIjGjQpIu3OP/4B3/++O/n5wANwzDFJJ4pMwQ6arK6ubryxtjtJUVHLdtKjhzsb/tlnsGFDJLnmzJkTyX6iplx+lMtPqLkg7GyFJGObHYCQj49QsymXn4LJ9cEH8KMfuWL7wgtbXWyH+nlFJe8K7srKysYbffpvAxgT+eI3N954YyT7iZpy+VEuP6HmgrCzFZKMbXYAQj4+Qs2mXH4KIteGDW5xm1Wr4Lvfhf/93zByBSjvCu6MfAtuUD9uERERkbb41a/ghRegb1+4++6WjaMrUIVbcGumEhEREZHWmTUL/vQn2GorN1iyd++kEwUt7wruPn36NN4YwBnuCRMmRLKfqCmXH+XyE2ouCDtbIcnYZgcg5OMj1GzK5add51qyBH76U3d7yhQYNqzNuwz184pK3hXcRZkGRgZwhrukpCSS/URNufwol59Qc0HY2QpJxjY7ACEfH6FmUy4/7TZXdTUceyysXw8nnwznnBNGrsDlXcG9bNmyxhsDOMNdVlYWyX6iplx+lMtPqLkg7Gy5Zoy5zRhTaYxZnLbt98aY940xC1OXw9Meu9gYU2GMedsYU5q2/RvGmEWpx/5kjDHNvXbGNjsAIR8foWZTLj/tMpe1cNpp8O9/w9e+5pZub74ZiD9XHsi7gjujAM5wi4gE7A5geIbt11lrB6UujwMYY74KnAQMSP3MNGNMx9Tz/wycCeydumTap4i0V1dfDQ8+CD17wv33Q9euSSfKG4VbcGuWEhEpENba54GWNnYjgbuttTXW2neBCmCoMaYP0MNa+6J1K6bNAI6OJbCIhOeZZ+Dii93tGTOgnXcBiVreFdw9evRovDGAM9ylpaXNPykByuVHufyEmgvCzhaQccaY11NdTrZLbesLLE97zorUtr6p2w23Nyljmx2AkI+PULMpl592levdd+HEE2HzZrjkEjjqqDBy5ZHYJkw0xtwGHAlUWmv3y/D4IcDDwLupTQ9Yayc1t9/emaadCeAM97hx4yLZT9SUy49y+Qk1F4SdLRB/BiYDNnV9DfATIFOHTNvE9kaMMWfiup7Qq1cvRowYUffYddddB8D5559ft23UqFGMHj2asrIyqlJtcr9+/ZgyZQpTp06lvLy87rnTp0+noqKCyZMn120bO3Ysw4cPr/c6Q4YMYeLEiUyaNIl58+bVbZ89ezZz5syhvLy8br8TJkygpKSkXh/S0tJSxo0bx/jx41m6dCkAxcXFTJ8+nZkzZzJr1qxY31NJSYn3e0pfOCSu9wTk9PeUi/cU17FXVVVVb3tI76n2fbXkPRVt3MjUV19lp6oq/r3PPvxq4UI2px6L+j2Vl5cH00ZEvhCPtTaWC/Bt4ABgcZbHDwEe9d1v7969bSM33GAtWDtmTOPHslmwwP3M177W8p9pwnnnnRfJfqKmXH6Uy0+ouawNNxsw38bU7jZ1AfZooj2uewy4GLg47bFy4CCgD/BW2vZRwM3NvW7GNjsAoR4f1oabTbn8tItcmzdbe+KJrl7ae29rV6+OK1awn1dUbXZsXUqsX5/BFqupqWm8MYAz3LV/sYVGufwol59Qc0HY2UKQ6pNd6xigdgaTR4CTjDFFxpg9cYMj51prVwJrjTHDUrOTnIr7lrJJGdvsAIR8fISaTbn8tItcV18N994L3bvDww/DttuGkSsPJb0G50HGmNeAD4ALrLVLWrWXAPpwi4iEyhgzC/etYi9jzArgd8AhxphBuG4hy4CzAKy1S4wx9wJvABuBsdbaTaldnYOb8aQL8ETqIiLt0Zw5cNFF7vadd0L//snmyXNJFtwLgN2ttZ+l5n99CHcmpZH0/oBFRUWN+gMWf/ABxcDf/v537nn99Zb19bGWBzp0YKvPP2fe888z6eqr6/bZmr4+CxcurHt+SP3ManOF1ndu4cKFzJw5M9E+WZneU2VlJUBw/QGLi4uD7ONYWVlZL2tI/TZrj/122x/Qg7V2VIbNtzbx/EuBSzNsnw80GpPTlE6dkj6vk1lx7UmXAIWaTbn85HWuigoYNcrNu/3738PIkWHkymPGdU+JaefG7IHrp91sA22MWQYMttZ+3NTzBg8ebOfPn19/4wUXwDXXwFVXwa9+1fKAffrAhx/CihXQt9nB9iIibWaMecVaOzjpHLmSsc0WkXCtXQsHHeSWbx85Eh54ADrk3aR2kYmqzU7sEzTG7FS7SpkxZmgqS7P9O6oy9bmu7SPo06UEIu3HPXPmzDbvIw7K5Ue5/ISaC8LOVkgyttkBCPn4CDWbcvnJy1ybN0NZmSu2+/d3823nqNgO9fOKSmyfYqrP4IvAvsaYFcaYM4wxZxtjzk495XhgcaoP95+Ak2wLTrdnbLxb04cbIu3Hnf61fEiUy49y+Qk1F4SdrZCEWnCHfHyEmk25/ORlrssu27KS5EMPQQ7n0Q/184pKbJ3rsvQZTH98KjA1khdrbcHdq5e7/uijSGKIiIiI5KXZs2HiRDAGZs6EffZJOlG7EuZoFl+tLbh32sldf/hhtHlERERE8sXixTB6tBskedllcPjhSSdqd/KuF/yuu+7aeGNrC+4+qWloV65sWyi2zO4QGuXyo1x+Qs0FYWcrJBnb7ACEfHyEmk25/ORNro8+ghEj4LPP4Ec/2jIVYNK52pm8K7gzamvB/cEH0eYRERERCd2GDXD88bBsGQweDLff7rqUSOTyruBevnx5442tLbh33tldR3CGO32O35Aolx/l8hNqLgg7WyHJ2GYHIOTjI9RsyuUn+FzWwpgx8Pzzrh56+GHo0iX5XO1U3hXcGQXQpUREREQkb0yZArfe6orshx/echJSYqGCG1Rwi4iISOF44gm3aCDAHXe47iQSq7wruDMu/dnagnuHHdyE7h9/7PoxtcGoUU3OgpgY5fKjXH5CzQVhZyskoS7XHPLxEWo25fITaq6zvvUtOOkkt8jN734HJ56YdCQg3M8rKrEu7R6HjMsE77qrW579v/91t33svLM7w/3ee7DbbtEFFRHJQEu7i0hiPv4YDjwQ3nkHTjgB7r67oJdtb4m8X9q9tZYtW9Z4Y+0Z7qIi/x1G1K2krKysTT8fF+Xyo1x+Qs0FYWcrJBnb7ACEfHyEmk25/ASXq3ZGknfegW98w3UlCajYDu7ziljeLXyzcePGxhtb26UE3BnuBQv8C+6KCrjzTujUCcrKgl2+WLn8KJefUHNB2NkKScY2OwAhHx+hZlMuP0HlshbGjYN//INPiorY/uGHoWvXpFPVE9TnFYO8K7gzakvB3Zoz3E8/DUcdBZ9/7u5fdRVf228//9cWERERidt118Ett0Dnzlw6eDDX9u2bdKKCE853CS1U1LDbyMaN7mIMbLWV/w59F79Ztcr1e/r8czj6aDjiCPjsMy559VVYutT/9WPWr1+/pCNkpFx+lMtfyNkKSaM2OxAhHx+hZlMuP8HkevDBLTOSTJ/O5m98I9k8WQTzecUk/wdNrlsH22zj5pFcv95/h7feCj/9KZxyCsyY0fzzTzsNpk+H0lJ4/HG37dhj3RyW3/0uPPOMVmkSkaw0aFJEcmbuXDjkEHeS8PLLE1u2PZ8V7KDJysrK+hva0p0EYM893fW77zb/3KVLXb/trbaCG290gw06dIDbbuOLrl3h2WfhgQdalyMmU6dOTTpCRsrlR7n8hZyttYwxxxpj/mOM+dQYU22MWWuMqU46V1MatdmBCPn4CDWbcvlJPNeyZTBihCu2zzgDLrwwjFxZhJorKnlXcFdXN/i/JZcF93XXuXkrTz4Z0r/6KC7mtr32crcvvNB1cQlEeXl50hEyUi4/yuUv5GxtcBVwlLW2p7W2h7W2u7W2R9KhmtKozQ5EyMdHqNmUy0+iudascV1eKyvh0EPhz3+u+/Zdn1cy8q7gbqStBfeuu0LHjq4Pd01N069z553u9vnnN3q4fLfdXBG+dCncc0/rsoiING2VtfbNpEOISMBqp/974w0YMAD+/vfWjXGTSKng7tTJFd3WuoVzsnniCaiuhq9/HQYObPTw5g4d4OKL3Z3LLnNnwkVEojXfGHOPMWZUqnvJscaYY5MOJSKBsBbOOceNJ9txR3jsMejZM+lUgsegSWNMN2vtupjzNGvQoEF24cKFWza88goMHuwK4QULWrfT733P9b8uL4cf/jDzc048Ee67D666Cn71q0YPV1VVUbzNNlBSAsuXw/33u8GUCauqqgpyaWXl8qNc/kLN1pYBOMaY2zNsttban7QxVmwatdmBCPX4gHCzKZefRHJddhlccombSOIf/4AhQ8LI1QKh5srZoEljzMHGmDeAN1P39zfGTGvrC7dWTcNuH209ww2wxx7u+p13Mj9eXQ2zZ7vbJ52U8SkVFRWw9dZbivErr3R/aSasoqIi6QgZKZcf5fIXcrbWstaenuESbLENGdrsQIR8fISaTbn85DzX3Xe7YtsYmDkzY7GdSK4WCjVXVFrSpeQ6oBT4BMBa+xrw7ThDNWVlwwVqoii4993XXb+ZpWvkww+71/nWt1z3kwwmT57sbpxxBvTq5abiee651meKSF2uwCiXH+XyF3K21jLG7GKMedAYU2mMWWWMud8Ys0vSuZrSqM0ORMjHR6jZlMtPTnP9619u2mKAa65x64Rkoc8rGS3qw22tXd5g06YYsrROFAV37SqRixdnfnzWLHc9enTz++raFc49192+/PLWZxIRaex24BFgZ6AvMDu1TUQK1VtvudWva2pg7FgYPz7pRJJBSwru5caYgwFrjNnaGHMBqe4lQYii4K4dBLloUePHPvoInnzSDa48/viW7W/sWOjWDZ56yvUxFxGJxg7W2tuttRtTlzuAHZIOJSIJWbkShg+Hqio48kiYMkWL7wWqJQX32cBY3NmUFcCg1P1E9O7du/6GKAruXXeF7t1dcd1wkYa//x02bXKDKXv1yrqLsWPTPpLiYjjrLHf7yitbnysC9XIFRLn8KJe/kLO1wcfGmJONMR1Tl5NJdfcLVaM2OxAhHx+hZlMuP7Hnqq6Gww6D996DAw90UxJ36pR8rlYKNVdU8n9p99tuc/2mTz/d3W6tgw+GF1+Ep5+G739/y/ZvfhNeeMHNwX3yyS3f3/vvu0V1Nm6Et9+GvfdufTYRaTfaOEvJbsBU4CDAAi8A51lr34swYqS0tLtIDDZsgMMPd9P/7b23q1OaOCkorRf7LCXGmBuMMX/KdmnrC7dWo1GsUZzhBjetIMDLL2/ZtnSpO4i7dYNjjmnyx0eMGFF/Q9++cOqpbqaSyy5rW7Y2aJQrEMrlR7n8hZyttay1/7XWHmWt3cFa29tae3TIxTaEO/NAyMdHqNmUy09suTZvdicZa+faLi/3KrYL7vMKRFNdSuYDrzRxCUNUBfe3vuWun39+y7a77nLXxxzjim5fF17ovt6ZPr31c4SLSMEzxvw6dZ3xREjS+UQkhy66yE37t8028Pjj7tt0CV7Wzj7W2um5DNJqURXc307NdPjPf8Lnn7tC+a9/ddtOOaV1+9x7bzdjybXXwnnnuWJegxlExF/tQHX1zRApZNdfD1df7WqUBx6AAw5IOpG0ULO9640xs3F9BdN9imv4b7bWfhFHsGy6NTzTHFXBvfPObsXK+fPdIjdffOFWjOzfHw49tNkfH5JlgnkmTHD9v//5T7jpJrfkag5lzZUw5fKjXP5CzubLWptaeYv11tr70h8zxpyQQKQWa9RmByLk4yPUbMrlJ/Jc994L55/vbt92G/zgB63aTcF8XoFpdtCkMeZ63LRTqcmo+RHwIdAF6GGtbeXp39ZpNADn1792f+1deaW73RbXX+/mr9x9d1dwr1rlDurTT2/bfu+9F370I/dHwQsvbOkvLiIFp42DJhdYaw9obltINGhSJAL/+IebLW3DBrjiCtdlVXIiZ0u7A1+31o621s5OXU4GhlprxwI5b+RjWWmy1s9+5pZ5f+89V2z/z/+4gY8tMGnSpOwPnniim0nliy/cP5hM833HpMlcCVIuP8rlL+RsvowxhxljbgD6Nui/fQewMeF4TQp1pcmQj49QsymXn8hyvfYajBzpiu1x49p8crHdf16BaknBvUNqKiqgblqq2uGwG2JJ1YR169bV3xBlwd21Kzz2mBskefbZ8Mgj0LFji3503rx5TT/hxhvd5PQff+zmy7z2Wvjss7ZnbmuuhCiXH+XyF3K2VvgA143vC+oPXn8EKE0wV7MatdmBCPn4CDWbcvmJJFdFBZSWwqefusX3IljYpl1/XgFrfoZ0+CXwT2PMUsAAewJjjDHdgOQHVkZZcAN89atuIELUiorcfs8+G2bMgF/+En7zG9h/f9hhB5e/Y0fo0MFdd+wIW20FX/kKDB0KBx3k7otIwbHWvga8Zoy5y1ob9BltEYnI+++7ftqrVrmxZH/7W4tPAkp4mi24rbWPG2P2Br6CK7jfShsoOSXGbC0TdcEdpy5d3BSBxx/v+pz/618wd27Lfnb77eHoo13XlGHDNNuJSAExxtxrrT0ReNUYkz7wxgDWWvu1hKKJSByqqtyZ7WXL3Em3Bx90J+4kb7VopUljzMHAHqQV6NbaGfHFyq7RAJyjjnKzijz8sLudT6qq4I033HVNjVtCfvNmd71pk5ue8LXX4Lnn3GqVtQ44AH7+czcQs0uXxOKLiL/WDMAxxvSx1q40xuye6fGQF7/RoEkRT+vWuTPaL73kvnV//nl30k0SkbNBk8aYO4E/Av8DDEld2vzCrVVdXV1/Q+0Z7oT/8pszZ47/DxUXu4GZRx0FJ5wAJ50Eo0e7eb9PO81NIXjTTfDmm7B4sRuVvP32bhGd00+HXXaBX/3KrYgZZa4cUC4/yuUv5Gy+rLW1Iw8/BpanCuwiYH9c/+5gNWqzAxHy8RFqNuXy06pcNTVw7LGu2N59d3jyyciL7Xb1eeWRlgyaHAx801o7xlr789Tl3LiDZVNZWVl/QyBdSm688cb4dm4MDBjgpgJavhxuv92d5a6qgj/+EUpK3IDMO+9023KVqw2Uy49y+Qs5Wxs8D3Q2xvQFngFOB+5INFEzGrXZgQj5+Ag1m3L58c61aZM74fbkk25s11NPQd++yefKkVBzRaUlBfdiYKe4g7RaIAV3znTp4s5+z58PL7/sbnfuDOXlbgrD3r3dMvUXXQQPPcTOn30GX36ZdGoRiYax1q4HjgVusNYeA3w14Uwi0lbWwpgxcN990KOH+z99772TTiURasksJb2AN4wxc4Ga2o3W2jA6TBdawV3LGDeQYuhQuOYamDXLDap47jm3quU//wnAzeCmO9xzT9htN9h1V9cVJf16112hZ08NxBQJnzHGHAT8GDgjta0l7biIhOySS+Avf3G1zOzZWiCvHWrJSpPfybTdWvuPWBI1o3///vbNN9/csmGffeA//4G33oJ9900iEgBz585l6NChib1+ndWr3ewnL74I8+ZRs3gxRR9+6P56bkq3bo2L8IaFec+ekcUM5vNqQLn8hJoLws3WxpUmv4ObqvVf1torjTF7AeOT7ObXnEZtdiBCPT4g3GzK5afFua680n0r3bGjO3E2YkQYuXIs1FxRDZpsybSA9QprY8w3gdFAIgV3UcPBkYGc4S4pKUn09etstx0ceaS7AOuqqijq0gXeecf1/16+HFasqH+9fLkbFf3WW+6STbdursvKDju4S+3t4mL3FVj37luuG97u2tXNMZ4SzOfVgHL5CTUXhJ2ttVLt8T+MMd2NMdtYa98Bgi22IUObHYiQj49QsymXnxbluuEGV2wbA3fcEXux3eJcCQg1V1RaOi3gIFyRfSLwLnC/tXZqvNEy23bbbe2aNWu2bOjdGz76CD78EHbcMYlIAIwYMYLZs2cn9vrZtCiXtW4Vq4ZFeMPCfP36toUpKnKFd5cufLBmDTv361d3P+t1UZFb8Kc1l06dtiwklH6daVvq+uRTT+Vvs2Y1/XPG5Lz7TV4fXwkJNVsbz3APBGYAxbg5uD8CTrXWLokwYqQatdmBCPX4gHCzKZefZnP99a/ws5+52zffDGeeGUauhISaK/Yz3MaYfYCTgFHAJ8A9uAL9u2190UgFcoY7rxkD227rLgMHZn6OtVBd7f64qax017W3V6+GtWvdpbo68+316910RzU1sHo1OwMsWpS799hCfwPo1avlP1BbeNcW4THd/9vate6Py5b+vI/W/vFgDLd+9JGbuirO12rlz9yyahXstVfsr5NjNwO/sNY+617eHALcAhzc1A8ZY24DjgQqrbX7pbYV49r1PYBlwInW2tWpxy7G9RHfBJxrrS1Pbf8GblaULsDjwHm2JWdtRKS+u+7aUmBPmZKzYluS01SXkreA/wNGWGsrAIwx5+cklQ8V3LlhjOvD3bOnm4bQl7Xud/X557B+PWeecgp/mTKl7n7ddcPbGza4WVZac6ldQKh2MaH060zbNm3i0zVr6LnNNk3/XHp9UXs75pqjJ7g/cALTG+C//006RkY7Abz7btIxotatttgGsNY+Z4zp1oKfuwOYijs7Xusi4Blr7RXGmItS9y80xnwVd7JlALAz8LQxZh9r7Sbgz8CZwEu4gns48ETb35ZIAXngASgrc/9vXHYZnHde0okkB5oquI/DNbrPGmPmAHfjvsJMVI8ePbbc2bRpy5R3W2+dTKCU0tLSRF8/m2ByGeO6iHTpAsXFfO2442D//ZNO1cidU6cybty45p9obf1iO+b7t/71r5zxk5+07Pk+WvuHQurnpk+fTllZWXyv1YafmTFjBqeeemrsr+P9M/vs4/9zW7xjjJkA3Jm6fzKum18zL2ufN8bs0WDzSOCQ1O3pwHPAhantd1tra4B3jTEVwFBjzDKgh7X2RQBjzAzgaJopuOu12QEJpm3MINRsyuUnY67HH3eL3G3aBL/9LVx8cRi5AhBqrqi0ZJaSbrhGdRTwPVzD/KC19slmfq7RV5gNHjfA9cDhwHrgNGvtguYC11smeP16N5Cvc2d3VlREJHBt7MO9HfAH3Mq/4BbC+UNtV5BmfnYP4NG0LiVrrLXbpj2+2lq7nTFmKvCStfZvqe234orqZcAV1tpDU9u/BVxorT2yqdfV0u4iKc88A0cc4bpW/uIXbuE6TccbvFzOUrIOuAu4K9Xn7wTcV49NFtxk/goz3WHA3qnLgbivKg9sLs/y5cu33AmoO8n48eOZMmVK0jEaUS4/yuUn1FwQdjZfxpjOwNlACbAI+KW1Nq4VrTJVALaJ7Y13YMyZuK4ndO7cmRFpMy9cd911AJx//pYeiqNGjWL06NGUlZVRlVott1+/fkyZMoWpU6dSXl5e99zp06dTUVHB5MmT67aNHTuW4cOH13udIUOGMHHiRCZNmsS8efPqts+ePZs5c+YwZswYBgwYAMCECRMoKSmp921NaWkp48aNY/z48SxduhSA4uJipk+fzsyZM5k1a1Zs72nJkiVMmzbN+z2lr9QXx3v68ssvmTNnTk5/Ty15TxUVFQA5/z01954OOuggeqXGBPWvquKqhQuhpobHd9+dP7/9Nhx1VCLH3rXXXlt37Ofy99Tce1qyZAkDBgwIpo2IfOVLa21sF9xgnMVZHrsZGJV2/22gT3P77Nmzp63z/vvui/SddrJJO/LII5OOkJFy+VEuP6HmsjbcbMB869+W3oMb03sW8BAwpRX7qNcep7e5QB/g7dTti4GL055XDhyUes5badtHATc397r12uyAhHp8WBtuNuXyU5dr7lxre/Rw9UpZmbWbNoWRKzCh5mpNm53p0pKl3ePSF0g7Xc2K1LaWC+gMt4hIjL5qrT3ZWnszcDzw7Qj2+QhQewqqDHg4bftJxpgiY8yeuG8h51prVwJrjTHDUl0CT037GRHJZP58+MEP3KxdJ57opgLskGTpJUlJckngVn09WVRUVPeVwK5r1zIN2NChA8elfU2QxFcPCxcurHt+SF9P1uZK6iuibO9p4cKFzJw5M9GviDK9p8rKSoBEvkZu6j0VFxcn9jVyU++psrKyXtaQvp6sPfbbydeTdd1HrLUbjWe/T2PMLNwAyV7GmBXA74ArgHuNMWcA/8V1F8Rau8QYcy/wBrARGGvdDCUA57BlWsAnaMEMJZ06hbnyfHFxcdIRsgo1m3L5OWDzZldsf/opHHss/O1vbn2IhIX6eYWaKyotWvim1TtvMEinwWM3A89Za2el7r8NHJI6i5JVvQE4r74KBxwAgwa52yIigWvNABxjzCZgXe1dXMG7PnXbWmvDnAoEDZqUAvXKK3DoobBmjSu2777bLcgmeSeqQZNJfq/xCHCqcYYBnzZXbAN1Z6OAoLqUzJw5M+kIGSmXH+XyE2ouCDubL2ttR2ttj9Slu7W2U9rtYIttaNBmByTk4yPUbMrVQunF9jHHBFdsB/d5pYSaKyqxFdyprzBfBPY1xqwwxpxhjDnbGHN26imPA+8AFbiV0sa0ZL+hFtzpX8uHRLn8KJefUHNB2NkKSagFd8jHR6jZlKsFFixw3UjWrOHFHXcMrtiGwD6vNKHmikpsnYmstaOaedwCY9v0IgEV3CIiIlLAFixwZ7ZXr4ajj+aqDRt4MOFF+SQc+T1UVgW3iIiIJC292B45Eu65h42ajUTSxDpoMg4DBw60ixYtcndmzYLRo90yqQl/FVFRUUFJSUmiGTJRLj/K5SfUXBButqgG4OSLem12QEI9PiDcbMqVRcNi+957Yeutk8+VhXL5aQ+DJttOZ7hFREQkKS+9BN/7niu2jzqqrtgWaSjvCu5Ql3ZPn+M3JMrlR7n8hJoLws5WSOq12QEJ+fgINZtyNfD881vm2T7uOLjvvnrFtj4vP6HmikreFdz1BFRwi4iISIF46ikYPhw++8x1bb37bp3Zliap4BYRERFpqUcfhREj4PPP4YwzYMaMIFaQlLDlXcFdb+nPgAruUaOanAUxMcrlR7n8hJoLws5WSEJdrjnk4yPUbMoF3H+/WzmypgbGjIG//AU6dkw+lwflSkbezVJSb5ngCy+Eq66CK65wt0VEAldos5RoaXdpN2bOhFNPhU2b4Je/hKuvBmOSTiUxK9hZSpYtW7blTu0Z7qKiRLKkKysrSzpCRsrlR7n8hJoLws5WSOq12QEJ+fgINVtB57rtNjj5ZFdsT5jQomK7oD+vVgg1V1TyruDeuHHjljsBdSkJdfli5fKjXH5CzQVhZysk9drsgIR8fISarWBz/elPrq+2tXDppTBpUovObBfs59VKoeaKSt4V3PUEVHCLiIhIO2It/O53cN557v6118JvfpNsJslbeTestii9+0hABXe/fv2SjpCRcvlRLj+h5oKwsxWSogC6/GUS8vERaraCyrV5syu0p06FDh3gr3+F009PPlcElCsZ+T1ocuRIeOQRePBBOProRHOJiLSEBk2KBO7LL+G009wgya23dnNsH3NM0qkkIQU7aLKysnLLnYDOcE+dOjXpCBkplx/l8hNqLgg7WyGp12YHJOTjI9RsBZFr/XpXXM+cCdtsA0880epiuyA+rwiFmisqeVdwV1dXb7kTUMFdXl6edISMlMuPcvkJNReEna2Q1GuzAxLy8RFqtnafa80aKC2Fxx6D7beH//f/4HvfSz5XxJQrGXnXh7uegApuERERyVOrVrml2hcuhL593dLt/fsnnUrakbw7w12PCm4RERFpi3fegW99yxXb++wD//qXim2JXN4Nmhw0aJBduHChu7PvvvDvf8MbbyT+j6OqqirIJYyVy49y+Qk1F4SbrdAGTdZrswMS6vEB4WZrl7leeQUOPxwqK+HrX4c5c6B37+RzxUi5/BTsoMmampotd2rPcHfpkkyYNBUVFUlHyEi5/CiXn1BzQdjZCkm9NjsgIR8foWZrd7nmzIHvfMcV24ceCs89F1mx3aZcMVOuZORdwb1y5cotdwIquCdPnpx0hIyUy49y+Qk1F4SdrZDUa7MDEvLxEWq2dpXrjjtgxAhYt84t2f7YY9CjR/K5ckC5kpF3BXc9n3/urtWHW0RERJpTuzz76afDxo1w0UUwY4abb1skRpqlRERERNq/jRth3Di4+WYwBm64AcaOTTqVFIi8K7h71/av2rTJrQZlTBB/mY4N9B+tcvlRLj+h5oKwsxWS3hH2iY1SyMdHqNnyOtf69TBqlFudunNnt7BNzKtH5vXnlYBQc0Ul72YpqVsmeP166NbN9d9evz7pWCIiLVJos5RoaXdJXGUljBwJL70E220Hs2fDN7+ZdCrJEwU7S0ndKNbA+m+PGDEi6QgZKZcf5fITai4IO1shCXXmgZCPj1Cz5WWuN96AYcNcsb3bbm6O7RwV23n5eSUo1FxRybuCu476b4uIiEg2Tz8NBx8M774LQ4bAyy8nvmaHFC4V3CIiItK+3HKLW6r900/huOPcHNs77ZR0KilgeVdwd+vWzd0IaA5ugCFDhiQdISPl8qNcfkLNBWFnKyR1bXZgQj4+Qs2WF7k2bYJf/QrOPNPdvugiuPde6No12VwBUa5k5O+gyfnz3VdEBxzglmYVEckDGjQpEpN16+DHP4aHH4ZOndz0fz/5SdKpJM8V7KDJulXLAutSMmnSpKQjZKRcfpTLT6i5IOxshSTUlSZDPj5CzRZ0rg8+gG9/2xXb224LTz6ZeLEd9OcVoFBzRSXvCu5169a5G4EV3PPmzUs6QkbK5Ue5/ISaC8LOVkjq2uzAhHx8hJot1FxrnnzSfeO9YAH06+dmJPnud5OOFeznpVzJyLuFb+oE1odbREREcuzOO7nixRdh82Z3hvv++6FXr6RTiTSSd2e46wQ2D7eIiIjkyMaNcMEFcOqpbL15M5x9Njz1lIptCVb+Dpq880449VQ3QOJvf0s6lohIi2jQpEgbrV4NJ53k+ml36gQ33OAKbpEYFOygyerqancjsD7cc+bMSTpCRsrlR7n8hJoLws5WSOra7MCEfHyEmi2IXG++CUOHumJ7hx3gmWeYs8ceSafKKIjPKwPlSkbeFdyVlZXuRmB9uG+88cakI2SkXH6Uy0+ouSDsbIWkrs0OTMjHR6jZEs/16KNw4IFQUQGDBsG8efDtbyefKwvl8hNqrqjkXcFdR324RURE2r/Nm2HSJDjqKFi7Fk48Ef75T9h996STibRY/s9SooJbRESkfVq9Gk4+GR5/HIyBSy+Fiy92t0XySN4V3H369HE3Aiu4J0yYkHSEjJTLj3L5CTUXhJ2tkNS12YEJ+fgINVvOc736Khx3HLz7LhQXw8yZUFqafK4WUi4/oeaKSt51KSkqKnI3AuvDXVJSknSEjJTLj3L5CTUXhJ2tkNS12YEJ+fgINVtOc91+Oxx8sCu2Bw92i9pkKLZznsuDcvkJNVdU8q7gXrZsmbsRWB/usrKypCNkpFx+lMtPqLkg7GyFpK7NDkzIx0eo2XKS64sv4Kyz3LLsX3wBZ54J//d/TfbXLujPqxWUKxl516WkTmBdSkRERKQN3nsPjj8e5s93/7dPmwann550KpFI5N0Z7joquEVE2swYs8wYs8gYs9AYMz+1rdgY85Qx5j+p6+3Snn+xMabCGPO2MSbzd/wivmbPhgMOcMX2nnvCCy+o2JZ2JdaC2xgzPNUoVxhjLsrw+CHGmE9TDf1CY8zE5vbZo0cPdyOwPtylWfqWJU25/CiXn1BzQdjZAvRda+2gtNXULgKesdbuDTyTuo8x5qvAScAAYDgwzRjTsakd17XZgQn5+Ag1Wyy5NmyAX/zCTflXVQVHHAGvvAJf/3qyuSKgXH5CzRWV2JZ2TzXC/wZ+AKwA5gGjrLVvpD3nEOACa+2RLd1v3TLBw4dDebmbKuiww6INLyISk9CWdjfGLAMGW2s/Ttv2NnCItXalMaYP8Jy1dl9jzMUA1trLU88rB35vrX0x2/61tLtk9c478KMfubPanTrBFVfA+edDh/z98l3an6ja7Dj7cA8FKqy17wAYY+4GRgJvNPlTzVi+fLm7EViXkvHjxzNlypSkYzSiXH6Uy0+ouSDsbIGxwJPGGAvcbK39C7CjtXYlQKro7p16bl/gpbSfXZHaVo8x5kzgTIDOnTszYsSIuseuu+46AM4///y6baNGjWL06NGUlZVRVVUFQL9+/ZgyZQpTp06lvLy87rnTp0+noqKCyZMn120bO3Ysw4cPr/c6Q4YMYeLEiUyaNIl58+bVbZ89ezZz5sxhzJgxDBgwAHDTkZWUlNQbtFVaWsq4ceMYP348S5cuBaC4uJjp06czc+ZMZs2aFdt7WrJkCdOmTfN+T+kr9cXxnr788kvmzJkTye/pmx98wM9ff51uGzfC7rtz6w9/yEPPPQfPPef9nioqKgBy/ntq+J4a/p4OOuggevXqlfPfU3Pv6dprr6079qP69xTFe1qyZAkDBgwIpo2IfOVLa20sF+B44K9p908BpjZ4ziHAJ8BrwBPAgOb227NnT2uttfbAA60Fa194wYbgyCOPTDpCRsrlR7n8hJrL2nCzAfNtTO1uay7Azqnr3qm2+NvAmgbPWZ26vhE4OW37rcBxTe2/rs0OTKjHh7XhZosk1+efW3vOOe7/b7D22GOtrapKPlcMlMtPqLmiarPjPMOdaRmohv1XFgC7W2s/M8YcDjwE7N1oR2lnS7baaitGjBjB9W++yV7Afz/6iA0VFYn/JTR37ty654d0tqQ2V1J/sWZ7T3PnzmXmzJmJ/sWa6T29++67AImc1WrqPQGJndVq6j29++679bKGdLak9thvt2dLImKt/SB1XWmMeRD37eQqY0wfu6VLSWXq6SuAXdN+fBfgg5wGlvz11luuC8nrr8PWW8O118KYMVo1UgpDFFV7pgtwEFCedv9i4OJmfmYZ0Kup52y//fbuT4599nF/Hb/5Ztv/fInAqaeemnSEjJTLj3L5CTWXteFmI6Az3EA3oHva7RdwgyGvBi5Kbb8IuCp1ewDuLHgRsCfwDtCxqdeoa7MDE+rxYW242Vqda/Nma6dNs7ZLF/f/dkmJtQsWJJ8rZsrlJ9RcUbXZcQ6a7IQbNPl94H3coMnR1tolac/ZCVhlrbXGmKHA33FnvLOGqhuAs/vu8N//ulWo9tgjlvcgIhK1kAZNGmP2Ah5M3e0EzLTWXmqM2R64F9gN+C9wgrW2KvUzlwA/ATYC4621TzT1Gho0WeBWrYIzzoDHHnP3TzkFbrwRundPNpdIC0XVZsc2FNhauxEYB5QDbwL3WmuXGGPONsacnXra8cBiY8xrwJ+Ak5oqtoG6r39DGzQ5c+bMpCNkpFx+lMtPqLkg7GyhsNa+Y63dP3UZYK29NLX9E2vt9621e6euq9J+5lJrbT9r7b7NFduQ1mYHJuTjI9Rs3rlmz4aBA12xvd12cM89MGNG5MV2u/m8ckS5khHr3DvW2settfukGufahvwma+1NqdtTU438/tbaYdbaF5rbZ6OCO5B5uNP7wYZEufwol59Qc0HY2QpJqAV3yMdHqNlanGvdOjj7bDe39kcfwfe+5/ptn3hisrlyTLn8hJorKvk72eXnn7vrQM5wi4iIFLz5892KkTff7AZG/vGP8NRTsMsuSScTSVScs5TEZ9Mm+PJLd3vrrZPNIiIiUug2bIDLLoNLL4WNG2HAALjrLth//6STiQQhtkGTcRk4cKBd9PLL0K2bO7tde6Y7YRUVFZSUlCQdoxHl8qNcfkLNBeFmC2nQZC4MHDjQLlq0KOkYjYR6fEC42bLmWrgQTjsNXnvN3T/vPLdqZI6+gc67zythyuUn+EGTsVJ3EhERkWR9+SX84Q8wZIgrtvfay60UOWWK/n8WaSDvCu7ly5fD+vXuTrduyYZJk76oRkiUy49y+Qk1F4SdrZAsX7486QgZhXx8hJqtXq7XXoOhQ+H3v3ddSMaNcwMjv/OdZHMFRLn8hJorKnlXcANuBDQEVXCLiIi0e19+CZMmweDBrivJnnvCs8/CDTfo/2SRJuTnoMnaM9xduyabQ0REpEDsvXq1K7Rff91tGDvW9dXeZptkg4nkgbwruIuLi4MsuEeNGpV0hIyUy49y+Qk1F4SdrZAUFxcnHSGjkI+P4LKtXQu//S3XvPACWOvOat96K3z3u0knAwL8vFKUy0+ouaKSd7OUDB482M6/7DIoLYVDD3Xze4qI5IlCm6VES7vnudmzYcwYWLECOnaECy6AiRODOuElEqeCnaVk2bJlQZ7hLisrSzpCRsrlR7n8hJoLws5WSJYtW5Z0hIxCPj6CyLZypVsZ8qijXLE9eDATDjvMdSEJ6P9eCOTzykC5/ISaKyp5V3Bv3LgxyII71OWLlcuPcvkJNReEna2QbNy4MekIGYV8fCSabdMmt0pk//5w331uIOSUKfDSSyxMLlWTQv1dKpefUHNFJe/6cAOapURERCRqL7/sBkK+8oq7f8QRMG0a7LZbsrlE2oG8O8NdVFQU5Bnufv36JR0hI+Xyo1x+Qs0FYWcrJEVFRUlHyCjk4yPn2T76CH76Uxg2zBXbu+wC997r+m+nFduhfmbK5Ue5kpGfgyaPOw5+8xu48ELXn0xEJE9o0KQEo7b7yG9/C6tXw1ZbwS9/CZdcoqn+RFIKdtBkZWVlkGe4p06dmnSEjJTLj3L5CTUXhJ2tkFRWViYdIaOQj4+cZHvxRbck+9ixrtj+4Q9h0SK4/PKsxXaon5ly+VGuZORdwV1dXR1kwV1eXp50hIyUy49y+Qk1F4SdrZBUV1cnHSGjkI+PWLO99x6MHg0HHwyvvgq77gr33w9z5sC++yaXqw2Uy49yJUODJkVERNq76mrXBfPaa6GmBoqKXPeR3/xG/5eK5EB+FtwBnuEWEREJzsaNcNttMGEC1HbvGT0aLrsMdt892WwiBSTvBk0OGjTILiwpcV+B3XsvnHBC0pEAN39kiEsYK5cf5fITai4IN1uhDZocNGiQXbhwYdIxGgn1+IAIsz35pDuLvXixu3/QQe4M97BhyeaKmHL5US4/BTtosqamJsgz3BUVFUlHyEi5/CiXn1BzQdjZCklNTU3SETIK+fhoc7Z589wgyNJSV2zvsQfccw/861+tLrYjyRUT5fKjXMnIu4J75cqVQRbckydPTjpCRsrlR7n8hJoLws5WSFauXJl0hIxCPj5ane3NN+H442HoUHjqKejZE6680m0/8UQwJplcMVMuP8qVjPzsw61BkyIiIs5778Ef/gDTp8PmzdC5M5x7rlurIsCv6EUKUX4W3AGe4RYREcmpDz90M4/8+c+wYQN06gRnnukGSO68c9LpRCRN3hXcvXv3DrLgHjt2bNIRMlIuP8rlJ9RcEHa2QtK7d++kI2QU8vHRbLYPPoCrrnKrRH7xhesqMnq0O8tdUpJcroQolx/lSkbezVIyePBgO3/5cje90cqVsNNOSUcSEWmxQpulREu7R2jFCtcn+5Zb3FzaACNHukJ7//2TzSbSThXsLCUVFRVBnuEeMWJE0hEyUi4/yuUn1FwQdrZCEurMAyEfH42yvfcenHMO9OsHU6e6Yvv442HhQnjooZwV26F+ZsrlR7mSkXddSoAtgyYDKrhFREQitXgx/PGPcNddbgEbY+Ckk+CSS2C//ZJOJyIe8q7g7mAtWOuK7U55F19ERCQ7axn48cdw+OHwxBNuW4cO8OMfu0K7f/9k84lIq+RdxdqtSxeoroYePZKOUs+QIUOSjpCRcvlRLj+h5oKwsxWSboFO3xrc8bFxIzzwAFx9NZfV9nnv0gV+8hP4xS9gr72SzUeAn1mKcvlRrmTk36DJ/faz85csgX32gbffTjqOiIgXDZqUetasgTvugBtugHfecdt69YKf/xzGjHG3RSQxBTto8uNVq9yN7t2TDdLApEmTko6QkXL5US4/oeaCsLMVklBXmkz8+Fi8GM4+G/r2hfPPd8V2v34wbRqXnXUWTJwYXLGd+GeWhXL5Ua5k5F2XkprPP3c3AutSMm/evKQjZKRcfpTLT6i5IOxshWRd7SD3wCRyfGzcCI884s5mP/fclu3f/z6MHQtHHQUdO/JioLM1hPpvSrn8KFcy8q7g7lDbBSawgltERCSj//4Xbr8dbr0Vli9327p1g7IyV2h/9avJ5hOR2OVvwR1YlxIREZE6NTXubPatt8KTT7rZtQD23hvGjXPFds+eyWYUkZzJv0GTu+3mVpocMwZuvDHpOCIiXjRosp1bssQV2XfeCR9/7LZtvTUceyyccQZ873tumj8RyQsFO2jyi9pVJgPrUjJnzpykI2SkXH6Uy0+ouSDsbIWkuro66QgZRXp8LF8OV18NX/+6W5Dmuutcsf21r8H118MHH8CsWXDooS0qtkM9dpXLj3L5CTVXVPKu4P78s8/cjcAK7hsDPduuXH6Uy0+ouSDsbIWksrIy6QgZtfn4qKqCv/wFDjkEdt8dfv1rt9T6ttvCWWfBvHnu/rnnwvbb5zZbTJTLj3L5CTVXVNSHW0REpCU+/tj1y37wQSgvhy+/dNs7d4YRI2D0aDjsMCgqSjaniAQn7wrujrUF9w47JBtERETav+XL4aGH3CqQzz8Pmze77R06wA9/6JZcP/ro4L51FZHAWGvz6jKoSxdrwdpnnrEhefnll5OOkJFy+VEuP6HmsjbcbMB8G0BbmqvLV77ylYg+uWhlPT42bbJ27lxr//AHa4cMcf/f1F622sra0lJrb77Z2g8/zH22hCmXH+XyE2quqNrs/DvDXXt2IbAz3CUlJUlHyEi5/CiXn1BzQdjZCklRoN0r6h0fH33kuojMmeOua2cXAeja1XUTOeYYOOII10c7l9kColx+lMtPqLmikneDJjdt2OBuBFZwl5WVJR0hI+Xyo1x+Qs0FYWcrJMuWLUs6QmNr1jB1+HD45S9h8GDYcUc45RS46y5XbO++u1t2/eGHXTH+97+7riM5KLYh3GNXufwol59Qc0Ul/85w1/bh7tUr2SAiIhI+a91Kj/Pmwb/+Bf/4ByxcyERr3TZwgxy/8x0YPtydzd53XzAm2dwi0q7kXcENQHExdMrP6CIiEhNrYdUqWLDAFdPz5sHcue4sdbqttmJJ9+4MGDPGFdoHH+y6joiIxCTWqtUYMxy4HugI/NVae0WDx03q8cOB9cBp1toFze44sO4kAKWlpUlHyEi5/CiXn1BzQdjZ8l1zbXu6HnHN3rF2rVvVcdEid1m82F2n97+uVVwMQ4e6yyGHwLBhPHvrrQwYNy6ebG0U6rGrXH6Uy0+ouaIS29LuxpiOwL+BHwArgHnAKGvtG2nPORz4Oa7gPhC43lp7YFP7HWyMnV9a6ga3iIjkmXxf2r0lbXu6Ni3t/tlnsGIF/Oc/Wy4VFe76vfcy/0yPHrD//jBkiCuwhwyBPfdUFxERaZV8WNp9KFBhrX3HWrsBuBsY2eA5I4EZqZlXXgK2Ncb0aXbPX/lK5GHbavz48UlHyEi5/CiXn1BzQdjZ8lxL2vY66xcvhpNOghNPhOOPh+OOczN+jBwJRx0FRx4Jhx/u+k+XlrouHvvs4xY3694d+vd3z/vlL+Gmm+Dpp12xvfXWrrA++WS44gp47DHXV3vNGjdf9jXXwI9+BHvtlbHYDvn4CDWbcvlRLj+h5opKnF1K+gLL0+6vwJ3Fbu45fYGV6U8yxpwJnAnwDWDiSy/x6ogRAFx33XUAnH/++XXPHzVqFKNHj6asrIyqqioA+vXrx5QpU5g6dSrl5eV1z50+fToVFRVMnjy5btvYsWMZPnw4I1KvATBkyBAmTpzIpEmTmFc70AaYPXs2c+bMYdasWSxduhSACRMmUFJSUm/EbWlpKePGjWP8+PF1zysuLmb69OnMnDmTWbNm1T03yvc0d+5cli5d2qr3lL7MatTvae7cuQwdOjTnv6fm3tO7774LkPPfU3PvaenSpYn8npp7T08//XS9rLn6PbXkPdUe+6G0Ee1o2eJm2/aGbTb33NO6V+rcmTXbbMNSY1jZtSsfdOvGqAkTeLdTJybefjubOnSANWsYu//+7abNBtc+Dh8+PKg2G+D9998HCKrNLi0tZenSpcG12RMnTuSJJ56oy+T7nuI89tKPfbXZuWuz4+xScgJQaq39aer+KcBQa+3P057zGHC5tfafqfvPAL+21r6Sbb8lRUW24osvgvt6cMSIEcyePTvpGI0olx/l8hNqLgg3WzvoUtJs255u165d7fLbbnMrMxrjrhveTr+/9dZumr4+fdw0fDG19aEeHxBuNuXyo1x+Qs0VVZsd5xnuFcCuafd3AT5oxXPqWdO9e3DFNri/1EKkXH6Uy0+ouSDsbHnOq93+vGtX16UkMCEfH6FmUy4/yuUn1FxRifMMdyfcwJrvA+/jBtaMttYuSXvOEcA4tgya/JO1dmhT+23TABwRkYS1gzPczbbt6dRmi0g+C37QpLV2I66YLgfeBO611i4xxpxtjDk79bTHgXeACuAWYExz+63tuxOamTNnJh0hI+Xyo1x+Qs0FYWfLZ9na9mzPV5vtL9RsyuVHufyEmisqsS7tbq193Fq7j7W2n7X20tS2m6y1N6VuW2vt2NTjA621zZ4GCbXxTh/EEBLl8qNcfkLNBWFny3eZ2vZs1Gb7CzWbcvlRLj+h5opKrAW3iIiIiEihU8EtIiIiIhKj2AZNxmXgwIF20aJFScdopKKigpKSkqRjNKJcfpTLT6i5INxs+T5o0pfabH+hZlMuP8rlJ9RcwQ+aFBERERGRPCy4ly9f3vyTEpC+IlJIlMuPcvkJNReEna2QqM32F2o25fKjXH5CzRWVvCu4RURERETyiQpuEREREZEY5d2gSWPMWuDtpHNk0Av4OOkQGSiXH+XyE2ouCDfbvtba7kmHyBW12a0Sajbl8qNcfkLNFUmb3SmKJDn2dogj/I0x85Wr5ZTLj3L5CzWbMabQ1jlXm+0p1GzK5Ue5/IScK4r9qEuJiIiIiEiMVHCLiIiIiMQoHwvuvyQdIAvl8qNcfpTLX6jZQs0Vl1Dfb6i5INxsyuVHufy061x5N2hSRERERCSf5OMZbhERERGRvBFswW2MGW6MedsYU2GMuSjD48YY86fU468bYw7IQaZdjTHPGmPeNMYsMcacl+E5hxhjPjXGLExdJsadK/W6y4wxi1Kv2WhEbUKf175pn8NCY0y1MWZ8g+fk5PMyxtxmjKk0xixO21ZsjHnKGPOf1PV2WX62yWMxhlxXG2PeSv2eHjTGbJvlZ5v8nceQ6/fGmPfTfleHZ/nZXH9e96RlWmaMWZjlZ+P8vDK2DSEcY7miNts7m9rs5vOo3W57LrXbmfeb+zbbWhvcBegILAX2ArYGXgO+2uA5hwNPAAYYBrycg1x9gANSt7sD/86Q6xDg0QQ+s2VAryYez/nnleF3+iGwexKfF/Bt4ABgcdq2q4CLUrcvAq5szbEYQ64fAp1St6/MlKslv/MYcv0euKAFv+ecfl4NHr8GmJjA55WxbQjhGMvFRW12q7KpzW4+g9rttudSu515vzlvs0M9wz0UqLDWvmOt3QDcDYxs8JyRwAzrvARsa4zpE2coa+1Ka+2C1O21wJtA3zhfM0I5/7wa+D6w1Fr7Xg5fs4619nmgqsHmkcD01O3pwNEZfrQlx2Kkuay1T1prN6buvgTsEtXrtSVXC+X886pljDHAicCsqF6vpZpoGxI/xnJEbXb0CrrNBrXbUeRqoYJrt5Nos0MtuPsCy9Pur6BxI9mS58TGGLMH8HXg5QwPH2SMec0Y84QxZkCOIlngSWPMK8aYMzM8nujnBZxE9n9QSXxeADtaa1eC+8cH9M7wnKQ/t5/gznJl0tzvPA7jUl+Z3pblq7YkP69vAaustf/J8nhOPq8GbUM+HGNRUJvtT2126+TDvym12y2XeLudqzY71ILbZNjWcDqVljwnFsaYbYD7gfHW2uoGDy/AfQW3P3AD8FAuMgHftNYeABwGjDXGfLvB40l+XlsDRwH3ZXg4qc+rpZL83C4BNgJ3ZXlKc7/zqP0Z6AcMAlbivgZsKLHPCxhF02dJYv+8mmkbsv5Yhm35Nn2U2mx/arPjo3Z7C7XbTchlmx1qwb0C2DXt/i7AB614TuSMMVvhfjl3WWsfaPi4tbbaWvtZ6vbjwFbGmF5x57LWfpC6rgQexH3lkS6RzyvlMGCBtXZVwweS+rxSVtV+RZu6rszwnKSOszLgSODHNtVprKEW/M4jZa1dZa3dZK3dDNyS5fWS+rw6AccC92R7TtyfV5a2IdhjLGJqsz2pzW61YP9Nqd32k3S7nes2O9SCex6wtzFmz9Rf2icBjzR4ziPAqcYZBnxa+zVAXFJ9jW4F3rTWXpvlOTulnocxZijuM/4k5lzdjDHda2/jBm8sbvC0nH9eabL+BZvE55XmEaAsdbsMeDjDc1pyLEbKGDMcuBA4ylq7PstzWvI7jzpXev/RY7K8Xs4/r5RDgbestSsyPRj359VE2xDkMRYDtdl+udRmt16Q/6bUbrdKYu12Im22jWGkbBQX3Ajtf+NGgl6S2nY2cHbqtgFuTD2+CBicg0z/g/va4HVgYepyeINc44AluFGrLwEH5yDXXqnXey312kF8XqnX7YprjHumbcv554X7z2Ml8CXur9MzgO2BZ4D/pK6LU8/dGXi8qWMx5lwVuP5htcfYTQ1zZfudx5zrztSx8zqucekTwueV2n5H7TGV9txcfl7Z2obEj7FcXTK9h6TboCZ+L2qzs2cLos1OvZba7bbnUrudOVPO22ytNCkiIiIiEqNQu5SIiIiIiLQLKrhFRERERGKkgltEREREJEYquEVEREREYqSCW0REREQkRiq4pd0zxmxrjBmTur2zMebvSWcSEZHs1G5Le6NpAaXdM8bsATxqrd0v6SwiItI8tdvS3nRKOoBIDlwB9DPGLMRNZt/fWrufMeY04GigI7AfcA2wNXAKUAMcbq2tMsb0wy1AsQOwHviZtfatXL8JEZEConZb2hV1KZFCcBGw1Fo7CPhVg8f2A0YDQ4FLgfXW2q8DLwKnpp7zF+Dn1tpvABcA03IRWkSkgKndlnZFZ7il0D1rrV0LrDXGfArMTm1fBHzNGLMNcDBwnzGm9meKch9TRERS1G5L3lHBLYWuJu325rT7m3H/PjoAa1JnWUREJHlqtyXvqEuJFIK1QPfW/KC1thp41xhzAoBx9o8ynIiINKJ2W9oVFdzS7llrPwH+ZYxZDFzdil38GDjDGPMasAQYGWU+ERGpT+22tDeaFlBEREREJEY6wy0iIiIiEiMV3CIiIiIiMVLBLSIiIiISIxXcIiIiIiIxUsEtIiIiIhIjFdwiIiIiIjFSwS0iIiIiEiMV3CIiIiIiMfr/Y4WnQ89p6f8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#PID controller\n",
    "def StateSpacePID_ode(x, t, A, B, C, u, y_ref, k_p, k_i, k_d):\n",
    "    state = x[:-1]\n",
    "    y = np.dot(C, state)\n",
    "    er = y_ref - y\n",
    "    d_int_e = er\n",
    "    u_pid = k_p * er + k_i * x[-1] - k_d * np.dot(C, np.dot(A, state))\n",
    "    dstate = StateSpace(state, t, A, B, u_pid)\n",
    "    return np.concatenate((dstate, np.array([d_int_e])))\n",
    "\n",
    "def StateSpacePID(x0, t, A, B, C, u, y_ref, k_p, k_i, k_d):\n",
    "    x0_ext = np.concatenate((x0, [0])) # add aditional variable for integral term\n",
    "    solution = odeint(StateSpacePID_ode, x0_ext, t, args=(MatA, MatB, MatC, F, 1, k_p, k_i, k_d))\n",
    "    return solution\n",
    "\n",
    "A = np.array(jacobian_A.subs({M:2.5, m:0.2, b:0.1,l:0.7, I:0.006, g:9.81}).evalf()).astype(float)\n",
    "B = np.array(jacobian_B.subs({M:2.5, m:0.2, b:0.1,l:0.7, I:0.006, g:9.81}).evalf()).astype(float).reshape(4,1)\n",
    "C = np.array([[1,0,0,0],[0,0,1,0]]) \n",
    "D = np.array([0])\n",
    "\n",
    "x0 = np.array([0,\n",
    "               0,\n",
    "               0,\n",
    "               0])  # initial state\n",
    "\n",
    "MatC = [0,0,1,0] # Let us assume that we measure only theta\n",
    "y_ref = 1 # and we want to track constant reference equal to 1 rad. The latter correspond a state responce if the system\n",
    "\n",
    "k_p = -60;\n",
    "k_i = -20;\n",
    "k_d = -6;\n",
    "\n",
    "t0 = 0 # Initial time \n",
    "tf = 20 # Final time\n",
    "t = np.linspace(t0, tf, 1000) \n",
    "\n",
    "solution = StateSpacePID(x0, t, MatA, MatB, MatC, F, 1, k_p, k_i, k_d)\n",
    "\n",
    "figure(figsize=(12, 5))\n",
    "\n",
    "y = solution[:,2]\n",
    "\n",
    "subplot(1, 2, 1)\n",
    "plot(t, y, linewidth=2.0, color = 'red')\n",
    "grid(color='black', linestyle='--', linewidth=1.0, alpha = 0.7)\n",
    "grid(True)\n",
    "xlabel('time')\n",
    "xlim([t0, tf])\n",
    "ylabel(r'Angle ')\n",
    "\n",
    "y = solution[:,0]\n",
    "subplot(1, 2, 2)\n",
    "plot(t, y, linewidth=2.0, color = 'red')\n",
    "grid(color='black', linestyle='--', linewidth=1.0, alpha = 0.7)\n",
    "grid(True)\n",
    "xlim([t0, tf])\n",
    "xlabel('time')\n",
    "ylabel(r'Position ')\n",
    "\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "0585be4b",
   "metadata": {},
   "source": [
    "## EX 2:  Mass-spring damper system\n",
    "\n",
    "Let us consider mass-spring damper system\n",
    "\n",
    "![TP1_SDSystem.png](attachment:TP1_SDSystem.png)\n",
    "\n",
    "with the following system parameters:\n",
    "\n",
    "      mass  m = 1.0 kg\n",
    "\n",
    "      spring constant k = 5.0 N/m\n",
    "\n",
    "      damping constant $\\rho$ = 2 Ns/m\n",
    "\n",
    "Let us suppose that measured output of the system is a position of the mass.\n",
    "\n",
    "## TODO\n",
    "\n",
    "Design a state feedback controller which satisfies the following time domain specifications for a step response ($y_{ref}(t) = 1, t>0$ $y_{ref}(t) = 0, t<0$)\n",
    "\n",
    "    Rise time < 5 s\n",
    "    \n",
    "    Overshoot < 10%\n",
    "    \n",
    "    Steady-state error < 2%\n",
    "    \n",
    "1) Is system controllable?    \n",
    "    \n",
    "2) Choose eigenvalues which ensure the disered behavior of closed-loop system\n",
    "\n",
    "3) Place them\n",
    "\n",
    "4) Add integral term, if needed\n",
    "\n",
    "Remark: if you want to avoid manual tuning (for chosing eigenvalues) you might read this additional material about time domain perfomance of second order system (which is our system)\n",
    "https://www.tutorialspoint.com/control_systems/control_systems_time_domain_specifications.htm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e5d3df0d",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}