INSE 6640: Smart Grids and Control System

Security

Lecture 7 - Mathematical Modeling of Physical Systems

Prof. Walter Lucia

Fall 2020

Fall 2020 1 / 54
Outline

1 Motivations

2 Introduction to Dynamical Systems

3 Linear Systems

4 Classes of Dynamical Systems

5 Equilibrium Points and Stability

Fall 2020 2 / 54
What are we going to study in the second part of the
class?

• We have seen that control systems play a major role in Smart Grid
operations;
• A smart grid is a “cyber-physical system” (CPS), where we have
physical systems, cyber components and communication
infrastructures. Such components are highly coupled;
• Security of Smart Grid is not only a cyber problem. Cyber-attacks
have physical consequences.
• In the first part of the class (Lectures 1-6) we have seen cyber
models and cyber security measures for Smart Grid. In the
second part of the class (Lectures 8-11), we will investigate
control systems security.

Fall 2020 3 / 54
Networked Control Systems in Smart Grid

• From an hight level point of view, control systems in smart grids

are networked control systems.

• It consists of 4 main ingredients:

• The physical plant, the controller logic, the state estimator and
the network infrastructure;
• To study the effects of attacks on control systems we need to know
how all these subsystems behave.

Fall 2020 4 / 54
Security Questions in Networked Control Systems

• We would like to answer the following questions:

• How the output y changes if a False Data Injection (FDI) attack
affect the communication channels?
• Can we design the control module to physically detect the presence
of attacks?
• Can we design the control logic to make the system resilient to
cyber-attacks?

Fall 2020 5 / 54
Security Questions in Networked Control Systems

Can we answer the previous questions? We will proceed as follows:

• Background on control theory: 2 lectures (Lectures 8-9),
• Detection of FDI attacks: 1 lecture (Lecture 10)
• Overview of existing solutions for control system security
(Lecture 11)

Fall 2020 6 / 54
Physical Plant (System) in CPS

• A physical system is an entity (or a set of entities) that evolves

over time, possibly under external excitation

Input Output

Input Output
Input Output Input Output

Input

Output
Input Output Input

Output

Fall 2020 7 / 54
Physical Plant Dynamical Model

• The way the system evolves over time is called the dynamics of
the system
• A dynamical model of a system is a set of mathematical
equations that explain how the system evolves over time (y(t)),
under the effect of an external excitation u(t) or autonomously (if
u(t) ≡ 0)

Fall 2020 8 / 54
Why is a Dynamical Model Useful?

• A mathematical model is useful in several ways:

1 It help to understand how inputs influence outputs
2 It is useful to design the control strategy in such a way that the plant
could reach the desired performance
3 It help us to understand if the system is behaving how we expect
Do we have a physical fault in the system?
Do we have a cyber-attack affecting the command inputs u(t) or the
measurements y(t) ?

Fall 2020 9 / 54
How can we derive a dynamical model of a system?

• This dynamical model can be derived in two main ways:

1 Knowledge-based methods: A mathematical model is derived by
resorting to the physic laws governing the considered process, e.g.,
Kirchhoff’s voltage and current laws
2 Identification based methods: The physical laws governing the
process dynamics are unknown (black box) or partially known (gray
box)

Fall 2020 10 / 54
Introduction to Linear Systems

Fall 2020 11 / 54
What is it a linear system?

• A linear system is a dynamical system which mathematical model

y(t) = f (x(t), u(t)) satisfies the superposition principle:

Recap: Linear Function

• A function y = f (u) is linear if it has two properties:
(
f (αu) = αf (u)
f (u1 + u2 ) = f (u1 ) + f (u2 )

• These two conditions could be combined together (superposition

principle)
f (αu1 + βu2 ) = αf (u1 ) + βf (u2 )

• Is y(t)=3u(t) a linear function?

Fall 2020 12 / 54
Dynamical Model of a Linear System

A linear dynamical system is represented by a system of n ordi-

nary differential equations of 1st order.

• Example: If the system has 1 single-input u(t), one single-output

y(t) and state vector x(t) ∈ IRn , then it can be represented as
follows:

ẋ1 (t)

 = a11 x1 (t) + · · · + a1n xn (t) + b1 u(t)
..
 .

ẋ (t)
n = an1 x1 (t) + · · · + ann xn (t) + bn u(t)
y(t) = c1 x1 (t) + · · · + cn xn (t) + d1 u(t)
x1 (0) = x10 , · · · , xn (0) = xn0
Fall 2020 13 / 54
Linear System: State-Space Representation

A linear dynamical system can always be represented in the following

matrix form (state-space form):
(
ẋ(t) = Ax(t) + Bu(t)
y(t) = Cx(t) + Du(t)

where:
• x(t) := [x1 , x2 , · · · , xn ]T ∈ IRn is state vector of the system
• u(t) := [u1 , u2 , · · · , ur ]T ∈ IRr is input vector of the system
• y(t) := [y1 , y2 , · · · , ym ]T ∈ IRm is output vector of the system
• A ∈ IRn×n , B ∈ IRn×r , C ∈ IRm×n , and D ∈ IRm×r

Fall 2020 14 / 54
Linear Systems: block scheme representation

(
ẋ(t) = Ax(t) + Bu(t)
y(t) = Cx(t) + Du(t)

Fall 2020 15 / 54
Example: mass-spring-damper system (1/4)

• Find the state space representation of the mass-spring-damper

system

Fall 2020 16 / 54
Example: mass-spring-damper system (2/4)

• To solve this example, we use Newton’s Second Law

(knowledge-based method):
X
F~i (t) = M z̈(t)
i
⇒ u(t) = M z̈(t) + β ż(t) + Kz(t)
• We define the states of this system as follows: x1 (t) = z(t) and
x2 (t) = ż(t) (other equivalent choices for x1 and x2 are possible)
( (
ẋ1 (t) = ż(t) x˙1 (t) = x2 (t)
⇒ β K 1
ẋ2 (t) = z̈(t) x˙2 (t) = − M x2 (t) − M x1 (t) + M u(t)
Fall 2020 17 / 54
Example: mass-spring-damper system (3/4)

( (
ẋ1 (t) = ż(t) x˙1 (t) = x2 (t)
⇒ β K 1
ẋ2 (t) = z̈(t) x˙2 (t) = − M x2 (t) − M x1 (t) + M u(t)

• Now we are able to rewrite the equation in a matrix form

• Consider x(t) = [x1 (t), x2 (t)]T
   
0 1 0
ẋ(t) = k β x(t) + 1 u(t)
−M −M M
| {z } | {z }
 A
 B
y(t) = 1 0 x(t)
| {z }
C

Fall 2020 18 / 54
Example: mass-spring-damper system (4/4)

1
u=1 xdot s x y
Input x0=[1;0] Measurements
B C

Fall 2020 19 / 54
How can we obtain the states evolution of the system?
(
ẋ(t) = Ax(t) + Bu(t)
y(t) = Cx(t) + Du(t)

• We need to solve the above system of ODEs;

• The solution is given by the following Lagrange’s Formula

Proposition
Given the continuous-time linear system ẋ = Ax + Bu, with initial
condition x(0) = x0 ∈ IRn , there exist a unique solution x(t)

forced response
natural response z }| {
z }| { Z t
At A(t−τ )
x(t) = e x0 + e Bu(τ )dτ
Z t 0 
At A(t−τ )
y(t) = C(e x(0)) + C e Bu(τ )dτ + Du(t)
0

Fall 2020 20 / 54
Example: Lagrange’s Formula

Evolution of the system = Evolution due to the initial condition +

Evolution due to the input
Z t 
At A(t−τ )
y(t) = C(e x(0)) + C e Bu(τ )dτ + Du(t)
0

1
xdot s x ytotal
Input u=1 x0=[1;0] y
B C

1
u=0 xdot s x
x0=[1; 0] y due to the initial condition x0
B1 C1
Comparison

A1 xdot

1 Lagrange Formula
u=1 xdot s x
Input u =1 x0=[0;0] y due to the external input
B2 C2

Fall 2020 21 / 54
Math Recap: Eigenvalues and Eigenvectors

Fall 2020 22 / 54
Eigenvalues and Eigenvectors

• Given a square matrix A ∈ IRn×n , the eigenvalues are the roots

of its characteristic polynomial

det(λI − A) = λn + βn−1 λn−1 + · · · + β1 λ + a0 = 0

• An eigenvector of A is any vector vi ∈ IRn such that

Avi = λi vi

• The algebraic multiplicity of λ1 is equal to the number of

coincident roots λi
• The geometric multiplicity of λi is the number of linear
independent eigenvectors associated to λi

Fall 2020 23 / 54
Diagonalizable Matrix

• If for each eigenvalue of A the algebraic multiplicity is equal to the

geometric multiplicity, then matrix A is said diagonalizable

A = T ∆T −1

 
λ1 0 . . . 0
 0 λ2 . . . 0 
∆=  = T −1 AT, T = [v1 |v2 . . . |vn ]
 
.. .. . . .
 . . . .. 
0 0 . . . λn

Fall 2020 24 / 54
Example - Mass-Spring-Damper
• Find the eigenvalues and eigenvectors if M = 25, K = 24, and
β=8

   
0 1 0
ẋ(t) = k β x(t) + 1 u(t)
−M −M M
|{z } | {z }
A B
y(t) = 1 0 x(t)
| {z }
C

Fall 2020 25 / 54
Example - Mass-Spring-Damper: Eigenvalues
• By using the MATLAB command “eig(A)” we obtain
• V, reading by column, contains the eigenvectors
• D, on the diagonal, contains the eigenvalues

Fall 2020 26 / 54
Why eigenvalues and eigenvectors are so important?

• Let us consider the following autonomous system (u ≡ 0 )

ẋ = Ax

• By resorting to the Lagrange’s formula, the evolution of the system

is given by the natural response:

x(t) = eAt x0

Fall 2020 27 / 54
Why eigenvalues and eigenvectors are so important?
• If A is diagonalizable we can rewrite the response as
 λt 
e 1 0 ... 0
 0 e λ2 t . . . 0 
x(t) = eAt x0 = T e∆t T −1 x0 = [v1 . . . vn ]  . α
 
.. .. ..
| {z }  .. . . . 
α
0 0 . . . eλn t
 
α1 n
λn t  .. 
X
λ1 t
= [v1 e . . . vn e ]  .  = αi eλi t vi
αn i=1
x(t)
• The evolution depends upon the eigenvalues λi of A
• The exponentials eλi (t) are called the modes of the system
• E.g. What will happen over time if x(t) = 2e−3t + 2e5t ?

Fall 2020 28 / 54
 Classes of Dynamical Systems

• Beside linear systems, we can find other classes of dynamical

systems:

Fall 2020 29 / 54
Other classes of dynamical systems (1/2)

• Let us consider the following linear continuous-time state-space

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

• If A, B, C, D are all constant, the system is said linear-time

invariant (LTI)
• If any matrix A, B, C, D is not constant, the system is said
linear-time varying (LTV)

Fall 2020 30 / 54
Other classes of dynamical systems (2/2)

• The system
ẋ(t) = f (x(t), u(t))
y(t) = g(x(t), u(t))
with

f : IRn+m → IRn , g : IRn+m → IRp , arbitrary nonlinear functions

defines a time-invariant nonlinear system.

• The system
ẋ(t) = f (t, x(t), u(t))
y(t) = g(t, x(t), u(t))
defines a time-varying nonlinear system.

Fall 2020 31 / 54
Equilibrium Points and Stability

Fall 2020 32 / 54
Equilibrium points of generic nonlinear systems

• What is an equilibrium point?

• Let us consider the continuous-time
nonlinear system
ẋ(t) = f (x(t), u(t))
y(t) = g(x(t), u(t))

Definition
A state xe ∈ IRn and an input ue ∈ IRm are an equilibrium pair if for an
initial condition x(0) = xe and constant input u(t) ≡ ue , the state
remains constant, i.e.
x(t) ≡ xe , ∀t ≥ 0

Fall 2020 33 / 54
How to find equilibrium points?
ẋ(t) = f (x(t), u(t))
y(t) = g(x(t), u(t))

Definition
A state xe ∈ IRn and an input ue ∈ IRm are an equilibrium pair if

f (xe , ue ) = 0

Definition
If xe ∈ IRn and an input ue ∈ IRm are an equilibrium pair the
• xe is called equilibrium state
• ue is called equilibrium input

• For nonlinear system we can have 0,1,..many...infinitely many

equilibrium pairs

Fall 2020 34 / 54
Stability of Equilibrium Points

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

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

• How does the system behave when we start from an initial

condition non exactly on the equilibrium point?
• What will happen if the system is in an equilibrium configuration
and an attack create a little perturbation?

Fall 2020 35 / 54
Stability of Equilibrium Points

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

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

• Different equilibria can have different stability

• Stability is not in general a system property but it refers to a
specific equilibrium pair (local)
• The equilibrium pair (xe , ue ) can be
• Stable
• Asymptotically Stable
• Unstable

Fall 2020 36 / 54
Stable and asymptotically stable

Figure: Stable

Figure: Asymptotically stable

1
1
https://www.math24.net/stability-theory-basic-concepts/
Fall 2020 37 / 54
Equilibrium Points and Stability of Linear Systems

Fall 2020 38 / 54
Equilibrium points of linear systems

(
ẋ(t) = Ax(t) + Bu(t)
y(t) = Cx(t) + Du(t)

• Equilibrium points:
• Pairs (xe , ue ) such that Axe + Bue = 0
• (xe , ue ) = (0, 0) is a obvious equilibrium, valid for any linear system

Fall 2020 39 / 54
Stability of linear systems
(
ẋ(t) = Ax(t) + Bu(t)
y(t) = Cx(t) + Du(t)

• Stability of Equilibrium Points:

• All the equilibria share the same type of stability (system property
not local)
• Stability/Instability depends only on the eigenvalues of A

Theorem
Let λ1 , . . . , λm , m ≤ n be the eingenvalues of A ∈ IRn×n .
• Asymptotically stable if and only if (iff) <[λi ] < 0, ∀i
• Unstable if ∃ i such that <[λi ] > 0 or ∃ <[λi ] = 0 with algebraic
multiplicity different from the geometric multiplicity.
• Stable or Marginally Stable if <[λi ] ≤ 0, ∀i and the eigenvalues
with <[λi ] = 0 have equal algebraic and geometric multiplicity
Fall 2020 40 / 54
Example - Mass-Spring-Damper

• By using MATLAB we have found the eigenvalues:

λ1 = −0.16 + 0.96i, λ2 = −0.16 − 0.96i

• Since all the eigenvalues have a negative real part we can

conclude that the system is asymptotically stable

Fall 2020 41 / 54
Mass-Spring-Damper: Simulation

For any input (u = 1 in the simulation below), the mass-spring-damper

system will asymptotically settle down (after a transitory period) to an
equilibrium point

Fall 2020 42 / 54
Stability of equilibrium points of Nonlinear Systems

Fall 2020 43 / 54
Stability of equilibrium points of Nonlinear Systems

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

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

• Stability is hard to prove. For different equilibrium points I can

have different stability.
• Methods to prove stability of equilibrium points:
• Lyapunov indirect method
• Lyapunov direct method (we will not see this)
Fall 2020 44 / 54
Lyapunov Indirect Method

Fall 2020 45 / 54
Lyapunov indirect method: Idea
ẋ(t) = f (x(t), u(t))

• We know how to check stability for Linear Systems (i.e.

eigenvalues of A)
• We can approximate the non-linear system as a linear system
around the equilibrium pair of interest by using the Taylor
expansion formula truncated to the first order
∂f (x, u)  ∂f (x, u) 
x̃˙ = f (xe , ue ) + (x − xe ) + (u − ue )
∂x ∂u
 
(xe ,ue ) (xe ,ue )

Fall 2020 46 / 54
Linearization around an equilibrium point
• Given a NLTI system ẋ(t) = f (x(t), u(t))
• Linearizing the system, around (xe , ue ), using the Taylor series of
f (x, u) approximated to the first order derivative
   
∂f (x, u)  ∂f (x, u) 
f (x, u) = f (xe , ue ) + x̃ + ũ
∂x ∂u
 
| {z } (xe ,ue ) (xe ,ue )
=0
we obtain

˙
x̃(t) = Ax̃(t) + B ũ(t)

 x̃ = x − xeq, u − ueq
ũ = 
 
A = ∂f ∂x
(x,u) 
 , B= ∂f (x,u) 
∂u (x ,u )
(xe ,ue ) e e

Can we understand the stability of (xe , ue ) by studying the

properties of the linearized model?

Fall 2020 47 / 54
Lyapunov indirect method: Theorem

Theorem
Consider a NLTI system ẋ = f (x) with f differentiable, xe = 0 an
equilibrium point. Given the LTI approximation of ẋ = f (x)

˙ ∂f
x̃(t) = Ax̃(t), A= (x − 0)
∂x

• If x̃˙ = Ax̃ is asymptotically stable, then x = 0 is an asymptotically

stable equilibrium point for the NLTI system
• If x̃˙ = Ax̃ is unstable, then x = 0 is an unstable equilibrium for the
NLTI system
• If x̃˙ = Ax̃ is marginally stable , nothing can be said about the
stability of x = 0 for the NLTI system

Fall 2020 48 / 54
Example - Pendulum

Next Lecture
The non-linear state-space representation of the pendulum is:

ẋ(t) = f (x(t))

where
   
ẋ1 (t) x2 (r) h
ẋ(t) = , f (x(t)) = g , H :=
ẋ2 (t) − l sin x1 (t) − Hx2 (t), ml2

and x1 (t) = y(t) and x2 (t) = ẏ(t).

Fall 2020 49 / 54
Pendulum - Equilibrium Points

• The equilibrium points of the system are obtained by the following

condition (f (xe ) = 0):
    
x2 0 x2e = 0
= ⇒
− gl sin x1 − Hx2 0 x1e = ±kπ, k = 0, 1, . . .

Next Lecture

Fall 2020 50 / 54
Pendulum - Stability of first equilibrium points

First case x2e = 0 and

x1e = kπ, k = 0, 2, 4, . . . ⇒
Linearized Model:
 
˙x̃(t) = 0 1
x̃(t)
− gl −H
Next Lecture

• The eigenvalues of the linear system matrix A are

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

• Re[λ1,2 ] < 0 ⇒ the equilibrium is asymptotically stable

Fall 2020 51 / 54
Pendulum - Stability of first equilibrium points

Second case x2e = 0 and

x1e = kπ, k = 1, 3, 5, . . . ⇒
Linearized Model:
 
˙x̃(t) = 0g 1
x̃(t)
l −H
Next Lecture

• The eigenvalues of the linear system matrix A are

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

• Re[λ1 ] < 0 and Re[λ2 ] > 0 ⇒ the equilibrium is unstable

Fall 2020 52 / 54
Thank you!

Fall 2020 53 / 54
References

Slotine, Jean-Jacques E and Li, Weiping and others

Applied nonlinear control
Prentice hall Englewood Cliffs, NJ 199(1), 1991.

Antsaklis, Panos J and Michel, Anthony N

A linear systems primer
Birkhäuser Boston, 2007.
Emilio Frazzoli, and Munther Dahleh
Dynamic Systems and Control
Massachusetts Institute of Technology: MIT OpenCourseWare, spring 2011.

A. Bemporad
Lecture Slides: Identification, Analysis and Control of Dynamical Systems
IMT School for Advanced Studies Lucca, 2018.

Fall 2020 54 / 54