MID_tutorial_interractive.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutoriel MID Générique\n",
    "<u>Objectif</u> : Illustrer les propriétés de la MID.\n",
    "\n",
    "Concrètement, nous considérons un système qui a un controle par une loi de commande retardée (expression mathématique) et l'on prend compte et tirons avantage de ce retard dans le calcul des gains de la loi de commande retardée pour stabiliser le système.\n",
    "\n",
    "<u> Informations importantes </u> :\n",
    "\n",
    "- Stabilité : Lorsque l'on a l'équation décrivant un système dans l'espace de Laplace (transformée de Laplace), la stabilité est traduite par l'ensemble des racines (complexes) de l'équation ayant leurs parties réelles négatives. On dit que les racines sont à gauche (de l'axe des ordonnées).\n",
    "- Transformée de Laplace (juste pour info) : Dans notre cas, la partie la plus importante dans la transformée de Laplace est la dérivation, c'est ce qui nous permet de passer d'une équation différentielle à une équation polynomiale. L'idée très grossière est qu'une dérivée d'ordre $n$ dans l'espace \"classique\"/\"temporel\" correspond à une multiplication par un terme $s^n$ dans l'espace de Laplace (fréquentiel). Donc, par exemple, $4\\times k \\times y''(t)\\rightarrow 4ks^2y(s)$\n",
    "\n",
    "On traite une équation différentielle du type :\n",
    "$$\n",
    "\\frac{d^n y(t)}{dt^n} + \\sum\\limits_{i=0}^{n-1}a_i\\frac{d^i y(t)}{dt^i} + \\sum\\limits_{j=0}^m \\alpha_j\\frac{d^j y(t-\\tau)}{dt^j}=0\n",
    "$$\n",
    "\n",
    "- $n, m\\in\\mathbb{N}, n>m$\n",
    "- $a_i, \\alpha_j \\in \\mathbb{R}$\n",
    "- $s\\in\\mathbb{C}$,\n",
    "- $\\tau \\in \\mathbb{R}$\n",
    "\n",
    "Par une transformée de Laplace :\n",
    "$$\n",
    "Q(s) = s^n+a_{n-1}s^{n-1}+...+a_0 + \\left[\\alpha_ms^m+...+\\alpha_0\\right]e^{-s\\tau}\n",
    "$$\n",
    "\n",
    "Cette équation est polynomiale et non différentielle ce qui la rend beaucoup plus simple à traiter.\n",
    "\n",
    "On souhaite réaliser ce que l'on appelle un placement de pôle. Cela signifie dans ce cas que l'on veut trouver les coefficients ($a_i$ et $\\alpha_j$ parfois notés $b_j$) soient de telle sorte que $Q(s_0)=0$. \n",
    "\n",
    "On connaît donc une valeur de racine à \"forcer\" pour le quasipolynôme en jouant sur les coefficients.\n",
    "\n",
    "Le code suivant a été développé sur la base du pseudo-code de la MID Générique, décrit par Franck, puis adapté en version interactive.\n",
    "\n",
    "Version : 1.5 (01/06/2021) ayrton\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>div.output_scroll { height: 55em; }</style>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from IPython.core.display import display, HTML\n",
    "display(HTML(\"<style>div.output_scroll { height: 55em; }</style>\"))\n",
    "%matplotlib notebook\n",
    "import numpy as np\n",
    "import sympy as sp\n",
    "import cxroots as cx\n",
    "import scipy.special as spsp\n",
    "import matplotlib.pyplot as plt\n",
    "import ipywidgets as widgets\n",
    "from IPython.display import display\n",
    "from IPython.display import Code\n",
    "from IPython.display import Markdown\n",
    "\n",
    "\n",
    "def MID_function(n,m,s0,value_tau):\n",
    "    \n",
    "    s = sp.symbols('s')  # define variable s for our problem to be solved\n",
    "    tau = sp.symbols('tau')  # define variable tau : delay\n",
    "    a = sp.symbols([\"a{:d}\".format(i) for i in range(n)], real=True)\n",
    "    alpha = sp.symbols([\"alpha{:d}\".format(i) for i in range(m + 1)], real=True)\n",
    "    Polynomial = s**n + np.array(a).dot([s**i for i in range(n)]) # Revient à faire s^n + a_{n-1}^{n-1}...\n",
    "    Delayed = np.array(alpha).dot([s**i for i in range(m+1)])*sp.exp(-s*tau) # Revient à faire \n",
    "    Q = Polynomial + Delayed \n",
    "    SysDerivatif = [Q]\n",
    "    for i in range(n+m+1):\n",
    "        DerniereDerivee = SysDerivatif[-1]\n",
    "        SysDerivatif.append(DerniereDerivee.diff(s)) # Dérivée par rapport à s\n",
    "    sol = sp.linsolve(SysDerivatif[:-1], alpha + a).args[0] # Solveur selon les alpha et les a\n",
    "    solNum = sol.subs({s : s0})\n",
    "    solNum = solNum.subs({tau : value_tau})\n",
    "    a_num = list(solNum[m + 1:])\n",
    "    alpha_num = list(solNum[:m + 1])\n",
    "    QNumerique = s**n + np.array(a_num).dot([s**i for i in range(n)])+\\\n",
    "    np.array(alpha_num).dot([s**i for i in range(m+1)])*sp.exp(-s*tau)\n",
    "    QNumerique = QNumerique.subs(tau, value_tau)\n",
    "    sysRootFinding = [QNumerique, QNumerique.diff(s)]\n",
    "    sysFunc = [sp.lambdify(s, i) for i in sysRootFinding]\n",
    "    rect = cx.Rectangle([-100, 10], [-100, 100])\n",
    "    roots = rect.roots(sysFunc[0], sysFunc[1], rootErrTol=1e-5, absTol=1e-5, M = n + m + 1)\n",
    "    xroot = np.real(roots[0])\n",
    "    yroot = np.imag(roots[0])\n",
    "    return xroot,yroot,QNumerique\n",
    "\n",
    "def factorization_integral_sympy(n, m, s0, value_tau):\n",
    "    s = sp.symbols(\"s\", complex = True)\n",
    "    t = sp.symbols(\"t\", real = True)\n",
    "    return value_tau**(m+1) * (s - s0)**(n + m + 1) / sp.factorial(m) * sp.Integral(t**m * (1-t)**n * sp.exp(-t*value_tau*(s - s0)), (t, 0, 1))\n",
    "\n",
    "def factorization_integral_latex(n, m, s0, tau):\n",
    "    factor = str(tau**(m+1) / spsp.factorial(m))\n",
    "    parenthesis = \"(s + \" + str(-s0) + \")\"\n",
    "    power = \"^{\" + str(n + m + 1) + \"}\"\n",
    "    \n",
    "    return r\"\\$\" + \"\\Delta(s) = \" + factor + parenthesis + power + r\"\\int_0^1 t^{\" + str(m) + r\"} (1 - t)^{\" + str(n) + \"} e^{-\" + str(tau) + \"t\" + parenthesis + \"} \\mathrm{d}t\" + \"$\" "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_slider = widgets.Output()\n",
    "output_equation = widgets.Output()\n",
    "output_factorisation_integrale = widgets.Output()\n",
    "\n",
    "slider_n = widgets.IntSlider(min=0,max=3,step=1,description='n :',value=2)\n",
    "slider_m = widgets.IntSlider(min=0,max=3,step=1,description='m :',value=1)\n",
    "FloatText_s0 = widgets.BoundedFloatText(min=-10.0,max=0.0,step=0.1,description='s0 :',value=-3.0,disabled=False)\n",
    "FloatText_tau = widgets.BoundedFloatText(min=0.50,max=2.00,step=0.01,description='tau:',value=0.73,disabled=False)\n",
    "\n",
    "      \n",
    "def common_slider(n,m,s0,value_tau):\n",
    "    output_slider.clear_output()\n",
    "    output_equation.clear_output()\n",
    "    output_factorisation_integrale.clear_output()\n",
    "    \n",
    "    with output_slider:\n",
    "        fig, ax = plt.subplots(figsize=(6, 4))\n",
    "        ax.set_title(\"Root spectrum\")\n",
    "        ax.set_xlabel(r\"Re$(s)$\")\n",
    "        ax.set_ylabel(r\"Im$(s)$\")\n",
    "        ax.axhline(linewidth=2, color='black', zorder = 2)\n",
    "        ax.axvline(linewidth=2, color='black', zorder = 2)\n",
    "        ax.axvspan(0, 1, alpha=0.5, color='red', zorder = 2)\n",
    "        ax.grid()\n",
    "        line, = ax.plot([], [], 'o', color='steelblue', zorder = 3)\n",
    "        \n",
    "        xr,yr,eq = MID_function(n,m,s0,value_tau)\n",
    "        line.set_data(xr, yr)\n",
    "        ax.relim()\n",
    "        ax.autoscale_view()\n",
    "        fig.canvas.draw()\n",
    "        \n",
    "    with output_equation :\n",
    "        print (\"Fp = \",eq)\n",
    "    \n",
    "    with output_factorisation_integrale :\n",
    "        factorization_integral_latex_equation = factorization_integral_latex(n, m, s0, value_tau)\n",
    "        display(Markdown(factorization_integral_latex_equation))\n",
    "        \n",
    "\n",
    "def slider_n_eventhandler(change):\n",
    "    common_slider(change.new, slider_m.value, FloatText_s0.value, FloatText_tau.value)\n",
    "def slider_m_eventhandler(change):\n",
    "    common_slider(slider_n.value, change.new, FloatText_s0.value, FloatText_tau.value)\n",
    "def FloatText_s0_eventhandler(change):\n",
    "    common_slider(slider_n.value,slider_m.value , change.new, FloatText_tau.value)\n",
    "def FloatText_tau_eventhandler(change):\n",
    "    common_slider(slider_n.value,slider_m.value , FloatText_s0.value, change.new)\n",
    "\n",
    "slider_n.observe(slider_n_eventhandler, names='value')\n",
    "slider_m.observe(slider_m_eventhandler, names='value')\n",
    "FloatText_s0.observe(FloatText_s0_eventhandler, names='value')\n",
    "FloatText_tau.observe(FloatText_tau_eventhandler, names='value')\n",
    "\n",
    "input_widgets = widgets.HBox([slider_n, slider_m, FloatText_s0, FloatText_tau])\n",
    "\n",
    "subtab_factorisation_integrale = widgets.VBox([output_factorisation_integrale])\n",
    "\n",
    "tab = widgets.Tab([output_slider, output_equation, subtab_factorisation_integrale])\n",
    "tab.set_title(0, 'Root plot ')\n",
    "tab.set_title(1, 'Output equation')\n",
    "tab.set_title(2, 'Factorized integral equation')\n",
    "\n",
    "description_label = widgets.Label('Insert degree of polynomial n, degree of delay polynomial m, s0 and tau : ')\n",
    "\n",
    "titre_label = widgets.HTML(\"<h3 class = 'text-center'><font color='blue'>Generic MID</font><size='10'>\")\n",
    "\n",
    "dashboard = widgets.VBox([titre_label,description_label,input_widgets, tab])   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ec592383fa6747298432ecb823635475",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(HTML(value=\"<h3 class = 'text-center'><font color='blue'>Generic MID</font><size='10'>\"), Label…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(dashboard)"
   ]
  },
  {
   "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
}