alchemyst/Skogestad-Python

{
  "cells": [
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "from __future__ import print_function\n",
        "\n",
        "import numpy as np\n",
        "import scipy as sp\n",
        "import matplotlib.pyplot as plt\n",
        "%matplotlib inline\n",
        "\n",
        "\n",
        "# TODO: This should be reworked!\n",
        "# This version is not really a rework, but provides clearer output\n",
        "# Calculation reworked for transfer functions of any order\n",
        "\n",
        "# Input transfer function and Proportional controller range\n",
        "Gk = 3.0\n",
        "Gz = [[-2, 1]]\n",
        "Gp = [[10, 1], [5, 1]]\n",
        "Krange = np.linspace(-2, 3, 1000)\n",
        "\n",
        "# Check for stability over given range\n",
        "poly_zeros = 1\n",
        "for i in Gz:\n",
        "    poly_zeros *= sp.poly1d(i)\n",
        "poly_poles = 1\n",
        "for i in Gp:\n",
        "    poly_poles *= sp.poly1d(i)\n",
        "\n",
        "unstable_vec = []\n",
        "for K in Krange:\n",
        "    poly_char = float(K)*Gk*poly_zeros + poly_poles\n",
        "    unstable_vec.append(\n",
        "        any(np.real(r) > 0 for r in np.roots(poly_char.coeffs)))\n",
        "\n",
        "# Output stability margin\n",
        "trigger = 0\n",
        "for i in range(0, len(Krange) - 1):\n",
        "    if unstable_vec[i]:\n",
        "        if unstable_vec[i] != unstable_vec[i - 1]:\n",
        "            print('Change from stable to unstable at Kc = ' +\n",
        "                  '%.2f' % Krange[i])\n",
        "            Limit1 = Krange[i]\n",
        "            trigger = 1\n",
        "    else:\n",
        "        if unstable_vec[i] != unstable_vec[i - 1]:\n",
        "            print('Change from unstable to stable at Kc = ' +\n",
        "                  '%.2f' % Krange[i])\n",
        "            Limit2 = Krange[i]\n",
        "            trigger = 1\n",
        "\n",
        "if trigger == 0:\n",
        "    print('No stability margins could be found')\n",
        "else:\n",
        "    print('Stable between Kc = ' + '%.2f' % Limit1 +\n",
        "          ' and Kc = ' + '%.2f' % Limit2)\n",
        "\n",
        "plt.plot(Krange, unstable_vec, 'rD', label='0=unstable\\n1=stable')\n",
        "plt.xlabel('Kc')\n",
        "plt.ylabel('Stability')\n",
        "plt.legend()\n",
        "plt.show()\n"
      ],
      "outputs": [],
      "execution_count": null
    }
  ],
  "metadata": {
    "anaconda-cloud": {},
    "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.6.1"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 1
}