486133.txt

Computer-Controlled Systems
Solutions Manual
Karl J. Åström
Björn Wittenmark
Department of Automatic Control
Lund Institute of Technology
October 1997
Preface
This Solutions Manual contains solutions to most of the problems in the fourth
edition of
Åström, K. J. and B. Wittenmark (1997): Computer controlled Systems –
heory and Applications, Prentice Hall Inc., Englewood Cliffs, N. J.
Many of the problems are intentionally made such that the students have to
use a simulation program to verify the analytical solutions. This is important
since it gives a feeling for the relation between the pulse transfer function and
the time domain. In the book and the solutions we have used Matlab/Simulink.
Information about macros used to generate the illustrations can be obtained by
It is also important that a course in digital control includes laboratory exercises.
The contents in the laboratory experiments are of course dependent on the available equipment. Examples of experiments are
Illustration of aliasing
Comparison between continuous time and discrete time controllers
State feedback control. Redesign of continuous time controllers as well as
controllers based on discrete time synthesis
Controllers based on input-output design
Control of systems subject to stochastic disturbances.
Finally we would like to thank collegues and students who have helped us to test
the book and the solutions.
Karl J. Åström
Björn Wittenmark
Department of Automatic Control
Lund Institute of Technology
S-220 00 Lund, Sweden
Solutions to Chapter 2
Problem 2.1
−ax + bu
The system is described by
Sampling the system using (2.4) and (2.5) gives
e−ah x( kh) + 1
x( kh + h)
y( kh)
cx( kh)
− e−
u( kh)
a
ah
The pole of the sampled system is exp( ah). The pole is real. For small values
of h the pole is close to 1. If a > 0 then the pole moves towards the origin when
h increases. If a < 0 then the pole moves along the positive real axis.
Problem 2.2
a. Using the Laplace transform method we find that
e
e B ds
b. The system has the transfer function
G ( s)
( s + 1)( s + 2)
Using the Table 2.1 gives
H ( q)
− e− − 1 − e−
− e− 2 q − e−
One state space realization of the system is
then
e Ah
e B ds
and
H ( q)
C ( qI
(q
h ( q + 1)/2 h( q 1) ( q 1)
h3 ( q2 + 4q + 1)
6( q 1)3
Problem 2.3
a.
e Ah
y( k)
− 0.5 y(k − 1)
6u( k
b as
e
a
eas b ds
6a
ln 0.5 ⇒ a
eah
ax( t) + bu( t)
as
e b ds
x( t)
y( t)
x( t)
(continuous time system)
ah
0.5x( kh) + 6u( kh)
y( kh)
x( kh)
(discrete-time system)
eah
x( kh + h)
b ah
e
a
 x( kh) + 
 u( kh)
 x( kh + h) 
 y( kh)  1 1  x( kh)
Eigenvalue to Φ :
det( sI Φ) 
0 ⇔ ( s + 0.5)( s + 0.3)
Both eigenvalues of Φ on the negative real axis ⇒ No corresponding continuous system exists.
y( k) + 0.5 y( k
−0.5 y(k) + 6u(t)
y( k + 1)
H ( q)
6u( k
one pole on the negative real axis.
⇒ No equivalent continuous system exists.
Problem 2.4
Harmonic oscillator (cf. A.3 and 3.2a).
Φ x( k) + Γ u( k)
C x( k)
x( k + 1)
y( k)
Φ( h)
Γ( h)
a.
Pulse transfer operator
H ( q)
C ( qI
Y ( z)
y( k)
where θ ( k
H ( z) U ( z)
(k
− 1)θ (k − 1) + θ (k − 1)
1) is a step at k
G ( s)
Y ( s)
s( s2 + 1)
y( kh)
[see Probl. 2.2]
y( t)
θ (k
The same way as a. ⇒ y( t) 1 cos t, and y( kh) 1 cos π k Notice that
the step responses of the continuous time and the zero-order hold sampled
systems are the same in the sampling points.
Problem 2.5
Do a partial fraction decomposition and sample each term using Table 2.1:
G ( s)
H ( q)
− 36 q − 1 + 1 1 − e− − 27 q − e−
12 ( q − 1)
8 q− e
−e
s2 ( s + 2)( s + 3)
Problem 2.6
Integrating the system equation
gives
x( kh + h)
e
e
e As B δ ( s
x( kh) +
− kh)u(kh)ds
x( kh) + Bu( kh)
Problem 2.7
The representation (2.7) is
x( kh + h)
y( kh)
Φ x( kh) + Γ u( kh)
C x( kh)
and the controllable form realization is
z( kh + h)
y( kh)
where
Φ z( kh) + Γ u( kh)
C z( kh)
From Section 2.5 we get
This gives the following relations
(1)
(2)
(3)
(4)
Equations (1)–(4) now give
1/(2h)
Problem 2.8
The pulse transfer function is given by
z( z 0.5) 0
C ( zI − Φ)−1 Γ
H ( z)
z( z
2( z
z( z
Problem 2.9
The system is described by
 a
The eigenvalues of A is obtained from
(λ + a)(λ + d)
which gives
λ 2 + ( a + d)λ + ad
−a + d ±
(a
The condition that a b c, and d are nonnegative implies that the eigenvalues, λ 1
and λ 2 , are real. There are multiple eigenvalues if both a d and bc 0. Using
the result in Appendix B we find that
e Ah
eλ 1 h
and
e
This gives
λ 1 eλ 2 h
−λ e
eλ 1 h eλ 2 h
(λ 1 λ 2 ) h
− α ah
To compute Γ we notice that α 0 and α 1 depend on h. Introduce
α 0 ( s) ds
e
sα 1 ( s) ds
then
eλ 1 h
eλ 1 h
eλ 2 h
Problem 2.10
a.
Using the result from Problem 2.9 gives
eλ 1 h
eλ 2 h
Further
and
(β 0 I + β 1 A) B
The pulse transfer operator is now given by
H ( q)
C ( qI − Φ)−1 Γ
which agrees with the pulse transfer operator given in the problem formulation.
Problem 2.11
The motor has the transfer function G ( s)
tation is given in Example A.2, where
exp( Ah)
e−h
1 − e−h
e B ds
( sI − A)−1
s( s + 1) 1
− e− 
h + e− − 1
1/( s( s + 1)). A state space represen-
e−s
e−s
This gives the sampled representation given in Example A.2.
a.
The pulse transfer function is
H ( z)
C ( zI
  z e−h
0 −1  1 e−h 
1 + e−h z 1
h + e−h 1
0   1 e−h 
( z e−h)( z 1) 1 e−h z e−h
h + e−h 1
− − 1)z + (1 − e− − he− )
(h + e
( z − e− )( z − 1)
( h + e− − 1) z + (1 − e− − he− )
z − (1 + e− ) z + e−
The pulse response is
h( k )
e Ah
e Akh
This gives
− e−
h − e− −
− e− −
1 − e− −
( k 1) h
e−( k−1) h
0   1 e−h 
h + e−h 1
( k 1) h
( k 1) h
+ h + e−h
+ e−kh
An alternative way to find h( k) is to take the inverse z-transform of H ( z).
A difference equation is obtained from H ( q)
y( kh + 2h)
− (1 + e− ) y(kh + h) + e− y(kh)
( h + e− − 1)u( kh + h) + (1 − e− − he− )u( kh)
1 and z exp(−h). The second pole will move from 1 to
The poles are in z
the origin when h goes from zero to infinity.
The zero is in
1 e−h he−h
f ( h)
h + e−h 1
The function f ( h) goes from
−1 to 0 when h goes from zero to infinity. See Fig. 2.1.
Problem 2.12
Consider the following transfer operators of infinity order.
G ( s)
e
(τ < h)
u( t
0 ⋅ x + 1 ⋅ u( t
(infinite order)
Pole and zero
Sampling period h
Figure 2.1 The pole (solid) and zero (dashed) in Problem 2.11.d as function of the
sampling period.
a.
Discrete time system is given by (cf. CCS p. 38–40).
Φ x( k) + Γ 0 u( k) + Γ 1 u( k
x( k + 1)
e Ah
e0
e ds B
e
A( h τ )
e ds B
⇒ x( k + 1)
− 1) (Notice that τ < h)
x( k) + 0.5u( k) + 0.5u( k
State space model:
 x( k + 1) 
u( k )
 1 0.5   x( k)   0.5 
 u( k )
u( k 1 )
  x( k) 
y( k)  1 0  
u( k 1 )
Φ x( k)+ Γ u( k)
The system is of second order.
The pulse-transfer function is
H ( z)
C ( zI
z( z 1 ) 0 z 1
z( z 1 )
0.5( z + 1)
z( z 1 )
Invers transform of H ( z) using Table 2.3 gives
0.5( z + 1)
z( z 1 )
H ( z)
The poles are z
⇒ h( kh)
1 and the zeros: z
Problem 2.13
a.
This is the same system as in Example 2.6 but with τ
have d 2 and τ ′ 0.5 and we get
Φ x( k) + Γ 0 u( k
x( k + 1)
where
e−1
1.5. In this case we
− 1 ) + Γ u( k − 2 )
1 − e−0.5
e−0.5 − e−1
A state representation is obtained from (2.12)
 x( k + 1) 
 u( k 1 ) 
u( k )
Γ 0   x( k)   0 
1   u( k 2 )  +  0  u( k )
u( k 1 )
  x( k) 
 1 0 0   u( k 2 ) 
u( k 1 )
y( k)
The pulse transfer function is
H ( z)
z2 ( z Φ)
z2 ( z 0.37)
Some calculations give
h( k )
Γ 0 e−( k−2) + Γ 1 e−( k−3) k ≥ 3
Using Theorem 2.4 and the definition of the z-transform (2.27) we can also
get the pulse response by expanding H ( z) in z−1 , i.e.,
H ( z)
The pulse response is now given by the coefficients in the series expansion.
There are three poles p1
0.37 and one zero z1
0 and p3
Problem 2.14
Sampling the given differential equation using the same procedure as in Problem 2.13 gives
a e−α h
Using the definition of τ ′ and with h
1 it follows that d
Further
4 and
e−α s b ds
− e−
α ( h τ ′)
e−α ( h−τ )
e−α s b ds
e−α ( h−τ
− e−
Straightforward calculations give
e−α ( h−τ
ab3 + b4
and
ab3 + b4
− ln1a ln
where it has been used that
ab3 + b4
− ln a
Problem 2.15
− 0.5 y(k − 1) u(k − 9) + 0.2u(k − 10) ⇔
y( k + 10) − 0.5 y( k + 9) u( k + 1) + 0.2u( k)
( q − 0.5q ) y( k) ( q + 0.2)u( k)
 A( q) q − 0.5q
y( k)
Pole excess
1 + 0.2q−1 u( k
y( k)
System order
B ( q)
deg A( q)
A∗ q−1 y( k)
− deg B (q)
deg A( q)
B ∗ q−1 u( k
d (cf. CCS p. 49).
Remark
q−10( q + 0.2)
B ( q)
A( q)
Problem 2.16
FIR filter:
y( k)
+ b n q − n u( k )
b 0 u( k ) + b 1 u( k
⇒ y( k + n)
+ b n u( k
H ( q)
Observable canonical form:
 a1 1 0
 a2 0 1
an 0
Problem 2.17
+ b n u( k )
+ b n u( k )
n:th order system
b 0 u( k + n) + b 1 u( k + n
⇒ q n y( k)
⇒ H ( q)
B ( q)
A( q)
− 1.5 y(k + 1) + 0.5 y(k)
q y( k) − 1.5qy( k) + 0.5 y( k)
u( k + 1 ) ⇔
y( k + 2)
qu( k)
Use the z-transform to find the output sequence when
y(0)
u( k )
y( 1)
Table 2.2 (page 57):
F1( z)
f ( jh) z− j
− y(0) − y(1)z− − 1.5z Y − y(0)
− u( 0 )
Compute
y(1)
− 0.5 y(−1) + 1.5 y(0) 1 − 0.5 + 0.75 1.25
⇒ z Y − 0.5 − 1.25z− − 1.5z( Y − 0.5) + 0.5Y z( U − 1)
U ( z)
0.5z( z − 1)
U ( z)
( z − 1)( z − 0.5) ( z − 1)( z − 0.5)
(step) ⇒
z − 0.5 ( z − 1) ( z − 0.5)
z − 0.5 ( z − 1)
u( 0 )
Y ( z)
Y ( z)
Table 2.3 (page 59) gives via inverse transformation.
e−1/T
e−k/T
e−1/T
e−( k−1)⋅ln2
(z
⇒ y( k)
0.5e−k⋅ln2
+ 2( k − 1) + e− −
(z
 y(1)
y(0)
y( 1)
(h
( k 1) ln2
Checking the result:
y(2)
0.5e−2 ln2 + 2 + e− ln 2
y(2)
u(1) + 1.5 y(1)
− 0.5 y(0)
Problem 2.18:
Verify that
h2 z( z + 1 )
2( z 1)3
( kh)2
F ( z)
(cf. Table 2.3)
( kh)2 z−k
differentiate twice
⇒(multiply by z)
f ( z−1)
( k2 + k) z−k−2
( z 1)3
( z 1)3
− − ( z − 1)
(z
multiplication by z2 gives ⇒ Σ( k2 + k) z−k
⇒ F ( z)
( z 1) z
( z 1)2
( z 1)2
− k ( − k − 1 ) z− −
2( z − 1)
( z − 1)
( z − 1)
f ( z−1)
z( z + 1 )
( z 1)3
h2 z( z + 1 )
2 ( z 1)3
Remark A necessary condition for convergence of f ( z−1) is that f is analytic
for z > 1. In that case there exists a power series expansion and termwise
differentiation is allowed.
Double integrator: The first step is to translate G ( s) to the corresponding pulse
transfer operator H ( z). Use the method of page 58.
The sampled system.
u( k )
Y ( s)
G ( s) ⋅
step ⇒ U ( z)
⇒ Y ( z)
h2 z( z + 1 )
2 ( z 1)3
Y ( z)
Y ( z)
U ( z)
cf. example A1: H ( z)
C ( zI
h2 z( z + 1 ) z 1
2 ( z 1)3
( Table 3.3)
H ( z) U ( z) ⇒
H ( z)
− 1 ⇒ U (s)
h2 ( z + 1 )
2 ( z 1)2
Problem 2.19
a.
The transfer function of the continuous time system is
( s + a)
G ( s)
Equation (2.30) gives
Ress
e−ah
H ( z)
−a
esh
e−ah 1
−ah
e
a
a z
− e−
− e−
s+a
ah
ah
This is the same result as obtained by using Table 2.1.
The normalized motor has the transfer function
s( s + 1)
G ( s)
and we get
esh 1
Res
z e
s( s + 1)
H ( z)
H ( z)
( z − 1)( e− − 1) + h( z − e− )
( z − e− )( z − 1)
( e− − 1 + h) z + (1 − e− − he− )
( z − e− )( z − 1)
e−h 1
z e−h ( 1)2
Compare Problem 2.11.
Problem 2.21
Consider the discrete time system
is lead if β < α
z+ a
arg
eiω h + b
eiω h + a
arctan
arg
cos ω h + a + i sin ω h
− arctan
a + cos ω h
Phase lead if
arctan
> arctan
a + cos ω h
a + cos ω h
We thus get phase lead if b < a.
0 we can use the series expansion of ( esh 1)/s
To compute the residue for s
and get
Time
Figure 2.2 Simulation of the step response of the system in Problem 2.22 for b −0.9
(upper solid), −0.75 (upper dashed), −0.50 (dash-dotted), 0 (dotted), 0.5 (lower solid),
and 1 (lower dashed).
Problem 2.22
A state space representation of the system is
x( k + 1) 
y( k)
1 b  x( k)
−0.9 −0.75 −0.5 0 0.5 and 1 is shown in Fig. 2.2.
Simulation of the system for b
Problem 2.23
A state space representation of G ( s) is
−ax + (b − a)u
The assumption that the system is stable implies that a > 0. Sampling this
system gives
b a
u( kh)
x( kh + h) e−ah x( kh) + (1 e−ah)
a
y( kh)
x( kh) + u( kh)
The pulse transfer operator is
H ( q)
− a)(1 − e− )/a + 1
q − e−
q − e−
b− a
1 − e−
− e−
(b
ah
ah
ah
where
ah
a
a
ah
− e− − 1
ah
The inverse is stable if
a
− e− − 1 < 1
(1
a
Since a > 0 and (1
− e−
ah
− e−
ah
) > 0 then
ah
2a
1 e−ah
For b > 0 we have the cases
b ≤ 2a
The inverse is stable independent of h.
b > 2a
Stable inverse if ah < ln( b/( b
Problem 2.28
y( k)
y( k
− 1) + y(k − 2)
y(0)
− 2a)).
y(1)
The equation has the characteristic equation
which has the solution
y( k)
Using the initial values give
The solution has the form
√5 (√5 + 1)
√ (√5 − 1)
Problem 2.29
The system is given by
2q−10 + q−11 u( k)
y( k)
Multiply by q11. This gives
The system is of order 11 (
(2q + 1) u( k)
deg A( q)) and has the poles
(multiplicity 9)
and the zero
y( k)
Problem 2.30
The system H1 has a pole on the positive real axis and can be obtained by sampling
a system with the transfer function
s+a
G ( s)
H2 has a single pole on the negative real axis and has no continuous time equivalent.
The only way to get two poles on the negative real axis as in H3 is when the
continuous time system have complex conjugate poles given by
π /ω where
Further we must sample such that h
We have two possible cases
G1 ( s)
G2 ( s)
Sampling G1 with h
π /ω gives, (see Table 2.1)
(1 + α )( q + α )
( q + α )2
H ( z)
where
e−ζ ω 0 h
i.e., we get a pole zero concellation. Sampling G2 gives H ( z)
0. This implies
that H3 cannot be obtained by sampling a continuous time system. H4 can be
rewritten as
H4 ( q ) 2 +
q( q 0.8)
The second part can be obtained by sampling a first order system with a time
delay. Compare Example 2.8.
Problem 2.31
We can rewrite G ( s) as
G ( s)
Using Table 2.1 gives
H ( q)
0.02 we get
H ( q)
(q
− e−
− e−
− e−
− e−
− 0.9802)(q − 0.9418)
Problem 2.32
 x +   u( t
eh−τ (1
 e 
se
 e
he
eh
e
−1 + e − 
1 + ( h − τ − 1) e −
e −
0   −1 + e
(h − τ )e − e −
1 + (τ − 1) e
eh−τ ( 1 + eτ )
e − −1
( h − τ ) eh−τ − eh−τ + 1
 e 
se 0
 e 
se
e A( h−τ )
( s 1)2
− h + τ ) + ( h − 1) e
The pulse transfer operator is given by
H ( q)
C ( qI
− Φ)− (Γ
where
Solutions to Chapter 3
Problem 3.1
Jury’s criterion is used to test if the roots are inside the unit circle.
a.
The roots are inside the unit circle since the underlined elements are positive.
(The roots are 0.75 ± 0.58i.)
One of the underlined elements is negative and there is thus one root outside the
unit circle. (The roots are 2.19 and 0.40 ± 0.25i.)
There are two roots outside the unit circle. (The roots are 0.35 and 0.82 ± 0.86i.)
One root is outside the unit circle. (The roots are
e.
−5 −0.5, and 0.5.)
The table breaks down since we can not compute α 1 . There is one of the underlined
elements that is zero which indicates that there is at least one root on the stability
boundary. (The roots are 0.7 and 0.5 ± 0.866i. The complex conjugate roots are on
the unit circle.)
Problem 3.2
The characteristic equation of the closed loop system is
z( z
− 0.2)(z − 0.4) + K
The stability can be determined using root locus. The starting points are z 0 0.2
and 0.4. The asymptotes have the directions ±π /3 and π . The crossing point
of the asymptotes is 0.2. To find where the roots will cross the unit circle let
z a + ib, where a2 + b2 1. Then
( a + ib)( a + ib
− 0.2)(a + ib − 0.4) − K
− ib and use a + b 1.
a − 0.6a − b + 0.08 + i(2ab − 0.6b) − K ( a − ib)
Multiply with a
Equate real and imaginary parts.
a2
− 0.6a − b
b(2a
−K a
− 0.6a − 1 − a + 0.08 −a(2a − 0.6)
4a − 1.2a − 0.92 0
The solution is
a
This gives K 2a − 0.6 0.70 and −1.30. The root locus may also cross the unit
circle for b 0, i.e. a ±1. A root at z −1 is obtained when
−1(−1 − 0.2)(−1 − 0.4) + K 0
0 then
a2
Figure 3.1
Re
The root locus for the system in Problem 3.2.
There is a root at z
1 when
1(1
− 0.2)(1 − 0.4) + K
The closed loop system is thus stable for
The root locus for K > 0 is shown in Fig. 3.1.
Problem 3.3
Sampling the system G ( s) gives the pulse transfer operator
H ( q)
The regulator is
u( kh)
where K > 0 and τ
a.
When τ
K uc ( kh
− τ ) − y(kh − τ )
K e( kh
0 then regulator is
u( kh)
K e( kh)
and the characteristic equation of the closed loop system becomes
The system is thus stable if
When there is a delay of one sample (τ
h) then the regulator is
e( kh)
u( kh)
and the characteristic equation is
K h + z( z
The constant term is the product of the roots and it will be unity when the
poles are on the unit circle. The system is thus stable if
Consider the continuous system G ( s) in series with a time delay of τ seconds.
The transfer function for the open loop system is
e
Go ( s)
The phase function is
arg Go ( iω )
and the gain is
Go ( iω )
The system is stable if the phase lag is less than π at the cross over frequency
Go ( iω c )
The system is thus stable if
The continuous time system will be stable for all values of K if τ
0 and
h. This value is about 50% larger than the value
for K < π /(2h) when τ
obtained for the sampled system in b.
Problem 3.4
The Nyquist curve is Ho ( eiω ) for ω
Ho ( eiω )
[0 π ]. In this case
eiω
Re
Figure 3.2
0 then Ho (1)
The Nyquist curve for the system in Problem 3.4.
2 and for ω
arg Ho
π then Ho ( 1)
−arctg
The argument is π /2 for ω
π /3(cos ω
values for the real and imaginary parts.
−2/3. The argument is
− 0.5). The following table gives some
ReHo
The Nyquist curve is shown in Fig. 3.2.
Problem 3.5
Consider the system
 x( k) +   u( k)
 0 1  x( k)
x( k + 1)
y( k)
We have
y(1)
x2 (1)
y(2)
x2 (2)
x1 (1) + x2 (1)
x1 (1)
Further
x(2)
x(3)
   +   ⋅ ( 1)  
Φ x(1) + Γ u(1)
The possibility to determine x(1) from y(1), y(2) and u(1) is called observability.
Problem 3.6
a.
The observability matrix is
The system is not observable since det Wo
The controllability matrix is
det Wc
2 and the system is reachable.
Problem 3.7
The controllability matrix is
For instance, the first two colums are linearly independent. This implies that Wc
has rank 2 and the system is thus reachable. From the input u′ we get the system
x( k + 1) 
 x( k) +   u′ ( k)
In this case
rankWc
1 and the system is not reachable from u′ .
Problem 3.8
a.
x( k + 1)
x(1)
x(2)
u( 0 )
u( 1 )
 x( k) +  1  u( k)
  1  +  u( 0 )   3 + u( 0 ) 
  0   3 + u( 0 ) 
3   3 + u( 0 )  +  u( 1 )   u( 1 ) 
⇒ x(2)
Two steps, in general it would take 3 steps since it is a third order system.
not full rank ⇒ not reachable, but may be controlled.
 1 1 1  is not in the column space of Wc and can therefore not be reached
from the origin. It is easily seen from the state space description that x3 will be
0 for all k > 0.
Problem 3.11
The closed loop system is
y( k)
a.
Hcl ( q)uc( k)
u c ( k)
K where K > 0 we get
y( k)
− 0.5q + K u (k)
Using the conditions in Example 3.2 we find that the closed loop system is
stable if K < 1. The steady state gain is
Hcl (1)
With an integral controller we get
Hcl ( q)
q( q
1)( q 0.5) + K q
q( q2
The system is stable if
and
and we get the condition
The steady state gain is Hcl (1)
Problem 3.12
The z-transform for a ramp is given in Table 2.3, also compare Example 3.13, and
we get
U c ( z)
( z 1)2
Using the pulse transfer functions from Problem 3.11 and the final value theorem
gives
lim e( k)
− y(k)
( 1 − H ( z)
lim uc ( k)
U c ( z)
if K is chosen such that the closed loop system is stable.
Time
Figure 3.3 The continuous time (solid) and sampled step (dots) responses for the system
in Problem 3.13.
a.
To use the final value theorem in Table 2.3 it is required that (1 z−1 ) F ( z)
does not have any roots on or outside the unit circle. For the proportional
controller we must then look at the derivative of the error, i.e.
lim e′ ( k)
z − 0.5z + K ( z − 1)
The derivative is positive in steady-state and the reference signal and the
output are thus diverging.
For the integral controller we get
lim e( k)
− 1 z(z − 1.5z + 0.5 + K − K ) z
z( z − 1.5z + 0.5 + K ) ( z − 1)
Problem 3.13
Consider the system
G ( s)
0.1 and the undamped natural frequency is
The damping of the system is ζ
1. The step response of the system will have an oscillation with the frequency
− ζ √0.99 rad/s
The sampled system will not detect this oscillation if h k2π /ω . The frequency
Fig. 3.3 shows the continuous time and the sampled step responses when h
6.3148. To formalize the analysis we can sample the system with h
2π /ω . The pulse transfer function is (Table 2.1)
ω will then have the alias ω ′
H ( q)
(1
− α )(q − α )
(q − α )
where α
exp( 0.1h). There is a pole zero cancellation and only the first order
exponential mode is seen at the sampling points.
Re
Figure 3.4 The root locus when the tank system in Problem 3.14 is controlled with an
integral controller.
Problem 3.14
a.
When Hr
K then the pulse transfer function for the closed loop system is
H c ( q)
K (0.030q + 0.026)
1.65q + 0.68 + K (0.030q + 0.026)
The characteristic equation is
z2 + (0.030K
The closed loop system is stable if
This gives the stability region
The steady state gain is
H c (1)
The steady state error when the input is a unit step is
e(
− H (1)
The minimum error is about 4%.
If the closed loop system is stable then the steady state error is zero if the
integral controller is used. The characteristic equation is then
( z2
− 1.65z + 0.68)(z − 1) + K z(0.03z + 0.026)
A scetch of the root locus is shown in Fig. 3.4.
Time
Figure 3.5 Step response of the system in Problem 3.14 when the proportional controller
is used with K 0.5 (solid), 1 (dashed), and 2 (dash-dotted).
Time
Figure 3.6 Step response of the system in Problem 3.14 when the integral controller is
used with K 0.01 (solid), 0.025 (dashed), and 0.05 (dash-dotted).
Fig. 3.5 and 3.6 show the step responses when the proportional controller
with K 0.5, 1, and 2 and the integral controller with K
0.01, 0.025, and
0.05 are used on the tank system.
Problem 3.16
The pulse transfer function for the closed loop system is
H c ( z)
H o ( z)
1 + H o ( z)
and the characteristic equation is
z2 + z + 0.16 + K (4z + 1)
z2 + (1 + 4K ) z + 0.16 + K
Using the conditions in Example 3.2 give
This implies
Since it is assumed that K > 0 we get the condition
for stability.
Problem 3.17
Using (3.17) we find that the transformation is given by
where Wc is the controllability matrix of the original system and Wc is the controllability of the transformed system.
The pulse transfer function of the system is
H ( z)
Thus the transformed system is given by
and
The transformation to give the controllable form is thus
In the same way the tranformation to the observable form is given by
Probem 3.18
a)
e)
Poles are mapped as z esh . This mapping maps the left half
plane on the unit circle.
see a)
When a system with relative degree > 1, “new” zeros appear
as described on p. 63 CCS.
These zeros may very well be outside the unit circle.
See Example 2.18 p. 65 CCC
Sample the harmonic oscillator, Example A.3 p. 530 CCS.
 1 0  not controllable
See e)
See p. 63, CCS
Problem 3.19
a.
The controllability matrix is
Since one column is zero we find that the system is not reachable (det Wc 0),
see Theorem 3.7.
The system may be controllable if the matrix Φ is such that Φ n x(0) 0. In
this case
and the origin can be reached from any initial condition x(0) by the control
sequence u(0) u(1) u(2) 0.
Probem 3.20
a.
The system is stable if, see p. 82,
e( k)
E ( z)
tot ( z)
U c ( z)
0.5) the
If K is chosen such that the closed loop system is stable (e.g. K
final-value theorem can be used.
lim e( k)
− 1 E ( z)
Re
Figure 3.7
The root locus for Problem 3.21b.
Problem 3.21
y( k)
a.
u( k )
(0.4q + 0.8) K
− 1.2q + 0.5 + K (0.4q + 0.8)
(0.4q + 0.8) K
q2 + q(0.4K 1.2) + 0.5 + 0.8K
The system is stable if, see p. 82 CCS,
q( q2
(0.4q + 0.8) K
1.2q + 0.5) + K (0.4q + 0.8)
(0.4q + 0.8) K
+ q(0.5 + 0.4K ) + 0.8K
Using root locus, Fig. 3.7, we can determine that the closed loop system is
stable if
Solutions to Chapter 4
Problem 4.1
The characteristic equation of the closed loop system is
det ( zI
− (Φ − Γ L))
( a11 + a22 b2 2 b1 1 ) z + a11 a22
+ ( a12 b2 a22 b1 ) 1 + ( a21 b1 a11 b2 ) 2
−a
12 a21
Identifying with the desired characteristic equation gives
a12 b2 a22 b1 a21 b1 a11 b2
p2 det Φ
where tr Φ
a11 + a22 and det Φ a11 a22 a12 a21 . The solution is
1  a21 b1 a11 b2
a12 b2 + a22 b1 b1
p2 det Φ
where
To check when ∆
a21 b2
−a
+ b1 b2 ( a22
−a
0 we form the controllability matrix of the system
 b a b + a b 
b2 a21 b1 + a22 b2
and we find that
det Wc
There exists a solution to the system of equations above if the system is control1 and dead beat response we have
lable. For the double integrator with h
a11 a12 a22 b2 1, a21 0, b1 0.5, and p1 p2 0.
This gives
  ( 1) 
This is the same as given in Example 4.4.
Problem 4.2
In this case the desired characteristic equation is
(z
− 0.1)(z − 0.25)
Using the result from Problem 4.1 we find that
and L is obtained from
To check the result we form
The matrix Φ − Γ L thus has the desired eigenvalues 0.25 and 0.1.
Output u(0)
Sampling period h
u(0) s function of h in Problem 4.3.
Figure 4.1
Problem 4.3
For the motor in Example A.2 we have
where α
e−h. The dead beat controller is characterized by
From Problem 4.1 it follows that
1  (1 α )2 α ( h 1 + α ) 1 h α   1 + α 
1  (1 α )2 (1 + α ) + α 2 (1 h α ) 
where
If x(0)
∆ (1 α )3 + (1
 1 1  then
u( 0 )
− α ) (h − 1 + α )
h( 1
(
This function is shown in Fig. 4.1. We thus want to find h such that
h( 1 − α )
2(1
− e− )
− 2e− + 1
2(1
2e−2h
Time
Figure 4.2
Deadbeat control of the motor when xT (0)
 1 1  and h
This is a nonlinear equation of the form
f ( h)
One way to find h is to use the iterative scheme
f ( h k)
Starting with h0
2 gives
 0.49 0.51 . Fig. 4.2 shows the response of the
motor controlled with a deadbeat controller when h
2.21. We see that the
system settles after two samples and that the constraint on the control signal is
fulfilled.
2.21 we get L
Problem 4.4
a.
Using the results in Problem 4.1 gives
The problem is easily solved using Matlab by giving the commands
ga=[0.01; 0.16];
P=roots([1 -0.63 0.21])
L=place(fi,ga,P)
Time
Figure 4.3
x( 0 )
The response and the control signal of the system in Problem 4.4 when
 1 0  and L  9.22 3.11 .
From Example 4.4 we find that the continuous-time characteristic polynomial
s2 + 2ζ ω s + ω 2 corresponds to z2 + p1 z + p2 with
−2e−
e−2ζ ω h
cos(ω h
This has the solution ζ
0.7 and ω
5.6, so the characteristic polynomial
becomes s2 + 7.8s + 31.7. In Matlab you can do
rd=roots([1-0.63 0.21])
rc=Log(rd)/h
Ac=poly(rc)
The chosen sampling interval is higher than recommended by the rule of
thumb, since
The closed loop system when using L  9.22 3.11  is shown in Fig. 4.3
when x(0)  1 0  .
Problem 4.5
a.
In this case
We get
−  y(k 1)  + Ψ u(k 1)
y( k)
 3.55 3.55   y( k 1)   0.114 
 u( k
y( k)
x( k)
The dynamical observer (4.28) has the form
− K C )xˆ(k k − 1) + Γu(k) + K y(k).
In this case we choose K such that the eigenvalues of Φ − K C are in the
x ( k + 1 k)
(Φ
origin. Using the results from Problem 4.1 but with Φ T and C T instead of Φ
and Γ gives
The reduced order observer (4.32) has the form
x ( k k)
− K C )(Φ xˆ(k − 1 k − 1) + Γu(k − 1)) + K y(k).
(I
In this case we want to find K such that
( I K C )Φ has eigenvalues in the origin. The first condition implies that
k2 1. Further
( I K C )Φ 
The eigenvalues will be in the origin if
The observer is then
 x( k
x ( k k) 
Since x2 ( k k)
x ( k k)
− 1 k − 1) +  0.114  u(k − 1) +  3.1  y(k)
y( k) we get
 3.55 3.55   y( k 1)   0.114 
 u( k
y( k)
which is the same as the observer obtained by direct calculation.
Problem 4.6
From Problem 2.10 we get for h 12
 x( kh) + 
 u( kh)
x( kh + h) 
y( kh)  0 1  x( kh)
The continuous time poles of the system are 0.0197 and 0.0129. The observer
should be twice as fast as the fastest mode of the open loop system. We thus
choose the poles of the observer in
e−0.0394⋅12
x1 and xe1
x2 and xe2
Time
Figure 4.4
The states (solid) and their estimates (dots) for the tank system in Problem 4.6
The desired characteristic equation of Φ
Using the results from Problem 4.1 gives
Fig. 4.4 shows the states and the estimated states when x(0)
when u( kh) is zero up to t
 1 1  and
250 and one thereafter.
Problem 4.7
The observer and the controller are described by
x( k k)
( I K C )Φ x( k
L x ( k k).
u( k )
− 1 k − 1) + (I − K C )Γu(k − 1) + K y(k)
In the state equation both x and y have the time argument k. Introduce
ξ ( k)
then
ξ ( k)
x( k k)
− K y(k)
− K C )Φ[ξ (k − 1) + K y(k − 1)] + (I − K C )Γu(k − 1)
− K C )Φξ (k − 1) + (I − K C )Φ K y(k − 1)
− (I − K C )Γ L[ξ (k − 1) + K y(k − 1)]
( I − K C )(Φ − Γ L)ξ ( k − 1) + ( I − K C )(Φ − Γ L) K y( k − 1)
Φ ξ ( k − 1) + Γ y( k − 1)
(I
(I
The output of the regulator can be written
u( k )
− Lξ (k) − LK y(k)
Coξ ( k) + Do y( k).
The observer and the regulator can thus be written in the form given in the
formulation of the problem.
Problem 4.8
The constant disturbance v( k) can be described by the dynamical system
w( k + 1)
w( k)
v( k)
w( k)
The process can thus be described on the form given in (4.43) with
a.
If the state and v can be measured then we can use the controller
u( k )
− Lx(k) − L w(k).
This gives the closed loop system
Φ x( k) + Φ xw w( k) Γ Lx( k) Γ Lw w( k)
(Φ Γ L) x( k) + (Φ xw Γ Lw )w( k)
C x( k)
x( k + 1)
y( k)
In general it is not possible to totally eliminate the influence of w( k). This is
only possible if Φ xw Γ Lw is the zero matrix. We will therefore only consider
the situation at the output in steady state
y(
− (Φ − Γ L)]− (Φ − Γ L )w(∞)
Hw (1)w(
The influence of w (or v) can be zero in steady state if
Hw ( 1 )
This will be the case if
is the ( i j ) element of Φ − Γ L and γ
γ 1 (1
where ϕ cij
i is the i:th element of Γ .
Assume that L is determined to give a dead beat regulator then
and
and
In this case is the state but not the disturbance measurable. The disturbance
can now be calculated from the state equation
Φ xw w( k
x( k)
− Φ x(k − 1) − Γu(k − 1).
The first element in this vector equation gives
w( k
[1 0]( x( k)
− Φ x(k − 1) − Γu(k − 1))
Since w( k) is constant and x( k) is measurable it is possible to calculate w( k)
w( k 1). The following control law can now be used
u( k )
− Lx(k) − L w(k)
where Lw is the same as in (a). Compared with the controller in (a) there is
a delay in the detection of the disturbance.
Estimate ve
Time
Time
Figure 4.5 The output of the system in Problem 4.8 for the regulators in a) (upper left),
b) (upper right) and c) (lower left and right). The estimate of v is also shown for case c).
Notice the difference in scale in the upper left curve.
If only the output is measurable then the state and the disturbance can be
estimated using an observer of the form (4.41)
 x ( k + 1) 
w( k + 1)
ε ( k)
 Φ Φ xw   x( k)   Γ 
 ε ( k)
 +   u( k ) + 
w( k)
y( k) C x( k)
The gain vector can now be determined such that the error goes to zero provided the augmented system is observable. The error equation is
 x( k + 1) 
w( k + 1)
Φ xw   x( k) 
w( k)
The characteristic equation of the system matrix for the error is
z3 + ( k 1
+ (1.05
The eigenvalues are in the origin if
The controller now has to be
u( k )
− Lxˆ (k) − L w(k)
where L and Lv are the same as in (a). The solutions above have the drawback that there may be an error in the output due to the disturbance if there
are small errors in the model. Fig. 4.5 show that the output when the controllers in (a), (b) and (c) are used.
Problem 4.9
a.
The state equation for the tank system when h 12 was given in the solution
to Problem 4.6. The desired characteristic equation is
Using the result in Problem 9.1 give
An integrator can be incorporated as shown in Section 4.5 by augmenting the
system with
x3 ( kh + h) x3 ( kh) + uc ( kh) C x( kh)
and using the control law
− Lx(kh) −
u( kh)
3 x3 ( kh) + c u c ( k)
The closed loop system is then
 x( k + 1) 
x3 ( k + 1)
+  0.0296 c  uc ( k)
  x( k) 
 x3 ( k) 
The characteristic equation is
+ (2.3240 − 0.5218 − 0.0035 − 0.0296 ) z+
+ (−0.6770 + 0.2408 − 0.0261 − 0.0261 )
z3 + ( 2.647 + 0.281
Assume that two poles are placed in the desired location and that the third
pole is in p. We get the following system of equations to determine the state
feedback vector.
This gives
The parameter c will not influence the characteristic equation, but it is a
feedforward term from the reference signal, see Fig. 4.11 CCS. Fig. 4.6 shows
the step response for some values of p. Fig. 4.7 shows the influence of c when
Time
Figure 4.6 The stepresponse for the tank process when p
and when c 0.
0 (solid) and 0.5 (dashed),
Time
Figure 4.7 The step response for the tank process when p
3 (dashed) and 6 (dash-dotted).
0 and when
0 (solid),
Problem 4.10
The process is
x( kh + h)
 x( kh) + 
 u( k ) + 
 v( k)
where v( k) is a sinusoidal. I.e. it can be described by
w( kh + h)
v( k)
 w( kh)
 1 0  w( kh)
The augmented system (9.33) is now
 x( kh + h) 
w( kh + h)
where α
cos ω o h and β
  x( kh)   h 
 0  u( kh)
β  w( kh)
sin ω o h. Assume first that the control law is
− Lx(kh) − L w(kh)
u( kh)
where
 1 0  we would
totally eliminate v since v and u are influencing the system in the same way.
Compare the discussion in the solution of Problem 4.8.
Since x( kh) and ξ ( kh) cannot be measured we use the observer of the structure
(9.35) where
i.e., a dead beat controller for the states. Further if Lw
and
The error equation is then
 x( k + 1) 
w( k + 1)
− K C Φ   x˜(k) 
− K C 1   w(k) 
0   x( k) 
−β  w(k) 
Let the desired characteristic equation of the error be
(z
1 and ω o
0.1π then for γ
which has the solution
0.5 we get the following system of equations
Fig. 9.8 shows the states of the double integrator and their estimates when v( t)
sin(ω o t). It is seen that the controller is able to eliminate the disturbance.
Problem 4.11
We have to determine the feedback vector L such that
x1 and xe1
x2 and xe2
Time
The states of the double integrator 
(solid) and their estimates (dots) when
Figure 4.8
v( t)
 1 1.5  and Lw
sin(ω o t). The controller is defined by L
a.
has all eigenvalues in the origin. This gives the condition
det(λ I
λ 2 + ( 1.6 +
I.e.,
The stationary gain of the closed loop system is given by stationary gain of
the
m ⋅ C ( I Φ + Γ L)−1 Γ m
To get unit steady state gain we choose m
The closed loop characteristic equation is stable if (See Example 3.2)
This gives
−1 + (−1.6 +
> −1 − (−1.6 +
Deadbeat
Figure 4.9
The stability area for L in Problem 4.11.
The stability area is shown in Fig. 4.9. Assume that the deadbeat control
in a. is used. The closed loop system will be unstable if the feedback from
0), but the system will remain stable if x2 is
x1 is disconnected (i.e., if 1
disconnected (if 2 0).
Problem 4.12
a.
The deadbeat requirement implies that the characteristic equation of the
closed loop system is
det(λ I Φ + Γ L) det 
λ2 + λ(
There are infinitely many solutions, one is
The controllability matrix is
det Wc
0 implies that the system is not reachable and arbitrary states
 2 8  is in the
cannot be reached from the origin. However, x( k)
column space of Wc and the point can thus be reached.
Φ 2 x(0) + ΦΓ u(0) + Γ u(1)
 x(0) + 
 u( 0 ) +   u( 1 )
 2 8  and x(0)  0 0  we get
x(2)
With x(2)
2.25u(0) + u(1)
9u(0) + 4u(1)
One solution is
u( 1 )
u( 0 )
The observer should have the characteristic equation
(λ
det(λ I
det 
λ + ( k − 2.25)λ + 0.5k − 2k
Identifying coefficients give the system
which has the solution
Solutions to Chapter 5
Problem 5.1
Euclid’s algorithm defines a sequence of polynomials A0
B and
An where A0
If the algorithm terminates with An+1 then the greatest common factor is An . For
the polynomials in the problem we get
This gives
−0.4z + 0.42z − 0.11 −0.4(z − 1.05z + 0.275)
The next step of the algorithm gives
( z 0.95)/( 0.4)
Finally
(z
0.1775( z
The greatest common factor of A and B is thus z
Problem 5.2
a.
To use Algorithm 5.1 we must know the pulse-transfer function B ( z)/ A( z)
and the desired closed-loop characteristic polynomial Acl ( z). Since the desired
pulse-transfer function from uc to y is Hm ( z) (1 + α )/( z + α ), we know at
least that ( z + α ) must be a factor in Acl ( z).
Step 1. You easily see that, with A and Acl being first order polynomials and
B a scalar, you can solve the equation for the closed loop system using scalars
R ( z) r0 and S ( z) s0 :
A( z) R ( z) + B ( z) S ( z)
( z + a) ⋅ r0 + 1 ⋅ s0
Acl ( z)
Identification of coefficients gives the following equation system:
ar0 + s0
with the solution r0
1 and s0
− a.
Step 2. Factor Acl ( z) as Ac ( z) Ao( z) where deg Ao ≤ deg R
and Ac z + α . Choose
T ( z)
Ac (1)
A o ( z)
B (1)
t0 Ao ( z)
With this choice of T , the static gain from uc to y is set to 1 ( Hm (1) 1), and
the observer dynamics are cancelled in the pulse-transfer function from uc to
y. In this case, there are no observer dynamics, though, since deg Ao 0.
The resulting control law becomes
R ( q ) u( k )
u( k )
− S(q) y(k)
(1 + α )u ( k) − (α − a) y( k)
T ( q)u c ( k)
i.e., a (static) proportional controller.
Solution with higher order observer: The solution above is not the only
one solving the original problem. We can, for example, decide to have another
closed loop pole in z
β , say.
Step 1. To solve the equation for the closed loop characteristic polynomial we
must increase the order of R by one. This gives
( z + α )( z + β )
( z + a)( r0 z + r1 ) + 1 ⋅ s0
and the equation system becomes
ar0 + r1 α + β
ar1 + s0 α β
Step 2. Splitting Acl into factors Ao
T ( z)
t0 Ao ( z)
Ac (1)
A o ( z)
B (1)
u( k )
T ( q)u c ( k)
− α + β − a u( k − 1 ) +
α +β a
α β a(α + β
z + β and Ac
The resulting control law becomes
R ( q ) u( k )
− a)
z + α gives
(1 + α )( z + β )
− S(q) y(k)
− α β − a(α + β − a)
u c ( k) + β u c ( k
y( k
The controller thus is a dynamical system. In this case there is a delay of
one sample from the measurements y to the control signal u. This could have
been avoided by choosing deg S 1.
The closed loop characteristic polynomial is given by AR + B S, i.e. ( z + α ) in
the first solution, and ( z + α )( z + β ) in the second one. In both cases we get
the same closed loop pulse-transfer function from uc to y since the observer
polynomial is cancelled by T ( z):
H m ( z)
B ( z) T ( z)
A( z) R ( z) + B ( z) S ( z)
t0 B ( z)
A c ( z)
Fig. 5.1 shows the response when the two different controllers are used. It is
0.99, and that the design parameters
assumed in the simulations that a
0.7 and β
0.5. You can see the effect of the observer polynomial
are α
when regulating a nonzero initial state, but not in the response to a set point
change.
Time
Figure 5.1 The output of the process with y(0) 1, and uc ( k) is 0 for k < 25 and −1
for k > 25. The dots corresponds to the zero order controller, and the crosses to the first
order controller.
Problem 5.3
a.
The desired closed loop pulse-transfer function is
H m ( z)
B m ( z)
In this case, the process zero is cancelled by the controller, so ( z + 0.7) is not a
factor of Bm . By choosing Bm 0.2z we make the steady state gain Hm (1)
1, and the pole excess in Hm equals the pole excess in H. Cancellation of
process poles and zeros is handled by Algorithm 5.3 or through the following
B + B − , where B + is the part
First, the process numerator is factored as B
of the numerator which should be cancelled by the controller, i.e., B + z + 0.7
and B −
1. B + must be a part of the R polynomial as well as Acl . This
gives the Diophantine equation
A( z) B + ( z) R( z) + B + ( z) B −( z) S ( z)
( z2
− 1.8z + 0.81)R(z) + S(z)
B + ( z) Acl ( z)
Acl ( z)
If R ( z) is a constant r0 , the left hand side is of second order, and so must Acl
¯ the causality condition (deg S ≤ deg R 1) leads
be. With this choice of R,
us to set S ( z) s0 z + s1 . Now, we can solve the Diophantine equation above,
since we have 3 indeterminates ( r0 , s0 and s1 ) and 3 coefficients to set:
( z2 1.8z + 0.81) ⋅ r0 + ( s0 z + s1 )
B + ( z) R ( z)
Thus, we have R ( z)
the desired
B + ( z2 1.5z + 0.7)
H m ( z)
z + 0.7 and S ( z)
− 0.11. To obtain
we must select
T ( z)
The controller is now
−0.7u(k − 1) + 0.2u (k) − 0.3 y(k) + 0.11 y(k − 1).
u( k )
In this case we do not want to cancel the process zero, so
H m ( z)
in order to get Hm (1)
given by the identity
( z2
1.7 ( z + 0.7)
B m ( z)
1. The closed loop characteristic polynomial is now
− 1.8z + 0.81)R(z) + (z + 0.7)S(z)
Acl ( z)
The simplest choice, a zero order controller, will not suffice in this case since
it would only give 2 parameters r0 and s0 to select the 3 parameters in the
second order polynomial z2 1.5z + 0.7. Thus, we must increase the order of
the controller by one and, consequently, add an observer pole which is placed
at the origin, i.e. Ao z and Acl z3 1.5z2 + 0.7z. Letting
the identity then gives the system of equations
Further T
u( k )
The controller is thus
−0.0875u(k − 1) + 0.1176u (k) − 0.2125 y(k) + 0.1012 y(k − 1)
Fig. 5.2 shows the output and the control signal for the controllers in Case a
and Case b. Case a should probably be avoided because of the ringing in the
control signal.
Problem 5.4
a.
Using the controller
u( k )
S ( q)
(u c ( k)
R ( q)
− y(k))
gives the closed loop system
Time
Time
Figure 5.2 The output and the control signal for the controllers in Case a (left) and
Case b (right) in Problem 5.3. The ringing in the control signal in Case a is due to the
cancellation of the process zero on the negative real axis.
Section 5.10 gives one solution to the problem
B ( Am
With the given system and model we get
1( z + α
( z + a)(1 + α )
The controller contains an integrator. Further the pole of the process is cancelled.
The characteristic polynomial of the closed loop system is
( z + α )( z + a)
The closed loop system will contain an unstable mode if a > 1. The controller
can be written
From this we can conclude that in order to get a stable closed loop we must
fulfill the following constraints.
Bm must contain the zeros of B that are outside the unit circle.
Bm must contain the poles of the process that are outside the unit
circle. The first constraint is the same as for the polynomial design discussed
in Chapter 5.
Problem 5.5
a.
Equation (5.33) gives the pulse transfer operator from uc and v to y:
y( k)
v( k)
u c ( k) +
The design in Problem 5.2 gave
−a
We thus get
If v( k) is a step there will thus be a steady state error 1/(1 + α ) in the output.
By inspection of the transfer function from v to y we see that we must make
R (1) 0 in order to remove the steady state error after a load disturbance
step. By forcing the factor ( z 1) into R ( z) we thus have obtained integral
action in the controller. The design problem is solved by using the general
Algorithm 5.3 or through a discussion like the one below.
With R ( z)
(z
− 1)R(z) the closed loop characteristic equation becomes
A( z)( z − 1) R( z) + B ( z) S ( z) A ( z)
( z + a)( z − 1) R( z) + 1 ⋅ S ( z) A ( z)
If R ( z) is a constant r0 , the left hand side is of second order, and so must Acl
be. With this choice of R, the causality condition (deg S ≤ deg R 1) leads
us to set S ( z) s0 z + s1 . Now, we can solve the Diophantine equation above,
since we have 3 indeterminates ( r0 , s0 and s1 ) and 3 coefficients to set:
(a
( z + a)( z
−ar
+ ( s0 z + s1 )
( z + α )( z + β ) ⇒
s0 α + β a + 1
s1 α β + a
To obtain the desired
H m ( z)
B ( z) T ( z)
A( z) R ( z) + B ( z) S ( z)
we must select
T ( z)
T ( z)
( z + α )( z + β )
(1 + α )( z + β )
The controller is now
u( k )
− 1) − (α + β − a + 1) y(k) − (a + α β ) y(k − 1)
+ (1 + α )u ( k) + β (1 + α )u ( k − 1)
u( k
Fig. 5.3 shows the controllers in Problem 5.2 and the controller with an
integrator. The reference value is zero and there is an initial value of the
state in the process. At t 25 a constant load disturbance is introduced. It is
assumed that a
0.99, and the design parameters are chosen as α
and β
Problem 5.6
B / A while the true
It is assumed that the design is based on the model H
model is H 0 B 0 / A0 . The pulse transfer operator of the closed loop system is
The design gives
and
Time
Figure 5.3 The output of the process in Problem 5.5 when the controller does not contain
an integrator (dots) and when an integrator is introduced (crosses).
This gives
Problem 5.7
a.
The design in Problem 5.2 gives
−a
Assume that the true process is
z + a0
Equation (5.41) gives
z + a0 + α
−a
Frequency
Figure 5.4 The left hand side of the inequality in Problem 5.7 when z eiω , 0 < ω < π
for a0
−0.955 (dashed), −0.9 (dash-dotted) and 0.6 (dotted). The right hand side of
the inequality is also shown (solid).
The closed loop system is stable if
a0 + α
−a <1
With the numerical values in the problem formulation we get
−1.4 < a
Equation (5.40) gives the inequality
H ( z)
− H ( z) <
H ( z)
H m ( z)
(1 + α )( z + a)
H f f ( z)
H f b ( z)
− 0.9 − z + a
α a
eiω the left hand side of the inequality for different
values of a . The right hand side is also shown.
Problem 5.8
Section 5.6 shows that the control signal is given by (5.52)
u( k )
H m ( q)
u c ( k)
H ( q)
(1 + α )( q + a)
u c ( k)
We may assume that both the process and the model have a continuous time
correspondence. This implies that a and α are less than zero. Further the desired
model is stable, i.e. α < 1. The control signal is now obtained by studying the
step response of Hm / H, which is a stable first order system. The largest value
is then either at the first step or the final value. The magnitude at the first step
can be determined either through the initial value theorem or by using series
expansion and the value is 1 + α . The final value is 1 + a. If α < a then
the closed loop system is faster than the open loop system and the control signal
is largest at the first step. If the desired response is slower than the open loop
system then the final value is the largest one.
Problem 5.14
a.
The rule of thumb on p. 130 gives
Identifying with
gives ω
0.1. Thus an appropriate sampling interval is
Using Example 2.16 we get sampled data characteristic equation
z2 + a1 z + a2
where
a1
a2
−2e−
e−2ζ ω h
cos(
The poles are in 0.66 ± 0.25i.
Problem 5.15
This solution demonstrates how to use Algorithm 5.3.
Data: The process is given by A
q2 1.6q + 0.65 and B
will at least contain Ac
q2 0.7q + 0.25, other factors may be added later
1 since no given factors are forced into the controller. The
desired response to command signals is assumed to be Hm
0.55/( q2 0.7q + 0.25) (cancelled process zero, Hm (1) 1).
Pole excess condition:
deg Am
− deg B ≥ deg A − deg B
Remark: The fact that we cancel one zero and do not introduce any other zero
in Bm causes the delay from the command signal to be one time unit more
than the delay of the process.
Model following condition:
Degree condition:
deg Acl
2 deg A + deg Am + deg R d + deg Sd
A+ B + Am Acl and Acl
Step 1. A+ 1, A−
A q2 1.6q + 0.65, B + q + 0.75 and B − 0.4 achives
cancellation of the process zero, but no cancellation of process poles.
Step 2. Using the degree condition above we may conclude that
− deg A − deg B − deg A − deg A
deg Ao
deg Acl
The Diophantine equation to solve thus becomes
(q
Since this is of second order, R must be a constant, r0 , say. In order to solve the
identity we must have two more parameters, so we let S s0 q + s1 :
( q2
− 1.6q + 0.65)r + (s q + s ) ⋅ 0.4
This gives the system of equations
Step 3. The controller polynomials are now given by (5.45):
Am ( q + 0.75)
Am (2.25q
Since, in this case, Am
Ac , this factor can (and should) of course be cancelled
in all controller polynomials, giving
The corresponding degree of the closed-loop polynomial AR + B S will thus be 3
instead of 5.
Problem 5.16
In this case we want to have an integrator in the controller, i.e., R d ( q 1). This
will increase the degree of the closed loop by one compared to Problem 5.15 (see
(5.42)), which is done by having Ao
( q + ao ), say. This gives the Diophantine
equation
(q
1.6q + 0.65)( q 1) R + 0.4 S ( q
0.7q + 0.25)( q + ao )
R must still be a constant (which as usual will be 1) and S must be of second
order:
( q2
− 1.6q + 0.65)(q − 1) + (s q
This gives
( q2
− 0.7q + 0.25)(q + a )
−2.6 + 0.4s a − 0.7
−0.7a + 0.25
−0.65 + 0.4s 0.25a
and
Using a0
S ( q)
2.5 ( a0 + 1.9) q2
− (0.7a
+ 2) q + 0.25a0 + 0.65
−0.25, (5.45) gives (after cancelling the common factor A
R B ( q − 1) q − 0.25q − 0.75
Problem 5.17
The minimum degree solution has deg A0 1 and gives a unique solution to the
Diophantine equation. Let us instead use deg A0 2 and deg S deg R 1. This
gives the equation
( z + 1)( z + 2)( z2 + r1 z + r2 ) + z ⋅ ( s1 z + s0 )
with the solution
The controller is causal. Using Theorem 5.1 we also have the solutions
− Q(z − 1)(z − 2)
where Q is an arbitrary polynomial. Choose for instance Q
7z + 6 + ( z
−1. This gives
This is also a causal controller. The closed loop systems when using R 0 S0 , T0
S0 and R S, T S respectively are
The number of zeros are different.
z(7z + 6)
z( z2 + 4z + 8)
Solutions to Chapter 6
Problem 6.1
In the first case it is assumed that we have a control structure as in Fig. 6.1.
There are three subsystems each with the transfer function
Gi ( s)
and the total transfer function from uc to y is
s( s + K 1 )( s + K 2 )( s + K 3 ) + K 1 K 2 K 3 K 4
If either of the gains K i is increased sufficiently much the closed system will
become unstable. Fig. 6.2 shows the response when uc is an impulse and when
K 1 K 2 K 3 1 and K 4 0.1, 0.25, and 0.75.
A disturbance in the process will propagate in the direction of the flow. In the
case of control in the direction opposite to the flow each of the subprocesses has
Raw
material
Final
Figure 6.1
Block diagram for the control in the direction of the flow in Problem 6.1.
Time
Figure 6.2 Impulse response for the control in the direction of the flow when K4
(solid), 0.25 (dashed), and 0.75 (dash-dotted).
Figure 6.3
Block diagram for the control in the direction opposite of the flow in Problem 6.1.
a transfer function of the type Gi . The system is then represented with the block
diagram in Fig. 6.3. Notice the order of the states. The system will remain stable
for all positive values of K i . A disturbance will now propagate in the direction
opposite the flow. A disturbance in uc will now only influence the first subprocess
and will not propagate along with the flow. The reader is strongly recommended to
compare with the case where the disturbance appears at the final product storage
instead.
Problem 6.2
Fig. 6.3 in CCS contains several examples of couplings of simple control loops.
a.
Cascade control loops are found for the cooling media flow and for the output
Feedforward is used for the level control loop where the input flow is used as
a measurable disturbance. The input flow is also used as feedforward for the
cooling of the jacket.
Nonlinear elements are used in the flow control loops of the product output
and the coolant flow. The flow is probably measured using differential pressure which is proportional to the square of the flow. The square root device
is thus used to remove the nonlinearity of the measurement device. An intentional nonlinearity is introduced in the selector. Either the temperature
or the pressure is used to control the coolant flow depending on the status of
the process.
Solutions to Chapter 7
Problem 7.1
Which frequencies will the signal
a1 sin 2π t + a2 sin 20t
f ( t)
give rise to when sampled with h 0.2?
Since sampling is a linear operation we consider each component of f ( t) separately.
The sampled signal has the Fourier transform, see (7.3)
F (ω + kω s )
where F (ω ) is the Fourier transform of the time continuous signal. The Fourier
transform of sin ω 0 t has its support (i.e., the set where it is 0) in the two points
±ω 0 . More precisely, it equals π i δ ω + ω 0
δ ω ω 0 . Thus, if the signal
sin ω 0 t is sampled with the sample interval h its Fourier transform will be 0 in
the points
2π and ω s
10π we get the angular frequencies
20 gives rise to
The output of the sampler is composed of the frequencies
Problem 7.2
We have the following specifications on the choice of sampling period and presampling filter:
All frequencies in the interval ( f 1 f 1 ) should be possible to reproduce from
the samples of the continuous time signal.
We want to eliminate the disturbance with the known and fixed frequency
The sampling theorem states that the first specification will be satisfied if and
only if the sample frequency f s is chosen such that
Moreover, for the disturbance f 2 not to fold on the data signal
( f s /2
Two cases:
a) h=2π/10
Figure 7.1
Folding for different frequencies in Problem 7.5.
Filter out the disturbance using an antialiasing filter. Then sample with f s >
2 f 1 . Suppose the disturbance should be attenuated 20 dB without effecting
the datasignal. A n:th order filter gives maximally n ⋅ 20 dB/decade. So to
achieve 20 dB in log f2
0.699 decades takes n 2.
If f s > 6 f 1 , the disturbance does not mix with the data signal. It can instead
be removed using digital filters.
Problems 7.5 and 7.6
The magnitude of the spectrum of the sampled signal can be obtained by folding
the spectrum of the time continuous signal around the angular frequency ω N
ω s /2 π /h. See Fig. 7.1 and Fig. 7.2.
Problem 7.7
The rotation frequency of the wheel ω r 2π r.
The frequency of the camera shutter ω s 2π /h.
The picture will not move if ω r n ⋅ ω s ; for integer values n.
A correct picture will be seen, if ω s > 2ω r according to the sampling theorem.
(The eye acts like a low pass filter).
The wheel will appear to rotate with a frequency lower than r if ω s < 2ω r . See
Fig. 7.3. For instance let ω s 4/3 ω r . Aliasing will give a frequency ω
The wheel then appears to rotate three times slower and in the wrong direction.
If ω s ω r the wheel will appear to stand still. Compare the stroboscope.
a) ω s =120
Figure 7.2
Folding for different frequencies in Problem 7.6.
Figure 7.3
Folding in Problem 7.7 when ω s < 2ω r .
Problem 7.9
The signal is
u( t)
[sin(6ω 0 t) + sin(2ω 0 t)]
sin(4ω 0 t) cos(2ω 0 t)
Sampling the signal with
gives the Nyquist frequency ω N
3ω 0 . Sampling the signal u( t) gives the alias
of sin(6ω 0 t) in ω 0. We thus get the frequencies
in the sampled signal.
Solutions to Chapter 8
Problem 8.1
The three transformations Euler’s method (forward difference (8.4)), backward
difference (8.5) and Tustin’s approximation (8.6) have different stability properties. This can be seen by finding how the left half s-plane is transformed into the
z-plane. For Euler’s method we have
This implies that the stability boundary in the sh-plane (the imaginary axis) is
translated one unit to the right, see Fig. 8.1a. When the backward difference is
used then
iω we get
This represents a circle with radius 0.5 and going through the points 0 and 1, see
Finally for Tustin’s approximation with s iω
arg z
2 arctan ω h
The imaginary axis is thus transformed into the unit circle in the z-plane. If a
transfer function is stable in the s-plane it will be translated into a stable discrete
time system if using the backward difference or Tustin’s approximation.
Problem 8.2
G ( s)
Forward differences
a
s+a
a>0
Backward differences
Figure 8.1 Transformation of the left half s-plane when using a. Eulers method, b.
Backward difference and c. Tustin’s approximation.
a.
Using Euler’s method we get
a
H ( z)
(z
ah
− 1)/h + a
− 1 + ah
This corresponds to the difference equation
− 1) y(kh)
y( kh + h) + ( ah
ahu( kh).
The difference equation is stable if
ah
0 < h < 2/a.
The approximation may, however, be poor even if the difference equation is
stable.
For Tustin’s approximation we get
H ( z)
( z + 1) ah/2
(1 + ah/2) z + ( ah/2
a
−1 +a
ah/2
1 + ah/2
ah/2
ah/2 + 1
The pole of the discrete time system will vary from 1 to 1 when h vary from
0 to infinity. The discrete time approximation is always stable when a > 0.
Using Tustin’s approximation with prewarping gives
a/α
1 + a/α
a
H ( z)
where
−1 +a
a/α
a/α + 1
a
tan( ah/2)
tan( ah/2)
1 + tan( ah/2)
H ( z)
tan( ah/2)
tan( ah/2) + 1
Problem 8.3
The lead network
Gk ( s)
should be approximated using different methods. At ω
argument 19○ and the gain 2.95.
a.
1.6 rad/s it has the
Euler’s method gives
H E ( z)
(z
(z
Backward differences
H B ( z)
(z
(z
− 1)/(zh) + 1
− 1)/(zh) + 2
z( 1 + h) 1
z(1 + 2h) 1
Tustin’s approximation
H T ( z)
− ( 1 − h/2 )
− ( 1 − h)
− (1 − 1/α )
− (1 − 2/α )
z( 1 + h/2 )
z( 1 + h)
Tustin’s approximation with prewarping
HTW ( z)
z(1 + 1/α )
z(1 + 2/α )
where
tan(ω 1 h/2)
Within two decimals this is the same as in (c).
e.
Zero order hold sampling gives
H Z O H ( z)
− e−
− 4 ⋅ 1 1 − e−
− e− − (1 − e−
z − e−
All five approximations have all the form
H ( z)
The gain and the phase at ω
H ( eiω h )
arg H ( eiω h )
H ( eiω h)
z+ a
1.6 are obtained from
( eiω h + a)( e−iω h + b)
eiω h + a
eiω h + b
( e + b)( e−iω h + b)
1 + ab + ( a + b) cos(ω h) + i( b a) sin(ω h)
1 + b2 + 2b cos(ω h)
( b a) sin(ω h)
arctan
1 + ab + ( a + b) cos(ω h)
1 + a2 + 2a cos(ω h)
1 + b2 + 2b cos(ω h)
The different approximations give at ω
1.6 rad/s.
H ( eiω )
Euler
Backward
Tustin with prewarping
Zero order hold
arg H ( eiω )
Gain
Frequency (rad/s)
Phase (deg)
Figure 8.2 The Bode diagrams for the filter in Example 8.3 when h 0.25 continuous
time filter (full); Euler’s method (dashed); backward difference (dotted); Tustin (dashdotted).
Gain
Frequency (rad/s)
The same as Fig. 8.2 but when h
Phase (deg)
Figure 8.3
Fig. 8.2 shows the Bode diagrams for the continuous time system and for the
Euler, backward and Tustin approximations. Fig. 8.3 is the same as Fig. 8.2 but
with h 0.05. The shorter sampling period gives a better approximation.
Problem 8.4
It is assumed that the sample and hold circuit can be approximated by a delay
of h/2 seconds. Further we will allow a decrease of the phase margin of 5○
This approximately corresponds to a decrease of the damping by 0.05 0.15. A
time delay of h/2 seconds gives at the crossover frequency a decrease of
ω c h/2[rad]
This gives
or approximately
Problem 8.5
The transfer function of the integral part of the PID-controller (8.22) is
GI ( s)
Using Euler’s approximation (8.4) gives
Ti ( z 1)
H I ( z)
which is the same integral part in (8.23). The derivative part of (8.22) has the
transfer function
G D ( s)
Using backward difference gives
T ( z − 1)
K Td( z
H D ( z)
K Td( z 1)
z( h + Td / N ) Td / N
which is the same as the derivative part on page 308.
a.
Approximation of the integral part with backward difference gives
Ti( z 1)
H I ( z)
An error will then directly influence the computation of the integral part.
Euler’s approximation gives a delay of one sampling interval before an error
will influence the integral part. The sampling interval is, however, usually
short for digital PID-algorithms.
Euler’s approximation for the derivative part gives
H d ( z)
K N ( z − 1)
A small value of Td can make Hd unstable. Since the D-part of a PIDcontroller sometimes is not used it is necessary that the regulator remains
stable when Td 0.
Problem 8.6
Using the bilinear transformation gives
H T ( z)
(2Ti + h)( z 1)
This is of the same form as (8.24) with
Problem 8.7
The tank process in Problem 2.10 has the transfer function
G ( s)
a.
( s + 0.0197)( s + 0.0129)
At the desired cross over frequency we have
G ( iω c )
arg G ( iω c )
We will use a PI controller of the form
K ( Ts + 1)
Gr ( s)
and we want the gain 1/0.523 and the phase
−15 degrees at ω . This gives
The characteristic equation of the closed loop system is
The roots are s1 2
have a damping ζ
−0.0135 ± 0.0281i and s −0.006. The complex poles
0.43. The zero of the closed loop system is −0.0062.
Tustin’s approximation with warping gives with α
H r ( z)
ω c / tan(ω c h/2)
1.85(α + 0.0067)
(α + 0.0067)( z
Using the rule of thumb from Section 8.2 for choosing the sampling period
gives
h 6 20 seconds
The choice h
12 seems to be reasonable. This gives α
H r ( z)
0.165 and
Time
Figure 8.4 Step response of the tank process when controlled with a continuous time
(solid) and a discrete time PI controller. The sampling interval is 6 (dash-dotted) and 12
seconds (dashed).
Fig. 8.4 shows simulations of the step response of the system controlled with
the continuous time and the approximate discrete time PI-controller when
h 6 and 12 seconds.
Problem 8.9
a.
The continuous time controller is
u( t)
Muc ( t)
− Lx(t).
A discretization is obtained by sampling uc and x and letting u be constant
between the sampling period points i.e. we get
u( kh)
L( I + ( A B L) h/2)
( I − LB h/2) M
4(1
Fig. 8.5 shows the stepresponse of the system when using the continuous controller and the controllers in a) and b) when h 0.25. It is possible to calculate backwards to find out the corresponding damping and natural frequency
for the controllers in a) and b). A discrete time state space representation of
the motor is given in (A.6). Using L  1 2  gives
− Lx(kh)
Using (8.24) and (8.25) give
Muc ( kh)
1− e
e−h
− e− )
− (1 − e− ) 
( h − 1 + e− ) 1 − ( h − 1 + e− )
1 (1
Time
Figure 8.5 Stepresponses for the motor in Problem 8.9 when a continuous time (solid),
a discretized (dash-dotted) and a modified discretized state (dashed) feedback controller
is used when h 0.25.
0.25 and L
 2 4  we get
and for h
 1.75 3  we get
0.25 and L
These two matrices have the characteristic equations
and
From the equations given in Example 2.16 we can calculate the corresponding
 2 4 ) we
continuous time systems. For the discretized controller ( L
get
and for the modified controller ( L
− 4h ) we get
The change in the damping is smaller when the modified controller is used
and the change in the undamped natural frequency is also smaller.
Problem 8.10
a.
We first want to compute a state feedback such that A
teristic equation
Assume L  1 2  then
s2 + (5 +
This gives
Modifying L using (8.16) gives
The characteristic equation of A
( s + 3)( s + 2 +
− B L has the charac-
L( I + ( A B L) h/2)
 17(1 3h) 3 + h 
Fig. 8.6 shows the output when using the discrete time controller in a) for
different values of h. The response when using the modified discrete time
controller from b) is shown in Fig. 8.7.
Problem 8.12
a.
Using (8.16) and (8.17) give
I + (A
Using the backward difference approximation gives
− q− I xˆ(k)
(I
− Ah + K C h)xˆ(k)
(A
(I
This gives
x ( k)
− K C )xˆ(k) + Bu(k) + K y(k)
q−1 x( k) + B hu( k) + K hy( k)
Introduce
Φ 0 x( k 1) + Φ 0 B hu( k) + Φ 0 K hy( k)
 x ( k 1) + 
 u( k ) + 
 y( k)
Time
Figure 8.6 The response of the system in Problem 8.10 when the state feedback con0.1 (dashed) , 0.2 (dash-dotted) and 0.3 (dotted). The
troller in a) is used with h
response for the continuous-time controller is also shown (solid).
Time
Figure 8.7 The response of the system in Problem 8.10 when the modified state feedback
controller in b) is used with h
0.1 (dashed), 0.2 (dash-dotted) and 0.3 (dotted). The
response for the continuous-time controller is also shown (solid).
Solutions to Chapter 10
Problem 10.1
a.
Better sensors, for instance a tachometer with less noise.
Flow control with local feedback.
Temperature control in houses.
Problem 10.2
a.
The time function y( t)
sin(ω t) has the z-transform
Y ( z)
See Table 2.1. Consider a system with
H d ( z)
Y ( z)
The impulse response of Hd will thus be sin( khω ). That this is the correct
answer is easily seen by making long division.
H d ( z)
sin( hω ) z−1 + sin(2hω ) z−2 + sin(3hω ) z−3 + ⋅ ⋅ ⋅
The time function t ⋅ e−t has the z-transform
he−hz
( z e−h)2
Y ( z)
This can be found by looking in a table of z-transforms. The desired system
thus has the z-transform
he−hz
( z e−h)2
H d ( z)
Long division gives
H d ( z)
he−h z−1 + 2he−2hz−2 + ⋅ ⋅ ⋅
Problem 10.3
Using the model of the disturbance gives
y( k + m)
C ( q)
w( k + m).
A( q)
Introduce the identity
q m−1 C ( q)
where deg F
y( k + m)
− 1 and deg G
F ( q)w( k + 1) +
A( q) F ( q) + G ( q)
− 1. Then
qG ( q)
w( k)
A( q)
F ( q)w( k + 1) +
qG ( q)
y( k)
C ( q)
If w( k + 1)
y( k + m) is
w( k + m) are assumed to be zero then the best prediction of
y( k + m)
qG ( q)
y( k).
C ( q)
The operator qG ( q)/ C ( q) is casual since deg C
Let A( q) q 0.5, C ( q) q and m 3 then
(q
deg G + 1.
− 0.5)(q + f q + f ) + g
This gives the system of equations
The predictor at k + 3 given data up to and including k is thus
y( k)
y( k + 3 k)
Let w( k) be zero except at k
y( k)
0.125 y( k)
0 and 5 when it is assumed to be one then
0.5 y( k
− 1) + w(k)
y( k k
Problem 10.4
Using (10.11) we find that the stationary variance of x fulfils (Φ stable)
The stationary covariance exists since the system is stable. Since R 1 is symmetric
P is also symmetric. Introduce
then
This gives
The solution is
The stationary covariance function is
Ex( k + τ ) x( k) T
It remains to compute Φτ . The eigenvalues of Φ are λ 1
the results on matrix functions in Appendix B we get
0.4 and λ 2
where α 0 and α 1 are obtained from
−λ λ (λ − − λ − )
The solution is
and
−3(0.4 − 0.2 )
Finally
rx (τ )
Problem 10.5
From the state space description we get the input-output description of the process
y( k) + a1 y( k
− 1) + ⋅ ⋅ ⋅ + a y(k − n)
c1 v( k
− 1) + ⋅ ⋅ ⋅ + c v(k − n)
where ai
n are the coefficients of the characteristic polynomial of
the matrix Φ . Multiply the equation above by y( k τ ) and take the mathematical
expectation. y( k τ ) is independent of all the terms on the right hand if τ > n + 1.
r y (τ ) + a1 r y (τ 1) + ⋅ ⋅ ⋅ + an r y (τ n) 0.
This is called the Yule-Walker equation.
Problem 10.6
There are two noise sources v1 and v2 that is used to generate y( k). Using Theorem
10.2 we find that the spectral density of y is
H ( z)φ v H T ( z−1)
where z
H ( z)
eiω
z + a
C ( zI − Φ)−1 Γ
and
z+a
z+a
z−1 +a
−1 + a) + (z + b)(z−1 + b)
( z + a)( z
σ 1 ( z + b)( z−1 + b) + σ 2 ( z + a)( z−1 + a)
−1 + a)(z−1 + b)
( z + a)( z + b)( z
Using the spectral factorization theorem (Theorem 10.3) we find that we can
generate the same spectral density by sending white noise through
H1 ( z)
( z + a)( z + b)
this gives the spectral density
( z + c )( z−1 + c )
( z + a)( z−1 + a)( z + b)( z−1 + b)
Identification with φ y gives the relationship
λ 2 (1 + c 2 )
σ 1 (1 + b2) + σ 2 (1 + a2 )
σ 1 b + σ 2a
Problem 10.7
The process is
x( k + 1)
ax( k) + v( k)
y( k)
x( k) + e( k)
This is the same as in Example 10.3 with the exception that Ev( k) e( s) r12δ ( k
s). The covariance function for x will thus be the same but the covariance of y will
contain an additional term due to the correlation between v and e. From Example
10.3 we know that
rx (τ ) a τ
1 a2
The covariance of y is
r y (τ )
E y( k + τ ) y( k)
where it has been used that rex (τ )
rxe (τ + 1)
E [ x( k + τ ) + e( k + τ )] [ x( k) + e( k)]
rx (τ ) + rxe (τ ) + rex (τ ) + re (τ )
rxe ( τ ) in stationarity.
E x( k + τ + 1) e( k)
arxe (τ ) + rve (τ )
rx (τ ) + rxe (τ ) + rxe ( τ ) + re (τ )
E [ ax( k + τ ) + v( k + τ )] e( k)
The last term is zero except for τ
0, where it is r12 . The first term is zero for
τ ≤ 0. It is natural that white noise in the future is uncorrelated with the present
value of x. This gives
aτ −1 r12
rxe (τ )
aτ −1 r12
 a 1−1a2 + 0 + a−τ −1 r12 + 0 τ < 0
 1− a r
a 1−a2 + a
and
r y (τ )
aτ
+ aτ
1 a2
1 a2
The definition of spectral density gives
φ y (ω )
r y (τ ) e−iωτ
− ra (e − 1a)(ea− − a) − 1 − a − ra + 1 − a + r
1 r a + r (1 − a ) − r ( e − a)( e− − a) + r a( e − a)( e− − a)
a( e − a)( e− − a)
a2
(1)
where it has been used that
a τ e−iωτ
( eiω
The spectral density for
y( k)
1 a2
a)( e−iω
− a)
− a ε ( k)
is (see Theorem 10.2)
λ 2 ( eiω
2π ( eiω
φ y(ω )
− c)(e− − c)
− a)(e− − a)
(2)
Identification of (1) with (2) gives the relation
−a −r
a
− r a −λ c
1 + a2
1 + a2
+ r2 a
a
a
λ 2 (1 + c 2 )
r2 (1 + a2 ) + r1
r2 a
A more elegant solution
λ 2 (1 + c 2 )
The ouput can be written as
y( k)
− 2ar
 H ( q) 1   v( k) 
e( k)
where H ( z)
z a
The spectral density of y is given by
r12   H ( z−1 ) 
 H ( z) 1   1
r1 H ( z) H ( z−1) + r12 H ( z) + r12 H ( z−1) + r2
2π ( z a)( z−1 a)
z a
a
1 z( r12
− r a) + r − 2ar + r (1 + a ) + z− (r − r a)
( z − a)( z− − a)
which gives the same equations as in the previous method
− 2ar
ar − r
λ 2 (1 + c 2 )
r2 (1 + a2 ) + r1
Problem 10.8
The process can be written as
y( k)
− c e(k)
−a
a
− c e(k) + e(k)
−a
where
x( k + 1)
ax( k) + ( a
x( k) + e( k)
− c)e(k)
Using Problem 10.7 we find that
( a c )2
a c
r y (τ )
with a
0.5 we get for τ
0.7 and c
r y (τ )
a τ r1
+ aτ
1 a2
(0.7) τ
− 0.49 + 0.7 (0.7)
0.36(0.7) τ
( a c )2
1 a2
1 a2
The variance can also be obtained from Theorem 10.4.
r y (0)
r y (0)
1 + c 2 2ac
1 a2
Further (10.17) gives
where z
1 (z
2π ( z
− c)(z− − c) )
− a)(z− − a)
eiω .
Problem 10.9
The variance of the process
y( k)
e( k)
e( k)
a0 q2 + a1 q + a2
can be determined using Theorem 10.4. The formula for I2 gives
1 + (0.2)2 + 0
2(0.2 + 0)
e1
and
r y (0)
(1
0.49)1.7 1.5 ⋅ 1.5(1
The recursive scheme in Section 10.4 can also be used to get
This gives
These calculations have been performed in high precision and thereafter round-ed.
Problem 10.10
The process is
y( k)
r y (τ )
This gives
r y (0)
e( k)
− 2e(k − 1) + 3e(k − 2) − 4e(k − 3)
Ey( k + τ ) y( k)
− 2e(k − 1) + 3e(k − 2) − 4e(k − 3))
e( k) + 4e( k − 1) + 9e( k − 2) + 16e(6 − 3)
E ( e( k)
+ crossterms
The mean value of the crossterms are zero since
E e( k + τ ) e( k)
In the same way
r y (1)
r y (2)
1 ⋅ 3 + ( 2) ⋅ ( 4)
r y (3)
1 ⋅ ( 4)
r y (4)
1 ⋅ ( 2) + ( 2) ⋅ 3 + 3 ⋅ ( 4)
Problem 10.11
a.
H ( z)
Using Theorem 10.2 we get
φ y(ω )
H ( eiω )φ u (ω ) H ( e−iω )
−a
eiω
e−iω
H ( eiω )
−a
1 + a2
2a cos ω
a2
−a
H ( e−iω )
eiω + e−iω
a⋅2
1 + a2
2a
Identifying with the desired spectral density gives
2π (1 + a2 )
2π (2a)
Divide the two equations
1 + a2
2a
a2 +
a+1
a
The desired filter is
H ( z)
Theorem 10.4 CCS
Var y
Problem 10.12
x( k + 1)
2a
 x( k) +  0  u( k) + v( k)
The stationary covariance is given by
Problem 10.13
Example 10.4 shows that the filter
H ( z)
−a
gives the spectral density
Φ(ω )
In this case r1
1. Identify with the desired spectral density gives
This gives a
2π 1 + a2 2a cos ω
0.9 and b
Solutions to Chapter 11
Problem 11.1
For the system we have
−aΦ(t kh)
Φ( t kh)
Φ( kh kh)
This differential equation has the solution
e−a( t−kh)
Φ( t kh)
and
e−a( t−s)b ds
Γ( t kh)
(1
a
− e−
a( t kh)
The discrete time loss function is obtained from CCS (11.6)-(11.8), which gives
e−2a( s−kh) ds
e−a( s−kh) (1
a
− e−
− e−
2ah
(1
2a2
a( s kh)
− e−
ah 2
− e− − ) + ρ ds
3 − 4e− + e−
2a
(1
a2
(1
2a
a2
a( s kh) 2
ah
2ah
Notice that there will be a Q12 term even if Q12c
Problem 11.2
Sample the system
e
e
Sample the loss function Q1c
0. Using (11.6)-(11.8) we get
Problem 11.3
The Riccati equation(11.17) gives
s( k)
a2 s( k + 1) + 1
− a b bs(sk(k + )1)
and (11.19) gives
b−2 ba
L( k)
a
which gives the controller
u( k )
− L(k)x(k) − ab x(k) − a x(k)
The minimum loss is given by
x(0)2 s(0)
x(0)2
and the closed loop system is
x( k + 1)
The state is zero after one step and the resulting controller is thus a dead-beat
controller.
Problem 11.5
a.
The loss function is Σ y2 + ρ u, i.e.
The steady state value of the Riccati equation is obtained from
Let
S Γ( Q2 + Γ T S Γ)−1 Γ T S Φ
For the inventory model we get
Pole
Figure 11.1
The pole in Problem 11.5 as a function of the control weighting ρ .
The solution is
The feedback vector is
K (ρ )  1
where
K (ρ )
Control weighting rho
The dynamics of the closed loop system is
− K (ρ ) − K (ρ )
The poles are obtained from
(λ
− 1)(λ + K (ρ )) + K (ρ )
λ (λ
− 1 + K (ρ ))
There is one pole in the origin and one in 1 K (ρ ). For ρ
0 then 1 K (ρ )
0 and as ρ
then 1 K (ρ )
1. Fig 11.1 shows how the pole varies as
a function of ρ .
The poles of the closed loop system can also be determined from (11.40). For the
inventory system
A( z) q( q 1)
B ( z)
and we get
ρ ( z−2
− z− )(z − z) + 1
r( z2 + p1 z + p2 )( z−2 + p1 z−1 + p2 )
a)
y(k), u(k)
y(k), u(k)
Time
y(k), u(k)
y(k), u(k)
Time
Figure 11.2 The output (solid) and the control signal (dashed) for the system in Prob0, b) ρ
5 and d) ρ
lem 11.5 when a) ρ
r( p1 + p1 p2 )
Since r
r(1 + p2 + p2 )
0 then the first equation implies that p2
0 and we get
r(1 + p2 )
which has the solution
The poles of the closed loop system are thus one at the origin and one at
It is easily seen that this is the pole in 1
− K (ρ ).
Problem 11.6
The system has the transfer function
H ( z)
Only the output and the control signals are penalized in the loss function. The
closed loop system has poles that are determined by the stable roots of
ρ + H ( z) H ( z−1)
This gives
Figure 11.3
Re
The closed loop poles in Problem 11.6 when ρ is varying.
ρ (0.68z4
Fig. 11.3 shows the closed loop poles when ρ is varying.
Problem 11.9
Solving LQ-problem for system of higher order than one by hand causes much
work. With numerical packages, like Control toolbox for Matlab, the design is
significantly simplified for the control engineer. The following Matlab-macro illustrates the solution of Problem 11.9, i.e. sampling of the system, sampling of
the loss function and solving the Riccati equation.
%Macro for Problem 11.9 CCS
alpha=0.001;
A=[0 1; 0 -alpha];
%Transform continuous-time LQG problem to the corresponding
%discrete-time LQG problem
[Phi,Gam,Q1,Q2,Q12,R1,Je] = lqgsamp(A,B,h,Q1c,Q2c,Q12c,zeros(2,2));
Gain margin
Control weighting rho
Figure 11.4 The gain margin β min (dashed) and β max (solid) from (11.37).
%Linear quadratic regulator design for discrete-time system
[L,Lv,S] = lqrd(Phi,Gam,Q1,Q2,Q12);
The design gives
Problem 11.10
In Problem 11.5 we have determined r, then
Equation (11.37) may be used to get the exact values of the gain margin. With
A( z)
P ( z)
we get
− z + β (z
z( z + (β p1 + β
I.e., the system is stable if
β min and β max are also shown in Fig. 11.4.
Problem 11.11
The system is
x( k + 1)
y( k)
0.5x( k) + v( k)
x( k) + e( k)
β max
Pole
State weighting r1
Figure 11.5 The pole of the Kalman filter in Problem 11.11 for different values of r1
when r2 1.
Theorem 11.5 gives
 x( k + 1 k)
K ( k)
 P ( k + 1)
that the Kalman filter is defined by
0.5 x( k k
− 1) + K (k)( y(k) − xˆ(k k − 1))
0.5P ( k)
r2 + P ( k)
0.25P ( k) + r1
− 0.25PP((kk))
P (0)
x(0
The dynamics of the filter is determined by
− K ( k)
The steady state variance is given from
P 2 + (0.75r2
r2 + P ( k)
Consider three cases r1 > r2 , r1
r2 and r1 < r2 . In the first case P
r1 and
Φ K C 0. In the second case P 1.13r1 and Φ K C 0.23. Finally if r1 < r2
then P 1.33r1 and Φ K C 0.5. Fig. 11.5 shows Φ K C for different values
of r1 when r2 1.
Additional problem
Suppose that the system in Problem 11.11 has a control signal u( k), i.e. the
system is
x( k + 1) 0.5x( k) + v( k) + u( k)
y( k)
x( k) + e( k)
Determine a steady-state LQG-controller when Q1
0 and Q2
Solution to the additional problem
Equation (11.17) and (11.19) gives
which has the solution
(0.75ρ − 1) + 4ρ
(0.75ρ
Using the Kalman filter from Problem 11.11 and the LQ-controller gives
x ( k + 1 k)
Φ x( k k 1) + Γ u( k) + K ( y( k) C x( k k 1))
Φ x( k k 1) + Γ ( L x ( k k 1)) + K ( y( k) C x( k k
(Φ Γ L K C ) x( k k 1) + K y( k)
and thus
U ( q)
− L(qI − Φ + Γ L + K C )− K Y (q)
H ( q)
reg
Problem 11.12
a.
Equation (11.47) gives
Φ P ( k)Φ T + R 1
P ( k + 1)
− Φ P ( k) C
( R 2 + C P ( k) C T )−1 C P ( k)Φ T
we get
p11 ( k + 1)
p12 ( k + 1)
p22 ( k + 1)
p11 ( k) + 2p12 ( k) + p22 ( k)
− (p
11 ( k) +
p12 ( k))2
p11 ( k)
p11 ( k) p22( k) p12 ( k)2
p11 ( k)
p12 ( k)( p11( k) + p12 ( k))
p12 ( k) + p22 ( k)
p11 ( k)
p22 ( k) + 1
− p p ( k)
p11 ( k + 1)
1 + p11 ( k + 1)
Further
K ( k)
0 K (0)
 p11 ( k) + p12 ( k)  1
p11 ( k)
p12 ( k)
[1 0]T i.e. K is timevarying. The steady state value of P is
The poles of the filter are found from det(λ I
has a dead beat response.
− (Φ − K C ))
0 The filter
p12(k)
p11(k)
p22(k)
Time
Time
Figure 11.6
The elements of the variance matrix in Problem 11.12.
The initial values of the filter are
x(0
and assume that P (0) 3I. Fig. 11.5 shows the elements of the covariance
matrix as a function of time.
Problem 11.14
Introduce the notation
and
The state at time k + 3 can now be determined
x( k + 3)
Φ x( k + 2) + Γ 1 v( k + 2) + Γ 2
Φ 2 x( k + 1) + ΦΓ 1 v( k + 1) + ΦΓ 2 + Γ 1 v( k + 2) + Γ 2
Φ 3 x( k) + Φ 2 Γ 1 v( k) + ΦΓ 1 v( k + 1) + Γ 1 v( k + 2)
The best estimate of x( k) given y( k) is determined from (11.50). Since v( k),
v( k + 1) and v( k + 2) are independent of y( k) then the best estimate of x( k + 3)
is given by
x ( k + 3 k)
Φ 3 x( k k) + (Φ 2 + Φ + I )Γ 2
 x ( k k) + 
The variance of the estimation error is
P ( k + 3 k)
Φ 3 P ( k k)(Φ 3) T + 0.01(Φ 2 Γ 1 Γ 1 (Φ 2 ) T + ΦΓ 1 Γ 1 Φ T + Γ 1 Γ 1 )
 P ( k k) 
If x(0) is known then P (0 0)
0 and
y(3)
 1 3  x(0) + 4.5
and the variance of the prediction error is 0.05.
Problem 11.15
x( k + 1)
ax( k) + v( k)
x( k) + e( k)
y( k)
cov e
We use the exponential smoothing estimator
− 1 k − 1) + (1 − α ) y(k)
(1 − α ) q
y( k) − x( k)
(1 − α ) q
v( k) + e( k) −
v( k)
q−a
q−a
− (q − (aq)(−q1−) α ) v(k) + (1q− αα)q e(k)
α x( k
x( k k)
[ x( k k) x( k)]
Using Theorem 10.4 we get
(1 + a)(1 + α )(1
var x( k k)
− aα )
Minimize with respect to α , use Maple.
σ a(1 + a) + 1
σ (1 + a)2 + 1
2 (1 + a) + a
σa
Kalman with direct term
− 1 k − 1) + K ( y(k) − C Φ xˆ(k − 1 k − 1)
( I − K C )Φ x( k − 1 k − 1) + K y( k)
Φ x( k
x ( k k)
This will have the same form as the exponential smoothing estimator only when
a 1.
Kalman variance
a2 P 2
P a2 P + 1
(1
− a )P − 1
−a P
This gives
− σ (1 − a ) + 1
1 + 2σ (1 + a2 ) + σ 2 (1
−a )
The gain in the Kalman filter is
aP
var x( k k)
Numerical values: σ
1 a
a
var x( k k)
Kalman:
var x( k k)
Problem 11.16
The scalar state equations are
x( k + 1)
y( k)
ax( k) + u( k) + v( k)
v( k)
v1 ( k) + mv ;
x( k) + e( k)
e( k)
e1 ( k) + me ;
Ee1
Let
and
a 1
 X ( k + 1)
y( k)
me 
 X ( k ) +  0  u( k ) +  0  v ( k )
1  X ( k) + e1
The observability matrix is then
 a
a
a+1
rank Wo
This means that me and mv are not both observable and no Kalman-filter can be
designed. It is, however, possible to design a second order observer with reconstruction of a linear combination of me and mv . Redefining the state vector X
as
where
gives
 X ( k + 1)
y( k)
x + me
(a
e
a 1
 X ( k) +  1  u( k) +  1  v1 ( k)
  X ( k) + e1 ( k)
Reconstruction of x1 and m is possible if these states are observable. The observability matrix is
a 1
from which follows that rank Wo
Problem 11.17
The constants a1 and a2 are determined by the following two conditions
The correct mean value should be obtained.
2. The variance should be minimized.
Condition 1 gives that
a1 + a2
The variance of the estimation error is
E ( x( k)
− xˆ(k))
a2 ⋅ 1 + a2 ⋅ 9
− a x(k) − a e (k) − a x(k) − a e (k))
+ (1 − a ) ⋅ 9 10a + 9 − 18a
E ( x( k)
a2
Taking the derivative with respect to a1 gives the condition
20a1
The estimator is thus
x ( k)
a1
y1 ( k) +
y2 ( k)
The minimum value of V is
Using only y1 gives the variance 1 and only y2 the variance 9. Thus, a combination
of the measurements gives a better result than using only the best measurement.
Assume that the a priori estimate of x is zero and the variance of x is p, i.e.
x (0 0)
From (11.50) – (11.54) we get
and
P (1 0)
P (0 0)
and
K (1)
P (1 0) C T R 2 + C P (1 0) C T
This gives
x( k k)
y1 ( k) +
y2 ( k) +
x( k k
If p is large then the weights for y1 and y2 will be those calculated above. In the
example R 1
0 the steady state gain will then be zero since the estimate will
approach the true value of x.
Problem 11.20
x( k + 1)
y( k)
 x( k) +  1  u( k)
 0.5 0.38  x( k)
The steady state solution is obtained from (11.17). Using Matlab we get the
An alternative solution is to use (11.40)
ρ A( z) A( z−1) + B ( z) B ( z−1)
zB ( z)
P ( z)
2z(0.5z + 0.38)
Now we have to find L such that
(Φ
rP ( z) P ( z−1)
z( z + 0.76)
where controllable canonical form have been used. This gives
Problem 11.21
a.
First assume that η
0 and use linear quadratic design with
Theorem 11.1 gives the Riccati equation
− L(k)
S ( k)
L( k)
0.5S ( k)
ρ + S ( k + 1)
S ( k + 1)
This gives
S( N )
S (1)
L(1)
S (2)
L(0)
− L(k)x(k) − Ly(k). For different values of ρ we
The control law is u( k)
get
L(0)
L(1)
In this case η
1 and x( k) is reconstructed using a Kalman filter
x( k + 1)
0.5 x( k) + u( k) + K ( k) y( k)
− xˆ(k)
E x2 (0)
Theorem 11.5 gives
K ( h)
P ( k + 1)
with P (0)
0.5P ( k)
1 + P ( k)
0.25P ( k)
− 0.5P(k) K (k)
1. This gives
P ( k)
K ( k)
− L(u)xˆ(k).
The control law is u( k)
Problem 11.22
x( k + 1)
x( k) + v( k)
  x( k) + e( k)
y( k)
a.
x( k + 1 k)
x( k k
− 1) + K y(k)
(11.47) ⇒
+ (σ 1 + σ 2 ) p
σ 1 σ 2 + p(σ 1 + σ 2 )
Solutions to Chapter 12
Problem 12.1
The m-step ahead predictor is given by
qG ( q)
y( k)
C ( q)
y( k + m k)
where G is obtained from the identity (12.17) and the variance of the prediction
error is given by (12.19). For m 1 we get
which gives
G ( q)
The predictor is then given by
−0.2q + 0.1q y(k)
y( k + 1 k)
and the variance of the prediction error
E y2 ( k + 1 k)
q( q2
This gives
( q2
F ( q)
− 1.2q + 0.4)(q + f ) + g q + g
G ( q)
and
σ 2 (1 + f 1 )
E y2 ( k + 2 k)
3 we get
q2 ( q2
( q2
− 1.2q + 0.4)(q
which gives
F ( q)
G ( q)
E y ( k + 3 k)
σ 2 (1 + f 1 + f 2 )
Using Theorem 10.4 we can compute the variance of y to
var( y)
(1+1.42+0.52 )(1+0.4)+2( 1.4 1.4 ⋅ 0.5)⋅ 1.2+ 2 ⋅ 0.5(1.22 0.4(1+0.4))
(1 0.42 )(1 + 0.4) + ( 1.2 + 1.2 ⋅ 0.4) ⋅ 1.2
This implies that the prediction variance is almost the same as the variance of y
when m ≥ 3.
Problem 12.2
The identity (12.17) is
q m−1 ( q + c )
( q + a)( q m−1 + f 1 q m−2 + ... + f m−1) + g0
This gives
a + f1
a f1 + f2
a f m−2 + f m−1
a f m−1 + g0
−a
(−a)( c − a)
(−a) ( c − a)
The solution is
(−a)
( a) m−2 ( c
The m-step ahead predictor is
− a)
− (c − a)
( a) m−1 ( c
y( k + m k)
− a)q y(k)
and the variance of the prediction error is
E y( k + m k)
− a) + a (c − a)
1−a −
1 + ( c − a)
1−a
σ 2 (1 + ( c
+ ⋅ ⋅ ⋅ + a2( m−2)( c
2( m 1)
− a) )
Problem 12.3
a.
The C-polynomial has a zero outside the unit circle. Example 12.1 shows how
the C-polynomial should be changed. It is replaced by
C ∗ ( z)
5( z + 0.2)
The equivalent process is thus
y( k)
− 0.9 y(k − 1)
5( e( k) + 0.2e( k
The two-step-ahead predictor is obtained from
q( q + 0.2)
This gives
(q
− 0.9)(q + f ) + g
F ( q)
G( g )
This predictor is
y( k + 2 k)
and
E y2 ( k + 2 k)
y( k)
25(1 + 1.12 )
Problem 12.4
Using the data given up to time k 7 it is possible to calculate y( k) and zd ( k)
z( k) y( k). zd is the deterministic part of z.
z( k )
y( k)
z d ( k)
The prediction of the demand for August through November is
z(8 7)
zd (8) + y(8 7)
z(11 7)
zd (11) + y(11 7)
We thus need the 1, 2, 3 and 4 step ahead predictors of y. Those are given by
solving the identity (12.17) and give
F ( q)
G( g )
The prediction is
qG ( q)
y( k)
C ( q)
y( k + m k)
g0 y( k) + g1 y( k
which gives the predicted values and their standard deviation σ .
y(7 + m 7)
zd(7 + m)
z( 7 + m 7 ) σ
Problem 12.5
The polynomials are
It is easily seen that C is a stable polynomial, e.g. by the stability triangle. This
is a necessary condition for the minimum variance design method.
The pole excess is
d deg A deg B 2
The Diophantine equation is
q d−1 C ( q)
q( q + 0.8q + 0.25q)
A( q) F ( q) + G ( q)
( q3
+ 0.5q)( q + f 1 ) + ( g0 q2 + g1 q + g2 )
Identifying coefficients of powers of q gives
Solving these equations
The minimum variance regulator is obtained from (12.27)
u( k )
− B (G)(F)(q) y(k) (q +.55q)(q + .1.8) y(k)
The loss function is
0.52 (1 + 1.82 )
Problem 12.6
The noise sequence has a non zero mean value. Introduce
2 + ε ( k)
u + u( k )
e( k)
u( k )
where ε ( k) is zero mean white noise. The process is now
y( k)
Choose u
gives
− 0.5 y(k − 1)
u + u( k 2) + 2 + ε ( k) 0.7(2 + ε ( k 1))
u( k 2) + ε ( k) 0.7ε ( k 1) + u + 0.6
−0.6 and the problem is reduced to the standard problem. The identity
y( k)
and
u( k )
y( k)
Problem 12.7
a.
The identity gives
F ( q)
G ( q)
−a
a( a − c ) q
and the minimum variance controller is
u( k )
The expression above gives the optimal controller u( k)
0 if a
0. The
process is then a moving average process of first order. This implies that
u( k 2) cannot be used to decrease the variance of y( k).
(
− qa+a(−−)a) y(k)
Figure 12.1
Time
The output of the open loop system in Problem 12.8.
Problem 12.8
a.
The identity gives for d
F ( q)
and for d
G ( q)
F ( q)
G ( q)
The minimum variance controller is
u( k )
− B (G)(F)(q) y(k)
and the minimum variance in the two cases are
Fig. 12.1 shows the output of the open loop system for one realization of the
noise sequence. The output is drifting since the A-polynomial contains an
integrator. Fig. 12.2 and Fig. 12.3 shows the output and the control signal
for the same noise realization when the minimum variance controller is used
with d 1 and d 2.
Problem 12.9
a.
Assume that
H ( z)
z+ a
Sending white noise e( k) through this filter gives the spectral density (see
Theorem 10.2)
φ (ω )
2π 1 + a2 + 2a cos ω
Time
Figure 12.2 The output and the control signal when d
controller is used for the process in Problem 12.8.
1 and the minimum variance
Time
Figure 12.3
Same as Fig. 12.2 but when d
This implies that λ
1 and a
process is now described by
y( k)
( q2 + 0.1q
0.6 gives the desired spectral density. The
− 0.3) y(k)
Use the controller
u( k )
e( k) + u( k)
( q + 0.6) u( k) + q e( k)
− K y(k)
This gives
− 0.3) y(k) −(q + 0.6) K y(k) + q e(k)
y( k)
e( k)
q + (0.1 + K ) q + (0.6K − 0.3)
The system is stable if
( q2 + 0.1q
Theorem 10.4 gives an expression for the variance of a second order process.
I2 ( K )
1 we get I2
(1.3
− 0.6K )((0.7 + 0.6K ) − (0.1 + K ) )
The minimum value of I2 is obtained from
This gives the third order equation
2.009, 1.839 and 0.004. Only the root K
which has the roots K
gives a stable closed loop system. The value of the variance is
I2 (0.004)
The minimum variance is Ey2
1 since d
Problem 12.10
a.
With the proportional controller
u( k )
− K y(k)
we get the closed loop system
e( k)
q2 + ( K 0.25) q + 0.5
y( k)
Using the results in Theorem 10.4 gives
2(1 ⋅ 0.5 + 0.5 ⋅ 0)
e1
and
1.25 ⋅ 1.5 ( K 0.25)
(1 0.25) ⋅ 1.5 ( K 0.25)2(1 0.5)
0.5(1.75 K )(1.25 + K )
Taking the derivative of I2 and putting the derivative equal to zero leads to
the equation
I2 ( K )
with the solutions K
1 and K
I2 (1)
1 gives the variance
This is minimal variance for the present control law. With minimum variance
control we would get E y2 1.
From Example 3.2 we find that the closed loop system is stable if
Both K 3.25 (the second root in a.) and K 2.125 give an unstable closed
loop system and the calculation of I2 is no longer valid.
Problem 12.11
y( k)
gives
− 1.5 y(k − 1) + 0.7 y(k − 2) u(k − 2) − 0.5u(k − 3) + v(k)
A( q) q − 1.5q + 0.7q
B ( q) q − 0.5
Note that the process zero (0.5) is stable and well damped. This means that the
process zero can be cancelled without problem. The degree conditions gives
deg Am
a.
v( k)
− deg B
≥ deg A deg B 2
deg Ao ≥ 2 deg A deg Am deg B +
0; Deadbeat Control
Am ( q)
B + ( q)
B − ( q)
B m ( q)
Ao ( q)
The polynomials R 1 and S are obtained from the Diophantine equation
A( z) R 1 ( z) + B − ( z) S ( z)
Recalling the condition deg S
( z3
deg A
A m ( z) A o ( z)
− 1 the equation becomes
+ 0.7z)( z + r1 ) + s0 z2 + s1 z + s2
This gives the regulator
u( k )
where
R ( q)
T ( q)
Assuming that uc ( k)
u( k )
T ( q)
u c ( k)
R ( q)
R 1 ( q) B + ( q)
− R((q)) y(k)
( q + 1.5)( q
0 gives
− R((q)) y(k) − (q1+ 1.5)(q1−05q5)
v( k)
e( k)
− 0.2e(k − 1); Minimum variance The polynomial C is given by
C ( q) q − 0.2q
The minimum variance control law can be written as
− B (G)(F)(q) y(k)
u( k )
where the polynomials F and G are obtained from
q d−1 C ( q)
A( q) F ( q) + G ( q)
deg F
deg G
which in this case is
q( q3
( q3
+ 0.7q)( q + f 1 ) + g0 q2 + g1 q + g2
which yields the equations
The minimum variance controller is thus given by
− (q1− 0.5)(q0+91q3) y(k)
u( k )
The output is in the deadbeat case given by
y( k)
e( k)
e( k)
q−4 C ( q) R 1 ( q) e( k) C ∗ ( q−1) R ∗ ( q−1) e( k)
(1
− 0.2q− )(1 + 1.5q− )e(k)
(1 + 1.3q−1
− 0.3q− )e(k)
which is a moving average (MA) process. The variance of the output is then
simply calculated as
E ( y2 )
(1 + 1.32 + 0.32 )σ 2
In the minimum variance case the output is
y( k)
q−( d−1) F ( q) e( k)
F ∗ ( q−1) e( k)
(1 + 1.3q−1) e( k)
which also is an MA process. The output variance is
E ( y2 )
(1 + 1.32 )σ 2
Time
Figure 12.4 The output and the sum of the square of the output for the two regulators
in Problem 12.11. The deadbeat controller is used in the upper diagram and the minimum
variance controller in the middle diagram. The accumulated loss functions are shown in
the lower diagram, deeadbeat (dashed) and minimum variance (solid).
Fig. 12.5 shows the output and the sum of the square of the output for the
regulators in a and b.
Problem 12.12
Introduce the polynomials
A( q) D ( q)
B1 ( q)
B ( q) D ( q)
C1 ( q)
and the noise e1
A1 ( q)
A( q) C ( q)
λ e. The system can then be written as
y( k)
u( k ) + λ
e1 ( k)
(12.1)
Since A, C and D are monic, A1
AD and C1
AC will also bo monic. Let
d. The minimum-variance controller for
d1 deg A1 deg B1 deg A deg B
the rewritten system (12.1) is then calculated as
u( k )
where
q d1−1 C1 ( q)
− B (Gq)(F )(q) y(k)
(12.2)
A1 ( q) F1( q) + G1 ( q)
(12.3)
Equation (12.3) is equivalent to
q d−1 A( q) C ( q)
A( q) D ( q) F1( q) + G1 ( q)
which implies that A must divide G1 , i.e. G1 ( q)
gives the minimum-variance control law
u( k )
A( q) G ( q). Putting F
− B (G)(D)(A()q)(q) y(k) − B (qG− ()qD ()qA− ()qF ()q− ) y(k)
where F and G satisfy the equation
q d−1 C ( q)
D ( q) F ( q) + G ( q)
with deg F
d 1 and deg G
deg D 1. We see that if A
reduces to the normal minimum variance controller.
D the controller
Problem 12.13
In this case
A( q)
q+a
B ( q)
C ( q)
D ( q)
The identity in Problem 12.12 gives
The minimum variance controller is thus
u( k )
− c(qbq a) y(k) − b y(k) − ac y(k − 1)
Problem 12.15
We have
y1 ( k)
y2 ( k)
u( k 1 ) +
e( k)
y1 ( k)
u( k 2) + ∗ e( k 1)
u( k − 2 ) +
e( k − 1)
(1 + α q− )(1 − 0.7q− )
(1 + α q− )(1 − 0.7q− )
To normalize the notations it is convenient to introduce a new noise ε ( k)
a.
Assume that y1 can be measured. The minimum variance controller for y1 is
then
−0.2 y (k)
u( k )
e( k 1).
The variances of y1 and y2 are
The minimum variance controller for y2 when y2 is measurable is obtained
by using the identity
(1 + α q−1 )(1
This gives
− 0.7q− )(1 + f q− ) + q− ( g
F ∗ ( q−1)
This gives
y2 ( k)
(1 + q−1 )ε ( k)
and
u( k )
G ∗ ( q−1)
and the controller
ε (k
y2 ( k)
− 1) + e(k − 2)
(1 + q−1)(1 + α q−1 ) e( k + 1)
y1 ( k)
The variances of the two signals are
1 + (α + 1)2 + α 2
In this case both y1 and y2 are measurable and y1 ( k) will contain more recent
information about the noise process than y2 ( k). It is now only necessary to
predict one step ahead. Introduce the identity
This gives
y2 ( k + 2)
e( k + 2) + ∗ ∗ u( k)
e( k + 2) +
e( k + 1) + ∗ ∗ u( k)
y1 ( k)
ε ( k + 1)
u( k
This gives
y2 ( k + 2)
y ( k)
u( k
ε ( k + 2) + ∗ 1 ∗ y1 ( k) + ∗ ∗ ∗ C ∗
ε ( k + 2) +
ε ( k + 2) +
− 1) + ABA u(k)
− G q − u( k )
y ( k ) + B ∗ u( k )
The variance is thus minimized when
u( k )
y1 ( k)
which is an admissible control law. For the specific system we get
u( k )
y1 ( k)
With this control law we get
ε ( k)
y2 ( k)
e( k
and
y1 ( k)
(1 + α q−1 ) y2 ( k + 1)
(1 + α q−1 )ε ( k + 1)
(1 + α q−1 ) e( k)
The variances of the two output signals are
and
We can thus conclude that the extra measurement will greatly improve the
control of the output of the system.
Problem 12.16
The same arguments can be used in this case as when deriving the normal minimum variance controller. Assume that the system has a stable inverse. Also
assume that deg A deg B deg D i.e. that the process can be written as
A∗ ( q−1 ) y( k)
B ∗ ( q−1)u( k
− d) + C (q− )e(k) + D (q− )v(k − d)
The identity (12.17) can be used here also and gives
y( k + d)
e( k + d) + ∗ u( k) + ∗ v( k)
F ∗ e( k + d) + ∗ e( k) + ∗ u( k) + ∗ v( k)
F ∗ e( k + d) + ∗
y( k)
u( k d)
v( k
+ ∗ u( k) + ∗ v( k)
F ∗ e( k + d) + ∗ y( k) + ∗ ∗ C ∗ q−d G ∗ u( k)
+ ∗ ∗ C ∗ q−d G ∗ v( k)
∗ e( k + d) + 1 ( G ∗ y( k) + B ∗ F ∗ u( k) + D ∗ F ∗ v( k))
The minimum variance controller is thus
u( k )
y( k)
v( k)
Problem 12.17
A( q) y( k)
A( q)
B ( q)u( k) + C ( q) e( k)
B ( q)
C ( q)
LQG-Control: Minimize E ( y2 + ρ u2). Let P ( z) be the closed loop system characteristic equation (12.45)
rP ( z) P ( z−1)
ρ A( z) A( z−1) + B ( z) B ( z−1 )
P ( z) contains stable zeros of the right hand expression (Lemma 12.1)
− 0.9)(z− − 0.9) + 1 ⋅ 1
r( z + p1 )( z−1 + p1 )
ρ (z
r(1 + p2 ) + rp1 z + rp1 z−1
This gives the system of equations
r(1 +
which has two solutions, one of which gives a stable P ( z)
(
Determine the LQG-regulator by means of pole placement:
P ( z)
C ( z)
Control law:
u( k )
− R((q)) y(k)
where S (0) 0. See computational procedure in Section 12.5.
The Diophantine equation to be solved is
− 0.5) (z − 0.9)(z + r ) + zs
( z + p1 )( z
This gives
The solution is given by
which results in the controller
u( k )
y( k)
For the closed loop system we get
y( k)
y( k)
u( k ) +
e( k)
e( k)
Theorem 10.4 gives
Var y
e( k)
− R y(k)
e( k)
e( k)
e( k)
The input u is
u( k )
− R y(k) − R ⋅ R e(k) − S e(k) − q s+qp
e( k)
Theorem 10.4 gives
Var u
In the following table the calculations are summarized for ρ
Var y
0.1 1 and 10.
Var u
Problem 12.24
From Example 12.16 we get a polynomial description of the process
A( q) y( k)
where
B ( q)u( k) + C ( q) e( k)
A( q)
B ( q)
C ( q)
The minimum variance regulator is given by the Diophantine
where
q d−1 B ( q)
P ( q)
The solution is
and the minimum variance regulator is
− c + 1 y(k)
u( k )
− design we use a state-space description of the process (see Example
x( kh + h) x( kh) + hu( kh) + v( kh + h) − v( kh)
y( kh)
x( kh) + ε ( kh)
To obtain the LQ-controller we have to sample the loss function.
From (11.6) - (11.8) we get
The Riccati equation is
Φ T S Φ + Q1 LT ( Q2 + Γ T S Γ) L
( Q2 + Γ T S Γ)−1 (Γ T S Φ + Q12 )
and the solution is
The closed loop system has a pole in
and
P ( z)
To get the regulator we solve the Diophantine (see p. 481)
( q + p1 )( q + c )
(q
S (0)
− 1)(q + r ) + hs q
( p1 + c + 1 + p1 c )
which has the solution
and thus
u( k )
h ( p1
y( k)
Problem 12.29
The process is described by
B ( z)
C ( z)
A( z)
a.
To get the minimum variance controller we use the identity
A( z) F ( z) + G ( z)
This gives
The dead-beat controller is obtained from the identity
A( z) F ( z) + G ( z)
which gives
G ( z)
G ( z)
The control law is
zd−1 C ( z)
Minimum variance control
var( y)
Dead-beat control
(1 + 0.42 )σ 2
var( y)
Problem 12.30
a.
Assume that C
0 and use the identity
This gives
−a
and the controller
ay
The closed loop system becomes
y( k)
e( k) + ce( k
and the variance
var( y)
Assume that C ( z)
is given by
z + c. The minimum variance controller for this model
−a
The closed loop system is now
e( k)
y( k)
which has the variance, see Theorem 6.4
var( y)
(1 + c 2 ) 2c c
It is better to use a guessed value c if
Problem 12.31
The C polynomial has a pole outside the unit circle. Replace C by a new polynomial obtained from spectral factorization
C ( z) C ( z−1)
(z
− 1.25)(z− − 1.25)
The new process is
y( k)
1.252 ( z
− 0.8)(z− − 0.8)
ε ( k)
where ε has the standard deviation 1.25.
a.
The 2-step ahead predictor is given by
zm−1 C ( z)
z( z
A( z) F ( z) + G ( z)
( z2
− 1.1z + 0.3)(z + f ) + g z + g
which gives
F ( z)
G ( z)
and
0.03( q 3)
y( k)
qG ( q)
y( k)
C ( q)
y ( k + 2 k)
The prediction error variance is
1.252(1 + 0.32 )
Problem 12.32
a.
2. This gives the identity
zC ( z)
where deg F
z3 ( z
1 and deg G
( z3
The solution is
A( z) F ( z) + G ( z)
2, i.e.,
− 1.7z + 0.8z − 0.1)(z + f ) + g z
F ( z)
G ( z)
and the controller is
u( k )
− BGF y(k) − 1.292q−−.1.)(q ++106.16 y(k)
(q 0 9
The output variance when using the minimum variance controller in a. is
given by
var( y) 1 + f 1
Since B ( z) has a root outside the unit circle we use the procedure in Theorem
12.3. Factor B ( z) as
B ( z) B + ( z) B − ( z)
with B −∗ ( q) monic, so
B + ( z)
B − ( z)
The Diophantine to solve is
q d−1 C ( q) B −∗( q)
q3 ( q
− 0.1)(q − 0.9)
A( q) F ( q) + B − ( z) G ( q)
− 1.7q + 0.8q − 0.1)( f q + f q + f )
+ (−0.9q + 1)( g q + g q + g )
( q3
This gives the system of equations
which has the solution
Equation (12.31) gives
u( k )
G ( q)
+ ( q) F ( q)
− −1.245q 0.7q.931q7210.136 y(k)
y( k)
The closed loop system is
Ay( k)
Bu( k) + C e( k)
y( k)
− BGF y(k)
+ C e( k)
y( k) + C e( k)
− e(k) qd−1C B −∗ e(k)
e( k)
q( q 0.9)
e( k)
Theorem 10.4 gives the output variance
Var y
Problem 12.33
y( k)
q( q 0.7)
u( k ) +
e( k)
( q 1)( q 0.7)
( q 1)( q 0.7)
Determine the controller from
(q
− 1)(q − 0.7)(q + r ) + (0.9q + 1)(s q + s )
q2 ( q
This gives the system of equations
with the solution
y( k)
e( k)
var y( k)
(1 + 0.526q−1) e( k)
(1 + 0.5262)σ 2
Compare to Example 12.9 p. 468
var y( k)