Simulating a sub-system in a Python (non real-time) simulation environment through a transfer function

Question

I am coding a dynamic system simulation (fixed step, non real-time; it runs on my desktop) and I want to model some of the system's components (e.g. filters...) through the tools made available by scipy.signal (e.g. dlsim). For those components, I know their representation in the classical form of transfer functions.

With scipy.signal, it is pretty easy and straightforward to simulate the output of the transfer functions "statically", that is once the time and input vectors are already known; on the other hand, I couldn't find a way to compute it within each simulation step. My simulator also include some closed-loop controllers, thus the outputs change dinamically as the simulation moves forward.

Any ideas?

PS I found this thread which seems to be quite similar, but I must admit that I do not understand the solution given by the author...: How to simulate one step to a transfer function in python

Answer 1

The described solution of the question How to simulate one step to a transfer function in python works as follows. You generate only two simulation steps for the inputs U (array-like) and T (also array-like). With both variables and your system you call the function scipy.signal.lsim (or scipy.signal.dsim for discrete systems) and you set also the initial values for the system states X. As result you get the output values and the new states Xn+1 which you store in the state variable for X.

In the next loop you take the last values of U and T and add the next input and time step. Now you call the lsim again but this time with the states X of the last iteration and so on.

Here some sample code of a Two-Mass-System:

(sorry, it's not beautiful but it works.)

import math
import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt

class TwoMassSystem:
    def __init__(self):
        # Source: Feiler2003 
        JM = 0.166 # kgm^2   moment of inertia (drive-mass)
        JL = 0.333 # kgm^2   moment of inertia (load-mass)
        d  = 0.025 # Nms/rad damping coefficient (elastic shaft)
        c  = 410.0 # NM/rad  stiffness (elastic shaft)

        self.A  = np.array([[-d/JM, -c/JM,  d/JM],
                            [  1.0,   0.0,  -1.0],
                            [ d/JL,  c/JL, -d/JL]] )
        self.B  = np.array([[ 1/JM, 0.0, 0.0],
                            [  0.0, 0.0, 0.0],
                            [  0.0, 0.0,-1/JL]]    )                                                    
        self.C  = np.array([ 0.0, 0.0, 1.0 ])
        self.D  = np.array([ 0.0, 0.0, 0.0 ]       )  
        self.X  = np.array([ 0.0, 0.0, 0.0 ]       ) 
        self.T  = np.array([ 0.0])
        self.U  = np.array([[0.0, 0.0, 0.0]])
        self.sys1 = signal.StateSpace(self.A, self.B, self.C, self.D)

    def resetStates(self):
        self.X  = np.array([ 0.0, 0.0, 0.0 ]       ) 
        self.T  = np.array([ 0.0])
        self.U  = np.array([[0.0, 0.0, 0.0]])      

    def test_sim(self):
        self.resetStates()   
        h = 0.1
        ts = np.arange(0,10,h)
        u = []
        t = []
        for i in ts:
            uM = 1.0
            if i > 1:
                uL = 1.0
            else:
                uL = 0.0
            u.append([ uM,
                      0.0,
                       uL])
            t.append(i)
        
        tout, y, x = signal.lsim(self.sys1, u, t, self.X)
        return t, y

    def test_step(self, uM, uL, tn):
        """
        test_step(uM, uL, tn)

        The call of the object instance simulates the two mass system with
        the given input values for a discrete time step. 
        Parameters
        ----------
        uM : float
            input drive torque 
        uL : float
            input load torque 
        tn : float
            time step
        Returns
        -------
        nM : float
            angular velocity of the drive
        nL : float
            angular velocity of the load
        """
        u_new = [ uM,
                  0.0,
                  uL]

        self.T = np.array([self.T[-1], tn])
        self.U = np.array([self.U[-1], u_new])

        tout, y, x = signal.lsim(self.sys1, self.U, self.T, self.X)
        # x and y contains 2 simulation points the newer one is the correct one.
        self.X = x[-1] # update state      
        return y[-1] #        

if __name__ == "__main__":
    a = TwoMassSystem()
    tsim, ysim = a.test_sim()
    h = 0.1
    ts = np.arange(h,10,h)
    ys = []
    a.resetStates()
    for i in ts:
        uM = 1.0
        if i > 1:
            uL = 1.0
        else:
            uL = 0.0     
        ys.append(a.test_step(uM, uL, i))

    plt.plot(tsim, ysim, ts, ys)
    plt.show()        

However:

However, as you can see in the plot the result isn't identical and here is the problem, because I'dont know why. Therefore, I created a new question: Why does the result of a simulated step differ from the complete simulation?

Sources:

@InProceedings{Feiler2003,
  author    = {Feiler, M. and Westermaier, C. and Schroder, D.},
  booktitle = {Proceedings of 2003 IEEE Conference on Control Applications, 2003. CCA 2003.},
  title     = {Adaptive speed control of a two-mass system},
  year      = {2003},
  pages     = {1112-1117 vol.2},
  volume    = {2},
  doi       = {10.1109/CCA.2003.1223166}
}