ELP7100: Autonomous Systems & Control Lab

Introduction to Scientific Computing Softwares

Course Instructors: Subashish Datta subashish@ee.iitd.ac.in

Shaunak Sen shaunak@ee.iitd.ac.in
Course TAs: Shuvam Routray, Khushall Gourav, Anvisha Chaturvedi

1. Similarity Transformation (Controller Canonical Form):

Consider a system:
ẋ = Ax + Bu
where x ∈ Rn is the state vector and u ∈ Rm is the control input. Assume that the pair (A, B) is control-
lable. Transform this system to an equivalent one:

x̄˙ = Āx̄ + B̄u

such that Ā = PAP−1 , B̄ = PB are in controller canonical form.

Steps to obtain Ā & B̄:

(a) Compute the controllability matrix:

C = B AB A2 B · · · An−1 B
 

= b1 b2 · · · bm Ab1 Ab2 · · · Abm · · · An−1 b1 An−1 b2 · · · An−1 bm

 

(b) Starting from the left and moving towards the right, select n independent columns of C . This is
always possible since the pair (A, B) is controllable.
(c) Reorder these columns by taking first b1 , Ab1 .... untill all colums involving b1 have been used.
Then take b2 , Ab2 ... till all columns involving b2 are used and then lastly use bm , Abm .
(d) Construct the following matrix:

C ′ = b1 Ab1 · · · Aµ1 −1 b1 b2 Ab2 · · · Aµ2 −1 b2 · · · bm Abm · · · Aµm −1 bm .

 

The positive integers µi are called the controllability indices of the system. Further,

µ := max µi
i

is called the controllability index of the system.

(e) Define
k
σk = ∑ µi k = 1, 2, . . . , m
i=1

i.e.
σ1 = µ1
σ2 = µ1 + µ2
..
.
σm = µ1 + µ2 + · · · + µm = n

(f) Compute [C¯]−1 , and denote qk as the σkth row of [C¯]−1 .

1
 (g) Construct the transformation matrix P as:
 
... q1 ...
. . . q 1A . . .
 .. .. .. 
 
 . . . 
. . . q1 Aµ1 −1 . . .
 
 
. . . q2 . . .
 
. . . q2 A . . .
 . .. .. 
 .
 . . . 

. . . q2 Aµ2 −1 . . .
 
 . .. .. 
 .. . . 
 
. . . qm . . .
 
. . . qm A . . .
 
 . .. .. 
 .. . . 
. . . qm Aµm −1 ...

(h) Finally, compute:

Ā = PAP−1 , B̄ = PB

Programming Task: Write code in MATLAB, Scilab, and Python to compute the controller
canonical form. Your code should work for any (A, B) controllable pair for any n and m. You
should not write this code for a particular example. Observe the structure of Ā and B̄.
• Input: Matrices A, B
• Output:
(a) Canonical matrices Ā, B̄
(b) Controllability indices µ1 , . . . , µm

2. Feedback Control Design

Given a controllable system:

ẋ = Ax + Bu, x(0) = x0 , x ∈ Rn , u ∈ R, (m = 1)

Choose any system of size n ≥ 4, m = 1. First, convert the system into controller canonical form:
x̄˙ = Āx̄ + B̄u using the previous algorithm (Problem 1). Then, apply state feedback control:

u = F̄ x̄

where F̄ ∈ R1×n is the feedback gain vector. The closed-loop system becomes:


x̄˙ = Ā + B̄F̄ x̄

Design Objective:


• Choose F̄ such that the eigenvalues of Ā + B̄F̄ satisfy:
(a) Damping ratio ζ ≥ 0.7,
(b) Decay rate λ ≤ 1.
• Compute F = F̄P

2


• Compute the eigenvalues of (A + BF), and then compare them with the eigenvalues Ā + B̄F̄ .
Comments on your observations.
• Simulate the closed-loop system for an initial condition x0 ∈ Rn .

Simulation Task: Implement the controller and simulate the closed-loop response in:
• MATLAB
• Python
• Scilab
Submit:
• Code file (.m, .py, .sce)
• Plot showing system response (state trajectories over time)
• Choice of closed-loop poles and justification.

3. Numerical Solution of Non-linear Systems

An equivalent circuit of the Wien-Bridge oscillator is shown below:

The non-linear voltage controlled voltage source is given by:

g(v2 ) = 3.23v2 − 2.195v32 + 0.666v52

Let x1 = v1 and x2 = v2 . Then, a state-space model of the circuit is:

1
ẋ1 = (−x1 + x2 − g(x2 ))
C1 R1
1 1
ẋ2 = − (−x1 + x2 − g(x2 )) − x2
C2 R1 C2 R2

Assume that: C1 = C2 = R1 = R2 = 1. Then:

(i) Solve the above set of equations to obtain x1 and x2 using numerical integration methods. Do not
use any existing ODE solver command. Let the solutions be x̄1 and x̄2 .
(ii) Solve these equations for x1 and x2 using existing commands (e.g., built-in ODE solvers). Plot x1
and x2 . Also plot x̄1 and x̄2 (obtained from part (i)). Compare these plots.

Programming Task: For both parts (i) and (ii), use MATLAB, Scilab, and Python.

[Example taken from Book: Nonlinear Systems, by Khalil, Chapter 2]

3
4. Modeling and Simulation with OpenModelica:
Using OpenModelica, perform modeling and simulation of the following systems:

• Simple electric circuits

• Mass-Spring Damper system
• An electro-mechanical system

Analysis Task: Comment on the type of model you have obtained in comparison to the ordinary state-
space model.

5. Consider the following model of a simple pendulum θ̈ + sin θ = 0. Define state variables as x1 = θ and
x2 = θ̇ . Rewrite the pendulum model as a set of two first-order differential equations.
We will explore numerical integration of these equations. The simplest method is the Euler method
that approximates the derivative as ẋ1 ≈ (x1 (t + h) − x1 (t))/h, where h is a small number. Using this
approximation and a similar one for ẋ2 , rewrite the two first-order differential equations as corresponding
difference equations.
The two first-order difference equations can be numerically solved in a computer for any initial condition
x1 (0) and x2 (0) and a finite time interval [0, T ]. Divide the time interval in T /h steps of h. Denote
xi (nh) ≡ xn for i = 1, 2 and n = 0, . . . T /h. Calculate the approximate trajectory for x1 (0) = 1, x2 (0) = 0,
T = 10, h = 0.1, 0.01, 0.001. What are your observations for different values of h? Compare this with a
standard ODE solver in MATLAB or Python.

6. Many control applications are to heat transfer. The heat equation in one dimension is

∂T ∂ 2T
= . (1)
∂t ∂ x2
For the x domain 0 ≤ x ≤ 1, time interval 0 ≤ t ≤ 1, the initial condition T (x, 0) = 37 and boundary
conditions T (0,t) = 25 = T (1,t), we wish to obtain the spatiotemporal profile of the temperature T .
Consider a discretization
1 1
(Ti, j+1 − Ti, j ) = 2 (Ti+1, j − 2Ti, j + Ti−1, j ), (2)
k h
where h is the step size in the x-direction, k is the step size in the t-direction, and T (ih, jk) ≡ Ti, j . Solve
for the temperature profile for different values of h and k. What are your observations?

7. If one puts a bacterial cell in a medium, it grows and divides. This process of growth continues until the
resources are exhausted. Suppose we have measurements of the growth (Ok ) at discrete instants of time
(tk , attached data file). The broad question is to estimate the growth rate.
(a) Calculate the growth rate by taking a discrete approximation γk = (Ok+1 − Ok )/(tk+1 − tk ).
(b) Treat the growth rate as a parameter to be estimated from a model Ṅ = γ(t)N, where N(t) is the
number of cells at time t. Use an Extended Kalman Filter, to estimate γ from output measurements
y(tk ) = N(tk ) (second column in Excel File).
The “Extension” in the Kalman Filter arises through a redefinition of the parameter γ as a state. So
the two states are N and γ and the corresponding equations are Ṅ = γN, γ̇ = 0. The equation may be
linearized around a point (N ∗ , γ ∗ ) as ∆N = N ∗ ∆γ + γ ∗ ∆N, ∆γ = 0, where ∆N = N − N ∗ and ∆γ = γ − γ ∗ .
This way the equation is in a standard linear form ẋ = Ax and y(tk ) = Cx(tk ). The Kalman Filter for this
is:

4
 • Between Observations (tk < t < tk+1 )

x̂˙tt = Axtt , (3)

Ṗtt = APtt + Ptt A′ + Q. (4)

• At an Observation (t = tk )
+ − −
t t t
x̂˙tkk = xtkk + K(tk )(y(tk ) −Cxtkk ), (5)
tk+ tk− tk−
Ṗ = P − K(tk )CP .
tk tk tk (6)

Here, the x̂’ s represent the state estimates, P’s represent the estimate covariances:
t− t−
K(tk ) = Ptkk C′ (CPtkk C′ + R)−1

is the Kalman Gain, Q ∈ R2×1 is the process noise covariance, and R ∈ R1×1 is the measurement noise
covariance.
(c) Compare the growth rate estimated in (a) and (b) above.