16.333 Lecture # 10
State Space Control
? Basic state space control approaches
? ?
? ?
?
?
Fall 2004 16.333 9–1
State Space Basics
? State space models are of the form
x˙(t) = Ax(t) + Bu(t)
y(t) = Cx(t) + Du(t)
with associated transfer function
G(s) = C(sI ? A)
?1
B + D
Note: must form symbolic inverse of matrix (sI ? A), which is hard.
? Time response: Homogeneous part x˙ = Ax, x(0) known
– Take Laplace transform
X(s) = (sI ? A)
?1
x(0) ? x(t) = L
?1
(sI ? A)
?1
x(0)
I A A
2
– But can show (sI ? A)
?1
= +
s
2
+
s
3
+ . . .
s
1
so L
?1
(sI ? A)
?1
= I + At +
2!
(At)
2
+ . . . = e
At
– Gives x(t) = e
At
x(0) where e
At
is Matrix Exponential
1
3 Calculate in MATLAB
R?
using expm.m and not exp.m
? Time response: Forced Solution – Matrix case x˙ = Ax + Bu
where x is an n-vector and u is a m-vector. Cam show
t
x(t) = e
At
x(0) + e
A(t?τ )
Bu(τ )dτ
0
t
y(t) = Ce
At
x(0) + Ce
A(t?τ )
Bu(τ )dτ + Du(t)
0
– Ce
At
x(0) is the initial response
– Ce
A(t)
B is the impulse response of the system.
1
MATLAB
R?
is a trademark of the Mathworks Inc.
?
?
?
Fall 2004 16.333 9–2
Dynamic Interpretation
? Since A = T ΛT
?1
, then
? ?? ?? ?
T
e
At
| |
λ
1
t
.
.
.
? w
.
.
1
?
?
e = T e
Λt
T
?1
=
?
v
1
?? ??
.··· v
n
| |
λ
n
t T
e ? w
n
?
where we have written
? ?
T
? w
.
1
?
?.
T
?1
=
?
.
T
? w
n
?
which is a column of rows.
? Multiply this expression out and we get that
n
At λ
i
t T
e = e v
i
w
i
i=1
? Assume A diagonalizable, then x˙ = Ax, x(0) given, has solution
x(t) = e
At
x(0) = T e
Λt
T
?1
x(0)
n
= e
λ
i
t
v
i
{w
i
T
x(0)}
i=1
n
λ
i
t
= e v
i
β
i
i=1
? State solution is a linear combination of the system modes v
i
e
λ
i
e
λ
i
t
– Determines the nature of the time response
v
i
– Determines extent to which each state contributes to that mode
β
i
– Determines extent to which the initial condition excites the mode
? ?
? ?
Fall 2004 16.333 9–3
? Note that the v
i
give the relative sizing of the response of each part
of the state vector to the response.
1
v
1
(t) = e
?t
mode 1
0
0.5
v
2
(t) = e
?3t
mode 2
0.5
? Clearly e
λ
i
t
gives the time modulation
– λ
i
real – growing/decaying exponential response
– λ
i
complex – growing/decaying exponential damped sinusoidal
? Bottom line: The locations of the eigenvalues determine the pole
locations for the system, thus:
– They determine the stability and/or performance & transient be-
havior of the system.
– It is their locations that we will want to modify with the
controllers.
Fall 2004 16.333 9–4
Full-state Feedback Controller
? Assume that the single-input system dynamics are given by
x˙ = Ax + Bu
y = Cx
so that D = 0.
– The multi-actuator case is quite a bit more complicated as we
would have many extra degrees of freedom.
? Recall that the system poles are given by the eigenvalues of A.
– Want to use the input u(t) to modify the eigenvalues of A to
change the system dynamics.
r u
y
?
OO //
A, B, C
//
x
K
Assume a full-state feedback of the form:?
u = r ? Kx
where r is some reference input and the gain K is R
1×n
– If r = 0, we call this controller a regulator
? ?
? ?
? ?
? ?
Fall 2004 16.333 9–5
? Find the closed-loop dynamics:
x˙ = Ax + B(r ? Kx)
= (A ? BK)x + Br
= A
cl
x + Br
y = Cx
? Objective: Pick K so that A
cl
has the desired properties, e.g.,
– A unstable, want A
cl
stable
– Put 2 poles at ?2± 2j
? Note that there are n parameters in K and n eigenvalues in A, so it
looks promising, but what can we achieve?
? Example #1: Consider:
1 1
x˙ =
1 2
1
x + u
0
– Then
det(sI ? A) = (s ? 1)(s ? 2)? 1 = s
2
s + 1 = 0 ? 3
so the system is unstable.
– De?ne u = k
1
k
2
x = ?Kx, then ?
? ?
? ?
1 1
1
? ?
1? k
1
1? k
2
A
cl
= A?BK = k
1
k
2
=
1 2
?
0
1 2
– So then we have that
det(sI ? A
cl
) = s
2
+ (k
1
? 3)s + (1? 2k
1
+ k
2
) = 0
?
? ?
? ?
? ?
? ?
? ?
? ?
Fall 2004 16.333 9–6
– Thus, by choosing k
1
and k
2
, we can put λ
i
(A
cl
) anywhere in the
complex plane (assuming complex conjugate pairs of poles).
? To put the poles at s = ?5, ? 6, compare the desired characteristic
equation
(s + 5)(s + 6) = s
2
+ 11s + 30 = 0
with the closed-loop one
2
s + (k
1
? 3)x + (1? 2k
1
+ k
2
) = 0
to conclude that
k
1
? 3 = 11 k
1
= 14
1? 2k
1
+ k
2
= 30 k
2
= 57
so that K = 14 57 , which is called Pole Placement.
? Of course, it is not always this easy, as the issue of controllability
must be addressed.
? Example #2: Consider this system:
1 1
x˙ =
0 2
1
x + u
0
with the same control approach
1 1
1
? ?
1? k
1
1? k
2
A
cl
= A ? BK = k
1
k
2
=
0 2
?
0
0 2
so that det(sI ? A
cl
) = (s ? 1 + k
1
)(s ? 2) = 0
The feedback control can modify the pole at s = 1, but it cannot
move the pole at s = 2.
? This system cannot be stabilized with full-state feedback
control.
? What is the reason for this problem?
? ?
?
? ?
Fall 2004 16.333 9–7
– It is associated with loss of controllability of the e
2t
mode.
Basic test for controllability: rank M
c
=
n?
? ?
? ?
?
B AB
1 1
=
1
0
1
0
M
c
=
0 2
So that rank M
c
1
=
< 2.
? Must assume that the pair (A, B) are controllable.
? ?
? ?
?
?
? ?
?? ??
Fall 2004 16.333 9–8
Ackermann’s Formula
? The previous outlined a design procedure and showed how to do it
by hand for second-order systems.
– Extends to higher order (controllable) systems, but tedious.
? Ackermann’s Formula gives us a method of doing this entire design
process is one easy step.
K = 0 ... 0 1 M
?1
Φ
d
(A)
c
– M
c
= B AB ... A
n?1
B
– Φ
d
(s) is the characteristic equation for the closed-loop poles,
which we then evaluate for s = A.
– It is explicit that the system must be controllable because we
are inverting the controllability matrix.
? Revisit Example #1: Φ
d
(s) = s
2
+ 11s + 30
? ?
? ?
? ? ?
1 1
1
1 1
?
1
= B AB = =M
c
0
1 2
0
0 1
So
?
?
?
1 1
?
?1
?
? ?
2
? ?
?
1 1 1 1
K = 0 1
?
+ 11 + 30I
?
0 1 1 2 1 2
? ?
43 14
? ?
= 0 1 = 14 57
14 57
? Automated in Matlab: place.m & acker.m (see polyvalm.m too)
? ?
? ?
Fall 2004 16.333 9–9
? Origins? For simplicity, consider a third-order system (case #2), but
this extends to any order.
? ? ? ?
?a
1
?a
2
?a
3
1
A =
?
1 0 0
?
B =
?
0
?
C = b
1
b
2
b
3
0 1 0 0
– This form is useful because the characteristic equation for the
2
system is obvious ? det(sI ? A) = s
3
+ a
1
s + a
2
s + a
3
= 0
Can show that ?
? ? ? ?
?a
1
?a
2
?a
3
1
A
cl
= A ? BK =
?
1 0 0
?
?
?
0
?
k
1
k
2
k
3
0 1 0 0
? ?
?a
1
? k
1
?a
2
? k
2
?a
3
? k
3
=
?
1 0 0
?
0 1 0
so that the characteristic equation for the system is still obvious:
2
Φ
cl
(s) = det(sI ? A
cl
) = s
3
+ (a
1
+ k
1
)s + (a
2
+ k
2
)s + (a
3
+ k
3
) = 0
? We then compare this with the desired characteristic equation devel-
oped from the desired closed-loop pole locations:
2
Φ
d
(s) = s
3
+ (α
1
)s + (α
2
)s + (α
3
) = 0
to get that
?
a
1
+ k
1
= α
1 ?
k
1
= α
.
1
? a
1
.
. .
. .
?
a
n
+ k
n
= α
n
k
n
= α
n
? a
n
? Pole placement is a very powerful tool and we will be using it for
most of our state space work.
? ?
?
?
? ?
Fall 2004 16.333 9–10
Aircraft State Space Control
? Can now design a full state feedback controller for the dynamics:
x˙
sp
= A
sp
x
sp
+ B
sp
δ
e
with desired poles being at ω
n
= 3 and ζ = 0.6 ? s = ?1.8 ± 2.4i
φ
d
(s) = s
2
+ 3.6s + 9
Ksp=place(Asp,Bsp,[roots([1 2*0.6*3 3^2])’])
w
? Design controller u = ?0.0264 ?2.3463
q
? With full model, could arrange it so phugoid poles remain in the same
place, just move the ones associated with the short period mode
s = ?1.8 ± 2.4i, ?0.0033 ± 0.0672i
ev=eig(A);
% damp short period, but leave the phugoid where it is
Plist=[roots([1 2*.6*3 3^2])’ ev([3 4],1)’];
K1=place(A,B(:,1),Plist)
?
?
?
?
?
?
u
w
?
?
?
?
? u = 0.0026 ?0.0265 ?2.3428 0.0363
q
θ
?
?
?
?? ?
?
?
? ?
Fall 2004 16.333 9–11
? Can also add the lag dynamics to short period model with θ included
x˙
sp
= A
?
sp
x
sp
+
?
δ
a
δ
a
=
4
δ
c
B
sp
e
;
e
s + 4
e
→ x˙
δ
= ?4x
δ
+ 4δ
e
c
, δ
a
= x
δ
e
x˙
sp
=
?
A
?
sp
B
sp
x
sp
0
+ δ
c
e
?
x˙
δ
0 ?4 x
δ
4
? Add s = ?3 to desired pole list
Plist=[roots([1 2*.6*3 3^2])’,-.25,-3];
At2=[Asp2 Bsp2(:,1);zeros(1,3) -4];Bt2=[zeros(3,1);4];
Kt=place(At2,Bt2,Plist);
step(ss(At2-Bt2*Kt2,Bt2,[0 0 1 0],0),35)
?
?
?
?
?
?
w
q
θ
x
δ
?
?
?
?
u = 0.0011 ?3.4617 ?4.9124 0.5273
0 5 10 15 20 25 30 35?0.25
?0.2
?0.15
?0.1
?0.05
0 Step Response
Time (sec)
q
rads
? No problem working with larger systems with state space tools
? Main control issue is ?nding “good” locations for closed-loop poles
Fall 2004 16.333 9–12
Estimators/Observers
Problem: So far we have assumed that we have full access to the ?
state x(t) when we designed our controllers.
– Most often all of this information is not available.
? Usually can only feedback information that is developed from the
sensors measurements.
– Could try “output feedback”
?
u = Kx ? u = Ky
– Same as the proportional feedback we looked at at the beginning
of the root locus work.
– This type of control is very di?cult to design in general.
? Alternative approach: Develop a replica of the dynamic system
that provides an “estimate” of the system states based on the mea-
sured output of the system.
? New plan:
1. Develop estimate of x(t) that will be called x?(t).
2. Then switch from u = ?Kx(t) to u = ?Kx?(t).
? Two key questions:
– How do we ?nd x?(t)?
– Will this new plan work?
?
Fall 2004 16.333 9–13
Estimation Schemes
? Assume that the system model is of the form:
x˙ = Ax + Bu , x(0) unknown
y = Cx
where
1. A, B, and C are known.
2. u(t) is known
3. Measurable outputs are y(t) from C = I
? Goal: Develop a dynamic system whose state
x?(t) = x(t)
for all time t ≥ 0. Two primary approaches:
– Open-loop.
– Closed-loop.
Fall 2004 16.333 9–14
Open-loop Estimator
? Given that we know the plant matrices and the inputs, we can just
perform a simulation that runs in parallel with the system
x?
˙
(t) = Ax? + Bu(t)
? Then x?(t) ≡ x(t) ? t provided that x?(0) = x(0)
? Major Problem: We do not know x(0)
System A,B,C
? x(t)
y(t)
//
u(t)
//
//
Observer A,B,C
? ?x(t)
?y(t)
//
? Analysis of this case. Start with:
x˙(t) = Ax + Bu(t)
x?
˙
(t) = Ax? + Bu(t)
? De?ne the estimation error: x?(t) = x(t) ? x?(t).
– Now want x?(t) = 0 ? t.
– But is this realistic?
?
Fall 2004 16.333 9–15
? Subtract to get:
d
(x ?x?) = A(x ?x?) ? x?
˙
(t) = Ax?
dt
which has the solution
x?(t) = e
At
x?(0)
– Gives the estimation error in terms of the initial error.
? Does this guarantee that x? = 0 ? t?
Or even that x? → 0 as t →∞? (which is a more realistic goal).
– Response is ?ne if x?(0) = 0. But what if x?(0) = 0?
? If A stable, then x? → 0 as t →∞, but the dynamics of the estima-
tion error are completely determined by the open-loop dynamics of
the system (eigenvalues of A).
– Could be very slow.
– No obvious way to modify the estimation error dynamics.
? Open-loop estimation does not seem to be a very good idea.
Fall 2004 16.333 9–16
Closed-loop Estimator
? Obvious way to ?x the problem is to use the additional information
available:
– How well does the estimated output match the measured output?
Compare: y = Cx with y? = Cx?
– Then form y? = y ? y? ≡ Cx?
System A,B,C
→ x(t)
y(t)
//
+
u(t)
//
//
L
Observer A,B,C
→ ?x(t)
?y(t)
//
?
OO
? Approach: Feedback y? to improve our estimate of the state. Basic
form of the estimator is:
x?
˙
(t) = Ax?(t) + Bu(t) + Ly?(t)
y?(t) = Cx?(t)
where L is a user selectable gain matrix.
? Analysis:
˙
x? = x˙ ? x?
˙
= [Ax + Bu] [Ax? + Bu + L(y ? y?)]?
= A(x ? x?) ? L(Cx ? Cx?)
= Ax? ? LCx? = (A ? LC)?x
Fall 2004 16.333 9–17
? So the closed-loop estimation error dynamics are now
x?
˙
= (A ?LC)?x with solution x?(t) = e
(A?LC)t
x?(0)
? Bottom line: Can select the gain L to attempt to improve the
convergence of the estimation error (and/or speed it up).
– But now must worry about observability of the system model.
? Note the similarity:
– Regulator Problem: pick K for A ?BK
3 Choose K ∈R
1×n
(SISO) such that the closed-loop poles
det(sI ?A + BK) = Φ
c
(s)
are in the desired locations.
– Estimator Problem: pick L for A ?LC
3 Choose L ∈R
n×1
(SISO) such that the closed-loop poles
det(sI ?A + LC) = Φ
o
(s)
are in the desired locations.
? These problems are obviously very similar – in fact they are called
dual problems.
Fall 2004 16.333 9–18
Estimation Gain Selection
? For regulation, were concerned with controllability of (A, B)
For a controllable system we can place the eigenvalues of
A ? BK arbitrarily.
? For estimation, were concerned with observability of pair (A, C).
For an observable system we can place the eigenvalues of
A ? LC arbitrarily.
? Test using the observability matrix:
??
rank M
o
? rank
?
?
?
?
?
?
C
CA
CA
2
.
.
.
CA
n?1
?
?
?
?
?
?
= n
? The procedure for selecting L is very similar to that used for the
regulator design process.
Fall 2004 16.333 9–19
? One approach:
– Note that the poles of (A ? LC) and (A ? LC)
T
are identical.
– Also we have that (A ? LC)
T
= A
T
? C
T
L
T
– So designing L
T
for this transposed system looks like a standard
regulator problem (A ? BK) where
A A
T
?
B C
T
?
K L
T
?
So we can use
K
e
= acker(A
T
, C
T
, P ) , L ≡ K
T
e
? Note that the estimator equivalent of Ackermann’s formula is that
?
?
L = Φ
e
(s)M
?1
o
?
?
?
?
0
.
.
.
0
1
?
?
?
?
? ?
? ?
? ?
? ?
Fall 2004 16.333 9–20
Simple Estimator Example
? Simple system
? ?
? ? ? ?
?1 1.5
1 .5
A = , B = , x(0) =
?0
1 ?2
0 ?1
C = 1 0 , D = 0
– Assume that the initial conditions are not well known.
– System stable, but λ
max
(A) = ?0.18
– Test observability:
C
rank
CA
1 0
= rank
?1 1.5
? Use open and closed-loop estimators. Since the initial conditions are
0
not well known, use x?(0) =
0
? Open-loop estimator:
x?
˙
= Ax? + Bu
y? = Cx?
? Closed-loop estimator:
x?
˙
= Ax? + Bu + Ly? = Ax? + Bu + L(y ? y?)
= (A ? LC)?x + Bu + Ly
y? = Cx?
– Which is a dynamic system with poles given by λ
i
(A ? LC) and
which takes the measured plant outputs as an input and generates
an estimate of x.
? ?
? ?
? ? ? ? ? ?
?
?
?
? ?
? ?
? ?
? ?
?
?
Fall 2004 16.333 9–21
? Typically simulate both systems together for simplicity
? Open-loop case:
x˙ = Ax + Bu
y = Cx
x?
˙
= Ax? + Bu
y? = Cx?
?
?
?0.5
?1
x˙
A 0
x
?
?
?
?
?
?
?
?
B x(0)
= + u , =
x?(0)
?
x?
˙
0 A
x? B
y
?
C
0
x
=
y?
0 C
x?
? Closed-loop case:
x˙ = Ax + Bu
x?
˙
= (A ? LC)?x + Bu + LCx
x˙
A 0
x B
= u
B
?
x?
˙
LC A ? LC
x?
+
? Example uses a strong u(t) to shake things up
0
0
Fall 2004 16.333 9–22
0 0.5 1 1.5 2 2.5 3 3.5 4?1
?0.5
0
0.5
1 Open loop estimator
states
time
x1 x2
0 0.5 1 1.5 2 2.5 3 3.5 4?1
?0.5
0
0.5
1
estimation error
time
Figure 1: Open-loop estimator. Estimation error converges to zero, but very slowly.
0 0.5 1 1.5 2 2.5 3 3.5 4?1
?0.5
0
0.5
1 Closed?loop estimator
states
time
x1 x2
0 0.5 1 1.5 2 2.5 3 3.5 4?1
?0.5
0
0.5
1
estimation error
time
Figure 2: Closed-loop estimator. Convergence looks much better.
?
Fall 2004 16.333 9–23
Aircraft Estimation Example
? Take Short period model and assume that we can measure q. Can
we estimate the motion associated with the short period mode?
x˙
sp
= A
sp
x
sp
?
+ B
sp
u
y = 0 1 x
sp
– Take x
sp
(0) = [?0.5; ?0.05]
T
? System stable, so could use an open loop estimator
? For closed-loop estimator, put desired poles at ?3, ?4
? For the various dynamics models as before
Csp=[0 1]; % sense q
Ke=place(Asp’,Csp’,[-3 -4]);Le=Ke’;
0 2 4 6 8 10?20
?15
?10
?5
0
5
10
x1
Open?loop estimator
time 0 2 4 6 8 10
?0.1
?0.05
0
0.05
0.1
x1
time
0 2 4 6 8 10?8
?6
?4
?2
0
2
4
x1 error
time 0 2 4 6 8 10
?0.05
?0.04
?0.03
?0.02
?0.01
0
0.01
0.02
x2 error
time
0 2 4 6 8 10?20
?15
?10
?5
0
5
10
15
x1
Closed?loop estimator
time 0 2 4 6 8 10
?0.1
?0.05
0
0.05
0.1
x1
time
0 2 4 6 8 10?20
?15
?10
?5
0
x1 error
time 0 2 4 6 8 10
?0.05
?0.04
?0.03
?0.02
?0.01
0
0.01
x2 error
time
Figure 3: Closed-loop estimator. Convergence looks much better.
? As expected, the OL estimator does not do well, but the closed-loop
one converges nicely
Fall 2004 16.333 9–24
Where to put the Estimator Poles?
? Location heuristics for poles still apply – use Bessel, ITAE, ...
– Main di?erence: probably want to make the estimator faster than
you intend to make the regulator – should enhance the control,
which is based on x?(t).
– ROT: Factor of 2–3 in the time constant ζω
n
associated with the
regulator poles.
? Note: When designing a regulator, were concerned with “band-
width” of the control getting too high ? often results in control
commands that saturate the actuators and/or change rapidly.
Di?erent concerns for the estimator: ?
– Loop closed inside computer, so saturation not a problem.
– However, the measurements y are often “noisy”, and we need to
be careful how we use them to develop our state estimates.
? High bandwidth estimators tend to accentuate the e?ect of sens-
ing noise in the estimate.
– State estimates tend to “track” the measurements, which are ?uc-
tuating randomly due to the noise.
? Low bandwidth estimators have lower gains and tend to rely more
heavily on the plant model
– Essentially an open-loop estimator – tends to ignore the measure-
ments and just uses the plant model.
Fall 2004 16.333 9–25
? Can also develop an optimal estimator for this type of system.
– Which is apparently what Kalman did one evening in 1958 while
taking the train from Princeton to Baltimore...
– Balances e?ect of the various types of random noise in the
system on the estimator:
x˙ = Ax + Bu + B
w
w
y = Cx + v
where:
3 w: “process noise” – models uncertainty in the system model.
3 v: “sensor noise” – models uncertainty in the measurements.
Final Thoughts
? Note that the feedback gain L in the estimator only stabilizes the
estimation error.
– If the system is unstable, then the state estimates will also go to
∞, with zero error from the actual states.
? Estimation is an important concept of its own.
– Not always just “part of the control system”
– Critical issue for guidance and navigation system
? More complete discussion requires that we study stochastic processes
and optimization theory.
? Estimation is all about which do you trust more: your mea-
surements or your model.
? ?
? ? ? ?
? ?
Fall 2004 16.333 9–26
Combined Regulator and Estimator
? As advertised, we can change the previous control u = ?Kx to the
new control u = ?Kx? (same K). We now have
x˙ = Ax + Bu
y = Cx
x?
˙
= Ax? + Bu + L(y ? y?)
y? = Cx?
with closed-loop dynamics
? ? ? ?? ?
x˙ A ?BK x
x?
˙
=
LC A ? BK ? LC x?
? x˙
cl
= A
cl
x
cl
? Not obvious that this system will even be stable: λ
i
(A
cl
) < 0?
? To analyze, introduce x? = x ? x?, and the similarity transform
I 0
T = = T
?1
I ?I
x x
? Rewrite the dynamics in terms of the state = T
x? x?
A
cl
? T
?1
A
cl
T ≡ A
cl
and when you work through the math, you get
A
cl
=
A ? BK BK
!!!
0 A ? LC
Fall 2004 16.333 9–27
? Absolutely key points:
1. λ
i
(A
cl
) ≡ λ
i
(A
cl
) why?
2. A
cl
is block upper triangular, so can ?nd poles by inspection:
det(sI ? A
cl
) = det(sI ? (A ? BK)) · det(sI ? (A ? LC))
The closed-loop poles of the system consist of the
union of the regulator and estimator poles
? So we can design the estimator and regulator separately with con?-
dence that combination of the two will work VERY well.
? Compensator is a combination of the estimator and regulator.
x?
˙
= Ax? + Bu + L(y ? y?)
= (A ? BK ? LC)?x + Ly
u = ?Kx?
? x˙
c
= A
c
x
c
+ B
c
y
u = ?C
c
x
c
– Keep track of this minus sign. We need one in the feedback path,
but we can move it around to suit our needs.
Fall 2004 16.333 9–28
? Let G
c
(s) be the compensator transfer function where
u
= ?C
c
(sI ? A
c
)
?1
B
c
= ?G
c
(s)
y
= ?K(sI ? (A ? BK ? LC))
?1
L
so by my de?nition, ? u = ?G
c
y ≡ G
c
(?y)
? Reason for making the de?nition is that when we implement the
controller, we often do not just feedback ?y(t), but instead have to
include a reference command r(t)
– Use servo approach and feed back e(t) = r(t)? y(t) instead
r
//
e
//
G
c
(s)
u
//
G(s)
y
//
?
OO
– So now u = G
c
e = G
c
(r ? y).
– And if r = 0, then we still have u = G
c
(?y)
? Important points:
– Closed-loop system will be stable, but the compensator dynamics
need not be.
– Often very simple and useful to provide classical interpretations of
the compensator dynamics G
c
(s).
? ?
? ?
Fall 2004 16.333 9–29
? Mechanics of closing the loop
G(s) : ˙x = Ax + Bu
y = Cx
G
c
(s) : ˙x
c
= A
c
x
c
+ B
c
e
u = C
c
x
c
and e = r ? y, u = G
c
e, y = Gu.
? Loop dynamics L = GG
c
? y = L(s)e
x˙ = Ax + Bu = Ax + BC
c
x
c
x˙
c
= A
c
x
c
+ B
c
e
? ? ? ?? ? ? ?
x˙ A BC
c
x 0
= + e
x˙
c
0 A
c
x
c
B
c
? ?
x
y = C 0
x
c
? Now form the closed-loop dynamics by inserting e = r ? y
? ? ? ?? ? ? ?? ? ??
x˙ A BC
c
x 0
? ?
x
=
x˙
c
0 A
c
x
c
+
B
c
r ? C 0
x
c
? ?? ? ? ?
A BC
c
x 0
= + r
?B
c
C A
c
x
c
B
c
? ?
x
y = C 0
x
c
?
?
? ? ? ?
?
?
Fall 2004 16.333 9–30
Performance Issue
? Often ?nd with state space controllers that the DC gain of the closed
loop system is not 1. So y = r in steady state.
? Relatively simple ?x is to modify the original controller with scalar N
u = r ? Kx ? u = N r ? Kx
? Closed-loop system on page 5 becomes
x˙ = Ax + B(N r ? Kx) = A
cl
x + BN r
G
cl
(s) = C(sI?A
cl
)
?1
BN
y = Cx
– Analyze steady state step response ? y
ss
= G
cl
(0)r
step
G
cl
(0) = C(?A
cl
)
?1
BN
– And pick N so that G
cl
(0) = 1 ? N =
1
(C(?A
cl
)
?1
B)
? A bit more complicated with a combined estimator and regulator
– One simple way (not the best) of achieving a similar goal is to add
N to r and force G
cl
(0) = 1
– Now the closed-loop dynamics on page 29 become:
?
?
?
?
?
x˙ x
= A
cl
+ B
cl
N r
x˙
c
1
N = →
(C
cl
(?A
cl
)
?1
B
cl
)
x
c
?
?
?
?
x
y = C
cl
x
c
? Note that this ?xes the steady state tracking error problems, but in
my experience can create strange transients (often NMP).
? ?
? ? ? ?
Fall 2004 16.333 9–31
Example: Compensator Design
1 x˙ = Ax + Bu
G(s) =
s
2
+ s + 1
?
y = Cx
where
?
0 1
? ?
0
?
? ?
A = B = C = 1 0
?1 ?1 1
? Regulator: Want regulator poles to have a time constant of τ
c
=
1/(ζω
n
) = 0.25 sec ? λ(A ? BK
r
) = ?4 ± 4j which can be found
using place or acker
K_r=acker(a,b,[-4+4*j;-4-4*j]);
to give K
r
= [ 31 7 ]
? Estimator: want the estimator poles to be faster, so use
τ
e
= 1/(ζω
n
) = 0.1 sec. Use real poles, ? λ(A ? L
e
C) = ?10
L_e=acker(a’,c’,[-10 -10]’)’;
19
which gives L
e
=
80
? Form compensator G
c
(s)
ac=a-b*K_r-L_e*c;bc=L_e;cc=K_r;dc=0;
A
c
=
?19 1 19
? ?
B
c
= C
c
= 31 7
?112 80?8
(s + 2.5553) u
G
c
(s) = 1149 =
s
2
+ 27s + 264 e
Low frequency zero, with higher frequency poles (like a lead)
Fall 2004 16.333 9–32
10?1 100 101 102
100
101
102
Freq (rad/sec)
Mag
Plant G
Compensator Gc
10?1 100 101 102
?200
?150
?100
?50
0
50
Freq (rad/sec)
Phase (deg)
Plant G
Compensator Gc
Figure 4: The compensator does indeed look like a high frequency lead (ampli?cation
from 2–16 rad/sec). Plant pretty simple looking.
Fall 2004 16.333 9–33
10?1 100 101 102
100
101
102
Freq (rad/sec)
Mag
Loop L
10?1 100 101 102
?250
?200
?150
?100
?50
0
Freq (rad/sec)
Phase (deg)
Figure 5: The loop transfer function L = G
c
G shows a slope change around ω
c
= 5
rad/sec due to the e?ect of the compensator. Signi?cant gain and phase margins.
Fall 2004 16.333 9–34
?150
?100
?50
0
50
Magnitude (dB)
10?2 10?1 100 101 102 103
?270
?225
?180
?135
?90
?45
0
Phase (deg)
Bode Diagram
Gm = 14.3 dB (at 14.9 rad/sec) , Pm = 45.8 deg (at 4.84 rad/sec)
Frequency (rad/sec)
Figure 6: Quite signi?cant gain and phase margins.
Fall 2004 16.333 9–35
?40 ?35 ?30 ?25 ?20 ?15 ?10 ?5 0 5 10
?25
?20
?15
?10
?5
0
5
10
15
20
25
Root Locus
Real Axis
Imaginary Axis
Figure 7: Freeze the compensator poles and zeros and draw a root locus versus an
?
additional plant gain α, G(s) ? G(s) =
(s
2
+
α
s+1)
. Note location of the closed-loop
poles!!
Fall 2004 16.333 9–36
10?1 100 101 102
10?2
10?1
100
101
102
Freq (rad/sec)
Mag
Factor of N=1.0899 applied to Closed loop
Plant G
closed?loop Gcl
Figure 8: Closed-loop transfer – system bandwidth has increased substantially.
Fall 2004 16.333 9–37
Estimator Design (est1.m)
1 clear all
2 close all
3 figure(1);clf
4 set(gcf,’DefaultLineLineWidth’,2)
5 set(gcf,’DefaultlineMarkerSize’,10)
6 figure(2);clf
7 set(gcf,’DefaultLineLineWidth’,2)
8 set(gcf,’DefaultlineMarkerSize’,10)
9
10 load b747 % get A B Asp Bsp
11 Csp=[0 1]; % sense q
12
13 Ke=place(Asp’,Csp’,[-3 -4]);Le=Ke’;
14
15 xo=[-.5;-.05]; % start somewhere
16
17 t=[0:.01:10];N=floor(.15*length(t));
18 % hit on the system with an input
19 %u=0;u=[ones(15,1);-ones(15,1);ones(15,1)/2;-ones(15,1)/2;zeros(41,1)]/5;
20 u=0;u=[ones(N,1);-ones(N,1);ones(N,1)/2;-ones(N,1)/2]/20;
21 u(length(t))=0;
22
23 [y,x]=lsim(Asp,Bsp,Csp,0,u,t,xo);
24 plot(t,y)
25
26 % closed-loop estimator
27 % hook both up so that we can simulate them at the same time
28 % bigger state = state of the system then state of the estimator
29 A_cl=[Asp zeros(size(Asp));Le*Csp Asp-Le*Csp];
30 B_cl=[Bsp;Bsp];
31 C_cl=[Csp zeros(size(Csp));zeros(size(Csp)) Csp];
32 D_cl=zeros(2,1);
33
34 % note that we start the estimators at zero, since that is
35 % our current best guess of what is going on (i.e. we have no clue :-) )
36 %
37 [y_cl,x_cl]=lsim(A_cl,B_cl,C_cl,D_cl,u,t,[xo;0;0]);
38 figure(1)
39 subplot(221)
40 plot(t,x_cl(:,[1]),t,x_cl(:,[3]),’--’)
41 ylabel(’x1’);title(’Closed-loop estimator’);xlabel(’time’);grid
42 subplot(222)
43 plot(t,x_cl(:,[2]),t,x_cl(:,[4]),’--’)
44 ylabel(’x1’);xlabel(’time’);grid
45 subplot(223)
46 plot(t,x_cl(:,[1])-x_cl(:,[3]))
47 ylabel(’x1 error’);xlabel(’time’);grid
48 subplot(224)
49 plot(t,x_cl(:,[2])-x_cl(:,[4]))
50 ylabel(’x2 error’);xlabel(’time’);grid
51 print -depsc spest_cl.eps
52 jpdf(’spest_cl’)
53
54 % open-loop estimator
55 % hook both up so that we can simulate them at the same time
56 % bigger state = state of the system then state of the estimator
57 A_ol=[Asp zeros(size(Asp));zeros(size(Asp)) Asp];
58 B_ol=[Bsp;Bsp];
59 C_ol=[Csp zeros(size(Csp));zeros(size(Csp)) Csp];
60 D_ol=zeros(2,1);
61
62 [y_ol,x_ol]=lsim(A_ol,B_ol,C_ol,D_ol,u,t,[xo;0;0]);
63 figure(2)
64 subplot(221)
65 plot(t,x_ol(:,[1]),t,x_ol(:,[3]),’--’)
66 ylabel(’x1’);title(’Open-loop estimator’);xlabel(’time’);grid
Fall 2004 16.333 9–38
67 subplot(222)
68 plot(t,x_ol(:,[2]),t,x_ol(:,[4]),’--’)
69 ylabel(’x1’);xlabel(’time’);grid
70 subplot(223)
71 plot(t,x_ol(:,[1])-x_ol(:,[3]))
72 ylabel(’x1 error’);xlabel(’time’);grid
73 subplot(224)
74 plot(t,x_ol(:,[2])-x_ol(:,[4]))
75 ylabel(’x2 error’);xlabel(’time’);grid
76 print -depsc spest_ol.eps
77 jpdf(’spest_ol’)
78
Fall 2004 16.333 9–39
Regular/Estimator Design (reg est.m)
1 % Combined estimator/regulator design for a simple system
2 % G= 1/(s^2+s+1)
3 %
4 % Jonathan How
5 % Fall 2004
6 %
7 close all;clear all
8 for ii=1:5
9 figure(ii);clf;set(gcf,’DefaultLineLineWidth’,2);set(gcf,’DefaultlineMarkerSize’,10)
10 end
11
12 a=[0 1;-1 -1];b=[0 1]’;c=[1 0];d=0;
13 k=acker(a,b,[-4+4*j;-4-4*j]);
14 l=acker(a’,c’,[-10 -10]’)’;
15 %
16 % For state space for G_c(s)
17 %
18 ac=a-b*k-l*c;bc=l;cc=k;dc=0;
19
20 G=ss(a,b,c,d);
21 Gc=ss(ac,bc,cc,dc);
22
23 f=logspace(-1,2,400);
24 g=freqresp(G,f*j);g=squeeze(g);
25 gc=freqresp(Gc,f*j);gc=squeeze(gc);
26
27 figure(1);clf
28 subplot(211)
29 loglog(f,abs(g),f,abs(gc),’--’);axis([.1 1e2 .2 1e2])
30 xlabel(’Freq (rad/sec)’);ylabel(’Mag’)
31 legend(’Plant G’,’Compensator Gc’);grid
32 subplot(212)
33 semilogx(f,180/pi*angle(g),f,180/pi*angle(gc),’--’);
34 axis([.1 1e2 -200 50])
35 xlabel(’Freq (rad/sec)’);ylabel(’Phase (deg)’);grid
36 legend(’Plant G’,’Compensator Gc’)
37
38 L=g.*gc;
39
40 figure(2);clf
41 subplot(211)
42 loglog(f,abs(L),[.1 1e2],[1 1]);axis([.1 1e2 .2 1e2])
43 xlabel(’Freq (rad/sec)’);ylabel(’Mag’)
44 legend(’Loop L’);
45 grid
46 subplot(212)
47 semilogx(f,180/pi*phase(L.’),[.1 1e2],-180*[1 1]);
48 axis([.1 1e2 -290 0])
49 xlabel(’Freq (rad/sec)’);ylabel(’Phase (deg)’);grid
50 %
51 % loop dynamics L = G Gc
52 %
53 al=[a b*cc;zeros(2) ac];
54 bl=[zeros(2,1);bc];
55 cl=[c zeros(1,2)];
56 dl=0;
57 figure(3)
58 rlocus(al,bl,cl,dl)
59 %
60 % closed-loop dynamics
61 % unity gain wrapped around loop L
62 %
63 acl=al-bl*cl;bcl=bl;ccl=cl;dcl=d;
64
65 N=inv(ccl*inv(-acl)*bcl)
66
Fall 2004 16.333 9–40
67 hold on;plot(eig(acl),’d’);hold off
68 grid
69 %
70 % closed-loop freq response
71 %
72 Gcl=ss(acl,bcl*N,ccl,dcl);
73 gcl=freqresp(Gcl,f*j);gcl=squeeze(gcl);
74
75 figure(4);clf
76 loglog(f,abs(g),f,abs(gcl),’--’);
77 axis([.1 1e2 .01 1e2])
78 xlabel(’Freq (rad/sec)’);ylabel(’Mag’)
79 legend(’Plant G’,’closed-loop G_{cl}’);grid
80 title([’Factor of N=’,num2str(N),’ applied to Closed loop’])
81
82 figure(5);clf
83 margin(al,bl,cl,dl)
84
85 figure(1);orient tall;print -depsc reg_est1.eps
86 jpdf(’reg_est1’)
87 figure(2);orient tall;print -depsc reg_est2.eps
88 jpdf(’reg_est2’)
89 figure(3);print -depsc reg_est3.eps
90 jpdf(’reg_est3’)
91 figure(4);print -depsc reg_est4.eps
92 jpdf(’reg_est4’)
93 figure(5);print -depsc reg_est5.eps
94 jpdf(’reg_est5’)