How to correctly plot the step response of a MIMO system with python control package

I would need to plot the step responses of a MIMO system with the python control package.

I've tried so far by using the function step_response, that however converts the system into a SISO before computing the step response, so that only one set of output is computed.

I then tried using the function forced_response with different setup for the input (i.e. constant unity value, numpy array of ones etc..., just for the sake of trying). I get different step responses, so related to other output, but not all the responses (i.e. number of input x number of output).

Here is a minimum sample code that implements a simple 2nd order model with 2 input and 4 output and dummy data. In attachment a plot of the responses I get.

stepResponses

In my test I first run the step_response function, yout results to be of size 4 x size_time (so only the first 4 output are excited).

Then I run the forced_response function, and youtForced still results of size 4 x size_time, instead of size 4 x size_time x 2 (or similar) as I expected (in the hypothesis forced_response treats the system as a MIMO).

Is there a way to have full control of the step response via the forced_response function (similarly to what the MATLAB step function does)?

Unfortunately there is poor documentation and very few practical examples about this.

Many thanks to who can help.

from control import ss, step_response, forced_response
import numpy as np
import matplotlib.pyplot as plt

sz = 2

f1 = 1*2*np.pi
f2 = 1.5*2*np.pi
OM2 = [-f1**2, -f2**2]
ZI = [-2*f1*0.01, -2*f2*0.01]

A11 = np.zeros((sz, sz))
A12 = np.eye(sz)
A21 = np.diag(OM2)
A22 = np.diag(ZI)

A = np.vstack((np.concatenate((A11, A12), axis=1), np.concatenate((A21, A22), axis=1)))

B1 = np.zeros((sz, sz))    
B2 = [[1e-6, 1e-7],[2e-6, 2e-7]]
B = np.vstack((B1, B2))

C1 = np.zeros((sz, sz*2))
C1[0] = [1e-4, 2*1e-4, 3*1e-4, 5*1e-5]
C1[1] = [2e-4, 3.5*1e-4, 1.5*1e-4, 2*1e-5]
C2 = np.zeros((sz*2, sz))
C = np.concatenate((C1.T, C2), axis=1)

D = np.zeros((sz*2, sz))

sys = ss(A, B, C, D)

tEnd = 1
time = np.arange(0, tEnd, 1e-3)
tout, youtStep = step_response(sys, T=time)
tout, youtForced, xout = forced_response(sys, T=time, U=1.0)
plt.figure()
for k, y in enumerate(youtStep):
    plt.subplot(4,1,k+1)
    plt.grid(True)
    plt.plot(tout, y,label='step')
    plt.plot(tout, youtForced[k], '--r',label='forced')
    if k == 0:
        plt.legend()
plt.xlabel('Time [s]')

Solution

OK the step response is easily manageable via the function control.matlab.step which actually allows the selection of the different input of the MIMO system, something I initially ignored, but was well reported in the official documentation:

https://python-control.readthedocs.io/en/0.8.1/generated/control.matlab.step.html

Here's the output [MIMO step response output]

Luckily it was an easy fix :)

from control import ss
import control.matlab as ctl
import numpy as np
import matplotlib.pyplot as plt

sz = 2

f1 = 1*2*np.pi
f2 = 1.5*2*np.pi
OM2 = [-f1**2, -f2**2]
ZI = [-2*f1*0.01, -2*f2*0.01]

A11 = np.zeros((sz, sz))
A12 = np.eye(sz)
A21 = np.diag(OM2)
A22 = np.diag(ZI)

A = np.vstack((np.concatenate((A11, A12), axis=1), np.concatenate((A21, A22), axis=1)))

B1 = np.zeros((sz, sz))    
B2 = [[1e-6, 1e-7],[2e-6, 2e-7]]
B = np.vstack((B1, B2))

C1 = np.zeros((sz, sz*2))
C1[0] = [1e-4, 2*1e-4, 3*1e-4, 5*1e-5]
C1[1] = [2e-4, 3.5*1e-4, 1.5*1e-4, 2*1e-5]
C2 = np.zeros((sz*2, sz))
C = np.concatenate((C1.T, C2), axis=1)

D = np.zeros((sz*2, sz))

sys = ss(A, B, C, D)

tEnd = 100
time = np.arange(0, tEnd, 1e-3)
yy1, tt1 = ctl.step(sys, T=time, input=0)
yy2, tt2 = ctl.step(sys, T=time, input=1)
plt.figure()
for k in range(0, len(yy1[1,:])):
    plt.subplot(4,1,k+1)
    plt.grid(True)
    plt.plot(tt1, yy1[:,k], label='input=0')
    plt.plot(tt2, yy2[:,k], label='input=1')
    if k == 0:
        plt.legend()
plt.xlabel('Time [s]')