Digital Control Lecture Notes

Download as pdf or txt
Download as pdf or txt
You are on page 1of 110

ELC 514E STATE SPACE DESIGN AND DIGITAL

CONTROL

CONTENTS:

• STATE SPACE DESIGN


• LYAPUNOV STABILITY
• DIGITAL CONTROL

Instructor: Dr.-Ing. S. I. Kamau


Department of Electrical and Communications Engineering,
Moi University
Contents

1 STATE SPACE DESIGN 1


1.1 State feedback design - a review . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Observer Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.2.2 Derivation of the observer equations . . . . . . . . . . . . . . . . . . . . . 6
1.2.3 Determination of matrix Ke . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.2.4 Selection of observer poles . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.2.5 Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.3 Combined state feedback controller and observer . . . . . . . . . . . . . . . . . . 14
1.3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.3.2 The independence of state feedback and observer design . . . . . . . . . . 14
1.3.3 Transfer function of the controller-observer . . . . . . . . . . . . . . . . . 16
1.4 State feedback design for systems with external inputs . . . . . . . . . . . . . . . 18
1.4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
1.4.2 Case 1: System has an integrator . . . . . . . . . . . . . . . . . . . . . . . 18
1.4.3 Case 2: System has no integrator (integral control) . . . . . . . . . . . . . 23
1.5 Selection of poles for a good design . . . . . . . . . . . . . . . . . . . . . . . . . . 27
1.5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
1.5.2 Dominant second order . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
1.5.3 Prototype design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
1.5.4 Integral of Time Multiplied by Absolute Error (ITAE) prototype . . . . . . 29
1.5.5 The Bessel prototype . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
1.5.6 Example: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

2 LYAPUNOV STABILITY 34
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.2 Equilibrium points . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

i
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau ii

2.3 Stability in the sense of Lyapunov . . . . . . . . . . . . . . . . . . . . . . . . . . . 37


2.4 Linearization and local stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
2.5 Lyapunov functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
2.5.1 Positive define and negative definite functions . . . . . . . . . . . . . . . . 41
2.5.2 The quadratic form . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
2.5.3 Lyapunov functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
2.6 Lyapunov stability analysis for LTI systems . . . . . . . . . . . . . . . . . . . . . . 44

3 DIGITAL CONTROL 47
3.1 Introduction to digital control systems . . . . . . . . . . . . . . . . . . . . . . . . 47
3.1.1 System configuration and terminology . . . . . . . . . . . . . . . . . . . . 47
3.1.2 Advantages and disadvantages of digital control . . . . . . . . . . . . . . . 48
3.2 The z-transform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
3.2.1 Definition of the z-transform . . . . . . . . . . . . . . . . . . . . . . . . . 49
3.2.2 z-transforms of common functions . . . . . . . . . . . . . . . . . . . . . . 52
3.2.3 Theorems of z-transform . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
3.2.4 Inverse z-transform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
3.3 Difference equations and pulse transfer functions . . . . . . . . . . . . . . . . . . 59
3.3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
3.3.2 Solving difference equations by recursion . . . . . . . . . . . . . . . . . . 60
3.3.3 Pulse transfer functions and weighting sequence . . . . . . . . . . . . . . 61
3.4 Mapping between the s-plane and the z-plane . . . . . . . . . . . . . . . . . . . . 62
3.4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
3.4.2 Mapping of the jω axis to the z-plane . . . . . . . . . . . . . . . . . . . . 63
3.4.3 Lines of constant attenuation . . . . . . . . . . . . . . . . . . . . . . . . . 63
3.4.4 Lines of constant frequency . . . . . . . . . . . . . . . . . . . . . . . . . . 64
3.4.5 Important note . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
3.5 Data holds . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
3.5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
3.5.2 The Zero-Order-Hold (ZOH) . . . . . . . . . . . . . . . . . . . . . . . . . . 67
1 − e−T s
 
3.5.3 z-transform of functions involving the term . . . . . . . . . . 69
s
3.6 Block diagram analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
3.6.1 Open-loop systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
3.6.2 Closed-loop systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau iii

3.6.3 Closed-loop systems with a Zero-Order-Hold . . . . . . . . . . . . . . . . . 76


3.7 Reconstruction of signals from sampled signals . . . . . . . . . . . . . . . . . . . 79
3.7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
3.7.2 Shannon’s Sampling Theorem . . . . . . . . . . . . . . . . . . . . . . . . . 79
3.7.3 Proving Shannon’s theorem . . . . . . . . . . . . . . . . . . . . . . . . . . 80
3.7.4 The ideal low pass filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
3.7.5 Frequency response of the Zero-Order-Hold . . . . . . . . . . . . . . . . . 83
3.7.6 Aliasing or frequency folding . . . . . . . . . . . . . . . . . . . . . . . . . 84
3.8 Discretization of controllers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
3.8.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
3.8.2 Discretization methods based on numerical integration . . . . . . . . . . . 87
3.8.3 Impulse and Step invariance methods . . . . . . . . . . . . . . . . . . . . 91
3.8.4 The matched pole-zero mapping method . . . . . . . . . . . . . . . . . . . 92
3.9 State space analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
3.9.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
3.9.2 Canonical forms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
3.9.3 Solving discrete-time state space equations . . . . . . . . . . . . . . . . . 99
3.9.4 Obtaining pulse transfer function from state space description . . . . . . . 102
3.9.5 Discretization of continuous-time state space equations . . . . . . . . . . . 103
Chapter 1

STATE SPACE DESIGN

1.1 State feedback design - a review

In the previous course ELC 412, we learnt that if we can measure the n internal state
variables x1 , x2 · · · xn of a system, then it is possible to feed back each state variable
xi through a gain ki to the input of the system. By choosing the feedback gains K =
[k1 k2 · · · kn ], we can arbitrarily locate the poles of the closed-loop system. However,
the system must be completely state controllable, and all the states must be accessible for
feedback.
It is assumed without loss of generality that the system is a Single Input Single Output
(SISO) system with a state space representation:

ẋ(t) = Ax(t) + Bu(t) (1.1a)


y(t) = Cx(t). (1.1b)

The control law is


u(t) = −Kx(t) (1.2)
(it is assumed that the closed-loop system has no external input.)
Substituting (1.2) into (1.1a) gives

ẋ(t) = Ax(t) + B (−Kx(t)) = (A − BK) x(t). (1.3)

The state matrix of the closed-loop system is therefore (A − BK). The block diagram of
the closed-loop system is shown in Figure 1.1.
When designing by pole placement, it is important to keep in mind that the control effort
required is proportional to how far the system poles are moved by the state feedback
- the further the poles move, the greater the control effort required (the K matrix will
have high gains).

1
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 2

rag replacements
ẋ(t) R x(t) y(t)
B C
+
+

A
SYSTEM

u(t) = −Kx(t)
−K STATE
FEEDBACK

Figure 1.1: State feedback - case when there is no external input

There are three methods of determining the feedback gain matrix K and they are out-
lined below without proofs – a reader interested in the proofs can refer to the handout
for ELC 412.

Method 1: The characteristic equation of system (1.1) is given by

sI − A = sn + an−1 sn−1 + an−2 sn−2 + · · · + a1 s + a0 = 0. (1.4)

Suppose the desired closed-loop poles are µ1 , µ2 , . . . µn . Then the desired char-
acteristic equation is

(s − µ1 ) (s − µ2 ) . . . (s − µn ) = sn +αn−1 sn−1 +αn−2 sn−2 +· · ·+α1 s+α0 = 0. (1.5)

Then
h i −1
K = α0 − a 0 α1 − a 1 α2 − a 2 ··· αn−1 − an−1 MW ,
(1.6)
where M is the controllability matrix:

M = B AB A2 B An−1 B
 
··· (1.7)

and  
a1 a2 a3 · · · an−1 1
a2 a3 a4 ··· 1 0
 
 

W= .. .. .. .. .. ..
. (1.8)
 . . . . . .
 an−1 1 0 ··· 0 0 
 

1 0 0 ··· 0 0

Method 2: Ackermann’s formula


−1
A2 B An−1 B

K = [0 0 ··· 0 1] B AB ··· φ (A) ,
(1.9)
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 3

where
φ (A) = An + αn−1 An−1 + αn−2 An−2 + · · · + α1 A + α0 I
(compare this with the desired characteristic equation (1.5)). This method is
known as Ackermann’s formula.

Method 3: Direct substitution method From equation (1.3), the characteristic equa-
tion of the closed loop system is

sI − (A − BK) = sI − A + BK = 0. (1.10)

For the direct substitution method, the above equation is equated to the desired
characteristic equation (1.5), i.e.

sI − A + BK = sn + αn−1 sn−1 + · · · + α1 s + α0 = 0. (1.11)

This is known as the Direct Substitution Method. The elements of the matrix K
are obtained by equating coefficients of equal powers of s in (1.11) above. Note
that this method becomes tedious if the number of states n is large - it is mostly
used for cases where n ≤ 3.

Example
Given is a continuous-time system:
" # " #" # " #
ẋ1 (t) 0 1 x1 (t) 0
= 2
+ u(t).
ẋ2 (t) −ω0 0 x2 (t) 1

This system has open-loop poles at s = ±jω0 . We would like to use state feedback law
u = −Kx to locate both of the closed-loop poles at s = −2ω0 (i.e. the closed-loop
system should have double poles at s = −2ω0 ). Determine K.
Solution:
The system is first tested for controllability. The controllability matrix is
" #
0 1
M = [B AB] = .
1 0

|M| = −1 6= 0,
hence the system is completely state controllable and design by pole placement is pos-
sible.

(i) Determination of K – Method 1:

|sI − A| = s2 + ω02 = s2 + a1 s + a0 ,
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 4

hence a1 = 0 and a0 = ω02 .


The controllability matrix M is obtained above as
" #
0 1
M= .
1 0

For a second order system,


" # " #
a1 1 0 1
W= =
1 0 1 0
(compare with (1.8) above.)
" #" # " #
0 1 0 1 1 0
MW = = = I,
1 0 1 0 0 1
hence MW = (MW)−1 = I.
The desired characteristic equation is

(s + 2ω0 ) (s + 2ω0 ) = s2 + 4ω0 s + 4ω02 = s2 + α1 s + α0 ,

hence α1 = 4ω0 and α0 = 4ω02 .


Substituting all the values in (1.6) gives

α1 − a1 ] (MW)−1 = 3ω02
 
K = [α0 − a0 4ω0 .

(ii) Determination of K – Method 2: For a second-order system, Ackermann’s formula


is given by
K = [0 1] [B AB]−1 φ (A) , (1.12)
where

φ (A) = A2 + α1 A + α0 I = A2 + 4ω0 A + 4ω02 I


" #" # " # " #
0 1 0 1 0 1 1 0
= 2 2
+ 4ω0 2
+ 4ω02
−ω0 0 −ω0 0 −ω0 0 0 1
" #
3ω02 4ω0
= .
−4ω03 3ω02
From the working above,
" #−1 " #
0 1 0 1
[B AB]−1 = = .
1 0 1 0
It follows that
" #" #
0 1 3ω02 4ω0
= 3ω02
 
K = [0 1] 4ω0 .
1 0 −4ω03 3ω02
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 5

(iii) Determination of K – Method 3: Let K = [k1 k2 ]. Then


" #
s −1
[sI − A + BK] = .
ω02 + k1 s + k2

s −1
= s2 + sk2 + ω02 + k1 = s2 + 4ω0 s + 4ω02 .

sI − A + BK =
ω02 + k1 s + k2
Equating coefficients of equal powers of s gives k1 = 3ω02 and k2 = 4ω0 , hence

K = 3ω02
 
4ω0 .

1.2 Observer Design

1.2.1 Introduction

For design by state variable feedback (design by pole placement), it is assumed that all
the states are available for feedback. However, there are situations where not all states
are available for feedback for the following reasons:

• some states may not be measurable

• it may be too expensive to measure all the states

• the state measurements may be noisy.

When some or all the states are not available, we need to estimate the unavailable state
variables. This is done using a state observer (also known as a state estimator). A state
observer is an electronic circuit or a computer program which estimates the states of
the system from the measured input and output signals.
An observer which estimates all the state variables of a system regardless of whether
some of the state variables are available for direct measurement is known as a full-order
state observer.
An observer which estimates fewer than n state variables, where n is the order of the sys-
tem is called a reduced order state observer. If the reduced order observer only estimates
those states which cannot be measured, it is called a minimum-order state observer.
A state observer can only be designed if the system is observable. The design also
requires an accurate model of the system.
A block diagram of a system connected to a full-order state observer is shown in Fig-
ure 1.2.
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 6

SYSTEM
ẋ(t) x(t) y(t)
C
R
B
+
+

replacements
A

u(t)
x̂(t)

−+
˙
x̂(t) x̂(t) ŷ(t)
C
R
B
++ +
+

−K 
Ke y(t) − ŷ(t) y(t) − ŷ(t)
Ke
FULL ORDER OBSERVER

Figure 1.2: A system connected to a full-order state observer

Note that the full-order observer has two inputs, u(t) and y(t). The feedback term
Ke (y(t) − ŷ(t)) is needed to ensure that
lim x̂(t) −→ x(t), (1.13a)
t→∞
lim ŷ(t) −→ y(t), (1.13b)
t→∞

i.e.
 
lim x̂(t) − x(t) = 0, (1.14a)
t→∞
 
lim ŷ(t) − y(t) = 0. (1.14b)
t→∞

The convergence is achieved by an appropriate selection of the feedback gain matrix


Ke , and this is discussed in more detail in the next section.

1.2.2 Derivation of the observer equations

From Figure 1.2, we can write the state equation of the full-order observer (the lower
part of the block diagram) as:
 
˙
x̂(t) = Ax̂(t) + Bu(t) + Ke y(t) − ŷ(t) . (1.15)
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 7

But y(t) = Cx(t) and ŷ(t) = Cx̂(t), hence

 
˙
x̂(t) = Ax̂(t) + Bu(t) + Ke Cx(t) − Cx̂(t) . (1.16)
For the system (the upper part of the block diagram), the state equation is given by

ẋ(t) = Ax(t) + Bu(t). (1.17)

Subtracting equation (1.16) from equation (1.17) gives the observer error equation, i.e.
the equation which gives the error between the actual state x(t) and the estimated state
x̂(t):
 
˙
ẋ(t) − x̂(t) = Ax(t) + Bu(t) − Ax̂(t) − Bu(t) − Ke Cx(t) − Cx̂(t)
      
= A x(t) − x̂(t) − Ke C x(t) − x̂(t) = A − Ke C x(t) − x̂(t) . (1.18)

Defining the error vector e(t) = (x(t) − x̂(t)), then the above equation can be written as
 
ė(t) = A − Ke C e(t). (1.19)

From (1.19) above, it can be seen that the behaviour of the error vector depends on
the eigenvalues of the matrix (A − Ke C). If the matrix (A − Ke C) is a stable matrix,
then the error vector e(t) will converge to zero for any initial error vector e(0), i.e. x̂(t)
will converge to x(t) regardless of the x̂(0) and x(0). The eigenvalues of the matrix
(A − Ke C) also determine how fast the estimated state x̂(t) converge to the state x(t).
The characteristic equation for the error vector system is

sI − A + Ke C = sn + αn−1 sn−1 + αn−2 sn−2 + · · · + α1 s + α0 = 0. (1.20)

Now recall dual systems discussed in ELC 412. For a system given by the state space
representation

ẋ(t) = Ax(t) + Bu(t) (1.21a)


y(t) = Cx(t), (1.21b)

the dual system is given by

ż(t) = AT z(t) + CT v(t) (1.22a)


w(t) = BT z(t), (1.22b)

where z(t) is the state vector of the dual system, v(t) the input and w(t) the output.
Consider system (1.21). For the pole placement problem, the control law is given by
equation (1.2), reproduced here for convenience:

u(t) = −Kx(t). (1.23)


ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 8

The closed-loop system is given by


   
ẋ(t) = Ax(t) + B − Kx(t) = A − BK x(t). (1.24)
The state matrix of the closed-loop system is (A − BK) and this matrix determines the
eigenvalues of the closed-loop system. The characteristic equation of the closed-loop
system is
sI − A + BK = 0. (1.25)

Now consider the dual system (1.22). For the pole-placement problem, the correspond-
ing control law is:
v(t) = −Kz(t). (1.26)
The closed-loop system is therefore given by
   
ż(t) = AT z(t) + CT − Kz(t) = A − CT K z(t). (1.27)

The state matrix of the closed-loop system is (AT −CT K) and this matrix determines the
eigenvalues of the closed-loop dual system. The characteristic equation of the closed-
loop dual system is
sI − AT + CT K = 0. (1.28)

From matrix theory, the eigenvalues of a matrix F are the same as the eigenvalues of the
transpose of that matrix FT . It follows that the eigenvalues of the matrix (AT − CT K)
are the same as those of (AT − CT K)T .
But  T
AT − C T K = A − KT C.
From the above equation, the characteristic equation of the closed-loop dual system
given in (1.28) can be written as
sI − A + KT C = 0. (1.29)

Comparing equations (1.20) and (1.29), it can be seen that the matrix Ke can
be determined by pole placement approach for the dual of the original system
then using Ke = KT .
For pole placement of the dual system, the dual system (1.22) must be completely state
controllable. The controllability matrix of the dual system (1.22) is
h i
T T T T 2 T n−1
 T  T
N= C A C A C ··· A C . (1.30)
We recognize (1.30) as the observability matrix of the original system (1.21), which
leads us to the necessary condition for observer design:
A system (1.21) must be observable for it to be possible to design an observer.
Conversely, if the system (1.21) is observable, then it is possible to choose a
matrix Ke such that the matrix (A − Ke C) has the desired eigenvalues.
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 9

1.2.3 Determination of matrix Ke

We now look at the methods of determining the matrix Ke . There are three methods
which are derived from the methods given in Section 1.1.

Determination of Ke : Method 1

This method is based on equation (1.6).


Suppose we would like to design an observer for the system

ẋ(t) = Ax(t) + Bu(t) (1.31a)


y(t) = Cx(t). (1.31b)

As explained above, observer design (determination of Ke ) can be carried out by pole


placement approach for the dual of the original system then using Ke = KT .
The dual of the above system is

ż(t) = AT z(t) + CT v(t) (1.32a)


w(t) = BT z(t). (1.32b)

The characteristic equation of the dual system is

sI − AT = sI − A = sn + an−1 sn−1 + an−2 sn−2 + · · · + a1 s + a0 = 0, (1.33)

while the controllability matrix of the dual system is


h 2 i
T n−1
N = CT A T CT A T CT T

··· A C . (1.34)

Since the characteristic equations of the original system and its dual are identical, then
 
a1 a2 a3 · · · an−1 1
 a2 a3 a4 · · · 1 0 
 
 . .. .. . . .. .. 
W= .. . . . . .  (1.35)
 
 an−1 1 0 · · · 0 0 
 

1 0 0 ··· 0 0

(compare with equation (1.8)).


Suppose the desired characteristic equation for the observer is

sI − A + Ke C = sn + αn−1 sn−1 + αn−2 sn−2 + · · · + α1 s + α0 = 0. (1.36)

Then analogous to equation (1.6) we can write


h i
K = α0 − a 0 α1 − a 1 α2 − a 2 ··· αn−1 − an−1 (NW)−1 . (1.37)
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 10

But Ke = KT hence
h i T
−1
Ke = α0 − a 0 α1 − a 1 α2 − a 2 ··· αn−1 − an−1 (NW) , (1.38)

or
 
α0 − a 0
α1 − a 1
 
−1 



Ke = WNT  α2 − a 2 . (1.39)
..
 
.
 
 
αn−1 − an−1

Determination of Ke : Method 2

This method is based on equation (1.9) – Ackermann’s formula.


For pole placement of the dual system, we can write analogous to equation (1.9)
h i−1
T T T T 2 T T n−1 T
φ AT .
  
K = [0 0 · · · 0 1] C A C A C ··· A C (1.40)

But Ke = KT hence
h h 2 n−1 T i−1  iT
Ke = [0 0 · · · 0 1] CT AT CT A T CT · · · A T C φ AT

T h T i−1 T
T 2 T n−1
T T T
[0 0 · · · 0 1]T .
 T  T
= φ A C A C A C ··· A C

(1.41)

T
From matrix theory, φ AT = φ (A), hence

 −1  
C 0
CA 0
   
   
CA2
   
Ke = φ (A)    0 , (1.42)
.. ..
   
. .
   
   
CAn−1 1

where
φ (A) = An + αn−1 An−1 + αn−2 An−2 + · · · + α1 A + α0 I
(compare this with the desired characteristic equation (1.36)). This method is known
as Ackermann’s formula for observer design.
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 11

Determination of Ke : Method 3

Suppose the desired observer poles are s = µ1 , s = µ2 , . . . s = µn . Then the desired


characteristic equation is

sI − A + Ke C = (s − µ1 )(s − µ2 ) · · · (s − µn ) = sn + αn−1 sn−1 + · · · + α1 s + α0 = 0.


(1.43)
This is known as the Direct Substitution Method for observer design. The elements
of the matrix Ke are obtained by equating coefficients of equal powers of s in (1.43)
above. Note that this method becomes tedious if the number of states n is large - it is
mostly used for cases where n ≤ 3.

1.2.4 Selection of observer poles

The location of the observer poles is a trade-off between the time it takes for the esti-
mated state x̂(t) to converge with the actual state x(t), the control effort required and
the effects of noise and disturbances.
As a rule of thumb, the observer poles (eigenvalues of the matrix (A−Ke C)) are chosen
to be a factor of 2 to 6 times faster than the controller poles (eigenvalues of the matrix
(A − BK) from pole placement).
If the observer poles are chosen faster than 6 times the controller poles, the estimated
state x̂(t) converges with the actual state x(t) in a short time but a large control effort
is required and the effects of noise and disturbances are highly amplified.
If the observer poles are chosen slower than 2 times the controller poles, the estimated
state x̂(t) takes a long time to converge with the actual state x(t), but the control effort
required is low and the effect of noise and disturbances is minimal.

1.2.5 Example

Given a continuous-time system:


      
ẋ1 (t) 1 2 0 x1 (t) 2
 ẋ2 (t)  =  3 −1 1   x2 (t)  +  1  u(t)
      

ẋ3 (t) 0 2 0 x3 (t) 1


 
h i x1 (t)
y(t) = 0 0 1  x2 (t)  ,
 

x3 (t)
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 12

suppose that the states of the system are not available for feedback. Design a full order
state observer with the eigenvalues of the matrix (A − Ke C) are located at s = −4,
s = −3 ± j.
Solution:
Observability matrix    
C 0 0 1
N =  CA  =  0 2 0  .
   

CA2 6 −2 2
N = 12 6= 0
hence the system is observable and observer design is possible.
Method 1:
Since the given system is a third order system,
 
 −1 α0 − a0
Ke = WNT  α1 − a 1  .
 

α2 − a 2
The desired characteristic equation is
(s + 4)(s + 3 + j)(s + 3 − j) = s3 + 10s2 + 34s + 40 = s3 + α2 s2 + α1 s + α0 = 0,
hence α2 = 10, α1 = 34 and α0 = 40.
Also
sI − A = s3 − 9s + 2 = s3 + a2 s2 + a1 s + a0 ,
hence a2 = 0, a1 = −9 and a0 = 2.

 
h i 0 0 6
2
N = CT A T CT AT CT =  0 2 −2  .
 

1 0 2
For a third order system,
   
a1 a2 1 −9 0 1
W =  a2 1 0  =  0 1 0  .
   

1 0 0 1 0 0
From the above working,
    −1  
 −1 α 0 − a 0 −9 0 1 0 0 1 38
T
Ke = WN  α1 − a 1  =   0 1 0   0 2 0   43 
      

α2 − a 2 1 0 0 6 −2 2 10
    
1 1 7
6 6 6
38 25 61
=  0 12 0   43  =  21 12  .
    

0 0 1 10 10
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 13

Method 2:
Ackermann’s formula for observer design for a third order system is
 −1  
C 0
Ke = φ (A)  CA   0 
   

CA2 1
where
φ (A) = A3 + α2 A2 + α1 A + α0 I.
From the working above, α2 = 10, α1 = 34 and α0 = 40, hence
 
151 86 20
φ (A) = A3 + 10A2 + 34A + 40I =  129 85 33  .
 

60 66 58
Also,    
C 0 0 1
 CA  =  0 2 0  .
   

CA2 6 −2 2
It follows that
  −1  
151 86 20 0 0 1 0
Ke =  129 85 33  0 2 0   0 
    

60 66 58 6 −2 2 1
     
151 86 20 − 31 16 1
6
0 25 1
6
=  129 85 33   0 12 0   0  =  21 12  .
     

60 66 58 1 0 0 1 10

Method 3:
For the direct substitution method for observer design,

sI − A + Ke C = (s + 4)(s + 3 + j)(s + 3 − j) = s3 + 10s2 + 34s + 40.

Let  
k e1
Ke =  k e 2  .
 

k e3
Then
     
s 0 0 1 2 0 0 0 k e1
sI − A + Ke C =  0 s 0  −  3 −1 1  +  0 0 k e2 
     

0 0 s 0 2 0 0 0 k e3
 
s − 1 −2 k e1
=  −3 s + 1 −1 + ke2 .
 

0 −2 s + k e3
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 14

sI − A + Ke C = s3 + (ke3 ) s2 + (2ke2 −9 ) s + (−7ke3 − 2ke2 + 6ke1 + 2)


= s3 + 10s2 + 34s + 40 = 0.

Equating coefficients of equal powers of s in the above equation gives ke1 = 25 16 ,


ke2 = 21 12 and ke3 = 10 hence
 
25 16
Ke =  21 12  .
 

10

1.3 Combined state feedback controller and observer

1.3.1 Introduction

A combined state feedback controller and observer is shown in Figure 1.3. The dif-
ference between Figures 1.2 and 1.3 is that in Figure 1.2, the estimated state is not
connected to anything but in Figure 1.3, the estimated state is used for state feedback
and there is no external input.

1.3.2 The independence of state feedback and observer design

From Figure 1.3, we can write for the system (the upper part of the block diagram)

ẋ(t) = Ax(t) + Bu(t) (1.44a)


y(t) = Cx(t). (1.44b)

This time, the control law is


u(t) = −Kx̂(t). (1.45)
Substituting (1.45) into (1.44a) gives

ẋ(t) = Ax(t) − BKx̂(t) = Ax(t) − BKx̂(t) + BKx(t) − BKx(t)


   
= A − BK x(t) + BK x(t) − x̂(t) . (1.46)

Defining the error vector e(t) = (x(t) − x̂(t)), then the above equation can be written as
 
ẋ(t) = A − BK x(t) + BKe(t). (1.47)

For the differential equation for the observer (the lower part of Figure 1.3), recall the
observer error equation (1.19), reproduced here for convenience:
 
ė(t) = A − Ke C e(t). (1.48)
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 15

SYSTEM
ẋ(t) x(t) y(t)
C
R
B
+
+

replacements A

u(t) = −Kx̂(t)
−K

˙
x̂(t) x̂(t) ŷ(t) +

C
R
B
+ +
+ +

y(t) − ŷ(t)
 Ke
Ke y(t) − ŷ(t)
FULL ORDER OBSERVER

Figure 1.3: Combined state feedback controller and observer

Equations (1.47) and (1.48) can be written in vector-matrix form as


" # " #" #
ẋ(t) A − BK BK x(t)
= , (1.49)
ė(t) 0 A − Ke C e(t)

which is the differential equation of the combined controller-observer. From (1.49)


above, the characteristic equation is

sI − A + BK −BK
= 0. (1.50)
0 sI − A + Ke C

From matrix theory,


P Q
= P R. (1.51)
0 R
It follows that the characteristic equation of the combined controller-observer is

sI − A + BK sI − A + Ke C = 0. (1.52)

The equation |sI − A + BK| is the characteristic equation of the closed-loop system
with state feedback (compare with equation (1.11)). On the other hand, |sI − A + Ke C|
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 16

is the characteristic equation of the observer (compare with equation (1.20)). This leads
us to the following conclusion:
Pole placement and observer design are independent of each
other, i.e. K and Ke are determined independently.
Note also that the poles of the closed-loop system shown in Figure 1.3 are the poles of
the controller (eigenvalues of (A − BK)) and the poles of the observer (eigenvalues of
(A − Ke C)). These poles can also be obtained by solving equation (1.52).

1.3.3 Transfer function of the controller-observer

From Figure 1.3, the differential equation for the observer (the lower part of the block
diagram) is  
˙
x̂(t) = Ax̂(t) + Bu(t) + Ke y(t) − ŷ(t) . (1.53)
But
ŷ(t) = Cx̂(t) (1.54)
and
u(t) = −Kx̂(t). (1.55)
Substituting (1.54) and (1.55) into (1.53) gives
 
˙
x̂(t) = A − BK − Ke C x̂(t) + Ke y(t). (1.56)

To obtain the transfer function of the controller-observer, we take the Laplace trans-
forms of (1.56) and (1.55) and eliminate the state variable.
The Laplace transform of (1.56) is
 
sX̂(s) = A − BK − Ke C X̂(s) + Ke Y (s), (1.57)

which can be rewritten as


 
sI − A + BK + Ke C X̂(s) = Ke Y (s). (1.58)

We make X̂(s) the subject of the formula by pre-multiplying both sides of (1.58) by
(sI − A + BK + Ke C)−1 , resulting in
 −1
X̂(s) = sI − A + BK + Ke C Ke Y (s). (1.59)

Taking the Laplace transform of (1.55) gives

U (s) = −KX̂(s). (1.60)


ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 17

Substituting (1.59) into (1.60) gives


 −1
U (s) = −K sI − A + BK + Ke C Ke Y (s). (1.61)

The transfer function of the controller-observer is obtained from equation (1.61) as

U (s)  −1
= K sI − A + BK + Ke C Ke . (1.62)
−Y (s)

This can be represented in a block diagram as shown in Figure 1.4.


replacements
Controller−Observer System
R(s) = 0 −Y (s)  −1 U (s) Y (s)
K sI−A+BK+KeC Ke C(sI − A)−1 B + D
+

Figure 1.4: Alternative representation of the controller-observer

Exercise:
Given a continuous-time system
" # " #" # " #
ẋ1 (t) 0 1 x1 (t) 0
= + u(t)
ẋ2 (t) 0 0 x2 (t) 1
" #
h i x (t)
1
y(t) = 1 0 + [0]u(t),
x2 (t)

(i) Determine K such that the poles of the state feedback controller (roots of the
equation |sI − A + BK|) are located at s = −2 ± j2.

(ii) Determine Ke such that the observer poles (roots of the equation
|sI − A + Ke C|) are located at s = −4 ± j4.

(iii) Determine the transfer function of the controller-observer.

Solutions:

(i) K = [8 4]
" #
8
(ii) Ke =
32
 −1 64(3s + 4)
(iii) K sI − A + BK + Ke C Ke = .
s2+ 12s + 72
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 18

1.4 State feedback design for systems with external in-


puts

1.4.1 Introduction

So far, we have considered systems which do not have external inputs. Such systems
are known as regulator systems. A regulator control system has no external input and
the purpose of the system is to return all state variables to zero when the states have
been perturbed.
However, many systems require that the output track an input signal. Such a control
signal, in which the output must follow command signals, is known as a servo system.
In a servo system it is generally required that the system have one or more integrators
within the closed-loop. If the plant to be controlled does not have an integrator, it is
usually necessary to add one or more integrators within the loop to eliminate/minimize
steady state errors.

1.4.2 Case 1: System has an integrator


replacements
Suppose the system shown in Figure 1.1 was modified to incorporate an external input.
The resulting system is shown in Figure 1.5.

r(t) u(t) ẋ(t) R x(t) y(t)


kr B C
+ +
EXTERNAL +
INPUT

A
SYSTEM

STATE
K FEEDBACK

Figure 1.5: State feedback - case when there is an external input

In this case, we assume that the system has at least one integrator. Since we are feeding
back the states and not the output, the input gain required for the output to track
the input is not necessarily equal to 1 so the input gain is set to kr . To simplify the
derivations, we assume that the external input is a unit step.
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 19

The state equation of the system (enclosed in the dotted lines) is

ẋ(t) = Ax(t) + Bu(t) (1.63)

while the control law is

u(t) = −Kx(t) + kr r(t) = −k1 x1 (t) − k2 x2 (t) − · · · − kn xn (t) + kr r(t). (1.64)

The system shown in Figure 1.5 is an n-th order system, but equation (1.64) has (n + 1)
unknowns k1 , k2 , . . . kn and kr , so we cannot uniquely determine all the unknowns in
equation (1.64).
To be able to solve for all the gains, we reconfigure the closed-loop system so that it
behaves like a regulator. From the design point of view, we desire

lim y(t) = y(∞) = r(t) = 1, (1.65)


t→∞

hence  
lim y(t) − r(t) = 0. (1.66)
t→∞

Equations (1.65) and (1.66) are used to help us convert the system shown in Figure 1.5
to a system which behaves like a regulator.
To convert the system into a regulator system, we let the output y(t) be equal to one of
the states, usually state x1 (t), and we also let kr be equal to one of the feedback gains
(the gain corresponding to the state equal to the output).
Let

x1 (t) = y(t) (1.67a)


kr = k 1 . (1.67b)

Substituting (1.67) into (1.64) gives


 
u(t) = −k1 y(t)−k2 x2 (t)−· · ·−kn xn (t)+k1 r(t) = −k1 y(t)−r(t) −k2 x2 (t)−· · ·−kn xn (t).
(1.68)
From (1.68), we can see that

lim u(t) = u(∞) = 0, (1.69)


t→∞

hence the system behaves like a regulator system, and the methods used to determine
the feedback gains for regulator systems can now be applied to this system.
Equation (1.68) can be rewritten as

u(t) = −Kx(t) + k1 r(t). (1.70)

The new configuration of the closed-loop system is shown in Figure 1.6.


ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 20

PLANT
replacements x1 (t)
r(t) u(t) x2 (t)
y = x1 (t)
k1 ẋ(t) = Ax(t) + Bu(t) x3 (t) C
+ +
xn (t)

k2

k3

kn

Figure 1.6: State feedback where there is an external input – alternative configuration

From Figure 1.6,


ẋ(t) = Ax(t) + Bu(t) (1.71)
and the control law is
u(t) = −Kx(t) + k1 r(t). (1.72)
Substituting (1.72) into (1.71) gives
 
ẋ(t) = A − BK x(t) + Bk1 r(t), (1.73)

which is the state equation of the closed-loop system. The characteristic equation of the
closed-loop system is
sI − A + BK = 0 (1.74)
(compare this with the characteristic equation of the regulator system, equation (1.10)).

If the system (1.63) is completely state controllable, then the matrix K in


(1.73) above can be determined using the pole placement methods reviewed
in Section 1.1.
Example
A plant is given by the transfer function
Y (s) 1
= 2.
U (s) s

For the system configuration shown in Figure 1.6, determine the feedback gain matrix
K so that the control law u(t) = −Kx(t) + k1 r(t), where r(t) is a unit step input, results
in a closed-loop system with poles located at s = −4 ± j4.
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 21

Solution

Y (s) 1
= 2
U (s) s
Hence
s2 Y (s) = U (s).

Converting the above equation into time domain gives

ÿ(t) = u(t).

Since this is a system with an external input, the states are selected as x1 (t) = y(t) and
x2 (t) = ẏ(t).
The state equations can then be written as:

ẋ1 (t) = ẏ(t) = x2 (t)


ẋ2 (t) = ÿ(t) = u(t)

The above equations can be written in vector-matrix notation as:


" # " #" # " #
ẋ1 (t) 0 1 x1 (t) 0
= + u(t) = Ax(t) + Bu(t)
ẋ2 (t) 0 0 x2 (t) 1
" #
x1 (t)
y(t) = [1 0] = Cx(t)
x2 (t)

The controllability matrix is


" #
h i 0 1
M= B AB = .
1 0

|M| = −1 6= 0,
hence the system is completely state controllable.
Since no pole placement method is specified, we can use any of the three methods in
Section 1.1 – in this case, Ackermann’s formula is chosen.
The desired characteristic equation is

(s + 4 + j4)(s + 4 − j4) = s2 + 8s + 32 = 0 = s2 + α1 s + α0 = 0,

hence
" #" # " # " # " #
0 1 0 1 0 1 1 0 32 8
φ(A) = A2 + 8A + 32I = +8 + 32 = .
0 0 0 0 0 0 0 1 0 32
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 22

" #−1 " #


h i h i 0 1 32 8 h i
K= k1 k2 = 0 1 = 32 8 .
1 0 0 32

A block diagram of the closed-loop system is shown in Figure 1.7 (compare this with
replacementsthe block diagram shown in Figure 1.6).

PLANT
r(t) u(t) 1 x2 (t) 1 y(t) = x1 (t)
32
+ + s s
k1

8
k2

Figure 1.7: Example

Note that the state space representation of the closed-loop system is


" #" # " #
  0 1 x1 (t) 0
ẋ(t) = A − BK x(t) + Bk1 r(t) = + r(t)
−32 −8 x2 (t) 32
" #
x1 (t)
y(t) = Cx(t) = [1 0] ,
x2 (t)

and the transfer function of the closed-loop system can be obtained using the equation

Y (s)
= C (sI − A + BK)−1 Bk1 .
R(s)
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 23

Exercise:
A plant is given by the transfer function

Y (s) 1
= .
U (s) s(s + 1)(s + 2)

For the system configuration shown in Figure 1.6, use the direct substitution method
to determine the feedback gain matrix K so that the control law

u(t) = −Kx(t) + k1 r(t),

where r(t) is a unit step input, results in a closed-loop system with poles located at

s = −2 ± j2 3, s = −10, and draw a block diagram of the resulting system.
Solution:
h i
K= 160 54 11 .

1.4.3 Case 2: System has no integrator (integral control)

State space design does not produce integral control unless special steps are taken in
the design process. If the system does not have an integrator (i.e. is type 0), the system
will have finite steady-state errors for step inputs and infinite steady-state errors for
ramp and parabolic inputs.
replacements
To simplify the mathematical analysis, we restrict our system to step inputs. To eliminate
steady-state errors for step inputs, we need to add an integrator to a type 0 system to
convert it to a type 1 system. The integrator is added to the closed-loop system as shown
in Figure 1.8.
SYSTEM
r(t) ζ̇(t) R ζ(t) u(t) ẋ(t) x(t) y(t) = x1 (t)
kI
R
+ + B +
C
INTEGRATOR +

STATE
K FEEDBACK

Figure 1.8: State feedback - Integral control

The closed-loop system shown in Figure 1.8 is of the order (n + 1). From the figure, we
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 24

can write the following equations:

ẋ(t) = Ax(t) + Bu(t) (1.75a)


y(t) = Cx(t) (1.75b)
u(t) = −Kx(t) + kI ζ(t) (1.75c)
ζ̇(t) = r(t) − y(t) = r(t) − Cx(t). (1.75d)

Equations (1.75a) and (1.75d) can be written in vector-matrix notation:


" # " #" # " # " #
ẋ(t) A 0 x(t) B 0
= + u(t) + r(t), (1.76)
ζ̇(t) −C 0 ζ(t) 0 1

while the control law is given by equation (1.75c), i.e.

u(t) = −Kx(t) + kI ζ(t). (1.77)

At t = ∞,
" # " #" # " # " #
ẋ(∞) A 0 x(∞) B 0
= + u(∞) + r(∞), (1.78)
ζ̇(∞) −C 0 ζ(∞) 0 1

and
u(∞) = −Kx(∞) + kI ζ(∞). (1.79)
Subtracting (1.78) from (1.76), keeping in mind that r(t) = r(∞) = 1, gives
    " #    " #
ẋ(t) − ẋ(∞) A 0 x(t) − x(∞) B  
   =    + u(t) − u(∞) . (1.80)
ζ̇(t) − ζ̇(∞) −C 0 ζ(t) − ζ(∞) 0

Similarly, subtracting (1.79) from (1.77) gives


     
u(t) − u(∞) = −K x(t) − x(∞) + kI ζ(t) − ζ(∞) . (1.81)

Let

x(t) − x(∞) = xe (t) (1.82a)


ζ(t) − ζ(∞) = ζe (t) (1.82b)
u(t) − u(∞) = ue (t). (1.82c)

Substituting equations (1.82) into ((1.80)) and (1.81) gives


" # " #" # " #
ẋe (t) A 0 xe (t) B
= + ue (t), (1.83)
ζ̇e (t) −C 0 ζe (t) 0

and " #
h i xe (t)
ue (t) = −Kxe (t) + kI ζe (t) = − K −kI (1.84)
ζe (t)
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 25

respectively.

" # " # " #


A 0 B xe (t) h i
Letting  = , B̂ = , e(t) = and K̂ = K −kI ,
−C 0 0 ζe (t)

equations (1.83) and (1.84) can be rewritten as

ė(t) = Âe(t) + B̂ue (t) (1.85)

and
ue (t) = −K̂e(t) (1.86)
respectively.
Equations (1.85) and (1.86) are similar to the regulator equations (1.1a) and (1.2),
which leads us to the following conclusion:

If the system (1.85) is completely state controllable, then the matrix K̂ in


(1.86) above can be determined using the pole placement methods reviewed
in Section 1.1.
Note that since the closed-loop system shown in Figure 1.8 is an (n + 1)-th order system,
matrix K̂ is also an (n + 1) dimension vector, i.e.
h i h i
K̂ = K −kI = k1 k2 k3 · · · kn −kI .

Example
A plant is given by the transfer function

Y (s) 1
= .
U (s) s+3

Use Ackermann’s formula to design a type 1 servo system with closed-loop poles at

s = −2.5 ± j2.5 3.
Solution
The given system is a first order system, and after adding an integrator, the final system
should be a type 1 second order system which gives zero steady-state errors to step
inputs.
Since the system has an external input, the output should be chosen to be equal to one
of the states.
From the given transfer function above, we can write

(s + 3)Y (s) = U (s).


ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 26

Converting to time domain gives

ẏ(t) + 3y(t) = u(t).

This is a first order system so there is only one state, which we can name x(t).
Letting
x(t) = y(t),
then
ẋ(t) = ẏ(t) = −3y(t) + u(t) = −3x(t) + u(t).
The state space representation of the system is

ẋ(t) = −3x(t) + u(t)


y(t) = x(t)

hence A = −3, B = 1, C = 1 and D = 0.


It follows that
" # " # " # " #
A 0 −3 0 B 1 h i
 = = , B̂ = = , and K̂ = k1 −kI .
−C 0 −1 0 0 0

The controllability matrix is


" #
h i 1 −3
M̂ = B̂ ÂB̂ = .
0 1

M̂ = −1 6= 0,
hence the system is completely state controllable.
Also " #
1 −3
M̂−1 = .
0 −1
The desired characteristic equation is
 √  √ 
s + 2.5 + j2.5 3 s + 2.5 − j2.5 3 = s2 + 5s + 25.

Hence " #
19 0
φ(Â) = Â2 + 5Â + 25I = .
−2 25
From the working above,
" #" #
h i h i 1 −3 19 0 h i
K̂ = k1 −kI = 0 1 = 2 −25 ,
0 −1 −2 25
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 27

PLANT
replacements
−3

+
r(t) 1 u(t) + 1 y(t) = x(t)
25
+ s + s
kI

2
k1

Figure 1.9: Example

hence k1 = 2 and kI = 25.


A block diagram of the closed-loop system is shown in Figure 1.9 (compare this with
the block diagram shown in Figure 1.8).
Exercise:
A plant is given by the transfer function

Y (s) 1
= .
U (s) (s + 1)(s + 2)

Use the direct substitution method to design a type 1 servo system with closed-loop

poles at s = −10, s = −2 ± j2 3, and draw a block diagram of the resulting system.
Solution:
h i h i
K̂ = k1 k2 −kI = 54 11 −160 .

1.5 Selection of poles for a good design

1.5.1 Introduction

A system with state feedback control is superior to a PID controller since it allows the
placing of all roots of the closed-loop system at the desired locations.
The roots of the closed-loop system are selected to satisfy the performance specifica-
tions. It thus follows that the first task in state-feedback design is to translate all per-
formance specifications into desired root locations. The control law is then obtained by
feeding back all the state variables with constant gains. The values of the gains depend
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 28

on the root locations.


Some common specifications include:

 Stable operation with adequate stability margin

 Transient response specifications (rise time, settling time, overshoot, etc.)

 Frequency response specifications (phase margin, gain margin, bandwidth)

 Desired system type (type 0, type 1, type 2, etc.)

 Sensitivity of system to parameter variations (robustness).

When selecting poles, it is useful to remember that the control effort required is related
to how far the open-loop poles are moved by the state feedback – the further the poles
are moved, the higher the control effort required.
Two commonly-used methods of pole selection criteria are

 dominant second order

 prototype design

1.5.2 Dominant second order

In this case, it is assumed that the system has a dominant second order transient which
has complex poles at a radius ωn and a damping factor ζ. The rest of the poles are
selected to have real parts corresponding to sufficiently damped modes of the system so
that the system approximates to a second order system (see Figure 1.10). This method
is commonly used when the performance specifications are given in terms of transient
response parameters e.g. rise time, overshoot, settling time.

1.5.3 Prototype design

This method is used when the system cannot be usefully approximated as a second
order system.
There are different prototypes and they have different characteristics. One category
of prototypes is based on performance indices. A performance index is a number that
indicates the “goodness” of a system performance.
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 29

j ω axis
s−plane

ωn sin θ = ζ
θ
real axis

Figure 1.10: Dominant second order poles

1.5.4 Integral of Time Multiplied by Absolute Error (ITAE) proto-


type

One commonly used index is the Integral of Time Multiplied by Absolute Error (ITAE)
criterion. This method minimizes the performance index
Z ∞
J= t|e(t)|dt. (1.87)
0

The ITAE criterion has been applied to a closed-loop transfer function of the form
Y (s) a0
= n (1.88)
U (s) s + an−1 sn−1 + · · · + a1 s + a0
and the coefficients which minimize the performance index J in equation (1.87) have
been determined. The ITAE characteristic equations for systems of order 1 to order 6
are given below (in each case, a0 = (ω0 )n ):

1st order: s + ω0 (1.89a)


2nd order: s2 + 1.4ω0 s + ω02 (1.89b)
3rd order: s3 + 1.75ω0 s2 + 2.15ω02 s + ω03 (1.89c)
4th order: s4 + 2.1ω0 s3 + 3.4ω02 s2 + 2.7ω03 s + ω04 (1.89d)
5 4
5th order: s + 2.8ω0 s + 5.0ω02 s3 + 5.5ω03 s2 + 3.4ω04 s+ω05 (1.89e)
6th order: s6 +3.25ω0 s5 +6.6ω02 s4 +8.6ω03 s3 +7.45ω04 s2 +3.95ω05 s+ω06 (1.89f)

Design of a system to meet the ITAE criterion results in a system with a small rise time and
zero steady-state error to a step input.
The designer tries different values of ω0 and selects the value which gives the best
response. In general, increasing the value of ω0 makes the system faster (smaller rise
time and smaller settling time).
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 30

1.5.5 The Bessel prototype

The Bessel prototype is based on Bessel functions. For a closed-loop transfer function
of the form of equation (1.88), the Bessel characteristic equations for systems of order
1 to order 6 are given below (in each case, a0 = (ω0 )n ):

1st order: s + ω0 (1.90a)


2nd order: s2 + 1.73ω0 s + ω02 (1.90b)
3rd order: s3 + 2.43ω0 s2 + 2.47ω02 s + ω03 (1.90c)
4th order: s4 + 3.12ω0 s3 + 4.39ω02 s2 + 3.20ω03 s + ω04 (1.90d)
5th order: s5 + 3.81ω0 s4 + 6.78ω02 s3 + 6.89ω03 s2 + 3.94ω04 s+ω05 (1.90e)
6th order: s6 +4.5ω0 s5 +9.62ω02 s4 +12.36ω03 s3 +9.92ω04s2 +4.67ω05 s+ω06 (1.90f)

Again, the designer tries different values of ω0 and selects the value which gives the best
response. In general, increasing the value of ω0 makes the system faster (smaller rise
time and smaller settling time).
Systems designed using the Bessel prototypes has no overshoots and zero steady-state error
for step inputs. In addition, for a given value of ω0 , the step response of a Bessel prototype
is slower than the step response of an ITAE prototype of the same order.
Figure 1.11 shows the step responses of the ITAE prototypes given by equations (1.89a)
to (1.89f) while Figure 1.12 shows the step responses of the Bessel prototypes given by
equations (1.90a) to (1.90f). In all cases, a value of ω0 = 1 is used. As an example, to
simulate the 6th order ITAE step response, the (closed-loop) transfer function used is

Y (s) ω06
= 6
U (s) s + 3.25ω0 s5 + 6.6ω02 s4 + 8.6ω03 s3 + 7.45ω04 s2 + 3.95ω05 s + ω06

As can be seen from Figure 1.12, the Bessel prototypes result in a closed-loop system
with no overshoots.

1.5.6 Example:

Given a second order system:


Y (s) 10
=
U (s) s(s + 2)
determine the feedback matrix K so that the control law u(t) = −Kx(t) + k1 r(t), where
r(t) is a unit step input, results in closed-loop poles at the second order ITAE locations
for ω0 = 2. (The second order ITAE polynomial is given as s2 + 1.4ω0 s + ω02 ).
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 31

1.2 1.2

1st order system 1st order system


1 1

0.8 0.8

0.6 0.6

6th order system 6th order system

0.4 0.4

0.2 0.2

0 0
0 2 4 6 8 10 12 0 2 4 6 8 10 12
time −−> (sec) time −−> (sec)

Figure 1.11: Step response – ITAE prototypes Figure 1.12: Step response – Bessel prototypes

Solution:

Y (s) 10
= 2 (1.91)
U (s) s + 2s
Hence
s2 Y (s) + 2sY (s) = 10U (s) (1.92)

Converting equation (1.92)into time domain gives

ÿ(t) + 2ẏ(t) = 10u(t) (1.93)

Since this is a system with an external input, the states are selected as x1 (t) = y(t) and
x2 (t) = ẏ(t).
The state equations can then be written as:

ẋ1 (t) = ẏ(t) = x2 (t) (1.94a)


ẋ2 (t) = ÿ(t) = −2ẏ(t) + 10u(t) = −2x2 (t) + 10u(t) (1.94b)

Equations (1.94a) and (1.94b) can be written in vector-matrix notation as:


" # " #" # " #
ẋ1 (t) 0 1 x1 (t) 0
= + u(t) = Ax(t) + Bu(t) (1.95a)
ẋ2 (t) 0 −2 x2 (t) 10
" #
x1 (t)
y(t) = [1 0] = Cx(t). (1.95b)
x2 (t)

The controllability matrix


" #
0 10
M = [B AB] = . (1.96)
10 −20
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 32

|M| = −100 6= 0, (1.97)


hence the system is completely state controllable and design by pole placement is pos-
sible.
For ω0 = 2, a second order ITAE polynomial is given by

s2 + 2.8s + 4 (1.98)

which is the desired characteristic equation.


Hence
" #" # " # " # " #
0 1 0 1 0 1 1 0 4 0.8
φ(A) = A2 +2.8A+4I = +2.8 +4 = .
0 −2 0 −2 0 −2 0 1 0 2.4
(1.99)
Using Ackermann’s formula,
" #−1 " #
0 10 4 0.8
K = [0 1] [B AB]−1 φ(A) = [0 1] = [0.4 0.08]. (1.100)
10 −20 0 2.4

The state equation of the original (open-loop) system is

ẋ = Ax(t) + Bu(t), (1.101)

while the control law is given as

u(t) = −Kx(t) + k1 r(t). (1.102)

Substituting (1.102) into (1.101) gives

ẋ = (A − BK)x(t) + Bk1 r(t), (1.103)

which is the state equation of the closed-loop system. The output equation of the
closed-loop system is the same as the output equation of the open-loop system (equation
(1.95b).
Note that the transfer function of the closed-loop system is given by

Y (s) 4
= C (sI − A + BK)−1 Bk1 = 2 (1.104)
R(s) s + 2.8s + 4

which is the second order ITAE prototype for ω0 = 2.


The step response of the closed-loop system can be simulated using MATLAB or Simulink.
The MATLAB code is given below:
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 33

A = [0 1;0 -2]; %state matrix of the open-loop system


B = [0;10]; %input matrix of the open-loop system
ITAEpoly2 = [1 2.8 4]; %second order ITAE polynomial for w0 = 2
p = roots(ITAEpoly2); %roots of the polynomial
K = acker(A,B,p) %(the command K = place(A,B,p) can also be used)
A1 = A-B*K; %state matrix of the closed-loop system
B1 = B*K(1); %input matrix of the closed-loop system
C = [1 0]; %output matrix
D = 0; %direct transmission matrix
sys1 = ss(A1,B1,C,D); %define a state space object named sys1
tmax = 20; %upper time limit for the step response
figure,step(sys1,tmax) %step response

A Simulink model to simulate resulting system is shown in Figure 1.13 (the transfer
4
function block s2 +2.8s+4 is included for comparison of the system response with the
expected response).

Figure 1.13: Simulink model for given example

Exercise
Repeat the above example, this time designing for a closed-loop system with poles
at the second order Bessel locations for ω0 = 2. Show that the closed-loop transfer
function of the resulting system is given by

Y (s) 4
= 2 .
R(s) s + 3.46s + 4
Chapter 2

LYAPUNOV STABILITY

2.1 Introduction

Stability plays a major role in control system design. There are many methods for
determining the stability of a given system. In this course, it is assumed that the reader
is familiar with the following methods of determining the stability of a system:

Routh stability criterion: based on the coefficients of the characteristic equation

Nyquist criterion: based on the frequency response of the system.

The above methods are known as classical methods, and they are only applicable if the
system being considered is Linear-Time-Invariant (LTI).
However, if the system is nonlinear, or linear but time varying, the methods above
cannot be used to determine the stability of such systems.
The Lyapunov stability methods are more general methods of determining the stabil-
ity of systems, and they are applicable to both linear and non-linear, as well as time-
invariant and time-varying systems. These methods are named after a Russian scientist,
A. M. Lyapunov, who studied dynamic systems developed the methods in the 19th cen-
tury.
Lyapunov developed 2 methods of determining the stability of systems: the first method
(also known as the indirect method) and the second method (a.k.a. the direct method).

First method (indirect method): This method involves solving of the differential equa-
tions (linear or nonlinear) describing the system to determine the stability of the
system. This method is not commonly used because it is usually difficult and
time-consuming to solve nonlinear differential equations.

34
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 35

Second method (the direct method): Does not require the solution of the differential
equations describing the system hence by using this method, the stability of the
system can be determined without having to solve the state equations.
The second method of Lyapunov is the most general method for determination of
the stability of nonlinear and/or time-varying systems as well as hybrid/switched
systems. This method also applies to linear-time-invariant systems.

2.2 Equilibrium points

Consider an autonomous system (i.e. a system without external inputs)

ẋ = f (x, t) (2.1)

where x is the state vector with components x1 , x2 , . . . , xn , and f (x, t) is an n-state


vector whose elements are functions of x1 , x2 , . . ., xn and t.
We assume that the system has a unique solution starting at a given initial state x(0).
Suppose we have an autonomous system in which all states have settled down to con-
stant values (not necessarily zero). Such a system is said to be in equilibrium.
The state xe is an equilibrium state of the system (2.1) if and only if

ẋe = f (xe , t) = 0 for all t, (2.2)

i.e. since the states have settled down to constant values, ẋe = 0.
Example
Find the equilibrium state of the autonomous LTI system described by the state equation
" # " #" #
ẋ1 (t) −1 2 x1 (t)
= (2.3)
ẋ2 (t) 1 −1 x2 (t)

Solution
At the equilibrium state xe , equation (2.2) applies hence

0 = −x1 + 2x2 (2.4a)


0 = x 1 − x2 (2.4b)

which gives a solution x1 = x2 = 0, hence the equilibrium state for this system is at the
origin, i.e. (x1e , x2e ) = (0, 0).
In general, a linear time-invariant autonomous system

ẋ = Ax (2.5)
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 36

has a single equilibrium point at the origin if the matrix A is non-singular. If A is


singular, system (2.5) has infinity equilibrium points.
Example
Find the equilibrium state of the nonlinear system described by the state equations

ẋ1 = x2 − x1 x21 + x22



(2.6a)
ẋ2 = −x1 − x2 x21 + x2 2

(2.6b)

Solution
Let the equilibrium point be (x1 , x2 ) = (x1e , x2e ). At the equilibrium point, ẋ1e = ẋ2e = 0
hence we can write

0 = x2e − x1e x21e + x22e



(2.7a)
0 = −x1e − x2e x21e + x22e

(2.7b)

Multiplying equation (2.7a) by x1e and (2.7b) by x2e gives

0 = x1e x2e − x21e x21e + x22e



(2.8a)
x22e x21e x22e

0 = −x1e x2e − + (2.8b)

Adding equation (2.8a) to equation (2.8b) gives


2
0 = x21e + x22e , (2.9)

whose only solution is x1e = x2e = 0.


Example
Find the equilibrium state(s) of the nonlinear system described by the state equations

ẋ1 = x2 (2.10a)
ẋ2 = −x1 − x21 − x2 (2.10b)

Solution
At equilibrium,

0 = x2e (2.11a)
0 = −x1e − x21e − x2e (2.11b)

hence x2e = 0 and x1e =0 or -1, hence the equilibrium states for this system are (x1e , x2e ) =
(0, 0) and (x1e , x2e ) = (−1, 0).
In general, a non-linear system may have more than one equilibrium points.
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 37

Fixed
surface

R
θ

Figure 2.1: Pendulum

Example
Consider the pendulum shown in Figure 2.1. It is assumed that the weight at the bottom
of the pendulum is joined to the fixed surface using a rigid body (not a string!) of
negligible mass.
The dynamics of the pendulum are given by the nonlinear differential equation

M R2 θ̈ + bθ̇ + M gR sin θ = 0 (2.12)

where R is the pendulum’s length, M its mass, b the frictional coefficient at the hinge,
and g the gravitational acceleration. The angle θ is in radians.
Defining two states x1 and x2 such that x1 = θ and x2 = θ̇, the state equations of the
system are given by

ẋ1 = x2 (2.13a)
g b
ẋ2 = − sin x1 − x2 (2.13b)
R M R2
At equilibrium, x2e = 0, and sin x1e = 0. This gives two equilibrium states (x1e , x2e ) =
(0, 0) and (x1e , x2e ) = (π, 0). Physically, these points correspond to the pendulum resting
in the bottom position and the upright position respectively.

2.3 Stability in the sense of Lyapunov

Consider the two-dimensional state space shown in Figure 2.2.


Suppose the system described by the state equation (2.1) has an initial state x0 . The
path traced by the state x(t) of the system for t ≥ 0, is known as the system trajectory
(see Figure 2.2).
For a two-dimensional state space such as the one shown in Figure 2.2, the Euclidean
norm is defined as q
kxk = x21 + x22 . (2.14)
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 38

x2 2−dimensional
state space
x0

0 x1
PSfrag replacements

system
trajectory

Figure 2.2: 2-dimensional state space

For an n-dimensional state space, the Euclidean norm is defined as


q
kxk = x21 + x22 + . . . + x2n . (2.15)

Suppose that the state vector x lies within a hyperspherical region of radius R, i.e.
kxk < R. (In a 2-dimensional system, the hyperspherical region would be a circle of
radius R, for a 3-dimensional system, the hyperspherical region would be a sphere of
radius R, etc.)
Define two regions with radius δ and ε within R such that δ < ε < R (see Figure 2.3 for
a 2-dimensional state-space).
x2
2−dimensional
state space

R
ε

0 δ x1

PSfrag replacements

Figure 2.3: 2-dimensional hyperspherical region

Definition: Lyapunov stability


A system is said to be stable in the sense of Lyapunov (or Lyapunov stable) if for every
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 39

radius ε, there exists a radius δ < ε so that every trajectory starting within the region of
radius δ will always remain within radius ε with continually increasing time.
,→ If an equilibrium state is Lyapunov stable, the state trajectory will remain within a
given neighbourhood if the initial state is chosen close to the equilibrium point.
Definition: Asymptotic stability
If in addition to satisfying the Lyapunov stability condition the trajectory starting within a
radius δ converges to the origin with continuing time, the system is said to be asymptotically
stable.
Definition: Asymptotic stability in the large
If a system is Lyapunov stable and all trajectories starting anywhere in the state space
converge to the origin with continuing time, the system is said to be asymptotically stable
in the large. Sometimes this referred to as global stability.
Note that for linear-time-invariant systems, asymptotic stability and asymptotic stability
in the large are the same.
Figure 2.4 shows three trajectories: trajectory 1 is a Lyapunov-stable trajectory, trajec-
tory 2 is asymptotically stable, while trajectory 3 is unstable.

x2
3 2−dimensional
state space

1 δ x1

R
PSfrag replacements

Figure 2.4: stable and unstable trajectories

2.4 Linearization and local stability

Lyapunov’s linearization method is concerned with the local stability of a system. To use
this method, it is assumed that the equilibrium state of the system is at the origin. The
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 40

nonlinear system is then approximated by a linear model (i.e. the system is linearized)
and the behaviour of the linear system analyzed.
Consider a scalar function f (h). For small x, Taylor’s approximation states that
∂f (h) 1 ∂ 2 f (h) 2
f (h + x) ≈ f (h) + x+ x +... (2.16a)
∂h 2! ∂h2
∂f (h)
≈ f (h) + x + higher order terms. (2.16b)
∂h
If h = 0, equation (2.16b) can be written as
∂f (h)
f (x) ≈ f (0) + x + higher order terms. (2.17)
∂h h=0

Now consider the system


ẋ = f (x) . (2.18)
Using Taylor’s theorem for small changes of state x around the origin, we can approxi-
mate this system as
∂f (x)
ẋ = f (x) ≈ f (0) + x + higher order terms. (2.19)
∂x x=0

In practice, linearization (or linear approximation) of the system is done by neglecting


the higher order terms. Furthermore, since it is assumed that the equilibrium state of
the system is at the origin, then f (0) = 0, hence equation (2.19) becomes
∂f (x)
ẋ ≈ x. (2.20)
∂x x=0

Equation (2.20) is used to approximate the nonlinear system with a linear one. Note
that since Taylor’s theorem is only valid for small variations around a point, the model
obtained using equation (2.20) is only useful in analyzing the stability in the neighbour-
hood of the equilibrium point at the origin.
Note that equation (2.20) is a vector matrix equation, and
∂f1 ∂f1 ∂f1
 
 ∂x1 ∂x2 · · · ∂xn 
 ∂f2 ∂f2 ∂f2 
 
∂f (x)  ···
∂x1 ∂x2 ∂xn  ,

= (2.21)
∂x  ... .. .. .. 

 . . . 

 ∂fn ∂fn ∂fn 
···
∂x1 ∂x2 ∂xn
which is known as a Jacobian.
For a second order system, the approximation can therefore be written as
∂f1 ∂f1
 
" # " #
ẋ1  ∂x1 ∂x2  x1
≈  ∂f ∂f2  (2.22)
ẋ2 2 x2
∂x1 ∂x2 x1 =x2 =0
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 41

The stability of the linearized system is analyzed as follows:

 If the eigenvalues of the linearized system lie strictly on the left half plane, the
equilibrium point of the nonlinear system is asymptotically stable.

 If at least one of the eigenvalues of the linearized system lies on the right half
plane, the equilibrium point of the nonlinear system is unstable.

 If the eigenvalues of the linearized system lie on the left half plane, but at least
one of them lies on the imaginary axis, no conclusion can be drawn about the
stability of the equilibrium point of the nonlinear system.

Example
Consider the nonlinear system

ẋ1 = x22 + x1 cos x2 (2.23a)


ẋ2 = x2 + x21 + x1 + x1 sin x2 . (2.23b)

The linearized approximation is given by


" # " #" #
ẋ1 1 0 x1
= . (2.24)
ẋ2 1 1 x2

The linearized system has eigenvalues on the right half plane so the equilibrium point
of the nonlinear system is unstable.

2.5 Lyapunov functions

2.5.1 Positive define and negative definite functions

Let us define a scalar function V (x) (this is a scalar function of a vector x).
Definition: Positive definite function
A scalar function V (x) is positive definite in a region if

V (x) > 0 if x 6= 0 (2.25a)


V (x) = 0 only when x = 0. (2.25b)

Definition: Negative definite function


A scalar function V (x) is negative definite if −V (x) is positive definite
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 42

Definition: Positive semi-definite function


A scalar function V (x) is positive definite in a region if

V (x) ≥ 0 if x 6= 0 (2.26a)
V (x) = 0 only when x = 0. (2.26b)

Definition: Negative semi-definite function


A scalar function V (x) is negative semi-definite if −V (x) is positive semi-definite
Examples

V (x) = x21 + x22 : positive definite (2.27a)


V (x) = (x1 + x2 )2 : positive semi-definite (2.27b)
V (x) = −x21 − (3x1 + 2x2 )2 : negative definite (2.27c)
V (x) = x1 x2 + x22 : indefinite (2.27d)
2 2x22
V (x) = x1 + : positive definite (2.27e)
1 + x22

The examples are fairly simple, and it is fairly straight-forward to tell which functions
are positive/negative definite/semi-definite. However, in general, there is no general
way to tell if a function is positive or negative definite. However, if V (x) is in quadratic
form, we can use Sylvester’s criterion to determine if V (x) is positive definite.

2.5.2 The quadratic form

The quadratic form is given by


  
p11 p12 p13 · · · p1n x1
p12 p22 p23 · · · p2n x2
  
  
T
  
V (x) = x Px = [x1 x2 x3 . . . xn ]  p13 p23 p33 · · · p3n  x3 , (2.28)
.. .. .. . . .. ..
  
. . . . . .
  
  
p1n p2n p3n · · · pnn xn

where P is a real, symmetric matrix.


Sylvester’s criterion states that the necessary and sufficient condition that the quadratic
form V (x) = xT Px is positive definite is that all successive principal minors of P be
positive, i.e.

p11 p12 p13


p11 p12
p11 > 0, > 0, p12 p22 p23 > 0, ··· P > 0. (2.29)
p12 p22
p13 p23 p33
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 43

Example
Use Sylvester’s criterion to show that the function below is positive definite:

V (x) = 10x21 + 4x22 + x33 + 2x1 x2 − 2x2 x3 − 4x1 x3 . (2.30)

Solution
For a 3rd order system, the quadratic form is given by
  
p11 p12 p13 x1
V (x) = xT Px = [x1 x2 x3 ]  p12 p22 p23   x2 
  

p13 p23 p33 x3


= p11 x21 + p22 x22 + p33 x23 + 2p12 x1 x2 + 2p13 x1 x3 + 2p23 x2 x3 .

It follows that  
10 1 −2
P= 1 4 −1  ,
 

−2 −1 1

and this matrix satisfies Sylvester’s criterion for positive-definiteness so the function
given in equation (2.30) is positive definite.

2.5.3 Lyapunov functions

Suppose the function V (x) has the following properties:

1. V (0) = 0.

2. V (x) > 0 for x 6= 0.

3. V (x) is continuous and has continuous derivatives with respect to all components
of x.

4. V̇ (x) ≤ 0 (i.e. the time derivative of V (x) is less or equal to zero) along the
trajectories of the system.

Such a function is known as a Lyapunov function.


If there exists such a function V (x) and it has continuous partial derivatives and V̇ (x)
along the trajectories of the system, then the equilibrium state at the origin is asymp-
totically stable.
Conversely, if the origin for a particular system is asymptotically stable, then Lyapunov
functions with the required properties will always exist.
Please note that
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 44

• Finding a Lyapunov function may be quite difficult.

• Failure to find a Lyapunov function for a particular system does not mean that a
Lyapunov function for that system does not exist, and in such cases, no conclusions
about the stability of the system can be made.

• The Lyapunov function for a system is not unique.

2.6 Lyapunov stability analysis for LTI systems

Consider the linear-time-invariant (LTI) system

ẋ = Ax. (2.31)

Let V (x) be in quadratic form, i.e. V (x) = xT Px. By use of the chain rule, we can write
that
V̇ (x) = ẋT Px + xT Pẋ. (2.32)
Substituting (2.31) in (2.32) gives

V̇ (x) = (Ax)T Px + xT P (Ax) = xT AT P + PA x = −xT Qx.



(2.33)

where
Q = − AT P + PA .

(2.34)

To satisfy the condition V̇ (x) < 0, then xT Qx must be positive definite.


Instead of specifying a positive-definite matrix P and examining whether or not Q is
positive definite, it is more convenient to specify a positive definite matrix Q first and

then work back through the equation Q = − AT P + PA to determine P, and then
check if the matrix P is positive definite. If P is positive definite, then the system is
asymptotically stable. Usually, Q is chosen as the identity matrix.
Note that any P and Q matrices where P is real-symmetric and both P and Q are

positive definite satisfying Q = − AT P + PA give sufficient proof of stability.

To determine matrix P, we equate AT P + PA to −Q element by element.
If the system is stable, then xT Px is the Lyapunov function for the system.
Example
Use Lyapunov’s direct method to determine if the system below is stable:
" # " #" #
ẋ1 −1 2 x1
= . (2.35)
ẋ2 1 −1 x2
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 45

Solution
Let
" # " #
p11 p12 1 0
P= (a real-symmetric matrix), and Q = . (2.36)
p12 p22 0 1

Then
" #" # " #" # " #
T
 −1 1 p11 p12 p11 p12 −1 2 −1 0
A P + PA = + = −Q =
2 −1 p12 p22 p12 p22 1 −1 0 −1
(2.37)
or " # " #
−2p11 + 2p12 2p11 − 2p12 + p22 −1 0
= (2.38)
2p11 − 2p12 + p22 4p12 − 2p22 0 −1
which gives the simultaneous equations

−2p11 + 2p12 = −1 (2.39a)


2p11 − 2p12 + p22 = 0 (2.39b)
4p12 − 2p22 = −1 (2.39c)

Solving these equations gives p11 = −0.25, p12 = −0.75 and p22 = −1, hence
" #
−0.25 −0.75
P= . (2.40)
−0.75 −1

The matrix P above does not meet Sylvester’s criteria for positive-definiteness hence we
conclude that the given system is not stable, and a Lyapunov function for this system
does not exist.
Example
Repeating the above example for the system
" # " #" #
ẋ1 −5 −4 x1
= (2.41)
ẋ2 2 1 x2

gives " #
1 7
3 12
P= 7 11
. (2.42)
12 6

In this case, p11 > 0 and |P| > 0 (satisfying Sylvester’s criterion for positive-definiteness),
hence xT Px is a positive definite function and the system is asymptotically stable. The
Lyapunov function is given by
" #" #
1 7
x 1 1 2 7 11 2
V (x) = xT Px = [x1 x2 ] 37 12 11
= x 1 + x 1 x 2 + x2 . (2.43)
12 6
x 2 3 6 6
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 46

Also note that


V̇ (x) = −xT Qx = − x21 + x22 ,

(2.44)
which is a negative-definite function.
Note that this method for linear systems does not work if the system has poles along the
imaginary axis – in such a case, equation (2.34) gives inconsistent simultaneous equations
and no conclusions can be drawn about the stability of such a system.
Chapter 3

DIGITAL CONTROL

3.1 Introduction to digital control systems

Digital control refers to the use of digital devices or digital computers to control dynamic
systems. The theory of digital control is very wide, and cannot be covered in a single
course. This chapter covers the basic materials of digital control systems. For more
details, the reader is referred to the references at the end of the handout.
PSfrag replacements

3.1.1 System configuration and terminology

Consider the closed-loop system shown in Figure 3.1.


controller plant
R(s) E(s) U (s) Y (s)
C(s) G(s)
r(t) + e(t) u(t) y(t)

C(s)
G(s)

Figure 3.1: continuous-time closed-loop system

The controller C(s) is a continuous-time controller and it can be implemented using an


analog electronic circuit.
The last two decades have seen an exponential increase in the computing power and
speed of digital computers as well as the dramatic fall in their prices, factors which have
led to their widespread use in process control. As a result, many controllers are nowa-
days implemented in digital form rather than analog form. Figure 3.2 shows a closed-
loop system similar to the one shown in Figure 3.1, but with the continuous controller

47
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 48

replaced with a digital controller. Since the digital controller works in discrete-time,
replacementssuch a system is also referred to as a discrete-time control system.
plant
r(t) e(t) sample e(kT ) u(kT )
D/A
u(t) y(t)
and A/D digital smoothing G(s)
+ hold converter controller converter filter

G(s)

Figure 3.2: closed-loop system with a digital controller

Sample-and-Hold Circuit: used to hold the analog input constant when the Analog-
to-Digital converter is performing the conversion

Analog-to-Digital Converter: Converts the analog input into a binary (digital) code.
Note that many Analog-to-Digital converters have built-in sample-and-hold cir-
cuits.

Digital Controller: Computes the control signal u(kT ).

Digital-to-Analog converter: Converts the control signal given in the form of a binary
code to an analog signal.

Smoothing filter: Used to smooth the output of the Digital-to-Analog converter.

3.1.2 Advantages and disadvantages of digital control

The advantages of digital controllers over analog controllers are:

• The control law can easily be modified by changing the program.

• Digital controllers have high repeatability since digital devices do not suffer the
drift associated with analog devices.

• It is easier and more economical to handle large control laws (i.e. high order
controllers).

Disadvantages:

• In general, converting a continuous-time control system to a digital control system


introduces time delays in the loop which degrades the system stability margin.
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 49

• Mathematical analysis and design of digital control systems is generally more com-
plicated than the analysis of continuous-time control systems.

3.2 The z-transform

3.2.1 Definition of the z-transform

The z-transform serves the same role for discrete-time control systems that the s-transform
(Laplace transform) serves for continuous-time control systems.
Consider the continuous-time signal f (t) shown in Figure 3.3

f (t)

PSfrag replacements

0 t

Figure 3.3: A continuous-time signal

Suppose the signal f (t) is sampled every T seconds using an ideal sampler. The ideal
sampler is represented as a switch as shown in Figure 3.4.

PSfrag replacements f (t) f ∗ (t)

ideal sampler
interval T

Figure 3.4: An ideal sampler

The output of the ideal sampler can be thought of as a train of impulses whose respective
impulse areas are equal to the magnitudes of f (t) at the corresponding sampling instants
t = kT , k = 0, 1, 2, . . . (see Figure 3.5).
The output of the ideal sampler at various time instants is tabulated below:
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 50
PSfrag replacements

f ∗ (t)

0 T 2T 3T 4T 5T 6T t

Figure 3.5: The sampled signal

Time Sampler output


t=0 f (0)δ(t)
t=T f (T )δ(t − T )
t = 2T f (2T )δ(t − 2T )
t = 3T f (3T )δ(t − 3T )
.. ..
. .
t = kT f (kT )δ(t − Kt)
.. ..
. .
The output of the sampler can thus be written as

X
f (t) =

f (kT )δ(t − kT ). (3.1)
k=0

The Laplace transform of the sampled signal f ∗ (t) is given by


Z ∞X∞
L (f (t)) = F (s) =
∗ ∗
f (kT )δ(t − kT )e−st dt. (3.2)
0 k=0

The integration and the summation can be interchanged, so we can write


∞ Z ∞
X
F (s) =

f (kT )δ(t − kT )e−st dt (3.3a)
k=0 0

X∞ Z ∞
= f (kT ) δ(t − kT )e−st dt. (3.3b)
k=0 0

From signals theory,


Z ∞
δ(t − kT )g(t)dt = g(kT ) = g(t)|t=kT , (3.4)
0

so equation (3.3) reduces to



X
F ∗ (s) = f (kT )e−skT . (3.5)
k=0
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 51

The exponential term e−skT in equation (3.5 introduces nonlinear functions in the
Laplace transform of discrete-time systems. This problem is resolved by the change
of variable

1
z = esT =⇒ s= ln z. (3.6)
T
With this change of variable, equation (3.5) can then be written as

1 X
F ∗ (s = ln z) = F (z) = f (kT )z −k (3.7)
T k=0

which is the definition of the z transform.


The z-transform is usually denoted using the symbol z, i.e.

X
z (f (t)) = z (f (kT )) = z (f (k)) = f (kT )z −k (3.8a)
k=0
= f (0) + f (T )z −1 + f (2T )z −2
+ · · · + f (kT )z −k + · · · . (3.8b)

The z-transform is therefore an infinite power series in z −1 .


It is assumed that f (t) = 0 for t < 0. This is called a one-sided z-transform (or unilateral
z-transform) because it assumes f (t) = 0 for t < 0.
If f (t) 6= 0 for t < 0, then

X
F (z) = f (kT )z −k (3.9)
k=−∞

which involves positive and negative powers of z. Equation (3.9) is known as the two-
sided or bilateral z-transform.
In general, any continuous functions which have a Laplace transform also have a z-
transform. It follows that when working with discrete-time systems, we can solve the
system equations using Laplace transform techniques but with a change of variable from
s to z so that the computations remain algebraic.
For a function f (t), the functions f (kT ) and F (z) are said to be a z-transform pair, i.e.
z
f (kT ) ←→ F (z). (3.10)
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 52

3.2.2 z-transforms of common functions

An important formula

The following formula from the theory of geometric series is useful in obtaining the z
transform of many functions:
N −1
2 N −1
X A − AxN
A + Ax + Ax + · · · + Ax = Axn = . (3.11)
n=0
1−x

If |x| < 1, then xN −→ 0 as N −→ ∞, so



X A
Axn = for |x| < 1. (3.12)
n=0
1−x

Step function
(
1 t≥0
f (t) = (3.13)
0 t<0
hence f (kT ) = 1. It follows that

F (z) = 1 + z −1 + z −2 + · · · + z −k + · · · (3.14a)
1 z
= = for |z −1 | < 1 or |z| > 1. (3.14b)
1−z −1 z−1
In other words,
z
z
1 ←→ for |z| > 1. (3.15)
z−1
The region of convergence is the outside of a circle of unit radius centered at the origin.

Geometric sequence ak

The z-transform of a geometric sequence ak where k is an arbitrary constant is given by

∞ ∞
X
k −k
X k
F (z) = a z = az −1 (3.16a)
k=0 k=0
2 −2
= 1 + az + a z + · · · + ak z −k + · · ·
−1
(3.16b)
1 z
= = for |az −1 | < 1 (or z > a) (3.16c)
1 − az −1 z−a
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 53

The exponential function e−at


(
e−at t ≥ 0
f (t) = . (3.17)
0 t<0

∞ ∞ ∞
X X X k
F (z) = f (kT )z −k
= e z
−akT −k
= e−aT z −1 (3.18a)
k=0 k=0 k=0
= 1 + e−aT z −1 + e z +···+e z +···
−2aT −2 −akT −k
(3.18b)
1 z
= = for |e−aT z −1 | < 1 (or |z| > e−aT ).(3.18c)
1−e z−aT −1 z−e −aT

The sinusoidal function


(
sin ωt t ≥ 0
f (t) = . (3.19)
0 t<0
But
1 jωt 
sin ωt = e − e−jωt . (3.20)
2j

z
z ejωt |z| > ejωT ,

= for (3.21a)
z − ejωT
 z
z e−jωt = for |z| > e−jωT , (3.21b)
z − e−jωT
hence
 
1 z z z sin ωT
z (sin ωt) = − = for |z| > 1. (3.22)
2j z − ejωT z − e−jωT z 2 − 2z cos ωT + 1

Similarly
z 2 − z cos ωT
z (cos ωt) = for |z| > 1. (3.23)
z 2 − 2z cos ωT + 1

The impulse function δ(t)



X
z (δ(t)) = δ(kT )z −k = δ(0) + δ(T )z −1 + δ(2T )z −2 + · · · = 1. (3.24)
k=0

3.2.3 Theorems of z-transform

Multiplication by a constant

If
z (x(t)) = X(z), (3.25)
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 54

then
z (ax(t)) = aX(z). (3.26)
Proof: ∞ ∞
X X
z (ax(t)) = ax(k)z −k
=a x(k)z −k = aX(z). (3.27)
k=0 k=0

Linearity

If
x(k) = αf (k) + βg(k), (3.28)
then
X(z) = αF (z) + βG(z). (3.29)

Multiplication by ak

If
z (x(k)) = X(z), (3.30)
then
z ak x(k) = X(a−1 z).

(3.31)
Proof:

X ∞
X z
k
 k
−k
z a x(k) = a x(k)z −k
= x(k) a z
−1
= X(a z) = X
−1
. (3.32)
k=0 k=0
a

Real translation theorem (shifting theorem)

If
z (x(t)) = X(z), (3.33)
then
z (x(t − nT )) = z −n X(z). (3.34)
Proof: ∞ ∞
X X
z (x(t − nT )) = x(kT − nT )z −k
= x((k − n)T )z −k . (3.35)
k=0 k=0

Let k − n = m. Then equation (3.35) can be rewritten as



X ∞
X
z (x(t − nT )) = x(mT )z −m−n
=z −n
x(mT )z −m . (3.36)
m=−n m=−n
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 55

A one sided z-transform of x(mT ) is zero for m < 0, hence the limit can be changed
from m = −n to m = 0, giving

X ∞
X
z (x(t − nT )) = z −n x(mT )z −m = z −n x(mT )z −m = z −n X(z). (3.37)
m=−n m=0

Example
Find the z-transform of the signal f (t) shown in Figure 3.6.

f(t)

0 2T t

Figure 3.6:

Solution
Figure 3.6 shows a step signal which has been shifted to the right by two sample inter-
vals.
1
A step signal starting at t = 0, denoted as 1(t), has a z-transform of 1−z −1
.
Using the real translation theorem, it follows that the z-transform of f (t) = 1(t − 2T )
shown in Figure 3.6 is given by
z −2
.
1 − z −1
Example
Find the z-transform of the signal f (t) shown in Figure 3.7.

f(t)

0 2T 4T t

Figure 3.7:

Solution
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 56

The signal shown in Figure 3.7 can be written as

f (t) = 1(t − 2T ) − 1(t − 4T ).

Using the real translation theorem,

z −2 z −4
F (z) = −
1 − z −1 1 − z −1

Complex translation theorem

If
z (x(t)) = X(z), (3.38)
then
z e−at x(t) = X zeaT .
 
(3.39)
Proof:

X ∞
X ∞
X −k
x(kT ) eaT z = X(zeaT ).
 
z e −at
x(t) = e −akT
x(kT )z −k
= x(kT ) e z
−akT −k
=
k=0 k=0 k=0
(3.40)

Initial value theorem

If
z (x(t)) = X(z), (3.41)
then
x(0) = lim X(z). (3.42)
z→∞

Proof:

X
X(z) = x(kT )z −k = x(0) + x(T )z −1 + x(2T )z −2 + · · · + x(kT )z −k + · · · . (3.43)
k=0

Letting z −→ ∞ in equation (3.43) yields equation (3.42).

Final value theorem

If
z (x(k)) = X(z) (3.44)
where x(k) = 0 for k < 0, then
  
lim x(k) = lim 1 − z −1 X(z) . (3.45)
k→∞ z→1
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 57

Proof:

X
z (x(k)) = X(z) = x(k)z −k , (3.46a)
k=0

X
z (x(k − 1)) = z −1
X(z) = x(k − 1)z −k , (3.46b)
k=0

hence ∞ ∞
X X
x(k)z −k − x(k − 1)z −k = X(z) − z −1 X(z). (3.47)
k=0 k=0

Taking the limit as z tends to 1,


"∞ ∞
#
X X   
lim x(k)z −k − x(k − 1)z −k = lim 1 − z −1 X(z) (3.48)
z→1 z→1
k=0 k=0

But " #

X ∞
X ∞
X ∞
X
lim x(k)z −k
− x(k − 1)z −k
= x(k) − x(k − 1). (3.49)
z→1
k=0 k=0 k=0 k=0

Recalling that x(k) = 0 for k < 0, we can expand the right hand side of equation (3.49)
to get

X ∞
X
x(k) − x(k − 1) = [x(0) − x(−1)] + [x(1) − x(0)] + [x(2) − x(1)] + · · · = x(∞),
k=0 k=0
(3.50)
hence
  
lim x(k) = lim 1 − z −1 X(z) . (3.51)
k→∞ z→1

3.2.4 Inverse z-transform

Introduction

There are 3 methods for working out the inverse z-transform:

(i) The direct division method.

(ii) Partial fraction expansion method

(iii) Inversion integral method.

In this course, we are going to look at the first two methods.


ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 58

Direct division method

This method involves the direct division of the numerator polynomial of X(z) by the
denominator polynomial of X(z). Before the division is done, both the numerator and
the denominator polynomials must be written in increasing powers of z −1 . By the divi-
sion, X(z) is expanded into an infinite power series in z −1 , and then the values of x(kT )
can be determined by inspection.
Example
Given
10z + 5
X(z) = , (3.52)
(z − 1)(z − 0.2)
find x(kT ) for k = 0, 1, 2, 3, 4.
Solution
First, X(z) is rewritten as a ratio of polynomials in z −1 :
10z −1 + 5z −2
X(z) = . (3.53)
1 − 2z −1 + 0.2z −2
Dividing the numerator by the denominator using long division yields

X(z) = 10z −1 + 17z −2 + 18.4z −3 + 18.68z −4 + · · · (3.54a)


= x(0) + x(T )z −1 + x(2T )z −2 + x(3T )z −3 + x(4T )z −4 + · · · . (3.54b)

It follows that x(0) = 0, x(T ) = 10, x(2T ) = 17, x(3T ) = 18.4 and x(4T ) = 18.68.
Note that in general, the direct division method does not yield a solution in closed form.

Partial fraction expansion method

The linearity property of the z-transform enables us to apply the partial fraction expan-
sion method to obtain the inverse z-transform of X(z).
Note that if X(z) has one or more zeros at the origin, it is easier to obtain the inverse
z-transform of X(z) by partial fraction expansion of X(z)
z
instead of the partial fraction
of X(z).
Example
Obtain the inverse z-transform of
10z
X(z) = . (3.55)
(z − 1)(z − 0.2)

X(z) has one zero at the origin. To show the difference between expanding X(z) and
X(z)
z
into partial fractions, both functions are expanded into partial fractions, starting
with X(z).
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 59

Partial fraction expansion of X(z) gives


12.5 2.5 12.5z −1 2.5z −1
X(z) = − = − . (3.56)
z − 1 z − 0.2 1 − z −1 1 − 0.2z −1

   
z −1 1
z −1
= z −1
z −1
(3.57a)
1 − z −1 1 − z −1
(
1 k = 1, 2, 3, . . .
= (3.57b)
0 k<1

Similarly, (
0.2k−1 k = 1, 2, 3, . . .
 
z −1
z−1 = (3.58)
1 − 0.2z −1 0 k<1
From (3.57) and (3.58) above,

x(k) = 12.5 − 2.5 (0.2)k−1 for k ≥ 1 (3.59a)


= 12.5 1 − 0.2 × (0.2)k−1

for k ≥ 1 (3.59b)
k

= 12.5 1 − (0.2) for k ≥ 1 (3.59c)
k

= 12.5 1 − (0.2) for k ≥ 0. (3.59d)

X(z)
Let us now expand z
into partial fractions:
 
X(z) 10 1 1
= = 12.5 − . (3.60)
z (z − 1)(z − 0.2) z − 1 z − 0.2
Hence
   
z z 1 1
X(z) = 12.5 − = 12.5 − . (3.61)
z − 1 z − 0.2 1−z −1 1 − 0.2z −1
x(k) = 12.5 1 − (0.2)k

for k ≥ 0. (3.62)

It can be seen from above that if a function has one or more finite zeros, it is easier to
obtain the inverse z-transform by expanding X(z)
z
into partial fractions than expanding
X(z) into partial fractions.

3.3 Difference equations and pulse transfer functions

3.3.1 Introduction

In time domain, dynamic systems in continuous-time are described using differential


equations. On the other hand, dynamic discrete-time systems are described using dif-
ference equations.
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 60

PSfrag replacements u(kT ) y(kT )


Discrete−time
input system output

Figure 3.8: discrete-time system

u(kT ) (sometimes written simply as u(k)) is a known input sequence, while y(kT ) (or
y(k)) is the output.
Consider a system described by the following linear difference equation:

y(kT ) + a1 y(kT − T ) + · · ·+ an y(kT − nT ) = b0 u(kT ) + b1 u(kT − T ) + · · ·+ bm u(kT − mT ).


(3.63)
This equation can be written more simply as

y(k) + a1 y(k − 1) + · · · + an y(k − n) = b0 u(k) + b1 u(k − 1) + · · · + bm u(k − m). (3.64)

Equation (3.63) (or (3.64)) is an nth order constant coefficient difference equation or an
nth order recursion equation.
y(kT ) (or y(k))is the value of the output at the current time, y(kT − T ) is the value of
the output at time (KT − T ) (the immediate past sampling instant), y(kT − 2T ) is the
value of the output at (kT −2T ) (the second-last sampling instant), and so on. Similarly,
u(kT ) is the value of the input at the current time, u(kT − T ) is value of the input at the
immediate past sampling instant, and so on.

3.3.2 Solving difference equations by recursion

Equation (3.63) (or (3.64)) can easily be solved using a digital computer by a procedure
known as recursion. To illustrate this procedure, suppose u(k) = 0 for k < 0 (i.e. the
input is zero for t < 0) as is usually the case.
Substituting k = 0 in (3.64), we can write

y(0) = −a1 y(−1) − a2 y(−2) − · · · − an y(−n) + b0 u(0). (3.65)

Recall that the input sequence is known, so, u(0) in the above equation is known. But
to determine y(0), the values y(−1), y(−2) · · · y(−n) must be known. These values are
known as the initial conditions of the difference equation (3.64).
Now substituting k = 1 into equation (3.64),

y(1) = −a1 y(0) − a2 y(−1) − · · · − an y(−n + 1) + b0 u(1) + b1 u(0). (3.66)

(Note that the value of y(0) determined from (3.65) is used in (3.66) above).
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 61

Similarly, substituting k = 2 into equation (3.64),

y(2) = −a1 y(1) − a2 y(0) − · · · − an y(−n + 2) + b0 u(2) + b1 u(1) + b2 u(0), (3.67)

and so on. This procedure of solving a difference equation is known as recursion. The
procedure can easily be programmed into a digital computer.
Note that in the special case where u(k) = 0 for all k (input is identically zero), the
difference equation (3.64) is known as a homogeneous difference equation, and the
solution to the equation is known as the zero input response.
Also note that if the initial conditions y(−1), y(−2) · · · y(−n) are not known, we can
only obtain a general solution to the difference equation (3.64), not a particular solu-
tion.

3.3.3 Pulse transfer functions and weighting sequence

Taking z-transform of (3.64), we can write

Y (z) + a1 z −1 Y (z) + · · · + an z −n Y (z) = b0 U (z) + b1 z −1 U (z) + · · · + bm z −m U (z), (3.68)

or
 
1 + a1 z −1 + · · · + an z −n Y (z) = b0 + b1 z −1 + · · · + bm z −m U (z). (3.69)
Hence
Y (z) b0 + b1 z −1 + · · · + bm z −m
= G(z) = . (3.70)
U (z) 1 + a1 z −1 + · · · + an z −n

G(z) is known as the pulse transfer function of the discrete-time system. The pulse
transfer function relates the input u(kT ) to the output y(kT ) at the sampling instants
only. It does not give any information about the nature of output between the sampling
instants.
The inverse z-transform of G(z), i.e. z−1 [G(z)] = g(k), is known as the weighting se-
quence of the system.
Example
Obtain the weighting sequence of the system described by the difference equation:

y(k) − ay(k − 1) = u(k). (3.71)

Solution:
Taking the z-transform of (3.71) above gives

Y (z) − az −1 Y (z) = U (z), (3.72)


ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 62

giving a pulse transfer function


Y (z) 1
= G(z) = . (3.73)
U (z) 1 − az −1
It follows that
g(k) = z−1 [G(z)] = ak . (3.74)
Exercise:
Obtain the weighting sequence of the system described by the difference equation:

y(k) + 3y(k − 1) + 2y(k − 2) = u(k). (3.75)

3.4 Mapping between the s-plane and the z-plane

3.4.1 Introduction

The stability of linear-time-invariant (LTI) continuous-time systems depends on the lo-


cations of the poles in the s-plane. For discrete-time systems, the stability of a system
depends on the location of its poles in the z-plane.

s = σ + jω
s−plane

0 σ
PSfrag replacements

Figure 3.9: s-plane

Consider the s-plane shown in Figure 3.9. Every point in the s-plane can be expressed
as a sum of a real number σ and an imaginary number jω, i.e.

s = σ + jω. (3.76)

In Section 3.2.1, it was shown that


z = esT , (3.77)
where T is the sampling time (the position of point in the z-plane depends on the its
position in the s-plane and the sampling time).
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 63

Substituting (3.76) into (3.77) gives

z = e(σ+jω)T = eσT ejωT . (3.78)

But ejωT = 1∠ωT , hence


z = eσT ∠ωT, (3.79)
i.e. z is a vector of magnitude eσT at an angle of ωT in the z-plane. Both the magnitude
and angle depend on the sampling time T .

3.4.2 Mapping of the jω axis to the z-plane

Along the jω axis (imaginary axis), σ = 0, hence

z = 1∠ωT, (3.80)

Equation (3.80) is the equation of a unit circle centered at the origin. The imaginary
axis on the s-plane is therefore mapped onto a unit circle centered at the origin in the
z-plane.
Note that the imaginary axis on the s-plane winds itself around the unit circle in the z-
plane infinitely many times - the origin in the s-plane is mapped to the point (1,0) in the
z-plane - then the positive part of the jω axis winds itself infinitely many times around
the unit circle in the anti-clockwise direction, while the negative part of the jω axis
winds itself around the unit circle in the clockwise direction. Points along the jω axis
in the s-plane whose frequencies differ by integral multiples of 2π T
are mapped to the
jπ j3π j5π jπ
same point in the z-plane, e.g. points T , T , T and − T are all mapped onto the point
1∠π = (−1, 0) in the z-plane.

3.4.3 Lines of constant attenuation

Figure3.10 shows two lines of constant, attenuation in the s-plane. The line to the right
of the jω axis is mapped onto a circle of radius eσ1 T > 1 centered at the origin in the
z-plane. The line to the left of the jω axis is mapped onto a circle of radius eσ2 T < 1
centered at the origin in the z-plane as shown in Figure 3.11.
From the above, we can deduce that the left half plane in the s-plane is mapped onto the
inside of the unit circle centered at the origin in the z-plane, while the right half-plane
in the s-plane is mapped onto the outside of the unit circle in the z-plane
,→ A linear-time-invariant discrete-time system is stable if all its poles lie inside of the unit
circle centered at the origin in the z-plane.
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 64

s−plane
PSfrag replacements
σ2 0 σ1 σ

Figure 3.10: lines of constant attenuation

z−plane unit circle

e σ1 T
0

PSfrag replacements
e σ2 T

Figure 3.11: mapping of lines of constant attenuation in the z-plane

3.4.4 Lines of constant frequency

Lines of constant frequency are shown in Figure 3.12.


The line at frequency jω1 is mapped to

z = e(σ+jω1 )T = eσT ejω1 T (3.81)

in the z-plane. This is a straight line starting at the origin in the z-plane at an angle ω1 T
radians to the horizontal axis as shown in Figure 3.13. Similarly, the line at a frequency
jω2 is mapped to z = eσT ejω2 T in the z-plane, while the line at the frequency −jω1 is
mapped to z = eσT e−jω1 T in the z-plane.
Example:
Consider the region in the s-plane that is shaded in Figure 3.14.
The shaded region is mapped into the z-plane as shown in Figure 3.15.
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 65


s−plane

PSfrag replacements jω2


jω1

0 σ
−jω1

Figure 3.12: lines of constant frequency

eσT

eσT
z−plane
PSfrag replacements
ω2 T
ω1 T
0 −ω1 T

eσT

Figure 3.13: mapping of lines of constant frequency in the z-plane


s−plane


  jω
PSfrag replacements 

 
  
2


 



 
  
 
 
 
 
 
 
 
0  σ
σ

 
 σ

 
2 1

jω1

Figure 3.14:
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 66

eσT  
                  
                  
                  
z−plane                   
                  
        ω   T         
       2
           
      0     ω   T       
         1
         
PSfrag replacements                   
eσT
unit circle

Figure 3.15:

3.4.5 Important note

Note that points in the s-plane with the same attenuation factors but whose frequencies
differ by integral multiples of 2π
T
are mapped onto the same locations in the z-plane.
This means that there are infinitely many values of s for each value of z, although a single
point in the s-plane corresponds to a single point in the z-plane.

3.5 Data holds

3.5.1 Introduction

?
x(t) x(kT )

PSfrag replacements t t

x(t) x(kT ) Data h(t)

Ideal Hold
sampler

Figure 3.16: Ideal sampler and a data hold


ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 67

 The ideal sampler converts a continuous time signal to a train of pulses occurring
at the sampling instants.

 The data hold is a device for generating a continuous-time signal from a discrete-
time sequence x(kT )

The signal h(t) during the time interval kT ≤ t ≤ (k + 1)T may be approximated by the
polynomial in τ as

h(kT + τ ) = an τ n + an−1 τ n−1 + . . . + a1 τ + a0 (3.82)

where 0 ≤ τ < T . Equation (3.82) is a polynomial extrapolator.


If τ = 0,
h(kT ) = a0 = x(kT ). (3.83)
Hence we can write that

h(kT + τ ) = an τ n + an−1 τ n−1 + . . . + a1 τ + x(kT ) (3.84)

If the data hold is an n-th order extrapolator, it is called an n-th order hold:

• n = 0: Zero order hold

• n = 1: First order hold

An n-th order hold uses the past (n + 1) discrete data samples x(kT − nT ), x(kT − nT +
T ), x(kT − nT + 2T ), . . . x(kT ) to generate the signal h(kT + τ ).
Because the higher data order hold uses past samples to extrapolate a continuous-time
signal between the present sampling instant and the next sampling instant, the accu-
racy of the approximation improves as the number of past samples used is increased.
However, this better accuracy is obtained at the cost of greater time delays since a large
number of past samples must be stored and the computation time increases with an
increase in the number of samples and this leads to increased time delays. In closed-
loop control systems, any added time delay in the loop will decrease the stability of the
system, and in some cases, it may lead to instability.

3.5.2 The Zero-Order-Hold (ZOH)

The simplest data hold is obtained when n = 0 and hence

h(kT + τ ) = x(kT ) for 0 ≤ τ < T, k = 0, 1, 2, . . . (3.85)


ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 68

x(t) x(kT ) h(t)

PSfrag replacements t t t

x(t) x(kT ) h(t)


ZOH
Ideal
sampler

Figure 3.17: Output of a Zero-Order Hold

This implies that the data hold circuit holds the amplitude of the sample constant from
one sampling instant to the next. Such a data hold is known as a ZERO ORDER HOLD
(ZOH).
The ZOH holds the discretized input at a constant value throughout the sample period
until a new sample is obtained, at which time the output changes. This results in a
discontinuous signal with steps at each sampling instant. Sometimes, the output of a
ZOH is passed through a low-pass filter to remove the discontinuities.
For the ZOH,

h(kT + τ ) = h(t) = x(kT ) for kT ≤ τ < (k + 1)T, (3.86)

hence we can write

h(t) = [1(t) − 1(t − T )] x(0) + [1(t − T ) − 1(t − 2T )] x(T ) + · · ·


+ [1(t − kT ) − 1(t − (k + 1)T )] x(kT ) + · · · (3.87)

which can be rewritten as



X
h(t) = [1(t − kT ) − 1(t − (k + 1)T )] x(kT ). (3.88)
k=0

Taking the Laplace transform of (3.88) gives


Z ∞
∞X
L(h(t)) = H(s) = [1(t − kT ) − 1(t − (k + 1)T )] x(kT )e−st dt. (3.89)
0 k=0

It is known that
1
L [1(t)] = (3.90)
s
and from the properties of Laplace transform, we can write that

e−kT s
L [1(t − kT )] = (3.91)
s
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 69

hence
∞ ∞
e−kT s − e−(k+1)T s 1 − e−T s −kT s
X   X  
L(h(t)) = H(s) = x(kT ) = x(kT ) e (3.92)
k=0
s k=0
s

which can be rewritten as


 ∞
1 − e−T s X

H(s) = x(kT )e−kT s (3.93)
s k=0

which is a product of 2 terms.


−kT s
P
Recall that the term ∞ k=0 x(kT )e is the Laplace transform of the sampled signal,
which was earlier defined as X (s) (see Section 3.2.1). It follows that the Laplace

transform of the ZOH is given by:


1 − e−T s
 
Gh (s) = .
s

1 − e−T s
 
In other words, a ZOH block can be replaced by the transfer function .
s

1 − e−T s
 
3.5.3 z-transform of functions involving the term
s
Suppose a transfer function G(s) follows a ZOH as shown in Figure 3.18.
PSfrag replacements ZOH System

1 − e−T s
 
G(s)
s

Figure 3.18: Transfer function G(s) following a ZOH

The overall transfer function is the product of the two transfer functions, and can be
written as:
1 − e−T s  G(s)
X(s) = G(s) = 1 − e−T s (3.94)
s s
The z-transform of the overall transfer function X(s), denoted by X(z), is given by
 
 G(s)
X(z) = z [X(s)] = 1 − z −1
z . (3.95)
s

Proof:
The overall transfer function is given by
1 − e−T s  G(s)
G(s) = 1 − e−T s = 1 − e−T s G1 (s)

X(s) = (3.96)
s s
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 70

G(s)
where G1 (s) = .
s
Let X1 (s) = e−T s G1 (s). Since X1 (s) is a product of two Laplace-transformed functions,
the inverse Laplace transform is the following convolution integral:
Z ∞
x1 (t) = g0 (t − τ )g1 (τ )dτ (3.97)
0

−T s

where g0 (t) = L−1 e = δ(t − T ) and g1 (t) = L−1 (G1 (s)), hence
Z ∞
x1 (t) = δ(t − T − τ )g1 (τ )dτ = g1 (t − T ). (3.98)
0

The z-transform of g1 (t − T ) is z −1 G1 (z), hence the z-transform of X(s) is


  
X(z) = z G1 (s) − e−sT G1 (s) = G1 (z) − z −1 G1 (z) = 1 − z −1 G1 (z) (3.99)

or  
 G(s)
X(z) = z [X(s)] = 1 − z −1
z (3.100)
s
End of Proof

PSfrag replacements
3.6 Block diagram analysis

3.6.1 Open-loop systems

Consider the system shown in Figure 3.19.

x(t) x∗ (t) y(t)


G(s)
δT
X(s) X ∗ (s) Y (s)

Figure 3.19:

δT is an ideal sampler with a sampling interval T .


From the block diagram, we can write that

Y (s) = G(s)X ∗ (s). (3.101)

Taking the starred Laplace transform on both sides of equation (3.101),

Y ∗ (s) = [G(s)X ∗ (s)]∗ = G∗ (s)H ∗ (s). (3.102)


ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 71

Since the z-transform is understood to be as the starred Laplace transform with eT s


replaced by z (see Section 3.2.1), we can then write

Y (z) = G(z)X(z). (3.103)

Proof:
Converting equation (3.101) into time domain,

y(t) = L−1 [G(s)X ∗ (s)] . (3.104)

Multiplication in the frequency domain is equivalent to convolution in the time domain,


hence
Z ∞ Z ∞ ∞
X
y(t) = g(t − τ )x (τ ) dτ =

g(t − τ ) x(nT )δ(τ − nT ) dτ (3.105a)
0 0 n=0
∞ Z ∞
X
= g(t − τ )x(nT )δ(τ − nT ) dτ (3.105b)
n=0 0

From the properties of the Dirac function, the above equation reduces to

X
y(t) = g(t − nT )x(nT ). (3.106)
n=0

At time t = kT ,

X
y(kT ) = g(kT − nT )x(nT ). (3.107)
n=0

The z-transform of y(kT ) is Y (z). Writing this mathematically,


∞ X
X ∞ ∞ X
X ∞
Y (z) = z [y(kT )] = g(kT − nT )x(nT )z −k = g((k − n)T )x(nT )z −k .
k=0 n=0 k=0 n=0
(3.108)
Letting k − n = m, the above equation becomes
∞ X
X ∞
Y (z) = g(mT )x(nT )z −(m+n) . (3.109)
m=−n n=0

Taking a unilateral z-transform,


∞ X
X ∞ ∞
X ∞
X
Y (z) = g(mT )x(nT )z −m z −n = g(mT )z −m x(nT )z −n = G(z)X(z).
m=0 n=0 m=0 n=0
(3.110)
End of Proof
Now consider the system shown in Figure 3.20.
For this system,
Y (s) = G(s)X(s). (3.111)
PSfrag replacements
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 72

x(t) y(t)
G(s)
X(s) Y (s)

Figure 3.20:

Taking the starred Laplace transform of (3.111) gives

Y ∗ (s) = [G(s)X(s)]∗ = [GX(s)]∗ , (3.112)

hence
Y (z) = z [GX(s)] . (3.113)
Note that the z-transform for Y (s) in Figure 3.19 (equation (3.103))is different from the
z-transform of Y (s) in Figures 3.20 (equation (3.113)): for the system in Figure 3.19,
Y (z) is the product of the z-transform of X(s) and the z-transform of G(s), while for
the system in Figure 3.20, Y (z) is the z-transform of the product G(s)X(s) (abbreviated
as GX(s)).
Example:
Suppose that in both Figures 3.19 and 3.20,
1 1
X(s) = (step input), and G(s) = .
s s+a
For the system in Figure 3.19, Y (z) = G(z)X(z).
   
1 1 1 1
G(z) = z = and X(z) = z = ,
s+a 1 − z e−aT
−1 s 1 − z −1
hence
1
Y (z) = .
(1 − z −1 e−aT ) (1 − z −1 )

For the system in Figure 3.20,


 
1
Y (z) = z [GX(s)] = z .
s(s + a)
Applying partial-fraction expansion to the expression on the right,
    
1 1 1 1 1 1
Y (z) = z − = −
a s s+a a 1 − z −1 1 − z −1 e−aT
 !
1 z −1 1 − e−aT
=
a (1 − z −1 ) (1 − z −1 e−aT )

Consider the system shown in Figure 3.21. It is assumed that all the samplers are
synchronized.
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 73

x(t) x∗ (t) u(t) u∗ (t) y(t) y ∗ (t)


δT G(s) δT H(s) δT
X(s) X ∗ (s) U (s) U (s)

Y (s) Y ∗ (s)

Figure 3.21:

From the figure, we can write

Y (s) = H(s)U ∗ (s). (3.114)

Taking the starred Laplace transform,

Y ∗ (s) = [H(s)U ∗ (s)]∗ = H ∗ (s)U ∗ (s). (3.115)

But
U (s) = G(s)X ∗ (s) =⇒ U ∗ (s) = [G(s)X ∗ (s)]∗ = G∗ (s)X ∗ (s), (3.116)
hence
PSfrag replacements
Y ∗ (s) = H ∗ (s)G∗ (s)X ∗ (s) (3.117)

or
Y (z) = H(z)G(z)X(z). (3.118)
The pulse transfer function of the block diagram shown in Figure 3.21 is

Y (z)
= H(z)G(z). (3.119)
X(z)

Consider the system shown in Figure 3.22. Again, it is assumed that all the samplers
are synchronized.

x(t) x∗ (t) u(t) y(t) y ∗ (t)


δT G(s) H(s) δT
X(s) X (s)∗
U (s) Y (s) Y ∗ (s)

Figure 3.22:

For this system,


Y (s) = G(s)H(s)X ∗ (s), (3.120)
hence

Y ∗ (s) = [G(s)H(s)X ∗ (s)]∗ = [G(s)H(s)]∗ X ∗ (s) = [GH(s)]∗ X ∗ (s) (3.121)

where GH(s) = G(s)H(s).


ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 74

It follows that
Y (z) = GH(z)X(z), (3.122)
and the pulse transfer function of the system is

Y (z)
= GH(z). (3.123)
X(z)

replacementswhere GH(z) = z [GH(s)] = z [G(s)H(s)].

3.6.2 Closed-loop systems

Consider the closed-loop system shown in Figure 3.23.

r(t) e(t) e∗ (t) y(t) y ∗ (t)


δT G(s) δT
R(s) + E(s) E (s) ∗ Y (s) Y ∗ (s)

H(s)

Figure 3.23:

From the figure,


E(s) = R(s) − H(s)Y (s). (3.124)
But
Y (s) = G(s)E ∗ (s), (3.125)
hence
E(s) = R(s) − H(s)G(s)E ∗ (s). (3.126)
Taking the starred Laplace transform,

E ∗ (s) = R∗ (s) − (GH(s))∗ E ∗ (s) (3.127)

hence
R∗ (s)
E ∗ (s) = . (3.128)
1 + (GH(s))∗
But Y (s) = G(s)E ∗ (s), so
Y ∗ (s) = G∗ (s)E ∗ (s). (3.129)
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 75

Substituting (3.128) into (3.129) gives


G∗ (s)R∗ (s) Y ∗ (s) G∗ (s)
Y ∗ (s) = =⇒ = . (3.130)
1 + (GH(s))∗ R∗ (s) 1 + (GH(s))∗
The pulse transfer function of the closed-loop system is therefore
Y (z) G(z)
rag replacements = (3.131)
R(z) 1 + GH(z)
where GH(z) = z [GH(s)] = z [G(s)H(s)]
Now consider the system shown in Figure 3.24. It is assumed that the samplers are
synchronized.

r(t) e(t) e∗ (t) y(t) y ∗ (t)


δT G(s) δT
R(s) + E(s) E (s)

Y (s) Y ∗ (s)

H(s)

Figure 3.24:

From the figure,


E(s) = R(s) − H(s)Y ∗ (s). (3.132)
But
Y (s) = G(s)E ∗ (s) =⇒ Y ∗ (s) = G∗ (s)E ∗ (s). (3.133)
hence

E(s) = R(s) − H(s)G∗ (s)E ∗ (s) =⇒ E ∗ (s) = R∗ (s) − H ∗ (s)G∗ (s)E ∗ (s). (3.134)

Making E ∗ (s) the subject of the formula gives


R∗ (s)
E ∗ (s) = . (3.135)
1 + G∗ (s)H ∗ (s)
Substituting (3.135) into (3.133) gives
G∗ (s)R∗ (s) Y ∗ (s) G∗ (s)
Y ∗ (s) = =⇒ = . (3.136)
1 + G∗ (s)H ∗ (s) R∗ (s) 1 + G∗ (s)H ∗ (s)
The pulse transfer function of the closed-loop system is therefore
Y (z) G(z)
= . (3.137)
R(z) 1 + G(z)H(z)
rag replacements

ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 76

3.6.3 Closed-loop systems with a Zero-Order-Hold

Consider the system shown in Figure 3.25.

r(t) e(t) e∗ (t) u(t) y(t) y ∗ (t)


rag replacements δT ZOH G(s) δT
R(s) + E(s) E ∗ (s) U (s) Y (s) Y ∗ (s)

Figure 3.25:

To obtain the pulse transfer function of the closed-loop system, we first replace the ZOH
block with its transfer function as shown in Figure 3.26.

r(t) e(t) e∗ (t) u(t) y(t) y ∗ (t)


1 − eT s G(s)
δT δT
R(s) + E(s) E (s)
∗ s U (s) Y (s) Y ∗ (s)

G1 (s)

Figure 3.26:

G1 (s) is the combined transfer function of the ZOH and G(s), and the z-transform of
G1 (s) can then obtained as described in Section 3.5.3, i.e.
 
 G(s)
G1 (z) = z [G1 (s)] = 1 − z −1
z (3.138)
s

The block-diagram of the closed-loop system reduces to what is shown in Figure 3.27.
The pulse transfer function of the closed-loop system is

Y (z) G1 (z)
= . (3.139)
R(z) 1 + G1 (z)
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 77

r(t) e(t) e∗ (t) y(t) y ∗ (t)


δT G1 (s) δT
R(s) + E(s) E ∗ (s) Y (s) Y ∗ (s)

Figure 3.27:

Exercise
For the system shown in Figure 3.25, suppose
K
G(s) =
s+1
where K is a constant.

(i) Show that 


Kz −1 1 − e−T
G1 (z) =
(1 − z −1 e−T )

(ii) Show that the pulse transfer function of the closed-loop system is

Y (z) Kz −1 1 − e−T
=
R(z) (1 − z −1 e−T ) + Kz −1 (1 − e−T )
replacements
Figure 3.28 shows a typical digital control system.
plant
r(t) e(t) sample e(kT ) u(kT )
D/A u(t) y(t)
and A/D digital
G(s)
+ hold converter controller converter

G(s)

Figure 3.28: typical digital control system

To obtain the pulse transfer function of the closed-loop system, the sample-and-hold-
A/D-converter block is replaced with an ideal sampler of sampling interval T , the digital
controller is represented using a transfer function C(z) (or C ∗ (s)). Unless otherwise
specified, the D/A converter is assumed to be a Zero-Order-Hold (ZOH).
The new block diagram of the system is shown in Figure 3.29.
The block diagram shown in Figure 3.29 can be further reduced as shown in Figure 3.30,
rag replacements

ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 78

R(s) E(s) E ∗ (s) 1 − eT s Y (s)


C ∗ (s) G(s)
input + δT s output

G1 (s)

Figure 3.29:

where
1 − eT s
G1 (s) = G(s).
s
PSfrag replacements

R(s) E(s) E ∗ (s) Y (s)


C (s)

G1 (s)
input + δT output

Figure 3.30:

From Figure 3.30,

Y (s) = E ∗ (s)C ∗ (s)G1 (s) =⇒ Y ∗ (s) = E ∗ (s)C ∗ (s)G∗1 (s). (3.140)

Also,

E(s) = R(s)−Y (s) =⇒ E ∗ (s) = R∗ (s)−Y ∗ (s) = R∗ (s)−E ∗ (s)C ∗ (s)G∗1 (s) (3.141)

hence
R∗ (s)
E ∗ (s) = . (3.142)
1 + G∗ (s)C ∗ (s)
Substituting (3.142) into (3.140) gives

R∗ (s)C ∗ (s)G∗1 (s) Y ∗ (s) C ∗ (s)G∗1 (s)


Y ∗ (s) = =⇒ = (3.143)
1 + G∗ (s)C ∗ (s) R∗ (s) 1 + G∗ (s)C ∗ (s)

hence the pulse transfer function of the closed-loop system is

Y (z) C(z)G1 (z)


= (3.144)
R(z) 1 + C(z)G1 (z)
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 79

Exercise
For the system shown in Figure 3.29, suppose
K
G(s) =
s+1
where K is a constant. Show that the pulse transfer function of the closed-loop system
is 
Y (z) Kz −1 1 − e−T
=
R(z) 1 − z −1 (K − Ke−T − e−T − 1) − z −2 e−T

3.7 Reconstruction of signals from sampled signals

3.7.1 Introduction

Sampling is made necessary because of the discrete nature of digital systems. Sampling
means that the continuous-time signal is replaced by a train of impulses. This handout
examines if it is possible to make an exact reconstruction the original signal from its
samples.

3.7.2 Shannon’s Sampling Theorem

Consider a signal which has a Fourier transform that is zero outside the interval (−B, B)
Hz (or (−2πB, 2πB) radians per second). Such a signal is said to be bandlimited to B
(see Figure 3.31).
PSfrag replacements

−B 0 B f (Hz)

Figure 3.31: Spectrum of a bandlimited signal

Shannon’s sampling theorem states that a signal bandlimited to B can be reconstructed


from its samples taken uniformly at a rate of not less than 2B samples per second.
To put it in another way, the theorem states that a bandlimited continuous-time signal
can be reconstructed from its samples provided that the sampling frequency is at least
twice as high as the highest frequency component in the continuous-time signal.
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 80

3.7.3 Proving Shannon’s theorem

Consider the continuous-time signal g(t) shown in Figure 3.32. Assume that this signal
is bandlimited.

g(t)

PSfrag replacements

0 t
PSfrag replacements
Figure 3.32: A bandlimited signal

Also consider a train of unit impulses shown in Figure 3.33. It can be written that

X
s(t) = δ(t − kT ). (3.145)
k=−∞

s(t)

−5T −4T −3T −2T −T 0 T 2T 3T 4T 5T t

Figure 3.33: A train of unit impulses

Multiplication of g(t) by the train of unit impulses s(t) yields the sampled signal gs (t)
shown in Figure 3.34. It follows that

X
gs (t) = g(t) · δ(t − kT ). (3.146)
k=−∞

We now obtain the Fourier transforms of the signals g(t), s(t) and gs (t) above. It was
assumed that the signal g(t) is bandlimited and we can assume without loss of gen-
erality that its Fourier transform G(ω) is bandlimited to 2πB radians per second, with
a maximum amplitude A, and its Fourier transform G(ω) can be drawn as shown in
Figure 3.35.
PSfrag replacements
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 81

gs (t)

0 T 2T 3T 4T 5T 6T t
PSfrag replacements
Figure 3.34: The sampled signal

−2πB 0 2πB ω(rad/sec)

Figure 3.35: Fourier transform of the bandlimited signal g(t)

The Fourier
PSfrag transform of the train of unit impulses (s(t)) is given by
replacements
∞ ∞
2π X 2πk 2π X
F (s(t)) = S(ω) = δ(ω − )= δ(ω − kωs ) (3.147)
T k=−∞ T T k=−∞

where ωs = 2π T
(this is derived using the theory of Fourier transforms for periodic sig-
nals). This signal is illustrated in Figure 3.36

S(ω)

−4ωs −3ωs −2ωs −ωs 0 ωs 2ωs 3ωs 4ωs ω

Figure 3.36: Fourier transform of a train of unit impulses

Recall that multiplication in the time domain is equivalent to convolution in the fre-
quency domain. Formally, this can written as
1 1
Z ∞
F
gs (t) = s(t) · g(t) ←→ Gs (ω) = [S(ω) ∗ G(ω)] = S(ω − τ )G(τ )dτ (3.148)
2π 2π −∞

The signal Gs (ω) is shown in Figure 3.37


As can be seen from Figure 3.37, the frequency spectrum of the impulse sampled signal
is reproduced an infinite number of times and is attenuated by a factor of T1 . It can
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 82

complementary primary complementary


components component components
A
T

−2ωs −ωs −2πB 0 2πB ωs 2ωs ω(rad/sec)

Figure 3.37: Fourier transform of gs (t)

written that ∞
A X
Gs (ω) = G(ω − kωs ) (3.149)
T k=−∞

where ωs = T
.
If ws > 2(2πB), there is no overlap in the shifted replicas of G(ω) and hence g(t) can
recovered exactly from gs (t) by means of a low pass filter with a gain T and a cut-
off frequency ωc greater than 2πB and less than (ωs − 2πB). This proves Shannon’s
sampling theorem.
1
The maximum allowable sampling interval Ts = 2B is known as the Nyquist interval.
The corresponding sampling rate (2B samples per second) is known as the Nyquist
sampling rate.

3.7.4 The
PSfrag ideal low pass filter
replacements

The ideal low pass filter required for Shannon reconstruction is shown in Figure 3.38.

H(ω)
T primary
component

−ωc 2πB 0 −2πB ωc ω

Figure 3.38: Frequency spectrum of an ideal low pass filter

This ideal low pass filter shown in Figure 3.38 has a transfer function
(
T |ω| ≤ ωc
H(ω) = (3.150)
0 |ω| > ωc
where ωc is the cut-off frequency of the filter. This filter will perform exact interpolation
of the signal between the sampling instants.
Performing the inverse Fourier transform on H(ω) gives
Z ωc
1 1
Z ∞
jωt
T ejωt dω

h(t) = H(ω)e dω = (3.151)
2π −∞ 2π −ωc
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 83

T ejωc t e−jωc t T ejωc t − e−jωc t


   
T T ωc sin(ωc t)
= − = = sin(ωc t) = (3.152)
2π jt jt πt 2j πt π ωc t
From equation (3.152), it can be seen that h(t) is a sinc function, and its profile is shown
in Figure 3.39.
1.2

1
T ωc
π
0.8

0.6

0.4

PSfrag replacements 0.2

0 t
−0.2

−0.4
−5 −4 −3 −2 −1 0 1 2 3 4 5

Figure 3.39: Impulse response of the ideal low pass filter

This signal is non-causal and it is physically non-realizable.


As stated in the previous section, the cut-off frequency ωc should be greater than 2πB
and less than (ωs − 2πB). In practice, ωc is usually selected at the midpoint between
these two frequencies, i.e. ωc = ω2s = Tπ .

3.7.5 Frequency response of the Zero-Order-Hold

From Section 3.5.2, it was shown that the transfer function of the ZOH is given by
1 − e−sT
Gh (s) =.
s
Replacing s with jω in the above equation gives
jωT
! jωT jωT
!
1 − e−T s
−jωT −jωT
1 − e−jωT jωT e 2 −e 2 2e 2 e 2 −e 2
Gh (jω) = = = 2e 2 =
s s=jω jω 2jω ω 2j
  
ωT
jωT   sin
2e 2 ωT jωT  2 
= sin = Te 2   . (3.153)
ω 2  ωT 
2
The magnitude response of the ZOH is given by
    
ωT ωT

sin  sin 2
jωT  2  
|Gh (jω)| = T e 2   =T   . (3.154)
 ωT   ωT 
2 2
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 84

The magnitude curve of the ZOH is shown in Figure 3.40. The frequency response of
the ideal low pass filter is drawn in dotted lines for comparison.

Gh (ω)
1 T

0.8

0.6
0.637T
PSfrag replacements

0.4

0.2 0.212T
0.127T
0

− 6π − 4π − 2π 0 2π 4π 6π ω
T T T T T T
−0.2
−20 −15 −10 −5 0 5 10 15 20

Figure 3.40: Frequency response of the Zero-Order-Hold

From the figure, it can be seen that the frequency response of the Zero-Order-Hold is
quite different from the frequency response of the ideal low pass filter. In particular,
the spectrum of the ZOH does not go to zero beyond the frequency ω = ω2s = Tπ , and
this causes distortion in the signal reconstructed using the ZOH. However, the ZOH is
physically realizable and is simple to implement and for these reasons it is commonly
used in digital control systems.

3.7.6 Aliasing or frequency folding

When sampling is done at frequencies higher than 2B samples per second, the frequency
spectrum of a sampled bandlimited signal is as shown in Figure 3.37. However, if the
sampling frequency is less than 2B samples per second, the frequency spectrum of the
sampled signal looks as shown in Figure 3.41.
In this case, the replicated spectra of G(ω) overlap and g(t) is no longer recoverable by a
low pass filter. This overlapping effect is known as aliasing. (One can think of aliasing as
follows: from Figure 3.41, it can be seen that the high frequency tail of the component
2 overlaps with a lower frequency part of component 1, hence we say that the higher
frequencies of component 2 take on an alias of a lower frequency of component 1).
PSfrag replacements

ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 85

A
T

2 1

−2ωs −ωs 0 ωs 2ωs ω(rad/sec)

Figure 3.41: Aliasing

When aliasing takes place, the reconstructed signal will have distortion for two reasons:

• Loss of the high frequency tail of Gs (ω) beyond |ω| > ωs /2.

• The same high frequency tail appears inverted or folded onto the spectrum at
frequencies |ω| < ωs /2.

This tail inversion is known as spectral folding and hence aliasing is sometimes referred
to as frequency folding.
In practice, aliasing always occurs regardless of the sampling frequency since bandlim-
ited signals do not exist in reality. It can be shown that if a signal is bandlimited, it
cannot be time-limited and vice-versa. All physical signals exist over a finite time dura-
tion hence they are not bandlimited.
To reduce the effects of aliasing, the following measures are usually adopted:

• Sampling at a frequency much higher than 2B samples per second.

• Filtering the signals before sampling (this reduces the high frequency tail).

In practice, the rule of thumb is to choose the sampling time T such that, either
1
= 30 × closed-loop bandwidth (3.155)
T
or
smallest time constant of system
T = . (3.156)
20
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 86

3.8 Discretization of controllers

3.8.1 Introduction

Discrete-time controller design methods can be broadly divided into two categories:

Direct methods: The discrete-time controller is designed directly in the z-plane.

Indirect methods: A continuous-time controller is first designed and then it is con-


verted to discrete-time form (i.e. the controller is discretized).

Indirect methods are more popular for the following reasons:

 There are many well-established methods of continuous-time controller design


and the designer is more likely to have more experience in these methods.

 Direct discrete-time controller design requires the pre-determination of the sam-


pling time T . To determine the sampling time, one needs to know the closed-loop
bandwidth, which in turn can only be determined if the continuous controller is
designed first.

 Discrete-time controller design only considers the system behaviour at the sam-
pling instants. Problems of inter-sampling behaviour can only be solved if the
continuous-time controller is designed first and then discretized.

 Physical insight is more easily retained when a continuous-time controller is first


designed and then discretized.

In the discretization, one may wish to preserve some of the following properties:

 impulse, step or ramp response.

 number of poles and zeros.

 d.c. gain.

 gain and phase margin.

 closed loop transfer function.

 closed loop stability.

Not all these properties can be preserved at the same time after discretization and the
choice of the discretization method is normally based on what the designer wants to
preserve.
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 87

3.8.2 Discretization methods based on numerical integration

Introduction

The fundamental concept in these methods is to represent the continuous-time con-


troller as a differential equation and then to derive a difference equation whose solution
is an approximation to the differential equation.
Consider a continuous controller described by the transfer function
U (s) b
= (3.157)
E(s) s+a

Cross-multiplying gives
sU (s) + aU (s) = bE(s). (3.158)
Converting (3.158) into time domain gives the differential equation:

u̇(t) + au(t) = be(t), (3.159)

hence we can write Z t


u(t) = (−au(t) + be(t)) dt. (3.160)
0
At time t = kT , Z kT
u(kT ) = (−au(t) + be(t)) dt. (3.161)
0

Similarly, at time t = kT − T ,
Z kT −T
u(kT − T ) = (−au(t) + be(t)) dt. (3.162)
0

Subtracting (3.162) from (3.161) gives


Z kT
u(kT ) = u(kT − T ) + (−au(t) + be(t)) dt. (3.163)
KT −T

The numerical integration methods depend on how the integral on the right of equation
(3.163) is approximated, i.e. how the area under the curve −au(t) + be(t) in the region
between t = (kT − T ) and t = kT is approximated.

Forward rectangular rule

By the forward rectangular rule, the area between t = (kT − T ) and t = kT is approxi-
mated by
Z kT
(−au(t) + be(t)) dt ≈ T (−au(kT − T ) + be(kT − T )) . (3.164)
KT −T
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 88


 
−au(t)+be(t)
  
      
         
         
        
   
kT−T kT t

Figure 3.42: Forward rectangular rule

(see Figure 3.42). Substituting this approximation in (3.163) gives

u(kT ) = u(kT −T )+T (−au(kT − T ) + be(kT − T )) = (1−aT )u(kT −T )+bT e(kT +T ).


(3.165)
Taking z-transforms,

U (z) = (1 − aT ) z −1 U (z) + bT z −1 E(z) (3.166)

hence
U (z) bT z −1 b
= = . (3.167)
E(z) 1 − z + aT z
−1 −1 1−z −1
+a
T z −1
Comparing (3.167) with (3.157) shows that

1 − z −1 z−1
s= = . (3.168)
Tz −1 T

Equation (3.168) is known as the forward rectangular rule.


From (3.168),
z = sT + 1, (3.169)
hence the jω axis on the s-plane is mapped to the line z = jωT + 1. The stable region in
the s-plane is mapped to the left of the line z = jωT + 1, as shown in Figure 3.43. It is
thus evident that a stable s-plane controller can be mapped outside the unit circle in the
z-plane, resulting in an unstable system. For this reason, this method is not commonly
used.
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 89
       
       z−plane

       z=jwT+1

       


    0  
       unit
circle

      


Figure 3.43:

−au(t)+be(t) &&&&'' #"#" #"#" #"#"


 
  
   
& & ' 
! 
! 
! ! #"#" #"#" #"#"
             !! !!!!!! "# "# "#
            !! !!!!!! #"#" #"#" #"#"
            !! !!!!!! #"#" #"#" #"#"
           !! !!!!!! "##" "##" "##"
%$ %$%$ %$ %$ %$  %$  %$  %$   %$  %$  %$  %$ %$ %$ #" %$ #" %$ #" %$ %$ %$t
kT−T kT

Figure 3.44: Backward rectangular rule

Backward rectangular rule

By the backward rectangular rule, the area between t = (kT − T ) and t = kT is approx-
imated by
Z kT
(−au(t) + be(t)) dt ≈ T (−au(kT ) + be(kT )) . (3.170)
KT −T

(see Figure 3.44). Substituting this approximation in (3.163) gives

u(kT ) = u(kT − T ) + T (−au(kT ) + be(kT )) . (3.171)

Taking z-transforms,

U (z) = z −1 U (z) − aT U (z) + bT E(z) = z −1 − aT U (z) + bT E(z) (3.172)

hence
U (z) bT b
= = . (3.173)
E(z) 1 − z −1 + aT 1−z −1
+a
T
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 90

Comparing (3.173) with (3.157) shows that


1 − z −1 z−1
s= = . (3.174)
T Tz
Equation (3.174) is known as the backward rectangular rule.
From (3.174),
   
1 1 1 1 1 1 1 + sT
z= = + − = − . (3.175)
1 − sT 2 1 − sT 2 2 2 1 − sT
Substituting s = jω in (3.175) gives
 
1 1 1 + jωT
z= − , (3.176)
2 2 1 − jωT
hence
1 1
z−= . (3.177)
2 2
Equation (3.177) describes a circle of radius 12 , centered at ( 12 , 0).

z−plane

*)(*)( *)(*)( *)(*)( *)(*)( *)(*)( *(*(


*()*()*)(*)( *)(*)( *)(*)( *)(*)( *(*(
0 *()*)( *)(*)( *)(*)( *)(*)( *)(*)( *(*( 1

Figure 3.45:

The entire left half of the s-plane is mapped into a circle of radius 21 , centered at ( 12 , 0) in
the z-plane, (see Figure 3.45), hence stable s-plane controllers produce stable z-plane
controllers (it is even possible for an unstable s-plane controller to be mapped into
the stable region in the z-plane). However, since the entire left half of the s-plane is
squeezed into a circle within the unit circle in the z-plane, there is a lot of distortion in
the transient and frequency response of the controller obtained using this method. To
minimize the distortion, a small sampling time should be used.

The trapezoidal rule

By the trapezoidal rule, the area between t = (kT − T ) and t = kT is approximated by


Z kT
T
(−au(t) + be(t)) dt ≈ (−au(kT − T ) + be(kT − T ) − au(kT ) + be(kT )) .
KT −T 2
(3.178)
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 91

+,,+ 8778 8778 8778


0/0/ 0/0/ ,+,+ 212121 0/0/ 212121 212121 43 212121 43 43 5665 43 5665 5665 6565
−au(t)+be(t)
87878787 8787
.-.- .-.- .-.- 0/0/0/0/ 2121 0/0/ 2121 2121 4343 2121 4343 4343 6565 4343 6565 6565 6565 87878787 8787
.-.-.-.- .-.- 0/0/0/0/ 2121 0/0/ 2121 2121 4343 2121 4343 4343 6565 4343 6565 6565 6565 7887 7887 7887
.-.-.-.- .-.- 0/0/ 0/0/ 2121 0/0/ 2121 2121 4343 2121 4343 4343 6565 4343 6565 6565 6565 87878787 8787
.-.- .-.- .-.- 0/0/ 21 0/ 21 21 43 21 43 43 65 43 65 65 65 8778 8778 8778
kT−T kT t

Figure 3.46: The trapezoidal rule

(see Figure 3.46). Substituting this approximation in (3.163) gives


T
u(kT ) = u(kT − T ) + (−au(kT − T ) + be(kT − T ) − au(kT ) + be(kT )) . (3.179)
2
By taking z-transforms and rearranging, it can be shown that

U (z) b
=   , (3.180)
E(z) 2 1 − z −1
+a
T 1 + z −1

hence    
2 1 − z −1 2 z−1
s= = (3.181)
T 1 + z −1 T z+1
Equation (3.181) is known by several names: trapezoidal rule, Tustin’s method, Tustin’s
transformation or bilinear transformation.
For this method,
sT
1+
z= 2 , (3.182)
Ts
1−
2
hence for s = jω, |z| = 1. The left half s-plane is mapped into the unit circle in the
z-plane, hence stable s-plane controllers will produce stable z-plane controllers.

3.8.3 Impulse and Step invariance methods

Impulse invariance method

For this method, the impulse response of the continuous-time controller is equal to the
impulse response of the discrete-time controller at the sampling instants.
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 92

Given a continuous controller G(s), the discrete-time equivalent using the impulse in-
variance method is given by

G(z) = T z(G(s)). (3.183)

Step invariance method

For the step-invariance method, the step responses of the continuous and the discrete-
time controllers agree at the sampling instants.
Given a continuous-time controller G(s), the discrete-equivalent using the step-invariance
method is given by

1 − e−T s
   
 G(s)
G(z) = z G(s) = 1−z −1
z (3.184)
s s

(compare equation (3.184) with equation (3.95) – the controller discretized using this
method is sometimes referred to as the ZOH equivalent).

3.8.4 The matched pole-zero mapping method

Using this method, the poles and finite zeros of the continuous controller are mapped
onto the z-plane using the rule z = esT . This method is easy to apply if the transfer
function of the continuous controller is given in poles and zeros.
The complete procedure of discretizing a continuous-time controller using the matched
pole-zero mapping method is as follows:

1. The controller G(s) must be factored before this method can be applied.

2. The poles of G(s) are mapped onto the z-plane according to the relation z = esT ,
e.g. a pole at s = −a is mapped to z = e−aT .

3. All finite zeros are mapped according to the relation z = esT .

4. Each of the zeros at infinity are mapped at z = −1.

5. The gain of the discrete-time controller is adjusted to match the gain of the con-
tinuous controller.
If the controller is a low-pass filter, the gain of the discrete-time controller at z = 1
should be the same as the gain of the continuous controller at s = 0.
For a high pass filter, match the gains at z = −1 and s = ∞.
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 93

Example
Use the matched pole-zero mapping method to obtain the discrete-time equivalent of
the continuous-time controller
a
G(s) = . (3.185)
s+a
Solution
The given controller is a low-pass filter with one pole at s = −a and one zero at infinity.
The pole and zero are mapped to z = e−aT and z = −1 respectively, giving a discrete
time equivalent
z+1
G(z) = K , (3.186)
z − e−aT
where K is the gain.
Since this is a low-pass filter,
G(z)|z=1 = G(s)|s=0 , (3.187)
which gives
1 − e−aT
K= , (3.188)
2
hence 
1 − e−aT (z + 1)
G(z) = . (3.189)
2 (z − e−aT )

Exercise
Use the matched pole-zero mapping method to show that the discrete-time equivalent
of the continuous-time controller
s+b
G(s) = , (3.190)
s+a
where our interest is in the low-frequency range, is given by
 
b 1 − e−aT z − e−bT
G(z) = (3.191)
a (1 − e−bT ) (z − e−aT )
Example
Use the matched pole-zero mapping method to obtain the discrete-time equivalent of
the high-pass filter
s
G(s) = . (3.192)
s+a
Solution
The pole at s = −a and the zero at s = 0 are mapped to z = e−aT and z = 1 respectively,
giving a discrete time equivalent
z−1
G(z) = K , (3.193)
z − e−aT
where K is the gain.
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 94

For a high-pass filter,


G(z)|z=−1 = G(s)|s=∞ , (3.194)
which gives
1 + e−aT
K= , (3.195)
2
hence  
1 + e−aT (z − 1) 1 + e−aT (1 − z −1 )
G(z) = = . (3.196)
2 (z − e−aT ) 2 (1 − e−aT z −1 )

3.9 State space analysis

3.9.1 Introduction

The dynamics of a discrete-time system can be written as

x(k + 1) = f (x(k), u(k), k) (3.197a)


y(k) = g (x(k), u(k), k) (3.197b)

where (3.197a) is the state equation while (3.197b) is the output equation.
If these equations can be linearized, then it can be written that

x(k + 1) = G(k)x(k) + H(k)u(k) (3.198a)


y(k) = C(k)x(k) + D(k)u(k) (3.198b)

which is linear-time varying system.


If the system is linear-time-invariant, the above equations become

x(k + 1) = Gx(k) + Hu(k) (3.199a)


y(k) = Cx(k) + Du(k) (3.199b)

where G, H, C and D are constant matrices. For a Multiple-Input-Multiple-Output


(MIMO) system with n states, r inputs and m outputs, the dimensions of these matrices
are as follows:

G: State matrix, dimensions n × n

H: Input matrix, dimensions n × r

C: Output matrix, dimensions m × n

D: Direct transmission matrix, also referred to as feed-through matrix, dimensions


m×r
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 95

Strictly speaking, equation (3.199) should be written as

x(kT + T ) = Gx(kT ) + Hu(kT ) (3.200a)


y(kT ) = Cx(kT ) + Du(kT ). (3.200b)

The state equation (3.200a) gives the next state x(kT + T ) (or x(k + 1)) given the
current state x(kT ) and input u(kT ), while the output equation (3.200b) gives the
current output y(kT ) given the current state and input.
In some books, equation (3.200) is written as

xk+1 = Gxk + Huk (3.201a)


yk = Cxk + Duk (3.201b)
rag replacements
The system given by equation (3.199) is illustrated in Figure 3.47.

x(k) + y(k)
u(k) x(k+1) +
H z I
−1 C
+
+

Figure 3.47: Discrete-time LTI state space model

Figure 3.47 is similar to the one for a continuous-time LTI state-space model encoun-
tered earlier, the only difference being that the integrator in the continuous-time model
is replaced by the z −1 operator in the discrete-time model. The z −1 operator is known as
the unit delay operator and it is analogous to the integrator in continuous-time models.

3.9.2 Canonical forms

Introduction

Just as in continuous-time systems, discrete time models may be expressed in the canon-
ical forms, namely

• the control canonical form,

• the observable canonical form, and


ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 96

• the diagonal (Jordan) canonical form.

Consider a discrete-time system described by the difference equation:

y(k) + a1 y(k − 1) + a0 y(k − 2) = b2 u(k) + b1 u(k − 1) + b0 u(k − 2) (3.202)

Taking the z transform of (3.202) gives

Y (z) + a1 z −1 Y (z) + a0 z −2 Y (z) = b2 U (z) + b1 z −1 U (z) + b0 z −2 U (z), (3.203)

which gives a pulse transfer function

Y (z) b2 + b1 z −1 + b0 z −2
= . (3.204)
U (z) 1 + a1 z −1 + a0 z −2

Control canonical form


rag replacements

The pulse transfer function (3.204) can be represented as shown in Figure 3.48

U (z) 1 X(z) Y (z)


1 + a1 z −1 + a0 z −2 b2 + b 1 z −1
+ b0 z −2

Figure 3.48:

From Figure 3.48,


X(z) 1
= (3.205)
U (z) 1 + a1 z + a0 z −2
−1

hence
X(z) = U (z) − a1 z −1 X(z) + a0 z −2 X(z). (3.206)
Also,
Y (z) = b2 X(z) + b1 z −1 X(z) + b0 z −2 X(z). (3.207)
Equations (3.206) and (3.207) can be modelled as shown in Figure 3.49.
With the states labelled as shown in Figure 3.49, the state-space representation in con-
trol canonical form is:
" # " #" # " #
x1 (k + 1) 0 1 x1 (k) 0
= + u(k) (3.208a)
x2 (k + 1) −a0 −a1 x2 (k) 1
" #
h i x (k)
1
y(k) = b0 − a 0 b0 b1 − a 1 b2 + [b2 ] u(k) (3.208b)
x2 (k)
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 97

PSfrag replacements b2

b1

+
U (z) X(z) x2 x1 +
+ Y (z)
+ + z −1 z −1 b0
+

−a1

−a0

Figure 3.49: Model in control canonical form

Observable canonical form

The pulse transfer function (3.204) can be written as


Y (z) b2 z 2 + b 1 z + b 0
= 2 . (3.209)
U (z) z + a1 z + a 0
Cross-multiplication gives
z 2 Y (z) + a1 zY (z) + a0 Y (z) = b2 z 2 U (z) + b1 zU (z) + b0 U (z). (3.210)
The observable canonical model is shown in Figure 3.50 where
f1 (z) = z 2 Y (z) + a1 zY (z) − b2 z 2 U (z) − b1 zU (z)
and
f2 (z) = zY (z) − b2 zU (z).

With the states labelled as shown in Figure 3.50, the state-space representation in ob-
servable canonical form is:
" # " #" # " #
x1 (k + 1) 0 −a0 x1 (k) b0 − a 0 b0
= + u(k) (3.211a)
x2 (k + 1) 1 −a1 x2 (k) b1 − a 1 b2
" #
h i x (k)
1
y(k) = 0 1 + [b2 ] u(k) (3.211b)
x2 (k)

Diagonal canonical form

Consider the pulse transfer function


Y (z) bm z m + bm−1 z m−1 + bm−2 z m−2 + · · · + b1 z + b0
= (3.212)
U (z) z n + an−1 z n−1 + an−2 z n−2 + · · · + a1 z + a0
PSfrag replacements

ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 98

U (z)

b0 b1 b2

+ f1 (z) + f2 (z) + Y (z)


x1 + x2 +
z −1
z −1
+ +

−a0 −a1

Figure 3.50: Model in observable canonical form

where n > m (a strictly proper pulse transfer function). Assuming that the characteristic
equation of the pulse transfer function (3.212) has n distinct poles p1 , p2 , . . . , pn , the
pulse transfer function can be expanded using partial fraction expansion as

Y (z) c1 c2 cn
= + +···+ . (3.213)
U (z) z − p1 z − p2 z − pn
PSfrag replacements
A block diagram of the system is shown in Figure 3.51.

1 X1 (z)
c1
z − p1

U (z) X2 (z) + Y (z)


1 c2
z − p2 +
+

1 Xn (z)
cn
z − pn

Figure 3.51: Model in diagonal canonical form

From Figure 3.51, we can write that

X1 (z) 1
= . (3.214)
U (z) z − p1

Cross-multiplying gives
zX1 (z) − p1 X1 (z) = U (z). (3.215)
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 99

Converting (3.215) into time domain gives

x1 (k + 1) = p1 x1 (k) + u(k). (3.216)

Similarly
x2 (k + 1) = p2 x2 (k) + u(k)
x3 (k + 1) = p3 x3 (k) + u(k)
.. . (3.217)
.
xn (k + 1) = pn xn (k) + u(k)
The output equation is

Y (z) = c1 X1 (z) + c2 X2 (z) + · · · + cn Xn (z) (3.218)

which is written in time-domain as

y(k) = c1 x1 (k) + c2 x2 (k) + · · · + cn xn (k). (3.219)

Equations (3.216), (3.217) and (3.219) can be written in vector matrix as


      
x1 (k + 1) p1 0 · · · 0 x1 (k) 1
 x2 (k + 1)   0 p · · · 0   x2 (k)   1 
      
 =  . .2
..  . . . . . ..   ..  +  ..  u(k) (3.220a)
    
.  . . .  .   . 
 
 
xn (k + 1) 0 0 · · · pn xn (k) 1
 
x1 (k)
i x2 (k) 
h 

y(k) = c1 c2 · · · cn  ..  (3.220b)
. 
 

xn (k)

which is the system representation in diagonal canonical form.

3.9.3 Solving discrete-time state space equations

Time domain approach

In general, discrete-time equations are easier to solve than differential equations be-
cause difference equations can be solved easily by means of a recursion procedure. The
recursion procedure can easily be programmed for solution with a digital computer.
Consider the LTI discrete-time system

x(k + 1) = Gx(k) + Hu(k) (3.221a)


y(k) = Cx(k) + Du(k). (3.221b)
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 100

Equation (3.221a) can be solved for any positive integer k using a recursion procedure
by substituting k = 0, k = 1, k = 2 and so on.
Substituting k = 0 in equation (3.221a) gives

x(1) = Gx(0) + Hu(0). (3.222)

For k = 1,
x(2) = Gx(1) + Hu(1). (3.223)
Substituting (3.222) in (3.223) gives

x(2) = G2 x(0) + GHu(0) + Hu(1). (3.224)

Similarly
x(3) = G3 x(0) + G2 Hu(0) + GHu(1) + Hu(2) (3.225)
and so on. By repeating this procedure, we obtain
k−1
X
k
x(k) = G x(0) + Gk−j−1 Hu(j). (3.226)
j=0

x(k) consists of two parts, the term Gk x(0) which is the contribution by the initial state
Pk−1 k−j−1
x(0), and j=0 G Hu(j) which is the contribution by the input u(j). The term
G x(0) is also referred to as the zero input response while the term k−1
k k−j−1
P
j=0 G Hu(j)
is referred to as the zero initial state response. The continuous-time analogy to equation
(3.226) is Z t
x(t) = eAt x(0) + eA(t−τ ) Bu(τ )dτ (3.227)
0
Rt
where eAt x(0) is the zero input response and 0
eA(t−τ ) Bu(τ )dτ is the zero initial state
response.
Substituting (3.226) in (3.221b) gives
k−1
X
k
y(k) = CG x(0) + C Gk−j−1 Hu(j) + Du(k). (3.228)
j=0

State transition matrix

For discrete-time systems, the state transition matrix is defined as

Ψ = Gk . (3.229)

(This is analogous to the matrix exponential eAt for continuous-time systems).


ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 101

z transform approach

Consider the state equation of an LTI discrete-time system

x(k + 1) = Gx(k) + Hu(k). (3.230)

Taking the z transform on both sides of (3.230) gives

zX(z) − zx(0) = GX(z) + HU(z). (3.231)

Rearranging,
(zI − G) X(z) = zx(0) + HU(z). (3.232)
Pre-multiplying both sides of (3.232) by (zI − G)−1 gives

X(z) = (zI − G)−1 zx(0) + (zI − G)−1 HU(z). (3.233)

Comparing equations (3.233) and (3.226), we obtain

Gk = z−1 (zI − G)−1 z


 
(3.234)

and
k−1
X
Gk−j−1 Hu(j) = z−1 (zI − G)−1 HU(z) .
 
(3.235)
j=0

Equation (3.234) is used to compute the state transition matrix Gk .


Example
Obtain the state transition matrix for the system described by the state equation
" # " #" # " #
x1 (k + 1) 0 1 x1 (k) 1
= + u(k).
x2 (k + 1) −0.16 −1 x2 (k) 1

Solution
Gk = z−1 (zI − G)−1 z .
 

Obtaining Gk involves finding the inverse z-transform of (zI − G)−1 z. Recall that if
a function X(z) has one or more zeros at the origin, it is easier to obtain the inverse
z-transform by first expanding X(z)z
into partial fractions then multiplying the result by
z (see Section 3.2.4). For that reason, we first expand (zI − G)−1 into partial fractions.
" #−1
z −1
(zI − G)−1 =
0.16 z + 1
z+1 1
 
 (z + 0.2)(z + 0.8) (z + 0.2)(z + 0.8) 
 
= 


 −0.16 z 
(z + 0.2)(z + 0.8) (z + 0.2)(z + 0.8)
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 102

4 1 5 5
 
3 3 3 3
 z + 0.2 − z + 0.8 −
z + 0.2 z + 0.8 
 
=
.

 0.8 0.8 1 4 
− 3 + 3
− 3 + 3
z + 0.2 z + 0.8 z + 0.2 z + 0.8
Multiplying the above result by z gives
4 1 5 5
 
3
z 3
z 3
z 3
z

 z + 0.2 z + 0.8 −
 z + 0.2 z + 0.8 

−1
(zI − G) z =  
,

0.8 0.8 1 4
 z z z z 
− 3 + 3 − 3 + 3
z + 0.2 z + 0.8 z + 0.2 z + 0.8
hence
4 1 5 5
 
(−0.2)k − (−0.8)k k
(−0.2) − (−0.8) k
  3 3 3 3 
Gk = z−1 (zI − G)−1 z = 
  

 0.8 0.8 1 4 
− (−0.2)k + (−0.8)k − (−0.2)k + (−0.8)k
3 3 3 3
which is the state transition matrix.

3.9.4 Obtaining pulse transfer function from state space descrip-


tion

Consider a SISO LTI system described by the equations

x(k + 1) = Gx(k) + Hu(k) (3.236a)


y(k) = Cx(k) + Du(k). (3.236b)

Taking the z-transform of these equations gives

zX(z) − zx(0) = GX(z) + HU (z) (3.237a)


Y (z) = CX(z) + DU (z). (3.237b)

Recall that the transfer function relates the inputs to the outputs for zero initial condi-
tions hence we can ignore x(0) in equation (3.237a) to write

X(z) = (zI − G)−1 HU (z). (3.238)

Substituting (3.238) into (3.237b) gives


 
Y (z) = C (zI − G)−1 H + D U (z). (3.239)

Since this is a single-input-single-output system, the pulse transfer function is given by


Y (z)
= C (zI − G)−1 H + D. (3.240)
U (z)
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 103

3.9.5 Discretization of continuous-time state space equations

In this section, we obtain the state space representation of a continuous-time system


preceded by a Zero-Order-Hold as shown in Figure 3.52.
PSfrag replacements ZOH System
−T s
 
1−e
G(s)
s

Figure 3.52: Plant preceded by a ZOH

Assume the plant is a SISO LTI continuous-time system

ẋ(t) = Ax(t) + Bu(t) (3.241a)


y(t) = Cx(t) + Du(t). (3.241b)

The solution to the state equation (3.241a) is


Z t
At
x(t) = e x(0) + eA(t−τ ) Bu(τ )dτ (3.242)
0

At time t = kT + T , equation (3.242) becomes


Z kT +T
A(kT +T )
x(kT + T ) = e x(0) + eA(kT +T −τ ) Bu(τ )dτ. (3.243)
0

Similarly for time t = kT ,


Z kT
x(kT ) = e AkT
x(0) + eA(kT −τ ) Bu(τ )dτ. (3.244)
0

Pre-multiplying (3.244) by eAT gives


Z kT
A(kT +T )
e AT
x(kT ) = e x(0) + eA(kT +T −τ ) Bu(τ )dτ. (3.245)
0

Subtracting (3.245) from (3.243) gives


Z kT +T
x(kT + T ) = e AT
x(kT ) + eA(kT +T −τ ) Bu(τ )dτ. (3.246)
kT

The input to the plant comes from the ZOH. Recall that the output of the ZOH is constant
between sampling instants. Between the sampling instant kT and kT + T , the output of
the ZOH is u(kT ). Substituting this in (3.246) above gives
Z kT +T
AT
x(kT + T ) = e x(kT ) + eA(kT +T −τ ) Bu(kT )dτ. (3.247)
kT
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 104

Let kT + T − τ = λ. Then dτ = −dλ. If τ = kT , λ = T , and for τ = kT + T , λ = 0.


Making these substitutions in (3.247) above gives
Z 0
AT
x(kT + T ) = e x(kT ) − eAλ Bu(kT )dλ. (3.248)
T

or Z T
AT
x(kT + T ) = e x(kT ) + eAλ Bu(kT )dλ. (3.249)
0

At time t = kT , the output equation can be written as

y(kT ) = Cx(kT ) + Du(kT ). (3.250)

Comparing (3.249) with equation (3.200a), it follows that

G = eAT (3.251)

and Z T
H= eAλ Bdλ. (3.252)
0

Equation (3.251) is used to compute the state matrix of the discrete-time system shown
in Figure 3.52, while equation (3.252) is used to determine the input matrix. As can
be seen from equation (3.250), the output and the direct transmission matrices C and
D respectively remain unchanged. This model is known as the ZOH equivalent of the
continuous-time plant. The discrete-time model given by equations (3.249) and (3.250)
is also the step-invariant discrete-time equivalent of the continuous model (3.241).
In MATLAB, the ZOH equivalent of a continuous-time plant is obtained using the com-
1
mand c2d. Figure 3.53 shows the step response of the continuous-time system s2 +0.4s+1
(plotted using the continuous line), and the step response for the ZOH equivalent dis-
crete model (plotted with the asterisks). As can be seen from the figure, the two plots
agree at the sampling instants.
ELC 514E: STATE SPACE DESIGN AND DIGITAL CONTROL by Dr.-Ing. S. I. Kamau 105

1.6

1.4

1.2

0.8

0.6

0.4

0.2

0
0 2 4 6 8 10 12 14 16 18 20

Figure 3.53: Step response

The MATLAB program used to generate Figure 3.53 is listed below:


num = 1; %numerator polynomial
den = [1 0.4 1]; %denominator polynomial
[A,B,C,D] = tf2ss(num,den); %conversion to state-space form
sys1 = ss(A,B,C,D); %define a state-space object
t = 0:0.1:20; %time scale for continuous model
y = step(sys1,t); %step response for continuous model

Ts = 1; %sampling time
sysd = c2d(sys1,Ts); %discretization - ZOH equivalent
td =0:Ts:20; %sampling times for discrete model
N = length(td); %total number of sample points
AD = sysd.A; %state matrix for discrete model
BD = sysd.B; %input matrix for discrete model
CD = sysd.C; %output matrix for discrete model
DD = sysd.D; %feedthrough matrix for discrete model
yd = dstep(AD,BD,CD,DD,1,N); %step response for discrete model

figure
plot(t,y,td,yd,’*’),grid, %plotting the results
=== × === × === × === × === × === × === × === × === × === × ===
Bibliography

[1] K. Ogata, Discrete-Time Control Systems, Prentice Hall, 1987.

[2] G. F. Franklin, J. D. Powell, M. Workman, Digital Control of Dynamic Systems,


Third Edition, Pearson Education, 2000.

[3] G. Franklin, J. D. Powell, A. Emami-Naeini, Feedback Control of Dynamic Sys-


tems, Addison-Wesley, 1997.

[4] N. S. Nise, Control Systems Engineering, Fourth Edition, John Wiley & Sons
Inc., 2004.

SK ’07
T HIS DOCU MEN T WAS PREPARED U SIN G T HE LATEX
WORD − PROCESSOR RU N N IN G ON A LIN U X PLAT F ORM.
PICT U RES GEN ERAT ED U SIN G T HE Xfig PROGRAM.

106

You might also like