ERA example, mass spring damper system

Code

 1# era_msd.py
 2# Johannes Kaisinger, 4 July 2024
 3#
 4# Demonstrate estimation of State Space model from impulse response.
 5# SISO, SIMO, MISO, MIMO case
 6
 7import numpy as np
 8import matplotlib.pyplot as plt
 9import os
10
11import control as ct
12
13# set up a mass spring damper system (2dof, MIMO case)
14# Mechanical Vibrations: Theory and Application, SI Edition, 1st ed.
15# Figure 6.5 / Example 6.7
16# m q_dd + c q_d + k q = f
17m1, k1, c1 = 1., 4., 1.
18m2, k2, c2 = 2., 2., 1.
19k3, c3 = 6., 2.
20
21A = np.array([
22    [0., 0., 1., 0.],
23    [0., 0., 0., 1.],
24    [-(k1+k2)/m1, (k2)/m1, -(c1+c2)/m1, c2/m1],
25    [(k2)/m2, -(k2+k3)/m2, c2/m2, -(c2+c3)/m2]
26])
27B = np.array([[0.,0.],[0.,0.],[1/m1,0.],[0.,1/m2]])
28C = np.array([[1.0, 0.0, 0.0, 0.0],[0.0, 1.0, 0.0, 0.0]])
29D = np.zeros((2,2))
30
31xixo_list = ["SISO","SIMO","MISO","MIMO"]
32xixo = xixo_list[3] # choose a system for estimation
33match xixo:
34    case "SISO":
35        sys = ct.StateSpace(A, B[:,0], C[0,:], D[0,0])
36    case "SIMO":
37        sys = ct.StateSpace(A, B[:,:1], C, D[:,:1])
38    case "MISO":
39        sys = ct.StateSpace(A, B, C[:1,:], D[:1,:])
40    case "MIMO":
41        sys = ct.StateSpace(A, B, C, D)
42
43
44dt = 0.1
45sysd = sys.sample(dt, method='zoh')
46response = ct.impulse_response(sysd)
47response.plot()
48plt.show()
49
50sysd_est, _ = ct.eigensys_realization(response,r=4,dt=dt)
51
52step_true = ct.step_response(sysd)
53step_true.sysname="H_true"
54step_est = ct.step_response(sysd_est)
55step_est.sysname="H_est"
56
57step_true.plot(title=xixo)
58step_est.plot(color='orange',linestyle='dashed')
59
60plt.show()
61
62if 'PYCONTROL_TEST_EXAMPLES' not in os.environ:
63    plt.show()

Notes

1. The environment variable PYCONTROL_TEST_EXAMPLES is used for testing to turn off plotting of the outputs.0