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.