Identification, Analysis and

Control of Dynamical Systems

Alberto Bemporad

[Link]/~bemporad/intro_control_course.html

Academic year 2020-2021

Course contents

1. Linear dynamical systems in continuous and discrete-time

2. Linearization and discretization, stability analysis

3. Controllability and observability analysis

4. Synthesis of feedback controllers and state estimators

5. System identification (=learn dynamical models from data)

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 2

Dynamical systems
Dynamical systems
• A dynamical system is an object (or a set of objects) that evolves over time,
possibly under external excitations.

• Examples: an engine, a satellite, a tank reactor, a human transporter, ...

throttle engine
position torque

thrusters attitude

reactant
concentration/
temperature
motor vertical
vessel voltage position
concentration/
temperature

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 3

Dynamical systems
• ... a supply chain, a portfolio, a computer server

asset portfolio
quantities wealth task quality
allocation of service

raw items
materials sold

factory warehouse distributor/retailer

• The way the system evolves over time is called the dynamics of the system.

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 4

Dynamical models

• A dynamical model of a system is a set of mathematical laws that explain how

the system evolves over time, usually under the effect of external excitations, in
a quantitative way.

• What is the purpose of a dynamical model ?

1. Understand the system (“How does X influence Y ?”)

2. Simulation (“What happens if I apply action Z on the system ?”)

3. Estimate (“How to estimate variable X from measuring Y ?”)

4. Control (“How to make the system behave autonomously the way I want ?”)

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 5

Linear systems
Continuous-time linear systems
• System of n first-order linear ordinary differential equations (ODEs) with inputs



ẋ1 (t) = a11 x1 (t) + . . . + a1n xn (t) +b1 u(t)



 ẋ2 (t) = a21 x1 (t) + . . . + a2n xn (t) +b2 u(t) [ ]
 dx
. . .
. . . ẋ =

 . . . dt



 ẋn (t) = an1 x1 (t) + . . . + ann xn (t) +bn u(t)

x1 (0) = x10 , ... xn (0) = xn0

• Set x = [x1 . . . xn ]′ ∈ Rn . The equivalent matrix form of the linear ODE

system is the so-called linear system

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

with initial condition x(0) = x0 , with vector x0 = [x10 . . . xn0 ]′ ∈ Rn

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 6

Example: Mass-spring-damper system
x!(t), x"(t)
K
u(t)
! M

(
ẋ1 (t) = x2 (t) velocity = derivative of traveled space

M ẋ2 (t) = u − βx2 (t) − Kx1 (t) Newton's law

Rewrite as the 2nd order linear system

( dx1 (t)
dt = x2 (t)
dx2 (t)
dt = −M
β
x2 (t) − K
M x1 (t) + 1
M u(t)
or in matrix form
" # " #
0 1 0
ẋ(t) = x(t) + u(t)
−M
K
−M
β 1
M
| {z } | {z }
A B
©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 7
nth -order linear ODE with input
dy (n) (t) dy (n−1) (t)
n
+ an−1 + · · · + a1 ẏ(t) + a0 y(t)
dt dtn−1
du (n−1) (t) du(n−2) (t)
= bn−1 + bn−2 + · · · + b1 u̇(t) + b0 u(t)
dt dt

By inspection, the nth -order ODE = 1st -order linear system of ODEs

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




ẋ1 (t) = x2 (t) y(t) = Cx(t) + Du(t)



 ẋ2 (t) = x3 (t)
  
. . 0 1 0 ... 0 0
. . 0 0 1 ... 0

 . .  . 
0


 −a0 x1 (t) + . . . − an−1 xn (t) + u(t) A= . .. . ,B = 
 .
.
ẋ
 n
(t) =  .. .
. . .
.  .
y(t) = b0 x1 (t) + . . . + bn−1 xn (t) 0 0 0 ... 1 0
−a0 −a1 −a2 ... −an−1 1
C= [ b0 b1 b2 ... bn−1 ], D = 0
st
The linear system of 1 -order ODEs is called the state-space realization of the
nth -order ODE. There are infinitely many state-space realizations.

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 8

Lagrange's formula
• Starting from the initial condition x(0) = x0 , the continuous-time linear system
ẋ = Ax + Bu has the unique solution x(t)

Z t
At
x(t) = e x0 + eA(t−τ ) Bu(τ )dτ
| {z }
natural response | {z }
0

forced response

• The exponential matrix is defined as

MATLAB
A2 t 2 An t n E=expm(A*t)
e At
≜ I + At + + ... + + ...
2 n!
(Moler, Van Loan, 2003) Python
from [Link] import expm
E=expm(A*t)
©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 9
State vector

• Given x(0) and u(t), ∀t ∈ [0, T ], Lagrange’s formula allows us to compute x(t)
and y(t), ∀t ∈ [0, T ]

• Generally speaking, the state of a dynamical system is a set of variables that

completely summarizes the past history of the system. It allows us to predict its
future motion

• Therefore, by knowing the initial state x(0) we can neglect all past history
u(−t), x(−t), ∀t ≥ 0

• The dimension n of the state x(t) ∈ Rn is called the order of the system

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 10

Eigenvalues and eigenvectors
• Let us recall some basic concepts of linear algebra:
 
a11 a12 ... a1n
 a21 a22 ... a2n 
A=  .  square matrix of order n, A ∈ Rn×n
 .. .
.
. ...
.
.
.

an1 an2 ... ann
 
1 0 ... 0
0 1 ... 0
I= . .  identity matrix of order n
 .. .
.
. ... . 
.
0 0 ... 1

• Characteristic equation of A:

det(λI − A) = 0

• Characteristic polynomial of A:

P (λ) = det(λI − A) = λn + an−1 λn−1 + . . . + a1 λ + a0

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 11

Eigenvalues and eigenvectors
• The eigenvalues of A ∈ Rn×n are the roots λ1 , …, λn of its characteristic
polynomial
det(λi I − A) = 0, i = 1, 2, . . . , n
• An eigenvector of A is any vector vi ∈ Rn such that Avi = λi vi for some
i = 1, 2, . . . , n.
• The diagonalization of A is A = T ΛT −1 , where
 λ1 0 ... 0

0 λ2 ... 0
Λ =  .. .. . . ..  = T −1 AT, T = [v1 |v2 | . . . |vn ]
. . . .
0 0 ... λn

(not all matrices A are diagonalizable, see Jordan normal form)

• Algebraic multiplicity of λi = number of coincident roots λi of det(λI − A)

• Geometric multiplicity of λi = number of linearly independent eigenvectors vi

such that Avi = λi vi .

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 12

Eigenvalues and modes
• Let u(t) ≡ 0 and assume A diagonalizable

• The state trajectory is the natural response

" #
eλ1 t ... 0
At Λt −1 ..
x(t) = e x(0) = T e T x0 = [v1 . . . vn ] α
| {z } .
λn t
α 0 ... e
" #
h i α1 Xn
= v 1 e λ1 t ... v n eλn t .. = α i e λi t v i
.
αn i=1

where vi =eigenvector of A, λi =eigenvalue of A, α = T −1 x(0) ∈ Rn

• The evolution of the system depends on the eigenvalues λi of A, called modes

of the system (sometimes we also refer to eλi t as the i-th mode)

• A mode λi is called excited if αi 6= 0

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 13

Some classes of dynamical systems
• Causality: a dynamical system is causal if y(t) does not depend on future inputs
u(τ ) ∀τ > t (strictly causal if ∀τ ≥ t)

• A linear system is always causal, and strictly causal iff D = 0

• Linear time-varying (LTV) systems:

(
ẋ(t) = A(t)x(t) + B(t)u(t)
y(t) = C(t)x(t) + D(t)u(t)

• When A, B , C , D are constant, the system is said linear time-invariant (LTI)

• A generalization of LTV systems are linear parameter-varying (LPV) systems

(
ẋ(t) = A(p(t))x(t) + B(p(t))u(t)
y(t) = C(p(t))x(t) + D(p(t))u(t)

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 14

Some classes of dynamical systems
• Multivariable systems: more generally, a system can have m inputs
(u(t) ∈ Rm ) and p outputs (y(t) ∈ Rp ). For linear systems, we still have
(
ẋ(t) = Ax(t) + Bu(t)
y(t) = Cx(t) + Du(t)

with A ∈ Rn×n , B ∈ Rn×m , C ∈ Rp×n , D ∈ Rp×m

• Nonlinear systems (
ẋ(t) = f (x(t), u(t))
y(t) = g(x(t), u(t))
where f : Rn+m → Rn , g : Rn+m → Rp are (arbitrary) nonlinear functions

• Time-varying nonlinear systems

(
ẋ(t) = f (t, x(t), u(t))
y(t) = g(t, x(t), u(t))

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 15

Stability
Equilibrium
• Consider the continuous-time nonlinear system
(
ẋ(t) = f (x(t), u(t))
y(t) = g(x(t), u(t))

Definition
A state xr ∈ Rn and an input ur ∈ Rm are an equilibrium pair if for initial
condition x(0) = xr and constant input u(t) ≡ ur the state remains constant:
x(t) ≡ xr , ∀t ≥ 0.

• Equivalent definition: (xr , ur ) is an equilibrium pair if f (xr , ur ) = 0

• xr is called equilibrium state, ur equilibrium input
• The definition generalizes to time-varying nonlinear systems

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 16

Stability
• Consider the nonlinear system
(
ẋ(t) = f (x(t), ur )
y(t) = g(x(t), ur )

and let xr be an equilibrium state, f (xr , ur ) = 0

Definition
The equilibrium state xr is stable if for each initial conditions x(0) “close
enough” to xr , the corresponding trajectory x(t) remains near xr for all t ≥ 0.

Math definition: ∀ϵ > 0 ∃δ > 0 : ∀x(0) such that ∥x(0) − xr ∥ < δ ⇒ ∥x(t) − xr ∥ < ϵ, ∀t ≥ 0.

• The equilibrium point xr is called asymptotically stable if it is stable and

x(t) → xr for t → ∞

• Otherwise, the equilibrium point xr is called unstable

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 17

Stability of equilibria - Examples
1
1
1 0.8
0.8
0.8 0.6
0.6
0.6 0.4
0.4
0.4 0.2

x2(t)
0.2
0

x2(t)
0.2
x (t)

0
0 −0.2
2

−0.2
−0.2 −0.4
−0.4
−0.4 −0.6
−0.6
−0.6 −0.8
−0.8
−0.8 −1
−1
−1 −1 −0.8 −0.6 −0.4 −0.2 0 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 x1(t)
−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 x1(t)
x1(t)

unstable
stable equilibrium asymptotically
equilibrium
stable equilibrium
1
0.6
x1 2
0.5 0.4 x1
x2 0.2 0
0
x2
0 −2
−0.2
−4
−0.5
−0.4
−6 x1
−0.6
−1
0 0.5 1 1.5 2 2.5 3 3.5 4
−0.8
−8 x2
0 0.5 1 1.5 2 2.5 3 3.5 4
−10
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

   
dx −2x1 (t) − 4x2 (t) dx −x1 (t) − 2x2 (t)  
= = dx 2x1 (t) − 2x2 (t)
dt 2x1 (t) + 2x2 (t) dt 2x1 (t) − x2 (t) =
dt x1 (t)

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 18

Stability of first-order linear systems
• Consider the first-order linear system
ẋ(t) = ax(t) + bu(t)
• xr = 0, ur = 0 is an equilibrium pair
• For u(t) ≡ 0, ∀t ≥ 0, the solution is
x(t) = eat x0
• The origin xr = 0 is

– unstable if a > 0 4.5

– stable if a ≤ 0
4

3.5

– asymptotically stable if a < 0 3

a>0
x(t)
2.5

1.5
x0
1
a=0
0.5
xr a<0
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
t

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 19

Stability of continuous-time linear systems
Since the natural response of ẋ = Ax + Bu is x(t) = eAt x0 , the stability
properties depend only on A. We can therefore talk about system stability of a
linear system (A, B, C, D)
Theorem
Let λ1 , . . ., λm , m ≤ n be the eigenvalues of A ∈ Rn×n .
The system ẋ = Ax + Bu is
• asymptotically stable iff <λi < 0, ∀i = 1, . . . , m
• (marginally) stable if <λi ≤ 0, ∀i = 1, . . . , m, and the eigenvalues
with null real part have equal algebraic and geometric multiplicity
• unstable otherwise (in particular, if ∃ i such that <λi > 0).

The stability properties of a linear system only depend on the real part of
the eigenvalues of matrix A

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 20

Stability of continuous-time linear systems
Proof:
2 2 n n
• The natural response is x(t) = eAt x0 (eAt ≜ I + At + A 2t + . . . + An!t + . . . )
• If matrix A is diagonalizable1 , A = T ΛT −1 ,
λ   eλ 1 t 
1 0 ... 0 0 ... 0

0 λ2 ... 0
  0 eλ2 t ... 0  −1
. ⇒ e =T
At
Λ= . . . . . .. . T
. . .. . . . . .
. . . . . .
0 0 ... λn 0 0 ... eλn t

• Take any eigenvalue λ = a + jb:

|eλt | = eat |ejbt | = eat

• A is always diagonalizable if algebraic multiplicity = geometric multiplicity

□
1 If A is not diagonalizable, it can be transformed to Jordan form. In this case the natural response
x(t) contains modes tj eλt , j = 0, 1, . . . , alg. multiplicity - geom. multiplicity

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 21

Linearization of nonlinear systems
• Consider the nonlinear system
(
ẋ(t) = f (x(t), u(t))
y(t) = g(x(t), u(t))

• Let (xr , ur ) be an equilibrium, f (xr , ur ) = 0

• Objective: investigate the dynamic behaviour of the system for small
perturbations ∆u(t) ≜ u(t) − ur and ∆x(0) ≜ x(0) − xr .
• The evolution of ∆x(t) ≜ x(t) − xr is given by

˙
∆x(t) = ẋ(t) − ẋr = f (x(t), u(t))
= f (∆x(t) + xr , ∆u(t) + ur )
∂f ∂f
≈ (xr , ur ) ∆x(t) + (xr , ur ) ∆u(t)
∂x
| {z } |∂u {z }
A B

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 22

Linearization of nonlinear systems

• Similarly
∂g ∂g
∆y(t) ≈ (xr , ur ) ∆x(t) + (xr , ur ) ∆u(t)
|∂x {z } |∂u {z }
C D

where ∆y(t) ≜ y(t) − g(xr , ur ) is the perturbation of the output from its
equilibrium

• The perturbations ∆x(t), ∆y(t), and ∆u(t) are (approximately) ruled by the
linearized system
(
˙
∆x(t) = A∆x(t) + B∆u(t)
∆y(t) = C∆x(t) + D∆u(t)

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 23

Lyapunov's stability
Lyapunov's indirect method
• Consider the nonlinear system ẋ = f (x), with f differentiable, and assume
x = 0 is equilibrium point (f (0) = 0)
∂f
• Consider the linearized system ẋ = Ax, with A = ∂x
x=0

1. If x = 0 is an asymptotically stable equilibrium for ẋ = Ax ⇒ it is (locally)

asymptotically stable for the nonlinear system

2. If x = 0 is an unstable equilibrium for ẋ = Ax ⇒ it is unstable for the nonlinear

system

3. If x = 0 is marginally stable for ẋ = Ax ⇒ nothing can be said about its stability

for the nonlinear system

Aleksandr Mikhailovich Lyapunov

(1857-1918)

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 24

Example: Pendulum
y(t) = angular displacement
ẏ(t) = angular velocity
ÿ(t) = angular acceleration
l
u(t) = mg gravity force
y(t)
h hẏ(t) = viscous friction torque
m
l = pendulum length
u(t) = mg
ml2 = pendulum rotational inertia
• mathematical model

ml2 ÿ(t) = −lmg sin y(t) − hẏ(t)

• in state-space form (x1 = y , x2 = ẏ )

(
ẋ1 = x2
ẋ2 = − gl sin x1 − Hx2 , H≜ h
ml2

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 25

Example: Pendulum

Look for equilibrium states:

" # " # (
x2r 0 x2r = 0
= ⇒
− l sin x1r − Hx2r
g
0 x1r = ±kπ, k = 0, 1, . . .

u(t) = mg
h
m

h
m

u(t) = mg

x2r = 0, x1r = 0, ±2π, . . . x2r = 0, x1r = ±π, ±3π, . . .

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 26

Example: Pendulum
• Linearize the system around x1r = 0, x2r = 0
" #
0 1
∆ẋ(t) = ∆x(t)
− gl −H
| {z }
A

• find the eigenvalues of A

 r 
g 1 g
det(λI − A) = λ2 + Hλ + = 0 ⇒ λ1,2 = −H ± H2 − 4
l 2 l
1.5

1
• <λ1,2 < 0 ⇒ ẋ = Ax asymptotically
0.5
stable

y(t)
0

• by Lyapunov’s indirect method −0.5

xr = [ 00 ] is also an asymptotically −1

stable equilibrium for the pendulum −1.5

0 2 4 6 8 10
t

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 27

Example: Pendulum
• Linearize the system around x1r = π , x2r = 0
" #
0 1
∆ẋ(t) = ∆x(t)
g
l −H
| {z }
A

• find the eigenvalues of A

 r 
g 1 g
det(λI − A) = λ + Hλ − = 0 ⇒ λ1,2 =
2
−H ± H2 +4
l 2 l
4

2
• λ1 < 0, λ2 > 0 ⇒ ẋ = Ax unstable 1

y(t)
0
• by Lyapunov’s indirect method −1

xr = [ π0 ] is also an unstable −2

equilibrium for the pendulum −3

−4
0 2 4 6 8 10
t

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 28

Lyapunov's direct method
• A second method exists to analyze global stability of nonlinear systems, based
on the concept of Lyapunov functions
• Key idea: if the energy of a system dissipates over time, the system
asymptotically reaches a minimum-energy configuration
• Assumptions: consider the autonomous nonlinear system ẋ = f (x), with f (·)
differentiable, and let x = 0 be an equilibrium (f (0) = 0)
• Some definitions of positive definiteness of a function V : Rn 7→ R
– V is locally positive definite if V (0) = 0 and there exists a ball
Bϵ = {x : kxk2 ≤ ϵ} around the origin such that V (x) > 0 ∀x ∈ Bϵ \ 0
– V is globally positive definite if Bϵ = Rn (i.e. ϵ → ∞)
– V is negative definite if −V is positive definite
– V is positive semi-definite if V (x) ≥ 0 ∀x ∈ Bϵ
– V is negative semi-definite if −V is positive semi-definite

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 29

Lyapunov's direct method

• Example: let x = [x1 x2 ]′ , V : R2 → R

– V (x) = x21 + x22 is globally positive definite

– V (x) = x21 + x22 − x31 is locally positive definite

– V (x) = x41 + sin2 (x2 ) is locally positive definite and globally positive semi-definite

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 30

Lyapunov's direct method
Theorem
Given the nonlinear system ẋ = f (x), f (0) = 0, let V : Rn 7→ R be positive definite
in a ball Bϵ around the origin, ϵ > 0, V ∈ C 1 (R). If the function

V̇ (x) = ∇V (x)′ ẋ = ∇V (x)′ f (x)

is negative definite on Bϵ , then the origin is an asymptotically stable equilibrium point.

If V̇ (x) is only negative semi-definite on Bϵ , then the the origin is a stable equilibrium
point.

Such a function V : Rn 7→ R is called a Lyapunov function for the system

ẋ = f (x)
©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 31
Example of Lyapunov's direct method
• Consider the following nonlinear system ẋ = f (x) (Slotine, Li, 1991 - Example 3.8)
(
ẋ1 = x1 (x21 + x22 − 2) − 4x1 x22
ẋ2 = 4x21 x2 + x2 (x21 + x22 − 2)

• The state x = 0 is an equilibrium because ẋ = f (0) = 0

• Consider the candidate Lyapunov function

V (x1 , x2 ) = x21 + x22

which is globally positive definite. Its time derivative V̇ is

V̇ (x1 , x2 ) = 2(x21 + x22 )(x21 + x22 − 2)

• It is easy to check that V̇ (x1 , x2 ) is negative definite if kxk22 = x21 + x22 < 2
√
• Since for any Bϵ with 0 < ϵ < 2 the hypotheses of Lyapunov’s theorem are
satisfied, x = 0 is an asymptotically stable equilibrium
√
• Any Bϵ with 0 < ϵ < 2 is a also a domain of attraction

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 32

Example of Lyapunov's direct method (cont'd)
• Cf. Lyapunov’s indirect method: the linearization around x = 0 is
" # " #
∂f (0, 0) 3x21 − 3x22 − 2 −6x1 x2 −2 0
= =
∂x 10x1 x2 5x1 + 3x22 − 2
2
0 −2
x=0

which is an asymptotically stable matrix

• Lyapunov’s indirect method tells us that the origin is locally asymptotically
stable
• Lyapunov’s direct method also tells us that Bϵ is a domain of attraction for all
√
0<ϵ< 2

• Consider this other example: ẋ = −x3 . The origin as an equilibrium. But

df (0)
dx = −3 · 0 = 0, so Lyapunov indirect method is useless.
2

• Lyapunov’s direct method with V = x2 provides V̇ = −2x4 , and therefore we

can conclude that x = 0 is (globally) asymptotically stable

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 33

Case of continuous-time linear systems
• Let us apply Lyapunov’s direct method to linear systems ẋ = Ax and choose
V (x) = x′ P x, with P = P ′  0 (P =positive definite and symmetric matrix)
• The derivative V̇ (x) = ẋ′ P x + x′ P ẋ = x′ (A′ P + P A)x
• V̇ (x) is negative definite if and only if the Lyapunov equation

A′ P + P A = −Q

is satisfied for some Q  0 (for example, Q = I )

Theorem
The autonomous linear system ẋ = Ax is asymptotically stable ⇔ ∀Q  0 the
Lyapunov equation A′ P + P A = −Q has one and only one solution P  0

Python
MATLAB
import control as ctrl
P=lyap(A’,Q)
P=lyap([Link](),Q)
(note transposition of matrix A !)

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 34

Discrete-time systems
Discrete-time models
y(t), y(kTs) u(kTs)
3.5 4

3.5
3
3

2.5 2.5

2
2
1.5

1.5 1
0 1 2 3 4 5 0 1 2 3 4 5
time t time t

sampled continuous-time signal discrete-time signal

• Discrete-time models describe relationships between sampled variables

x(kTs ), u(kTs ), y(kTs ), k = 0, 1, . . .
• The value u(kTs ) is kept constant during the sampling interval [kTs , (k + 1)Ts )
• A discrete-time signal can either represent the sampling of a continuous-time
signal, or be an intrinsically discrete signal
• Discrete-time signals are at the basis of digital controllers (as well as of digital
filters in signal processing)

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 35

Difference equation
• Consider the first-order difference equation (autonomous system)
(
x(k + 1) = ax(k)
x(0) = x0

• The solution is x(k) = ak x0

3

2.5
a>1
2
x(k)

1.5

a=1
1

0.5

0<a<1
0
0 1 2 3 4 5 6 7 8 9 10
k

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 36

Linear discrete-time system
• Consider the set of n first-order linear difference equations forced by the input
u(k) ∈ R


 x1 (k + 1) = a11 x1 (k) + . . . + a1n xn (k) +b1 u(k)



 x 2 (k + 1) = a21 x1 (k) + . . . + a2n xn (k) +b2 u(k)

.. .. ..
 . . .



 xn (k + 1) = an1 x1 (k) + . . . + ann xn (k) +bn u(k)


x1 (0) = x10 , . . . xn (0) = xn0

• In compact matrix form:

(
x(k + 1) = Ax(k) + Bu(k)
x(0) = x0
" x1 #
where x = ... ∈ Rn .
xn

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 37

Linear discrete-time system

• The solution is
X
k−1
x(k) = Ak x 0 + Ai Bu(k − 1 − i)
| {z }
i=0
natural response | {z }
forced response

• If matrix A is diagonalizable, A = T ΛT −1
 λ1   λk 0 ... 0

0 ... 0 1
0 λ2 ... 0   0 λk2 ... 0
Λ =  .. .. . . ..  ⇒ Ak = T  . . . .  T −1
. . . . . . .
. . . ..
0 0 ... λn 0 0 ... λk
n

where T = [v1 . . . vn ] collects n independent eigenvectors.

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 38

Modal response
• Assume input u(k) = 0, ∀k ≥ 0
• Assume A is diagonalizable, A = T ΛT −1
• The state trajectory (natural response) is

X
n
x(k) = Ak x0 = T Λk T −1 x0 = αi λki vi
i=1

where
– λi = eigenvalues of A
– vi = eigenvectors of A
– αi = coefficients that depend on the initial condition x(0)
" α1 #
α= .. = T −1 x(0), T = [v1 . . . vn ]
.
αn

• The system modes depend on the eigenvalues of A, as in continuous-time

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 39

Example - Wealth of a bank account

• k = year counter
• ρ= interest rate
• x(k)= wealth at the beginning of year k
• u(k)= money saved at the end of year k
• x0 = initial wealth in bank account

(
Discrete-time model: x(k + 1) = (1 + ρ)x(k) + u(k)
x(0) = x0
Stored amount of money (keur)
x0 10 ke 50

45

u(k) 5 ke 40

35

ρ 10 % 30

x(k)
25

20

15

1 − (1.1) k 10

x(k) = (1.1) ·10+ k

5 = 60(1.1)k −50 5

1 − 1.1 0
0 0.5 1 1.5 2 2.5
k (years)
3 3.5 4 4.5 5

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 40

Example - Supply chain
!1x1(k) !2x2(k)
"2x2(k)
u(k) "1x1(k) y(k)
#3x3(k) $3x3(k)
S P R
x1(k) x2(k) x3(k)

• Problem statement:
– At each month k , S purchases the quantity u(k) of raw material
– A fraction δ1 of raw material is discarded, a fraction α1 is shipped to producer P
– A fraction α2 of product is sold by P to retailer R, a fraction δ2 is discarded
– Retailer R returns a fraction β3 of defective products every month
and sells a fraction γ3 to customers
• Mathematical model:
 k month counter

 x1 (k + 1) = (1 − α1 − δ1 )x1 (k) + u(k)

 x1 (k) raw material in S

 x2 (k + 1) = α1 x1 (k) + (1 − α2 − δ2 )x2 (k) x2 (k) products in P
+β3 x3 (k) x3 (k) products in R



 x3 (k + 1) = α2 x2 (k) + (1 − β3 − γ3 )x3 (k) u(k) raw material purchased by S


y(k) = γ3 x3 (k) y(k) products sold to customers

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 41

Example - Student population dynamics
• Problem statement:
– 3-years course
– percentage of promoted, repeaters, and dropouts are roughly constant
– direct enrollment in 2nd and 3rd academic year is not allowed
– students cannot enroll for more than 3 years

k Year
xi (k) Number of students enrolled in year i at year k, i = 1, 2, 3
u(k) Number of freshmen at year k
y(k) Number of graduates at year k
αi promotion rate during year i, 0 ≤ αi ≤ 1
βi failure rate during year i, 0 ≤ βi ≤ 1
γi dropout rate during year i, γi = 1 − αi − βi ≥ 0

• 3rd -order linear discrete-time system:



 x1 (k + 1) = x1 (k) − α1 x1 (k) − γ1 x1 (k) + u(k) = β1 x1 (k) + u(k)

 x (k + 1)
2 = α1 x1 (k) + β2 x2 (k)
 x3 (k + 1)
 = α2 x2 (k) + β3 x3 (k)


y(k) = α3 x3 (k)

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 42

Example - Student population dynamics
• In matrix form
    

 β1 0 0 1

    

  α1 0  x(k) +  0  u(k)
 x(k + 1) = β2
0 α2 β3 0



 h i


 y(k) = 0 0 α3 x(k)

• Simulation
80

α1 = .60 β1 = .20 60

40

α2 = .80 β2 = .15 20

α3 = .90 β3 = .08 0
2020 2022 2024 2026 2028 2030 2032 2034 2036 2038 2040

110

u(k) ≡ 100, k = 2020, . . . 105

100

lim y(k) ≈ 69.0537 95

k→∞ 90
2020 2022 2024 2026 2028 2030 2032 2034 2036 2038 2040

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 43

nth -order difference equation
• Consider the nth -order difference equation forced by u
an y(k − n) + an−1 y(k − n + 1) + · · · + a1 y(k − 1) + y(k)
= bn u(k − n) + · · · + b1 u(k − 1) + b0 u(k)

• Equivalent linear discrete-time system in canonical state matrix form

    


0 1 0 ... 0 0

   0

 0  0 1 ... 0   

 

 x(k + 1) = .  . .   . 
.  . .   x(k) +  .  u(k)
 . 
.  . .   

 0  0 0 ... 1  0



 −a −a −a −a

 [ n n−1 n−2 . . . ] 1 1


y(k) = (bn − b0 an ) . . . (b1 − b0 a1 ) x(k) + b0 u(k)

MATLAB Python
• There are infinitely many state-space realizations
tf2ss ctrl.tf2ss

• nth -order difference equations are very useful for digital filters, digital
controllers, and to reconstruct models from data (system identification)

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 44

Some state-space realization methods
• The following state-space realization is called controllable canonical form)
 0 1 0 ... 0
 0 MATLAB
0 0 1 ... 0 0 sysc=canon(ss(A,B,C,D),’companion’)
 . 
A= . .. . , B = 
 .
.
 .. .
. . .
.  . Python
0 0 0 ... 1 0
−an −an−1 −an−2 ... −a1 1 sysc,T=ctrl.canonical_form(
C = [ bn bn−1 ... b1 ], D = 0 (b0 = 0) [Link](A,B,C,D), form=’reachable’)

• The following state-space realization is called observable canonical form

 −a1 1 0 0 ... 0   b1  MATLAB
−a2 0 1 ... 0 0 b2 sys=canon(ss(A’,C’,B’,D),’companion’)
 . .  
A=

.
.
. . .
. .. . . . ,B= .
.  syso=ss(sys.A’,sys.C’,sys.B’,D)
. . . . .  . 
−an−1 0 0 ... 0 1 bn−1
−an 0 0 ... 0 0 bn
Python
syso,T=ctrl.canonical_form(
C = [ 1 0 0 ... 0 0], D=0 (b0 = 0) [Link](A,B,C,D), form=’observable’)

• We will see later that (A, B) in controllable canonical form is a reachable pair,
(A, C) in observable canonical form is an observable pair

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 45

Mapping past I/O pairs to state vector
• The observable canonical form of the nth -order difference equation
X
n X
n
y(k) = − ai y(k − i) + bi u(k − i) (b0 = 0)
i=1 i=1

corresponds to the following definition of the state vector x(k):



 x1 (k) = y(k)

 P Pn

 x2 (k) = − n i=2 ai y(k + 1 − i) + i=2 bi u(k + 1 − i)



 .. ..

 . .
 Pn Pn
xj (k) = − ai y(k + j − 1 − i) + bi u(k + j − 1 − i)


i=j i=j

 .. ..

 . .



 −an−1 y(k − 1) − an y(k − 2) + bn−1 u(k − 1) + bn u(k − 2)

 x n−1 (k) =

xn (k) = −an y(k − 1) + bn u(k − 1)

• This is easy to verify by inspection, just compute x(k + 1) and check

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 46

Discrete-time linear system


 x(k + 1) = Ax(k) + Bu(k)
y(k) = Cx(k) + Du(k)

 x(0) = x0

• From a given initial condition x(0) and input sequence {u(k)}∞k=0 one can
predict the entire sequence of states x(k) and outputs y(k), ∀k ∈ N
• The state x(0) summarizes all the past history of the system
• The dimension n of the state x(k) ∈ Rn is called the order of the system
• The system is called proper (or strictly causal) if D = 0
• General multivariable case:
A ∈ Rn×n
x(k) ∈ Rn
B ∈ Rn×m
u(k) ∈ Rm
C ∈ Rp×n
y(k) ∈ Rp
D ∈ Rp×m

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 47

Equilibrium
• Consider the discrete-time nonlinear system
(
x(k + 1) = f (x(k), u(k))
y(k) = g(x(k), u(k))

Definition
A state xr ∈ Rn and an input ur ∈ Rm are an equilibrium pair if for initial
condition x(0) = xr and constant input u(k) ≡ ur , ∀k ∈ N, the state remains
constant: x(k) ≡ xr , ∀k ∈ N.

• Equivalent definition: (xr , ur ) is an equilibrium pair if f (xr , ur ) = xr

• xr is called equilibrium state, ur equilibrium input
• The definition generalizes to time-varying discrete-time nonlinear systems

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 48

Stability
• Consider the nonlinear system
(
x(k + 1) = f (x(k), ur )
y(k) = g(x(k), ur )
and let xr an equilibrium state, f (xr , ur ) = xr

Definition
The equilibrium state xr is stable if for each initial conditions x(0) “close
enough” to xr , the corresponding trajectory x(k) remains close to xr for all
k∈N
Analytic definition: ∀ϵ > 0 ∃δ > 0 : ∥x(0) − xr ∥ < δ ⇒ ∥x(k) − xr ∥ < ϵ, ∀k ∈ N.

• The equilibrium point xr is called asymptotically stable if it is stable and

x(k) → xr for k → ∞

• Otherwise, the equilibrium point xr is called unstable

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 49
Stability of first-order linear systems
• Consider the first-order linear system
x(k + 1) = ax(k) + bu(k)
• xr = 0, ur = 0 is an equilibrium pair

• For u(k) ≡ 0, ∀k = 0, 1, . . ., the solution is

x(k) = ak x0
• The origin xr = 0 is

3
– unstable if |a| > 1
2.5
– stable if |a| ≤ 1
2
– asymptotically stable if |a| < 1 a>1
x(k) 1.5
x0 a=1
1

0.5 0<a<1
xr a=0
0
0 2 4 6 8 10
k
©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 50
Stability of discrete-time linear systems
The natural response of x(k + 1) = Ax(k) + Bu(k) is x(k) = Ak x0 , so stability
only depend on A. We therefore talk about system stability
Theorem
Let λ1 , . . ., λm , m ≤ n be the eigenvalues of A ∈ Rn×n .
The system x(k + 1) = Ax(k) + Bu(k) is
• asymptotically stable iff |λi | < 1, ∀i = 1, . . . , m
• (marginally) stable if |λi | ≤ 1, ∀i = 1, . . . , m, and the eigenvalues with
unit modulus have equal algebraic and geometric multiplicity a
• unstable otherwise(in particular, if ∃ i such that |λi | > 1)
a Algebraic multiplicity of λ = number of coincident roots λ of det(λI − A). Geometric
i i
multiplicity of λi = number of linearly independent eigenvectors vi , Avi = λi vi

The stability properties of a discrete-time linear system only depend on

the modulus of the eigenvalues of matrix A

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 51

Stability of discrete-time linear systems
Proof:

• The natural response is x(k) = Ak x0

• If matrix A is diagonalizable2 , A = T ΛT −1 ,
 
λ 0 ... 0
 k
λ1 0 ... 0
1
0 λ2 ... 0  0 λk
2 ... 0 
    −1
. ⇒ A =T .
k
Λ= . . .
 . . . . T
.
.
. .. .
. . . .. . 
. . .
0 0 ... λn
0 0 ... λk
n

• Take any eigenvalue λ = ρejθ :

|λk | = ρk |ejkθ | = ρk

• A is always diagonalizable if algebraic multiplicity - geometric multiplicity □

• Lyapunov theorems also exist for nonlinear discrete-time systems (LaSalle, 1997)

2 If A is not diagonalizable, it can be transformed to Jordan form. In this case the natural response
x(t) contains modes kj λk , j = 0, 1, . . . , alg. multiplicity − geom. multiplicity

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 52

Zero eigenvalues
• Modes corresponding to λi =0 go to zero in finite-time

• This has no continuous-time counterpart, where instead all converging modes

tend to zero in infinite time (eλi t )

• Example: dynamics of a buffer

u(k)
x3(k) x2(k) x1(k) y(k)

     
 x1 (k + 1) = x2 (k) 
 0 1 0 0

 
    
 x (k + 1) = x3 (k)  x(k + 1) = 0 0 1  x(k) +  0  u(k)
2
  0 0 0 1
 x3 (k + 1)

= u(k) 
 [ ]
 

y(k) = x1 (k) y(k) = 1 0 0 x(k)
• Natural response: A3 x(0) = 0 for all x(0) ∈ R3

• For u(k) ≡ 0, the buffer deploys after at most 3 steps !

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 53
Exact sampling
• Consider the continuous-time system


 ẋ(t) = Ax(t) + Bu(t)
y(t) = Cx(t) + Du(t)

 x(0) = x0
• We want to characterize the value of x(t), y(t) at the time instants
t = 0, Ts , 2Ts , . . . , kTs , . . ., under the assumption that the input u(t) is
constant during each sampling interval (zero-order hold, ZOH)
y(t), y(kTs) u(t), u(kTs)
1.5 2

1.5

u(t) = ū(k), kTs ≤ t < (k + 1)Ts

1
1

0.5
0.5

• x̄(k) ≜ x(kTs ) and ȳ(k) ≜ y(kTs ) are the state 0

and the output samples at the k th sampling

0
−0.5

instant, respectively
−1
−0.5

−1.5

−1 −2
0 2 4 6 8 10 0 2 4 6 8 10
time t time t

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 54

Exact sampling
• Use Lagrange formula to get the response of the continuous-time system
between t0 = kTs and t = (k + 1)Ts from x(t0 ) = x(kTs ):
∫ t
x(t) = eA(t−t0 ) x(t0 ) + eA(t−σ) Bu(σ)dσ
t0
∫ (k+1)Ts
= eA((k+1)Ts −kTs ) x(kTs ) + eA((k+1)Ts −σ) Bu(σ)dσ
kTs

• Since the input u(t) is piecewise constant, u(σ) ≡ ū(k), kTs ≤ σ < (k + 1)Ts .
By setting τ = σ − kTs we get
Z !
Ts
x((k + 1)Ts ) = eATs x(kTs ) + eA(Ts −τ ) dτ Bu(kTs )
0

and hence
Z !
Ts
ATs A(Ts −τ )
x̄(k + 1) = e x̄(k) + e dτ B ū(k)
0

which is a linear difference relation between x̄(k) and ū(k) !

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 55

Exact sampling
• The discrete-time system
(
x̄(k + 1) = Āx̄(k) + B̄ ū(k)
ȳ(k) = C̄ x̄(k) + D̄ū(k)
depends on the original continuous-time system through the relations
Z !
Ts
A(Ts −τ )
Ā ≜ eATs , B̄ ≜ e dτ B, C̄ ≜ C, D̄ ≜ D
0

(if A is invertible then B̄ = (Ā − I)A−1 B )

• If u(t) is piecewise constant, (Ā, B̄, C̄, D̄) provides the exact evolution of state
and output samples at discrete times kTs y(t), y(kT )
0.4
s

0.2

MATLAB Python −0.2

−0.4
sys=ss(A,B,C,D); sys=[Link](A,B,C,D) 0 1 2 3 4 5
time t
6 7 8 9 10

sysd=c2d(sys,Ts); sysd=ctrl.c2d(sys,Ts) 1
u(t), u(kTs)

[Ab,Bb,Cb,Db]=ssdata(sysd); Ab,Bb,Cb,Db=[Link](sysd) 0.5

−0.5

−1

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 0 1 2 3 4 5

time t
6 7 8 9 10
56
Choice of sampling time

Rule of thumb: Ts ≈ 101

of rise time = time to move from 10% to
90% of the steady-state value, for input u(t) ≡ 1, x(0) = 0

rise time
y(t), y(kTs)
1 1
90 % 0.8
0.9

0.6
0.8
0.4
0.7
0.2

0.6 0
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
output

0.5
time t
u(t), u(kT )
0.4
s
2

0.3
1.5

0.2
1
10 %
0.1
0.5
rise time
0
0 1 2 3 4 5 6 7 8 9 10 0
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
time t
time t

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 57

Euler's forward method
x(t)
x((k + 1)Ts)
x((k+1)T
1)Ts ) − )¡x(kT
x((k +x((k+1)T x(kTs ) )
T
Ts

x((k + 1)Ts ) − x(kTs )

ẋ(kTs ) ≈ x(kT
x(kTs )
ẋ(kT
_
x(t) s)

Ts TTs
t
kT
kT s
(k+1)T
(k + 1)Ts Leonhard Paul Euler
(1707-1783)

• For nonlinear systems ẋ(t) = f (x(t), u(t)):

x̄(k + 1) = x̄(k) + Ts f (x̄(k), ū(k))

• For linear systems ẋ(t) = Ax(t) + Bu(t):

x((k + 1)Ts ) = (I + Ts A)x(kTs ) + Ts Bu(kTs )

Ā ≜ I + Ts A, B̄ ≜ Ts B, C̄ ≜ C, D̄ ≜ D
Tsn An
• e Ts A = I + T s A + . . . + n!
+ . . . Euler’s method ≈ exact sampling for Ts → 0

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 58

Eigenvalues mapping

• Let λi = eigenvalues of matrix A (continuous-time system), i = 1, . . . , n.

Assume A diagonalizable, A = T ΛT −1

• The eigenvalues of eTs A = T eTs Λ T −1 are eTs λi

<λi < 0 → |eTs λi | < 1

• The eigenvalues of I + Ts A = T (I + Ts Λ)T −1 are 1 + Ts λi

<λi < 0 6→ |1 + Ts λi | < 1 !

Euler’s forward method can make an asymptotically stable

continuous-time system unstable if Ts is not small enough!

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 59

Example - Hydraulic system
Continuous-time model
( √ p
d
dt h(t) = − a A2g h(t) + 1
A u(t)
√ p
q(t) = a 2g h(t)

Discrete-time model
( √ p
h̄(k + 1) = h̄(k) − Ts a 2g
h̄(k) + TAs ū(k)
√ p A
q̄(k) = a 2g h̄(k) (Torricelli's law)

level h(t) (m)

7

u 6

h 2

1
continuous time
q Euler approximation
0
0 5 10 15 20 25 30
time (s)

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 60

N -steps Euler method
• We can obtain the matrices A, B of the discrete-time linearized model while
integrating the nonlinear continuous-time dynamic equations ẋ = f (x, u)

• N -steps explicit forward Euler method: given x̄(k), ū(k), execute the following
steps
1. x = x̄(k), Ā = I , B̄ = 0
2. for n=1,2,…,N do
Ts ∂f
• Ā ← (I + N ∂x
(x, ū(k))Ā
Ts ∂f Ts ∂f
• B̄ ← (I + N ∂x
(x, ū(k))B̄ + N ∂u
(x, ū(k))
• x ← x + TNs f (x, ū(k))
3. end
4. return x̄(k + 1) ≈ x and matrices Ā, B̄ such that x̄(k + 1) ≈ Ax̄(k) + B ū(k).

• Property: the difference between the state x̄(k + 1) and its approximation x


computed by the above iterations satisfies kx̄(k + 1) − x)k = O TNs

• Explicit forward Runge-Kutta 4 method also available

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 61

Tustin's discretization method
• Assume u(t) constant within the sampling interval. Given the linear system
ẋ(t) = Ax(t) + Bu(t), apply the trapezoidal rule to approximate the integral
Z (k+1)Ts Z (k+1)Ts
x̄(k + 1) − x̄(k) = ẋ(t)dt = (Ax(t) + Bu(t))dt
kTs kTs
Ts
≈ (Ax̄(k) + B ū(k) + Ax̄(k + 1) + B ū(k)) (trapezoidal rule)
2
and therefore
Ts Ts
(I − A)x̄(k + 1) = (I + )x̄(k) + Ts B ū(k)
2 2
 −1    −1
Ts Ts Ts
x̄(k + 1) = I − A I+ A x̄(k) + I − A Ts B ū(k)
2 2 2

• Advantage: simpler to compute than exponential matrix, without too much loss
of approximation quality

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 62

Z-transform
Consider a function f (k), f : Z → R, f (k) = 0 for all k < 0

Definition
2

The unilateral Z-transform of f (k) is the function of 1.8

the complex variable z ∈ C defined by

1.6

1.4

1.2

f(k)
∞
X
1

−k
0.8

F (z) = f (k)z 0.6

0.4

k=0 0.2

0
−2 0 2 4 6 8 10

Once F (z) is computed using the series, it’s

extended to all z ∈ C for which F (z) makes sense

Z-transforms convert difference equations into

Witold Hurewicz
(1904-1956) algebraic equations.

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 63

Examples of Z-transforms

• Discrete impulse

0 6 0
if k =
f (k) = δ(k) ≜ ⇒ Z[δ] = F (z) = 1
1 if k = 0
• Discrete step

0 if k < 0 z
f (k) = 1I(k) ≜ ⇒ Z[1I] = F (z) =
1 if k ≥ 0 z−1

• Geometric sequence
z
f (k) = ak 1I(k) ⇒ Z[f ] = F (z) =
z−a

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 64

Properties of Z-transforms

• Linearity

Z[α1 f1 (k) + α2 f2 (k)] = α1 Z[f1 (k)] + α2 Z[f2 (k)]

Example: f (k) = 3δ(k) − 5

2k
1I(t) ⇒ Z[f ] = 3 − 5z
z− 12

• Forward shift3
Z[f (k + 1) 1I(k)] = zZ[f ] − zf (0)

Example: f (k) = ak+1 1I(k) ⇒ Z[f ] = z z−a

z
−z = az
z−a

3z is also called forward shift operator

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 65

Properties of Z-transforms

• Backward shift or unit delay 4

Z[f (k − 1) 1I(k)] = z −1 Z[f ]

Example: f (k) = 1I(k − 1) ⇒ Z[f ] = z

z(z−1)

• Multiplication by k
d
Z[kf (k)] = −z Z[f ]
dz

Example: f (k) = k 1I(k) ⇒ Z[f ] = z

(z−1)2

4 z −1 is also called backward shift operator

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 66

Discrete-time transfer function
Apply forward-shift & linearity rules to x(k + 1) = Ax(k) + Bu(k), and
linearity to y(k) = Cx(k) + Du(k):
X(z) = z(zI − A)−1 x0 + (zI − A)−1 BU (z)
Y (z) = zC(zI − A)−1 x0 + (C(zI − A)−1 B + D)U (z)
| {z } | {z }
Z-transform of natural response Z-transform of forced response

Definition
The transfer function of the discrete-time linear system (A, B, C, D) is

G(z) = C(zI − A)−1 B + D

that is the ratio between the Z-transform Y (z) of the output and the Z-
transform U (z) of the input signals for the initial state x0 = 0

MATLAB Python
sys=ss(A,B,C,D,Ts); sys=[Link](A,B,C,D,Ts)
G=tf(sys) G=[Link](sys)
©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 67
Discrete-time transfer function
u(k) y(k) U (z) Y (z)
A; B; C; D G(z)

x0 = 0
  0.5 
 x(k + 1) = 1
0 −0.5 x(k) + [ 01 ] u(k)
Example: The linear system
 y(k) = [ 1 −1 ] x(k)

with sampling time Ts = 0.1 s has the transfer function

MATLAB Python
Ts=0.1; import numpy as np
−z + 1.5 Ts=0.1
G(z) = 2 A=[0.5 1;0 -0.5]; A=[Link]([[0.5, 1],[0, -0.5]])
z − 0.25 B=[0;1];
C=[1 -1];
B=[Link]([[0],[1]])
C=[Link]([[1, -1]])
sys=ss(A,B,C,0,Ts); sys=[Link](A,B,C,0,Ts)
G=tf(sys) G=[Link](sys)

The transfer function does not depend Transfer function: Transfer function:
-z + 1.5 -z + 1.5
on the input u(k), it is only a property ————— —————
z^2 - 0.25 z^2 - 0.25
of the linear system.

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 68

Difference equations

• Consider the nth -order difference equation forced by u

an y(k − n) + an−1 y(k − n + 1) + · · · + a1 y(k − 1) + y(k)

= bn u(k − n) + · · · + b1 u(k − 1)

• For zero initial conditions we get the transfer function

bn z −n + bn−1 z −n+1 + · · · + b1 z −1
G(z) =
an z −n + an−1 z −n+1 + · · · + a1 z −1 + 1
b1 z n−1 + · · · + bn−1 z + bn
=
z n + a1 z n−1 + · · · + an−1 z + an

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 69

Impulse response

• Consider the impulsive input u(k) = δ(k), U (z) = 1. The corresponding output
y(k) is called impulse response

• The Z-transform of y(k) is Y (z) = G(z) · 1 = G(z)

• Therefore the impulse response coincides with the inverse Z-transform g(k) of
the transfer function G(z)

Example (integrator:)
2 2

1 1

u(k)

y(k)
u(k) = δ(k) h i
0 0

Z −1 1
= 1I(k − 1)
−1 −1
y(k) = z−1 −2 −2
−2 0 2 4 6 8 10 −2 0 2 4 6 8 10
k k

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 70

Poles, eigenvalues, modes
• Linear discrete-time system
(
x(k + 1) = Ax(k) + Bu(k)
NG (z)
y(k) = Cx(k) + Du(k) G(z) = C(zI−A)−1 B+D ≜
DG (z)
x(0) = 0
• Use the adjugate matrix to represent the inverse of zI − A

C Adj(zI − A)B
C(zI − A)−1 B + D = +D
det(zI − A)

• The denominator DG (z) = det(zI − A) !

The poles of G(z) coincide with the eigenvalues of A

• Well, not always ... There might be a zero/pole cancellation (we will see later)

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 71

Steady-state solution and DC gain
• Let A asymptotically stable (|λi | < 1). The natural response vanishes
asymptotically
• Assume constant u(k) ≡ ur , ∀k ∈ N. What is the asymptotic value
xr = limk→∞ x(k) ?

Impose xr (k + 1) = xr (k) = Axr + Bur and get xr = (I − A)−1 Bur

The corresponding steady-state output yr = Cxr + Dur is

yr = (C(I − A)−1 B + D) ur
| {z }
DC gain
• Cf. final value theorem in complex analysis:
yr = lim y(k) = lim (z − 1)Y (z) = lim (z − 1)G(z)U (z)
k→+∞ z→1 z→1
ur z
= lim (z − 1)G(z) = G(1)ur = (C(I − A)−1 B + D)ur
z→1 z−1
• G(1) is called the DC gain of the system

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 72

Example - Student population dynamics
• Recall student population dynamics
    

 .2 0 0 1

    
 x(k + 1) =  .6 .15 0  x(k) +  0  u(k)

 h 0 .8 i.08 0


 y(k) = 0 0 .9 x(k)

• DC gain:
h 1 0 0 i h .2
0 0 −1 1 i h i
[ 0 0 .9 ] 010 −.6 .15 0 0 ≈ 0.69
001 0 .8 .08 0
• Transfer function: G(z) = 0.432
z −0.43z +0.058z−0.0024
3 2 , G(1) ≈ 0.69
y(k)
35

30
MATLAB Python
A=[b1 0 0; a1 b2 0; 0 a2 b3]; A=[[b1, 0, 0],[a1, b2, 0],[0, a2, b3]]
25
B=[1;0;0]; B=[[1],[0],[0]]
20 C=[0 0 a3]; C=[0, 0, a3]
D=[0]; D=[0]
15
sys=ss(A,B,C,D,1); sys=[Link](A,B,C,D,1)
10 dcgain(sys) [Link](sys)
5
0.6905 0.6905
0
2006 2008 2010 2012 2014 2016
step k

• For u(k) ≡ 50 students enrolled steadily, y(k) → 0.6905 · 50 ≈ 34.5 graduates

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 73
Closed-loop control
Proportional integral derivative (PID) controllers
• PID (proportional integrative derivative) controllers are the most used
controllers in industrial automation since the ’30s
h Z
1 t
de(t) i
u(t) = Kp e(t) + e(τ )dτ + Td
Ti 0 dt
where e(t) = r(t) − y(t) is the tracking error

• Initially constructed by analog electronic components, today they are

implemented digitally
– ad hoc digital devices
– just few lines of C code included in the control unit

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 74

PID parameters

r + e Z t ++ u y
Process
1
Kp
-
e(⌧ )d⌧
Ti 0
+
de(t)
Td
dt
Controller

• Kp is the controller gain, determining the “aggressiveness” of the controller

• Ti is the reset time, determining the weight of the integral action. The integral
action guarantees that in steady-state y(t) = r(t)
de(t)
• Td is the derivative time. The term e(t) + Td dt provides a “prediction” of the
tracking error at time t + Td
• We call the controller P, PD, PI, or PID depending on the feedback terms
included in the control law

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 75

Structure of PID controller
• In practice one implements the following version of the PID controller
h Z i
1 t
u(t) = Kp br(t) − y(t) + (r(τ ) − y(τ ))dτ + d(t)
| {z } Ti 0 |{z}
| {z }
proportional derivative
integral
action action
action
Td ˙
d(t) + d(t) = −Td ẏ(t)
N
• the reference signal r(t) is not included in the derivative term (r(t) may have
abrupt changes)

• the proportional action Kp (br(t) − y(t) only uses a fraction b ≤ 1 of the

reference signal r(t)

• the derivative term d(t) is a filtered version of ẏ(t)

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 76

Digital implementation of PID controller

• In digital (=discrete-time) form with sampling time Ts , the PID controller takes
the following form

u(k) = P (k) + I(k) + D(k)

P (k) = Kp (br(k) − y(k))

Kp T s
I(k + 1) = I(k) + (r(k) − y(k)) forward differences
Ti

Td Kp Td N
D(k) = D(k − 1) − (y(k) − y(k − 1))
Td + N T s Td + N T s
backward differences

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 77

PID controller: pros and cons

• Very simple to implement, only 3 parameters to calibrate

• It only requires the measurement of the output signal y(t)

• The control law does not exploit the knowledge of the model of the process

• Achievable closed-loop performance is limited

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 78

State-feedback control
Reachability analysis
• Consider the linear discrete-time system

x(k + 1) = Ax(k) + Bu(k)

with x ∈ Rn , u ∈ Rm and initial condition x(0) = x0 ∈ Rn

X
k−1
• The solution is x(k) = Ak x0 + Aj Bu(k − 1 − j)
j=0
Definition
The system x(k + 1) = Ax(k) + Bu(k) is (completely) reachable if ∀x1 , x2 ∈
Rn there exist k ∈ N and u(0), u(1), . . ., u(k − 1) ∈ Rm such that

X
k−1
x 2 = Ak x 1 + Aj Bu(k − 1 − j)
j=0

• In simple words: a system is completely reachable if from any state x1 we can

reach any state x2 at some time k , by applying a suitable input sequence

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 79

Reachability
• Determine a sequence of n inputs transferring the state vector from x1 to x2
after n steps
 
u(n − 1)
 
  u(n − 2) 
x2 − A x1 = B AB . . . A
n n−1
B  
} 
| {z } | ..
{z  . 
X R
u(0)
| {z }
U

• This is equivalent to solve with respect to U the linear system of equations

RU = X

• Matrix R ∈ Rn×nm is called the reachability matrix of the system

• A solution U exists if and only if X ∈ Im(R)
(Rouché-Capelli theorem: a solution exists ⇔ rank([R X]) = rank(R))

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 80

Reachability
Theorem

The system (A, B) is completely reachable ⇔ rank(R) = n

Proof:
(⇒) Assume (A, B) reachable, choose x1 = 0 and x2 = x. Then ∃k ≥ 0 such
that
X
k−1
x= Aj Bu(k − 1 − j)
j=0

If k ≤ n, then clearly x ∈ Im(R). If k > n, by Cayley-Hamilton theorem we

have again x ∈ Im(R). Since x is arbitrary, Im(R) = Rn , so rank(R) = n.

(⇐) If rank(R) = n, then Im(R) = Rn . Let X = x2 − An x1 and

′
U = [ u(n − 1)′ . . . u(1)′ u(0)′ ] . The system X = RU can be solved with
respect to U , ∀X , so any state x1 can be transferred to x2 in k = n steps.
Therefore, the system (A, B) is completely reachable.

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 81

Minimum-energy control
• Let (A, B) reachable and consider steering the state from x(0) = x1 into
x(k) = x2 , k > n  
u(k − 1)
h i
 u(k − 2)


x2 − Ak x1 = B AB . . . Ak−1 B 
 .. 

| {z }
X | {z } . 
Rk
u(0)
| {z }
U

(Rk ∈ Rn×km is the reachability matrix for k steps)

• Since rank(Rk ) = rank(R) = n, ∀k > n (Cayley-Hamilton), we get
rank Rk = rank[Rk X] = n
• Hence the system X = Rk U admits solutions U
Problem
Determine the input sequence {u(j)}k−1
j=0 that brings the state from x(0) = x1
1X
k−1
1
to x(k) = x2 with minimum energy ku(j)k2 = U ′ U
2 j=0 2

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 82

Minimum-energy control
• The problem is equivalent to finding the solution U of the system of equations
X = Rk U
with minimum norm kU k
• We must solve the optimization problem
1
U ∗ = arg min
2
kU k subject to X = Rk U
2
• Let’s apply the method of Lagrange multipliers:
1
L(U, λ) = kU k + λ′ (X − Rk U ) Lagrangian function
2
2

∂L
∂U = U − Rk′ λ = 0 MATLAB
⇒ U ∗ = Rk′ (Rk Rk′ )−1 ·X U=pinv(Rk)*X
∂L
= X − Rk U = 0 | {z }
∂λ Python
Rk# = pseudoinverse from [Link] import pinv
U=pinv(Rk)@X

• Note that Rk Rk′ is invertible because rank(Rk ) = rank(R) = n, ∀k ≥ n

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 83
Controllability

• If the system is completely reachable, we have seen that we can bring the state
vector from any value x(0) = x1 to any other value x(n) = x2

• Let’s focus on the subproblem of determining a finite sequence of inputs that

brings the state to the final value x(n) = 0

Definition
A system x(k + 1) = Ax(k) + Bu(k) is controllable to the origin in k steps
if ∀x0 ∈ Rn there exists a sequence u(0), u(1), . . ., u(k − 1) ∈ Rm such that
Pk−1
0 = Ak x 0 + j=0 Aj Bu(k − 1 − j)

• Controllability is a weaker condition than reachability

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 84

Controllability, Stabilizability
• The linear system of equations
 
u(k − 1)
 
   u(k − 2) 
−An x0 = B AB . . . An−1 B  
} 
..
| {z  . 
R
u(0)
admits a solution if and only if An x0 ∈ Im(R), ∀x0 ∈ Rn
Theorem
The system is controllable to the origin (in n steps) if and only if

Im(An ) ⊆ Im(R)

Definition
A linear system x(k + 1) = Ax(k) + Bu(k) is called stabilizable if can be
driven asymptotically to the origin

• Stabilizability is a weaker condition than controllability

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 85
Canonical reachability decomposition
• Goal: Make a change of coordinates to separate reachable from unreachable
states

• Let rank(R) = nc < n and consider the change of coordinates

h i
T = wnc +1 ... wn v1 ... v nc

where {v1 , . . . , vnc } is a basis of Im(R), and {wnc +1 , . . . , wn } is a completion

to obtain a basis of Rn (i.e., a basis of ker(R′ ), R′ wi = 0)

• As Im(R) is A-invariant (Ax ∈ Im(R), ∀x ∈ Im(R), follows from

Cayley-Hamilton theorem), Avi has no components along the basis vectors
wnc +1 , …, wn

• Since T −1 Avi are the new coordinates of Avi , the first n − nc components of
T −1 Avi are zero

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 86

Canonical reachability decomposition

• The columns of B also have zero components along wnc +1 , …, wn , because

Im(B) ⊆ Im(R)

• In the new coordinates, the system has matrices à = T −1 AT , B̃ = T −1 B and

C̃ = CT in the canonical reachability form
(a.k.a. controllability staircase form)
" # " # MATLAB
Auc 0 0 h i
[At,Bt,Ct,Tinv]=
à = B̃ = C̃ = Cuc Cc ctrbf(A,B,C)
A21 Ac Bc

• Let x = [ xxuc
c
] be the coordinates of the state vector in the new coordinate
system, xuc ∈ Rn−nc , xc ∈ Rnc

• We have that xuc (k) = Akuc xuc (0), so xuc (k) does not depend on u(k)

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 87

Reachability and transfer function
Proposition
The eigenvalues of Auc are not poles of the transfer function
C(zI − A)−1 B + D

Proof: Let T transform (A, B) to canonical reachability decomposition (Ã, B̃).

The transfer function is
G(z) = C(zI − A)−1 B + D = C̃(zI − Ã)−1 B̃ + D
( [ ])−1 [ ]
[ ] Auc 0 0
= Cuc Cc zI − +D
A21 Ac Bc
[ ][ ]
[ ] (zI − A )−1 0 0
uc
= Cuc Cc +D
⋆ (zI − Ac )−1 Bc
= Cc (zI − Ac )−1 Bc + D

Clearly G(z) does not depend on the eigenvalues of Auc

Lack of reachability → zero/pole cancellations!

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 88

Reachability and transfer function
• Why are the eigenvalues of Auc not appearing in the transfer function G(z) ?
• Remember: G(z) explains the forced response, i.e., the response for x(0) = 0
• Expressed in canonical decomposition, the system evolution is


 xuc (k + 1) = Auc xuc (k)
xc (k + 1) = Ac xc (k) + Bc u(k) + A21 xuc (k)

 y(k) = Cuc xuc (k) + Cc xc (k) + Du(k)

• For xuc (0) = 0, xc (0) = 0, we get xuc (k) ≡ 0 and



 xc (k + 1) = Ac xc (k) + Bc u(k)
y(k) = Cc xc (k) + Du(k)

 xc (0) = 0

so the forced response does not depend at all on Auc !

• The input u(k) only affects the output y(k) through the reachable subsystem
(Ac , Bc , Cc , D), not through the unreachable part Auc

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 89

Reachability analysis of continuous-time systems
• Similar definitions of reachability, controllability, and stabilizability can be given
for continuous-time systems

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

• No distinction between controllability and reachability in continuous-time

(because no finite-time convergence of modal response exists)

• Reachability matrix and canonical reachability decomposition are identical to

discrete-time

• rank R = n is also a necessary and sufficient condition for reachability

• Auc asymptotically stable (all eigenvalues with negative real part) is also a
necessary and sufficient condition for stabilizability

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 90

Stabilization by state feedback
• Main idea: design a device that makes the process (A, B, C) asymptotically
stable by manipulating the input u to the process
!"#$%&'$()*+,'-..

v(k) + u(k) x(k) y(k)

+
A,B C

state feed
K back

• If measurements of the state vector are available, we can set

u(k) = k1 x1 (k) + k2 x2 (k) + . . . + kn xn (k) + v(k)

• v(k) is an exogenous signal exciting the closed-loop system

Problem
Find a feedback gain K = [k1 k2 . . . kn ] that makes the closed-loop system
asymptotically stable.

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 91

Stabilization by state feedback

v(k) + u(k) x(k) y(k)

+
A,B C

K
!"#$%&'"##()$*$+%,

• Let u(k) = Kx(k) + v(k). The overall system is

x(k + 1) = (A + BK)x(k) + Bv(k)

y(k) = (C + DK)x(k) + Dv(k)

Theorem
 
(A, B) ”reachable” (rank B AB . . . An−1 B = n) ⇒ the eigenvalues of
(A + BK) can be decided arbitrarily.

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 92

Eigenvalue assignment problem
Fact
(A, B) reachable ⇔ (A, B) is algebraically equivalent to a pair (Ã, B̃) in
controllable canonical form
   
0 0
 ..   .. 

à =  .
In−1   
, B̃ =  . 
 0  0
−a0 −a1 ... −an−1 1

The transformation matrix T such that à = T −1 AT , B̃ = T −1 B is

 
a1 a2 ... an−1 1
 a2 a3 ... 1 0
 
 . . . . . 
T = [B AB . . . An−1 B]
 .. .
.
.
.
.
.
. 
. 
 
 an−1 1 0 ... 0
1 0 0 ... 0

where a1 , a2 , . . ., an−1 are the coefficients of the characteristic polynomial

pA (λ) = λn + an−1 λn−1 + . . . + a1 λ + a0 = det(λI − A)
©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 93
• Let (A, B) reachable and assume m = 1 (single input)
• Characteristic polynomials:

pA (λ) = λn + an−1 λn−1 + . . . + a1 λ + a0 (open-loop eigenvalues)

n n−1
pd (λ) = λ + dn−1 λ + . . . + d1 λ + d0 (desired closed-loop eigenvalues)

• Let (A, B) be in controllable canonical form

   
0 0
 .
. In−1   .
. 
A=

. , B = 
 
. 

0 0
−a0 −a1 ... −an−1 1

• As K = [k1 . . . kn ], we have
 
0
 .
. In−1 
A + BK = 

. 

0
−(a0 − k1 ) −(a1 − k2 ) ... −(an−1 − kn )

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 94

• The characteristic polynomial of A + BK is therefore

λn + (an−1 − kn )λn−1 + . . . + (a1 − k2 )λ + (a0 − k1 )

• To match pd (λ) we impose

a0 − k1 = d0 , a1 − k2 = d1 , . . . , an−1 − kn = dn−1

Procedure
If (A, B) is in controllable canonical form, the feedback gain
h i
K = a0 − d0 a1 − d1 . . . an−1 − dn−1

makes pd (λ) the characteristic polynomial of (A + BK)

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 95

• If (A, B) is not in controllable canonical form we must set
h i
K̃ = a0 −d0 a1 −d1 . . . an−1 −dn−1
K = K̃T −1 ← don't invert T, solve instead T ′ K ′ = K̃ ′ w.r.t. K′ !

where  
a1 a2 ... an−1 1
 a2 a3 ... 1 0
 
T = R

.
.
.
.
.
.
.
.
.
.
.
.
. 
.
. 
a 1 0 ... 0

n−1
1 0 0 ... 0
−1
• Explanation: a matrix M and T M T have the same eigenvalues

det(λI − T −1 M T ) = det(T −1 T λ − T −1 M T ) = det(T −1 ) det(λI − M )

· det(T ) = det(λI − M )

• Since (Ã + B̃ K̃) = T −1 AT + T −1 BKT = T −1 (A + BK)T , it follows that

(Ã + B̃ K̃) and (A + BK) have the same eigenvalues

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 96

Ackermann's formula
• Let (A, B) reachable and assume m = 1 (single input)
• Characteristic polynomials:

pA (λ) = λn + an−1 λn−1 + . . . + a1 λ + a0 (open-loop eigenvalues)

n n−1
pd (λ) = λ + dn−1 λ + . . . + d1 λ + d0 (desired closed-loop eigenvalues)

• Let pd (A) = An + dn−1 An−1 + . . . + d1 A + d0 I ← (This is n × n matrix !)

MATLAB
Ackermann’s formula
K=-acker(A,B,P);
K = −[0 . . . 0 1][B AB . . . An−1 B]−1 pd (A) K=-place(A,B,P);

Python
where P = [λ1 λ2 . . . λn ] are the desired
K=-[Link](A,B,P)
closed-loop poles K=-[Link](A,B,P)

• Numerically robust methods to solve the pole assignment problem exist

(Tits, Yang, 1996)

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 97

Zeros of closed-loop system
Fact
In case of no feedthrough (D = 0) the zeros of the system are the same under
state feedback
NK (z) = N (z)

• Example for x ∈ R3 : change the coordinates to canonical reachability form

   
0 1 0 0 [ ]
   
A= 0 0 1  , B =  0  , K = k3 k2 k1
−a3 −a2 −a1 1
    
z 2 + a1 z + a2 z + a1 1 0 1
    
Adj(zI − A)B =  −a3 z(z + a1 ) z  0  =  z 
−a3 z −a2 z − a3 z 2 1 z 2

• Adj(zI − A)B does not depend on the coefficients a1 , a2 , a3

• Hence also Adj(zI − A − BK)B does not depend on a1 − k1 , a2 − k2 , a3 − k3

N (z) = C Adj(zI−A)B = C Adj(zI−A−BK)B = NK (z), ∀K ′ ∈ Rn

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 98

Example - Student population dynamics
• The open-loop poles are (0.8, 0.15, 0.2)
• Say we want to place the closed-loop poles in (0.1 ± 0.2j, 0.1) by setting

u(k) = Kx(k) + Hr(k)

where r(k) is the desired reference signal

• First, design K by pole placement:
MATLAB Python
K=-place(A,B,[.1+.2*j,.1-.2*j,.1]) K=-place(A,B,[.1+.2j,.1-.2j,.1])

• Then choose H such that the DC-gain from r to y is 1:

MATLAB Python
sys_cl=ss(A+B*K,B,C+D*K,D,1); sys_cl=[Link](A+B@K,B,C+D@K,D,1)
dc_cl=dcgain(sys_cl); dc_cl=[Link](sys_cl)
H=1/dc_cl; H=1/dc_cl

• We get K = [−0.1300 − 0.0698 0.0017], H = 1.7708

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 99

Example - Student population dynamics
• Compare open-loop vs. closed-loop response
60
open-loop
50 closed-loop
reference
40

30

20

10

0
2012 2013 2014 2015 2016 2017 2018 2019 2020 2021 2022

90

85

80

75

70
2012 2013 2014 2015 2016 2017 2018 2019 2020 2021 2022

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 100

State estimation
Observability
!"#$%&'$()*+,'-..

u(k) x(k) y(k)

A,B C

? ˆ
x(k) ./$/-)-./&%$/-

• Implementing a state feedback controller u(k) = Kx(k) requires the entire

state vector x(k)
• Problem: often sensors only provide the measurements of output y(k)
• Idea: is it possible to estimate the state x by measuring only the output y and
knowing the applied input u ?
• Observability analysis addresses this problem, telling us when and how the
state estimation problem can be solved

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 101

Observability
(
• Consider x(k + 1) = Ax(k) + Bu(k)
y(k) = Cx(k) + Du(k)

with5 x ∈ Rn , u ∈ R, y ∈ R and initial condition x(0) = x0 ∈ Rn

• The solution for the output is
Pk−1
y(k, x0 , u(·)) = CAk x0 + j=0 CAj Bu(k − 1 − j) + Du(k)

Definition
The pair of states x1 6= x2 ∈ Rn is called indistinguishable from the output
y(·) if for any input sequence u(·)

y(k, x1 , u(·)) = y(k, x2 , u(·)), ∀k ≥ 0

A linear system is called (completely) observable if no pair of states are indis-

tinguishable from the output
5 Everything here can be easily generalized to multivariable systems u ∈ Rm , y ∈ R p

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 102

Observability
• Consider the problem of reconstructing the initial condition x0 from n output
measurements, applying a known input sequence

y(0) = Cx0 + Du(0)

y(1) = CAx0 + CBu(0) + Du(1)
..
.
Pn−2
y(n − 1) = CAn−1 x0 + j=1 CAj Bu(n − 2 − j) + Du(n − 1)

• Define
 
  y(0) − Du(0)
C  
   y(1) − CBu(0) − Du(1) 
 CA   .. 
Θ=
 .. 
 Y =
 . 

 .   X
n−2 
 
CAn−1 y(n − 1) − CAj Bu(n − 2 − j) − Du(n − 1)
| {z } j=1
n × n matrix | {z }
n-th dimensional vector

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 103

Observability
• The initial state x0 is determined by solving the linear system

Y = Θx0

The matrix Θ ∈ Rn×n is called the observability matrix of the system

• If we assume perfect knowledge of the output (i.e., no noise on output

measurements), we can always solve the system Y = Θx0 . In particular:

– There is only one solution if rank(Θ) = n

– There are infinite solutions if rank(Θ) < n.

In this case, all solutions are given by x0 + ker(Θ), where x0 is any particular
solution of the system (e.g., the true initial state)

Pk−1
• Knowing x0 , we know x(k) = Ak x0 + i=0 Ai Bu(k − 1 − i) for all k ≥ 0

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 104

Observability

• The system of equations Θx0 = Y has a solution if and only if

rank(Θ) = rank([Θ Y ]) (Rouché-Capelli Theorem)

• Because we have Θ ∈ Rn×n , if rank(Θ) = n ⇒ rank([Θ Y ]) = n for each Y

• The solution is unique if and only if rank(Θ) = n

• The input u(k) only influences Y , not Θ

• Then, for linear systems the observability property only depends on A and C

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 105

Observability
Theorem

A linear system is observable if and only if rank(Θ) = n

• As the observability property of a system depends only on matrices A and C ,

we call a pair (A, C) observable if
 
C
 
 CA 
rank 

 ..  = n

 . 
CAn−1
• It can be proved that ker(Θ) is the set of states x ∈ Rn that are
indistinguishable from the origin x = 0
y(k, x, u(·)) = y(k, 0, u(·)), ∀k ≥ 0
for any input sequence u(·)
• Hence, since ker(Θ) = {0} if and only if rank(Θ) = n, a system is observable if
and only if there are no states that are indistinguishable from x = 0
©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 106
Reconstructability

• Under observability assumptions, we just saw that it is possible to determine

the initial condition x0 from n input/output measurements

x(0) = Θ−1 Y

• To close the control loop at time k it is enough to know the current x(k)

• If the initial condition x(0) is known, it is possible to calculate x(k) as

X
k−1
x(k) = Ak Θ−1 Y + Ai Bu(k − 1 − i)
i=0

• Question: Can we determine the current state x(k) even if the system is not
completely observable?

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 107

Reconstructability
Definition
A linear system x(k +1) = Ax(k)+Bu(k) is called reconstructable in k steps
if, for each initial condition x0 , x(k) is uniquely determined by {u(j), y(j)}k−1
j=0

The solutions of the system

 
y(0) − Du(0)  
 y(1) − CBu(0) − Du(1)  C
  
.   CA  
.   x
Yk ≜  . = . 
  .
 X
k−2   . 
 y(k − 1) − CAj Bu(k − 2 − j) + Du(k − 1)  CAk−1
j=1 | {z }
Θk

are given by x = x0 + ker(Θk ), where x0 is the “true” (unknown) initial state

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 108

Reconstructability
• Let x0 be the initial (unknown) “true” state, and x = x0 + x̄ be a generic initial
state, where x̄ ∈ ker(Θk ). An estimation x̂(k) of the current state x(k) is

X
k−1
k
x̂(k) = A x0 + A x̄ + k
Aj Bu(k − 1 − j)
j=1

• x̂(k) coincides with x(k) if and only if x̄ ∈ ker(Ak ). Because this must hold for
any x̄ ∈ ker(Θk ), we have the following

Lemma

A system is reconstructable in k steps if and only if ker(Θk ) ⊆ ker(Ak )

Definition

A system is detectable if it is reconstructable asymptotically for k → +∞

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 109

Canonical observability decomposition
Goal: Change coordinates to separate observable and unobservable states

• Let dim(ker(Θ)) = n − no ≥ 1 and consider the change of coordinates

h i
T = vno +1 ... vn w1 ... w no

where {vno +1 , . . . , vn } is a basis of ker(Θ), and {w1 , . . . , wno } is a completion

to obtain a basis of Rn
• By Cayley-Hamilton theorem, ker(Θ) is A-invariant (Ax ∈ ker(Θ), ∀x ∈ Θ),
and hence Avi has no components along the basis vector w1 , . . . , wno ,
∀i = no + 1, . . . , n
• Note also that Cvi = 0, because Θvi = 0, ∀i = no + 1, . . . , n
• In the new coordinates the system has matrices à = T −1 AT , B̃ = T −1 B and
C̃ = CT in the canonical observability form
" # " # MATLAB
Auo A12 Buo h i
[At,Bt,Ct,Tinv]=
à = B̃ = C̃ = 0 Co obsvf(A,B,C)
0 Ao Bo

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 110

Observability and transfer function
Proposition
The eigenvalues of Auo are not poles of the transfer function
C(zI − A)−1 B + D

Proof: Consider a matrix T changing the state coordinates to canonical

observability decomposition of (A, C). The transfer function is
G(z) = C(zI − A)−1 B + D = C̃(zI − Ã)−1 B̃ + D =
( [ ])−1 [ ]
[ ] Auo A12 Buo
0 Co zI − +D
0 Ao Bo
[ ][ ]
[ ] (zI − A )−1 ⋆ Bno
uo
= 0 Co +D
0 (zI − Ao )−1 Bo
= Co (zI − Ao )−1 Bo + D

Clearly G(z) does not depend on the eigenvalues of Auo

Lack of observability → zero/pole cancellations!

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 111

Observability and transfer function

• Why are the eigenvalues of Auo not appearing in the transfer function G(z) ?

• Expressed in canonical decomposition, the system evolution is



 xuo (k + 1) = Auo xuo (k) + A12 xo (k) + Buo u(k)
xo (k + 1) = Ao xo (k) + Bo u(k)

 y(k) = Co xo (k) + Du(k)

• The evolution of xo (k) is not affected by the unobservable states xuo (k)

X
k−1
xo (k) = Ako xo (0) + Aio Bo u(k − 1 − i)
i=0

so the output y(k) = Co xo (k) + Du(k) does not depend at all on Auo !

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 112

Canonical observability decomposition
Proposition
Ao ∈ Rno ×no and Co ∈ Rp×no are a completely observable pair
Proof:
   
• We have that CT C̃
 CT T −1 AT   C̃ Ã 
   
no = rank Θ = rank ΘT = rank 
 .  = rank 
  . 

 .   . 
. .
−1
CT T An−1 T C̃ Ãn−1
 
0 Co
0 
   
 . .  Co
 . . 
 . .   Co A o 
   
= rank 
0 Co A n
o
o −1  = rank 
  . 

0   . 
 Co A n o
o
 .
  no −1
 .. .
.  Co A o
 . . 
0 Co An−1o

• The last equality follows by Cayley-Hamilton theorem (the last n − no rows

[0 C0 Aio ] are a linear combination of the first no rows). Hence, (Ao , Co ) is
completely observable
©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 113
Duality
• Given a linear system (A, B, C, D), with x ∈ Rn , u ∈ Rm and y ∈ Rp , we call
dual system the system
(
x̃(k + 1) = A′ x̃(k) + C ′ ũ(k)
ỹ(k) = B ′ x̃(k) + D′ ũ(k)
where x̃ ∈ Rn , ũ ∈ Rp and ỹ ∈ Rm
• The reachability [observability] matrix of the dual system is equal to the
transpose of the observability [reachability] matrix of the original system
h i
R̃ = C′ A′ C ′ ... (A′ )n−1 C ′ = Θ′
 
B′
 
 B ′ A′ 
Θ̃ =   = R′
 .. 
 . 
′ ′ n−1
B (A )
• The system (A, B, C, D) is reachable [observable] if and only if its dual system
(A′ , C ′ , B ′ , D ′ ) is observable [reachable]
©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 114
State estimation
State estimation problem
At each time k construct an estimate x̂(k) of the state x(k), by only measuring
the output y(k) and input u(k).

• Open-loop observer: Build an artificial copy of the system, fed in parallel by

with the same input signal u(k)
!"#$%&'$()*+,'-..

u(k) x(k) y(k)

A,B /+0-)./$/-
C

A,B ˆ
x(k)
./$/-)-./&%$/-

• The “copy” is a numerical simulator x̂(k + 1) = Ax̂(k) + Bu(k) reproducing

the behavior of the real system

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 115

Open-loop observer
!"#$%&'$()*+,'-..

u(k) x(k) y(k)

A,B /+0-)./$/-
C

A,B ˆ
x(k)
./$/-)-./&%$/-

• The dynamics of the real system and of the numerical copy are

x(k + 1) = Ax(k) + Bu(k) True process

x̂(k + 1) = Ax̂(k) + Bu(k) Numerical copy

• The dynamics of the estimation error x̃(k) = x(k) − x̂(k) are

x̃(k + 1) = Ax(k) + Bu(k) − Ax̂(k) − Bu(k) = Ax̃(k)

and then x̃(k) = Ak (x(0) − x̂(0))

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 116

Open-loop observer
!"#$%&'$()*+,'-..

u(k) x(k) y(k)

A,B /+0-)./$/-
C

A,B ˆ
x(k)
./$/-)-./&%$/-

The estimation error is x̃(k) = Ak (x(0) − x̂(0)). This is not ideal, because

• The dynamics of the estimation error are fixed by the eigenvalues of A and
cannot be modified
• The estimation error vanishes asymptotically if and only if A is asymptotically
stable
• Note that we are not exploiting y(k) to compute the state estimate x̂(k) !

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 117

Luenberger observer
!"#$%&'$()*+,'-..

u(k) x(k) y(k)

A,B /+0-)./$/-
C
+
ˆ
x(k) ˆ -
y(k)
A,[B L] ./$/-)
C
-./&%$/-

./$/-),1.-+2-+

• Luenberger observer: Correct the estimation equation with a

feedback from the estimation error y(k) − ŷ(k)

x̂(k + 1) = Ax̂(k) + Bu(k) + L(y(k) − C x̂(k))

| {z }
feedback on estimation error David G. Luenberger
(1937–)

where L ∈ Rn×p is the observer gain

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 118
Luenberger observer
!"#$%&'$()*+,'-..

u(k) x(k) y(k)

A,B /+0-)./$/-
C
+
ˆ
x(k) ˆ -
y(k)
A,[B L] ./$/-)
C
-./&%$/-

./$/-),1.-+2-+

• The dynamics of the state estimation error x̃(k) = x(k) − x̂(k) is

x̃(k + 1) = Ax(k) + Bu(k) − Ax̂(k) − Bu(k) − L[y(k) − C x̂(k)]
= (A − LC)x̃(k)
and then x̃(k) = (A − LC)k (x(0) − x̂(0))
• Same idea for continuous-time systems ẋ(t) = Ax(t) + Bu(t)
dx̂(t)
= Ax̂(t) + Bu(t) + L[y(t) − C x̂(t) − Du(t)]
dt
dx̃(t)
The dynamics of the state estimation error are dt = (A − LC)x̃(t)
©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 119
Eigenvalue assignment of state observer
Theorem
If the pair (A, C) is “observable” (= (A′ , C ′ ) “reachable”), then the eigenvalues
of (A − LC) can be placed arbitrarily.

MATLAB Python
where P = [λ1 λ2 . . . λn ] =
L=acker(A’,C’,P)’; L=[Link](A.T,C.T,p).T
L=place(A’,C’,P)’; L=[Link](A.T,C.T,p).T desired observer eigenvalues

80
true state
estimator L1
60 estimator L2
estimator L3 response from initial conditions
 
40
x(0) = −1 0
1 , x̂(0) = [ 0 ] for
20 u(k) ≡ 0.1 for different choices of
the observer poles
0

−20
0 10 20 30 40
time (s)
©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 120
Dynamic compensators
Potential issues in state feedback control

• Measuring the entire state vector may be too expensive (many sensors)

• It may be even impossible (high temperature, high pressure, inaccessible

environment)

Can we use the estimate x̂(k) instead of x(k) to close the loop ?

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 121

Dynamic compensator
!"#$%&'$()*+,'-..

u(k) x(k) y(k)

A,B C

v(k) + + ˆ
x(k) .0$0-
K -.0&%$0,+

!"#$%&'),/0*/0)1--!2$'3)',#0+,((-+

• Assume the open-loop system is completely observable and reachable

• Construct the linear state observer

x̂(k + 1) = Ax̂(k) + Bu(k) + L(y(k) − C x̂(k))

• Set u(k) = K x̂(k) + v(k)
• The dynamics of the error estimate x̃(k) = x(k) − x̂(k) is

x̃(k+1) = Ax(k)+Bu(k)−Ax̂(k)−Bu(k)+L(Cx(k)−C x̂(k)) = (A−LC)x̃(k)

The error estimate does not depend on the feedback gain K !

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 122

Closed-loop dynamics
• Let’s combine the dynamics of the system, observer, and feedback gain


 x(k + 1) = Ax(k) + Bu(k)

 x̂(k + 1) = Ax̂(k) + Bu(k) + L(y(k) − C x̂(k))

 u(k) = K x̂(k) + v(k)


y(k) = Cx(k)
• Take x(k), x̃(k) as state components of the closed-loop system
" # " #" #
x(k) I 0 x(k)
= (it is indeed a change of coordinates)
x̃(k) I −I x̂(k)
• The closed-loop dynamics is
 " # " #" # " #

 x(k + 1) A + BK −BK x(k) B

 = + v(k)

 A − LC
 x̃(k + 1) 0 x̃(k) 0


 " #

 h i x(k)

 y(k) =
 C 0
x̃(k)
©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 123
Closed-loop dynamics
• The transfer function from v(k) to y(k) is
" #−1 " #
h i zI − A − BK BK B
G(z) = C 0
0 zI − A + LC 0
" #" #
h i (zI − A − BK)−1 ⋆ B
= C 0
0 (zI − A + LC)−1 0
N (z)
= C(zI − A − BK)−1 B =
DK (z)
• Even if we substituted x(k) with x̂(k), the input-output behavior of the
closed-loop system didn’t change !
The closed-loop poles can be assigned arbitrarily using dynamic output feed-
back, as in the state feedback case

The closed-loop transfer function does not depend on the observer gain L

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 124

Separation principle
Separation principle
The design of the control gain K and of the observer gain L can be done inde-
pendently

• Watch out ! G(z) = C(zI − A − BK)−1 B only represents the I/O

(=input/output) behavior of the closed-loop system

• The complete set of poles of the closed-loop system are given by

 A+BK −BK

det(zI− 0 A−LC ) = det(zI−A−BK) det(zI−A+LC) = DK (z)DL (z)

• A zero/pole cancellation of the observer poles has occurred:

" #
h i  A+BK  B N (z)DL (z)
G(z) = C 0 (zI − 0
−BK
A−LC )−1 =
0 DK (z)DL (z)

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 125

Transient effects of the estimator gain

• L has an effect on the natural response of the system !

h i
x(0)
• To see this, consider the effect of a nonzero initial condition x̃(0)
for v(k) ≡ 0

y(0) = Cx(0)
h i  h x(0) i
−BK
y(1) = C 0 A+BK 0 A−LC x̃(0)
h ih i
(A+BK)x(0)−BK x̃(0)
= C 0 (A−LC)x̃(0)
= C(A + BK)x(0) − CBK x̃(0)
h i  h i
−BK x(1)
y(2) = C 0 A+BK 0 A−LC x̃(1)

= C(A + BK)x(1) − CBK x̃(1)

= C(A + BK)2 x(0) − C(A + BK)BK x̃(0) − CBK(A − LC)x̃(0)

• If x̃(0) 6= 0, L has an effect during the transient !

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 126

Choosing the estimator gain

• Intuitively, if x̂(k) is a poor estimate of x(k) then the control action will also be
poor

Rule of thumb: place the observer poles ≈ 10 times faster than

the controller poles

• Optimal methods exist to choose the observer poles (Kalman filter)

• Fact: The choice of L is very important for determining the sensitivity of the
closed-loop system with respect to input and output noise

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 127

Example: Control of a DC Motor
d3 y d2 y dy u y
+β +α = Ku
dt dt dt

MATLAB MATLAB
K=1; beta=.3; alpha=1; % Closed-loop system, state=[x;xhat]
G=tf(K,[1 beta alpha 0]);
bigA=[A,B*K;L*C,A+B*K-L*C];
ts=0.5; % sampling time bigB=[B;B];
Gd=c2d(G,ts); bigC=[C,zeros(1,3)];
sysd=ss(Gd); bigD=0;
[A,B,C,D]=ssdata(sysd); clsys=ss(bigA,bigB,bigC,bigD,ts);

% Controller x0=[1 1 1]’; % Initial state

polesK=[-1,-0.5+0.6*j,-0.5-0.6*j]; xhat0=[0 0 0]’; % Initial estimate
polesKd=exp(ts*polesK); T=20;
K=-place(A,B,polesKd); initial(clsys, [x0;xhat0],T);
pause
% Observer
polesL=[-10, -9, -8]; t=(0:ts:T)’;
polesLd=exp(ts*polesL); v=ones(size(t));
L=place(A’,C’,polesLd)’; lsim(clsys,v);

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 128

Example: Control of a DC Motor
d3 y d2 y dy u y
+β +α = Ku
dt dt dt
Python
import numpy as np
Python
K,beta,alpha = 1,.3,1 # Closed-loop system, state=[x;xhat]
G=[Link](K,[1,beta,alpha,0])
bigA=[Link](([Link]((A,B@K)),
ts=0.5 # sampling time [Link]((L@C,A+B@K-L@C))))
Gd=ctrl.c2d(G,ts) bigB=[Link]((B,B))
sysd=[Link](Gd) bigC=[Link]((C,[Link]((1,3))))
A,B,C,D=[Link](sysd) bigD=0
clsys=[Link](bigA,bigB,bigC,bigD,ts)
# Controller
polesK=[Link]([-1,-0.5+0.6j,-0.5-0.6j]) x0=[Link]([[1],[1],[1]]) # Initial state
polesKd=[Link](ts*polesK) xhat0=[Link]((3,1)) # Initial estimate
K=-[Link](A,B,polesKd) t=[Link](0,20+ts,ts)
_,y=ctrl.initial_response(clsys,t,[Link]((x0,xhat0)))
# Observer
polesL=[Link]([-10, -9, -8]) v=[Link]([Link])
polesLd=[Link](ts*polesL) _,yf,xf=ctrl.forced_response(clsys,t,v)
L=[Link](A.T,C.T,polesLd).T

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 129

Example: Control of a DC Motor

y(k) y(k)
2 3

2.5
1.5

2
1
1.5
0.5
1

0
0.5

−0.5 0
0 5 10 15 20 0 5 10 15 20
time (s) time (s)

h1i h0i h0i

x(0) = 1 , x̂(0) = 0 , v(k) ≡ 0 x(0) = x̂(0) = 0 , v(k) ≡ 1
1 0 0

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 130

Linear quadratic regulation
Linear quadratic regulation (LQR)
• State-feedback control via pole placement requires one to assign the
closed-loop poles
• Any way to place closed-loop poles automatically and optimally ?
• The main control objectives are
1. Make the state x(k) “small” (to converge to the origin)
2. Use “small” input signals u(k) (to minimize actuators’ effort)
These are conflicting goals !

• LQR is a technique to place automatically and optimally the closed-loop poles

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 131

Finite-time optimal control

• Linear system x(k + 1) = Ax(k) + Bu(k) with initial condition x(0)

• We look for the optimal sequence of inputs

U = {u(0), u(1), . . . , u(N − 1)}

driving x(k) towards the origin while minimizing the performance index

X
N −1
′
J(x(0), U ) = x (N )QN x(N ) + x′ (k)Qx(k) + u′ (k)Ru(k) quadratic cost
k=0

where Q = Q′  0, R = R′  0, QN = Q′N  06

6 For a matrix Q ∈ Rn×n , Q ≻ 0 means that Q is a positive definite matrix, i.e., x′ Qx > 0 for all

x ̸= 0, x ∈ Rn . Q ⪰ 0 means positive semidefinite, x′ Qx ≥ 0, ∀x ∈ Rn .

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 132

Finite-time optimal control
• Example: Q diagonal Q = diag(q1 , . . . , qn ), single input, QN = 0
!
X
N −1 X
n
J(x(0), U ) = qi x2i (k) + Ru2 (k)
k=0 i=1

• Consider again the general linear quadratic (LQ) problem

X
N −1
J(x(0), U ) = x′ (N )QN x(N ) + x′ (k)Qx(k) + u′ (k)Ru(k)
k=0

– N is called the time horizon over which we optimize performance

– The first term x′ Qx penalizes the deviation of x from the desired target x = 0
– The second term u′ Ru penalizes actuator authority
– The third term x′ (N )QN x(N ) penalizes how much the final state x(N ) deviates
from the target x = 0
• Q, R, QN are the tuning parameters of optimal control design (cf. the
parameters of the PID controller Kp , Ti , Td )

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 133

Minimum-energy controllability
• Consider again the problem of controllability of the state to zero with minimum
energy input
  2
u(0)
 
 u(1) 
minU  
 .. 
 . 
u(N − 1) 2
s.t. x(N ) = 0

• The minimum-energy control problem can be seen as a particular case of the LQ

optimal control problem by setting
R = I, Q = 0, QN = ∞ · I

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 134

Solution to LQ optimal control problem
Pk−1
• By substituting x(k) = Ak x(0) + i=0 Ai Bu(k − 1 − i) in
X
N −1
J(x(0), U ) = x′ (k)Qx(k) + u′ (k)Ru(k) + x′ (N )QN x(N )
k=0

we obtain
1 ′ 1
J(x(0), U ) = U HU + x(0)′ F U + x(0)′ Y x(0)
2 2
where H = H ′  0 is a positive definite matrix
• The optimizer U ∗ is obtained by zeroing the gradient

0 = ∇U J(x(0), U ) = HU + F ′ x(0)
 
u∗ (0)
 
 u∗ (1) 
−→ U ∗ =   = −H −1 F ′ x(0)
 .. 
 . 
∗
u (N − 1)

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 135

[LQ problem matrix computation]
 ′   
x(1) Q 0 0 ... 0 x(1)
 x(2)   0 Q 0 ... 0  x(2) 
    
′  .   . . . .  . 
J(x(0), U ) = x (0)Qx(0) + 
 . 


 . .
.
. . . 
 . +

 .   . . . . .  . 
 x(N − 1)   0 ... 0 Q 0
  x(N − 1) 
x(N ) 0 0 ... 0 QN x(N )
| {z }
Q̄
  
R 0 ... 0 u(0)
 0 R ... 0  u(1) 
   
u′ (0) u′ (1) ... u′ (N − 1) 
 . . .   . 
.   
.
 . . .  . 
. . . . .
0 ... 0 R u(N − 1)
| {z }
R̄
S̄
z }| {  
   B 0 ... 0  A
x(1)  
 x(2)   AB B ... 0  u(0)  A2 
     
   . 
u(1)   
 .  =  . . .   +  .  x(0)
 .   . . . .  ...  . 
. . . . .  u(N − 1)  . 
x(N ) AN − 1 B A N −2 B ... B AN
| {z }
N̄
′ ′ ′
J(x(0), U ) = x (0)Qx(0) + (S̄U + N̄ x(0)) Q̄(S̄U + N̄ x(0)) + U R̄U
1 ′ ′ ′ ′ 1 ′ ′
= U 2(R̄ + S̄ Q̄S̄) U + x (0) 2N̄ Q̄S̄ U + x (0) 2(Q + N̄ Q̄N̄ ) x(0)
2 | {z } | {z } 2 | {z }
H F Y

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 136

Solution to LQ optimal control problem

• The solution  
u∗ (0)
 
 u∗ (1) 
U∗ = 
 ..  = −H −1 F ′ x(0)

 . 
u∗ (N − 1)
is an open-loop one: u∗ (k) = fk (x(0)), k = 0, 1, . . . , N − 1

• Moreover the dimensions of the H and F matrices is proportional to the time

horizon N

• We use optimality principles next to find a better solution (computationally

more efficient, and more elegant)

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 137

Dynamic programming
• Consider the following basic fact in optimization

V0 ≜ min f (z, y) = min{ min f (z, y) }

z,y z y
| {z }
this is a function of z

• In case f is separable in the sum of two functions

f (z, y) ≜ f0 (z) + f1 (z, y)

we get miny f (z, y) = f0 (z) + miny f1 (z, y)

• Therefore we can compute V0 in two steps:

V1 (z) = min f1 (z, y)

y
V0 = min{f0 (z) + V1 (z)}
z

• We apply the above reasoning to f = J(x(0), U ), z = [u′ (0) . . . u(k1 − 1)′ ]′ ,

y = [u′ (k1 ) . . . u(N − 1)′ ]′

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 138

Dynamic programming
• At a generic instant k1 and state x(k1 ) = x̄(z) consider the optimal cost-to-go
 
NX−1 
′ ′ ′
Vk1 (x̄(z)) = min x (k)Qx(k) + u (k)Ru(k) + x (N )QN x(N )
u(k1 ),...,u(N −1)  
k=k1

Principle of dynamic programming

V0 (x(0)) = min J(x(0), U )

U ≜{u(0),...,u(N −1)}
(k −1 )
X1
′ ′
= min x (k)Qx(k) + u (k)Ru(k) + Vk1 (x(k1 ))
u(0),...,u(k1 −1)
k=0

• Solving over [0, k1 ] with terminal weight equal to the optimal cost-to-go from k1
to N at x(k1 ) is the same as solving over [0, T ]

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 139

 Bellman's principle of optimality
Bellman’s principle
Given the optimal sequence U ∗ = [u∗ (0), . . . , u∗ (N − 1)]
(and the corresponding optimal trajectory x∗ (k)), the subsequence
[u∗ (k1 ), . . . , u∗ (N − 1)] is optimal for the problem on the horizon
[k1 , N ], starting from the optimal state x∗ (k1 ) Richard Bellman
(1920-1984)

optimal state x∗ (k)

• Given the state x∗ (k1 ), the optimal input trajectory u∗

on the remaining interval [k1 , N ] only depends on
time
0 k1 N
x∗ (k1 )
optimal input u (k)
∗
• Then each optimal move u∗ (k) of the optimal trajectory
on [0, N ] only depends on x∗ (k)
• The optimal control policy can be always expressed in
time state feedback form u∗ (k) = u∗k (x∗ (k)) !
0 k1 N
day, May 11, 2010

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 140

Bellman's principle of optimality

• The principle also applies to nonlinear systems

optimal state trajectories x∗
and/or non-quadratic cost functions: the optimal
control law can be always written in state-feedback
form

u∗ (k) = fk (x∗ (k)), ∀k = 0, . . . , N − 1

• Compared to the open-loop solution {u∗ (0), . . . , u∗ (N − 1)} = f (x(0)) the

feedback form u∗ (k) = fk (x∗ (k)) has the big advantage of being more robust
with respect to perturbations: at each time k we apply the best move on the
remaining period [k, N ]

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 141

Riccati iterations
By applying the dynamic programming principle, we can compute the optimal
inputs u∗ (k) recursively as a function of x∗ (k) (Riccati iterations):

1. Initialization: P (N ) = QN
2. For k = N, . . . , 1, compute recursively the following
matrix

P (k−1) = Q−A′ P (k)B(R+B ′ P (k)B)−1 B ′ P (k)A+A′ P (k)A

3. Define

K(k) = −(R + B ′ P (k + 1)B)−1 B ′ P (k + 1)A

The optimal input is

u∗ (k) = K(k)x∗ (k)

Jacopo Francesco Riccati
(1676–1754)

The optimal input policy u∗ (k) is a (linear time-varying) state feedback !

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 142
Linear quadratic regulation
• Consider the infinite-horizon optimal control problem
∞
X
∞
V (x(0)) = min x′ (k)Qx(k) + u′ (k)Ru(k)
u(0),u(1),...
k=0

Result
Let (A, B) be a stabilizable pair, R  0, Q  0. There exists a unique solution
P∞ of the algebraic Riccati equation (ARE)

P∞ = A′ P∞ A + Q − A′ P∞ B(B ′ P∞ B + R)−1 B ′ P∞ A

such that the optimal cost is V ∞ (x(0)) = x′ (0)P∞ x(0) and the optimal con-
trol law is the constant linear state feedback u(k) = KLQR x(k) with

KLQR = −(R + B ′ P∞ B)−1 B ′ P∞ A.

MATLAB
Python E = closed-loop poles
P∞ = dare(A,B,Q,R)
P∞ ,E,Km=[Link](A,B,Q,R)
[-K∞ ,P∞ ] = dlqr(A,B,Q,R) = eigenvalues of (A + BKLQR )
K∞ =-Km
[-K∞ ,P∞ ,E] = lqr(sysd,Q,R)
©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 143
Linear quadratic regulation
• Go back to Riccati iterations: starting from P (∞) = P∞ and going backwards
we get P (j) = P∞ , ∀j ≥ 0

• Accordingly, we get

K(j) = −(R + B ′ P∞ B)−1 B ′ P∞ A ≜ KLQR , ∀j = 0, 1, . . .

• The LQR control law is linear and time-invariant

• (A, B) stabilizable implies closed-loop asymptotic stability, ∀R, Q  07

• LQR is an automatic and optimal way of placing poles !

• A similar result holds for continuous-time linear systems

7 R ≻ 0, Q ⪰ 0 with Q = F ′ F , F ∈ Rn×nq , n = rank Q, and (A, F ) detectable also ensures
q
closed-loop asymptotic stability. Matrix F can be obtained for example by the LDLT
[D ] 1
decomposition Q = [L1 L2 ] 1 0
0 0
[L1 L2 ]′ , F = D12 L′1 , with D1 ∈ Rnq ×nq .

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 144

LQR with output weighting
• We often want to regulate only y(k) = Cx(k) to zero, so define
∞
X
V ∞ (x(0)) = min y ′ (k)Qy y(k) + u′ (k)Ru(k)
u(0),u(1),...
k=0

• The problem is again an LQR problem with equivalent state weight Q = C ′ Qy C

MATLAB Python
[-K∞ ,P∞ ,E] = dlqry(sysd,Qy,R) P∞ ,E,-K∞ = [Link](A,B,C.T@Qy@C),R)
Corollary
Let (A, B) stabilizable, (A, C) detectable, R > 0, Qy > 0. The LQR control
law u(k) = KLQR x(k) the asymptotically stabilizes the closed-loop system

lim x(t) = 0, lim u(t) = 0

t→∞ t→∞

Intuitively: the minimum cost x′ (0)P∞ x(0) is finite ⇒ y(k) → 0 and u(k) → 0.
y(k) → 0 implies that the observable part of the state → 0. As u(k) → 0, the unobservable states
remain undriven and go to zero spontaneously (=detectability condition)

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 145

LQR example
• Two-dimensional single input single output (SISO) dynamical system (double
integrator)
" # " #
1 1 0
x(k + 1) = x(k) + u(k)
0 1 1
h i
y(k) = 1 0 x(k)

• LQR (infinite horizon) controller defined on the performance index

∞
X 1 2
V ∞ (x(0)) = min y (k) + u2 (k), ρ > 0
u(0),u(1),... ρ
k=0
h1 i
• Weights: Q = [ 10 ] · 1
ρ · [1 0] = ρ 0 ,R = 1
0 0
• Note that only the ratio Q11 /R = ρ1
matters, as scaling the cost function does
not change the optimal control law

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 146

LQR Example
output y(k)
ρ = 0.1 (red line)
1

K = [−0.8166 − 1.7499]
0.5

−0.5
0 2 4 6 8 10 12 14 16 18 20
ρ = 10 (blue line)
input u(k) K = [−0.2114 − 0.7645]
1

0.5

−0.5
ρ = 1000 (green line)
K = [−0.0279 − 0.2505]
−1
0 2 4 6 8 10 12 14 16 18 20

Initial state: x(0) = [ 10 ]

∞
X 1 2
V ∞ (x(0)) = min y (k) + u2 (k)
u(0),u(1),... ρ
k=0

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 147

Kalman filtering
Kalman filtering - Introduction
• Problem: assign observer poles in an optimal way, that is to minimize the state
estimation error x̃ = x − x̂
• Information comes in two ways: from sensors measurements (a posteriori) and
from the model of the system (a priori)
• We need to mix the two information sources optimally, given a probabilistic
description of their reliability (sensor precision, model accuracy)

The Kalman filter solves this problem, and is now

the most used state observer in most engineering
fields (and beyond)

Rudolf E. Kalman∗
(1930–2016)
∗ R.E. Kalman receiving the Medal of Science from the President of the USA on October 7, 2009

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 148

Process model
• The process is modeled as the linear time-varying system with noise
x(k + 1) = A(k)x(k) + B(k)u(k) + G(k)ξ(k)
y(k) = C(k)x(k) + D(k)u(k) + ζ(k) x(0) = x0

• ξ(k) ∈ Rq = process noise, E[ξ(k)] = 0 (zero mean), E[ξ(k)ξ ′ (j)] = 0, ∀k 6= j ,

(white noise), E[ξ(k)ξ ′ (k)] = Q(k)  0 (covariance matrix)

• ζ(k) ∈ Rp = measurement noise, E[ζ(k)] = 0, E[ζ(k)ζ ′ (j)] = 0 ∀k 6= j ,

E[ζ(k)ζ ′ (k)] = R(k)  0
• x0 ∈ Rn is a random vector, E[x0 ] = x̄0 , P0 = E[(x0 − x̄0 )(x0 − x̄0 )′ ], P0  0

• Vectors ξ(k), ζ(k), x0 are uncorrelated: E[ξ(k)ζi (j)] = 0, E[ξ(k)x′0 ] = 0,

E[ζi (k)x′0 ] = 0, ∀k, j ∈ Z, ∀i = 1, . . . , p

• Probability distributions: we often assume normal (=Gaussian) distributions

ξ(k) ∼ N (0, Q(k)), ζ(k) ∼ N (0, R(k)), x0 ∼ N (x̄0 , P0 )

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 149

Kalman filter
In formulating the Kalman filter we use the following notation:

x̂(k|k − 1) state estimate at time k based on

data up to time k − 1
x̂(0| − 1) = x̄0 initial state estimate
x̃(k|k − 1) = x(k) − x̂(k|k − 1) state estimation error
′
P (k|k − 1) = E [x̃(k|k − 1)x̃(k|k − 1) ] covariance of state estimation error
x̂(k|k) state estimate at time k
based on data up to time k
x̃(k|k) = x(k) − x̂(k|k) state estimation error
′
P (k|k) = E [x̃(k|k)x̃(k|k) ] covariance of state estimation error
x̂(k + 1|k) state prediction at time k + 1
based on data up to time k

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 150

Kalman filter
• The Kalman filter provides the optimal estimate x̂(k|k) of x(k) given the
measurements up to time k

• Optimality means that the trace of the variance P (k + 1|k) is minimized

• The filter is based on two steps:

1. measurement update based on the most recent y(k)

M (k) = P (k|k − 1)C(k)′ [C(k)P (k|k − 1)C(k)′ + R(k)]−1

x̂(k|k) = x̂(k|k − 1) + M (k) (y(k) − C(k)x̂(k|k − 1) − D(k)u(k))
P (k|k) = (I − M (k)C(k))P (k|k − 1)
with initial conditions x̂(0| − 1) = x̄0 , P (0| − 1) = P0

2. time update based on the model of the system

x̂(k + 1|k) = A(k)x̂(k|k) + B(k)u(k)

P (k + 1|k) = A(k)P (k|k)A(k)′ + G(k)Q(k)G(k)′

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 151

Stationary Kalman filter
• Assume A, C , G, Q, R are constant (time-invariant case)
• Under suitable assumptions8 , P (k|k − 1), M (k) converge to constant matrices
−1
P∞ = AP∞ A′ + GQG′ − AP∞ C ′ [CP∞ C ′ + R] CP∞ A′
M∞ = P∞ C ′ (CP∞ C ′ + R)−1

• By setting L∞ = AM∞ the dynamics of the prediction x̂(k|k − 1) becomes the

Luenberger observer

x̂(k + 1|k) = Ax̂(k|k − 1) + B(k)u(k) + L∞ (y(k) − C x̂(k|k − 1) − D(k)u(k))

with all the eigenvalues of (A − L∞ C) inside the unit circle

Python
MATLAB
P∞ ,_,Lt=[Link](A.T,C.T,G.T@Q@G,R)
[~,L,P∞ ,M]=kalman(sys,Q,R)
L=Lt.T; M=[Link](A,L)[0]
8 (A, C) observable, and (A, GB ) stabilizable, where B is such that Q = B B ′ , cf. conditions
q q q q
for asymptotic stability of LQR

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 152

Tuning Kalman filters
• It is usually hard to quantify exactly the correct values of Q and R for a given
process

• The diagonal terms of R are related to how noisy are output sensors

• Q is harder to relate to physical noise, it mainly relates to how rough is the

(A, B) model

• After all, Q and R are the tuning knobs of the observer (similar to LQR)

• The “larger” is R with respect to Q the “slower” is the observer to converge (L,
M will be small)

• On the contrary, the “smaller” is R than Q, the more precise are considered the
measurments, and the “faster” observer will be to converge

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 153

Extended Kalman filter
• The Kalman filter can be extended to nonlinear systems
x(k + 1) = f (x(k), u(k), ξ(k))
y(k) = g(x(k), u(k)) + ζ(k)
1. Measurement update:
∂g
C(k) = (x̂(k|k − 1), u(k))
∂x
M (k) = P (k|k − 1)C(k)′ [C(k)P (k|k − 1)C(k)′ + R(k)]−1
x̂(k|k) = x̂(k|k − 1) + M (k) (y(k) − g(x̂(k|k − 1), u(k)))
P (k|k) = (I − M (k)C(k))P (k|k − 1)
2. Time update:

x̂(k + 1|k) = f (x̂(k|k), u(k), E[ξ(k)]), x̂(0| − 1) = x̂0

∂f ∂f
A(k) = (x̂(k|k), u(k), E[ξ(k)]), G(k) = (x̂(k|k), u(k), E[ξ(k)])
∂x ∂ξ
P (k + 1|k) = A(k)P (k|k)A(k)′ + G(k)Q(k)G(k)′ , P (0| − 1) = P0

• The EKF is in general not optimal and may even diverge, due to linearization.
But is the de-facto standard in nonlinear state estimation
©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 154
LQG control
• Linear Quadratic Gaussian (LQG) control combines an LQR control law and a
stationary Kalman predictor/filter
• Consider the stochastic dynamical system

x(k + 1) = Ax(k) + Bu(k) + ξ(k), ξ(k) ∼ N (0, QKF )

y(k) = Cx(k) + ζ(k), ζ(k) ∼ N (0, RKF )

with initial condition x(0) = x0 , x0 ∼ N (x̄0 , P0 ), P0 , QKF  0, RKF  0, and

ζ and ξ are independent and white noise terms.

• The objective is to minimize the cost function

" T #
1 X
J(x(0), U ) = lim E x′ (k)QLQ x(k) + u′ (k)RLQ u(k)
T →∞ T
k=0

when the state x is not measurable

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 155

LQG control
If we assume that all the assumptions for LQR control and Kalman
predictor/filter hold, i.e.

• the pair (A, B) is reachable and the pair (A, Cq ) with Cq such that
QLQ = Cq Cq′ is observable (here Q is the weight matrix of the LQ controller)
• the pair (A, Bq ), with Bq s.t. QKF = Bq Bq′ , is stabilizable, and the pair (A, C) is
observable (here Q is the covariance matrix of the Kalman predictor/filter)

Then, apply the following procedure:

1. Determine the optimal stationary Kalman predictor/filter, neglecting the fact

that the control variable u is generated through a closed-loop control scheme,
and find the optimal gain LKF

2. Determine the optimal LQR strategy assuming the state accessible, and find the
optimal gain KLQR

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 156

 LQG control
!"#$%&'$()*+,'-..

u(k) x(k) y(k)

A,B C

v(k) + + ˆ
x(k) 3$(%$#
KLQR 4(2-+

/01)',#2+,((-+

Analogously to the case of output feedback control using a Luenberger

observer, it is possible to show that the extended state [x′ x̃′ ]′ has eigenvalues
equal to the eigenvalues of (A + BKLQR ) plus those of (A − LKF C) (2n in
total)

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems''

Sunday, May 16, 2010
157
System identification
Model identification

• Designing a control system requires a dynamical model of the process

• Often a dynamical model can be difficult to obtain due to the complexity of the
process, whose dynamics may be even (partially or completely) unknown

• Even if we have a mathematical model, sometimes this is too complex to base a

controller design on it (large state dimensions, nonlinearities, etc.)

System identification is a procedure to build a mathematical model of

the dynamics of a system from measured data

Lecture based on
[1] L. Ljung, “System Identification,” Control Systems Handbook (W. Levine ed.), CRC Press, pp. 1033–1054, 1995
[2] L. Ljung, “System Identification: Theory for the User,” Prentice Hall, 1987

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 158

Model identification
input output
(measured)
u(t) ? (measured)
y(t)
estimated
dynamical process
identification model
algorithm
^
G(z)

Different types of identification:

• White box: model structure based on first principles (e.g., Newton’s law), model
parameters estimated from measured data

• Grey box: model structure partially known from first principles, the rest is
reconstructed from data

• Black box: model structure and its parameters completely unknown, they are
only estimated from I/O data

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 159

Step-response identification

u(t) ? y(t)
dynamical process
identification ^
G(z)
algorithm

• Excite the process with a step u(t) = 1I(t), record output response y(t)

• Observe the shape of y(t) and reconstruct G(z)

(1st-order response ? 2nd-order undamped response ? Any delay ? ...)

• Mostly used in process control: excitation experiment is easily done,

superposition of effects can be used in the multivariable case to identify each
entry Gij (z) of the transfer matrix G(z), one at the time

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 160

Black-box identification via linear regression
• Consider the black-box ARX (AutoRegressive eXogenous) model

y(k)+a1 y(k−1)+...+ana y(k−na ) = b1 u(k−nk )+...+bnb u(k−nk −nb +1)+e(k)

where e(k) is zero-mean white noise and y(k), u(k), e(k) ∈ R

• We can predict the next output value given previous observations

y(k) = −a1 y(k−1)−...−ana y(k−na )+b1 u(k−nk )+...+bnb u(k−nk −nb +1)+e(k)

• In more compact form

y(k) = φ′ (k)θ + e(k)
h i′
θ= a1 . . . a na b1 ... b nb unknown parameter vector
h i′
φ(k) = −y(k − 1) . . . − y(k − na ) u(k − nk ) ... u(k − nk − nb + 1)
regressor

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 161

Black-box identification via linear regression
• Let ŷ(k|θ) = φ′ (k)θ = output prediction, which is an estimate of y(k) based on
the parameter vector θ and past data (e(k) = 0 is the best we can assume)

• We don’t know θ, but we have collected a set Z N of measured data

Z N = {u(−n), y(−n), ..., u(N − 1), y(N − 1)} , n = max{na , nb + nk − 1}

• We solve a least-squares problem to estimate the vector θ∗ that best makes
ŷ(k|θ) fit y(k)

θ∗ = arg min V (θ, Z N )
θ

with

1 X
N −1
V (θ, Z N ) = (y(k) − ŷ(k|θ))2

y(k)
N
k=0
1 X
N −1

2
= y(k) − φ′ (k)θ
N
k=0 ϕ(k)

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 162

Black-box identification via linear regression
• V (θ, Z N ) is a quadratic function of θ. We find the minimum by zeroing the
derivative of V
N −1
d 2 X
0= VN (θ, Z N ) = − φ(k) (y(k) − φ′ (k)θ)
dθ N
k=0

or
X
N −1 X
N
φ(k)y(k) = φ(k)φ′ (k)θ
k=0 k=1

• The best parameter vector we can choose is therefore

"N −1 #−1 N −1
X X
θ∗ = φ(k)φ′ (k) φ(k)y(k)
k=0 k=0

MATLAB
θ ∗ =arx(Z N ,[na nb nk ])

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 163

Recursive linear regression

• Drawback of (batch) linear regression: if a new data pair u(N ), y(N ) is acquired
hP i−1
′ N
the new matrix k=0 φ(k)φ (k) is required to compute the new optimal
∗
parameter vector θ

• Computations become more and more expensive as N keeps growing

• Given the best estimate θ∗ (k − 1) obtained using k − 1 data points

 −1
X
k−1 X
k−1
θ∗ (k − 1) = P (k − 1) φ(j)y(j), P (k − 1) =  φ(j)φ′ (j)
j=0 j=0

we would like to get θ∗ (k) without solving the regression problem from scratch

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 164

Recursive linear regression
• Since
X
k−1
P −1 (k) = P −1 (k − 1) + φ(k)φ′ (k), P −1 (k − 1)θ∗ (k − 1) = φ(j)y(j)
j=0

we get
 
X
k X
k−1
θ∗ (k) = P (k) φ(j)y(j) = P (k)  φ(j)y(j) + φ(k)y(k)
j=0 j=0


= P (k) P −1 (k − 1)θ∗ (k − 1) + φ(k)y(k)


= P (k) (P −1 (k) − φ(k)φ′ (k))θ∗ (k − 1) + φ(k)y(k)

and therefore

θ∗ (k) = θ∗ (k − 1) + P (k)φ(k)(y(k) − θ∗ (k − 1)′ φ(k))

| {z }
estimation error

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 165

Recursive linear regression
• Since
 −1
X
k
 −1
P (k) =  φ(j)φ′ (j) = P −1 (k − 1) + φ(k)φ′ (k)
j=0

we can apply the Matrix Inversion Lemma to update P (k) recursively and get

P (k − 1)φ(k)φ′ (k)P (k − 1)
P (k) = P (k − 1) −
1 + φ′ (k)P (k − 1)φ(k)

X
k−1
• Let m = dim θ = na + nb . Matrix φ(j)φ′ (j) is not invertible for k < m
j=0
"m−1 #−1
X ′
• We can start the recursions after m steps with P (m − 1) = φ(j)φ (j)
j=0
X
m−1
(if the inverse exists) and θ∗ (m − 1) = P (m − 1)−1 φ(j)y(j)
j=0
©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 166
Recursive linear regression
• Alternatively, we can initialize P (−1) = ρ2 I , ρ > 0, and θ∗ (−1) = θ̄

• Interpretation: we are adding a regularization term on θ

1 Xk
′
2 1 X
m Xk

2
min 2
kθ − θ̄k 2
2 + y(j) − φ (j)θ = 2
(θ i − θ̄ i ) 2
+ y(j) − φ′ (j)θ
ρ j=0
ρ i=1 j=0
m 
X 2 X X
1 1
k

2 k

2
= θ̄i − e′i θ + y(j) − φ′ (j)θ = y(j) − φ′ (j)θ
i=1
ρ ρ j=0 j=−m

where we replaced y(−j) = ρ1 θ̄i , φ(−j) = ρ1 ej , ei = i-th column I , j = 1, . . . , m

 −1  −1
−1
X −1
X 1
• Therefore P (−1) =  φ(j)φ′ (j) = e e ′  = ρ2 I
2 −j −j
j=−m j=−m
ρ
−1
X Xm
1
and θ∗ (−1) = P (−1) φ(j)y(j) = ρ2 I θ̄i ei = θ̄
j=−m i=1
ρ2

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 167

Black-box identification: the general procedure

• Design of experiment: What kind of input excitation u(k) to apply ?

• Model structure: Which class of models do I choose to fit my data ?

• Fit criterion between data and model: How do I best choose the model within
that class (=the parameter vector) ?

• Validation criterion: Is the model that I have identified good enough to

reproduce the dynamics of the process ?

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 168

Design of experiment
• Collecting data is a very crucial (and most expensive) step

• Some theory is available, as well as some practical rules

• The data set Z N should be as informative as possible to fully identify the model9
1

0.8

• Pseudo-random binary signals (PRBS) randomly

0.6

0.4

switching between ±1 are a good choice

0.2

−0.2

MATLAB −0.4

−0.6

» u=idinput(N,’PRBS’); −0.8

−1
0 5 10 15 20 25 30

9 Sinusoidal signals u(t) = sin(ωt) are not good, as only G(jω) would be captured. The input

signal must at least contain as many different frequencies as the order of the chosen structure of
linear models. Step responses are not ideal but ok: |F [1I(t)]| = ω
1
(the Fourier transform of the
continuous-time signal 1I(t)) has in infinite number of frequencies, although decreasing in
amplitude.

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 169

Model structure
• A linear system with additive disturbance v(k) can be described as

y(k) = G(z)u(k) + v(k)

where G(z) is a transfer function (z −m x(k) = x(k − m))

B(z) b0 + b1 z −1 + . . . + bnb z −nb

G(z) = =
A(z) 1 + a1 z −1 + a2 z −2 + . . . + ana z −na

• Alternatively, it can be described by perturbing the difference equation:

y(k) = − a1 y(k − 1) − a2 y(k − 2) − . . . − ana y(k − na )

+ b0 u(k) + b1 u(k − 1) + . . . + bnb u(k − nb ) + w(k)

• The two models are equivalent if we set

w(k) = v(k) + a1 v(k − 1) + a2 v(k − 2) + . . . + ana v(k − na )

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 170

Model structure
• The disturbance v(k) is not necessarily white noise, but can be colored noise

v(k) = H(z)e(k)

where e(k) is white noise and H(z) is another transfer function

C(z) 1 + c1 z −1 + . . . + cnc z −nc

H(z) = =
D(z) 1 + d1 z −1 + d2 z −2 + . . . + dnd z −nd

and H(z) = 0 for z −1 = 0

e(k)
• The overall model is called Box-Jenkins (BJ) model C(z)
D(z)
B(z) C(z)
y(k) = u(k) + e(k) u(k) y(k)
A(z) D(z) B(z) +
+

A(z)

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' Tuesday, June 1, 2010
171
Model structure - Special cases

• Output Error (OE) model: v(k) is white noise (C(z) = D(z) = 1)

B(z)
y(k) = u(k) + e(k)
A(z)

• Auto-Regressive Moving-Average with eXogenous variable (ARMAX) model:

G(z) and H(z) have the same denominator (A(z) = D(z))

A(z)y(k) = B(z)u(k) + C(z)e(k)

• ARX models are a particular case of ARMAX models (C(z) = 1)

A(z)y(k) = B(z)u(k) + e(k)

• ARX and ARMAX models are the most used in practice

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 172

Model structure
• Differently from BJ models, in ARMAX models v(k) and u(k) are filtered by the
1
same dynamics D(z)

• This is justified if the source of disturbance enters early in the process, together
with the input
Example: in airplanes, the disturbances from wind blasts create the same kind
of forces on the airplane as the deflections of the control surfaces

• ARX models are the simplest to compute numerically

e(k)
e(k) C(z)
u(k) y(k) e(k)
B(z) +
+
u(k) +
1 y(k) u(k) +
1 y(k)
+ +

A(z) B(z) B(z)

A(z) A(z)

OE model ARMAX model ARX model

Tuesday, June 1, 2010 Tuesday, June 1, 2010

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 173

Fit criterion

• Let us consider BJ models, which is the most general structure

• Let θ collect all the parameters in the transfer functions G(z) and H(z) to be
estimated from data

y(k) = G(z, θ)u(k) + H(z, θ)e(k)

H −1 (z, θ)y(k) = H −1 (z, θ)G(z, θ)u(k) + e(k)
y(k) + H −1 (z, θ)y(k) = y(k) + H −1 (z, θ)G(z, θ)u(k) + e(k)

• Finally, we get


y(k) = 1 − H −1 (z, θ) y(k) + H −1 (z, θ)G(z, θ)u(k) + e(k)

• Note that 1 − H −1 (z, θ) = h1 z −1 + h2 z −2 + . . . for some coefficients {hi }i=1

∞

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 174

Fit criterion
• For e(k) = 0 (=the best estimate of e(k) we can make), the one-step ahead
prediction of y(k) based on previous measurements is


ŷ(k|θ) = 1 − H −1 (z, θ) y(k) + H −1 (z, θ)G(z, θ)u(k)

• Assuming we have enough data (N ≥ max(na , nb , nc , nd )), we compute the

residual
ϵ(k|θ) = y(k) − ŷ(k|θ)
• The most used fit criterion is
N −1
1 X 2
V (θ, Z N ) = ϵ (k|θ)
N
k=0

• The optimal vector θ∗ is determined by solving the optimization problem

θ∗ = arg min V (θ, Z N )

θ

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 175

Complexity

• Ideally ϵ should depend linearly on θ, so we can get the explicit solution of a

least-squares problem. This only happens for ARX models

• Besides choosing the model structure (ARX, ARMAX, etc.) we also need to
decide the order of the model, i.e., the number of free parameters to optimize

• A small number of parameters could make the model too simple, and not able to
explain the data

• A large number of parameters could make the model more complex than we
need and overfit the data in Z N , resulting in poor predictions on new data

• How to choose the right model complexity?

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 176

Model validation
• Usually to avoid being fooled by overfitting the data set we split Z N in two
subsets: estimation data Zest and validation data Zval :
– Zest is used to compute the optimal parameter vector θ ∗
– Zval is used to see how the estimated model behaves on fresh data
• A validation criterion is to look at one-step prediction errors
N −1
1 X nothing to optimize here,
V (θ∗ , Zval ) = (y(k) − ŷ(k|θ∗ ))2
N just substitute θ ∗ , Zval and evaluate
k=0

• Another validation criterion is to simulate the model completely in “open-loop”

ysim (k, θ∗ ) = G(z, θ∗ )u(k)

and to look at
N −1
1 X
(y(k) − ysim (k, θ∗ ))2
N
k=0

(or just observe how much the plots of y(k) and ysim (k, θ ∗ ) differ)

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 177

Model validation - Residual analysis

• Ideally the prediction error (or prediction residual) ϵ(k|θ) should be white
noise and uncorrelated with u(k)

• To test whiteness of ϵ(k|θ) we compute the auto-correlation function

−τ −1
NX
1
Rϵ (τ ) = ϵ(k + τ |θ)ϵ(k|θ)
N
k=0

• To test correlation between ϵ(k|θ) and u(k) we compute the sample covariance

N −1
1 X
Rϵu (τ ) = ϵ(k|θ)u(k − τ )
N
k=τ −1

• Both Rϵ (τ ) and Rϵu (τ ) should be small

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 178

Model selection
• Which model structure to choose (ARX, ARMAX, OE, etc. )? Which model
orders na , nb , nk , etc. ?

• Cross-validation is the procedure that compares the quality of fit of different

models, by validating them on a data set where neither of them was estimated

• Let θ1∗ , …, θs∗ a set of optimal parameters 0.2

V(θ,Zest)

for different model structures 0.15

0.1

• The best model θi∗ is the one for which 0.05

V (θi∗ , Zval ) is smallest 0

1 2 3 4 5 6 7 8 9 10

• Often Vi (θ∗ , Zest ) decreases as the model

4.5
V(θ,Zval)
4

complexity increases, while Vi (θ∗ , Zval ) 3.5

starts increasing when the model 2.5

complexity becomes excessive (=overfit of 1 2 3 4 5 6 7 8 9 10

model index i (growing complexity)

estimation data)

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 179

Model selection
• If fresh validation data are not available (=no cross-validation), we can use the
same performance figures, but in addition penalize overfit (we want a good
balance between simplicity and accuracy)
• Let di = number of elements of θi∗ (=model complexity)
• We look for the model that minimizes one of the following figures:

– Akaike’s Information theoretic Criterion (AIC):

 
1 X 2
N −1
2di
1+ ϵ (k|θi∗ )
N N
k=0

– Akaike’s Final Prediction Error (FPE):

!
1 X 2
di N −1
1+ N
ϵ (k|θi∗ )
1− di
N
N
k=0

– Rissanen’s Minimum Description Length (MDL):

1 X 2
N −1
di · ln N
ϵ (k|θi∗ ) +
N N
k=0

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 180

Identification example

• Assume that the real (unknown) process is

0.03726z −1 − 0.09676z −2 + 0.08355z −3 − 0.024z −4

G(z) =
1 − 3.464z −1 + 4.493z −2 − 2.586z −3 + 0.5577z −4
with sample time T = 0.04 s
y1
0.2

0.1

• Input excitation: PRBS sequence −0.1

−0.2
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5
Time

• We have 200 samples. The first 100 samples 1

u1

are used for estimation of θ∗ , the rest for 0.5

validation −0.5

−1
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5
Time

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 181

Identification example (cont'd)
We try five different ARX model structures ARX(na ,nb ,nk ):

A(z)y(t) = B(z)u(t) + e(t)

Nval −1
1 ∑ 1 ∑−1
Nest
ϵ2 (k|θi∗ ), ϵ2 (k|θi∗ )
Nval k=0
Nest k=0
−3
x 10
1.4

• i = 1: ARX(1,1,1) V(θ,Zest)
V(θ,Zval)

• i = 2: ARX(2,2,1)
1.35

• i = 3: ARX(3,3,1) 1.3

• i = 4: ARX(4,4,1) 1.25

• i = 5: ARX(5,5,1) 1.2

• i = 6: ARX(6,6,1) 1.15

1.1
1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6

model structure index i

The best model structures are i = 2, 3, 4

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 182
Identification example (cont'd)
Residual analysis:

Correlation function of residuals. Output y1 Correlation function of residuals. Output y1 Correlation function of residuals. Output y1
1 1 1

0.8 0.8 0.8

0.6 0.6 0.6

0.4 0.4 0.4

0.2 0.2 0.2

0 0 0

−0.2 −0.2 −0.2

−0.4 −0.4 −0.4

0 5 10 15 20 25 0 5 10 15 20 25 0 5 10 15 20 25
lag lag lag

Cross corr. function between input u1 and residuals from output y1 Cross corr. function between input u1 and residuals from output y1 Cross corr. function between input u1 and residuals from output y1
0.3 0.3 0.3

0.2 0.2 0.2

0.1 0.1 0.1

0 0 0

−0.1 −0.1 −0.1

−0.2 −0.2 −0.2

−0.3 −0.3 −0.3

−0.4 −0.4 −0.4

−25 −20 −15 −10 −5 0 5 10 15 20 25 −25 −20 −15 −10 −5 0 5 10 15 20 25 −25 −20 −15 −10 −5 0 5 10 15 20 25
lag lag lag

ARX(2,2,1) ARX(3,3,1) ARX(4,4,1)

A(z) = 1 − 0.2153z −1 − 0.5624z −2 A(z) = 1 + 0.1228z −1 − 0.3396z −2 A(z) = 1 + 0.1451z −1 − 0.319z −2

− 0.4444z −3 − 0.4258z −3 − 0.03208z −4

B(z) = 0.04041z −1 + 0.02456z −2

B(z) = 0.04014z −1 + 0.037z −2 B(z) = 0.03912z −1 + 0.03826z −2
+ 0.02247z −3 + 0.02476z −3 + 0.004177z −4

not much different from

ARX(3,3,1)

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 183

Concluding remarks on system identification
• System identification and control design are complementary: no controller
without a model, but identified model only useful for control and/or estimation

• If model parameters change on-line, one can use adaptive control, by

identifying the model and changing the controller accordingly in real-time
(caveat: closed-loop stability may be an issue )

• If linear model structures are not able to capture the model well, one should use
nonlinear models, like artificial neural networks, piecewise affine functions,
and other general function approximation methods available in machine
learning (Schoukens, Ljung, 2019)

• In general, the more a-priori knowledge of the process we can exploit (e.g., from
physical principles), the better. Sometimes black-box identification fails
because its very difficult to guess the right model structure

©2021 A. Bemporad - ``Identification, Analysis, and Control of Dynamical Systems'' 184