MTHM606 Applied AI and Control
Coursework Part 2: Control
For your second piece of coursework you will need to submit your answers to questions 1–3 and
either question 4 or question 5. Your answers must be submitted in pdf format to BART and
your original code files submitted to the module ELE page. You are strongly recommended to do
this by creating a MATLAB Live Script (.mlx) which can include both text, code and screenshots
of your Simulink block diagrams: submit the original .mlx file to ELE and generate a pdf file to
submit to BART. You are encouraged to make use of the code provided in lectures as a basis (it
needs to be appropriately referenced). This piece of work is worth 30% of the module mark.
1. Dynamical Systems and Introduction to Control. Consider the following model of
the spread of an infectious disease:
dS
dt = β(1 u)SI
dI
dt = β(1 u)SI γI δI,
in which S ≥ 0 and I ≥ 0 represent the proportion of susceptible and infected individuals
in the population, respectively; β, γ and δ are constant-valued parameters taking the values
β = 0.8, γ = 0.1, δ = 0.02; and u ≥ 0 is considered as an actuation input representing the
effects of social distancing measures. Assume that the proportion of infected individuals at
time t = 0 is equal to p = 0.01, i.e., S(0) = 1 p and I(0) = p.
(a) Simulate the model in Matlab or Simulink (Simulink is recommended) for the case in
which u steps up from a value of 0 at time t < 5 to the value of 0.8 for time t ≥ 5 (i.e.,
u(t) = 0 for t < 5 and u(t) = 0.8 for t ≥ 0). You should provide plots of S(t) and I(t)
in the interval 0 ≤ t ≤ 50, in addition to the details of your Matlab code and Simulink
block diagrams used to generate these plots.
(b) Verify that the system has an equilibrium point when I = 0, S = S and u = u for any
given S ≥ 0 and u ≥ 0.
(c) Find a linearised approximation to the system in the vicinity of the equilibrium point
in question 1b.
In preparing your answer, you are advised to refer to the Introduction to Control and Intro_xFFFE_duction to Dynamical Systems practicals.
2. Controllability and Stabilisability. The rotational motion of the Lego Gyroboy robot
is well approximated by the linear time-invariant system
dx
dt =
0 1
0
βW2
2J R2 + (mR2 + Jw)W2
x +
0
2RαW
2J R2 + (mR2 + Jw)W2
u,
1
where
x =
"
d
dt
#
,
and where
J =
M(W2 + D2
)
12
and Jw =
mR2
2
,
with
m = 0.031, M = 0.641, D = 0.1, W = 0.105,
R = 0.027, α = 0.000505, and β = 0.00358.
(a) Design an actuation input u for the rotational motion of the robot that transitions the
robot from a rotational angle of 0 at time t0 = 0 through a rotational angle of π
8
radians (corresponding to a sixteenth of a revolution) in 0.5 seconds, and where the
rate of change of rotational angle is zero at both ends of this interval, i.e.,
(0) = 0,
d
dt (0) = 0, (0.5) = π
8
, and d
dt (0.5) = 0.
You should plot your solution for u over the interval 0 ≤ t ≤ 0.5 and also provide your
Matlab code used to compute this. You are advised to use the function cont x0 to x1
and associated functions from the Controllability and Stabilisability practical.
(b) Simulate the state trajectory resulting from the application of the actuation input u
obtained in the first part of this question. This could be done either in Simulink using
the From Workspace block, or in Matlab using the function lti sys sim. You should
plot the resulting state trajectory (for both and d
dt ) over the time interval 0 ≤ t ≤ 0.5,
provide your Matlab code, and also provide a screenshot and brief description of your
Simulink block diagram if you have used Simulink.
3. Observability and Detectability. The rotational motion of the Lego Gyroboy robot is
well approximated by the linear time-invariant system
dx
dt =
0 1
0
βW2
2J R2 + (mR2 + Jw)W2
x +
0
2RαW
2J R2 + (mR2 + Jw)W2
u,
where the state x and the parameter values are as defined in question 2. It is desired to
estimate the state of this system using the meaurement
y =
1 0 x.
Design an observer for this syst
em using the eigenvalue placement technique, placing the
eigenvalues that determine the speed of convergence of state estimation error at 3 and 4.
You should provide your Matlab code or calculations used to design the observer, and also
state the differential equations governing the observer, specifying the numerical values of the
matrices in these differential equations.
2
4. Optimal Control. A model of infectious disease dynamics within a small closed com munity is described by the equations
dS
dt = b(S + E + I + R) dS cSI uS,
dE
dt = cSI (d + e)E,
dI
dt = eE (a + d + g)I,
dR
dt = gI dR + uS,
where S, E, I and R denote the numbers of susceptible, exposed, infected and recovered indi viduals (but won’t be constrained to integer values!). Here, u is a control input representing
the effects of vaccination, which is subject to the constraint
0 ≤ u(t) ≤ 0.9 for all 0 ≤ t ≤ T.
It is desired to design a vaccination strategy (i.e., a control input u) in order to minimise
the cost function
Z
0
T
2I(τ ) + u(τ )
2
dτ,
subject to the initial condition
S(0)
E(0)
I(0)
R(0)
=
1000
100
50
15
.
The parameter values are as follows:
b = 0.525, d = 0.5, c = 0.001, e = 0.5, g = 0.1, a = 0.2, and T = 20.
(a) Use Pontryagin’s maximum principle to show that the optimal input u necessarily
satisfies
u = min 0.9, max 0,
(λ1
2
λ4)S
,
where λ1 and λ4 denote the first and fourth components of the costate
λ =
λ1
λ2
λ3
λ4
.
(b) State the differential equation and boundary condition that must be satisfied by λ.
(c) Hence compute the optimal input and corresponding state and costate trajectories
using the forward-backward sweep method. You should provide details of your code
and plots of the optimal input, state and costate trajectories over the time interval
0 ≤ t ≤ T. You are advised to use the code from the Vaccination control example
in the Optimal Control practical as a starting point, and you will need to alter the
functions state eq.m and costate eq.m, the rule for updating the input u, and the
parameters defining the model (including the initial state x0).
3
5. Measurement feedback. The paper:
A. Kalbat, “Linear Quadratic Gaussian (LQG) Control of wind turbines,” 2013
3rd International Conference on Electric Power and Energy Conversion Systems,
2013, pp. 1-5, doi:10.1109/EPECS.2013.6713053
considers a (infinite horizon) Linear Quadratic Gaussian controller for a wind turbine con troller. The model of the dynamics of the wind turbine (the “plant”) is described by the
differential equation
dx
dt = Ax + Bu + Nw,
y = Cx + Mv.
Here, x =
x1 x2 x3
T
where x1 is proportional to the rotor speed, x2 to the drive train
torsion, and x3 to the generator speed. Also, u is proportional to the collective blade pitch
angle, w represents turbine system noise, and v represents measurement noise. The signals
w and v are considered to be uncorrelated zero-mean unit variance white noise, and the
parameters take the values
A =
0.1445 3.108 0
26.91 0 26.91
0 15.60 0
, B =
3.456
0
0
, C =
0 0 1 ,
N =
0.02496
0
0
, and M = 0.3162.
(a) Find a controller mapping y to u in order to minimise the power spectral density
E
lim
T→∞ T
1
Z
0
T
x(τ )
TQx(τ ) + u(τ )
TRu(τ )dτ
,
where
Q =
1 0 0
0 0.1 0
0 0 1
, and R = 1.
You should provide the solutions X and Y to the relevant Algebraic Riccati equations
and the equations describing your optimal controller.
(b) Simulate the system in Simulink with an end time of T = 1000 and an initial state of
x(0) = 0 0 0 T
, using band limited white noise blocks for the disturbances w
and v and setting the power level of these blocks appropriately. You should provide an
image of your Simulink block diagram, a short description of the block parameters, any
associated Matlab code, and plots of the signals u and y.
You are encouraged to consult the section “Infinite Horizon Linear Quadratic Gaussian
Controller” from the Measurement Feedback practical in completing this question.
4