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.