Monday, April 11, 2011

LTI Transient-Response Analysis using Python(numpy, scipy, matplotlib)

In undergraduate Engineering Physics, i learn how to modelling a system by using a mathematical model of the system, to analyze the character of the system.

This script can be used to generate plot as in the Book of Modern Control Engineering 4th edition, International Edition, by Katsuhiko Ogata, Prentice Hall

Once when taking a subject Automatic Control is often to make a plot of LTI system, that generated in s(signal) domain by transforming mathematical model of the system using Laplace Transformation.

This example i took from page 307 of that book.
To plot transfer function =
(6.3223 s^2 + 18 s +12.811) / (s^4 + 6 s^3 + 11.3223 s^2 + 18s + 12.811)

The plot generated by those script as below


  1. I want to learn PID design. For that ur blog is a good start point. I am trying to run ur code. I m getting following error:

    Traceback (most recent call last):
    File "D:\Project (All)\Projects\MagLev\PID_Analysis_Python\", line 15, in
    t, s = step(tf)
    File "C:\Python26\lib\site-packages\scipy\signal\", line 716, in step
    vals = lsim(sys, U, T, X0=X0)
    File "C:\Python26\lib\site-packages\scipy\signal\", line 490, in lsim
    GT = dot(dot(vti, diag(numpy.exp(dt * lam))), vt).astype(xout.dtype)
    File "C:\Python26\lib\", line 29, in _show_warning
    file.write(formatwarning(message, category, filename, lineno, line))
    TypeError: idle_formatwarning_subproc() takes exactly 4 arguments (5 given)

    I am new to python .. but took care of the dependencies. Need help.

    1. please refer to documentation your current version scipy.

    2. K corrected it. Thanks a lot. :)

  2. .. And also if I am not jumping how to start with a matrix of transfer functions, say:

    num1 = [-2.3833x10^3]
    den1 = [1,160.3460,-1.9620x10^3,-3.14598852x10^5]

    num2 = [31.9361,0,1.1132x10^5]
    den2 = [1,160.3460,-1.9620x10^3,-3.14598852x10^5]

    matrix0 = [0, lti(num1/den1); 1, lti(num2/den2)]

    Any clue will be helpful.

    1. because you put that on a list, not an array type.
      a = [1.923] # it's a list
      b = numpy.array([0.839]) # it's an array
      you can start by learning numpy.

    2. After some trail n errors seeing the numpy manual I tried as follows:

      num = np.array([[[0.0], [1.0]],[[-2383.3 ],[31.9361,0,111320.0 ]]])
      den = np.array([[[0.0], [1.0]],[[1.0,160.3460,-1962.0,-314598.852],[1.0,160.3460,-1962.0,-314598.852]]])
      tf = lti(num,den)

      But its giving error as

      Traceback (most recent call last):
      File "D:\Project (All)\Projects\MagLev\PID_Analysis_Python\", line 36, in
      tf = lti(num,den)
      File "C:\Python26\lib\site-packages\scipy\signal\", line 236, in __init__
      self.__dict__['num'], self.__dict__['den'] = normalize(*args)
      File "C:\Python26\lib\site-packages\scipy\signal\", line 276, in normalize
      raise ValueError("Denominator polynomial must be rank-1 array.")
      ValueError: Denominator polynomial must be rank-1 array.

      What can be done!

    3. please take a look example from Ogata Modern Control Engineering
      4th edition, International Edition page 307, and my code above is the way to implement it using python.