**Model Predictive Control**

**of Embedded Systems**

**Dissertation**

### zur Erlangung des akademischen Grades

**Doktoringenieur (Dr.-Ing.)**

### von

**Juan Pablo Menendez Zometa**

### geboren am 29. Mai 1980 in Santa Tecla

### genehmigt durch die Fakultät für Elektrotechnik und Informationstechnik

### der Otto-von-Guericke-Universität Magdeburg

### Gutachter:

### Prof. Dr.-Ing. Rolf Findeisen

### Prof. Dr. Daniel Limon

### eingereicht am 23. August 2016

### Promotionskolloquium am 12. Mai 2017

This thesis is the result of the work done as a research assistant at the Systems Theory and Automatic Control group of Prof. Rolf Findeisen, at the Otto-von-Guericke Univer-sity in Magdeburg, Germany.

This work would have not been possible without the support of innumerable people. Special thanks go to Prof. Findeisen for his supervision and support, and to friends, colleagues and family for their encouragement.

**1. Introduction** **1**

1.1. Linear model predictive control for embedded systems . . . 2

1.2. Contribution . . . 6

1.3. Outline . . . 8

**2. Model predictive control for embedded systems** **11**
2.1. Embedded digital control systems . . . 11

2.1.1. Digital control of linear dynamic systems . . . 11

2.1.2. Cyber-physical systems: real-time embedded systems . . . 15

2.2. On-line optimization . . . 19

2.2.1. Basics of convex optimization . . . 19

2.2.2. Quadratic programs . . . 20

2.2.3. Optimization theory for quadratic programs . . . 21

2.3. Model predictive control for linear systems . . . 23

2.3.1. A basic MPC setup . . . 25

2.3.2. MPC as a QP . . . 26

2.3.3. General moving horizon control formulation . . . 27

2.4. Summary . . . 28

**3. Tailored on-line optimization software tools for MPC** **29**
3.1. Exploiting the properties of the MPC algorithm . . . 29

3.1.1. Known-ahead maximum computation time . . . 29

3.1.2. Partial use of the solution by the controller . . . 30

3.1.3. Similarities between consecutive problems . . . 31

3.1.4. Special structure of the data . . . 33

3.2. Tailored optimization software tools . . . 35

3.2.1. Interior point methods . . . 36

3.2.2. Active set methods . . . 38

3.2.3. Gradient methods . . . 40

3.2.4. Comments on explicit methods . . . 41

3.3. A novel optimization algorithm for embedded MPC . . . 41

3.3.1. Fast gradient method . . . 42

3.3.2. Augmented Lagrangian method . . . 45

3.3.3. The ALM+FGM algorithm for MPC . . . 49

3.3.4. ALM+FGM for embedded MPC . . . 53

3.4. Summary . . . 55

**4. General description of multistage problems** **57**
4.1. General formulation of multistage problems . . . 57

4.1.1. Motivational example . . . 58

4.1.2. Abstract formulation of multistage problems . . . 59

4.2. Reformulation as a condensed optimization problem . . . 64

4.2.1. Reformulation as a general QCQP . . . 64

4.2.2. Reformulation as a QCQP in standard form . . . 67

4.2.3. Special case: condensed QP . . . 70

4.3. High-level multistage specification language . . . 72

4.3.1. Code generation of condensed formulation . . . 74

4.4. Summary . . . 75

**5. muAO-MPC: a free code generation tool for embedded MPC** **77**
5.1. Core features . . . 77

5.2. Automatic generation of C code . . . 79

5.2.1. Forming and solving the condensed QP . . . 80

5.2.2. Solving the QP with the ALM+FGM algorithm . . . 81

5.2.3. Further controller performance improvements . . . 82

5.3. Examples: code generation for a microcontroller . . . 83

5.3.1. Setup description . . . 83

5.3.2. Considered Embedded Hardware . . . 84

5.3.4. Discussion . . . 85

5.4. Summary . . . 88

**6. Application examples** **91**
6.1. Low-end example: A direct current motor . . . 91

6.1.1. System description . . . 91

6.1.2. Generating a fast embedded MPC controller . . . 93

6.1.3. Results . . . 93

6.2. High-performance example: An autonomous vehicle . . . 95

6.2.1. System description . . . 95
6.2.2. MPC Implementation . . . 97
6.2.3. Discussion . . . 99
6.3. Summary . . . 101
**7. Conclusions** **103**
7.1. Outlook . . . 104

**A. Forming a condensed parametric quadratic program** **119**
**B. System Matrices** **123**
B.1. Simulation examples . . . 123

The internet of things is a novel paradigm based on networked physical devices that interact with their environment. This paradigm enables the improvement of existing technologies, and is likely to serve as the basis for yet unforeseen technological ad-vancements. An internet-of-things device works on a computer network, and typically integrates sensors and actuators. At the core of any device is an embedded computer, commonly with very low computational capabilities. Model predictive control is an ad-vanced control method likely to play a role in future internet-of-things applications as a way to improve the performance of networked cyber-physical systems. Therefore, the easy implementation of model predictive control on embedded hardware is of interest.

Model predictive control is based on the repeated solution of an optimization problem that allows to naturally handle multi-input multi-output systems subject to constraints. At every sampling time, inputs to the system are determined that optimize the pre-dicted state evolution up to a certain future time instance. For this, an optimization problem must be solved in real time at each sampling period. This makes the implemen-tation of model predictive control on embedded systems challenging. The difficulty on the implementation arises due to the combination of two factors: the limited computa-tional capabilities of today’s embedded hardware and the relatively high computacomputa-tional demands of solving the necessary optimization problem. The efficient and reliable im-plementation of model predictive control requires a tailored imim-plementation.

A starting point for the efficient implementation is that for a given model predictive control formulation, the optimization problem belongs to a often narrow problem class (e.g. a quadratic program) with a particular structure. This structure does not change for many applications, and from one sampling time to the next only the problem data is different. The few available tools often exploit the special multistage structure of model predictive control problems. This leads to code whose computational complexity grows linearly with the horizon length by exploiting the sparsity of the problem. An alternative are, especially for short horizon lengths, condensed formulations.

We present a universal code generation tool that exploits condensed parametric formu-lations tailored for linear time-invariant discrete-time model predictive control. As base we introduce a simple modeling language that presents an easy to formulate model pre-dictive control problems. Based on this language, we outline how one can automatically generate C code that is tailored towards dense solution approaches, which are espe-cially suitable for low-performance and memory-constrained embedded processors. We furthermore implement an optimization algorithm tailored for model predictive control that explicitly takes into account the limitations of embedded hardware.

The presented modeling language is intuitive and easy to use, and offers great flexibil-ity in the problem formulation. We demonstrate via several examples that the outlined code generation approach allows for fast solutions even on low-cost embedded platforms.

Das Internet der Dinge ist ein neues Paradigma, welches auf vernetzten physikalischen Geräten basiert, die mit ihrer Umgebung interagieren. Dieses Paradigma ermöglicht die Verbesserung von existieren Technologien und wird wahrscheinlich als Grundlage für noch unvorhergesehene technologische Fortschritte dienen. Ein Internet-der-Dinge-Gerät arbeitet auf einem Computer und bezieht typischerweise Sensoren und Aktoren mit ein. Im Kern eines jeden Gerätes befindet sich ein eingebetteter Computer, welcher für gewöhnlich über eine geringe Rechenleistung verfügt. Modellprädiktive Regelung ist eine erweiterte Methode in der Regelungstechnik, die wahrscheinlich eine zukünftige Rolle in der Anwendung des Internets der Dinge spielen wird. Sie stellt einen Weg dar, die Leistung von vernetzten cyber-physikalischen Systemen zu verbessern. Deshalb ist die einfach umsetzbare Implementierung von modellprädiktiver Regelung auf eingebetteter Hardware von Interesse.

Modellprädiktive Regelung ist ein fortschrittliches Regelungsverfahren, basierend auf der wiederholten Lösung eines Optimierungsproblems, mit dem auf einfache Weise auch Mehrgrößensysteme und Systeme mit Beschränkungen behandelt werden können. Zu je-dem Abtastzeitpunkt wird ein Eingang für das System berechnet, welcher die prädizierte Systementwicklung optimiert. Somit muss zu jedem Abtastzeitpunkt wiederholt und in Echtzeit ein Optimierungsproblem gelöst werden. Die wiederholte Lösung dieser Probleme macht die Implementierung der Modellprädiktiven Regelung auf eingebet-teten Systemen zu einer Herausforderung. Die Schwierigkeiten, die bei der Umsetzung entstehen sind auf eine Kombination von zwei Faktoren zurückzuführen: die begren-zte Rechenleistung eingebetteter Systeme sowie die relativ hohen Rechenanforderungen zur Lösung der Optimierungsprobleme. Daraus folgt, dass eine effiziente Umsetzung der Modellprädiktiven Regelung eine maßgeschneiderte problemspezifische Implemen-tierung erfordert. Für ein gegebenes Modellprädiktives Regelungsproblem gehört das entsprechende Optimierungsproblem oftmals zu einer speziellen Klasse (z.B. Quadratis-ches Programm) mit spezieller Struktur. In vielen Anwendungen bleibt diese Struktur

von einem Abtastzeitpunkt zum nächsten konstant, es ändern sich lediglich die Prob-lemdaten. Die wenigen, verfügbaren Software-Tools nutzen die inhärente Struktur der Modellprädiktiven Regelung aus, was zu einem linearen Wachstum der Rechenkomplex-ität mit der Länge des Prediktionshorizonts führt. Diese Implementierungen nutzen oft dünnbesetzte Matrizen zur Beschleunigung der Lösung des Optimierungsproblems aus. In einigen Fällen, insbesondere für Probleme mit kurzen Prädiktionshorizonten auf eingebetteten Systemen, sind allerdings kondensierte Formulierungen besser geeignet. In dieser Arbeit wird ein Software-Tool zur Codegenerierung vorgestellt, welches auf einer parametrischen Formulierung der linearen zeitinvarianten diskreten Modellprädiktive Regelung beruht. Hierzu wird zunächst eine einfache Modellierungssprache vorgestellt, welche erlaubt modellprädiktive Regelungsprobleme einfach zu formulieren. Auf Basis dieser Sprache wird umrissen, wie automatisch C-Code generiert werden kann, der die Verwendung kondensierter Lösungsansätze ermöglicht. Darüber hinaus wird ein Op-timierungsalgorithmus implementiert, der die Einschränkungen eingebetteter Systeme berücksichtigt. Die vorgestellte Modellierungssprache ist intuitiv und bietet große Flex-ibilität bei der Problemformulierung. Anhand einiger Beispiele wird demonstriert, dass der vorgestellte Ansatz der Codegenerierung zur schnellen Lösung von Modellprädiktiven Regelungsproblemen, auch auf preisgünstigen eingebetteten Systemen, ermöglicht.

Despite many conceptual and theoretical advances, the systematic use of advanced con-trol approaches for constrained system in industrial applications is still limited [1]. Often simpler unconstrained methods like a proportional-integral control might be sufficient to stabilize a system with constraints. In such cases, a more complex algorithm that can ex-plicitly take into account the constraints may be, due to technical and economic reasons, hard to justify even if it delivers superior closed-loop performance [1]. A notable excep-tion is the process industry, where model predictive control (MPC) [2, 3, 4, 5, 6, 7, 8] has been consistently applied for nearly four decades [9]. This is due several factors, including that the processes are slow and that operation near state constraints offers great economic benefits.

In recent years, the use of MPC has been applied to a broad range of applications, e.g. path-following in robotics [10], building climate control [11], control of water canals [12], and biomedical applications such as artificial pancreases [13].

In the case of fast embedded systems subject to constraints, MPC has been tradition-ally considered difficult to apply due to the high computational demands of the MPC algorithm relative to the computational capabilities of embedded computers. Therefore, a widespread use of MPC on embedded applications has been limited in the past and has only gained momentum in recent years. The constantly decreasing cost/performance ratio of computers, and, most importantly, the recent improvement in MPC-tailored op-timization algorithms has made the implementation of MPC in embedded hardware possible [14, 15, 16, 17], albeit still challenging on low-cost devices. This work attempts to further simplify the application of MPC on low-cost embedded systems by introducing a new code generation software tool for MPC that explicitly considers the limitations of low-cost embedded hardware. The two main components of the generated code are an efficient optimization solver and data structures tailored for embedded applications.

**1.1. Linear model predictive control for embedded**

**systems**

Model predictive control is an advanced control method that takes into account system constraints, and is based on repeated optimization using information of a prediction model. According to [7, Ch. 1], model predictive control is the only generic control method that can routinely deal with constraints. Furthermore, [8] states that MPC is “perhaps, the most general way of posing the process control problem in the time domain”. The fact that it is generic implies that a broad range of applications that deal with the control of constrained systems can be considered within the MPC framework. A further appeal of MPC is that it is conceptually intuitive, although its inner operation rely on complex and less intuitive concepts. Many types of MPC formulations exist. What all of them have in common are the explicit use of a mathematical model to make predictions and a receding horizon optimization [18]. MPC has been developed at least in three major streams: the industrial or classical setting, the one based on adaptive control, and the synthesis approach [8]. We focus here on the synthesis approach. Two of the main characteristic of this approach are the use of state-space models and stability guarantees [18].

To this end, we consider the linear discrete-time dynamical system

*x*+*= Ax + Bu + h(w)*

subject to input and state constraints

*u ∈ U, x ∈ X,*

*with x the current state, x*+ _{the successor state, u the input to the system and h(w) a}*function affine in the parameter w. This parameter represents a reference trajectory or*
a known disturbances (or a combination of them), for example. The input constraints

*U* are often determined by actuator limits. Examples are the minimum and maximum

voltage that can be applied to an electric motor, or the maximum aperture of a valve.
These are usually hard constraints, i.e. they are enforced by the hardware and, under
*normal operation, cannot be violated. The state constraints X are usually determined*
by several factors, such as safety, comfort, and efficiency. These are commonly artificially

introduced by design. For example the acceleration of a vehicle can be artificially limited for comfort reasons, or the temperature or pressure of a chemical process can be limited for economic and safety reasons. These types of constraints are often considered as soft. In practice they may be violated, e.g. due to several sources of uncertainty in the control loop like disturbances and model-plant mismatch. In other words, soft constraints are desired limits that, in case of not being respected, the controller can relax.

To control this type of constrained system, we use the following linear MPC problem
formulation:
minimize_{u}*N −1*X
*j=0*
*`(xj, uj, wj) + `N(xN, wN*)
*subject to xj+1* *= Axj+ Buj* *+ h(wj*) *j* *= 0, ..., N − 1*
*uj* *∈ U* *j* *= 0, ..., N − 1*
*xj* *∈ X* *j* *= 0, ..., N − 1*
*xN* *∈ Xf* *⊆ X*
*x*0 *= x.*
(1.1)

*The stage cost `(xj, uj, wj), the terminal cost `N(xN, wN*) and the terminal constraint set

*Xf* *constitute the three ingredients that help establish stability of the closed-loop system,*

*[19]. The horizon length is given by N. The stage cost is commonly used to determine*
the desired behavior by way of weighting matrices. Both, the stage and terminal costs,
are commonly quadratic functions. The MPC formulation is an optimization problem
that needs to be solved at each sampling time. The solution to the problem for the state

*x* *at the current sampling time is used to compute the current control input u (see [3,*

Ch. 2] for further details).

*Note that the stage cost `(xj, uj, wj) is repeated for j = 0, · · · , N − 1. Thus, in this*

*work we would refer to MPC as belonging to the broader class of multistage problems.*
This classification also includes moving horizon estimation problems [3, Ch. 3].

We are particularly interested in MPC problems like (1.1) with a quadratic cost and affine inequality constraints and affine equality constraints (i.e. the linear dynamics). Such problems lead to quadratic programs (QP). To find a solution using iterative algo-rithms, the MPC problem is often brought into a special form. We distinguish between two types of problem formulations: sparse and condensed.

On the one hand, sparse formulations exploit the multi-stage nature of MPC to achieve a computational complexity that grows linearly with the horizon length, instead of cu-bically in the case of non-tailored interior point methods [20]. On the other hand, condensed problem formulations use dense matrices without a directly exploitable struc-ture. The computational complexity of a condensed solver typically grows cubically or quadratically with the horizon length. In general, condensed formulations are preferred in cases where the horizon length is short. In embedded applications, the horizon length is often selected to be as short as possible, while still delivering suitable performance. Thus, condensed formulations can be of computational advantage for embedded applica-tions. One specific problem that benefits from a condensed formulation is the case when the number of states is much larger than the number of inputs. For example, stochastic MPC formulations based on polynomial chaos expansion show this property [21].

One often cited drawback of MPC is that its computational requirements are sig-nificantly higher in comparison to simpler unconstrained controllers. This is partic-ularly relevant when implementing MPC on an embedded computer. An embedded computer in general provides only a small fraction of the computational performance of a modern desktop computer. Implementing MPC on an embedded computer is in general challenging. Although nowadays very powerful off-the-shelf microprocessors are cheaply available, in some applications, for technical or economical requirements, low-cost resource-limited embedded processors need to be used. Technical limitations may include low energy consumption, small footprint, etc. Furthermore, in cost-sensitive applications, the cheapest processor that can do the job is typically selected.

A great amount of research has been performed recently that aims at speeding up MPC implementations. Of particular interest are automatic code generation methods for MPC for embedded systems. [22, 23, 24, 25]. The current state of research on this field can be split into two branches: tailored solvers and automatic code generators.

*One way to speed up finding a solution is the use of explicit MPC [26]. In this*
approach, the solutions are computed off-line and stored in a lookup table. This approach
has high memory requirements and is in general limited to problems of small size [27].
In the following, we restrict our attention to MPC approaches that require solving an
optimization problem on-line. We focus on optimization algorithms that iteratively
search for a solution to the optimization problem. That is, the algorithm looks for a
point that minimizes the value of the cost function and that, at the same time, satisfies
all constraints.

A large variety of literature exist in the field of optimization [28, 29, 30]. Of particular
*interest in this work are convex optimization problems [31]. Although convex problems*
can be solved reliably, an off-the-shelf general purpose solver might not be fast enough
for an MPC application. Furthermore, such a solver might be difficult to use on an
embedded platform due to software dependencies and hardware limitations. This has
motivated the development of solvers that take into account the MPC characteristics
and limitations of embedded applications. Among the most commonly exploited
charac-teristics are the MPC structure and warm starting. The former commonly refers to the
exploitation of the banded structure of the problem matrices arising due to the
predic-tive nature of MPC. Warm start refers to the use of the solution to the MPC problem
at the previous sampling time as the initial guess for the solution of the problem at the
current sampling time.

One of the first reported on-line algorithms that fully exploited the MPC structure
for linear systems was presented in [20], which used a primal-dual interior point method
(IPM). This has been the base for recent code generation tools tailored for MPC like
the ones presented in [22, 23]. Primal-dual IPMs spends the first iterations finding
the so-called central path. Once in that path, the algorithm can quickly converge to a
solution. Hence, these methods do not typically make much progress on improving the
quality of the solution on the very first iterations. Furthermore, primal-dual IPM are
difficult to warm start [27]. This might restrict the application of tailored primal-dual
IPMs on embedded applications where acceptable performance can be attained by rough
solutions and only very limited time is available to compute the solutions. Alternatively,
a tailored primal barrier IPM that can be warm started was discussed in [32]. This
method is in general not as fast and accurate as primal-dual methods, but may offer a
good trade-off between speed an accuracy for embedded applications. Finally, active-set
methods have also been tailored to efficiently solve optimization problems arising from
*MPC, with perhaps the most popular being the online active set strategy [27], and its*
implementation, called qpOASES [33].

In general, the theoretical computational bounds are much larger than the ones ob-served in practice for interior point methods [31, 34] and active set methods [29, 35]. In safety-critical applications, error bound certification might be required. First-order methods in general have tighter bounds, and may be preferred for such applications. Some theoretical contribution in this direction have been made in [36, 37, 38].

generation tools. These tools commonly consist on a parser of a high-level domain-specific language, and a code generator of a tailored optimization solver. Typically the C programming language [39] is used for the generated code. Commonly, the generated code does not depend on external software libraries, i.e. an appropriate C compiler (e.g. [40]) suffices to generate executable binaries for a hardware target platform.

Several parsers for general purpose optimization problems exist, e.g. YALMIP [41],
AMPL [42], and CVX [43]. However, these are mainly used to feed data to general
purpose optimization routines, i.e. the solver does not take into account many of the
characteristics of each problem at hand, which makes them in general less efficient. In
*contrast, CVXGEN, [44] generates embedded solvers for convex quadratic programs. In*
the CVXGEN context, embedded means that, given a problem with certain
character-istics, solvers that exploit these characteristics are generated. CVXGEN reports several
orders of magnitude improvements in the solution time compared to CVX [44].
Neverthe-less, the C code generated by CVXGEN grows very quickly in size (see Subsection 3.2.1
for further details). While CVXGEN can handle MPC problems [45], the generated code
does not exploit the structure inherent in MPC. The combination of QCML and ECOS
as presented in [24], the former being a parser and the latter being a stand alone general
purpose second-order cone program solver [46], is based on similar ideas as CVXGEN.
Both CVXGEN and QCML do not eliminate the equality constraints due to the system
dynamics, leading to a sparse formulation.

In [47] two code generation tools are presented. One is FiOrdOs, which based on first-order methods that do not exploit structure inherent in MPC [48]. The other one is FORCES, a structure-exploiting primal-dual interior point method [22]. Parsers specific to nonlinear MPC have also been developed. The ACADO toolkit [49] is a C++ software tool that provides C-code generation capabilities for continuous-time nonlinear MPC problems [50]. Although it is possible to use the ACADO toolkit for discrete-time linear systems as well, the generated code is bulkier than code tailored for linear discrete-time MPC problems.

**1.2. Contribution**

This work addresses the computational drawback of MPC without compromising its
*advantages. The main contribution is a free MPC software tool called µAO-MPC that*

simplifies the deployment of MPC controllers, while at the same time taking into
*consid-eration the computational limitation inherent to embedded processors. Because *

*µAO-MPC relies on a flexible problem formulation, it can be used to deploy moving horizon*

estimators as well. In this sense, we retain the main advantages of the MPC framework, namely its intuitiveness and generality, while at the same time addressing one of its main cited drawbacks, i.e. prohibitive computational requirements for embedded computers.

*Broadly speaking, µAO-MPC consist of three parts: a parser for a simple MPC *
prob-lem specification, a C code generator and a tailored optimization algorithm. The typical
workflow is as follows: the user specifies an MPC problem using a language that closely
resembles the way MPC problems are expressed mathematically. This specification is
parsed and transformed into an optimization problem, in particular a quadratic program
(QP) with dense matrices. Finally, C code is generated that contains the QP together
with a first-order optimization algorithm that can approximately solve the optimization
problem efficiently.

In this work, we follow an approach similar to CVXGEN in the way of how MPC
problems are expressed (i.e. simple text-based mathematical representation). We focus
however on code generation of embedded MPC solvers for computationally limited
sys-tems using first-order methods, similar to FiOrdOs. At the moment, we have focused on
*solving QPs and our results show that µAO-MPC generates code of much smaller size*
than CVXGEN. One key difference from our approach and all other existing parser/code
generators for linear systems is that we eliminate the equality constraints to derived a
condensed formulation. Furthermore, the solver we use finds control inputs
appropri-ate for certain types of embedded control faster than CVXGEN. This control inputs
are based on rough approximate solutions, which are nevertheless in many applications
good enough for control. We do not compare our results against FiOrdOs, but we would
*expect that µAO-MPC and FiOrdOs deliver similar performance in many situation. A*
full comparison with many optimization algorithms is beyond the scope of this work.
The main contribution is therefore not in the field of convex optimization, but rather
on the implementation aspects of automatic code generation tools for model predictive
control on embedded systems.

In this work we extend the results presented in [25] by introducing a method and a parser that automatically generates easily usable code for the efficient on-line solution based on a condensed formulation. This approach, which is specially tailored for linear time-invariant MPC problems, allows (besides QP formulations) the consideration of

linear programs (arising e.g. in min-max robust MPC [51]), quadratically constrained
quadratic programs (arising in robust MPC by adding ellipsoidal constraints [52]) and
second-order cone programs (arising in stochastic MPC [53]). We show how to
system-atically go from the MPC description to the condensed formulation using the proposed
*method. Thus, we extend µAO-MPC capabilities two-fold: first, the presented approach*
allows to generate data tailored for linear programs, quadratically constrained programs
and second order cone programs in dense form. Second, we extend the family of
prob-lems that can be considered by introducing a multistage modeling language. This allows
not only to consider a broader range of multistage problems, but also moving horizon
estimation problems and special filtering problems. Note that we focus here on the
gen-eration of the dense problem data, and not on a particular solver. To the best of our
knowledge, this is the first tailored tool for linear time-invariant discrete-time MPC that
can automatically parse, condense the data, and generate easily usable code. Loosely
related, but tailored towards continuous-time nonlinear systems is the ACADO Toolkit
[49, 50].

**1.3. Outline**

This work is organized as follows. In Chapter 2 we describe the overall setup, namely digital control of embedded systems and the MPC algorithm. We discuss the properties of embedded systems most relevant to MPC applications. In particular, we discuss the computational limitations of embedded systems. Furthermore, we present the theoretical foundations of convex optimization applied in particular to convex quadratic programs. In Chapter 3 we discuss optimization algorithms tailored for embedded MPC. We briefly discuss some state-of-the-art approaches. The chapter focuses on a novel opti-mization algorithm tailored for embedded MPC applications. The algorithm is based on an augmented Lagrangian method combined with Nesterov’s fast gradient method (the ALM+FGM algorithm).

In Chapter 4 we introduce a general MPC language specification. We show how some types of MPC problems specified using this language can be transformed into condensed parametric quadratically constrained quadratic programs in a structured way.

*Chapter 5 presents µAO-MPC, an automatic code generation tool for MPC written*
in Python. The software generates portable C-code for MPC problems specified in

the language introduced in Chapter 4. By default the generated condensed QPs are
solved using the tailored ALM+FGM algorithm discussed in Chapter 3. We show with
*some simulation examples µAO-MPC’s suitability for embedded applications, taken into*
account the properties discussed in Chapter 2.

*In Chapter 6 we present two application examples that again demonstrate that *

*µAO-MPC is appropriate for a specific type of embedded applications, namely those that can*

be warm started, are well conditioned and where rather rough approximate solutions deliver good closed-loop performance.

**embedded systems**

Nowadays, most closed-loop control schemes are implemented on digital systems. This is
mainly due to the flexibility provided by software implementations, and the ubiquity of
microprocessors. Furthermore, there is a strong body of control theory and software tools
that covers the design of digital control systems. Particularly well developed is the theory
of linear unconstrained single-input single-output systems [54, 55]. Usually, the analysis
*is performed in the frequency domain (i.e. z-transform for discrete-time systems). In*
the case of unconstrained multi-input multi-output (MIMO) systems, modern control
methods exist that are well established, such as H2and H∞control [54, Ch. 40]. Optimal
control methods together with state-space representations are commonly used for the
*control synthesis of MIMO systems [56, 57]. In the case of constrained MIMO systems,*
model predictive control (MPC) is the most widespread technique [7, 8].

This chapter lays the theoretical foundations needed for a basic understanding of the main characteristics of MPC for embedded systems.

**2.1. Embedded digital control systems**

This section describes the most relevant characteristics of digitally implemented control systems.

**2.1.1. Digital control of linear dynamic systems**

We start this section with a few definitions motivated by [55]. We expand later on some concepts that are key for the further development of this topic.

Figure 2.1 shows the general scheme of a digital control system. We first focus on the right part of the figure, which represents the analog (or continuous) part of the control

CPU DAC
ADC
Plant
Sensor
*ˆuk* *u(t)* *y(t)*
*w(t)*
*v(t)*
*ˆyk*
Clock

Figure 2.1.: Block diagram of a basic digital control system, adapted from [55].
*system. These are represented with grey boxes. We consider the plant to be any physical*
process that needs to be controlled. Most of the physical plants of interest are dynamical
*systems which are continuous in nature. We want to modify the behavior y(t) of the*
*plant with a control input u(t). The variable w(t) represents disturbances. Finally, we*
*capture the magnitude of physical variables of the plant via sensors. We use v(t) to*
denote disturbances or noise in the sensor. All these variables change continuously in
time.

Let us now turn our attention to the middle part of Figure 2.1, which represents the
bridge between the digital and the analog part of the control system, represented by the
grey/white boxes. The analog-to-digital converter (ADC) samples a physical variable
and then converts this value into a digital number representation. We assume that these
samples are taken at constant time intervals, which are driven by the clock. We call
*the interval between any two samples the sampling period. The signals generated by the*
*ADC vary only at discrete times, more precisely at every t = kT , with T the sampling*
*period, and k = 0, 1, . . .. This type of signals are called discrete signals. For simplicity*
*of notation, we write any discrete signal y(kT ) in the form yk*. A system that has an

*interface to continuous signals (e.g. y(t)) and discrete signals (e.g. yk*) is known as a

as a digital number with a limited number of binary digits. This introduces the effect
*known as quantization, which refers to the effect of representing yk* *as a value ˆyk* that

has been approximated to a certain precision. Digital computers work, therefore, with

*digital signals, i.e. signals that are both quantized and discrete.*

In the middle part of Figure 2.1 the digital-to-analog converter (DAC) is also depicted,
which performs the opposite function of the ADC. The same concepts apply. The main
*difference is that the DAC usually implements a zero-order hold to keep the analog signal*
constant during a whole sampling period (the ADC also does this, but only for a fraction
of the sampling period).

Finally, on the left side of the figure we have the digital part of our control system,
represented by white boxes. The central processing unit (CPU) is in charge of computing
*a control action ˆuk* *based on the current sensor information ˆyk* via a control algorithm.

Two of the most important factors that affect the performance of sampled data systems are quantization and the sampling period. If the sampling period is sufficiently small, and the quantization effects are negligible, a sampled-data control system might behave in some cases almost identically to a continuous-time control system [54, Ch. 15]. Usually, in cost-sensitive embedded control systems the sample period and quantization effects are not negligible and need to be taken into account in the design of the control system. We further expand the concepts and effects of quantization and the sampling period in the following.

**Quantization**

Loosely speaking, quantization is an effect where a digital number with limited precision
*approximates a real number [54, Ch. 15]. Consider a signal z ∈ Qz* ⊆ R in the

*interval [z, z] that is represented by the quantized value ˆz. The set Qz* can represent

the values of a continuous signal (e.g. Q*z* *= R) or of a digital signal (i.e. z is itself*

*quantized). The bounds z and z represent, for example, actuator or sensor limits. In*
this context we define quantization Φ : Q*i* → Q*o* as a nonlinear map from the set

Q*i* *= {z ∈ Qz* *⊆ R | z ≤ z ≤ z} to the set Qo* *= {ˆz ∈ R | ˆz = z + sρ(z − z)}, where*

*ρ* = 2*−l, with l ∈ Z is the number of bits available to represent the digital number and*
*s ∈ Z, 0 ≤ s ≤ (2l*_{−}1) is the full range of raw binary values. Note that Q

*o* is defined

by an affine function of integer numbers. Other definitions are also possible, e.g. using logarithmic functions.

The quantization Φ is usually defined in either of the following two ways. The first is
truncation, which means mapping a real number in Q*i* to the largest previous number

in Q*o*. The second way is round-off, which is more commonly used in control computers.

Round-off simply maps any number in Q*i* to the nearest number in Q*o*. We assume in

the following that the quantization function uses round-off.
*We define the maximum quantization error β for Φ as*

*β*= 1

2*ρ(z − z).* (2.1)

*Furthermore, for the mapping Φ : z 7→ ˆz we have β ≤ |z − ˆz|. For example, consider*
the quantization happening in the ADC and the DAC. In the former case, we have
*Φ : y(kT ) 7→ ˆyk* with Q*z* = R. In the latter case, the main difference is that the input

set Q*i* is now also defined for a quantized signal, with Q*z* possibly defined similarly as

Q*o* (i.e. quantized).

Quantization have different effects in the CPU than in the ADC or the DAC. In many
applications, the round-off errors for floating-point operations performed by the CPU
*can be neglected. However, floating-point arithmetic can be sometimes replaced by *

*fixed-point arithmetic. The notation Qa.b denotes binary fixed-fixed-point numeric representation*

*using a + b + 1 bits, with a integer bits, b fractional bits and 1 sign bit. For example,*

*Q13.18 denotes 32-bit binary fixed-point numeric representation, with 13 integer bits*

plus 1 sign bit and 18 fractional bits. Fixed-point arithmetic operations are performed in general more efficiently (in terms of energy and speed) than floating-point operations. The main drawback of fixed-point arithmetic is an increase in round-off errors due to the lower precision when compared to floating-point arithmetic. See [58] for a discussion of fixed-point arithmetic used in model predictive control.

In contrast to the CPU, quantization effects play a major role on the DAC and ADC,
because these components use shorter word sizes (commonly between 8 and 16 bits). The
*maximum quantization error β may be large enough to noticeably degrade performance*
in cost-sensitive applications, where 8-bit or 10-bit DACs and ADCs are typically used.

**Sampling period**

The selection of the sampling rate is commonly based on a rule of thumb and a trade-off
*between cost and performance. The first thing to consider is that the sample rate ωs*

the Nyquist-Shannon sampling theorem, which states that a signal that is a bandlimited
*by ωb* *is uniquely characterized by its uniform samples taken with a frequency ωs* *>2ωb*,

see [54, Ch. 16]. In practice, higher ratios are required to have a smooth time response. The following rule of thumb is commonly used [55]:

*20 <* *ωs*

*ωb*

*<40.*

In [54, Ch. 16] other possible heuristic rules are also discussed.

The next consideration in determining the sampling rate is cost. In general, higher
sampling rates imply better closed-loop performance. If cost is no concern, a powerful
processor can be used that allows a very high sampling rate (e.g. *ωs*

*ωb* *>* 40). For

cost-sensitive applications, the sampling rate is selected as the lowest that meet all performance requirements: fast disturbance rejection, smooth time behavior, etc.

Although there exist applications that may benefit from a variable sampling period,
we focus here on systems with constant sampling period. Furthermore, we only deal with
continuous-time systems discretized via zero-order hold on the inputs, i.e. we consider
*discrete-time systems of the form xk+1* *= f(xk, uk*). As a side note, when using very

*fast sampling rates an incremental model of the form dxk* *= f(xk, uk, T*) [59] might be

necessary for the closed-loop control of embedded systems [58].

**2.1.2. Cyber-physical systems: real-time embedded systems**

*The term embedded system has different meanings in different contexts. In this work, we*
define an embedded system as a dedicated microprocessor-based system that controls a
*physical process. Recently, the term cyber-physical systems has been applied to describe*
embedded systems that control physical processes in real time [60, 61], possibly via
a computer network. This term has also different scopes in different disciplines, and
some authors emphasize the communication between embedded systems in the provided
definition [62]. We summarize the main characteristics of cyber-physical systems that
are relevant for later chapters.

**Hardware**

There exist a broad spectrum of microprocessors. We focus on devices at the bottom end of the spectrum in terms of computational capabilities, namely microcontrollers.

Even among microcontrollers there is a broad range of options. The CPU for example is driven by clock frequencies usually in the range of several tens to a few hundred megahertz. Another distinctive features is the word size, which is typically 8-bit, 16-bit or 32-bit. For computationally intensive applications, like the MPC algorithm, 32-bit microcontrollers are in general preferred. We briefly describe the most relevant hardware characteristics of microntrollers: the processor (or CPU) and the memory. In the case of the processor, we base our discussion on a popular family of 32-bit microprocessors: ARM Cortex-M. These state-of-the-art processors are nowadays in widespread use and are well suited for applications with complex data processing requirements like MPC [63, Ch. 1]. For a detailed description of this family of processors see [63].

We start by discussing the arithmetic capabilities of a microcontroller. Let us first consider the low-cost Cortex-M3 microprocessor. This CPU lacks a floating-point unit (FPU). This means that only integer arithmetic can be performed directly on hardware. Integer addition and multiplication can be executed in one clock cycle, whereas division needs several cycles. Furthermore, floating-point operations are emulated by software, with additions and multiplications being several times slower than their integer counter-parts. Empirical evidence shows that fixed-point operations (add, multiply) are around 4 times faster than floating-point operations for the Cortex-M3 [25] as well as for the classical ARM-7TDMI microprocessor [64]. In contrast, the more expensive Cortex-M4 includes digital signal processing (DSP) capabilities and, optionally, a single-precision FPU. This allows to execute 32-bit floating-point arithmetic operations directly on hard-ware.

The next big component of a microcontroller is memory. We distinguish between two types of memory: program memory, and data memory. The former is a non-volatile stor-age typically based on flash technology. It is commonly between ten kilobytes and a few megabytes, and is sometimes colloquially called ROM (read-only memory). The latter is a volatile type of storage typically based on SRAM (static random-access memory), and has a capacity of one kilobyte to a few hundred kilobytes.

**Real-time systems**

We now introduce some concepts related to applications that process data in real-time.
*We take most of our definitions from [60]. A real-time system is one whose correctness*
is based not only on the correctness of the results of the computations but also on the

physical time at which these results are produced.

*A deadline is the time by which a result must be produced. Deadlines can be classified*
in three types: firm, soft and hard. A firm deadline is characterized by the fact that a
result is no longer useful after the deadline has passed. If the deadline is soft, a result
that missed the deadline can still be used. A hard deadline is one that, if missed, failure
or catastrophic consequences may result. It follows that any hard deadline is a firm
deadline.

Similar as in the classification of deadlines, real-time systems are also separated in soft and hard real-time systems (also called safety-critical systems). The former is any real-time system with no hard deadlines, whereas the latter is a real-time system that must meet at least one hard deadline.

In this work, we consider digital control applications with firm deadlines. Depending
on the application, a control system may be soft or hard real-time. In the case of soft
real-time control systems, missing a deadline may result in degradation of performance.
*Another important concept is that of temporal determinism, which means that the*
*response time for any event is known. Computational delay is the time it takes the *
mi-croprocessor to execute the control algorithm. The variation on the computational delay
*is called the control jitter. We consider an algorithm to be (temporally) deterministic,*
if the control jitter is much smaller than the computational delay.

**Software**

Within software we make the distinction between application, operating system and compiler. We provide loose definitions of these terms, adapted from [65, Ch. 1]. An

*application performs an user-defined task. A digital control algorithm is considered an*

application. It may be a state observer, a model predictive control algorithm, etc.
*The operating system is in charge of interfacing the application with the hardware.*
*A real-time operating system (RTOS) is additionally in charge, among other things, of*
task scheduling. The control application is one of the many tasks that run concurrently
on the CPU. The RTOS scheduling algorithm guarantees, under certain conditions, that
our control algorithm is executed periodically with a constant sampling period [65, 66].
*In embedded applications the C programming language [39], or simply C, and C++*
are commonly used. However, C presents certain advantages over C++, namely: for
embedded hardware targets, object-oriented C++ may carry penalty on computational

performance and power consumption when compared to procedural C [67]. Furthermore, for embedded targets, C compilers are more common than C++ compilers.

*We use the term compiler as a broad term to denote the process of translating the user*
source code into binary code executable by the microcontroller. Hence, in this context a
compiler also includes a preprocessor, an assembler, and a linker. The source code refers
to both, the application and the operating system. In embedded systems, an application
is commonly written entirely in C or C++, whereas the OS is usually a mix of C/C++
and hardware specific assembler code.

**Memory management**

Although the topic of memory management is a complex one (see [65, Ch. 3]), here we reduce it to two simplified cases: dynamic and static memory allocation. A program-ming language typically offers only certain types of memory management functionality. Dynamic memory allocation is done during program execution. In C and C++ for in-stance, dynamic allocation is done with a call to the malloc function. In C++ the new operator can also be used to dynamically allocate memory. Furthermore, an RTOS is in charge of allocating memory for each task using the functionality provided by the base language.

The main advantage of dynamic allocation is that it offers flexibility. That is, the amount of memory required by certain types of data (e.g. arrays) does not need to be known at compile time. However, a dynamic memory allocation algorithm does not show temporal determinism [68, 69]. Therefore, its use is restricted in real-time systems. Moreover, the use of dynamic memory allocation is prohibited in critical systems accord-ing to the MISRA-C standard, [70]. Furthermore, some embedded C/C++ compilers do not implement dynamic memory allocation algorithms.

In contrast, static memory allocation is made at compile time. As its name implies, the amount of memory cannot be change during runtime. If the size of the data involved changes, the affected code needs to be recompiled. The main advantage is that there are no runtime penalties. Therefore, static memory allocation is commonly required in real-time applications.

**2.2. On-line optimization**

This section introduces some of the basic concepts of on-line optimization, which is fundamental for the implementation of MPC.

**2.2.1. Basics of convex optimization**

Before we proceed, let us state some required definitions. The following presentation is based on [31].

*We consider an optimization problem of the form*
minimize

*y∈Rny* *f*0*(y)*

*subject to fi(y) ≤ 0, for all i ∈ I,*

*fi(y) = 0, for all i ∈ E,*

(2.2)

*where the functions fi* are smooth real-valued functions in R*ny* and I and E are each a

*finite set of indices. The function f*0 : R*ny* *→ R is called the cost function or objective*

*function. The vector y ∈ Rny* *is the optimization variable. The functions f*

*i* : R*ny* →

*R, i ∈ I and fi* : R*ny* *→ R, i ∈ E are called inequality and equality constraint functions,*

*respectively. A point y is feasible if it satisfies all constraints. Problem (2.2) is said to*
*be feasible if there exist at least one feasible point. Otherwise it is called infeasible. The*
*set containing all feasible points is called the constraint set or feasible set. The optimal*

*value f*∗

0 of problem (2.2) is defined as:

*f*_{0}∗ *= inf{f*0*(x) | fi(x) ≤ 0, ∀ i ∈ I, fi(x) = 0, ∀ i ∈ E}.*

*A solution y*∗ _{to problem (2.2), also called an optimal point, is a point that is feasible}*and for which f*0*(y*∗*) = f*0∗ *holds. A feasible point y is called -suboptimal if satisfies*

*f*0*(y) ≤ f*0∗*+ , for > 0. We call the suboptimality level.*

*An optimization problem is called convex if the cost function is convex, the inequality*
constraint functions are convex, and the equality constraint functions are affine (i.e.

*fi(y) = Diy −di,for all i ∈ E). This implies that the feasible set of a convex optimization*

**2.2.2. Quadratic programs**

Of particular interest in this work are quadratic optimization problem (QP) of the form:
minimize
*y∈Rny*
1
2*y*
T
*Hy+ g*T
*y*
*subject to Cy ≤ c,*
*Dy= d.*
(2.3)

*With respect to (2.2) we have f*0*(y) =* 1_{2}*y*T*Hy* *+ g*T*y, fi(y) = Diy − di,* *for all i ∈ E,*

*and fi(y) = Ciy − ci, i ∈ I*. If we replace the affine inequality constraints in (2.3) with

quadratic constraints we obtain the problem:
minimize
*y∈Rny*
1
2*y*
T
*Hy+ g*T
*y*
subject to 1_{2}*y*T
*Piy+ qi*T*y ≤ ri,* *for all i ∈ I,*
*Dy* *= d,*

*with Pisymmetric positive semi-definite, the problem is called a quadratically constrained*

*quadratic program (QCQP). This type of problem arises in various MPC setups, see for*

*example [71, 72]. A related problem is the second-order cone program (SOCP) of the*
form:
minimize
*y∈Rny*
1
2*y*
T
*Hy+ g*T
*y*

subject to k ¯*Piy+ ¯pi*k2 ≤ *¯qi*T*y+ ¯ri,* *for all i ∈ I,*

*Dy* *= d.*

SOCPs arise in many engineering applications [73], including stochastic MPC problems [53, 74, 75]. We are mainly interested in QPs without equality constraints, which in standard form are represented as:

minimize
*y∈Rny*
1
2*y*
T
*Hy+ g*T
*y*
*subject to Cy ≤ c.* (2.4)

*In general, a QP is convex if and only if the Hessian matrix H is symmetric positive*
*semi-definite. Furthermore, we focus on strictly convex QPs. A QP is called strictly*
*convex if and only if H*T*= H > 0. Strict convexity implies that if a solution y*∗ _{exists}

then it is unique. Additionally, strict convexity is an important property exploited by many optimization algorithms, and in particular by the method discussed in Section 3.3.

**2.2.3. Optimization theory for quadratic programs**

One important concept is that of duality, since it commonly exploited by optimization
*algorithms. We start by defining the Lagrangian function of problem (2.2), which is*
*called the primal problem, as*

L*(y, λ) = f*0*(y) +*
X

*i∈E∪I*

*λifi(y),* (2.5)

*where λ denotes a Lagrange multiplier vector. Then, the Lagrange dual function is*
*defined as the minimum value of the Lagrangian over y for some fixed λ,*

*g*0*(λ) = inf*
*y∈Dy*
L*(y, λ) = f*0*(y) +*
X
*i∈E∪I*
*λifi(y).* (2.6)

One important characteristic of the dual function is that it is concave, even if problem
*(2.2) is not convex, [31, Ch. 5]. The Lagrange dual problem for problem (2.2) is defined*
as

*maximize g*0*(λ)*

*subject to λi* ≥*0, i ∈ I,*

(2.7)
*which is itself a convex optimization problem. We say that strong duality holds if the*
*optimal value of dual problem g*∗

0 *equals the optimal value of the primal problem f*0∗.
Duality is important since, for a strictly convex quadratic program (2.8) (for which strong
duality holds), the solution to the primal problem can be readily found by solving the
dual problem. Many optimization algorithms solve directly the primal problem, others
solve the dual problem instead to find the solution of the primal, or exploit information
of both primal and dual problems to find a solution.

We are interested in quadratic programs of the form (2.4). If the gradient and
*con-straint vector depend on a parameter w ∈ Rnw, we have the parametric problem P(w):*

minimize
*y∈Rny*
1
2*y*
T
*Hy+ g(w)*T
*y*
*subject to Cy ≤ c(w),* (2.8)

*with c(w) ∈ Rnc*. In MPC applications, this parameter is commonly given by (a

*com-bination of) the current state vector x, a reference trajectory, known disturbances. For*
notational simplicity and where no confusion may arise, the dependency on the
*param-eter w will be dropped.*

For our further discussion, we need a few further definitions.

**Definition 1** **(Feasible set). For the quadratic program P(w) the feasible set is given by**

F*(w) = {y ∈ Rny* * _{| Cy ≤ c}(w)}.* (2.9)

**Definition 2**

**(Active set). Consider problem P(w). A constraint C**iy ≤ ci(w), 1 ≤ i ≤*nc* *is called active at y ∈ F(w) 6= ∅ if Ciy* *= ci(w) holds, and inactive otherwise. The*

*set of indices*

A*(y; w) = {i ∈ {1, . . . , nc} | Ciy= ci(w)}* (2.10)

*is called set of active constraints or simply active set. If y*∗ _{is the solution to P(w) we}

*call A(y*∗_{; w) optimal active set.}

**Definition 3** **(LICQ). Given the pair (y, w), the constraints ˆc**i(y; w) = Ciy − ci(w),

*and the active set A(y; w), we say that the linear independence constraint qualification*

*(LICQ) holds if the set of active constraint gradients {∇ˆci(y; w), i ∈ A(y; w)}, is linearly*

*independent.*

Now we are ready to introduce the first-order optimality conditions, also known as Karush-Kuhn-Tucker (KKT) conditions (taken from [29, Theorem 12.1]), which charac-terize an optimal point.

**Theorem 1** * (First-order optimality conditions). Suppose that y*∗

*is a local solution to*

*(2.2), that the functions fi* *in (2.2) are continuously differentiable, and that the LICQ*

*holds at y*∗* _{. Then there is a Lagrange multiplier vector λ}*∗

*∗*

_{, with components λ}*i, i ∈ E ∪I,*

*such that the following conditions are satisfied at (y*∗* _{, λ}*∗

_{):}∇L*(y*∗*, λ*∗*) = 0,*
*fi(y*∗*) ≤ 0, ∀ i ∈ I,*
*fi(y*∗*) = 0, ∀ i ∈ E,*
*λ*∗* _{i}* ≥

*0, ∀ i ∈ I,*

*λ*∗

*∗*

_{i}fi(y*) = 0, ∀ i ∈ I ∪ E.*(2.11)

A proof can for example be found in [29, Ch. 12].

For convex problems like (2.8), the KKT conditions provide sufficient and necessary
*conditions for a pair (y*∗* _{, λ}*∗

*∗*

_{) to be optimal. The point y}_{is the unique minimizer if the}

*problem is strictly convex. Furthermore, if the LICQ holds, for any given y*∗

_{the optimal}

*Lagrange multiplier λ*∗

_{is also unique.}

For the discussion of the active set optimization algorithm and the explicit MPC method in Chapter 3, the following definitions are needed.

**Definition 4****(Set of feasible parameters). The set of feasible parameters of a parametric**

*quadratic program is given by*

W *= {w ∈ Rnp* _{| F}*(w) 6= ∅}.* (2.12)
**Definition 5****(Critical region). Let a strictly convex quadratic program P(w) with w ∈ W**

*be given. Furthermore, let y*∗* _{(w) denote the solution and A(y}*∗

_{(w); w) the corresponding}*optimal active set. Then, for every index set A ⊆ {1, . . . , nc*}*, the critical region of W*

*is given by*

*CR*A *= {w ∈ W | A(y*∗*(w); w) = A}.* (2.13)

**2.3. Model predictive control for linear systems**

We are ready to introduce the model predictive control algorithm. We consider the
linear MPC setup of the form (1.1). Figure 2.2 depicts the behavior of a system being
*controlled by a discrete MPC implementation. At time tk* *= kT , k ∈ Z and T the*

*sampling period, we estimate the state ˆxkfrom measurements ˆyk*, typically using a state

*observer. The MPC algorithm uses ˆxk* *and the linear model of the system to predict*

**the future trajectory of the state x***k* *for the following N steps. With this information,*

**a minimization problem is solved to find the input sequence u***k* that optimizes a given

*performance criterion for the given horizon N. In the figure, the continuous function*

*φ(t, ˆxk , uk) (shown in dashed red) represents the nominal state evolution starting at*

*ˆxk* **under the influence of the input u***k***. The predictive state trajectory x***k* coincides

*with φ(t, ˆxk , uk) at each sampling time tk*

*(represented by black dots over φ(tk, ·*)). The

*discrete-time MPC algorithm has no knowledge of φ(·), though. We have included this*
continuous time plot in the figure to simplify explanation of the concept. Furthermore,

*t*
*x, u*
*tk* *tk+1* *tk+2* *tk+N −1* *tk+N* *tk+N +1*
ˆ
*xk*
ˆ
*xk+1*
ˆ
*xk+2*
*φ(t,ˆxk ,uk*)

*φ(t,ˆxk+1*)

**,u**k+1**u**

*k*

**u**

*k+1*

*. . .*

**x**

*k*

**x**

*k+1*

*N*

*N*

Figure 2.2.: A visual representation of the MPC algorithm.

**u***k*is held constant at each sampling interval (a zero-order hold), as is common for digital

controllers [7].

Ideally, the state trajectory of the controlled system would follow the trajectory

*φ(t, ˆxk , uk*). However, due to several effects which are often not included in the MPC

model (disturbances, model-plant mismatch, errors in the state estimation,
quantiza-tion of the input, etc.), the real system follows a slightly different trajectory under the
**influence u***k* *during the interval from tk* *to tk+1* (shown in solid red).

*At the next sampling period tk+1*, we repeat the process. We first take a new

*mea-surement ˆyk+1* *and estimate ˆxk+1***. We then compute the input sequence u***k+1* (dashed

* blue) for the next following N steps, and apply the first part of uk+1* (solid blue) to the

real system during the current sampling interval, and so on for the subsequent sampling
*times tk, k = 2, 3 . . .*

To simplify our discussion, we focus on one particular moving horizon setup, namely regulation to the origin of a constrained linear system. We describe this sample setup in the following.

**2.3.1. A basic MPC setup**

We considered discrete-time linear time-invariant systems subject to input and state constraints described by

*x*+*=Ax + Bu,* (2.14)

*where x ∈ X ⊆ Rn, and u ∈ U ⊂ Rm* are the system state and input at the current

*sampling time, respectively. The successor state is denoted by x*+. The considered input
*constraint set U is a convex and compact box set containing the origin. Moreover, the*
*state constraint set X is a closed (not necessarily bounded) convex set with the origin in*
its interior. These type of constraints are commonly found in practice (see [3, Ch. 1]).

*The discrete-time system and input matrices are represented by A and B, respectively.*
*We assume the pair (A, B) to be stabilizable. The control objective is to bring the system*
to a desired equilibrium point while satisfying all constraints. For simplicity and without
loss of generality, in this section we only consider regulation to the origin.

MPC is a natural candidate to solve the constrained regulation problem. The MPC
*controller needs to solve at each sampling time k an optimization problem parameterized*
*by the current state x. This problem is, for the conditions given above, defined as follows:*

minimize** _{u}** 1

_{2}

*N −1*X

*j=0*

*(kxj*k2

*Q+ kuj*k2

*R*) + 1 2

*kxN*k2

*P*

*subject to xj+1*

*= Axj+ Buj,*

*j*

*= 0, ..., N − 1,*

*u ≤ uj*

*≤ u,*

*j*

*= 0, ..., N − 1,*

*e ≤ Exxj*

*+ Euuj*

*≤ e,*

*j*

*= 0, ..., N − 1,*

*f*

*≤ F xN*

*≤ f ,*

*x*0

*= x,*(2.15)

*where N ≥ 2 is an integer denoting the prediction horizon. We consider positive definite*
*quadratic forms kwk*2

*S* *= w*

T_{Sw >}*0 and kwk*2

*S* *= 0 only if w = 0, with S a symmetric*

*positive definite matrix. The matrices Q, R and P are symmetric positive definite and*
*denote the states, input and terminal weighting matrices, respectively. If Q is symmetric*
positive semi-definite, an additional requirement for stability is needed, namely that the
*pair (A, Q) is detectable [3, Ch. 2]. The matrices Q and R, loosely speaking, are a*

*specification of the desired controller behavior. The term `(xj, uj*) = 1_{2}*(kxj*k2*Q+ kuj*k2*R*)

*is called stage cost, whereas `N(xN*) = 1_{2}*kxN*k2*P* *is called terminal cost. Furthermore,*

*the box set U is defined by the lower and upper bounds u and u. Similarly, the state*
*constraints X are defined by the bounds e and e, and the matrices Ex* ∈ R*q×n* and

*Eu* ∈ R*q×m. The terminal state constraint set Xf* *⊆ X* is a polytope defined by the

*triplet f, f and F ∈ Rr×n*_{. This set, together with the cost function, is often used to}

guarantee closed-loop stability [3, Ch. 2].

* We optimize over the input sequence u = {u*0

*, . . . , uN −1*} ∈ U

*= UN*. The state

**constraints further impose the constraints u ∈ U***N (x) = {u ∈ U | x(x, u) ∈ X }.*

**The state sequence x(x, u) ∈ R**(N +1)n_{describes the trajectory followed by the states of}
* system (2.14) starting at x and under the input sequence u, i.e. x(x, u) = {x*0

*= x, x*1 =

*Ax*0 *+ Bu*0*, . . . , xN* *= AxN −1* *+ BuN −1*}**. The set X = {x ∈ R***(N +1)n* *| xk* *∈ X, k* =

*0, . . . , N − 1, xN* *∈ Xf*} denotes the state trajectories satisfying all constraints. Thus,

U*N(x) is the set of feasible input sequences of problem (2.15) for a given x.*

**In an MPC scheme, the first element of the optimal input sequence u**∗ _{is used as}
*feedback control action, i.e. during each sampling period u = u*∗

0 is applied to the plant.
In many applications, it is safe to assume that the input sequence found by optimization
**algorithms is indeed u**∗_{. However, in the case of embedded systems with fast sampling}
**rate, the MPC controller may only have enough time to compute u**∼ _{∈ U}

*N(x), a (possibly*

**rough) approximation to u**∗_{.}

**2.3.2. MPC as a QP**

Most implementations of optimization algorithms cannot directly solve problems like (2.15). Therefore, the MPC problem needs to be brought into a form that can be directly used together with a tailored solver.

The MPC optimization problem (2.15) can be brought into the following equivalent form (refer to Appendix A)

minimize
**u∈R***N m*
1
2**u**
T
* Hu+ g(x)*T

**u**

*(2.16)*

**subject to Cu ≤ c(x).***Compared to (2.4), we see that the gradient vector g(x) is now given by a function*

*particular c : Rn* _{→ R}*2((m+q)N +r) _{. Problem (2.16) is known as parametric QP, and}*

*is parameterized by the variable x. Note that this problem does not include equality*constraints (i.e. E = ∅).

There are several ways to guarantee strict convexity of (2.16). In the case of setpoint
*stabilization (2.15), selecting the (symmetric) weighting matrices such that Q ≥ 0, P ≥ 0*
*and R > 0 guarantees that the Hessian matrix is symmetric positive definite.*

A MPC problem like (2.16) can be reliably solved by many commercial (e.g. MAT-LAB, Gurobi) and non-commercial (e.g. CVXOPT [76], OpenOpt [77]) optimization packages [78]. Nevertheless, to efficiently apply MPC on embedded systems, a solver needs to exploit the structure of problems like (2.16) that are derived from (2.15). Next chapter addresses this topic.

**2.3.3. General moving horizon control formulation**

So far, we have only consider one particular type of MPC problem. However, one of the
strengths of moving horizon formulation is that is general and therefore can deal with
a great variety of problems. Nevertheless, we focus our attention on problems that can
**be equivalently expressed as a parametric quadratic program P(p) of the form:**

minimize** _{u}** 1

_{2}

**u**T

* Hu+ g(p)*T

**u**

* subject to Cu ≤ c(p),* (2.17)

**with H and C constant matrices, and g(p) and c(p) vectors affine in each of the ****param-eters of the parameter sequence p.**

Note that (2.16) is the standard QP representation the MPC setup (2.15). This particular problem (regulation to the origin subject to constraints) has the state vector

*x* **as its only parameter. That is, (2.16) is a particular case of (2.17) in which p consist**

* of a single parameter x, i.e. p = {x}. In general, p may consist of several parameters*
that depend on the MPC problem at hand. Thus, formulation (2.17) allows us to
consider a broad range of problems, including, but not limited to: setpoint stabilization,
trajectory tracking [79, 80], systems with dead time [81], and systems with known-ahead
disturbances [11, 82]. Notable systems that cannot be expressed in form (2.17) are
time-varying systems, nonlinear systems, and formulations with weighting matrices time-varying
with time.

**Require:** **At each sampling time k = 0, 1, . . . get the parameter sequence p**k

1: **Find the solution u**∗ **to problem (2.17) for p = p***k*

2: *Assign the MPC control action u _{k}*

**= u**∗

_{0}3:

**return u**k**Algorithm 1:** Model predictive control algorithm

The MPC scheme that was intuitively described in Figure 2.2 is formally presented in Algorithm algorithm 1.

**2.4. Summary**

We presented a general description of the characteristics of the problems we are inter-ested in. We introduced digital control of embedded systems and outlined important concepts that are relevant for MPC applications. We focused first on the limitations of embedded computers. We then introduced the model predictive control algorithm for linear systems, and showed that it is, under certain conditions, equivalent to a convex optimization problem. We presented the theoretical foundations of convex optimization which are needed in the following chapters. In particular, we discussed convex quadratic programs, a type of optimization problem commonly arising in linear MPC. In the next chapter we discuss how problems like (2.17) can be solved efficiently, while at the same time taking the limitations of embedded systems into account.

**software tools for MPC**

This chapter discusses how optimization algorithms can be tailored for efficient MPC implementations. We start by describing the main properties of the MPC algorithm that need to be taken into account. It follows the discussion on how some of these properties are exploited by state-of-the-art tailored optimization software tools. We then present a novel optimization algorithm that takes into account the properties of MPC and is well suited for embedded system. This allows the algorithm to achieve good controller performance with the limited computational resources available on embedded hardware platforms.

**3.1. Exploiting the properties of the MPC algorithm**

In this section we discuss the most relevant properties of the MPC algorithm that need to be taken into account in embedded systems. These properties can be used to develop efficient MPC implementations. We discuss some methods that work well in practice. From a theoretical perspective, however, they present some challenges. We, therefore, emphasize in the following the practical implications of these ideas, while we only briefly mention some theoretical aspects.

**3.1.1. Known-ahead maximum computation time**

As discussed in Section 2.1, one fundamental assumption that is commonly made in digital control systems is that the sampling period is constant. This imposes a limit on what the maximum computational time of the MPC algorithm can be. Because an MPC scheme runs alongside with several other software components on the same processor, the exact available computational time for the MPC can in most cases only determined