An experiment to see how `SymPy` could be used to analyze block diagrams in order to compute transfer functions (Laplace or frequency domain analysis).

sympy transfer function solving.ipynb

{
 "metadata": {
  "name": "sympy transfer function solving"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Solve for the Input-Output Transfer Function in a block diagram"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "An experiment to see how `SymPy` could be used to analyze block diagrams\n",
      "\n",
      "Pierre Haessig -- September 26th 2013"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from IPython.core.display import display\n",
      "import sympy\n",
      "from sympy import symbols, Eq, srepr"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 22
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "sympy.init_printing()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "**Example diagram:**\n",
      "\n",
      "A simple closed loop system (\"w0\" to \"w3\" the wire names):\n",
      "\n",
      "    U --w0--> (+/-) --w1--> |C(s)| --w2--+--> Y\n",
      "                ^                        |\n",
      "                +--w3------ |D(s)| <-----+\n",
      "\n",
      "The objective is to compute the Input-Output transfer $Y/U$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "**Create the symbols used in the block diagram**"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# All the wires (signal carriers)\n",
      "W = symbols('W_0:4', seq=True)\n",
      "# Inputs\n",
      "In = symbols('U(s)', seq=True)\n",
      "# Outputs\n",
      "Out = symbols('Y(s)', seq=True)\n",
      "# Tranfer function blocks\n",
      "C, D = symbols('C(s) D(s)')\n",
      "\n",
      "display(In[0], Out[0], [C,D], W)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "latex": [
        "$$U(s)$$"
       ],
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAACkAAAAWCAYAAABdTLWOAAAABHNCSVQICAgIfAhkiAAAAhpJREFU\nSInt1suLT2EYB/DPaFwbRW4l5JKNhdkwu9FYsJgk5Q9wWSALFnJJ2VlMuZWtNGooSblmpzHFZogV\nasgl99uQpDHD/Cze9+j0m3MOfr/fWOC7eU7f5/s877fn9D69/INoxxs0VFC7CCWsLxK1RNGVAs3s\nqHmUkWvEN+z6fX8/cBHPMC5P0KI6k+fxSWVTTNAc+2/LE7So3OQsYYrHK/cH6vAQPWlyRJVNE6yN\nvU7m5FfiMl7gC56jC5vLdCWcwnxhqjU1uSwecC0jtwHnsAAXcACXMBbrMvRXY1yeEPU1MDgai3Ef\n7zPyG9EvXKzXZbnJGfruGJckRC0mOV0w+rRA8xUDGfzbDO5V1M5NiFqYnBJjb07+hLBS7uAQVqVq\n8tCb1qRNDmZw5UhygymuFGNdTs1BrMFjbMEZYVqdwgIvOmcIGuOBtwtMNkXNzRQ3J3KdBXUJJqAV\nR4SV9U72VAfwJKvBGPRFwaScQ7ZGQ+0pblSsu/cLJtM4GnutLuOnRb4rIdJj7RP2XD32Gfr7ZmB7\n/D6W4vuFGzkPEzPMLM3oBVNj/FzGN8XYJQeThN9dwl0cxl504GPk2zLq9sTciozcB+Hmn8Z+YU92\nR/0NjCzTt8VcswI0YDeuR2MDeCks4tacmpnCmunIyG0SLssDYWq9uIUdGF+mrYu6HsOEs/7AA6Na\nLBRu7M4qevz0qVYLDPuj9z/+enwHTzOEP+dtjkMAAAAASUVORK5CYII=\n",
       "text": [
        "U(s)"
       ]
      },
      {
       "latex": [
        "$$Y(s)$$"
       ],
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAACcAAAAWCAYAAABDhYU9AAAABHNCSVQICAgIfAhkiAAAAhFJREFU\nSInt1juMTVEUBuDvTjAyQYhBRYhCFB6FO82YQhQSQSR6hUJQjgSVRkg0JNNIPKKg9IqIxmMy0Uw8\nIpEokDCJx0wIlRFmcBV7X27O7HPudS8K/M1J1v+ftf+z1l57H/4BnMYbTGvi3VWoYFue4EIU9BYk\nKWMczzCjJr4CX7CvCWNVXMFLdKTI2XiFj1iW4DvwKJpYneEu473mqlZFj1Cc3XmCdfiKB2jPcMfi\nywcz8QWC4bMtGIOS0JHHRaK+aOJITWx9jN3F5Ix+f+Q25OTbhBsYxiehOwPYldAejrl68sxNxUOh\ngmvRiRGMYklCfytqZyW47XGxYRzHIZzEbdxJ6DdG/YE8c7BS+MoXuBpf2JnQtQt7NK8V92KeuQmu\nMxGbF9caKDIHe6KwIkxSCosif7PA3Kh0VfMwhuf1RG1COypYmqPpivy5HL7Xj7YexWbMqbPuiNCN\nuhiKyRfm8OXIny/IsRWDwkRXhP3ZLxy8Kbz+Veaqbe1vINdMYepPRKNvpas4roG2NmJuivCVTxpJ\nVoNTMe+WTHzCQLT9ZOJajAnHwmLpTb9GOFyzqE7vh0y8Kz6/m5vUgjm4Jhya3SZO9UXhWhsUulCK\n2rIwydcz+u6anHUxpLitMB+fcSbB7YgGnwpVeof7wjE1PaMtRV3h9dUMLvkDF3+zWC5M4N4WchT+\nMrWK3/qz+R9/Fb4BipyLyhvaq/8AAAAASUVORK5CYII=\n",
       "text": [
        "Y(s)"
       ]
      },
      {
       "latex": [
        "$$\\begin{bmatrix}C(s), & D(s)\\end{bmatrix}$$"
       ],
       "output_type": "display_data",
       "text": [
        "[C(s), D(s)]"
       ]
      },
      {
       "latex": [
        "$$\\begin{pmatrix}W_{0}, & W_{1}, & W_{2}, & W_{3}\\end{pmatrix}$$"
       ],
       "output_type": "display_data",
       "text": [
        "(W\u2080, W\u2081, W\u2082, W\u2083)"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "1) Generate the exhaustive set of equations of the diagram"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The structure of the equations of a block diagram is *linear*:\n",
      "\n",
      "$$(Wires[:], Out[:])^T = M \\times (Wires[:], In[:])^T$$\n",
      "\n",
      "where the matrix $M$ is a very sparse matrix and contains $1, -1$ for summation blocks and Transfer functions ($H(s)$, ...)"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "For now, the equations are computed by hand:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Diagram_eqs = [\n",
      " Eq(W[0], In[0]),\n",
      " Eq(W[1], W[0]-W[3]),\n",
      " Eq(W[2], C*W[1]),\n",
      " Eq(W[3], D*W[2]),\n",
      " Eq(Out[0], W[2])\n",
      "]\n",
      "Diagram_eqs"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "latex": [
        "$$\\begin{bmatrix}W_{0} = U(s), & W_{1} = W_{0} - W_{3}, & W_{2} = C(s) W_{1}, & W_{3} = D(s) W_{2}, & Y(s) = W_{2}\\end{bmatrix}$$"
       ],
       "output_type": "pyout",
       "prompt_number": 4,
       "text": [
        "[W\u2080 = U(s), W\u2081 = W\u2080 - W\u2083, W\u2082 = C(s)\u22c5W\u2081, W\u2083 = D(s)\u22c5W\u2082, Y(s) = W\u2082]"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "**Next question** : generates this set of equations *automatically*..."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "2) Solving the equations"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Solving the equations involves solving for $Wires$ and $Out$ and thus eliminating all $Wires$ variables:\n",
      "\n",
      "$$(Wires[:], Out[:]) = f(In[:])$$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Concatenation of symbols:\n",
      "Out + W"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "latex": [
        "$$\\begin{pmatrix}Y(s), & W_{0}, & W_{1}, & W_{2}, & W_{3}\\end{pmatrix}$$"
       ],
       "output_type": "pyout",
       "prompt_number": 15,
       "text": [
        "(Y(s), W\u2080, W\u2081, W\u2082, W\u2083)"
       ]
      }
     ],
     "prompt_number": 15
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "sol = sympy.solve(Diagram_eqs, Out + W)\n",
      "sol"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 8,
       "text": [
        "\u23a7                   U(s)            C(s)\u22c5U(s)        C(s)\u22c5D(s)\u22c5U(s)          C\n",
        "\u23a8W\u2080: U(s), W\u2081: \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500, W\u2082: \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500, W\u2083: \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500, Y(s): \u2500\u2500\u2500\n",
        "\u23a9              C(s)\u22c5D(s) + 1      C(s)\u22c5D(s) + 1      C(s)\u22c5D(s) + 1         C(s\n",
        "\n",
        "(s)\u22c5U(s)  \u23ab\n",
        "\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u23ac\n",
        ")\u22c5D(s) + 1\u23ad"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The time it takes: about **50 ms** for this simple diagram. Question: how does it *scale* with the *size* of the diagram ?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit sympy.solve(Diagram_eqs, Out + W)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "10 loops, best of 3: 41.6 ms per loop\n"
       ]
      }
     ],
     "prompt_number": 25
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Out of all the resulting variables, we are only interested those of $Out$. We therefore select the interesting solution :"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Y = Out[0]\n",
      "Y_sol = sol[Y]\n",
      "Y_sol"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "latex": [
        "$$\\frac{C(s) U(s)}{C(s) D(s) + 1}$$"
       ],
       "output_type": "pyout",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAFEAAAAlCAYAAAA6PAXhAAAABHNCSVQICAgIfAhkiAAAA51JREFU\naIHt2UmIHUUcx/HPRB1Bw2TwEhfQKBEXgomICoqDgqKOIIrLxYvJTfSgOYqgghvqQVHBY+PBo2Bc\n4gI6qIgGRZC4IQoOJJJEJNGoaCLj4d9vXk+n3+tlXnqe0D9o+nVXdX2ra6r/9a/f0Ol/oTXYg/UV\n67+Ce1viLZc1Eq3Fs/gBf2MXtmM2U+cRvFyjzYuwFyem13N4vqDenThYcL8OL89qXevEoH2L23EO\nzsM9mE/rTIpOXlOz7a+wOf09p/ogNuFlWa3rTTGIqwvKptPzLH7HMbnyGXwiBuEAdmBDpvwxvJP+\nnlN9EIt4dViNtKrhcyfhOryg+JPan55n8AX+zZQdi1fxETbiUjyTq7MDl6d16yjPO5qsRTV9cD0m\n8E1JvbP1P+2epsRMfU3EUiIkZDWPE3BazX7leXVZP9XkoflMnKhYb8qRM/VXJHgbb2ArTs/V+S09\nr6nZrzzvaLIW1XQQv8eCWEiG6YDimLlZfFof4EZ8h2sz5VPpeb94yaIXnE7bL+PVYbWu7dht+MLy\nhOh8lbayacnN+EOEm6fECpqf/S/hvdy9KrxhrNZ1Fn4WMeY2keKci7v049INjlwtzxQvexnOwFVi\nlX8gU+dxvJvh/CVW6I0p5z4cwvW5PuV5dVkrolPwHH4UyfZukS7ckpZPYp+ledtasVPYlT4zjydx\nXKbOTmzJXF8s4toe8cl+ipsK+pPnNWGNpR5Vf8eyT3GYGDVvuazWNK3dvXMd3ljsnTt16jQ26uVe\nCyvai06dmmicnepR8JatNpzqRISPBbHr2Iv3cbelSfCoeMvRDLaJcVgQHuZQrdOOU52IbdbJwn7a\nJNyVX4R5mn35UfDKlOChAWWzwrS9FX+qMIhtOdUJXi9gbMA/eDj3EnWd6jyvTInBg5jVQQMGsWeF\njYNTvRNv6e+7i3hVWFV5I1MPMi5O9de4egivCivPa+RW11FvJo6LUz1hac7axKku492fttk77ii4\nd0VJP5eoN4htOtXDdL6w1YbxylhlvBfFYtY7thXc+6yknwPVllOdGLywHMKDmXtNnOo8r0yJES0s\nRJ42If4Keaf6y7TOh7hQuVN9gYhvPV2Cj3E4vT5epDinigViq/j/8ud4OvNcnleFVcRrqtX6s3OV\nCB2bFIeQRbXhVCf6yfZhkR/OiXx0MtefJk51nlemxOCZeGWmr9kjqdj2QI2zUz0KXisaZ6d6FLxO\nnTp16tRpZfUf62xhkNFqr2wAAAAASUVORK5CYII=\n",
       "prompt_number": 13,
       "text": [
        "  C(s)\u22c5U(s)  \n",
        "\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n",
        "C(s)\u22c5D(s) + 1"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Find the Transfer function $Y(s)/U(s)$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "U = In[0]\n",
      "H_Y = Y_sol/U\n",
      "H_Y"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "latex": [
        "$$\\frac{C(s)}{C(s) D(s) + 1}$$"
       ],
       "output_type": "pyout",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAFEAAAAlCAYAAAA6PAXhAAAABHNCSVQICAgIfAhkiAAAAzJJREFU\naIHt2U2IVWUcx/HPmI5QMoobzcCMjCyGmhAUDIeCophFFL0saqPuwhbi0o0FFRIulBBanlVLwell\nelkkFRFiBGFWREEDCqaU2iSkU9PiOad77unce86Md5656vnC4cB5/uf5/e+f5zwvv0vDdctynMH6\nmvGHsWv+0uk/VuEgfsJfOIUJjOViXsXbs+hzI37FLT3Ksa9ZJxTtezyHu3EPXsJkGjMoFOTRWfb9\nLbb3JMs+532hiMtK2lak9zH8gZsK7aP4ElO4gGMYzrW/jo96mWw/shL/YE9F3D58Wni2GL9jP+7E\nBjwvjOKMJ/FnGhuVmILrMYDvKuLu0vq0M4aEkfqOMJcSpoQ8k7gZt+GXq8p0liyKqDVQM25I+GTz\n/IYEH+I97MbaQszF9L58jvnNmZhF/BEz2j/BMi4onzO3Y7PwqT+BH/BYrn0ovZ+/ujT7nwmc1n1h\nKZsTO/WV3wY9ZYHmxJgjEXYKn/VxPCtscTbgRXyTxnyGB7SvzncIxd2C2/Ew7sPJXMwmfIHp+Uu/\nf7gVb+JnYbN9WtiaPJ22D+Ks9n3iKuFUcip9ZxJvYEku5gR2zGfi1xqvmf2J5azyaeKGZYXm7NzQ\n0BCf7BQxs6BZNDTMhdiu8zXncsdwnRNh+pjBlbTtE+E0s6Tk/YV0uUcxLtRhBtuqXlgnjuuc4GOs\nFmyrEcGVOSeYrvkfH8PlTvByh7Yxwex9BpfUKGIs1znBuyUaw7iMVwo/oqhXpVXUqyLRuYh5pnQo\nYmZArMTjOOT/Xh4te2kUX+PvXNtiHMHnuF+wqw4UYo7hQd0dlhP4QOsMXaZXR6uuXs/IRGK6zt04\niUe66NXRKurNu8udjcR+cZ0HtO9Zi3p1tKr09qR9ZtcLJc+2VuTZRlbEfnGd7xUssm56VVpVem8J\ni1l2jZc8O16RZ0diuc6JzgvLFezNPaujV9Qq6lWR6NHCQlzXeamwxVkjLBC7cRRfCX+LZhT16miV\n6c2VZVqjc5EwdYwon0L+I4brnGhttqeF/eFRYT86WMinqFdHq6hXRaLzSHwol2v+Smr23ZHYrvN1\n6XLHdp0bl7uhoaGh4cbiX/SvEsXWTwQ8AAAAAElFTkSuQmCC\n",
       "prompt_number": 14,
       "text": [
        "     C(s)    \n",
        "\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n",
        "C(s)\u22c5D(s) + 1"
       ]
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "repr(H_Y)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 16,
       "text": [
        "C(s)/(C(s)*D(s) + 1)"
       ]
      }
     ],
     "prompt_number": 16
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Look at the internals of the Tranfer Function:\n",
      "srepr(H_Y)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 23,
       "text": [
        "Mul(Symbol('C(s)'), Pow(Add(Mul(Symbol('C(s)'), Symbol('D(s)')), Integer(1)), \n",
        "Integer(-1)))"
       ]
      }
     ],
     "prompt_number": 23
    }
   ],
   "metadata": {}
  }
 ]
}