1

Stability analysis of Model Predictive Controllers
using Mixed Integer Linear Programming

arXiv:1604.00863v1 [math.OC] 4 Apr 2016

Daniel Simon and Johan Löfberg

Abstract—It is a well known fact that finite time optimal
controllers, such as MPC does not necessarily result in closed
loop stable systems. Within the MPC community it is common
practice to add a final state constraint and/or a final state
penalty in order to obtain guaranteed stability. However, for
more advanced controller structures it can be difficult to show
stability using these techniques. Additionally in some cases the
final state constraint set consists of so many inequalities that the
complexity of the MPC problem is too big for use in certain fast
and time critical applications. In this paper we instead focus on
deriving a tool for a-postiori analysis of the closed loop stability
for linear systems controlled with MPC controllers. We formulate
an optimisation problem that gives a sufficient condition for
stability of the closed loop system and we show that the problem
can be written as a Mixed Integer Linear Programming Problem
(MILP).

I. I NTRODUCTION
A linear Model Predictive Controller (MPC) solves, online
in each sample instant, a finite time horizon optimal control
problem of the form
Vk∗ =

min.

xk+i ,uk+i

N
−1
X

`(xk+i , uk+i ) + Ψ(xk+N )

(1a)

i=0

s.t. xk+i+1 = Axk+i + Buk+i
Exk+i ≤ f

(1b)

i = 1, . . . , N − 1

(1c)

T xk+N ≤ t
Guk+i ≤ h i = 1, . . . , N − 1

(1d)
(1e)

and implements the optimal solution, u∗k , as the input to the
system in a receding horizon fashion.
While the system dynamics (1b), the state and control constraints (1c) and (1e) reflects the system and the requirements
on it, the final state constraint (1d) and penalty, Ψ(xk+N )
are in a sense artificial and added merely to ensure recursive
feasibility and stability.
Stability of the MPC control law (1) is most often proven
a-priori by showing that the objective function (1a) is a valid
Lyapunov function for the closed loop system by designing
Ψ(xk+N ) and the terminal set such that
∗
Vk+1
− Vk∗

(2)

is guaranteed to be less that zero. See, e.g., [1] for details.
The authors are with the Departement of Electrical Engineering, Linköping
University, SE-581 83 Linköping, Sweden, e-mail: dansi@isy.liu.se, johanl@isy.liu.se
This work has been submitted to the 55th , IEEE Conference on Decision
and Control 2016

However there exist many MPC formulations, such as move
blocking [2] or soft constraints [3] (i.e., slack), for which it
can be difficult to show stability using this standard framework
presented in [1]. Even if it allows for the possibility to
guarantee stability, the addition of the terminal constraints can
add a significant complexity to the original problem which
might be unnecessary and limit the applicability of MPC
within certain fields. In reference tracking MPC the terminal
set can become very complex and some times not even finitely
determined, see e.g., [4] or [5]. Therefore it can often be
beneficial to analyse and verify the stability of a certain design
rather than building in the stability by adding extra constraints.
The problem of analysing stability of optimisation based
controllers is no new field. In [6] the authors derive a stability
test based on bounds of the cost function. The author of
[7] uses the KKT conditions for the MPC controller (1) to
derive, using the S-procedure, an LMI which gives a sufficient
condition for stability. Several papers have followed up on this
idea such as, e.g., [8]–[12]. These papers have extended the
idea to hold for more general cases or improved the complexity
of the resulting LMI. E.g., in [12] the authors extends the ideas
in [7] to hold for more general systems and show that given
that the system and its constraints are polynomial, the stability
can be analysed using sum-of-squares programming. In [9] the
authors derive an LMI of much lower dimension than that of
[7] which improve on the complexity of the problem.
In this paper we also exploit the KKT conditions of the MPC
problem (1) to analyse the difference of the value function
(2). But instead of deriving an LMI condition we formulate
the problem as an indefinite quadratic bilevel optimisation
problem. We then show that this problem can be rewritten
as a Mixed Integer Linear Program (MILP), which can have
major computational advantages compared to using the LMI
formulation.
II. T HE MILP S TABILITY TEST
In this section we will derive the proposed stability test.
We will formulate an optimisation problem that uses the MPC
problem’s objective function as a candidate Lyapunov function
and then minimises the difference of the Lyapunov function
between two consecutive time steps. Since the optimisation
problem only test the validity of a certain Lyapunov function
candidate the formulated problem can only verify stability, not
prove instability. Therefore it is a sufficient but not necessary
condition for stability. Nevertheless, in situations where the
value function is non-decreasing at some point, it might be an
indication of problems in the design and furtrher investigations
and simulations needs to be made.

2

The resulting optimisation problem is an indefinite quadratic
bilevel optimisation problem which are very difficult to solve.
We will finally show how this indefinite problem can be
rewritten as a mixed integer linear programming problem.
If the objective function (1a) is a quadratic or LPrepresentable function in x and u we can rewrite the MPC
problem (1), at time step k, in a more compact form as
minimize Vk

(3)

Uk

s.t. Exk + F Uk ≤ b

conditions
HUk + Gxk + F T λk = 0

(5a)

Exk + F Uk − b ≤ 0

(5b)

λk ≥ 0

(5c)

λ (Exk + F Uk − b) = 0

(5d)

T

and similar for time k + 1.
Note that equations (5d) are bilinear constraints but they
can for each row be modeled using a Big-M reformulation as
(i)
four linear constraints with a binary variable, zk , as
(i)

λk (eTi xk + fiT Uk − bi ) = 0 ⇒

where
Vk =

(i)

(i)

∗
minimize Vk∗ − Vk+1

s.t. xk+1 = Axk + Bu∗k
Uk∗

= arg min Vk

∗
Uk+1

(i)

0 ≤ bi − eTi xk − fiT Uk ≤ m2 (1 − zk ) (6)

and where Uk = [uTk , uTk+1 , . . . , uTk+N −1 ]T , xk is the current
measured state and the matrices, H, G, Q̄, E, F and b are
suitably defined. The objective function, Vk+1 , at time k + 1
is analogously defined.
Let us now formulate the stability test as the following
optimisation problem
∗
Uk∗ ,Uk+1
,xk

(i) (i)

0 ≤ λk ≤ m1 zk ,

1
1 T
U HUk + UkT Gxk + xTk Q̄xk
2 k
2

(4a)
(4b)

(i)

We can see that the binary variable zk forces either the
constraint, bi − eTi xk − fiT Uk , to be equal to zero or the
(i)
dual variable, λk , to be equal to zero, ensuring that the
complementarity constraint holds. In other words, we can view
(i)
this as that zk encodes whether the constraint is active or not,
(i)
i.e., zk = 1 ⇒ eTi xk + fiT Uk = bi .
Using the KKT conditions (5) and the binary reformulation (6) we can write the problem (4) using the notation,
T
T
y = [UkT , xTk , λTk , Uk+1
, xTk+1 , λTk+1 ]T and z̄ = [zkT , zk+1
]T as
1 T
y H̄y
2
s.t. Ēy = 0

(4c)

minimize

(7a)

y,z̄

= arg min Vk+1

(4d)

(7b)

¯
Āy ≤ b̄ + dz̄
Hence we want to find the state xk of the system that result
in the smallest possible difference in our candidate Lyapunov
function when controlled with the MPC controller.
If this difference is positive (or equal to zero) this means
that Vk is a valid Lyapunov function for the system and
hence it is stable. On the other hand if this difference is
negative this means that we have an increase in the Lyapunov
function candidate for som point xk and hence it is not a valid
Lyapunov function for the system.
Note that it is straightforward to extend the algorithm
to incorporate a condition of sufficient decrease and thus
asymptotic stability and robustness, e.g., by modifying the
objective function to

∗
Vk∗ − Vk+1
−  xTk Qxk − uTk Ruk ≥ 0
At this point it should be pointed out that it is assumed,
throughout the paper, that the MPC algorithms are recursively
feasible. This must of course be tested before one can apply the
proposed stability test. A method for performing the feasibility
test is presented in [13].
Note also that in the case of the MPC problem having
several non-unique solutions then the stability test derived here
is a pessimistic bound on the stability, i.e., it selects the worst
case combinations of optimal points.
Since the MPC problem (3) is a convex QP we can replace
it in (4c) and (4d) with the necessary and sufficient KKT

(7c)

where


H G
GT Q̄

 0
0
H̄ = 
 0
0

 0
0
0
0

B̄ A
Ē = H G
...

F
E
−F −E

. . .
0

. . .
0
Ā = 

...


...


...
...


b
M1 1 − b




0




0
,
b̄ = 


b


M1 1 − b




0
0

0
0
0
0
0
0
0 −H
0 −GT
0
0
0
FT
0

0
0
H

0
0
−I
I
0
0
0
0

...
0
0
0
F
−F
0
0


0
0
0
−G
−Q̄
0
−I
...
G


0
0

0

0

0
0

0

T
F


...
...
...
E
−E
0
0

0
−M1

 0

 M2
d¯ = 
 0

 0

 0
0







0

0

−I 
I

0
0 

0 

0 

0 

−M1 

0 
M2

3

Even though we have eliminated the bilinear constraints the
problem is still an indefinite QP, now with binary variables.
To make the final reformulation into a MILP one must
observe that


1
1
min. y T H̄y = min. min. y T H̄y
z̄
y
y,z̄ 2
2

is much faster to solve.
find xk
s.t.

and here we can replace the inner optimisation problem with
its KKT conditions in the minimisation over z̄ and thus we
arrive at
H̄y + ĀT η + Ē T µ = 0
¯ ≤0
Āy − b̄ − dz̄

(8a)
(8b)

Ēy = 0
T
¯ =0
η (Āy − b̄ − dz̄)

(8d)

η≥0

(8e)

(8c)

If we now multiply (8a) with 21 y T from the left we have

1 T
1
1
¯ Tη
y H̄y = − (y T ĀT η + y T Ē T µ) = − b̄ + dz̄
| {z }
2
2
2
=0

We see that in the optimum the objective can be equivalently
written as a linear term plus a bilinear term between a real
variable and a binary variable. We can see that the elements
of the bilinear term is either 0, if z̄ (i) = 0, or euqal to d¯Ti η
when z̄ (i) = 1. Hence we can introduce a new variable, w, as
z̄ T d¯T η = 1T w

(10b)

Ēy = 0

(10c)

− M3 (1 − z̄) ≤ w − d¯T η ≤ M3 (1 − z̄)

(10d)

− M3 z̄ ≤ w ≤ M3 z̄

(10e)

(10a)

¯ ≤0
− M4 (1 − q) ≤ (Āy − b̄ − dz̄)

(10f)

0 ≤ η ≤ M5 q

(10g)

Due to the special structure of the original problem there
exist a lot of structure in the problem (10) that should be
exploited to enhance the performance of the MILP representation.
III. E XPLOITING STRUCTURE IN THE MILP
First, let us observe that in MPC problems there are often
upper and lower bounds on the variables, e.g., umin ≤ uk+i ≤
umax . These constraints can not be fulfilled with equality
at the same time and hence the two corresponding binary
(i)
(j)
variables, zk and zk , in (6) can not be equal to one at the
same time. So we can introduce the constraint
(i)

(j)

zk + zk ≤ 1

(11)

for the appropriate indices i and j.
Let us now look at a single constraint eTi xk + fiT Uk ≤ bi
in the original MPC problem. This single constraint generates
through (6) the four constraints
eTi xk + fiT Uk − bi ≤ 0

where the elements wi can be modeled, yet again using the
Big-M formulation, with four linear constraints
−M3 (1 − z̄) ≤ w − d¯T η ≤ M3 (1 − z̄),

1
1
− b̄T η − 1T w < 0
2
2
H̄y + ĀT η + Ē T µ = 0

bi −

eTi xk

−M3 z̄ ≤ w ≤ M3 z̄

We can now combine all the pieces together and formulate
the stability problem (4) as the following MILP
1
1
minimize − b̄T η − 1T w
η,µ,y,w,z̄,q
2
2
s.t.

(9a)

H̄y + ĀT η + Ē T µ = 0

(9b)

Ēy = 0

(9c)

− M3 (1 − z̄) ≤ w − d¯T η ≤ M3 (1 − z̄)

(9d)

− M3 z̄ ≤ w ≤ M3 z̄

(9e)

¯ ≤0
− M4 (1 − q) ≤ (Āy − b̄ − dz̄)

(9f)

0 ≤ η ≤ M5 q

(9g)

Note that we are not really interested in finding the optimum
of problem (9), but rather find if there exist one point,
xk , where the objective is less than zero. Hence we can
reformulate the problem into a feasibility problem which often

−

fiT Uk
(i)
−λk
(i)
λk

≤

(12a)

(i)
m2 (1

−

(i)
zk )

≤0
≤

(12b)
(12c)

(i) (i)
m1 zk

(12d)

in (7c).
By formulating the KKT conditions for (7), for each of these
four constraints we have yet another binary variable, q (·) , in
(10f) and (10g). Since the binary variables forces a constraint
(i)
to be active we can see that if zk = 0 then (12a) and (12b)
can not both be active at the same time. This means that the
corresponding elements q (i) and q (nλ +i) can not both be equal
to one, i.e., we can constrain them as
(i)

q (i) + q (nλ +i) ≤ 1 + zk
(i)
Furthermore if zk
T
T
e i x k + f i U k = bi

(13)

= 1 we have from (12a) and (12b) that
which gives in (10f) and (10g) that

(i)

m4 (1 − q (i) ) ≤ eTi xk + fiT Uk − bi = 0
(nλ +i)

m4

(1 − q (nλ +i) ) ≤ bi − eTi xk − fiT Uk = 0

and hence we can in this case without any loss of generality
constrain both q (i) and q (nλ +i) to be equal to one. This is
done by adding the constraint
(i)

q (i) + q (nλ +i) ≥ 2zk

(14)

4

The corresponding argumentation can be used for (12c) and
(12d) to introduce the two additional constraints
(i)

q (3nλ +i) + q (4nλ +i) ≤ 1 + (1 − zk )
q (3nλ +i) + q (4nλ +i) ≥ 2(1 −

(i)
zk )

(15)
(16)

The constraints (13) -(16) concerns relations in one MPC
contraint, i, but as stated in (11) there exist relations also
between different MPC constraints. This gives additional relationships between the different binary variables which we
encode as the following constraints
(j)

(17)

(j)
1 − zk
(j)
2zk
(i)
1 − zk
(i)
1 − zk
(i)
2zk

(18)

q (nλ +i) ≤ 1 − zk
q

(2nλ +i)

≤

q (3nλ +i) + q (4nλ +i) ≥
q (nλ +j) ≤

q

(3nλ +j)

q

(2nλ +j)

≤

+q

(4nλ +j)

≥

(19)
(20)
(21)
(22)

The stability analysis optimisation problem to be solved
thus consist of the problem (10) with the additional constraints
(11) and (13) - (22). Note that the added binary constraints
are redundant in the original problem and does not effect
the optimal solution. They are only added to cut of binary
combinations in order to possibly increase the performance of
the solver.
Note that nowhere in the derivation of the algoritm do we
use the fact that it is the MPC controller objective function
that we have as candidate Lyapunov function, Vk . Hence we
can generalize the algorithm to the use of any other positive
definite Vk by appropriately modifying the matrix H̄.
IV. E XAMPLES
In this section we will look at three examples and try to
illustrate some properties and performance of the proposed
algorithm.
All examples have been implemented in Matlab using
YALMIP, [14]. The MILP problems have been solved using
the solver Gurobi 5.6.2, [15] and the LMI problems have been
solved using MOSEK 7.1, [16].
A. Sufficient but not necessary condition
The first example is taken from [17] where we consider the
following unstable system




1.216 −0.055
0.02763
xk+1 =
xk +
u
0.221 0.9947
0.002673 k
which we control using an unconstrained finite time MPC
controller with the objective function
Vk = minimize
uk+i

N
−1
X

xTk+i Qxk+i + Ru2k+i

i=0

and Q = 10I and R = 1. Since the controller is unconstrained
can we easily calculate for which prediction horizon length, N ,
the closed loop system is stable by looking at the eigenvalues
of the closed loop system matrix, A + BL. In the range

N = 1, . . . , 50 we calculate the eigenvalues of the closed loop
system and can conclude that it is stable for N ≥ 8. When we
use the algorithm derived in section II to test for stability we
get the result that the system is stable for N ≥ 21. Running
the LMI algorithm from [7] we obtain the same results.
Here we clearly see that the test is only a sufficient, but not
necessary, condition for stability since apparently there exist a
set of 8 ≤ N ≤ 20 where the system is stable but the objective
function does not constitute a valid Lyapunov function for the
closed loop system.
B. Complexity
To investigate the complexity of the proposed algorithm we
consider a very simple example. The system is the following
stable two state system




0.9744 0.0141
0.0106
xk+1 =
x +
u
−0.1023 0.9003 k
0.4878 k
This system is controlled with an input constrained MPC
controller with objective function
N
−1
X

xTk+i Qxk+i + Ru2k+i + xTk+N P xk+N

i=0

with Q = I, R = 10 and P =


10
1


1
and the constraints
2

−5 ≤ uk+i ≤ 5.
To investigate the computational complexity of the algorithm we have compared our developed MILP algorithm to
the LMI algorithm developed in [7].
Both algorithms verifies that the closed loop system is stable
for all tested prediction horizon lengths but the computational
time differs quite much between the two algorithms.
The following table shows the solvers computational time
in seconds for the two methods as a function of prediction
horizon.
N
LMI
MILP

2
0.03
0.06

4
0.52
0.17

6
5.36
0.28

8
22.42
0.42

10
80.06
0.61

As can be seen from the table the computation time grows
very rapidly with increasing prediction horizon for the LMI
algorithm while it remain relatively small for the MILP
approach.
C. Application to aircraft control
Let us now consider a more realistic example taken from
the aircraft industry.
We consider the stabilisation of the so called short period
dynamics, [18], of a fighter aircraft. The short period dynamics
can be approximated with a two state discrete time linear
system where xk = [αk qk ] is the angle of attack and
pitch rate, see Figure 1, and the input is the control surface
deflection.
The short period dynamics considered in this example are
those of the ADMIRE aircraft, [19], and we have linearised

5

Figure 1. The two states, angle of attack (α) and pitch rate (q) of the linear
short period dynamic’s approximation

the system at Mach 0.6 and altitude 4km. This result in the
following system




0.9798 0.0158
0.0106
xk+1 =
x +
u
0.1449 0.9787 k
0.4878 k
For this system we design a stabilising MPC controller of the
form (1) where the objective is a quadratic function
N
−1
X

xTk+i Qxk+i + Ru2k+i + xTk+N P xk+N

i=0



2 0
, R = 10 and P is the associated LQ cost.
0 0.1
The constraints are upper and lower bounds on the states and
control.


 
−10
10
≤ xk ≤
, −20 ≤ uk ≤ 20
−50
50

with Q =

Since the system is unstable and constrained it is ”generally
necessary” to have a final state constraint, [1]. We select the
final state constraint set to the invariant set of the associated
LQ controller, see e.g., [20] for more details.
Let us first consider the statement that it is necessary to
have a terminal state constraint. The MPC controller is by
construction stabilising and when we test the stability of the
closed loop system with the algorithm derived in previous
section we see that it indeed is stable for all N in the tested
range (N = 2, . . . , 10). The question is, is it still stable if we
would remove the terminal state constraint set?
Running the algorithm again, now without the terminal
constraint set in the controller, the test indicate that the
controller still is stable for all tested prediction horizons. To
verify this result we select 1000 random initial conditions and
simulate the closed loop system. As shown in Figure 2 all
initial conditions that are within the initially feasible set are
stable.
Let us now simplify the MPC controller even more by
introducing move blocking and analyse the stability of the
resulting closed loop system. We select the prediction horizon
N = 4 and introduce the move blocking structure


1 0
1 0

Uk = T Ûk , T = 
0 1
0 1

Figure 2. A simulation of the closed loop system from 1000 different initial
conditions for the MPC controller without terminal constraint. The blue circles
are the initial conditions, red lines are the trajectories and the black squares
are the final states. The cyan polytope is the set of initially feasible states.

where Ûk is the new reduced set of input signals. This blocking
structure gives a controller where the sampling time of the
controller output is half the internal sampling time used in the
predictions.
As pointed out in [2] the standard feasibility and stability
arguments can not be used to prove stability for this move
blocking strategy. Instead we aim to prove stability by using
our MILP test.
When applying the MILP algorithm it returns a sufficient
condition certificate that the move blocking MPC design is
stabilising. To verify this result we again simulate a set of
1000 random initial conditions throughout the state constraint
set, see Figure 3. It is clear from this figure that for all initially
feasible points the move blocking MPC controller remains
feasible and stabilises the system, as proven by the MILP test.

Figure 3. A simulation of the closed loop system from 1000 different initial
conditions for the move blocking MPC controller.

V. C ONCLUSIONS
In this paper we have derived an algorithm for a-postiori
stability analysis of MPC controllers.
The test allows the control engineer to verify more complex
MPC designs that normally is not possible to a-priori guarantee
stability for. We have shown that the stability test can be
written as a Mixed Integer Linear Program which has far less

6

computational complexity than other tests based on solving
LMIs.
As illustrated in one of the examples the test is only a
sufficient test for stability and if the test fails no conclusive
statements can be made about the lack of stability of the closed
loop system.
Future work consists of further analysing the structure of
the problem to reduce the computational complexity. Also
extending the algorithm to a robust stability test for systems
with disturbances and modelling errors.
Acknowledgement: This work has been done as a cooperation between Linköping University and Saab Aeronautics
and is funded by the Swedish Governmental Agency for
Innovation Systems (VINNOVA) and Centrum för industriell
informationsteknologi (CENIIT).
R EFERENCES
[1] D. Mayne, J. B. Rawlings, C. V. Rao, and P. Scokaert, “Constrained
model predictive control: Stability and optimality,” Automatica, vol. 36,
no. 6, pp. 789–814, jun 2000.
[2] R. Cagienard, P. Grieder, E. Kerrigan, and M. Morari, “Move blocking
strategies in receding horizon control,” Journal of Process Control,
vol. 17, no. 6, pp. 563–570, 2007.
[3] M. N. Zeilinger, M. Morari, and C. N. Jones, “Soft Constrained Model
Predictive Control With Robust Stability Guarantees,” IEEE Transactions on Automatic Control, vol. 59, no. 5, pp. 1190–1202, 2014.
[4] L. Chisci and G. Zappa, “Dual mode predictive tracking of piecewise
constant references for constrained linear systems,” International Journal of Control, vol. 76, no. 1, pp. 61–72, jan 2003.
[5] D. Limon, I. Alvarado, T. Alamo, and E. Camacho, “MPC for tracking
piecewise constant references for constrained linear systems,” Automatica, vol. 44, no. 9, pp. 2382–2387, sep 2008.
[6] J. A. Primbs and V. Nevistic, “A new approach to stability analysis for
constrained finite receding horizon control without end constraints,” Ieee
Transactions on Automatic Control, vol. 45, no. 8, pp. 1507–1512, 2000.
[7] J. A. Primbs, “The analysis of optimization based controllers,” Automatica, vol. 37, no. June, pp. 3297–3301, 2001.
[8] N. S. Ahmad, J. Carrasco, and W. Heath, “A Less Conservative LMI
Condition for Stability of Discrete-Time Systems with Slope-Restricted
Nonlinearities,” IEEE Transactions on Automatic Control, vol. PP,
no. 99, pp. 1–1, 2014.
[9] G. Li, W. P. Heath, and B. Lennox, “Min-max model predictive
control as a quadratic program,” IFAC Proceedings Volumes (IFACPapersOnline), vol. 16, no. 7, pp. 263–268, 2005.
[10] ——, “The stability analysis of systems with nonlinear feedback
expressed by a quadratic program,” Proceedings of the 45th IEEE
Conference on Decision and Control, vol. 2, no. 3, pp. 4247–4252, 2006.
[11] C. Løvaas, M. M. Seron, and G. C. Goodwin, “A dissipativity approach
to robustness in constrained model predictive control,” Proceedings of
the IEEE Conference on Decision and Control, vol. 2, pp. 1180–1185,
2007.
[12] M. Korda and C. N. Jones, “Stability and Performance Verification
of Optimization-based Controllers,” 2015. [Online]. Available: http:
//arxiv.org/abs/1501.03919
[13] J. Löfberg, “Oops! I cannot do it again: Testing for recursive feasibility
in MPC,” Automatica, vol. 48, no. 3, pp. 550–555, 2012.
[14] ——, “YALMIP: A Toolbox for Modeling and Optimization in MATLAB,” in Proceedings of the CACSD Conference, Taipei, Taiwan, 2004.
[15] G. O. Inc., “Gurobi Optimizer reference manual,” p. 572, 2014.
[Online]. Available: www.gurobi.com
[16] M. ApS, “The MOSEK optimization toolbox for MATLAB manual.
Version 7.1 (Revision 28).” 2015.
[17] J. Löfberg, “Minimax approaches to robust model predictive control,”
Ph.D. dissertation, Linköping University, 2003.
[18] B. L. Stevens and F. L. Lewis, Aircraft Control and Simulation, 2nd ed.
John Wiley & Sons, 2003.
[19] L. Forssell and U. Nilsson, “ADMIRE The Aero-Data Model In a Research Environment Version 4.0, Model Description,” Swedish Defence
Research Agency, Tech. Rep. December, 2005.

[20] J. B. Rawlings and D. Q. Mayne, Model Predictive Control: Theory and
Design, 1st ed. Nob Hill Publishing, 2009.

