-
Notifications
You must be signed in to change notification settings - Fork 411
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
SISO and MIMO mixed-sensitivity synthesis examples
Both examples taken from Skogestad and Postlethwaite's Multivariable Feedback Control.
- Loading branch information
Showing
2 changed files
with
282 additions
and
0 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,180 @@ | ||
"""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 numpy as np | ||
import matplotlib.pyplot as plt | ||
|
||
from control import tf, ss, mixsyn, feedback, 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 to 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 to 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 to 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.gamma | ||
|
||
|
||
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]') | ||
|
||
|
||
analysis() | ||
design() | ||
plt.show() |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,102 @@ | ||
"""Demonstrate mixed-sensitivity H-infinity design for a SISO plant. | ||
Based on Example 2.11 from Multivariable Feedback Control, Skogestad | ||
and Postlethwaite, 1st Edition. | ||
""" | ||
|
||
import numpy as np | ||
import matplotlib.pyplot as plt | ||
|
||
from control import tf, ss, mixsyn, feedback, step_response | ||
|
||
s = tf([1, 0], 1) | ||
# the plant | ||
g = 200/(10*s+1)/(0.05*s+1)**2 | ||
# disturbance plant | ||
gd = 100/(10*s+1) | ||
|
||
# first design | ||
# sensitivity weighting | ||
M = 1.5 | ||
wb = 10 | ||
A = 1e-4 | ||
ws1 = (s/M+wb)/(s+wb*A) | ||
# KS weighting | ||
wu = tf(1, 1) | ||
|
||
k1, cl1, info1 = mixsyn(g, ws1, wu) | ||
|
||
# sensitivity (S) and complementary sensitivity (T) functions for | ||
# design 1 | ||
s1 = feedback(1,g*k1) | ||
t1 = feedback(g*k1,1) | ||
|
||
# second design | ||
# this weighting differs from the text, where A**0.5 is used; if you use that, | ||
# the frequency response doesn't match the figure. The time responses | ||
# are similar, though. | ||
ws2 = (s/M**0.5+wb)**2/(s+wb*A)**2 | ||
# the KS weighting is the same as for the first design | ||
|
||
k2, cl2, info2 = mixsyn(g, ws2, wu) | ||
|
||
# S and T for design 2 | ||
s2 = feedback(1,g*k2) | ||
t2 = feedback(g*k2,1) | ||
|
||
# frequency response | ||
omega = np.logspace(-2,2,101) | ||
ws1mag,_,_ = ws1.freqresp(omega) | ||
s1mag,_,_ = s1.freqresp(omega) | ||
ws2mag,_,_ = ws2.freqresp(omega) | ||
s2mag,_,_ = s2.freqresp(omega) | ||
|
||
plt.figure(1) | ||
# text uses log-scaled absolute, but dB are probably more familiar to most control engineers | ||
plt.semilogx(omega,20*np.log10(s1mag.flat),label='$S_1$') | ||
plt.semilogx(omega,20*np.log10(s2mag.flat),label='$S_2$') | ||
# -1 in logspace is inverse | ||
plt.semilogx(omega,-20*np.log10(ws1mag.flat),label='$1/w_{P1}$') | ||
plt.semilogx(omega,-20*np.log10(ws2mag.flat),label='$1/w_{P2}$') | ||
|
||
plt.ylim([-80,10]) | ||
plt.xlim([1e-2,1e2]) | ||
plt.xlabel('freq [rad/s]') | ||
plt.ylabel('mag [dB]') | ||
plt.legend() | ||
plt.title('Sensitivity and sensitivity weighting frequency responses') | ||
|
||
# time response | ||
time = np.linspace(0,3,201) | ||
_,y1 = step_response(t1, time) | ||
_,y2 = step_response(t2, time) | ||
|
||
# gd injects into the output (that is, g and gd are summed), and the | ||
# closed loop mapping from output disturbance->output is S. | ||
_,y1d = step_response(s1*gd, time) | ||
_,y2d = step_response(s2*gd, time) | ||
|
||
plt.figure(2) | ||
plt.subplot(1,2,1) | ||
plt.plot(time,y1,label='$y_1(t)$') | ||
plt.plot(time,y2,label='$y_2(t)$') | ||
|
||
plt.ylim([-0.1,1.5]) | ||
plt.xlim([0,3]) | ||
plt.xlabel('time [s]') | ||
plt.ylabel('signal [1]') | ||
plt.legend() | ||
plt.title('Tracking response') | ||
|
||
plt.subplot(1,2,2) | ||
plt.plot(time,y1d,label='$y_1(t)$') | ||
plt.plot(time,y2d,label='$y_2(t)$') | ||
|
||
plt.ylim([-0.1,1.5]) | ||
plt.xlim([0,3]) | ||
plt.xlabel('time [s]') | ||
plt.ylabel('signal [1]') | ||
plt.legend() | ||
plt.title('Disturbance response') | ||
|
||
plt.show() |