Estimation of Makrov parameters

Code

  1# markov.py
  2# Johannes Kaisinger, 4 July 2024
  3#
  4# Demonstrate estimation of markov parameters.
  5# SISO, SIMO, MISO, MIMO case
  6
  7import numpy as np
  8import matplotlib.pyplot as plt
  9import os
 10
 11import control as ct
 12
 13def create_impulse_response(H, time, transpose, dt):
 14    """Helper function to use TimeResponseData type for plotting"""
 15    
 16    H = np.array(H, ndmin=3)
 17
 18    if transpose:
 19        H = np.transpose(H)
 20    
 21    q, p, m = H.shape
 22    inputs = np.zeros((p,p,m))
 23
 24    issiso = True if (q == 1 and p == 1) else False
 25    
 26    input_labels = []
 27    trace_labels, trace_types = [], []
 28    for i in range(p):
 29        inputs[i,i,0] = 1/dt # unit area impulse
 30        input_labels.append(f"u{[i]}")
 31        trace_labels.append(f"From u{[i]}")
 32        trace_types.append('impulse')
 33
 34    output_labels = []
 35    for i in range(q):
 36        output_labels.append(f"y{[i]}")
 37
 38    return ct.TimeResponseData(time=time[:m],
 39                            outputs=H,
 40                            output_labels=output_labels,
 41                            inputs=inputs,
 42                            input_labels=input_labels,
 43                            trace_labels=trace_labels,
 44                            trace_types=trace_types,
 45                            sysname="H_est",
 46                            transpose=transpose,
 47                            plot_inputs=False,
 48                            issiso=issiso)
 49
 50# set up a mass spring damper system (2dof, MIMO case)
 51# Mechanical Vibrations: Theory and Application, SI Edition, 1st ed.
 52# Figure 6.5 / Example 6.7
 53# m q_dd + c q_d + k q = f
 54m1, k1, c1 = 1., 4., 1.
 55m2, k2, c2 = 2., 2., 1.
 56k3, c3 = 6., 2.
 57
 58A = np.array([
 59    [0., 0., 1., 0.],
 60    [0., 0., 0., 1.],
 61    [-(k1+k2)/m1, (k2)/m1, -(c1+c2)/m1, c2/m1],
 62    [(k2)/m2, -(k2+k3)/m2, c2/m2, -(c2+c3)/m2]
 63])
 64B = np.array([[0.,0.],[0.,0.],[1/m1,0.],[0.,1/m2]])
 65C = np.array([[1.0, 0.0, 0.0, 0.0],[0.0, 1.0, 0.0, 0.0]])
 66D = np.zeros((2,2))
 67
 68
 69xixo_list = ["SISO","SIMO","MISO","MIMO"]
 70xixo = xixo_list[3] # choose a system for estimation
 71match xixo:
 72    case "SISO":
 73        sys = ct.StateSpace(A, B[:,0], C[0,:], D[0,0])
 74    case "SIMO":
 75        sys = ct.StateSpace(A, B[:,:1], C, D[:,:1])
 76    case "MISO":
 77        sys = ct.StateSpace(A, B, C[:1,:], D[:1,:])
 78    case "MIMO":
 79        sys = ct.StateSpace(A, B, C, D)
 80
 81dt = 0.25
 82sysd = sys.sample(dt, method='zoh')
 83sysd.name = "H_true"
 84
 85 # random forcing input
 86t = np.arange(0,100,dt)
 87u = np.random.randn(sysd.B.shape[-1], len(t))
 88
 89response = ct.forced_response(sysd, U=u)
 90response.plot()
 91plt.show()
 92
 93m = 50
 94ir_true = ct.impulse_response(sysd, T=dt*m)
 95
 96H_est = ct.markov(response, m, dt=dt)
 97# Helper function for plotting only
 98ir_est = create_impulse_response(H_est,
 99                                 ir_true.time,
100                                 ir_true.transpose,
101                                 dt)
102
103ir_true.plot(title=xixo)
104ir_est.plot(color='orange',linestyle='dashed')
105plt.show()
106
107if 'PYCONTROL_TEST_EXAMPLES' not in os.environ:
108    plt.show()

Notes

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