Fundamentals of Linear Control
Mauricio de Oliveira
Supplemental material for Chapter 8
Before You Start
In this script you will perform calculations with transfer-functions. In addition to the commands already used, the following MATLAB commands will be used:
- sim, to simulate a nonlinear model using simulink.
You will also use the following auxiliary commands:
- rectangle, to plot circles.
8.1 Closed-Loop Stability and Performance
Create tf objects to represent the pendulum linearized about the stable and the unstable equlibria:
% Pendulum parameters
m = 0.5
l = 0.3
r = l/2
b = 0
g = 9.8
J = m*l^2/12
Jr = (J+m*r^2)
% models linearized around equilibria
Gpi = tf(1/Jr, [1, b/Jr, -m*r*g/Jr])
G0 = tf(1/Jr, [1, b/Jr, m*r*g/Jr])
wn = 7;
Create tf objects to represent the controllers 
and 
from Sections 6.5 and 7.8:
% C6 controller
C6 = zpk([-6.855, -0.9573], [0, -21], 3.8321)
% C7 controlller
z1 = 1;
z2 = 5;
pc = 11;
C7 = tf(conv([1 z1],[1 z2]),conv([1 0],[1 pc]));
K = 2/abs(freqresp(C7,wn));
C7 = zpk(K * C7)
Calculate the closed-loop sensitivity transfer-functions with controller :
S6 = feedback(1,Gpi*C6)
damp(S6)
and with controller :
S7 = feedback(1,Gpi*C7)
damp(S7)
Calculate the loop transfer-function with :
L6 = Gpi*C6
and with :
L7 = Gpi*C7
Plot the combined frequency response:
xl = [-1, 2];
yl1 = [-30, 20];
yl2 = [0, 270];
[magS7,phsS7,w] = bode(S7, {10^xl(1), 10^xl(2)});
[magS6,phsS6] = bode(S6, w);
[magL6,phsL6] = bode(L6, w);
[magL7,phsL7] = bode(L7, w);
% Fix phases
phsL6 = phsL6 + 360;
phsL7 = phsL7 + 360;
xl = [-1 2];
yl1 = [-30 20];
yl2 = [0 270];
% Fig. 8.2 (a):
figure()
hline = semilogx(w, 20*log10([magS6(:)'; magS7(:)'; magL6(:)'; magL7(:)']));
grid on
xlim(10.^xl);
ylim(yl1);
set(gca, 'YTick', [-30 -20 -10 0 10 20]);
l = legend('$S_6$', '$S_7$', '$L_6$', '$L_7$', 'Location', 'NorthWest');
set(l, 'interpreter', 'latex');
ylabel('Magnitude (dB)');
xlabel('\omega (rad/s)');
% Fig. 8.2 (b):
figure()
hline = semilogx(w, ([phsS6(:)'; phsS7(:)'; phsL6(:)'; phsL7(:)']));
grid on
xlim(10.^xl);
ylim(yl2);
set(gca, 'YTick', [0 90 180 270 360]);
l = legend('$S_6$', '$S_7$', '$L_6$', '$L_7$', 'Location', 'SouthEast');
set(l, 'interpreter', 'latex');
ylabel('Phase (deg)');
xlabel('\omega (rad/s)');
8.4 Control of the Simple Pendulum - Part III
Create tf objects to represent the controller 
:
% C8 controller
z1 = 1;
z2 = 3;
pc = 10.5;
C8 = tf(conv([1 z1],[1 z2]),conv([1 0],[1 pc]))
K = 2/abs(freqresp(C8,wn));
C8 = zpk(K * C8)
Compare frequency response of the controllers:
% Bode plots
xl = [-1 2];
[magC6,phsC6,w] = bode(C6, {10^xl(1), 10^xl(2)});
[magC7,phsC7] = bode(C7, w);
[magC8,phsC8] = bode(C8, w);
% Fig. 8.8:
figure()
yl1 = [0 20];
hline = semilogx(w, 20*log10([magC6(:)'; magC7(:)'; magC8(:)']));
grid on
xlim(10.^xl);
ylim(yl1);
l = legend('$C_6$', '$C_7$', '$C_8$', ...
'Location', 'SouthWest');
set(l, 'interpreter', 'latex');
ylabel('Magnitude (dB)');
xlabel('\omega (rad/s)');
Calculate the loop transfer-function with :
L8 = Gpi*C8
Plot the combined frequency response:
[magL6,phsL6,w] = bode(L6, {10^xl(1), 10^xl(2)});
[magL7,phsL7] = bode(L7, w);
[magL8,phsL8] = bode(L8, w);
% Fix phases
phsL6 = phsL6 + 360;
phsL7 = phsL7 + 360;
phsL8 = phsL8 + 360;
% Fig. 8.9a:
figure()
yl2 = [-30 20];
hline = semilogx(w, 20*log10([magL6(:)'; magL7(:)'; magL8(:)']));
grid on
xlim(10.^xl);
ylim(yl2);
l = legend('$L_6$', '$L_7$', '$L_8$', ...
'Location', 'SouthWest');
set(l, 'interpreter', 'latex');
ylabel('Magnitude (dB)');
xlabel('\omega (rad/s)');
% Fig. 8.9b:
figure()
yl2 = [90 225];
hline = semilogx(w, ([phsL6(:)'; phsL7(:)'; phsL8(:)']));
grid on
xlim(10.^xl);
ylim(yl2);
set(gca, 'YTick', [0 90 180 270 360]);
l = legend('$L_6$', '$L_7$', '$L_8$', ...
'Location', 'NorthWest');
set(l, 'interpreter', 'latex');
ylabel('Phase (deg)');
xlabel('\omega (rad/s)');
Create auxiliary model for robust analysis:
a1 = b/Jr
a2 = m*g*r/Jr
b2 = 1/Jr
F = tf(1, [1 a1 0])
with controller :
G6 = -a2*feedback(F,b2*C6)
with controller :
G7 = -a2*feedback(F,b2*C7)
with controller :
G8 = -a2*feedback(F,b2*C8)
and evaluate the norm 
for each transfer-function:
norm(G6,inf)
norm(G7,inf)
norm(G8,inf)
which confirm that the small gain theorem can prove that the controllers and can stabilize the nonlinear model of the simple pendulum, but not .
Repeat the robustness study this time using the circle criterion by verifying if the Nyquist diagram of the transfer-functions 
, 
and 
are in the circle 
:
% Nyquist diagrams
eps0 = 0.02;
offset = 0.1;
[re,im] = nyquist(G6, {eps0,1e2});
b6 = re(:) + j*im(:);
[re,im] = nyquist(G7, {eps0,1e2});
b7 = re(:) + j*im(:);
eps0 = 0.1;
offset = 0.1;
[re,im] = nyquist(G8, {eps0,1e2});
b8 = re(:) + j*im(:);
% Unit circles C(-1,1)
t = 0 : 0.1 : 2*pi;
circ1 = cos(t) + j*sin(t);
% Circle C(-0.22,1)
alpha = -0.22;
beta = 1;
cc = -(1/alpha+1/beta)/2;
rr = abs(-1/alpha+1/beta)/2;
circ2 = cc + rr*(cos(t) + j*sin(t));
% Fig. 8.15:
figure()
xl = [-1.25 4.75];
yl = 1.2*[-1.25 1.25];
plot(real([b6; b6]), imag([b6; -b6]), ...
real([b7; b7]), imag([b7; -b7]), ...
real([b8; b8]), imag([b8; -b8]))
% The following commands plot the circles
% in spite of the name rectangle :)
h = rectangle('Position',[-1,-1,2,2], 'Curvature',[1,1], ...
'FaceColor', 0.1*[1 0 0] + .9*[1 1 1]);
uistack(h, 'bottom')
h = rectangle('Position',[-1,-rr,2*rr,2*rr], 'Curvature',[1,1], ...
'FaceColor', 0.1*[0 1 0] + .9*[1 1 1]);
uistack(h, 'bottom')
l = legend('$G_6$', '$G_7$', '$G_8$', 'Location', 'SouthEast');
set(l, 'interpreter', 'latex');
xlim(xl)
ylim(yl)
xlabel('Re');
ylabel('Im');
grid on
set(gca, 'layer', 'top')
pbaspect([1 diff(ylim)/diff(xlim) 1]);
Alternatively calculate the parameters (8.16):
aa = (alpha+beta)/2
bb = (beta-alpha)/2
and verify the algebraic condition (8.18):
norm((bb/aa)*feedback(aa*G6,1),inf)
norm((bb/aa)*feedback(aa*G7,1),inf)
norm((bb/aa)*feedback(aa*G8,1),inf)
which should all be less than one.
Be patient, this step may take a while!
Repeat this study varying the damping parameter 
:
% Damping study
bs = linspace(0, 0.5, 100);
Gs = zeros(3,length(bs));
i = 1;
for bb = bs
% model for uncertainty feedback
aa1 = bb/Jr;
FF = tf(1, [1 aa1 0]);
GG6 = -a2*feedback(FF,b2*C6);
GG7 = -a2*feedback(FF,b2*C7);
GG8 = -a2*feedback(FF,b2*C8);
Gs(:,i) = [norm(GG6,inf); norm(GG7,inf); norm(GG8,inf)];
i = i + 1;
end
xl = [bs(1) bs(end)];
yl = [.25 1.25];
% Fig. 8.11:
figure()
plot(bs, Gs, ...
xl, [1 1], 'k--');
grid on
xlim(xl);
ylim(yl);
set(gca, 'YTick', [.25 0.5 0.75 1 1.25]);
set(gca, 'XTick', [0 0.1 0.2 0.3 0.4 0.5]);
l = legend('$G_6$', '$G_7$', '$G_8$', 'Location', 'NorthEast');
set(l, 'interpreter', 'latex');
ylabel('||G||_\infty');
xlabel('b');
Select a controller for simulation:
ctr = tf(C7) % make sure it is a tf object
and simulate using the simulink model
SimplePendulumWithFeedback
This model takes some variables from the workspace:
thetaBar = pi/2 - pi/4;
theta0 = pi/2;
num = ctr.num{1};
den = ctr.den{1};
You can interact with the simulation on the GUI interface or using the command line interface. The following will simulate the closed-loop system for 10 seconds:
T = 10
sim('SimplePendulumWithFeedback', T)
The following will recreate Fig. 8.12:
T1 = 4
thetaBar1 = pi/2 - pi/4;
theta0 = pi/2 + pi/4;
bs = [0 .1 .5];
N = length(bs);
y1 = cell(1,N);
t1 = cell(1,N);
i = 1;
for b = bs
thetaBar = thetaBar1;
[ti,x,yi] = sim('SimplePendulumWithFeedback', T1, simset('OutputVariables', 'ty', 'RelTol', 1e-5));
t1{i} = ti;
y1{i} = yi;
i = i + 1;
end
leg = [ones(N, 1) * '$b = ' num2str(bs','%2.1f') ones(N, 1) * '$'];
thetaBar2 = pi/2 + pi/4;
T2 = T1;
N = length(bs);
y2 = cell(1,N);
t2 = cell(1,N);
i = 1;
for b = bs
thetaBar = thetaBar2;
theta0 = y1{i}(end);
[ti,x,yi] = sim('SimplePendulumWithFeedback', T2, ...
simset('OutputVariables', 'ty', 'RelTol', 1e-5));
t2{i} = ti;
y2{i} = yi;
i = i + 1;
end
Y1 = [t1; cellfun(@(x) (180/pi)*x, y1, 'UniformOutput', false)];
Y2 = [cellfun(@(x) x+T1, t2, 'UniformOutput', false); cellfun(@(x) (180/pi)*x, y2, 'UniformOutput', false)];
yl = [-45 180];
plot(Y1{:}, ...
Y2{:}, ...
[0 T1+T2], (180/pi)*thetaBar1*[1 1], 'k--', ...
[0 T1+T2], (180/pi)*thetaBar2*[1 1], 'k--', ...
[0 T1+T2], (180/pi)*0*[1 1], 'k-')
l = legend(leg, 'Location', 'SouthEast');
set(l, 'interpreter', 'latex');
ylim(yl)
xlabel('t (s)')
ylabel('\theta (deg)')
grid on
set(gca, 'YTick', [-45 0 45 90 135 180]);
b = 0 % set b back 0