ECPD-2023-2024 Lab Guide

Master in Electrical and Computer Engineering
Master in Aerospace Engineering

2023/2024
Estimação e Controlo Preditivo Distribuído
Distributed Predictive Control and Estimation
Predictive Control of a Thermal Plant
Picture source: https://apmonitor.com/heat.htm
Prepared by
Afonso Botelho
João Miranda Lemos

Instituto Superior Técnico
Departamento de Engenharia Eletrotécnica e de Computadores
Área Científica de Sistemas, Decisão e Controlo
IST – Estimação e Controlo Preditivo Distribuído – TCLab work 2023/2024 Page 1

Ethics statement
When delivering the report with the resolution of this work, the group of students that
signs it guarantees, by this act, that the text and all the software and results delivered
were entirely carried out by the elements of the group, with a significant participation
of all of them, and that no part of the work or software and results presented were
obtained from other persons or sources.
Objective
Design a model predictive controller and Kalman filter for a real thermal plant with one
input and one output. Perform system identification with real data, design controller
and observer in simulation, and then apply them to the real plant. The first two parts
of this laboratory work provide exercises on basic issues on function optimization and
receding horizon and predictive control.
Pedagogical concept
This laboratory work has been produced in the framework of the pedagogical project
MobilLab, financed by IST. According to the pedagogical concept in this project, the
students receive a lab kit to perform the experiments according to their own schedule.
The kit is then given back to IST at the end of the term. Laboratory classes are used to
present short progress reports by the student and to clarify doubts.
Underlying concepts
• Receding horizon control

• System identification
• Model predictive control
• Constrained optimization
• State estimation
Laboratory guide organization

This laboratory guide is organised in three parts:
1. Function optimization
2. Receding horizon and unconstrained MPC
3. Model identification
4. MPC and Kalman filter design
5. Application to a real system
Hereafter, Pi identifies part i of the work to be done, as listed above.

Kit used
This work uses the TCLab kit. See the document TCLab installation and setup
instructions in order to connect the kit to your personal computer.
The kit is physically connected to the computer by a USB cable. The kit consists of two
connected boards, one with the plant to control and the other an ARDUINO that is
used just for communication between the plant board and the personal computer.
The TCLab kit must also be connected to a power supply. See the manual TCLab
installation and setup instructions for instructions.
The plant consists of two power transistors that receive a command signal each
between 0 and 100% in order to dissipate heat on metallic plates, and two
temperature sensors (thermistors) that provide a temperature reading at each metallic
plate. In the present work, we will only use one heater and one sensor.
After interconnecting the TCLab kit with the personal computer, and installing the
interface software, it is possible to read each of the temperatures with a single
MATLAB command (in degree Celsius), as well as to send commands to each of the
transistor heaters (a number between 0 and 100, with 0 corresponding to no heating
and 100 corresponding to the maximum possible heating).
A typical MATLAB program to do temperature control consists therefore of a time

cycle in which the following activities are done:
1. Read the temperature from the plant (just one command line in MATLAB);
2. Compute the control variable using the predictive control algorithm (this task is
the core of the work to be done);
3. Write the value of the control variable to the plant (just one command line in
MATLAB);
4. Wait for the sampling time to end and go to step 1.
Safety precautions
Never touch the heating elements (see picture). They can be
very hot and injure your skin.
In order not to damage the kit, do not apply high values of

the control, such that the temperature is kept below 60 °C.

Software required
To perform this work, the following software is required:
• MATLAB, available through the IST campus license;

• MATLAB Optimization Toolbox, available through the IST campus license;
• MATLAB Control System Toolbox, available through the IST campus license;
• MATLAB System Identification Toolbox, available through the IST campus
license;
• TCLab MATLAB interface (function tclab.m), available at the Fénix page of the
course.
• The following MATLAB scripts: TCLab_openloop.m, TCLab_identification.m,
TCLab_simulator.m, TCLab_closedloop.m available at the Fénix page of the
course.
The information required to perform this work is described below.
Grades
The final grade of the laboratory is computed as
L=0.25 G1 + 0.15 G2 + 0.1 G3 + 0.3 G4 + 0.2 G5
where Gi is the grade obtained in the part Pi of the work.
Appendix 1 lists the good practices that must be followed when writing the report.
Failure to comply with these practices will imply a penalty on the grade. Appendix 1
can be used as a check list to verify the report and improve it.
Participation in laboratory classes

It is not mandatory for the students to perform the work during the laboratory
sessions. Whenever appropriate, this oral report can be accompanied by a
demonstration.
It is not mandatory for the students to attend the full class, but only during the period
required to present the progress report.
Not all the students must be present. At least 1 student must be present in each class.
During the whole laboratory classes, all the group members must have been present in
the lab at least once. Failure to comply with these guidelines imply a penalty of 1 point
in the lab grade.
Format of the report to deliver

The report to deliver must comply with the following format:

• The report can be written in either Portuguese or English, but only one
language can be used. It is advised that you use the language that you master
better. Poor language will be penalized.
• File format: pdf
• Each group of students delivers a report, signed by all the group members.
• Include in the 1st page:
o The name of the curricular unit and the academic year (ano letivo)
o Student number, name and email of all the students in the group
o An ethical commitment to originality with the following text:
The group of students identified above guarantees that the text of this
report and all the software and results delivered were entirely carried
out by the elements of the group, with a significant participation of all of
them, and that no part of the work or the software and results presented
was obtained from other people or sources.
• The answers must be given sequentially from page 2 of the report, indicating at
the beginning of each one the respective number (P3, P4, P5) in bold type.
• The maximum number of pages of the report to be delivered, including the
cover pages, is 7 for (P1+P2) and 10 for (P3+P4+P5), with character font size
12, in LATEX or another word processor.
• Pages must be numbered, with the numbers preferably placed in the centre of
the bottom line.
Software to deliver
Besides the report, you must deliver the following developed code and data:
• All the software to solve P1 and P2;

• The data files for the identification experiments openloop_data_1.mat and
openloop_data_2.mat;
• The identified model data file singleheater_model.mat;
• The mpc_solve function developed to solve the MPC optimization problem with
quadprog;
• The modified TCLab_simulation.m script configured to perform question 7 of
P4;
• The script TCLab_closedloop.m and data file closedloop_data_1.mat for the
closed-loop experiment of P5 in the real plant.

P1 – Basics on constrained optimization
Motivation
The key issue in MPC is that it is based on on-line minimization of a cost function to
obtain the value of the control variable to apply to the plant. In this way, MPC allows
to enforce constraints on several variables, a feature that is unique among on-line
design techniques (optimal control based on Pontryagin’s Principle may also embed
constraints, but this design technique is an off-line one).
Hence, it is important to get acquainted with algorithms and software for function
minimization. This work provides an exercise on the use of popular function
minimization Matlab tools. Although basic, this knowledge is very useful in writing our
own tools for MPC, as well as in using MPC toolboxes.
Objectives
Get acquainted to the basics of numerical function minimization using Matlab

Optimization Toolbox through the solution of simple two-dimensional problems,
namely
• fminunc, for unconstrained function minimization

• fmincon, for constrained function minimization
In two-dimensional problems the function to minimize can be represented graphically

either by using the Matlab functions
• surf, that draws a 3-dimensional plot of the function

• contour, that plots the level curves (the level curves of a function 𝑦 = 𝑓(𝑥1 , 𝑥2 )
are the curves defined by 𝑓(𝑥1 , 𝑥2 ) = 𝐶, for different values of the constant 𝐶).
Both surf and contour plot the function in a grid of points that can be defined using the
function
• meshgrid

In addition to the example provided below, it is advisable that you read the help and
doc (in the Matlab command window type, for instance, help fmincon or doc fmincon)
to get more information about these functions, as well as the syntax of their use.
A useful example
This example illustrates how to use the minimization software described above to
compute unconstrained and constrained minima for a function. This problem is a trivial
one, the aim being just to illustrate the use of the Matlab software.
The function to minimize is
𝑓(𝑥1 , 𝑥2 ) = (𝑥1 − 1)2 + (𝑥2 − 1)2 .
It is obvious that this function has an unconstrained minimum for 𝑥 = (1, 1).
This unconstrained minimization problem is simply formulated as
min 𝑓(𝑥)
𝑥
To use the Matlab solver fminunc to solve this problem, start by creating a Matlab
function script with the name of the function to minimize. For instance, calling the
function BasicFunction, you must create the script:
function f = BasicFunction(x)
f = (x(1)-1)^2 + (x(2)-1)^2;
end
To solve numerically the unconstrained optimization problem, define the initial guess
of the optimal value of 𝑥 (called x0 below) and run the following code
options = optimoptions('fminunc','Algorithm','quasi-newton');
xopt=fminunc(@BasicFunction,x0,options)
The first line defines the options of the solver (see Matlab help) and the second line
calls fminunc to actually solve the problem. The name of the script that contains the
definition of the function to minimize is given as the first argument of fminunc by the
name of the script preceded by the symbol @.

x1
Forbidden
1 region
Allowed
region
1 x2
Figure P1-1. Basic optimization example: definition of a linear inequality constraint.
Consider now the introduction of an inequality constraint, such as
𝑥2 ≤ 1 − 𝑥1 ,
that corresponds to the allowed values below the straight line in figure P1-1. This
constraint can be written in the form
𝐴𝑥 ≤ 𝐵, with 𝐴 = [1 1] and 𝐵 = 1.
This form of the constraint is the standard one to be specified in the Matlab solver
software function fmincon.
With the above constraint, the optimization problem is formulated as
min 𝑓(𝑥)
𝑥
𝑠𝑢𝑏𝑗𝑒𝑐𝑡 𝑡𝑜 𝐴𝑥 ≤ 𝐵
The numerical solution of this constrained problem using fmincon is coded in Matlab
by
Ac=[1 1];
Bc=1;
xoptconstr=fmincon(@BasicFunction,x0,Ac,Bc)

The first two lines define the constraint. The rest is similar to the use of fminunc
described above. Of course, one must also define the initial guess of the optimal value
of 𝑥 (called x0).
 Work to perform
0) The Matlab scripts BasicFunction.m and ProbBasic.m are available in Fénix.

These scripts allow to solve the two unconstrained, and constrained,
minimization problems described above. They also plot a 3-D view of the
function as well as its level curves. Study and run the script ProbBasic.m . You
can use it as prototype to the other activities in this set of questions. In this
question you don’t have to write anything on your report.
1) Compute the unconstrained minimum of the Rosenbrock function, given by
𝑓(𝑥) = 100(𝑥2 − 𝑥12 )2 + (1 − 𝑥1 )2 .
Compute also the minimum that satisfies the constraint
𝑥1 ≤ 0.5.
In both cases take as the initial estimate of the minimum
−1
𝑥0 = [ ].
1
In the same plot, show the level lines of the function, the initial estimate
(marked with “o”), the unconstrained minimum (marked with “x”), and the
constrained minimum (marked with “*”). In another figure, plot the function.

P2 – Basics on receding horizon control
Motivation
The receding horizon strategy is the key component of model predictive control (MPC).
It consists of designing the control action such as to optimize it over an extended
horizon, apply to the plant only the first control move of the optimized control
sequence, and then repeated the process at the next discrete time instant, starting
from the new state. It is the receding horizon strategy that allows to transform an
open-loop optimal control action into a feedback control action.
When enlarging the optimization horizon, the initial control move approximates the
one that would be obtained with an infinite horizon, which is stabilizing.
Objectives
The objective of the work to be done here is to compute the state feedback gain yield
by the receding horizon strategy in a problem with linear dynamics, quadratic cost and
without constraints, and to compare it with the infinite horizon gain computed using
the solution of an equivalent linear quadratic problem with infinite horizon, that
involves the solution of an algebraic Riccati equation.
1 Cost optimization
Infinite horizon LQ problem
Consider the infinite horizon linear quadratic (LQ) optimal control problem that
consists of minimizing the infinite horizon unconstrained quadratic cost
𝐽𝐿𝑄 (𝑢) = ∑ 𝑥 𝑇 (𝑡)𝑄𝑥(𝑡) + 𝑢𝑇 (𝑡)𝑅𝑢(𝑡)

𝑡=1
subject to
𝑥(𝑡 + 1) = 𝐴𝑥(𝑡) + 𝐵𝑢(𝑡).
The state 𝑥 has dimension 𝑛 and, in general, 𝑢 has dimension 𝑚.
1 This symbol means: study this part carefully at home.

The solution to this problem is obtained by solving, in order to the positive definite
𝑛 × 𝑛 matrix 𝑆, the discrete-time algebraic Riccati equation (ARE)
𝐴𝑇 𝑆𝐴 − 𝑆 − 𝐴𝑇 𝑆𝐵(𝐵 𝑇 𝑆𝐵 + 𝑅)−1 𝐵𝑇 𝑆𝐴 + 𝑄 = 0,
and then compute the optimal LQ gain by
𝐾𝐿𝑄 = (𝐵 𝑇 𝑆𝐵 + 𝑅)−1 𝐵𝑇 𝑆𝐴.
The optimal LQ control is given by the linear state feedback
𝑢(𝑡) = −𝐾𝐿𝑄 𝑥(𝑡).
The optimal LQ control exists if:
1. (𝐴, 𝐵) is stabilizable
2. 𝑅 ≻ 0 (read ≻ positive definite) and 𝑄 ≽ 0 (read ≽ positive semi-definite)
3. (𝑄, 𝐴) has no unobservable mode on the unit circle. If 𝑄 = 𝐶 𝑇 𝐶, with 𝐶 a
vector, this condition is verified if (𝐶, 𝐴) is observable.
Numerically, the discrete time LQ optimal gain can be computed using the MATLAB
Control Systems Toolbox function dlqr. The syntax to call this function is
[KLQ,S,lambda] = dlqr(A,B,Q,R)
where KLQ is the optimal gain, S is the positive definite solution of the Riccati equation
and lambda is the vector of eigenvalues of the closed-loop system dynamics 𝐴 − 𝐵𝐾𝐿𝑄 .
Finite horizon, receding horizon control
Consider now the finite horizon unconstrained quadratic cost, defined over a horizon
with 𝐻 steps, given by
𝐽𝑅𝐻 (𝑢; 𝑡) = ∑𝐻−1 𝑇 2

𝑖=0 𝑥 (𝑡 + 𝑖 + 1)𝑄𝑥(𝑡 + 𝑖 + 1) + 𝑅𝑢 (𝑡 + 𝑖).
It is assumed that the plant input is scalar, and hence 𝑅 is a positive scalar weight.

States influenced by future controls
over the horizon
t+0 t+1 t+2 t+H-1 t+H
Control sequence to optimize

Present time
u(t)=?
Figure P2-1. The control sequence to optimize in the RH cost and the corresponding
sequence of states that is influenced by it.
By selecting 𝑄 = 𝐶 𝑇 𝐶, where 𝐶 is a vector in the output equation
𝑦(𝑡) = 𝐶𝑥(𝑡),
the cost function penalizes the plant output and becomes
𝐽𝑅𝐻 (𝑢; 𝑡) = ∑𝐻−1 2 2

𝑖=0 𝑦 (𝑡 + 𝑖 + 1) + 𝑅𝑢 (𝑡 + 𝑖).
Figure P2-1 shows the time synchronization of the variables in this cost function. “We”
are at the present time 𝑡 and want to compute the control sample 𝑢(𝑡) to apply to the
plant. For that sake, we consider the outputs (or the states) between 𝑡 + 1 and 𝑡 + 𝐻.
Since the plant has a unit delay between applying a control sample and its influence on
the state and the output, this set of future outputs is influenced by the control inputs
between 𝑡 and 𝑡 + 𝐻 − 1.
The outcome of the optimization is, therefore, the sequence of future control inputs
between 𝑡 and 𝑡 + 𝐻 − 1. According to the receding horizon strategy, of this whole
control sequence, only the first element, 𝑢(𝑡) is applied to the plant, the whole
procedure being repeated at time 𝑡 + 1 (where the present time will be 𝑡 + 1, the
control sequence with respect to which the optimization is performed is 𝑢(𝑡 +
1), … , 𝑢(𝑡 + 𝐻), and the output sequence will be 𝑦(𝑡 + 2), … , 𝑢(𝑡 + 𝐻 + 1)) and so
on.
In order to simplify the notation, it is assumed, without loss of generality, that 𝑡 = 0.

Define the vectors
𝑦(1) 𝑢(0)
𝑌=[ ⋮ ] and 𝑈=[ ⋮ ].
𝑦(𝐻) 𝑢(𝐻 − 1)
With these definitions, the receding horizon cost can be written as
𝐽𝑅𝐻 = 𝑌 𝑇 𝑌 + 𝑅𝑈 𝑇 𝑈.
Using the state model, the vector of future output samples, 𝑌, can be expressed as a
function of the vector of future control samples, 𝑈, and of the initial state, thereby
resulting in a function of 𝑈. In order to follow this procedure, iterate the state model
to compute future values of the state, yielding (remember that time 0 actually
corresponds to present time 𝑡; 𝑏 is used instead of 𝐵 to emphasize that 𝑢 is scalar) the
following predictive models for the state
𝑥(1) = 𝐴𝑥(0) + 𝐵𝑢(0),
𝑥(2) = 𝐴2 𝑥(0) + 𝐴𝑏𝑢(0) + 𝑏𝑢(1),
𝑥(3) = 𝐴3 𝑥(0) + 𝐴2 𝑏𝑢(0) + 𝐴𝑏𝑢(1) + 𝑏𝑢(2),
or, in general,
𝑥(𝑖) = 𝐴𝑖 𝑥(0) + 𝐴𝑖−1 𝑏𝑢(0) + 𝐴𝑖−2 𝑏𝑢(1) + ⋯ + 𝑏𝑢(𝑖 − 1).
Predictive models for the output are obtained by multiplying by 𝐶, yielding
𝑦(𝑖) = 𝐶𝐴𝑖 𝑥(0) + 𝐶𝐴𝑖−1 𝑏𝑢(0) + 𝐶𝐴𝑖−2 𝑏𝑢(1) + ⋯ + 𝐶𝑏𝑢(𝑖 − 1).
This set of predictive models, for 𝑖 = 1, ⋯ , 𝐻, can be written in matrix form as
𝑦(1) 𝐶𝑏 0 0 ⋯ 0 𝑢(0) 𝐶𝐴
𝑦(2) 𝐶𝐴𝑏 𝐶𝑏 0 ⋯ 0 𝑢(1) 𝐶𝐴2
𝑦(3) = 𝐶𝐴2 𝑏 𝐶𝐴𝑏 𝐶𝑏 ⋱ ⋮ 𝑢(2) + 𝐶𝐴3 𝑥(0).
⋮ ⋮ ⋱ ⋱ ⋱ 0 ⋮ ⋮
[𝑦(𝐻)] [ 𝐶𝐴𝐻−1
𝑏 𝐶𝐴𝐻−2
𝑏 ⋯ 𝐶𝐴𝑏 ]
𝐶𝑏 [𝑢(𝐻 − 1)] [ 𝐶𝐴 𝐻]
Define the matrices

𝐶𝑏 0 0 ⋯ 0 𝐶𝐴
𝐶𝐴𝑏 𝐶𝑏 0 ⋯ 0 𝐶𝐴2
𝑊: = 𝐶𝐴2 𝑏 𝐶𝐴𝑏 𝐶𝑏 ⋱ ⋮ , Π ≔ 𝐶𝐴3 .
⋮ ⋱ ⋱ ⋱ 0 ⋮
[𝐶𝐴𝐻−1 𝑏 𝐶𝐴𝐻−2
𝑏 ⋯ 𝐶𝐴𝑏 𝐶𝑏] [𝐶𝐴𝐻 ]
With these definitions, the pencil of output predictive models can be written as
𝑌 = 𝑊𝑈 + Πx(0).
Inserting these predictive models in the cost
𝐽𝑅𝐻 = 𝑌 𝑇 𝑌 + 𝑅𝑈 𝑇 𝑈,
this same cost is written as
𝐽𝑅𝐻 = (𝑈 𝑇 𝑊 𝑇 + 𝑥 𝑇 (0)Π𝑇 )(𝑊𝑈 + Π𝑥(0)) + 𝑅𝑈 𝑇 𝑈
or
𝐽𝑅𝐻 = 𝑈 𝑇 (𝑊 𝑇 𝑊 + 𝑅𝐼)𝑈 + 2𝑥 𝑇 (0)Π𝑇 𝑊𝑈 + 𝑥 𝑇 (0)Π𝑇 Π𝑥(0),
which is a quadratic form in 𝑈. To simplify the notation, define the matrix
𝑀 ≔ 𝑊 𝑇 𝑊 + 𝑅𝐼.
With this definition, the finite horizon cost is written
𝐽𝑅𝐻 = 𝑈 𝑇 𝑀𝑈 + 2𝑥 𝑇 (0)Π𝑇 𝑊𝑈 + 𝑥 𝑇 (0)Π𝑇 Π𝑥(0).
The gradient of 𝐽𝑅𝐻 with respect to 𝑈, is2
∇𝑈 𝐽𝑅𝐻 = 2𝑈 𝑇 𝑀 + 2𝑥 𝑇 (0)Π𝑇 𝑊.
Equating the gradient to zero, yields the following equation, satisfied by the optimal
control sequence
𝑈 𝑇 𝑀 + 𝑥 𝑇 (0)Π𝑇 𝑊 = 0.
Transposing, and solving this equation with respect to 𝑈, yields the optimal control
sequence
2 The gradient with respect to 𝑥 of 𝑥 𝑇 𝐴𝑥 is 2𝑥 𝑇 𝐴; the gradient with respect to 𝑥 of 𝑏 𝑇 𝑥 is 𝑏 𝑇 .

𝑈 ∗ = −𝑀−1 𝑊 𝑇 Π𝑥(0).
According to the receding horizon strategy, only the first entry of this sequence is
applied to the plant. Recalling that 𝑥(0), actually, corresponds to 𝑥(𝑡), and that 𝑢(0)
corresponds to 𝑢(𝑡), the optimal receding horizon control is given by the state
feedback
𝑢(𝑡) = −𝐾𝑅𝐻 𝑥(𝑡),
with the optimal feedback gain given by
𝐾𝑅𝐻 = 𝑒1 𝑀−1 𝑊 𝑇 Π,
where
𝑒1 = [1 0 ⋯ 0]
is a row vector of dimension 𝐻.
 Work to perform
Consider the open-loop unstable 1st order plant
𝑥(𝑡 + 1) = 1,2 𝑥(𝑡) + 𝑢(𝑡).
For this plant, develop a Matlab script to answer the following questions:
1) Compute the optimal LQ state feedback gain;

2) Compute the optimal receding horizon gain for different values of the horizon
𝐻. Make a plot of the gain as a function of 𝐻 and superimpose it on a line that
corresponds to the optimal LQ gain.
3) Make a plot of the absolute value of the closed-loop eigenvalue (the
eigenvalue of 𝐴 − 𝑏𝐾) and compare it with the stability boundary.
4) Discuss the advantage of enlarging the horizon 𝐻 in relation to the above
example. Compare the RH and LQ gains.
5) Repeat the study for the open-loop stable 1st order plant
𝑥(𝑡 + 1) = 0,8 𝑥(𝑡) + 𝑢(𝑡).

In relation to stability, for which plant is it more advantageous to use enlarged
values of the horizon 𝐻?
In the above study, take 𝑄 = 1 and select values of 𝑅 that you consider that yield
representative results. Try with small and large values of 𝑅 (you must decide the actual
values. The maximum range of 𝐻 must also be decided by you.
P3 – Model identification
Motivation
MPC requires a dynamic model of the plant to control, in order to perform state
predictions. Although it is possible to build a very good model of the TCLab system
from first principles, via the laws of heat transfer and radiation, this may not be the
case for other systems, due to the high complexity of the system. Therefore, we
instead perform system identification, where a plant model is estimated from
experimental data.
Furthermore, in order to simplify the resulting MPC optimization problem, the model
to identify will be linear. This model will also be used later for the design of a Kalman
filter for state estimation, and for building a simulation environment that supports the
design of the controller and observer. Since we will only consider a single heater, the
identified model will be single-input/single-output (SISO).
Objectives
The objective for this part is the identification of a linear model for the TCLab system
from experimental data, considering only one heater and one temperature sensor.
Furthermore, we will validate the model using data obtained from a separate
experiment.
Linear model identification
We assume that the system can be described by a nonlinear difference equation such
as

𝑥(𝑘 + 1) = 𝑓(𝑥(𝑘), 𝑢(𝑘)),
where 𝑥 ∈ ℝ𝑛 is the vector of state variables, 𝑢 ∈ ℝ𝑚 is the vector of control

variables, 𝑓: ℝ𝑛 × ℝ𝑚 → ℝ𝑛 is a function that models the system dynamics, and 𝑘 is
the sample index (discrete time). Furthermore, we define the output equation
𝑦(𝑘) = 𝑔(𝑥(𝑘)),
where 𝑦 ∈ ℝ𝑝 is the vector of output variables and 𝑔: ℝ𝑛 → ℝ𝑝 is a function that

computes the output from the state variables.
For the TCLab system, the input variables are the two heater controls, and the output
variables are the measured temperatures of the two sensors. For this work, however,
we will only consider the control of one heater, in which case 𝑚 = 𝑝 = 1. We assume
that we do not know anything about the number of states 𝑛, their dynamics 𝑓, or
output function 𝑔. The system is essentially a “black box”, where we only have access
to its inputs and outputs. We must therefore obtain a model by exciting the system
sufficiently and appropriately.
Assume we want to control the output of the system around a constant equilibrium
value 𝑦̅, to which corresponds the equilibrium state 𝑥̅ , i.e.
𝑦̅ = 𝑔(𝑥̅ ),
and for which there exists a constant control 𝑢̅ such that the system is kept in
equilibrium, i.e.
𝑥̅ = 𝑓(𝑥̅ , 𝑢̅).
We want to find a linearised model of 𝑓 that holds well in the neighbourhood of 𝑥̅ and
𝑢̅, in order to design a controller that can regulate the plant output at 𝑦̅.
To this end, we can perform an experiment where the system is operated in this
neighbourhood, and decompose the commanded inputs 𝑢 and the measured outputs
𝑦 into the two components
𝑦 = 𝑦̅ + Δ𝑦,

𝑢 = 𝑢̅ + Δ𝑢,
where Δ𝑦 and Δ𝑢 are the deviations/increments from the equilibrium. We can then
estimate the 𝐴, 𝐵, 𝐶 and 𝐾𝑒 matrices that best explain the observed experimental
observations of the Δ𝑦 and Δ𝑢 dynamics with the discrete linear time-invariant model
Δ𝑥(𝑘 + 1) = 𝐴Δ𝑥(𝑘) + 𝐵Δ𝑢(𝑘) + 𝐾𝑒 𝑒(𝑘),
Δ𝑦(𝑘) = 𝐶Δ𝑥(𝑘) + 𝑒(𝑘),
where Δ𝑥 is the deviation from the steady-state 𝑥̅ such that
𝑥 = 𝑥̅ + ∆𝑥,
and 𝑒 is a Gaussian disturbance with a variance that is also estimated.
The controller will therefore consider this linear model for the incremental variables,
where we can achieve the desired output 𝑦̅ by driving Δ𝑦 to the origin, via a control
increment ∆𝑢 about 𝑢̅.
 Work to perform
The work to perform in this task consists of a series of experiments in which the plant
is operated in open-loop, with the control consisting of increments around the
equilibrium point.
0) Follow the TCLab installation and setup instructions from the guide available in
the Fénix course page. In the present work we will not use Simulink, so you may
skip the instructions related to that.
1) The MATLAB script TCLab_openloop.m allows for performing experiments in
open-loop, with a pre-defined control applied to the heaters. Modify this script
to use a sampling time of 5 seconds; this time interval will be the sampling time
of the discrete-time model to identify and of the MPC. Modify the input profile
of the first heater such that:

a. Initially, the input has a constant value for which the steady-state
temperature is between 40 and 50 degrees Celsius. After reaching this
equilibrium temperature, hold that input for some time.
b. Afterwards, raise and lower the input by small increments, e.g 5-10%,
allowing the temperature to reach equilibrium before changing the
input again (the transient duration may be in the order of 1000
seconds). Perform around five such input changes.
c. Change the name of the resulting .mat file to openloop_data_1.mat, to
be delivered with the final report. Show the plot of the experiment in
the report.
A simulated example of the desired profile is shown in Figure 1 for reference.
The idea is to drive the system to an equilibrium output 𝑦̅ for a certain input 𝑢̅,
and apply several step inputs Δ𝑢 around that equilibrium and observe the
output response Δ𝑦. We can then estimate the matrices for the linear
incremental dynamics that best fit the experimental data.
Figure 1 – Simulated example of experiment for model identification.
For the best results, try to minimize energy changes with the environment, e.g.
touching and moving the board, opening/closing windows, changing the room
air conditioning, etc.
2) Perform a second experiment that is shorter in duration (< 1000 seconds), but
which excites the system significantly more, with input changes of higher
amplitude and without necessarily waiting for the system to reach equilibrium

but keeping temperature approximately within the 40–50-degree Celsius
interval. Perform about 10 input changes. This second experiment will be used
to validate the identified model, given that this data will not be used at all in
the identification itself. A simulated example of the desired profile is shown in
Figure 2. Change the name of the resulting .mat file to openloop_data_2.mat,
and include it in the final delivery.
Figure 2 - Simulated example of experiment for model validation.
3) Run the script TCLab_identification.m. This script performs model identification

from the plant data, and has the following functionalities:
a. Loads the results of the first experiment and estimates the equilibrium
output 𝑦̅ (y_ss) corresponding to the initial input 𝑢̅ (u_ss) by
computing the mean value of 𝑦 during the initial equilibrium interval
(the time interval between the two red lines in Figure 1). Change the
k_ss_begin and k_ss_end variables to capture the correct initial and
final sample for that time interval.
b. Truncates the initial transient and computes the incremental output Δ𝑦
(Dy) and input Δ𝑢 (Du).
c. Uses the MATLAB ssest function to estimate the 𝐴, 𝐵, 𝐶 and 𝐾𝑒 matrices
that best fit the observations, as well as the variance of the disturbance
𝑒, with a default state dimension of 𝑛 = 1.
d. Estimates the initial state with the findstates function and propagates it
with the model, plotting the comparison between that result and the

experimental data. You must see a very close adherence of the output
propagated with the identified model to the experimental output data.
e. Loads the data from the second experiment, computes the incremental
variables, estimates the initial state and propagates it with the model,
again comparing it with the observations. Although the error between
the simulated output and the experimental data should now be greater,
the two data records must still match closely. The MSE error between
the two is printed in the command window.
Run the script multiple times for different values of the state dimension 𝑛 and
present a table with the validation errors for each case. Select the best order 𝑛
according to this metric and for that case present the plots for both tests in the
report. If two different state dimensions result in very similar MSE errors,
consider selecting the lower-order one. For the next parts, we will work with
the model generated for that state dimension, which is saved in the
singleheater_model.mat file.
P4 – MPC and Kalman filter design
Motivation
Given the model previously identified, we will now design a model predictive
controller for the system and test it in simulation with the same model, albeit with
some disturbances. We will initially consider an unconstrained MPC regulator, and
then incrementally complexify the controller to track a non-zero reference and include
control and output constraints. Given that the MPC requires knowledge of the state,
we will also design a state observer, namely a Kalman filter.
Objectives
The objectives of this part are to design a model predictive controller and a Kalman
filter, supported by a simulation environment, before applying them to the real system
in the next part.

Constrained MPC – practical aspects
We have seen in P2 that for an unconstrained problem with finite-horizon, quadratic

cost and linear dynamics, the optimal receding-horizon control is a simple state
feedback, where the feedback gain can be easily determined with algebraic
computations from the plant linear model and the cost weights. The real advantage of
using a receding-horizon method, however, is that we can directly include constraints
in the problem formulation, such as control limits and safety restrictions. On the other
hand, the inclusion of such constraints renders the computation of the solution more
complex, generally without a closed-form analytical solution being easily available, as it
is in the unconstrained case. In such cases, we need to use a numerical optimization
algorithm.
Given that the optimization problem must be solved in real-time, the computation
time of the optimizer is the main drawback of this control method. In this case, since
the optimization problems considered are quadratic programs, we will use the
MATLAB quadprog function, which implements an optimization algorithm specific for
that class of problem and is therefore more efficient at solving it than the previously
used fmincon (that must be used, for example, when the cost is not quadratic or the
dynamics nonlinear). Naturally, the bigger the problem – number of states, controls,
length of the prediction horizon, constraints – the higher the optimizer computational
load.
The quadprog function allows to minimize a quadratic program with the form
1 𝑇
minimize 𝑧 𝐹𝑧 + 𝑓 𝑇 𝑧
𝑧 2
subject to 𝐴𝑖𝑛𝑒𝑞 𝑥 ≤ 𝑏𝑖𝑛𝑒𝑞
𝐴𝑒𝑞 𝑧 = 𝑏𝑒𝑞
𝑙𝑏 ≤ 𝑧 ≤ 𝑢𝑏
where 𝑧 is the vector of optimization variables, 𝐹 is a positive semidefinite matrix, 𝑓 is

a vector, matrix 𝐴𝑖𝑛𝑒𝑞 and vector 𝑏𝑖𝑛𝑒𝑞 define the inequality constraints, matrix 𝐴𝑒𝑞

and vector 𝑏𝑒𝑞 define the equality constraints, and 𝑙𝑏 and 𝑢𝑏 define lower and upper-
bounds on 𝑧. The function may be called with the syntax
z = quadprog(F,f,Aineq,bineq,Aeq,beq,lb,ub),
where any unnecessary inputs may be replaced with an empty variable. To use this
function, we must therefore first write our problem in matrix form and determine the
parameters 𝐹, 𝑓, 𝐴𝑖𝑛𝑒𝑞 , 𝑏𝑖𝑛𝑒𝑞 , 𝐴𝑒𝑞 , 𝑏𝑒𝑞 , 𝑙𝑏 and 𝑢𝑏 that correspond to the specific
problem at hand.
Without loss of generality, we assume that we are at time 𝑘 = 0, and similarly to that
presented in P2 we consider the minimization of the cost function
𝐻−1
𝐽 = ∑ 𝑥̂ 𝑇 (𝑖 + 1)𝑄𝑥̂(𝑖 + 1) + 𝑅𝑢̂2 (𝑖),

𝑖=0
subject to the dynamics
𝑥̂(𝑖 + 1) = 𝐴𝑥̂(𝑖) + 𝐵𝑢̂(𝑖),
with the initial state 𝑥(0) and output equation
𝑦̂(𝑖) = 𝐶𝑥̂(𝑖).
The ‘hat’ notation in 𝑥̂, 𝑢̂ and 𝑦̂ denotes that these are virtual variables, in the sense
that they are predictions made within the MPC, and thus distinguish them from the
real variables 𝑥, 𝑢 and 𝑦 that have been or will be applied. From the minimization
yields the control profile (𝑢̂(0), … , 𝑢̂(𝐻 − 1)), out of which we only apply the first, i.e.
𝑢(0) = 𝑢̂(0),
and repeat the process in the next sampling interval.
Similarly to what was previously done in P2, we can define 𝑧 as the concatenation of all
the control variables

𝑢̂(0)
𝑧=𝑈=[ ⋮ ],
𝑢̂(𝐻 − 1)
which yields the cost terms
𝐹 = 2𝑀, 𝑓 = 2𝑥 𝑇 (0)Π𝑇 𝑊,
where 𝑥(0) is the initial state, and 𝑀, Π and 𝑊 are as previously defined in P2. If we
do not want to model any constraints, we can simply call the function with
U = quadprog(F,f),
and the solution will be the same as that obtained with the method shown in P2.
This is known as the dense formulation, since only the control variables are collocated
as optimization variables, yielding a problem that is smaller in size but where the cost
and constraint matrices are denser. This is in contrast to the sparse formulation, where
both the state and control variables are collocated, defining instead the optimization
vector
𝑋
𝑧 = [ ],
𝑈
where
𝑥̂(0)
𝑋 = [ ⋮ ].
𝑥̂(𝐻)
We can write the dynamics in matrix form as
𝑥̂(0) 0 ⋯ ⋯ 0 𝑥̂(0) 0 ⋯ 0 𝑢̂(0) 𝐼

𝑥̂(1) 𝐴 ⋯ 0 ⋮ 𝑥̂(1) 𝐵 ⋯ 0 𝑢̂(1) 0
[ ]=[ ][ ]+[ ][ ] + [ ] 𝑥(0),
⋮ ⋮ ⋱ ⋮ ⋮ ⋮ ⋮ ⋱ ⋮ ⋮ ⋮
𝑥̂(𝐻) 0 ⋯ 𝐴 0 𝑥̂(𝐻) 0 ⋯ 𝐵 𝑢̂(𝐻 − 1) 0
Defining the matrices
0 ⋯ ⋯ 0 0 ⋯ 0 𝐼
𝐴 ⋯ 0 ⋮ 𝐵 ⋯ 0 0
𝐴̃ ∶= [ ], 𝐵̃ ∶= [ ] , 𝐸 ∶= [ ],
⋮ ⋱ ⋮ ⋮ ⋮ ⋱ ⋮ ⋮
0 ⋯ 𝐴 0 0 ⋯ 𝐵 0

we can write
𝑋 = 𝐴̃𝑋 + 𝐵̃ 𝑈 + 𝐸𝑥(0).
In the sparse formulation, we must implicitly include the dynamics as an equality

constraint on 𝑧, yielding
𝐴𝑒𝑞 𝑧 = 𝑏𝑒𝑞 ,
where
𝐴𝑒𝑞 ∶= [𝐴̃ − 𝐼 𝐵̃ ], 𝑏𝑒𝑞 ∶= −𝐸𝑥̂(0).
Now for the cost function, we can write it in matrix form as
0 ⋯ ⋯ 0 𝑥̂(0)
⋮ 𝑄 ⋯ 0 𝑥̂(1)
𝐽(𝑋, 𝑈) = [𝑥̂ 𝑇 (0) 𝑥̂ 𝑇 (1) ⋯ 𝑥̂ 𝑇 (𝐻)] [ ][ ]
⋮ ⋮ ⋱ ⋮ ⋮
0 0 ⋯ 𝑄 𝑥̂(𝐻)
𝑅 ⋯ 0 𝑢̂(0)
+ [𝑢̂ 𝑇 (0) 𝑇 (𝐻
⋯ 𝑢̂ − 1)] [ ⋮ ⋱ ⋮][ ⋮ ].
0 ⋯ 𝑅 𝑢̂(𝐻 − 1)
Defining the augmented cost matrices
0 ⋯ ⋯ 0
⋮ 𝑄 ⋯ 0 𝑅 ⋯ 0
𝑄̃ ∶= [ ], 𝑅̃ ∶= [ ⋮ ⋱ ⋮]
⋮ ⋮ ⋱ ⋮
0 ⋯ 𝑅
0 0 ⋯ 𝑄
yields the quadprog cost matrix
𝑄̃ 0
𝐹 = 2[ ],
0 𝑅̃
and 𝑓 = 0.
The sparse and dense formulations are theoretically equivalent, however, depending
on the problem and on the optimization algorithm, one may be computationally faster
to solve than the other. Furthermore, it is easier in the sparse formulation to define

state or output constraints, since the former is directly available, which is not the case
for the dense formulation. Let us say we want to constrain the output as
𝐺𝑦̂(𝑖) ≤ 𝑔, 𝑖 = 1, … , 𝐻,
for some matrix 𝐺 and vector 𝑔, we can write the constraints in matrix form as
𝐺 ⋯ 0 𝑦̂(1) 𝑔
[⋮ ⋱ ⋮ ] [ ⋮ ] ≤ [ ⋮ ].
⏟0 ⋯ 𝐺 𝑦̂(𝐻) 𝑔
⏟
∶=𝐺̃ ∶=𝑔̃
For the sparse formulation, we can compute the output vector 𝑌 from 𝑋 with
0 ⋯ 0 𝑥̂(0)
𝑦̂(1)
𝐶 ⋯ 0 𝑥̂(1)
[ ⋮ ]=[ ][ ],
⋮ ⋱ ⋮ ⋮
𝑦̂(𝐻)
⏟0 ⋯ 𝐶 𝑥̂(𝐻)
∶=𝐶̃
and the constraints become
𝐺̃ 𝐶̃ 𝑋 ≤ 𝑔̃.
We can therefore implement these constraints with the function quadprog by defining
𝐴𝑖𝑛𝑒𝑞 = [𝐺̃ 𝐶̃ 0], 𝑏𝑖𝑛𝑒𝑞 = 𝑔̃.
For the dense formulation, the output vector is instead computed from 𝑈 and 𝑥(0)
with
𝑌 = 𝑊𝑈 + Π𝑥(0).
Therefore, the constraints become instead
𝐴𝑖𝑛𝑒𝑞 = 𝐺̃ 𝑊, 𝑏𝑖𝑛𝑒𝑞 = 𝑔̃ − 𝐺̃ Π𝑥(0).
The constraint
𝐺𝑦̂(𝑖) ≤ 𝑔, 𝑖 = 1, … , 𝐻

Is known as a hard constraint, because it must be strictly satisfied by the computed
solution. For practical reasons, it may be advantageous to instead consider the relaxed
constraint
𝐺𝑦̂(𝑖) ≤ 𝑔 + 𝜂𝑖 , 𝑖 = 1, … , 𝐻,
where the variables 𝜂𝑖 ≥ 0, 𝑖 = 1, … , 𝐻 are included as optimization variables, and

heavily penalized in the cost function with 𝛼𝜂𝑖2 where 𝛼 is very large compared to the
remaining terms. This constraint is known as a soft constraint, because although it is
strongly encouraged to be satisfied, it is not strictly required. In this case, the
optimization vector becomes
𝑋
𝑧 = [𝑈],
𝜂
where 𝜂 = [𝜂1𝑇 , … , 𝜂𝐻
𝑇 ]𝑇
. The quadprog cost function parameters 𝐹 and 𝑓 as well as the
constraint parameters 𝐴𝑒𝑞 , 𝑏𝑒𝑞 , 𝐴𝑖𝑛𝑒𝑞 , 𝑏𝑖𝑛𝑒𝑞 , 𝑙𝑏 and 𝑢𝑏 must be changed accordingly to
accommodate for the new optimization variables. In particular, the 𝑙𝑏 parameter must
ensure 𝜂𝑖 ≥ 0, 𝑖 = 1, … , 𝐻, in order to ensure the constraint can only be relaxed, not
tightened.
MPC regulator
Let us again consider the single-heater TCLab system. We want to design a model
predictive controller that drives the incremental output Δ𝑦 to 0, such that the absolute
output 𝑦 is driven to the equilibrium 𝑦̅. Hence, assuming without loss of generality that
we are at time 𝑘 = 0, we define the following quadratic cost function
𝐻−1
𝐽 = ∑ Δ𝑦̂ 2 (𝑖 + 1) + 𝑅Δ𝑢̂2 (𝑖),

𝑖=0
which we want to minimize given the system model previously identified
Δ𝑥̂(𝑖 + 1) = 𝐴Δ𝑥̂(𝑖) + 𝐵Δ𝑢̂(𝑖),

Δ𝑦̂(𝑖) = 𝐶Δ𝑥̂(𝑖),
for 𝑖 = 0, … , 𝐻 − 1, and given an initial state Δ𝑥(0).
Notice that we do not include the disturbance 𝑒 in the MPC model, since that is a
random variable that is obviously unknown a priori. Nonetheless, the MPC will reject
this disturbance by operating in closed loop, driving the system sufficiently close to the
origin despite it.
Tracking with feed-forward
Let us now assume we instead want to control the output to a reference 𝑟 = 𝑦̅ + Δ𝑟.
We can design a model predictive controller that tracks the incremental reference Δ𝑟,
rather than regulating Δ𝑦 to the origin. One possibility to achieve this is by penalizing
the reference tracking error in the cost function
𝐻−1
𝐽 = ∑(Δ𝑦̂(𝑖 + 1) − Δ𝑟)2 + 𝑅Δ𝑢̂2 (𝑖).

𝑖=0
This formulation, however, will have a steady-state tracking error unless the system
has integral action, due to the fact that tracking Δ𝑟 requires a non-zero control
increment Δ𝑢, and therefore the minimization of the tracking error is adversarial to
that of the control energy, and the result is something in-between. To overcome this,
we can determine the control increment Δ𝑢̅ that maintains the output Δ𝑦 at the value
Δ𝑟, and instead formulate the cost function
𝐻−1
𝐽 = ∑(Δ𝑦̂(𝑖 + 1) − Δ𝑟)2 + 𝑅(Δ𝑢̂(𝑖) − Δ𝑢̅)2 ,

𝑖=0
in which case the control will not be penalized for maintaining Δ𝑟. The control Δ𝑢̅ can
be determined by solving the steady-state equations
Δ𝑥̅ = 𝐴Δ𝑥̅ + 𝐵Δ𝑢̅,
Δ𝑟 = 𝐶Δ𝑥̅ ,

for Δ𝑥̅ and Δ𝑢̅. This tracking method is essentially a feedforward control, where the
MPC computes a control increment about Δ𝑢̅ that regulates the output to the
equilibrium Δ𝑟. This is in addition to the feedforward term 𝑢̅ that we are already using,
where the MPC computes a control increment Δ𝑢 about 𝑢̅.
We can obtain an equivalent result using the same MPC regulator as before, by first
applying the change of variables
δ𝑥̂ = Δ𝑥̂ − Δ𝑥̅ ,
δ𝑥̂ = Δ𝑥̂ − Δ𝑟,
δ𝑥̂ = Δ𝑥̂ − Δ𝑢̅,
thus considering δ𝑥̂ and δ𝑢̂ as the new optimization variables instead (or δ𝑢̂ only in
case of the dense formulation), which are subject to the model
δ𝑥̂(𝑖 + 1) = 𝐴δ𝑥̂(𝑖) + 𝐵δ𝑢̂(𝑖),
δ𝑦̂(𝑖) = 𝐶δ𝑥̂(𝑖).
Naturally, if control and output constraints are present, they must also consider this
change of variables.
Zero-offset tracking
The previous tracking formulation only guarantees no steady-state error if the real
plant matches the MPC model. In practice, however, this will not be the case, e.g. due
to changes in the ambient temperature. We can overcome this problem by admitting
that the system is affected by an input disturbance 𝑑, which results in the dynamics
Δ𝑥(𝑘 + 1) = 𝐴Δ𝑥(𝑘) + 𝐵(Δ𝑢(𝑘) + 𝑑(𝑘)) + 𝐾𝑒 Δ𝑒(𝑘).
This changes the steady-state equations to
Δ𝑥̅ = 𝐴Δ𝑥̅ + 𝐵(Δ𝑢̅ + 𝑑̅ ),
Δ𝑟 = 𝐶Δ𝑥̅ ,

which must still be solved for Δ𝑥̅ and Δ𝑢̅, and which requires knowledge of the steady-
state disturbance 𝑑̅. To avoid this assumption, we will estimate the disturbance 𝑑
together with the state 𝑥 in the state observer, as shown in the next section, yielding
the estimate 𝑑̂ . The best choice for 𝑑̅ will then be
𝑑̅ = 𝑑̂ .
Once we have computed Δ𝑥̅ and Δ𝑢̅, we can simply use the same change-of-variable
strategy as described in the previous section.
Kalman filter
The model predictive controller requires an initial condition on the state Δ𝑥 in order to
perform the predictions, whereas we only have measurements from the output 𝑦.
Furthermore, it requires an estimate of the disturbance 𝑑. Therefore, we need to
implement a state observer, which uses the output measurements, command inputs,
and a model of the system to estimate the state, which, given sufficient time, must
converge to the real state. For this purpose, we will use a Kalman filter, which designs
a linear state observer in an optimal way.
We assume the model
Δ𝑥(𝑘 + 1) = 𝐴Δ𝑥(𝑘) + 𝐵Δ𝑢(𝑘) + 𝐵𝑑(𝑘) + 𝐾𝑒 Δ𝑒(𝑘),
Δ𝑦(𝑘) = 𝐶Δ𝑥(𝑘) + 𝑒(𝑘),
which includes the input disturbance 𝑑 introduced in the previous section. The
Gaussian disturbances 𝑒 and 𝐾𝑒 𝑒 have the known covariances
𝑇
𝐸 [(𝐾𝑒 𝑒(𝑘))(𝐾𝑒 𝑒(𝑘)) ] = 𝑄𝐸 ,
𝐸[𝑒(𝑘)𝑇 𝑒(𝑘)] = 𝑅𝐸 .
Although in this case we have identified these covariances from experimental data,
typically this is very hard to do for more complex systems, and the covariances only
interpreted as tuning knobs for the filter, as discussed later.

In order to estimate the disturbance 𝑑 together with the state 𝑥, we will augment the
state with the disturbance, yielding the dynamics
Δ𝑥(𝑘 + 1) 𝐴 𝐵 Δ𝑥(𝑘) 𝐵
[ ]=[ ][ ] + [ ] Δ𝑢(𝑘),
𝑑(𝑘 + 1) 0 1 𝑑(𝑘) 0
and output equation
Δ𝑥(𝑘)
𝑦(𝑘) = [𝐶 0] [ ]
𝑑(𝑘)
Defining the augmented state
Δ𝑥
𝑥𝑑 = [ ],
𝑑
the above model can be written as
𝑥𝑑 (𝑘 + 1) = 𝐴𝑑 𝑥𝑑 (𝑘) + 𝐵𝑑 Δ𝑢(𝑘),
𝑦(𝑘) = 𝐶𝑑 𝑥𝑑 (𝑘),
where
𝐴 𝐵 𝐵
𝐴𝑑 = [ ], 𝐵𝑑 = [ ], 𝐶𝑑 = [𝐶 0].
0 1 0
Consider also the augmented state covariance matrix
𝑄𝐸 0
𝑄𝐸 𝑑 = [ ],
0 𝛿𝐸
where 𝛿𝐸 is a tuning parameter.
The estimation is made in two steps: prediction and correction. At the beginning of the
sampling interval 𝑘, we first propagate the previous state estimate with the model and
the previously commanded control action, i.e. we compute
𝑥̂𝑑 − (𝑘) = 𝐴𝑑 𝑥̂𝑑 (𝑘 − 1) + 𝐵𝑑 𝑢(𝑘 − 1),
which is the prediction step. Afterwards, we correct this estimate using the latest
measurement 𝑦(𝑘), thus yielding the final estimate 𝑥̂(𝑘)

𝑥̂𝑑 (𝑘) = 𝑥̂𝑑 − (𝑘) + 𝐿(𝑦(𝑘) − 𝐶𝑑 𝑥̂𝑑 − (𝑘)),
where 𝐿 is the observer gain. Note that here we use the ‘hat’ notation to denote an
estimate; this must not be confused with the ‘hat’ notation in the context of MPC,
which instead denotes a prediction.
The Kalman filter selects the observer gain that minimizes the variance of the
estimation error
𝐽𝐸 = 𝐸 [∑‖𝑥𝑑 (𝑘) − 𝑥̂𝑑 (𝑘)‖2 ],

𝑘=0
given the disturbance covariances 𝑄𝐸 𝑒 and 𝑅𝐸 . In practice, we can look at these

covariances as tuning knobs for the filter, where decreasing the state disturbance
covariance 𝑄𝐸 relative to the output disturbance variance 𝑅𝐸 increases our trust in the
model and thus the gain 𝐿 will be smaller, while increasing it instead increases our
trust in the measurements, resulting in a greater gain.
We can compute the Kalman gain 𝐿 with the MATLAB dlqe function, with the syntax
L = dlqe(Ad,eye(n+1),Cd,QEd,RE),
where n is the state dimension 𝑛.
 Work to perform
0) In this part you will use the TCLab_simulation.m MATLAB script provided, which
simulates the TCLab system with the model identified in the previous part and
an open-loop control profile, starting at the equilibrium temperature for 𝑢 = 0
(ambient temperature). Run the script to test that it runs properly. You do not
need to put anything in the report for this question.
Tip: During design, it can be advantageous to momentarily disable the random
disturbance 𝑒 by setting the standard deviation variable e_std to 0, in order to
test the controller and observer in ideal conditions.
1) Develop a MATLAB function named mpc_solve, in order to solve the
unconstrained MPC regulator for a given initial condition using quadprog, using
either the sparse or dense formulations. Note that a lot of the code developed

for solving the receding-horizon problem in P2 can likely be reused. As a sanity
check, you may want to apply it to one of the problems from P2 and ensure
that the solution obtained with quadprog is the same as before. As a
suggestion, your function may have the syntax
u0 = mpc_solve(x0,H,R,A,B,C),
where u0 is first control action, that is actually applied to the plant. This
function will be modified in the next sections, and the final version is to be
delivered with the report.
2) Implement in the TCLab_simulation script a closed-loop unconstrained model
predictive controller that regulates Δ𝑦 to 0 using the mpc_solve function.
Discuss the effect of changing the prediction horizon 𝐻 and control weight 𝑅,
and select a suitable horizon 𝐻, to be fixed onwards unless otherwise
necessary.
3) Choose a control weight 𝑅 such that the resulting control 𝑢 exceeds the
interval [0,100] %. Obviously, this control cannot be applied to the real plant
and will be saturated. Afterwards, implicitly include these control limits as
constraints in the model predictive controller, by modifying the mpc_solve
function appropriately. Show a plot of the result.
4) Track a reference Δ𝑟 of about 5 degrees with the feedforward tracking scheme,
by applying a suitable change of variables and using the same MPC regulator as
before. Plot and show the results. Then, perturb the simulation model constant
c1 with a 10% increase, which simulates a change in the ambient temperature.
Plot the results and discuss the observations. Leave this disturbance on for
questions 6 and 7 (turn it off for question 5).
5) Include a safety constraint for a maximum temperature 𝑦 of 55 degrees Celsius.
Define a reference above this value. Check whether the MPC optimization is
always successful, i.e. that the quadprog output exitflag is 1. If not, explain
why, and change the constraint to a soft constraint. Plot the results.
6) Design a Kalman filter with the dlqe function for the augmented state 𝑥𝑒 . Use
the identified disturbance variance e_var to initially select the 𝑄𝐸 and 𝑅𝐸
parameters. Implement the filter in the simulator, and test it without the MPC
in the loop, instead using the open-loop profile 𝑢 = 𝑢̅. Initialize the state

estimate by adding an error to the “real” initial state and tune the filter
parameters if necessary. Show a plot with comparison between the estimated
output 𝑦̂ and measured output 𝑦, and with the estimated disturbance 𝑑.
Discuss the effect of changing the 𝛿𝐸 parameter.
7) Add the MPC in the loop again, where it must now receive the state estimate
Δ𝑥̂ rather than the “real” state Δ𝑥. Compute the feedforward terms Δ𝑥̅ and Δ𝑢̅
using the disturbance estimate 𝑑̂ . Perform a simulation where the reference 𝑟
commutes to four different values: 50 degrees, 40 degrees, 60 degrees (which
cannot be reached due to the safety constraint) and 45 degrees. Allow
equilibrium to be maintained for a few minutes before changing the reference
again, to show the steady-state behaviour of the closed-loop system. Show a
plot and discuss the results. Discuss how this tracking scheme with the
estimation of 𝑑 and its compensation in the MPC resembles a controller with
integral action, e.g. PI controller.
One common rule of thumb when designing a MPC, is that the computation time is
less than 10% of the sampling time, since otherwise it introduces a non-negligible input
delay that may negatively affect performance. Although, most likely, this is not a
problem in this case, you must ensure that this rule is satisfied for your controller and,
if necessary, reduce the prediction horizon, choose a lower-order state dimension, or
select the dense/sparse formulation, in order to reduce the computational load.
P5 – Application to the real system
Motivation
In P3, a linearised model for the TCLab system was identified and, in P4, a model
predictive controller and Kalman filter were designed for it. The controller obtained by
coupling these two elements will now be applied to the real thermal plant.
Objectives
The objective is to apply the previously developed controller and observer to the real
system.

 Work to perform
In this section we will work on the script TCLab_closedloop, which is essentially a copy
of TCLab_openloop, with some placeholder suggestions of where to place the Kalman
filter and MPC. It should be straightforward to adapt your code from TCLab_simulation
to this script.
1) Implement the previously designed Kalman filter and MPC in the

TCLab_closedloop script. It is recommended that you implement the Kalman
filter first and test it with an open-loop profile, only afterwards adding the
MPC. Perform an experiment with the same specifications as Question 7) of P4
and save the resulting .mat file as closedloop_data_1.mat, to be delivered with
the report. Show a plot and discuss the results.
References
[1] Rawlings, James Blake, David Q. Mayne, and Moritz Diehl. Model predictive control:
theory, computation, and design. Vol. 2. Madison, WI: Nob Hill Publishing, 2017.
(Available at https://sites.engineering.ucsb.edu/~jbraw/mpc/)
[2] A. Botelho, B. Parreira, P. N. Rosa and J. M. Lemos. Predictive Control for Spacecraft
Rendezvous. Springer, 2021.

Appendix 1 – Good practice rules to write the report
The lab grades are given based on the quality evaluation of the reports, according to
the following good practice rules, as well as on the participation of the students on the
lab sessions and report elaboration.
1. Adequacy of the answers to the questions in the work guide

a) The answers must fully respond to the intended
b) Conciseness. Use concise text, but that fully expresses the answer
c) Clarity. The text and its articulation with the figures included in the report
must be clear. If variables not defined in the statement are used, they must
be defined in the report.
d) Technical correction. There must be neither technical errors in the answers,
nor errors in the use of concepts.
e) Critical analysis of results. Do the results make sense given what is known
and/or when compared to each other? What is the interpretation of the
results?
2. Imagination and originality in the answers and in the way they are presented.
3. Correction of the text and elegance of the style.
a) There must be no spelling or syntax errors.
b) Preferably use short and clear sentences
c) Well articulated text.
d) Proper punctuation, including correct punctuation of equations.
4. Readable and informative graphics and plots

a) Identify the variables that correspond to each axis and, if appropriate, the
units in which they are expressed.
b) Proper axes font (not too small or too large; type 14 is suggested).
c) Use appropriate scales, in order to show the important part of the evolution
of the curves and/or the global evolution of the variables, depending on the
situation.
d) Do not hide important parts of curves with labels.

e) Do not use a range of variables that is equal to the range of data. In
particular, if the variable to plot has periods during which it is constant, do
not select the range of the axis to coincide with this constant. For instance,
when plotting a square wave between 0 and 1, do note use vertical axes
between 0 and 1, but include a small margin, say, axis between -0.1 and 1.1.
f) Size of the graphics in the document suitable for what you want to show.
g) Concise graph legends that elucidate the main aspect of the graph.
5. Clear, well-structured and documented programs
1. Correction of programs.
2. Good structure of programs.
3. Adequate and well distributed comments throughout the program.
4. Good program intelligibility.
5. Adequate and short names of variables, explanation of variables with
comments.
6. Avoid using i and j as variables in your program because they are preassigned
to the unit imaginary number (type i in the Matlab command window and see
the answer; repeat for j). You can use ii and jj as counters, for instance.
– End of document –

ECPD-2023-2024 Lab Guide

Uploaded by

Copyright:

Available Formats

ECPD-2023-2024 Lab Guide

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

ECPD-2023-2024 Lab Guide

Uploaded by

Copyright:

Available Formats

Master in Electrical and Computer Engineering

Master in Aerospace Engineering

Estimação e Controlo Preditivo Distribuído

Distributed Predictive Control and Estimation

Predictive Control of a Thermal Plant

Picture source: https://apmonitor.com/heat.htm

João Miranda Lemos

Departamento de Engenharia Eletrotécnica e de Computadores

Área Científica de Sistemas, Decisão e Controlo

IST – Estimação e Controlo Preditivo Distribuído – TCLab work 2023/2024 Page 1

• Receding horizon control

Laboratory guide organization

Hereafter, Pi identifies part i of the work to be done, as listed above.

IST – Estimação e Controlo Preditivo Distribuído – TCLab work 2023/2024 Page 2

A typical MATLAB program to do temperature control consists therefore of a time

In order not to damage the kit, do not apply high values of

IST – Estimação e Controlo Preditivo Distribuído – TCLab work 2023/2024 Page 3

• MATLAB, available through the IST campus license;

The information required to perform this work is described below.

L=0.25 G1 + 0.15 G2 + 0.1 G3 + 0.3 G4 + 0.2 G5

where Gi is the grade obtained in the part Pi of the work.

Participation in laboratory classes

Format of the report to deliver

IST – Estimação e Controlo Preditivo Distribuído – TCLab work 2023/2024 Page 4

• All the software to solve P1 and P2;

IST – Estimação e Controlo Preditivo Distribuído – TCLab work 2023/2024 Page 5

Get acquainted to the basics of numerical function minimization using Matlab

• fminunc, for unconstrained function minimization

In two-dimensional problems the function to minimize can be represented graphically

• surf, that draws a 3-dimensional plot of the function

IST – Estimação e Controlo Preditivo Distribuído – TCLab work 2023/2024 Page 6

The function to minimize is

𝑓(𝑥1 , 𝑥2 ) = (𝑥1 − 1)2 + (𝑥2 − 1)2 .

This unconstrained minimization problem is simply formulated as

IST – Estimação e Controlo Preditivo Distribuído – TCLab work 2023/2024 Page 7

Figure P1-1. Basic optimization example: definition of a linear inequality constraint.

Consider now the introduction of an inequality constraint, such as

With the above constraint, the optimization problem is formulated as

IST – Estimação e Controlo Preditivo Distribuído – TCLab work 2023/2024 Page 8

0) The Matlab scripts BasicFunction.m and ProbBasic.m are available in Fénix.

IST – Estimação e Controlo Preditivo Distribuído – TCLab work 2023/2024 Page 9

Infinite horizon LQ problem

𝐽𝐿𝑄 (𝑢) = ∑ 𝑥 𝑇 (𝑡)𝑄𝑥(𝑡) + 𝑢𝑇 (𝑡)𝑅𝑢(𝑡)

𝑥(𝑡 + 1) = 𝐴𝑥(𝑡) + 𝐵𝑢(𝑡).

The state 𝑥 has dimension 𝑛 and, in general, 𝑢 has dimension 𝑚.

1 This symbol means: study this part carefully at home.

IST – Estimação e Controlo Preditivo Distribuído – TCLab work 2023/2024 Page 10

and then compute the optimal LQ gain by

𝐾𝐿𝑄 = (𝐵 𝑇 𝑆𝐵 + 𝑅)−1 𝐵𝑇 𝑆𝐴.

The optimal LQ control is given by the linear state feedback

𝑢(𝑡) = −𝐾𝐿𝑄 𝑥(𝑡).

The optimal LQ control exists if:

Finite horizon, receding horizon control

𝐽𝑅𝐻 (𝑢; 𝑡) = ∑𝐻−1 𝑇 2

IST – Estimação e Controlo Preditivo Distribuído – TCLab work 2023/2024 Page 11

t+0 t+1 t+2 t+H-1 t+H

Control sequence to optimize

By selecting 𝑄 = 𝐶 𝑇 𝐶, where 𝐶 is a vector in the output equation

the cost function penalizes the plant output and becomes