{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction to some Python control commands\n",
    "Requires `numpy`, `scipy`, `control`, `jupyter`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "from __future__ import division\n",
    "import numpy as np\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a href='https://python-control.readthedocs.io/en/0.8.2/'>python control tool box</a>\n",
    "\n",
    "<a href='https://python-control.readthedocs.io/en/0.7.0/matlab_strings.html'>conversion from Matlab to Python</a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "import control\n",
    "import control.matlab as mlab"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Let's randomly generate some systems and explore Chaps 1-4 in [JH]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A linear (time invariant) system in state-state space form looks like the following:\n",
    "\\begin{align} \n",
    "\\dot{x}(t)&=Ax(t)+Bu(t)\\\\\n",
    "y(t)&=Cx(t)+Du(t)\n",
    "\\end{align}\n",
    "where $x\\in \\mathbb{R}^{n}$, $y\\in \\mathbb{R}^p$, $u\\in \\mathbb{R}^m$, $A\\in \\mathbb{R}^{n\\times n}$, $B\\in \\mathbb{R}^{n\\times m}$, $C\\in \\mathbb{R}^{p\\times n}$, and $D\\in \\mathbb{R}^{p\\times m}$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A :\n",
      " [[0.77132064 0.02075195]\n",
      " [0.63364823 0.74880388]]\n",
      "B :\n",
      " [[0.49850701]\n",
      " [0.22479665]]\n",
      "C :\n",
      " [[0.19806286 0.76053071]]\n",
      "D :\n",
      " [[0.16911084]]\n"
     ]
    }
   ],
   "source": [
    "# Randomly generate a system\n",
    "n=2; m=1; p=1\n",
    "np.random.seed(10)\n",
    "A=np.random.rand(n,n)\n",
    "B=np.random.rand(n,m)\n",
    "C=np.random.rand(p,n)\n",
    "D=np.random.rand(p,m)\n",
    "sys1=control.StateSpace(A,B,C,D)\n",
    "print('A :\\n',sys1.A)\n",
    "print('B :\\n', sys1.B)\n",
    "print('C :\\n', sys1.C)\n",
    "print('D :\\n', sys1.D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A :\n",
      " [[0.22199317 0.87073231]\n",
      " [0.20671916 0.91861091]]\n",
      "B :\n",
      " [[0.48841119]\n",
      " [0.61174386]]\n",
      "C :\n",
      " [[0.76590786 0.51841799]]\n",
      "D :\n",
      " [[0.2968005]]\n"
     ]
    }
   ],
   "source": [
    "# randomly generate another system\n",
    "n=2; m=1; p=1 # system dimensions\n",
    "np.random.seed(5) # change the seed if you want a different system\n",
    "A=np.random.rand(n,n)\n",
    "B=np.random.rand(n,m)\n",
    "C=np.random.rand(p,n)\n",
    "D=np.random.rand(p,m)\n",
    "sys2=control.StateSpace(A,B,C,D)\n",
    "print('A :\\n',sys2.A)\n",
    "print('B :\\n', sys2.B)\n",
    "print('C :\\n', sys2.C)\n",
    "print('D :\\n', sys2.D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "A = [[0.77132064 0.02075195 0.         0.        ]\n",
       " [0.63364823 0.74880388 0.         0.        ]\n",
       " [0.09673612 0.37145171 0.22199317 0.87073231]\n",
       " [0.12116374 0.46525    0.20671916 0.91861091]]\n",
       "\n",
       "B = [[0.49850701]\n",
       " [0.22479665]\n",
       " [0.08259562]\n",
       " [0.10345252]]\n",
       "\n",
       "C = [[0.05878516 0.2257259  0.76590786 0.51841799]]\n",
       "\n",
       "D = [[0.05019218]]"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# creating a system in series\n",
    "sys=control.series(sys1,sys2) \n",
    "sys=sys2*sys1 # also creates a system in seiries\n",
    "sys"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "A = [[0.77132064 0.02075195 0.         0.        ]\n",
       " [0.63364823 0.74880388 0.         0.        ]\n",
       " [0.         0.         0.22199317 0.87073231]\n",
       " [0.         0.         0.20671916 0.91861091]]\n",
       "\n",
       "B = [[0.49850701]\n",
       " [0.22479665]\n",
       " [0.48841119]\n",
       " [0.61174386]]\n",
       "\n",
       "C = [[0.19806286 0.76053071 0.76590786 0.51841799]]\n",
       "\n",
       "D = [[0.46591134]]"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# system in parallel\n",
    "sys=control.parallel(sys1,sys2)\n",
    "# alt. you can use sys=sys1+sys2\n",
    "sys"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "A = [[0.77132064 0.02075195 0.         0.        ]\n",
       " [0.63364823 0.74880388 0.         0.        ]\n",
       " [0.         0.         0.22199317 0.87073231]\n",
       " [0.         0.         0.20671916 0.91861091]]\n",
       "\n",
       "B = [[0.49850701 0.        ]\n",
       " [0.22479665 0.        ]\n",
       " [0.         0.48841119]\n",
       " [0.         0.61174386]]\n",
       "\n",
       "C = [[0.19806286 0.76053071 0.         0.        ]\n",
       " [0.         0.         0.76590786 0.51841799]]\n",
       "\n",
       "D = [[0.16911084 0.        ]\n",
       " [0.         0.2968005 ]]"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# appends two systems\n",
    "sys=control.append(sys1,sys2)\n",
    "sys"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "A = [[0.80217406 0.13922428 0.40198704 0.27209189]\n",
       " [0.64756127 0.80222777 0.18127195 0.12269706]\n",
       " [0.1018481  0.39108091 0.28859679 0.91581412]\n",
       " [0.12756659 0.48983593 0.29014139 0.97507669]]\n",
       "\n",
       "B = [[0.5248504 ]\n",
       " [0.23667593]\n",
       " [0.08696035]\n",
       " [0.10891942]]\n",
       "\n",
       "C = [[0.20852941 0.80072063 0.13636792 0.09230299]]\n",
       "\n",
       "D = [[0.17804743]]"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# creates a feedback loop with two system, +1 positive feedback, default is -1\n",
    "sys=control.feedback(sys1,sys2,+1)\n",
    "sys"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a state space model\n",
    "ss_sys=control.ss(A,B,C,D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\frac{0.2968 s^2 + 0.3527 s + 0.05338}{s^2 - 1.141 s + 0.02393}$$"
      ],
      "text/plain": [
       "\n",
       "0.2968 s^2 + 0.3527 s + 0.05338\n",
       "-------------------------------\n",
       "    s^2 - 1.141 s + 0.02393"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# transfer function for state-space model\n",
    "control.tf(ss_sys)\n",
    "#mlab.tf(ss_sys) # also works"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Comments on discrete time transfer functions\n",
    "The discrete version has not been implemented in the control toolbox \n",
    "however there are many tools for discrete, linear constant coefficient \n",
    "difference equations in numpy and scipy\n",
    "\n",
    "e.g., `scipy.signal.TransferFunction` ([docs](https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.TransferFunction.html))\n",
    "\n",
    "[docs for `scipy.signal` package](https://docs.scipy.org/doc/scipy/reference/signal.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Characteristic polynomial\n",
    "[docs for `np.poly()`](https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.poly.html)\n",
    "returns coefficients for characteristic polynomial if array in input:\n",
    "\n",
    "`c[0] * x**(N) + c[1] * x**(N-1) + ... + c[N-1] * x + c[N]`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1.        , -1.14060408,  0.0239283 ])"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.poly(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "eigenvalues :\n",
      " [0.02137935 1.11922473]\n"
     ]
    }
   ],
   "source": [
    "# eigenvalues and eigenvectors of A\n",
    "e,v=np.linalg.eig(A)\n",
    "print('eigenvalues :\\n', e)"
   ]
  },
  {
   "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
