MIMO robust control example (SP96, Example 3.8)

Code

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
"""robust_mimo.py

Demonstrate mixed-sensitivity H-infinity design for a MIMO plant.

Based on Example 3.8 from Multivariable Feedback Control, Skogestad and Postlethwaite, 1st Edition.
"""

import os

import numpy as np
import matplotlib.pyplot as plt

from control import tf, ss, mixsyn, step_response


def weighting(wb, m, a):
    """weighting(wb,m,a) -> wf
    wb - design frequency (where |wf| is approximately 1)
    m - high frequency gain of 1/wf; should be > 1
    a - low frequency gain of 1/wf; should be < 1
    wf - SISO LTI object
    """
    s = tf([1, 0], [1])
    return (s/m + wb) / (s + wb*a)


def plant():
    """plant() -> g
    g - LTI object; 2x2 plant with a RHP zero, at s=0.5.
    """
    den = [0.2, 1.2, 1]
    gtf = tf([[[1], [1]],
              [[2, 1], [2]]],
             [[den, den],
              [den, den]])
    return ss(gtf)


# as of this writing (2017-07-01), python-control doesn't have an
# equivalent to Matlab's sigma function, so use a trivial stand-in.
def triv_sigma(g, w):
    """triv_sigma(g,w) -> s
    g - LTI object, order n
    w - frequencies, length m
    s - (m,n) array of singular values of g(1j*w)"""
    m, p, _ = g.freqresp(w)
    sjw = (m*np.exp(1j*p*np.pi/180)).transpose(2, 0, 1)
    sv = np.linalg.svd(sjw, compute_uv=False)
    return sv


def analysis():
    """Plot open-loop responses for various inputs"""
    g = plant()

    t = np.linspace(0, 10, 101)
    _, yu1 = step_response(g, t, input=0)
    _, yu2 = step_response(g, t, input=1)

    yu1 = yu1
    yu2 = yu2

    # linear system, so scale and sum previous results to get the
    # [1,-1] response
    yuz = yu1 - yu2

    plt.figure(1)
    plt.subplot(1, 3, 1)
    plt.plot(t, yu1[0], label='$y_1$')
    plt.plot(t, yu1[1], label='$y_2$')
    plt.xlabel('time')
    plt.ylabel('output')
    plt.ylim([-1.1, 2.1])
    plt.legend()
    plt.title('o/l response\nto input [1,0]')

    plt.subplot(1, 3, 2)
    plt.plot(t, yu2[0], label='$y_1$')
    plt.plot(t, yu2[1], label='$y_2$')
    plt.xlabel('time')
    plt.ylabel('output')
    plt.ylim([-1.1, 2.1])
    plt.legend()
    plt.title('o/l response\nto input [0,1]')

    plt.subplot(1, 3, 3)
    plt.plot(t, yuz[0], label='$y_1$')
    plt.plot(t, yuz[1], label='$y_2$')
    plt.xlabel('time')
    plt.ylabel('output')
    plt.ylim([-1.1, 2.1])
    plt.legend()
    plt.title('o/l response\nto input [1,-1]')


def synth(wb1, wb2):
    """synth(wb1,wb2) -> k,gamma
    wb1: S weighting frequency
    wb2: KS weighting frequency
    k: controller
    gamma: H-infinity norm of 'design', that is, of evaluation system
    with loop closed through design
    """
    g = plant()
    wu = ss([], [], [], np.eye(2))
    wp1 = ss(weighting(wb=wb1, m=1.5, a=1e-4))
    wp2 = ss(weighting(wb=wb2, m=1.5, a=1e-4))
    wp = wp1.append(wp2)
    k, _, info = mixsyn(g, wp, wu)
    return k, info[0]


def step_opposite(g, t):
    """reponse to step of [-1,1]"""
    _, yu1 = step_response(g, t, input=0)
    _, yu2 = step_response(g, t, input=1)
    return yu1 - yu2


def design():
    """Show results of designs"""
    # equal weighting on each output
    k1, gam1 = synth(0.25, 0.25)
    # increase "bandwidth" of output 2 by moving crossover weighting frequency 100 times higher
    k2, gam2 = synth(0.25, 25)
    # now weight output 1 more heavily
    # won't plot this one, just want gamma
    _, gam3 = synth(25, 0.25)

    print('design 1 gamma {:.3g} (Skogestad: 2.80)'.format(gam1))
    print('design 2 gamma {:.3g} (Skogestad: 2.92)'.format(gam2))
    print('design 3 gamma {:.3g} (Skogestad: 6.73)'.format(gam3))

    # do the designs
    g = plant()
    w = np.logspace(-2, 2, 101)
    I = ss([], [], [], np.eye(2))
    s1 = I.feedback(g*k1)
    s2 = I.feedback(g*k2)

    # frequency response
    sv1 = triv_sigma(s1, w)
    sv2 = triv_sigma(s2, w)

    plt.figure(2)

    plt.subplot(1, 2, 1)
    plt.semilogx(w, 20*np.log10(sv1[:, 0]), label=r'$\sigma_1(S_1)$')
    plt.semilogx(w, 20*np.log10(sv1[:, 1]), label=r'$\sigma_2(S_1)$')
    plt.semilogx(w, 20*np.log10(sv2[:, 0]), label=r'$\sigma_1(S_2)$')
    plt.semilogx(w, 20*np.log10(sv2[:, 1]), label=r'$\sigma_2(S_2)$')
    plt.ylim([-60, 10])
    plt.ylabel('magnitude [dB]')
    plt.xlim([1e-2, 1e2])
    plt.xlabel('freq [rad/s]')
    plt.legend()
    plt.title('Singular values of S')

    # time response

    # in design 1, both outputs have an inverse initial response; in
    # design 2, output 2 does not, and is very fast, while output 1
    # has a larger initial inverse response than in design 1
    time = np.linspace(0, 10, 301)
    t1 = (g*k1).feedback(I)
    t2 = (g*k2).feedback(I)

    y1 = step_opposite(t1, time)
    y2 = step_opposite(t2, time)

    plt.subplot(1, 2, 2)
    plt.plot(time, y1[0], label='des. 1 $y_1(t))$')
    plt.plot(time, y1[1], label='des. 1 $y_2(t))$')
    plt.plot(time, y2[0], label='des. 2 $y_1(t))$')
    plt.plot(time, y2[1], label='des. 2 $y_2(t))$')
    plt.xlabel('time [s]')
    plt.ylabel('response [1]')
    plt.legend()
    plt.title('c/l response to reference [1,-1]')


if __name__ == "__main__":
    analysis()
    design()
    if 'PYCONTROL_TEST_EXAMPLES' not in os.environ:
        plt.show()

Notes

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