# 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.