{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#A script to illustrate the incubator P/PI/PID control computations for Section 5.6.\n",
    "\n",
    "#System/Plant Model: The uncontrolled incubator temperature is governed by the\n",
    "#Newton cooling ODE y'(t) = -k*(y(t) - a(t)) where \"k\" is the cooling constant\n",
    "#and \"a(t)\" is ambient temperature. We let us take, for the moment,\n",
    "k = 0.05;\n",
    "var('t');\n",
    "a(t) = 0;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#The control function will be u(t), and the constant \"K\" in the controlled ODE\n",
    "#(equation (5.108) in the text) will be K = 1. The controlled ODE is thus\n",
    "#y'(t) = -k*(y(t)-a(t)) + K*u(t) with\n",
    "K = 1;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#The desired temperature (setpoint) will be 0 degrees for t < 20 and then 3\n",
    "#degrees for t > 20.\n",
    "r(t) = 3*heaviside(t-20);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Assume the initial condition is y(0) = y0 with\n",
    "y0 = 5;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#The Control: Choose the control gains for PID control. In this case we will use\n",
    "#PI control (so Kd = 0):\n",
    "Kp = 1; #Proportional gain\n",
    "Ki = 1/10; #Integral gain\n",
    "Kd = 0; #Derivative gain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#The plant transfer function Gp(s) and controller transfer Gc(s) are, from (5.111)\n",
    "#and (5.129) respectively,\n",
    "var('s');\n",
    "Gp(s) = K/(s+k); print(Gp(s))\n",
    "Gc(s) = Kp + Ki/s + Kd*s; print(Gc(s))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#From (5.118) the closed-loop transfer function is\n",
    "G(s) = Gp(s)*Gc(s)/(1+Gp(s)*Gc(s)); print(G(s))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#The System Response to a Disturbance: This system starts off at the wrong\n",
    "#temperature (5 degrees instead of the desired setpoint 0) and there is a\n",
    "#disturbance in the form of an abrupt setpoint change at time t = 20.\n",
    "\n",
    "#The governing ODE is y'(t) = -k(y(t) -a(t)) + K*u(t) where u(t) = r(t)  - y(t).\n",
    "#According to equation (5.133) the system response in the s domain can be\n",
    "#computed as Y(s) = G(s)*R(s) + y0*Gp(s)/(1+Gp(s)*Gc(s)). Using Sage:\n",
    "R = laplace(r(t),t,s,algorithm='sympy');\n",
    "Y = simplify(G(s)*R[0] + y0*Gp(s)/(1+Gp(s)*Gc(s)));\n",
    "print(Y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#The time domain response is then\n",
    "ysol(t) = inverse_laplace(Y,s,t,algorithm='sympy');\n",
    "print(ysol);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#A plot\n",
    "plt1 = plot(ysol(t), (t,0,100),color='red');\n",
    "plt2 = plot(r(t), (t,0,100),color='blue');\n",
    "show(plt1+plt2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Note the system starts off at 5 degrees (where the desired setpoint is zero),\n",
    "#so the control cools the system to 0 degrees, with a bit of overshoot. The\n",
    "#controller effectively responds to the change in the setpoint at time t = 20."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 9.2",
   "language": "sage",
   "name": "sagemath"
  },
  "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}