Analysis of discrete-time

systems
 Stability
• Stability of one solution (non-linear and/or time-varying
systems)
• System stability (global property of linear systems)
• Global stability vs. local stability (non-linear systems)
• (General) stability
• Asymptotic stability
• BIBO stability (Bounded Input - Bounded Output)

Stability of linear systems:

A linear, discrete, time-invariant system is
asymptotically stable, if and only if all the
eigenvalues of the system matrix Φ are inside the
unit circle.
 BIBO-stability vs. Lyapunov-stability

The general solution of the state equation x(t )  Ax(t )  Bu(t )

is

z
t

x(t )  e A ( t  t k ) x(t k )  e A ( t  s ') Bu( s' )ds'

tk

The behaviour of state x in the future depends on two terms:

autonomous part (initial conditions) and control inputs.

Lyapunov stability concerns the autonomous part. The initial state

is disturbed a bit, and it is investigated how this deviation behaves
in the future; no control inputs are used.
 BIBO-stability vs. Lyapunov-stability
BIBO-stability is related to the input/output-behaviour, and it is
connected to the second term of the solution. The system is BIBO-
stable (bounded input-bounded output), if a bounded input u leads
to a bounded output y.

E.g. a stock (integrator) and an ideal oscillator (harmonic oscillator)

are marginally stable (generally stable) and Lyapunov stable, but
not asymptotically stable nor BIBO-stable.

An inverted pendulum is unstable according to all definitions of

stability.

An ideally mixed vesssel (low-pass filter) is stable according to all

above definitions.
Marginal, asymptotic and general stability

System

Stable Unstable

Asymptotically stable Marginally stable

Often the definitions are simplified by stating that by stability
asymptotic stability is meant. (An asymptotically stable system is
always BIBO-stable.)
 Stability tests

• Direct calculation of the eigenvalues of Φ

• Study of the characteristic polynomial
• Root locus
• Nyquist criterion
• Lyapunov’s method
 The Jury stability criterion
Consider the characteristic polynomial
A( z )  a0 z n  a1z n 1    an
With Jury’s test it is easy to check, whether all the poles are
inside the unit circle i.e. whether the system is asymptotically
stable
a0 a1  an 1 an n  n
a
an an 1  a1 a0 a0
n 1
n 1 n 1 n 1 a
a0 a1  an 1  n 1  nn11
n 1 n 1 n 1 a0
an 1 an  2  a0
aik 1  aik   k akki

0 (a0n  a0 , a1n  a1 , , ann  an )
a0
 The Jury stability criterion
If a0 > 0, all the roots are inside the unit circle, if and only if

a0k  0, k  0,1,  , n  1
Consider a couple examples: The characteristic equation:

A( z )  3z 2  2 z  1 3 2 1 2  1
3
a0 a1 a2
1  ( 13 ) 2  ( 23 ) 3  (1)
a2 a1 a0
(3  13  ) 83 (2  23  ) 43 1  1
2
a01 a11 4
3  ( 23 ) 8
3  ( 43 )
a11 a01
( 83  23  )2
a00
 The Jury stability criterion

a0  0, a01  0, a00  0

so that all roots are inside the unit circle and the system is stable (in fact
the roots can easily be solved directly; – 0.3333 ± 0.4714i).

The Jury table is formed as follows: The last term of the first row is divided
by the first term. The second row is multiplied by this factor.
The second row is then subtracted from the first row, which gives the third
row. The fourth row is obtained by changing the row vector of the third row
upside down. Again a factor is formed by dividing the last term in the third
row by the first one. The procedure then continues as described above.
 The Jury stability criterion
The eigenvalues and roots of polynomials can most
conveniently be calculated by using numerical routines and
available software. The
Jury method on the other hand can easily be used in the case
of a small system, when a computer is not at hand.

The real power of the Jury criterion comes into action in

symbolic calculations. Stability can then be determined as a
function of one or even more parameters. Consider the
following example A( z )  z 2  a z  a
1 2

The characteristic polynomial:

 The Jury stability criterion
Form the Jury table
a0 a1 a2 1 a1 a2
 2  a2
a2 a1 a0 a2 a1 1
a1
1  ( a2 )2
a1 (1  a2 ) 1 
a01 a11 1  a2
a1 (1  a2 ) 1  ( a2 ) 2
a11 a01
( a ) 2
(1  a2 )
a00 1  ( a2 ) 2  1
1  a2

The system is stable,if R|1  (a )2

2
0
S|1  (a ) 2

( a1 ) 2
(1  a2 )
0
T 2
1  a2
 The Jury stability criterion
The Jury stability criterion
( a ) 2
(1  a ) ( a ) 2
(1  a2 )
1  ( a2 ) 
2 1 2
 (1  a2 )(1  a2 )  1
1  a2 1  a2
F
 (1  a )G (1  a ) 
(a ) I
J F
 (1  a )G
(1  a )
2

(a ) I
J
2 2

H 1 a K H 1 a 1 a K
1 2 1
2 2 2
2 2 2

F
 (1  a )G
(1  a )  ( a ) I 1  a
2

J  c (1  a )  ( a ) h
2

H 1 a K 1 a
2 1 2 2 2
2 2 1
2 2


b1  a gb1  a  a gb1  a  a g
2 2 1
0 2 1

1  a2
 The Jury stability criterion
The stability criteria thus become
R|b1  a gb1  a g  0
2 2

S| b1  a gb1  a  a gb1  a  a g  0
2 2 1 2 1

T 1 a 2

 S
R
| b1  a gb1  a g  0
2 2

|Tb1  a  a gb1  a  a g  0
2 1 2 1

R|1  a  1 R|1  a  1
2 2 The second group is not
 Sa   a  1  Sa   a  1 fulfilled with any values
|Ta  a  1
2

2 1
1
|Ta  a  1
2

2 1
1
of α1 and α2
 The Jury stability criterion
System with the characteristic polynomial A( z )  z 2  a1z  a2

is stable for the values α1 and α2 such that:

R|a
a2

2 1 1

S|a 2   a1  1 a1

Ta 2  a1  1 -2 -1 1 2

-1
 Stability in frequency domain
The frequency response of G(s) is G(iω), ω  [0, [ . It can
graphically be presented in the complex plane as the Nyquist
curve or as amplitude/phase curves as a function of frequency
(Bode diagram).

Correspondingly, for a discrete system H(z) the frequency

response is H(eiωh), ωh  [0, π] . This can also be presented
graphically as a discrete Nyquist or discrete Bode diagram.

The difference is, that in the discrete case only the frequency
interval ωh  [- π, π[ is considered.
 Stability in frequency domain
1
The continuous system G ( s)  2
. s 1
s  14
is sampled, (h = 0.4), giving the discrete ZOH-equivalent

0.066z  0.055
H ( z)  2
z  1450
. z  0.571
Compare the continuous frequency response G(iω) with the
discrete one H(eiωh), in the frequency range ω h  [0, π] .
 Stability in frequency domain
sys=tf(1,[1 1.4 1])
Bode Diagrams
Transfer function: N
1 0
--------------- continuous
s^2 + 1.4 s + 1 -20

Phase (deg); Magnitude (dB)

sysd=c2d(sys,0.4) -40
discrete
-60
Transfer function:
0.06609 z + 0.05481 0
---------------------
-50
z^2 - 1.45 z + 0.5712
-100 continuous

Sampling time: 0.4 -150

-200
discrete
w=logspace(-2,1,1000)'; -250
bode(sys,w) 10-1 100

hold Frequency  (rad/sec)

bode(sysd,w)
 Stability in frequency domain
Nyquist Diagrams

0.8

0.6

0.4
nyquist(sys,w) 0.2

Imaginary Axis
hold
0
nyquist(sysd,w)
-0.2

-0.4 continuous

-0.6

-0.8
discrete

-0.2 0 0.2 0.4 0.6 0.8 1

Real Axis
 Discrete Nyquist stability criterion
YREF (z) + Y(z)
Discrete control system H OL(z)
_
Characteristic equation 1  HOL ( z )  0

Stability can be determined by using the open loop HOL(z) Nyquist

diagram. HOL(eih) (open loop Nyquist curve) encircles the point -1
N times clockwise.
N ZP

in which Z is the number of zeros and P the number of poles of the

characteristic equation outside the unit circle.
 Discrete Nyquist stability criterion

This fact can be applied in stability analysis. The characteristic

equation (CE) has the form
numOL ( z ) denOL ( z )  numOL ( z ) numCE ( z )
1   0
denOL ( z ) denOL ( z ) denCE ( z )
The open loop (OL) poles are the same as the poles of the
characteristic equation. The zeros of the characteristic equation
determine stability so that if the characteristic equation has zeros
outside the unit circle, the closed loop system is unstable. The
stability criterion is thus obtained by setting Z= 0 and by demanding
that the Nyquist curve encircles point -1 P times counterclockwise.
(Z=N+P=0)
 Discrete Nyquist stability criterion
The criterion becomes simple, if the open loop pulse transfer
function has no poles outside the unit circle. Then the Nyquist
curve must not encircle the point –1 at all.
 Discrete Nyquist stability criterion

A process is controlled with a discrete P-controller, which has

the gain K (h=1)
YREF (z) + E(z) U(z) Y(z) 04.
K H(z) H(z) 
_ (z  05
. )(z  02
. )

The discrete Nyquist diagram is constructed with Matlab

» sysd=zpk([],[0.2 0.5],0.4,1);

» nyquist(sysd)
Discrete Nyquist stability criterion
Nyquist Diagrams
From: U(1)
1

0.8

0.6

0.4
Imaginary Axis

0.2
To: Y(1)

-0.2

-0.4

-0.6

-0.8

-1
-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1

Real Axis
 Discrete Nyquist stability criterion
The interception point in the x 10-4
Nyquist Diagrams
From: U(1)

real axis can be found e.g. by

5

using the zoom-command. By 3

inspection, the point is 2

Imaginary Axis
1

approximately

To: Y(1)
0

-0.4416 -1

-2

-3

The magnitude can thus be -4

-0.4421 -0.442 -0.4419 -0.4418 -0.4417 -0.4416 -0.4415 -0.4414 -0.4413 -0.4412

multiplied with (1/0.4416) to Real Axis

reach the critical point -1.

1
The controlled system is stable when K  0.4416  2.26
 Discrete Nyquist stability criterion
Stability can also be determined by direct calculus from the
pulse transfer function
0.4
H ( z) 
( z  0.5)( z  0.2)
Substitute z with eiωh = eiω = cos(ω)+i sin (ω) ,(Euler formula),
which gives the frequency response H (eiωh)
i 0.4
H ( e )  i
(e  0.5)( ei  0.2)
0.4

(cos   i sin   0.5)(cos   i sin   0.2)
 Discrete Nyquist stability criterion
0.4

(cos2   sin 2   0.7 cos   0.1)  ( 2 sin  cos   0.7 sin  )i
0.4

(2 cos2   0.7 cos   0.9)  ( 2 sin  cos   0.7 sin  )i

Setting the imaginary part 0 the interception point with the real
axis is obtained
b
2 sin  cos  0.7 sin   0  sin  0.7  2 cos  0 g
 (sin   0)  (0.7  2 cos  0)  (sin   0)  (cos  207 )
 (  0)  (  arccos 207 )
 Discrete Nyquist stability criterion

The frequency 0 describes the start point in the Nyquist curve

and the frequency arccos(7/20) the interception point with the
real axis. Substitute it to the frequency response function
0.4 4
) 7 2   0.444
7
i arccos 20
H (e
2( 20 )  0.7( 20 )  0.9
7
9
The interception point is -0.444. The gain of the controller K
can be multiplied by the factor (1/0.444) in order the crossing
at point –1 to take place. The controlled system is stable,when

K  94  2.25
Symbolic frequency response calculated with
the MATLAB symbolic toolbox

f='(4/10)/((z-(1/2))*(z-(2/10)))'
z='cos(w)+i*sin(w)'
g=subs(f,z)
g =

(4/10)/((cos(w)+i*sin(w)-(1/2))*(cos(w)+i*sin(w)-(2/10)))

g2=simplify(g)
g2 =

4/(20*cos(w)^2+20*i*cos(w)*sin(w)-7*cos(w)-9-7*i*sin(w))
 Gain and Phase margins

Definitions of gain and phase margins are identical to those of

the continuous time systems.

Open loop pulse transfer function H(z)

ωo is the lowest frequency, for which arg H (ei h )   o

ωc is the lowest frequency, for which H (ei h )  1 c

1
Gain marginAm arg  Phase margin  m arg    arg H (ei h )
c

H ( e i o h )
 Controllability, reachability, observability
and detectability
Principal questions:
* How can any state be transferred into any other state ?
* How can a state be determined from observations ?

RS
Consider the state-space x( k  1)   x( k )   u( k ), x(0)  x 0
realization T
y( k )  Cx( k )
The solution at the time instant n (n is the dimension of the
system or in other words the number of state components) is
x(n)   n x 0   n 1 u( 0)     u(n  1)   n x 0  Wc U
Wc      n 1
T
U  u ( n  1)
T
u (n  2)
T
 T
u ( 0)
 Controllability, reachability, observability
and detectability

If the rank of the matrix Wc is n, then n linear equations are

obtained, from which the controls U can be calculated, which
drive the system into any desired final state.

The system is controllable, if it is possible to find a control

sequence, which drives the system to the origin from any
state in a finite time interval.

The system is reachable, if it is possible to find a control

sequence, which drives the system from any state to any
state in a finite time interval.
 Controllability, reachability, observability
and detectability
Reachability is a stronger property than controllability. For
example, if Φn = 0, the state will go to the origin without any
controls, so that the system is controllable while not
necessarily reachable.
The system is reachable, if and only if the rank of Wc is n
(Wc is the controllability matrix)

Ex. System x( k  1)   x( k )   u( k ), x( 0)  2 1
T

is planned to be driven to i. origin x = [0 0]T

ii. state x = [0 3]T
iii. state x = [-1 2]T
 Controllability, reachability, observability
and detectability
L
M
0 0O
P L 1O
,  M P
N 1 0Q N 0Q
Consider the reachability

L L 1O L0 0O L1O O L L1O L0O O L1 0O

 MM P M P  MM P M P P  M
Wc   
0 1 0P M 0 P
NN Q N Q N Q Q NN Q N Q Q N Q
0 1 0 1P
For a square matrix the rank can be investigated by the
determinant det Wc  1  0 ,the rank is full,
so the system is reachable and thus controllable also
rank Wc  2  n
 Controllability, reachability, observability
and detectability
The original state can be driven to any other state with (n =
2) steps at the maximum.
x(0)  2 1
T

x(1) 
LM
0 0 OP LM OP LM OP LM OP
x(0) 
1
u(0) 
u(0)

u(0)
N
1 0 Q NQ N Q N Q
0 x1 (0) 2
L
x( 2)  M
0 0O L 1O L u(1) O L u(1) O
P0Q x(1)  MN0PQ u(1)  MN x (1)PQ  MNu(0)PQ
N1 1

i. x(2) = [0 0]T => U = [u(0) u(1)]T = [0 0]T

ii. x(2) = [0 3]T => U = [u(0) u(1)]T = [3 0]T
iii. x(2) = [-1 2]T => U = [u(0) u(1)]T = [2 -1]T
 Controllability, reachability, observability
and detectability

M
L 0 0O
P L
,  M P
0O
N
Consider reachability
1 0Q N 1Q

L L
W     M M P M
0O L0 0O L0O O L L0O L0O O L0 0O
P  MM P M P P  M
c
1 0P M 1 P
NN Q N Q N Q Q NN Q N Q Q N Q
1 1 0 1 0P
the rank is not full, so that the system is not reachable.
Nothing can be said about controllability by this analysis.
det Wc  0 rank Wc  n
 Controllability, reachability, observability
and detectability
Let us check how the state behaves x( 0)  2 1
T

L
x(1)  M
0 0 OP LM OP
x(0) 
0
u(0) 
0LM 
OP LM
0 OP
N1 0 Q NQ 1 N
x1 (0)  u(0) Q N
2  u(0) Q
L
x( 2)  M
0 0O L 0O L 0 O
P0Q x(1)  MN1PQ u(1)  MNu(1)PQ
N1
The first state cannot be influenced so that it is not
reachable. The origin can be reached and the system is
controllable. The last state (iii) cannot be reached by any
control sequence.
 Controllability, reachability, observability
and detectability
On the other hand, both the origin (i) and the second state (ii)
can be reached even with one control step

i. x(1) = [0 0]T => U = u(0) = -2

ii. x(1) = [0 3]T => U = u(0) = 1
iii. x(1) = [-1 2]T => Not reachable
If the aim is to reach the desired states only after two steps,
the first control step can be arbitrary

i. x(2) = [0 0]T => U = [u(0) u(1)]T = [* 0]T

ii. x(2) = [0 3]T => U = [u(0) u(1)]T = [* 3]T
iii. x(2) = [-1 2]T => Not reachable
 Controllability, reachability, observability
and detectability
RS
Consider x( k  1)   x( k )   u( k ), x(0)  x
Ty( k )  Cx(k ) 0

Let the output signal y and control signal u be known from

previous time instants. Based on this the aim is to find x0.
Consider the solution for y.
y(0)  Cx(0)  Cx0
b g
y(1)  Cx(1)  C x0   u(0)  Cx0  C u(0)
y(2)  Cx(2)  Cbx(1)   u(1)g  C x  C u(0)  C u(1)
2
0

 n 1
y(n  1)  C x0  C  n1i  u(i  1)
n 1

i 1
 Controllability, reachability, observability
and detectability
Because the control u is known at time instants k = 0 … n -
1, the determination of x0 does not depend on the weighted
sum of controls. The formula of y can be divided into two
parts; one depending on the initial condition yx and one
depending on controls yu.
n 1
y(n  1)  y x (n  1)  y u ( n  1)  C n 1x 0  C  n 1i  u(i  1)
i 1
The initial condition is found if yx can be solved at different
time instants.
Putting the equations at different time instants together gives
y x ( n  1)  C n 1x 0
 Controllability, reachability, observability
and detectability

LM y(0) OP LM C OP
Y 
x
M y(1)
P 
MM  PP MM  PPM C
P x 0  Wo x 0

Ny (n  1)Q NC Q
x
n 1

The initial condition x0 can be calculated,if Wo is of full rank.

The system is observable,if and only if the rank of Wo is n.

The system is detectable, if non-observable states are

stable.
 Canonical forms
The same system can be described by the difference
equation
y ( k  na )  a1 y ( k  na  1)  ana y ( k )  b0u( k  nb ) bnb u( k )

or by the pulse transfer operator or by the pulse transfer

function nb 1
nb 1 b z nb
 b z    bnb
b0q  b1q
nb
   bnb H ( z )  0 1
H (q )  na na 1 z na
 a z na 1
   ana
q  a1q    ana 1

or by the state-space representation. The last alternative is

not unique, since it is possible to form an indefinite number of
state-space representations, which give the same input-
output behaviour (e.g. diagonal form or Jordan form)
 Canonical forms
With respect to reachability and observability the most
important forms are the controllable canonical form and the
observable canonical form.
Controllable canonical form:
LM
 a1  a2   ana 1  ana OP LM OP
1

M
x( k  1)  M 0
1 0 

0 0
PP x(k )  MM00PP u(k )
MM  1
 
0

0

PP MM  PP
MN 0 0  1 0 PQ MN0PQ
y( k )  bnb  na 1 bnb  na  2  bnb 1 bnb x( k )
 Canonical forms
The observable canonical form:

LM a 1 1 0  0 OP LM
bnb  na 1 OP
MM a2 0 1  0
PP MM
bnb  na  2
PP
x( k  1)      x( k )  
MMana 1 0 0  1
PP MM
bnb 1
PP
u( k )

MN a
na 0 0  0 PQ bnb MN PQ
y( k )  1 0 0  0 x( k )
Both of these have an alternative version, in which the states
have been chosen in the reverse order.
 Canonical forms
The above forms are as such valid only for strictly proper
systems, but by small modifications they can be modified
also in the case that the D-matrix is non-zero. However, the
state-space representation must always be causal, i.e. na 
nb.

If nb is much smaller than na ,the formulas will point to

coefficients of the B-polynomial, which do not exist (e.g. b-1 ,
b-2 , …). These can always be set to zero.

0 0
  b2 u( k  nb  2)  b1 u( k  nb  1)  b0u( k  nb ) bnb u( k )
0 0
 Canonical forms

Develop controllable canonical forms for the given pulse

transfer functions:

nb  1 1 nb  0
2z  1 Hb ( z)  2 ,
Ha ( z)  2 , z  2z  1 na  2
z  2z  1 na  2
a1  2, a2  1, b0  2, b1  1 a1  2, a2  1, b0  1
L
x( k  1)  M
2 1O L1O L
x( k  1)  M
2 1O
P L1O
x( k )  M Pu( k )
Px( k )  M Pu( k )
N 1 0 Q N0Q N 1 0 Q N0Q
y( k )  2 1 x( k ) y( k )  0 1 x( k )
Non-reachable and/or non-observable systems
There may be several reasons, why a discrete system is not
reachable or observable:

* The original continuous system (which is then sampled) is

not reachable or observable.

* Hidden oscillations (sampling frequency too low)

* Pole-zero cancellation. Reachability is lost, if sampling

leads to a system with a common pole and zero. The
sampling interval must be changed.
 Analysis of simple control loops

Control problems:

* Regulator problem
Setpoint is constant

* Combination of regulator and servo problems

e.g. several but rare step changes in the setpoint.

* Servo problem
The changing setpoint trajectory must be followed.
 Analysis of simple control loops
Classification of disturbances:

* Load disturbance
Influence on control variables, often stepwise and
change the long term average (low frequencies)

* Measurement noise
Often high-frequency noise caused by measurement
devices.

* Parameter changes
System parameters change with time.
 Analysis of simple control loops
Typical disturbance models, which are used in system
analysis

* Impulse and pulse

* Step

* Ramp

* Sinusoid

* Noise
 Hidden oscillations
The following model is 2

simulated and the response 1

is sampled with different 0

0 5 10 15 20

sampling intervals (h = 3, 2, 2

1, 0.5) 1

0
0 5 10 15 20
2

0
0 5 10 15 20
2

0
0 5 10 15 20
 Hidden oscillations
2

1
What is really happening in
0
the system, can be seen from 2
0 5 10 15 20

the figures. 1

0
This is an example of hidden 2
0 5 10 15 20

oscillations or ripple. 1

0
0 5 10 15 20
2

0
0 5 10 15 20