MIT OpenCourseWare

http://ocw.mit.edu

16.323 Principles of Optimal Control

Spring 2008

For information about citing these materials or our Terms of Use, visit: http://ocw.mit.edu/terms.
 16.323 Lecture 6

Calculus of Variations applied to Optimal Control

ẋ = a(x, u, t)
ṗ = −HxT
Hu = 0
Spr 2008 16.323 6–1
Optimal Control Problems

• Are now ready to tackle the optimal control problem

– Start with simple terminal constraints
� tf
J = h(x(tf ), tf ) + g(x(t), u(t), t)dt
t0

with the system dynamics

ẋ(t) = a(x(t), u(t), t)
– t0, x(t0) fixed
– tf free
– x(tf ) are fixed or free by element
• Note that this looks a bit different because we have u(t) in the inte­
grand, but consider that with a simple substitution, we get
ẋ=a(x,u,t)
g̃(x, ẋ, t) → ĝ(x, u, t)
• Note that the differential equation of the dynamics acts as a constraint
that we must adjoin using a Lagrange multiplier, as before:
� tf
T
� �
Ja = h(x(tf ), tf )+ g(x(t), u(t), t) + p {a(x(t), u(t), t) − ẋ} dt
t0

• Find the variation:10

� tf
gxδx + guδu + (a − ẋ)T δp(t)
�
δJa = hxδxf + htf δtf +
t0
T T
� � �
+p (t){axδx + auδu − δẋ} dt + g + p (a − ẋ) (tf )δtf
• Clean this up by defining the Hamiltonian: (See 4–4)
H(x, u, p, t) = g(x(t), u(t), t) + pT (t)a(x(t), u(t), t)
10 Take partials wrt each of the variables that the integrand is a function of.

June 18, 2008

Spr 2008 16.323 6–2

• Then
� �
T
δJa = hxδxf + htf + g + p (a − ẋ) (tf )δtf
� tf
T T
� �
+ Hxδx + Huδu + (a − ẋ) δp(t) − p (t)δẋ dt
t0

11
• To proceed, note that by integrating by parts we get:
� tf � tf

− pT (t)δẋdt = − pT (t)dδx

t0 t0
� tf � �T
T
�tf dp(t)
= −p δx�t + δxdt
0
t0 dt
� tf
T
= −p (tf )δx(tf ) + ṗT (t)δxdt
t0
� tf
= −pT (tf ) (δxf − ẋ(tf )δtf ) + ṗT (t)δxdt
t0

• So now can rewrite the variation as:

� �
T
δJa = hxδxf + htf + g + p (a − ẋ) (tf )δtf
� tf � tf
Hxδx + Huδu + (a − ẋ)T δp(t) dt − pT (t)δẋdt
� �
+
t0 t0
� �
T T T
� �
= hx − p (tf ) δxf + htf + g + p (a − ẋ) + p ẋ (tf )δtf
� tf
T T
�� � �
+ Hx + ṗ δx + Huδu + (a − ẋ) δp(t) dt
t0

11
� �
udv ≡ uv − vdu

June 18, 2008

Spr 2008 16.323 6–3

• So necessary conditions for δJa = 0 are that for t ∈ [t0, tf ]

ẋ = a(x, u, t) (dim n)
ṗ = −HxT (dim n)
Hu = 0 (dim m)
– With the boundary condition (lost if tf is fixed) that

htf + g + pT a = htf + H(tf ) = 0

– Add the boundary constraints that x(t0) = x0 (dim n)

– If xi(tf ) is fixed, then xi(tf ) = xif

∂h

– If xi(tf ) is free, then pi(tf ) = (tf ) for a total (dim n)

∂xi

• These necessary conditions have 2n differential and m algebraic equa­

tions with 2n+1 unknowns (if tf free), found by imposing the (2n+1)
boundary conditions.

June 18, 2008

Spr 2008 16.323 6–4

• Note the symmetry in the differential equations:

� �T
∂H
ẋ = a(x, u, t) =
∂p
�T T
∂(g + pT a)
�
∂H
ṗ = − =−
∂x ∂x
� �T � �T
∂a ∂g
= − p−
∂x ∂x
– So the dynamics of p, called the costate, are linearized system
dynamics (negative transpose – dual)
⎡ ⎤
∂a1 ∂a1
∂x1 . . . ∂xn
� �
∂a ...
=⎣
⎢ ⎥
∂x
⎦
∂an ∂an
∂x1 . . . ∂xn

• These necessary conditions are extremely important, and we will be

using them for the rest of the term.

June 18, 2008

Spr 2008 16.323 6–5

Control with General Terminal Conditions

• Can develop similar conditions in the case of more general terminal
conditions with tf free and

m(x(tf ), tf ) = 0

• Follow the same procedure on 6–1 using the insights provided on 5–21
(using the ga form on 5–20) to form

w(x(tf ), ν, tf ) = h(x(tf ), tf ) + ν T m(x(tf ), tf )

• Work through the math, and get the necessary conditions are
ẋ = a(x, u, t) (dim n) (6.22)
ṗ = −HxT (dim n) (6.23)
Hu = 0 (dim m) (6.24)
– With the boundary condition (lost if tf fixed)
H(tf ) + wtf (tf ) = 0
– And m(x(tf ), tf ) = 0, with x(t0) and t0 given.
– With (since x(tf ) is not directly given)
� �T
∂w
p(tf ) = (tf )
∂x

• Collapses to form on 6–3 if m not present – i.e., does not constrain

x(tf )

June 18, 2008

Spr 2008 16.323 6–6
Example 6–1

• Simple double integrator system starting at y(0) = 10, ẏ(0) = 0,

must drive to origin y(tf ) = ẏ(tf ) = 0 to minimize the cost (b > 0)
1 2 1 tf 2
�
J = αtf + bu (t)dt
2 2 0

• Define the dynamics with x1 = y, x2 = ẏ so that

� � � �
0 1 0
ẋ(t) = Ax(t) + Bu(t) A = B=
0 0 1

• With p(t) = [p1(t) p2(t)]T , define the Hamiltonian

1
H = g + pT (t)a = bu2 + pT (t) (Ax(t) + Bu(t))
2

• The necessary conditions are then that:

∂H
ṗ1 = − ∂x 1
= 0 → p1(t) = c1
ṗ = −HxT , →
∂H
ṗ2 = − ∂x 2
= −p1 → p2(t) = −c1t + c2
p2 c2 c1
Hu = bu + p2 = 0 → u=− =− + t
b b b

• Now impose the boundary conditions:

1
H(tf ) + ht(tf ) = bu2(tf ) + p1(tf )x2(tf ) + p2(tf )u(tf ) + αtf = 0
2

= bu2(tf ) + (−bu(tf ))u(tf ) + αtf

2
1 1
= − bu2(tf ) + αtf = 0 → tf = (−c2 + c1tf )2
2 2bα

June 18, 2008

Spr 2008 16.323 6–7

• Now go back to the state equations:

c2 c1 c2 c1
ẋ2(t) = − + t → x2(t) = c3 − t + t2
b b b 2b

and since x2(0) = 0, c3 = 0, and

c2 c1
ẋ1(t) = x2(t) → x1(t) = c4 − t2 + t3
2b 6b
and since x1(0) = 10, c4 = 10

• Now note that

c2 c1
x2(tf ) = − tf + t2f = 0

b 2b
c2 2 c1 3
x1(tf ) = 10 − tf + tf = 0

2b 6b
c2 60b 120b
= 10 − t2f = 0 → c2 = 2 , c1 =
6b tf t3f
– But that gives us:
� �2
1 60b 120b (60b)2
tf = − 2 + 3 tf =
2bα tf tf 2bαt4f

so that t5f
= 1800b/α or tf ≈ 4.48(b/α)1/5, which makes sense
because tf goes down as α goes up.

– Finally, c2 = 2.99b3/5α2/5 and c1 = 1.33b2/5α3/5

Figure 6.1: Example 6–1

June 18, 2008

Spr 2008 16.323 6–8

Example 6–1

1 %

2 % Simple opt example showing impact of weight on t_f

3 % 16.323 Spring 2008

4 % Jonathan How

5 % opt1.m

6 %

7 clear all;close all;

8 set(0, ’DefaultAxesFontSize’, 14, ’DefaultAxesFontWeight’,’demi’)

9 set(0, ’DefaultTextFontSize’, 14, ’DefaultTextFontWeight’,’demi’)

10 %

11 A=[0 1;0 0];B=[0 1]’;C=eye(2);D=zeros(2,1);

12 G=ss(A,B,C,D);

13 X0=[10 0]’;

14 b=0.1;

15
16 alp=1;

17 tf=(1800*b/alp)^0.2;

18 c1=120*b/tf^3;

19 c2=60*b/tf^2;

20 time=[0:1e-2:tf];

21 u=(-c2+c1*time)/b;

22 [y1,t1]=lsim(G,u,time,X0);

23
24 figure(1);clg

25 plot(time,u,’k-’,’LineWidth’,2);hold on

26 alp=10;

27 tf=(1800*b/alp)^0.2;

28 c1=120*b/tf^3;

29 c2=60*b/tf^2;

30 time=[0:1e-2:tf];

31 u=(-c2+c1*time)/b;

32 [y2,t2]=lsim(G,u,time,X0);

33 plot(time,u,’b--’,’LineWidth’,2);

34
35 alp=0.10;

36 tf=(1800*b/alp)^0.2;

37 c1=120*b/tf^3;

38 c2=60*b/tf^2;

39 time=[0:1e-2:tf];

40 u=(-c2+c1*time)/b;

41 [y3,t3]=lsim(G,u,time,X0);

42 plot(time,u,’g-.’,’LineWidth’,2);hold off

43
44 legend(’\alpha=1’,’\alpha=10’,’\alpha=0.1’)

45 xlabel(’Time (sec)’)

46 ylabel(’u(t)’)

47 title([’b= ’,num2str(b)])

48
49 figure(2);clg

50 plot(t1,y1(:,1),’k-’,’LineWidth’,2);

51 hold on

52 plot(t2,y2(:,1),’b--’,’LineWidth’,2);

53 plot(t3,y3(:,1),’g-.’,’LineWidth’,2);

54 hold off

55 legend(’\alpha=1’,’\alpha=10’,’\alpha=0.1’)

56 xlabel(’Time (sec)’)

57 ylabel(’y(t)’)

58 title([’b= ’,num2str(b)])

59
60 print -dpng -r300 -f1 opt11.png

61 print -dpng -r300 -f2 opt12.png

June 18, 2008

Spr 2008 16.323 6–9
LQR Variational Solution

• Deterministic Linear Quadratic Regulator

Plant:
ẋ(t) = A(t)x(t) + Bu(t)u(t), x(t0) = x0

z(t) = Cz (t)x(t)

Cost:
� tf � T
z (t)Rzz(t)z(t) + u (t)Ruu(t)u(t) dt + x(tf )T Ptf x(tf )
T
�
2JLQR =
t0

– Where Ptf ≥ 0, Rzz(t) > 0 and Ruu(t) > 0

– Define Rxx = CzT RzzCz ≥ 0
– A(t) is a continuous function of time.
– Bu(t), Cz (t), Rzz(t), Ruu(t) are piecewise continuous functions of
time, and all are bounded.

• Problem Statement: Find input u(t) ∀t ∈ [t0, tf ] to min JLQR

– This is not necessarily specified to be a feedback controller.

• To optimize the cost, we follow the procedure of augmenting the con­

straints in the problem (the system dynamics) to the cost (integrand)
to form the Hamiltonian:
1� T
H = x (t)Rxxx(t) + u (t)Ruuu(t) + pT (t) (Ax(t) + Buu(t))
T
�
2
– p(t) ∈ Rn×1 is called the Adjoint variable or Costate
– It is the Lagrange multiplier in the problem.

June 18, 2008

Spr 2008 16.323 6–10

• The necessary conditions (see 6–3) for optimality are that:

T
1. ẋ(t) = ∂H = Ax(t) + B(t)u(t) with x(t0) = x0
∂p
T
2. ṗ(t) = − ∂H = −Rxxx(t) − AT p(t) with p(tf ) = Ptf x(tf )
∂x

3. ∂H = 0 ⇒ Ruuu + BuT p(t) = 0, so u� = −Ruu

−1 T
Bu p(t)

∂u
2
4. As before, we can check for a minimum by looking at H2 ≥ 0
∂
∂u
(need to check that Ruu ≥ 0)
• Note that p(t) plays the same role as Jx�(x(t), t)T in previous solutions
to the continuous LQR problem (see 4–8).
– Main difference is there is no need to guess a solution for J �(x(t), t)
• Now have:
−1 T
ẋ(t) = Ax(t) + Bu�(t) = Ax(t) − BuRuu Bu p(t)
which can be combined with equation for the adjoint variable
ṗ(t) = −Rxxx(t) − AT p(t) = −CzT RzzCz x(t) − AT p(t)
−1 T
� ��
A −BuRuu Bu
� � �
ẋ(t) x(t)
⇒ =
ṗ(t) −CzT RzzCz −AT p(t)
� �� �
H
where H is called the Hamiltonian Matrix.
– Matrix describes coupled closed loop dynamics for both x and p.
– Dynamics of x(t) and p(t) are coupled, but x(t) known initially
and p(t) known at terminal time, since p(tf ) = Ptf x(tf )
– Two point boundary value problem ⇒ typically hard to solve.

June 18, 2008

Spr 2008 16.323 6–11

• However, in this case, we can introduce a new matrix variable P (t)

and show that:
1. p(t) = P (t)x(t)

2. It is relatively easy to find P (t).

• How proceed?
1. For the 2n system
� � −1 T
��
A −B R B
� �
ẋ(t) u uu u x(t)
=
ṗ(t) −CzT RzzCz −AT p(t)
define a transition matrix
� �
F11(t1, t0) F12(t1, t0)
F (t1, t0) =
F21(t1, t0) F22(t1, t0)
and use this to relate x(t) to x(tf ) and p(tf )
� � ��
F11(t, tf ) F12(t, tf )
� �
x(t) x(tf )
=
p(t) F21(t, tf ) F22(t, tf ) p(tf )
so
x(t) = F11(t, tf )x(tf ) + F12(t, tf )p(tf )
� �
= F11(t, tf ) + F12(t, tf )Ptf x(tf )

2. Now find p(t) in terms of x(tf )

� �
p(t) = F21(t, tf ) + F22(t, tf )Ptf x(tf )

3. Eliminate x(tf ) to get:

� �� �−1
p(t) = F21(t, tf ) + F22(t, tf )Ptf F11(t, tf ) + F12(t, tf )Ptf x(t)
� P (t)x(t)
June 18, 2008
 Spr 2008 16.323 6–12

• Now have p(t) = P (t)x(t), must find the equation for P (t)
ṗ(t) = Ṗ (t)x(t) + P (t)ẋ(t)
⇒ − CzT RzzCz x(t) − AT p(t) =

−Ṗ (t)x(t) = CzT RzzCz x(t) + AT p(t) + P (t)ẋ(t)

−1 T
= CzT Rzz Cz x(t) + AT p(t) + P (t)(Ax(t) − BuRuu Bu p(t))

−1 T
= (CzT RzzCz + P (t)A)x(t) + (AT − P (t)BuRuu Bu )p(t)

−1 T
= (CzT RzzCz + P (t)A)x(t) + (AT − P (t)BuRuu Bu )P (t)x(t)

−1 T
� T T
�
= A P (t) + P (t)A + Cz RzzCz − P (t)BuRuu Bu P (t) x(t)

• This must be true for arbitrary x(t), so P (t) must satisfy

−1 T
−Ṗ (t) = AT P (t) + P (t)A + CzT RzzCz − P (t)BuRuu Bu P (t)
– Which, of course, is the matrix differential Riccati Equation.
– Optimal value of P (t) is found by solving this equation backwards
in time from tf with P (tf ) = Ptf

June 18, 2008

Spr 2008 16.323 6–13

• The control gains are then

−1 T −1 T
uopt = −Ruu Bu p(t) = −Ruu Bu P (t)x(t) = −K(t)x(t)

• Optimal control inputs can in fact be computed using linear

feedback on the full system state
– Find optimal steady state feedback gains u(t) = −Kx(t) using
K = lqr(A, B, CzT RzzCz , Ruu)

• Key point: This controller works equally well for MISO and MIMO
regulator designs.

June 18, 2008

Spr 2008 16.323 6–14
Alternate Derivation of DRE
• On 6-10 we showed that:
� �� �−1
P (t) = F21(t, tf ) + F22(t, tf )Ptf F11(t, tf ) + F12(t, tf )Ptf

• To find the Riccati equation, note that

d −1
M (t) = −M −1(t)Ṁ (t)M −1(t)
dt
which gives
� �� �−1
Ṗ (t) = Ḟ21(t, tf ) + Ḟ22(t, tf )Ptf F11(t, tf ) + F12(t, tf )Ptf
� �� �−1
− F21(t, tf ) + F22(t, tf )Ptf F11(t, tf ) + F12(t, tf )Ptf ·
� �� �−1
Ḟ11(t, tf ) + Ḟ12(t, tf )Ptf F11(t, tf ) + F12(t, tf )Ptf

12
• Since F is the transition matrix for the system (see 6–10), then
d
F (t, tf ) = HF (t, tf )
dt

F˙ 11 −1 T
� � � � � �
Ḟ12 A −BuRuu Bu F11 F12
(t, tf ) = (t, tf ) (t, tf )
Ḟ21 Ḟ22 −Rxx −AT F21 F22

12 Consider homogeneous system ẋ(t) = A(t)x(t) with initial condition x(t ) = x . The general solution to this differential
0 0
equation is given by x(t) = Φ(t, t0 )x(t0 ) where Φ(t1 , t1 ) = I. Can show the following properties of the state transition matrix Φ:
1. Φ(t2 , t0 ) = Φ(t2 , t1 )Φ(t1 , t0 ), regardless of the order of the ti

2. Φ(t, τ ) = Φ(τ, t)−1

d
3. dt
Φ(t, t0 ) = A(t)Φ(t, t0 )

June 18, 2008

Spr 2008 16.323 6–15

• Now substitute and re-arrange:

� �
Ṗ = [Ḟ21 + Ḟ22Ptf ] − P [Ḟ11 + Ḟ12Ptf ] [F11 + F12Ptf ]−1
F˙ 11 = −1 T
AF11 − BuRuu Bu F21

F˙ 12 = −1 T

AF12 − BuRuu Bu F22

Ḟ21 = −RxxF11 − AT F21

Ḟ22 = −RxxF12 − AT F22

�� �
T T
Ṗ = −RxxF11 − A F21 + (−RxxF12 − A F22)Ptf
� ��
−1 T −1 T
−P AF11 − BuRuu Bu F21 + (AF12 − BuRuu Bu F22)Ptf [F11 + F12Ptf ]−1

• There are four terms:

−Rxx(F11 + F12Ptf )[F11 + F12Ptf ]−1 = −Rxx

−AT (F21 + F22Ptf )[F11 + F12Ptf ]−1 = −AT P

−P A(F11 + F12Ptf )[F11 + F12Ptf ]−1 = −P A

−1 T
P BuRuu Bu (F21 + F22Ptf )[F11 + F12Ptf ]−1 = P BuRuu
−1 T
Bu P

• Which, as expected, gives that

−1 T
−Ṗ = AT P + P A + Rxx − P BuRuu Bu P

June 18, 2008

Spr 2008 16.323 6–16

CARE Solution Algorithm

• Recall from (6–10) that
� � −1 T
��
A −BuRuu Bu
� �
ẋ(t) x(t)
=
ṗ(t) −CzT RzzCz −AT p(t)

• Assuming that the eigenvalues of H are unique, the Hamiltonian can

be diagonalized into the form:
� � ��
−Λ 0
� �
ż1(t) z1(t)
=
ż2(t) 0 Λ z2(t)
where diagonal matrix Λ is comprised of RHP eigenvalues of H.

• A similarity transformation exists between the states z1, z2 and x, p:

� � � � � � � �
x(t) z1(t) z1(t) x(t)
=Ψ ⇔ = Ψ−1
p(t) z2(t) z2(t) p(t)
where
−1 −1
� � � �
Ψ11 Ψ12 (Ψ )11 (Ψ )12
Ψ= and Ψ−1 =
Ψ21 Ψ22 (Ψ−1)21 (Ψ−1)22
and the columns of Ψ are the eigenvectors of H.

• Solving for z2(t) gives

z2(t) = eΛtz2(0) = [(Ψ−1)21x(t) + (Ψ−1)22p(t)]
= [(Ψ−1)21 + (Ψ−1)22P (t)]x(t)
– For the cost to be finite, need limt→∞ x(t) = 0, so can show that
lim z2(t) = 0
t→∞
– But given that the Λ dynamics in the RHP, this can only be true
if z2(0) = 0, which means that z2(t) = 0 ∀t

June 18, 2008

Spr 2008 16.323 6–17

• With this fact, note that

x(t) = Ψ11z1(t)
p(t) = Ψ21z1(t)
which can be combined to give:
p(t) = Ψ21(Ψ11)−1x(t) ≡ Pssx(t)

• Summary of solution algorithm:

– Find the eigenvalues and eigenvectors of H
– Select the n eigenvectors associated with the n eigenvalues in the
LHP.
– Form Ψ11 and Ψ21.
– Compute the steady state solution of the Riccati equation using
Pss = Ψ21(Ψ11)−1

% alternative calc of Riccati solution

H=[A -B*inv(Ruu)*B’ ; -Rxx -A’];

[V,D]=eig(H); % check order of eigenvalues

Psi11=V(1:2,1:2);

Psi21=V(3:4,1:2);

Ptest=Psi21*inv(Psi11);

June 18, 2008

Spr 2008 16.323 6–18
Optimal Cost

• Showed in earlier derivations that the optimal cost-to-go from the

initial (or any state) is of the form
1
J = xT (t0)P (t0)x(t0)
2
– Relatively clean way to show it for this derivation as well.
• Start with the standard cost and add zero (Ax + Buu − ẋ = 0)
1 tf � T
�
T T
�
JLQR = x Rxxx + u Ruuu + p (Ax + Buu − ẋ) dt
2 t0
1
+ x(tf )T Ptf x(tf )
2
• Now use the results of the necessary conditions to get:
ṗ = −HxT ⇒ pT A = −ṗT − xT Rxx
Hu = 0 ⇒ pT Bu = −uT Ruu
with p(tf ) = Ptf x(tf )
• Substitute these terms to get
� tf
1 1
= x(tf )T Ptf x(tf ) −
� T T
�
JLQR ṗ x + p ẋ dt
2 2 t0

� tf � �
1 1 d T
= x(tf )T Ptf x(tf ) − (p x) dt
2 2 t0 dt

1 1 �
x(tf )T Ptf x(tf ) − pT (tf )x(tf ) − pT (t0)x(t0)
�
=
2 2

1 T 1 � T T
�
= x(tf ) Ptf x(tf ) − x (tf )Ptf x(tf ) − x (t0)P (t0)x(t0)
2 2
1
= xT (t0)P (t0)x(t0)
2

June 18, 2008

Spr 2008 16.323 6–19
Pole Locations
• The closed-loop dynamics couple x(t) and p(t) and are given by

� � −1 T
��
A −B R B
� �
ẋ(t) u uu u x(t)
=
ṗ(t) −CzT RzzCz −AT p(t)
with the appropriate boundary conditions.

• OK, so where are the closed-loop poles of the system?

– Answer: must be eigenvalues of Hamiltonian matrix for the system:
−1 T
� �
A −BuRuu Bu
H �
−CzT RzzCz −AT

so we must solve det(sI − H) = 0.

• Key point: For a SISO system, we can relate the closed-loop poles
to a Symmetric Root Locus (SRL) for the transfer function
N (s)
Gzu(s) = Cz (sI − A)−1Bu =
D(s)
– Poles and zeros of Gzu(s) play an integral role in determining SRL
– Note Gzu(s) is the transfer function from control inputs to perfor­
mance variable.

• In fact, the closed-loop poles are given by the LHP roots of

Rzz
Δ(s) = D(s)D(−s) + N (s)N (−s) = 0
Ruu
– D(s)D(−s) + RRuu zz
N (s)N (−s) is drawn using standard root locus
rules - but it is symmetric wrt to both the real and imaginary axes.
– For a stable system, we clearly just take the poles in the LHP.

June 18, 2008

 Spr 2008 16.323 6–20
Derivation of the SRL
• The closed-loop poles are given by the eigenvalues of
−1 T
� �
A −BuRuu Bu
H� T → det(sI − H) = 0
−Cz RzzCz −AT

• Note: if A is invertible:
� �
A B
det = det(A) det(D − CA−1B)
C D

⇒ det(sI − H) = det(sI − A) det (sI + AT ) − CzT Rzz Cz (sI − A)−1 Bu Ruu

−1 T
� �
Bu

= det(sI − A) det(sI + AT ) det I − Cz T Rzz Cz (sI − A)−1 Bu Ruu

−1 T
Bu (sI + AT )−1
� �

• Also: det(I + ABC) = det(I + CAB), and if D(s) = det(sI − A),

then D(−s) = det(−sI − AT ) = (−1)n det(sI + AT )

−1 T
det(sI−H) = (−1)n D(s)D(−s) det I + Ruu Bu (−sI − AT )−1 CzT Rzz Cz (sI − A)−1 Bu
� �

• If Gzu(s) = Cz (sI −A)−1Bu, then GTzu(−s) = Bu T

(−sI −AT )−1CzT ,
so for SISO systems
n −1 T
� �
det(sI − H) = (−1) D(s)D(−s) det I + Ruu Gzu(−s)RzzGzu(s)
� �
R zz
= (−1)nD(s)D(−s) I + Gzu(−s)Gzu(s)
Ruu
� �
R zz
= (−1)n D(s)D(−s) + N (s)N (−s) = 0
Ruu

June 18, 2008

Spr 2008 16.323 6–21

Example 6–2
• Simple example from 4–12: A scalar
� ∞system2 with ẋ = ax + bu with
cost (Rxx > 0 and Ruu > 0) J = 0 (Rzzx (t) + Ruuu2(t)) dt
2 2
• The steady-state
√ 2 P2 solves 2aP + Rzz − P b /Ruu = 0 which gives
a+ a +b Rzz /Ruu
that P = −1 b2
Ruu
>0
√2 2
−1 a+ a +b Rzz /Ruu
– So that u(t) = −Kx(t) where K = Ruu bP = b
– and the closed-loop dynamics are
� �
b �
ẋ = (a − bK)x = a − (a + a2 + b2Rzz/Ruu) x
b
�
= − a2 + b2Rzz/Ruu x = Acl x(t)
• In this case, Gzu(s) = b/(s−a) so that N (s) = b and D(s) = (s−a),
and the SRL is of the form:
Rzz 2
(s − a)(−s − a) + b =0
Ruu
Symmetric root locus

0.8

0.6

0.4

0.2
Imaginary Axis

−0.2

−0.4

−0.6

−0.8

−1
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2
Real Axis

• SRL is the same whether a < 0 (OL stable) or a > 0 (OL unstable)
– But the CLP is always the one in the LHP
– Explains result on 4–12 about why gain K �= 0 for OL unstable
systems, even for expensive control problem (Ruu → ∞)

June 18, 2008

Spr 2008 16.323 6–22
SRL Interpretations
• For SISO case, define Rzz/Ruu = 1/r.

• Consider what happens as r � ∞ – high control cost case

Δ(s) = D(s)D(−s) + r−1N (s)N (−s) = 0 ⇒ D(s)D(-s)=0
– So the n closed-loop poles are:
� Stable roots of the open-loop system (already in the LHP.)
� Reflection about the jω-axis of the unstable open-loop poles.

• Consider what happens as r � 0 – low control cost case

Δ(s) = D(s)D(−s) + r−1N (s)N (−s) = 0 ⇒ N(s)N(-s)=0
– Assume order of N (s)N (−s) is 2m < 2n
– So the n closed-loop poles go to:
� The m finite zeros of the system that are in the LHP (or the
reflections of the system zeros in the RHP).
� The system zeros at infinity (there are n − m of these).

• The poles tending to infinity do so along very specific paths so that

they form a Butterworth Pattern:
– At high frequency we can ignore all but the highest powers of s in
the expression for Δ(s) = 0
Δ(s) = 0 � (−1)ns2n + r−1(−1)m(bosm)2 = 0
2
2(n−m) n−m+1 bo
⇒ s = (−1)
r

June 18, 2008

Spr 2008 16.323 6–23

• The 2(n − m) solutions of this expression lie on a circle of radius

(b20/r)1/2(n−m)
at the intersection of the radial lines with phase from the negative
real axis:
lπ n−m−1
± , l = 0, 1, . . . , , (n-m) odd
n−m 2

(l + 1/2)π n−m
± , l = 0, 1, . . . , −1 , (n-m) even
n−m 2
n−m Phase
1 0
2 ±π/4
3 0, ±π/3
4 ±π/8, ±3π/8

• Note: Plot the SRL using the 180o rules (normal) if n − m is even
and the 0o rules if n − m is odd.
(s−2)(s−4)
Figure 6.2: G(s) = (s−1)(s−3)(s2 +0.8s+4)s2
Symmetric root locus
8

2
Imag Axis

−2

−4

−6

−8
−6 −4 −2 0 2 4 6
Real Axis

June 18, 2008

 Spr 2008 16.323 6–24

• As noted previously, we are free to pick the state weighting matrices

Cz to penalize the parts of the motion we are most concerned with.

• Simple example – consider oscillator with x = [ p , v ]T

� �
0 1
� �
0
A= , B=
−2 −0.5 1
but we choose two cases for z

� � � �
z=p= 1 0 x and z=v= 0 1 x
SRL with Position Penalty
SRL with Velocity Penalty

4
1.5

0.5
1
Imaginary Axis

Imaginary Axis

0 0

−1
−0.5

−2

−1
−3

−4 −1.5
−4 −3 −2 −1 0 1 2 3 4 −3 −2 −1 0 1 2 3
Real Axis Real Axis

Figure 6.3: SRL with position (left) and velocity penalties (right)

• Clearly, choosing a different Cz impacts the SRL because it completely

changes the zero-structure for the system.

June 18, 2008

Spr 2008 16.323 6–25
LQR Stability Margins
• LQR/SRL approach selects closed-loop poles that balance between
system errors and the control effort.
– Easy design iteration using r – poles move along the SRL.
– Sometimes difficult to relate the desired transient response to the
LQR cost function.
• Particularly nice thing about the LQR approach is that the designer
is focused on system performance issues

• Turns out that the news is even better than that, because LQR exhibits
very good stability margins
– Consider the LQR stability robustness.
� ∞
J = zT z + ρuT u dt
0
ẋ = Ax + Bu
z = Cz x, Rxx = CzT Cz

z
Cz �

u
� B (sI − A)−1 �
K �

–
x

• Study robustness in the frequency domain.

– Loop transfer function L(s) = K(sI − A)−1B
– Cost transfer function C(s) = Cz (sI − A)−1B

June 18, 2008

Spr 2008 16.323 6–26

• Can develop a relationship between the open-loop cost C(s) and the
closed-loop return difference I +L(s) called the Kalman Frequency
Domain Equality
1
[I + L(−s)]T [I + L(s)] = 1 + C T (−s)C(s)
ρ
• Sketch of Proof
– Start with u = −Kx, K = ρ1 B T P , where
1
0 = −AT P − P A − Rxx + P BB T P
ρ
– Introduce Laplace variable s using ±sP
1
0 = (−sI − AT )P + P (sI − A) − Rxx + P BB T P
ρ
– Pre-multiply by B T (−sI − AT )−1, post-multiply by (sI − A)−1B
– Complete the square . . .
1
[I + L(−s)]T [I + L(s)] = 1 + C T (−s)C(s)
ρ

• Can handle the MIMO case, but look at the SISO case to develop
further insights (s = jω)
[I + L(−s)]T [I + L(s)] = (I + Lr (ω) − jLi(ω))(I + Lr (ω) + jLi(ω))
≡ |1 + L(jω)|2
and
C T (−jω)C(jω) = Cr2 + Ci2 = |C(jω)|2 ≥ 0

• Thus the KFE becomes

1
|1 + L(jω)|2 = 1 + |C(jω)|2 ≥ 1
ρ

June 18, 2008

Spr 2008 16.323 6–27

• Implications: The Nyquist plot of L(jω) will always be outside the

unit circle centered at (-1,0).
4

3 |LN(jω)|

2 |1+LN(jω)|

1
Imag Part

0
(−1,0)
−1

−2

−3

−4
−7 −6 −5 −4 −3 −2 −1 0 1
Real Part

• Great, but why is this so significant? Recall the SISO form of the
Nyquist Stability Theorem:
If the loop transfer function L(s) has P poles in the RHP s-plane (and
lims→∞ L(s) is a constant), then for closed-loop stability, the locus
of L(jω) for ω : (−∞, ∞) must encircle the critical point (-1,0) P
times in the counterclockwise direction (Ogata528)
• So we can directly prove stability from the Nyquist plot of L(s).
But what if the model is wrong and it turns out that the actual loop
transfer function LA(s) is given by:
LA(s) = LN (s)[1 + Δ(s)], |Δ(jω)| ≤ 1, ∀ω

June 18, 2008

Spr 2008 16.323 6–28

• We need to determine whether these perturbations to the loop TF

will change the decision about closed-loop stability
⇒ can do this directly by determining if it is possible to change the
number of encirclements of the critical point
stable OL
3

|L|
1
|1+L|
Imag Part

0 ω=0

−1

ω
−2

−3
−2 −1 0 1 2 3 4
Real Part

Figure 6.4: Example of LTF for an open-loop stable system

• Claim is that “since the LTF L(jω) is guaranteed to be far from the
critical point for all frequencies, then LQR is VERY robust.”
– Can study this by introducing a modification to the system, where
nominally β = 1, but we would like to consider:
� The gain β ∈ R
� The phase β ∈ ejφ
K(sI − A)−1B � β �
– �

June 18, 2008

Spr 2008 16.323 6–29

• In fact, can be shown that:

– If open-loop system is stable, then any β ∈ (0, ∞) yields a stable
closed-loop system. For an unstable system, any β ∈ (1/2, ∞)
yields a stable closed-loop system ⇒ gain margins are (1/2, ∞)
– Phase margins of at least ±60◦

⇒ which are both huge.

Figure 6.5: Example loop transfer functions for open-loop stable system.

Figure 6.6: Example loop transfer functions for open-loop unstable system.

• While we have large margins, be careful because changes to some of

the parameters in A or B can have a very large change to L(s).

• Similar statements hold for the MIMO case, but it requires singular
value analysis tools.
June 18, 2008
Spr 2008 16.323 6–30

LTF for KDE

1 % Simple example showing LTF for KDE

2 % 16.323 Spring 2007
3 % Jonathan How
4 % rs2.m
5 %
6 clear all;close all;
7 set(0, ’DefaultAxesFontSize’, 14, ’DefaultAxesFontWeight’,’demi’)
8 set(0, ’DefaultTextFontSize’, 14, ’DefaultTextFontWeight’,’demi’)
9
10 a=diag([-.75 -.75 -1 -1])+diag([-2 0 -4],1)+diag([2 0 4],-1);
11 b=[
12 0.8180

13 0.6602

14 0.3420

15 0.2897];

16 cz=[ 0.3412 0.5341 0.7271 0.3093];

17 r=1e-2;

18 eig(a)

19 k=lqr(a,b,cz’*cz,r)

20 w=logspace(-2,2,200)’;w2=-w(length(w):-1:1);

21 ww=[w2;0;w];

22 G=freqresp(a,b,k,0,1,sqrt(-1)*ww);

23
24 p=plot(G);

25 tt=[0:.1:2*pi]’;Z=cos(tt)+sqrt(-1)*sin(tt);

26 hold on;plot(-1+Z,’r--’);plot(Z,’r:’,’LineWidth’,2);

27 plot(-1+1e-9*sqrt(-1),’x’)

28 plot([0 0]’,[-3 3]’,’k-’,’LineWidth’,1.5)

29 plot([-3 6],[0 0]’,’k-’,’LineWidth’,1.5)

30 plot([0 -2*cos(pi/3)],[0 -2*sin(pi/3)]’,’g-’,’LineWidth’,2)

31 plot([0 -2*cos(pi/3)],[0 2*sin(pi/3)]’,’g-’,’LineWidth’,2)

32 hold off

33 set(p,’LineWidth’,2);

34 axis(’square’)

35 axis([-2 4 -3 3])

36
37 ylabel(’Imag Part’);xlabel(’Real Part’);title(’Stable OL’)

38 text(.25,-.5,’\infty’)

39 print -dpng -r300 tf.png

40
41 %%%%%%%%%%%%%%%%%%%%%%

42
43 a=diag([-.75 -.75 1 1])+diag([-2 0 -4],1)+diag([2 0 4],-1);

44 r=1e-1;

45 eig(a)

46 k=lqr(a,b,cz’*cz,r)

47 G=freqresp(a,b,k,0,1,sqrt(-1)*ww);

48
49 p=plot(G);

50 hold on;plot(-1+Z,’r--’);plot(Z,’r:’,’LineWidth’,2);

51 plot(-1+1e-9*sqrt(-1),’x’)

52 plot([0 0]’,[-3 3]’,’k-’,’LineWidth’,1.5)

53 plot([-3 6],[0 0]’,’k-’,’LineWidth’,1.5)

54 plot([0 -2*cos(pi/3)],[0 -2*sin(pi/3)]’,’g-’,’LineWidth’,2)

55 plot([0 -2*cos(pi/3)],[0 2*sin(pi/3)]’,’g-’,’LineWidth’,2)

56 hold off

57 set(p,’LineWidth’,2)

58 axis(’square’)

59 axis([-3 3 -3 3])

60
61 ylabel(’Imag Part’);xlabel(’Real Part’);title(’Unstable OL’)

62 print -dpng -r300 tf1.png

June 18, 2008