MIMO robust control example (SP96, Example 3.8)

Code

  1"""robust_mimo.py
  2
  3Demonstrate mixed-sensitivity H-infinity design for a MIMO plant.
  4
  5Based on Example 3.8 from Multivariable Feedback Control, Skogestad and Postlethwaite, 1st Edition.
  6"""
  7
  8import os
  9
 10import numpy as np
 11import matplotlib.pyplot as plt
 12
 13from control import tf, ss, mixsyn, step_response
 14
 15
 16def weighting(wb, m, a):
 17    """weighting(wb,m,a) -> wf
 18    wb - design frequency (where |wf| is approximately 1)
 19    m - high frequency gain of 1/wf; should be > 1
 20    a - low frequency gain of 1/wf; should be < 1
 21    wf - SISO LTI object
 22    """
 23    s = tf([1, 0], [1])
 24    return (s/m + wb) / (s + wb*a)
 25
 26
 27def plant():
 28    """plant() -> g
 29    g - LTI object; 2x2 plant with a RHP zero, at s=0.5.
 30    """
 31    den = [0.2, 1.2, 1]
 32    gtf = tf([[[1], [1]],
 33              [[2, 1], [2]]],
 34             [[den, den],
 35              [den, den]])
 36    return ss(gtf)
 37
 38
 39# as of this writing (2017-07-01), python-control doesn't have an
 40# equivalent to Matlab's sigma function, so use a trivial stand-in.
 41def triv_sigma(g, w):
 42    """triv_sigma(g,w) -> s
 43    g - LTI object, order n
 44    w - frequencies, length m
 45    s - (m,n) array of singular values of g(1j*w)"""
 46    m, p, _ = g.frequency_response(w)
 47    sjw = (m*np.exp(1j*p)).transpose(2, 0, 1)
 48    sv = np.linalg.svd(sjw, compute_uv=False)
 49    return sv
 50
 51
 52def analysis():
 53    """Plot open-loop responses for various inputs"""
 54    g = plant()
 55
 56    t = np.linspace(0, 10, 101)
 57    _, yu1 = step_response(g, t, input=0, squeeze=True)
 58    _, yu2 = step_response(g, t, input=1, squeeze=True)
 59
 60    # linear system, so scale and sum previous results to get the
 61    # [1,-1] response
 62    yuz = yu1 - yu2
 63
 64    plt.figure(1)
 65    plt.subplot(1, 3, 1)
 66    plt.plot(t, yu1[0], label='$y_1$')
 67    plt.plot(t, yu1[1], label='$y_2$')
 68    plt.xlabel('time')
 69    plt.ylabel('output')
 70    plt.ylim([-1.1, 2.1])
 71    plt.legend()
 72    plt.title('o/l response\nto input [1,0]')
 73
 74    plt.subplot(1, 3, 2)
 75    plt.plot(t, yu2[0], label='$y_1$')
 76    plt.plot(t, yu2[1], label='$y_2$')
 77    plt.xlabel('time')
 78    plt.ylabel('output')
 79    plt.ylim([-1.1, 2.1])
 80    plt.legend()
 81    plt.title('o/l response\nto input [0,1]')
 82
 83    plt.subplot(1, 3, 3)
 84    plt.plot(t, yuz[0], label='$y_1$')
 85    plt.plot(t, yuz[1], label='$y_2$')
 86    plt.xlabel('time')
 87    plt.ylabel('output')
 88    plt.ylim([-1.1, 2.1])
 89    plt.legend()
 90    plt.title('o/l response\nto input [1,-1]')
 91
 92
 93def synth(wb1, wb2):
 94    """synth(wb1,wb2) -> k,gamma
 95    wb1: S weighting frequency
 96    wb2: KS weighting frequency
 97    k: controller
 98    gamma: H-infinity norm of 'design', that is, of evaluation system
 99    with loop closed through design
100    """
101    g = plant()
102    wu = ss([], [], [], np.eye(2))
103    wp1 = ss(weighting(wb=wb1, m=1.5, a=1e-4))
104    wp2 = ss(weighting(wb=wb2, m=1.5, a=1e-4))
105    wp = wp1.append(wp2)
106    k, _, info = mixsyn(g, wp, wu)
107    return k, info[0]
108
109
110def step_opposite(g, t):
111    """reponse to step of [-1,1]"""
112    _, yu1 = step_response(g, t, input=0, squeeze=True)
113    _, yu2 = step_response(g, t, input=1, squeeze=True)
114    return yu1 - yu2
115
116
117def design():
118    """Show results of designs"""
119    # equal weighting on each output
120    k1, gam1 = synth(0.25, 0.25)
121    # increase "bandwidth" of output 2 by moving crossover weighting frequency 100 times higher
122    k2, gam2 = synth(0.25, 25)
123    # now weight output 1 more heavily
124    # won't plot this one, just want gamma
125    _, gam3 = synth(25, 0.25)
126
127    print('design 1 gamma {:.3g} (Skogestad: 2.80)'.format(gam1))
128    print('design 2 gamma {:.3g} (Skogestad: 2.92)'.format(gam2))
129    print('design 3 gamma {:.3g} (Skogestad: 6.73)'.format(gam3))
130
131    # do the designs
132    g = plant()
133    w = np.logspace(-2, 2, 101)
134    I = ss([], [], [], np.eye(2))
135    s1 = I.feedback(g*k1)
136    s2 = I.feedback(g*k2)
137
138    # frequency response
139    sv1 = triv_sigma(s1, w)
140    sv2 = triv_sigma(s2, w)
141
142    plt.figure(2)
143
144    plt.subplot(1, 2, 1)
145    plt.semilogx(w, 20*np.log10(sv1[:, 0]), label=r'$\sigma_1(S_1)$')
146    plt.semilogx(w, 20*np.log10(sv1[:, 1]), label=r'$\sigma_2(S_1)$')
147    plt.semilogx(w, 20*np.log10(sv2[:, 0]), label=r'$\sigma_1(S_2)$')
148    plt.semilogx(w, 20*np.log10(sv2[:, 1]), label=r'$\sigma_2(S_2)$')
149    plt.ylim([-60, 10])
150    plt.ylabel('magnitude [dB]')
151    plt.xlim([1e-2, 1e2])
152    plt.xlabel('freq [rad/s]')
153    plt.legend()
154    plt.title('Singular values of S')
155
156    # time response
157
158    # in design 1, both outputs have an inverse initial response; in
159    # design 2, output 2 does not, and is very fast, while output 1
160    # has a larger initial inverse response than in design 1
161    time = np.linspace(0, 10, 301)
162    t1 = (g*k1).feedback(I)
163    t2 = (g*k2).feedback(I)
164
165    y1 = step_opposite(t1, time)
166    y2 = step_opposite(t2, time)
167
168    plt.subplot(1, 2, 2)
169    plt.plot(time, y1[0], label='des. 1 $y_1(t))$')
170    plt.plot(time, y1[1], label='des. 1 $y_2(t))$')
171    plt.plot(time, y2[0], label='des. 2 $y_1(t))$')
172    plt.plot(time, y2[1], label='des. 2 $y_2(t))$')
173    plt.xlabel('time [s]')
174    plt.ylabel('response [1]')
175    plt.legend()
176    plt.title('c/l response to reference [1,-1]')
177
178
179if __name__ == "__main__":
180    analysis()
181    design()
182    if 'PYCONTROL_TEST_EXAMPLES' not in os.environ:
183        plt.show()

Notes

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