SISO robust control example (SP96, Example 2.1)
Code
1"""robust_siso.py
2
3Demonstrate mixed-sensitivity H-infinity design for a SISO plant.
4
5Based on Example 2.11 from Multivariable Feedback Control, Skogestad
6and Postlethwaite, 1st Edition.
7"""
8
9import os
10
11import numpy as np
12import matplotlib.pyplot as plt
13
14from control import tf, mixsyn, feedback, step_response
15
16s = tf([1, 0], 1)
17# the plant
18g = 200/(10*s + 1) / (0.05*s + 1)**2
19# disturbance plant
20gd = 100/(10*s + 1)
21
22# first design
23# sensitivity weighting
24M = 1.5
25wb = 10
26A = 1e-4
27ws1 = (s/M + wb) / (s + wb*A)
28# KS weighting
29wu = tf(1, 1)
30
31k1, cl1, info1 = mixsyn(g, ws1, wu)
32
33# sensitivity (S) and complementary sensitivity (T) functions for
34# design 1
35s1 = feedback(1, g*k1)
36t1 = feedback(g*k1, 1)
37
38# second design
39# this weighting differs from the text, where A**0.5 is used; if you use that,
40# the frequency response doesn't match the figure. The time responses
41# are similar, though.
42ws2 = (s/M ** 0.5 + wb)**2 / (s + wb*A)**2
43# the KS weighting is the same as for the first design
44
45k2, cl2, info2 = mixsyn(g, ws2, wu)
46
47# S and T for design 2
48s2 = feedback(1, g*k2)
49t2 = feedback(g*k2, 1)
50
51# frequency response
52omega = np.logspace(-2, 2, 101)
53ws1mag, _, _ = ws1.frequency_response(omega)
54s1mag, _, _ = s1.frequency_response(omega)
55ws2mag, _, _ = ws2.frequency_response(omega)
56s2mag, _, _ = s2.frequency_response(omega)
57
58plt.figure(1)
59# text uses log-scaled absolute, but dB are probably more familiar to most control engineers
60plt.semilogx(omega, 20*np.log10(s1mag.flat), label='$S_1$')
61plt.semilogx(omega, 20*np.log10(s2mag.flat), label='$S_2$')
62# -1 in logspace is inverse
63plt.semilogx(omega, -20*np.log10(ws1mag.flat), label='$1/w_{P1}$')
64plt.semilogx(omega, -20*np.log10(ws2mag.flat), label='$1/w_{P2}$')
65
66plt.ylim([-80, 10])
67plt.xlim([1e-2, 1e2])
68plt.xlabel('freq [rad/s]')
69plt.ylabel('mag [dB]')
70plt.legend()
71plt.title('Sensitivity and sensitivity weighting frequency responses')
72
73# time response
74time = np.linspace(0, 3, 201)
75_, y1 = step_response(t1, time)
76_, y2 = step_response(t2, time)
77
78# gd injects into the output (that is, g and gd are summed), and the
79# closed loop mapping from output disturbance->output is S.
80_, y1d = step_response(s1*gd, time)
81_, y2d = step_response(s2*gd, time)
82
83plt.figure(2)
84plt.subplot(1, 2, 1)
85plt.plot(time, y1, label='$y_1(t)$')
86plt.plot(time, y2, label='$y_2(t)$')
87
88plt.ylim([-0.1, 1.5])
89plt.xlim([0, 3])
90plt.xlabel('time [s]')
91plt.ylabel('signal [1]')
92plt.legend()
93plt.title('Tracking response')
94
95plt.subplot(1, 2, 2)
96plt.plot(time, y1d, label='$y_1(t)$')
97plt.plot(time, y2d, label='$y_2(t)$')
98
99plt.ylim([-0.1, 1.5])
100plt.xlim([0, 3])
101plt.xlabel('time [s]')
102plt.ylabel('signal [1]')
103plt.legend()
104plt.title('Disturbance response')
105
106if 'PYCONTROL_TEST_EXAMPLES' not in os.environ:
107 plt.show()
Notes
1. The environment variable PYCONTROL_TEST_EXAMPLES is used for testing to turn off plotting of the outputs.