PID_designer_notebook.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__PID controller / compensator designer based on 2nd order approximation "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "from ipywidgets import interact, interactive, fixed, interact_manual\n",
    "import ipywidgets as widgets\n",
    "\n",
    "%matplotlib notebook\n",
    "# %pylab widgetsnbextension\n",
    "import matplotlib.pyplot as plt\n",
    "from numpy import log as ln, arange, roots\n",
    "from math import atan, sqrt, pi, degrees, tan, e\n",
    "# from sympy import symbols, simplify, fraction,solve, pprint, evalf, im, I\n",
    "from math import isclose\n",
    "from control import root_locus, tf, tfdata, feedback\n",
    "import sympy as S\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Run the code below, and you will get an interactive window below to get the Overshoot from the damping ratio and vice versa"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d07abfe8685e4c4aa869fbf5b928df3b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(Text(value='0', description='zeta'), Output()), _dom_classes=('widget-interact',))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "01b810cd45294c00a80fd53360761d3d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(Text(value='1', description='os'), Output()), _dom_classes=('widget-interact',))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.show_os(os)>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "def angle_from_damping_ratio(zeta):\n",
    "    '''the angle needed to meet the zeta condition, drawn from origin to + pole, and by symmetry can be extended to bottom angle\n",
    "    taken by drawing a line from the origin to the complex conjugate pole ie the -pole'''\n",
    "#     try:\n",
    "#         zeta = float(zeta)\n",
    "#     except:\n",
    "#         print(\"Invalid zeta\")\n",
    "#         return \n",
    "    if zeta == 1:\n",
    "        return 0\n",
    "    elif abs(zeta) == 0:\n",
    "        return 90\n",
    "    x = sqrt(1-zeta**2)/zeta\n",
    "    angle = atan(x)\n",
    "    return degrees(angle)\n",
    "\n",
    "def os_from_damping_ratio(zeta):\n",
    "    num = -pi * zeta\n",
    "    den = sqrt(1-zeta**2)\n",
    "    return 100 * e**num/den\n",
    "\n",
    "def damping_ratio_from_os(os):\n",
    "    os = os/100\n",
    "    zeta = -ln(os)/sqrt(pi**2+(ln(os))**2)\n",
    "    return zeta\n",
    "def show_zeta(zeta):\n",
    "    if zeta:\n",
    "        zeta = float(zeta)\n",
    "    else:\n",
    "        zeta = 0\n",
    "    if zeta>1 or zeta <0:\n",
    "        return \"Zeta is between 0 and 1\"\n",
    "    global_zeta = str(zeta)\n",
    "    os = os_from_damping_ratio(zeta)\n",
    "    ang = angle_from_damping_ratio(zeta)\n",
    "    out = {\"Zeta\":zeta, \"OS\":os, \"Angle deg\":ang}\n",
    "    print(out)\n",
    "\n",
    "def show_os(os):\n",
    "    if os:\n",
    "        os = float(os)\n",
    "    else:\n",
    "        os = 1\n",
    "    zeta = damping_ratio_from_os(os)\n",
    "    ang = angle_from_damping_ratio(zeta)\n",
    "    out = {\"Zeta\":zeta, \"OS\":os, \"Angle deg\":ang}\n",
    "    print(out)\n",
    "\n",
    "    \n",
    "interact(show_zeta, zeta = \"0\")\n",
    "interact(show_os, os = \"1\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Find PD Gain and zero location__ Give it an open loop transfer function, and the desired pole location. It will print the KD (deriviative gain), and the zero location in the s-plane. Run the cell right below then below that there's a cell with an example,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_zero_and_gain_for_PD_controller(oltf, pole, rounding_accuracy=3):\n",
    "    ''' takes open loop transfer func,oltf, and desired pole. returns the zero location, as well as\n",
    "     the gain of the PD controller, to meet that desired pole requirement. Assumes pole provided is in the 2nd quadrant, ie, left hand side of the JW axis or on it '''\n",
    "    #sanitize the pole to make sure it is on the top of the LHS ie 2nd quadrant\n",
    "    pole = ((-1* abs(pole.real)) + abs(pole.imag)*1j)\n",
    "    oltf_at_pole = oltf.horner(pole)\n",
    "    zero_location_angle = pi - angle(oltf_at_pole)\n",
    "    sigma = abs(pole.real)\n",
    "    wd = abs(pole.imag)\n",
    "    zero_loction = sigma + wd/tan(zero_location_angle)\n",
    "    kd = abs(1/( oltf_at_pole * (pole+zero_loction) ) )\n",
    "    print({\"kd\": round(kd, rounding_accuracy), \"zero\":round(zero_loction, rounding_accuracy) })"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Example\n",
    "# s = tf('s')\n",
    "#oltf = 1/(s**2 + 4)\n",
    "# desired_pole = (x,y*j)\n",
    "# find_zero_and_gain_for_PD_controller(oltf, desired_pole)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_routh_poly_state_space(a):\n",
    "    #extracts the polynomial from a state space\n",
    "    #det(sI-A) = the poly that u use for routh hurwitz table from state space\n",
    "    a = S.Matrix(a)\n",
    "    e = S.eye(a.rows)\n",
    "    x = S.symbols('s')\n",
    "    e = e*x\n",
    "    poly = S.expand(S.det(e - a))\n",
    "    print(poly)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "s**3 - 6*s**2 - 7*s - 52\n"
     ]
    }
   ],
   "source": [
    "a = [[0, 3, 1], [2, 8, 1], [-10, -5, -2]]\n",
    "find_routh_poly_state_space(a)"
   ]
  },
  {
   "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}