Towards Parallelizable Sampling–based
Nonlinear Model Predictive Control
R.V. Bobiti ∗ M. Lazar ∗

arXiv:1701.02660v2 [] 12 Jan 2017

∗

Department of Electrical Engineering, Eindhoven University of Technology,
The Netherlands (e–mails: r.v.bobiti@tue.nl, m.lazar@tue.nl)

Abstract: This paper proposes a new sampling–based nonlinear model predictive control (MPC)
algorithm, with a bound on complexity quadratic in the prediction horizon N and linear in the number of
samples. The idea of the proposed algorithm is to use the sequence of predicted inputs from the previous
time step as a warm start, and to iteratively update this sequence by changing its elements one by one,
starting from the last predicted input and ending with the first predicted input. This strategy, which
resembles the dynamic programming principle, allows for parallelization up to a certain level and yields
a suboptimal nonlinear MPC algorithm with guaranteed recursive feasibility, stability and improved
cost function at every iteration, which is suitable for real–time implementation. The complexity of the
algorithm per each time step in the prediction horizon depends only on the horizon, the number of
samples and parallel threads, and it is independent of the measured system state. Comparisons with
the fmincon nonlinear optimization solver on benchmark examples indicate that as the simulation time
progresses, the proposed algorithm converges rapidly to the “optimal” solution, even when using a small
number of samples.
Keywords: suboptimal nonlinear predictive control, sampling–based optimization, embedded model
predictive control, dynamic programming, constrained control
1. INTRODUCTION
A vast literature on nonlinear model predictive control (NMPC)
has proven both the theoretical (Grüne and Pannek, 2011), as
well as the practical advantage (Magni et al., 2009) of this
method in treating optimally the control of multi–variable nonlinear systems subject to constraints on state and inputs. Common research interests include methods for reducing the complexity of the NMPC algorithms, to make them applicable on
devices of low memory (ASIC, FPGA), such as explicit NMPC
(ENMPC), see, e.g., (Johansen, 2004). Much work has been
done on treating other limiting factors, such as the requirement
of NMPC of solving an optimization problem online. This is
not well achieved by common optimization tools, which have
no specific termination time, especially due to non–convexity
which involves multiple local–minima. Therefore, real time
requirements are not met, which limits the industrial impact of
NMPC. To alleviate these concerns, the literature has proposed
multiple solutions, such as approximate dynamic programming
(DP) (Bertsekas, 2005a), suboptimal MPC (Scokaert et al.,
1999), approximative DP and suboptimal control via rollout
algorithms (Bertsekas, 2005b), NMPC based on system approximation by neural models (Ławryńczuk, 2009).
An alternative strategy in NMPC is to draw samples from either
the state or input space, to design computationally feasible
NMPC methods, see for example, (Piovesan and Tanner, 2009)
which proposes a randomized approach to sampling the space
of predicted input sequences. More recently, in (Chakrabarty
et al., 2016), an ENMPC method was proposed based on
sampling of the state space for continuously differentiable
nonlinear systems. The method therein solves optimization
problems offline to find optimal control sequences, which are
used to construct the ENMPC strategy. While there are still

concerns in ENMPC related to robustness, feasibility of the
offline optimization and finding the neighbors in the sampled
grid for the off–grid states, ENMPC, when successful, reduces
significantly the computational load of MPC at the expense of
an acceptable cost degradation.
Input and state space sampling–methods for solving NMPC
via approximate DP have also been proposed, see (Bertsekas,
2005a), though they inherit the dimensionality issues of DP
(Lee and Lee, 2004).
Another relevant sampling–based strategy, the so–called sampling based MPC (SBMPC), was proposed in (Dunlap et al.,
2010). The method therein is applicable to nonlinear systems
in general, though, its performance is dependent on a user–
specified heuristic. A more ample discussion on sampling–
based DP and SBMPC, in the light of the method proposed in
this paper is reported in Section 2.3.
A common problem of sampling–guided methods for NMPC is
the sampling strategy. For example, with each input sample, a
tree is expanded. After the tree is built, the path of least cost in
the tree is selected from the initial state to the desired state.
If the sampling is performed over the input space, and each
sample is connected to all the samples in the input space for
the next time step in the control horizon, then the tree growth
is exponential with the horizon. Alternatively, as in randomized
MPC (Piovesan and Tanner, 2009), sampling randomly in the
input space, of dimension m, augmented to the horizon of
dimension N requires a large number of samples, in an mN
dimensional space, to achieve a significant accuracy.
In this paper we adopt a suboptimal formulation of NMPC,
as originally proposed in (Scokaert et al., 1999), where it
was shown that feasibility of a solution implies stability under

suitable conditions. This, together with the fact that suboptimal
NMPC has the same inherent robustness properties as optimal
NMPC, see (Pannocchia et al., 2011) and (Lazar and Heemels,
2009), suggest that suboptimal NMPC is a viable and in fact the
best one can hope for when a sampling–guided MPC strategy is
undertaken for the control of nonlinear systems. Furthermore,
we aim at a sampling method which provides a suboptimal
solution that yields good control performance, has a reasonable
computational complexity increase with the prediction horizon
and allows for parallel implementation up to some level.
In this paper, the main idea for achieving this goal is to use
the shifted sequence of predicted inputs from the previous time
step as a warm start, and to iteratively update this sequence
by changing its elements one by one, starting from the last
predicted input and ending with the first predicted input. This
strategy resembles the dynamic programming principle, especially the rollout algorithms, see (Bertsekas, 2005b), which improve a heuristic base policy for optimal control. Additionally,
in this paper, we sample the original input space, which is
typically represented by a proper set U ⊂ Rm . Sampling allows
for parallelization of the calculations performed for updating
each of the elements of the predicted sequence of inputs and
it enables limiting the computational time according to the
requirements of the considered application. An upper–bound
on the complexity of the overall algorithm is quadratic with the
prediction horizon N and linear with the number of samples in
U. This enables the usage of long prediction horizons or real–
time implementation on inexpensive computing devices such as
ASIC and FPGA. The suitability for real–time implementation
is also enhanced by the fact that the algorithm can be stopped
at any iteration performed within a sampling period, while the
complexity of the calculations per iteration depends only on
the horizon N , the number of samples and parallel threads,
and it does not depend on the measured state of the system.
Moreover, the updated predicted sequence of inputs obtained
at any iteration will guarantee recursive feasibility, stability
and an improved cost function under the same conditions as
in suboptimal NMPC (Scokaert et al., 1999).
The remainder of this paper is organized as follows. In Section 2, basic notation is introduced, together with the problem
formulation and a discussion on the relation with the existing
methods. Section 3 presents the main result as a prototype algorithm and its complexity analysis. In Section 4, three examples
illustrate the potential of the developed method, and Section 5
concludes the paper.
2. PRELIMINARIES
2.1 Notation
Let R, R+ , Z and Z+ denote the field of real numbers, the set
of non–negative reals, the set of integers and the set of non–
negative integers, respectively. For every c ∈ R and Π ⊆ R,
define Π≥c := {k ∈ Π | k ≥ c} and similarly Π≤c . Let int(S)
denote the interior of a set S. Let Sh := S × . . . × S for any
h ∈ Z≥1 denote the h–times Cartesian–product of S ⊆ Rn . A
set S ⊂ Rn is called proper if it is non–empty, compact and
0 ∈ int(S).
For a vector x ∈ Rn , the symbol kxk is used to denote
an arbitrary p–norm; it will be made clear when a specific
norm is considered. For a vector x ∈ Rn , define |x| :=
[|x1 | . . . |xn |]T . Also for a vector x ∈ Rn , by max |x| we denote

max{|x1 |, . . . , |xn |}. For a scalar x ∈ R, denote by dxe the
smallest integer number larger than x.
A function α : R+ → R+ is said to belong to class K, i.e.,
α ∈ K, if it is continuous, strictly increasing and α(0) = 0.
Furthermore, α ∈ K∞ if α ∈ K and lims→∞ α(s) = ∞.
2.2 Suboptimal MPC problem formulation
Let us consider the discrete–time system described by
xk+1 = f (xk , uk ),
n

(1)

m

where xk ∈ R is the state and uk ∈ R is the control vector
at discrete–time k ∈ Z+ . We assume that the map f : Rn ×
Rm → Rn satisfies f (0, 0) = 0, which is, the origin is an
equilibrium point for system (1).
The goal of MPC is to regulate the state to the origin while
satisfying control and state constraints, i.e., uk ∈ U ⊂ Rm and
xk ∈ X ⊂ Rn for all k ∈ Z+ , where U and X are proper
sets. MPC relies on a receding–horizon control law in order to
determine, for each k, a finite–sequence of control inputs
U (k) = {uk|k , uk+1|k , . . . , uk+N −1|k },

where N is the control and prediction horizon, which are
considered equal in this paper, for simplicity of exposition. If
the initial state is xk|k = xk and the control sequence is U (k),
the solution of system (1) in closed-loop with U (k) at time k +i
is denoted by φ(xk|k , U (k), i). The current control action uk , is
selected as the first control action in U (k), i.e., uk = uk|k .
To achieve this, optimal MPC minimizes, at each discrete–time
k, a cost function of the type
J(xk|k , U (k)) = Vf (xk+N |k ) +

N
−1
X

L(xk+i|k , uk+i|k ), (2)

i=0

where J : Rn × Rm → R+ is the total cost function, Vf :
Rn → R+ is a terminal cost and L : Rn × Rm → R+ is a stage
cost. The minimization is performed with respect to U (k) and
it is subject to
xk+j|k = f (xk+j−1|k , uk+j−1|k ),

∀j ∈ Z[1,N ] ,

(3)

and to the state and input constraints.
When the dynamics f is a nonlinear, possibly non-convex
function, the optimization of the cost (2) cannot be guaranteed
to converge to a global optimum, in general. Alternatively,
suboptimal MPC was proposed as a viable alternative, see, e.g.,
(Scokaert et al., 1999), to deal with this inherent shortcoming
of nonlinear global optimization.
Assume that the following constraints are required to hold at
each iteration of the MPC problem:
xk+i|k ∈ X, uk+i|k ∈ U,

∀i ∈ Z[0,N −1] ,

(4)

and
xk+N |k ∈ XT ,
(5)
where XT ⊆ X is a proper set which represents a terminal
constraint. Moreover, define by U(xk|k ) the set of control
sequences U (k) which, applied on xk|k , satisfy (3), (4) and (5).
Suboptimal MPC relies on an initial feasible solution, a warm
start sequence Uwarm (k) ∈ U(xk|k ) at each step k, which is
improved iteratively. The suboptimal MPC problem considered
in this paper is formulated as follows:

Problem 2.1. For each k ∈ Z≥1 , given Uwarm (k) find a
sequence U (k) ∈ U(xk|k ) such that
J(xk|k , Uwarm (k)) > J(xk|k , U (k)),
and the constraints (3), (4) and (5) are satisfied.

(6)

Consider now a locally stabilizing control law kf : XT → U.
Assume that XT is a sublevel set of Vf . For stability of the MPC
closed–loop system it is also required that:
• Vf (f (x, kf (x))) + L(x, kf (x)) ≤ Vf (x) for all x ∈ XT ;
• there exist α1 , α2 ∈ K∞ such that α1 (|x|) ≤ Vf (x) ≤
α2 (|x|) for all x ∈ XT ;
• there exist α3 ∈ K∞ such that L(x, u) ≥ α3 (|x|) for all
(x, u) ∈ X × U;
Remark 2.2. The first property above listed implies that XT is
positively invariant for the system xk+1 = f (xk , kf (xk )). The
second and third properties can be satisfied if, for example Vf
and L are positive definite quadratic functions.
These properties imply that the cost function J(·, ·) is a Lyapunov function, see (Mayne and Rawlings, 2009, Lemma 2.14).
As such, if F is the set of states in X for which there exists a control sequence U (k) which satisfies the constraints (3), (4) and
(5), then the solution to Problem 2.1 provides an asymptotically
stabilizing controller with a region of attraction F.
In (Mayne and Rawlings, 2009), an algorithm was proposed for
suboptimal MPC with stability guarantees. Given the current
state xk|k and the previous control sequence U (k − 1) as an
input, the steps of the algorithm therein can be summarized as
follows:
• If xk|k ∈
/ XT , use the warm start sequence:
Uwarm (k) ={uk|k−1 , uk+1|k−1 , . . . ,
uk+N −2|k−1 , kf (xk+N −1|k−1 )}.
(7)
Solve iteratively Problem 2.1 via optimization to improve
Uwarm (k) with a U (k) ∈ U(xk|k ). Apply control uk =
uk|k .
• If xk|k ∈ XT , set uk|k = kf (xk|k ), or, similarly to the previous case, use the warm start and solve an optimization
algorithm iteratively to find an improved control sequence
U (k) ∈ U(xk|k ).
Optimization solvers, in both optimal and suboptimal MPC,
present difficulties in terms of parallelization and a priori
known execution time independently of the current state xk .
To circumvent these drawbacks and enable a computationally
efficient and parallelizable nonlinear MPC algorithm, we propose a sampling based approach to solving Problem 2.1, as
illustrated in Section 3. To this end, we first review existing,
similar approaches within nonlinear MPC.
2.3 Existing approaches based on sampling
This section provides a brief in depth review of two main
existing approaches for NMPC based on sampling, namely
approximate DP and sampling–based MPC, which were also
mentioned in the Introduction.
An approximate version of DP, as a tool for solving optimization problems, was proposed in (Bertsekas, 2005a, Section
6.6.1). DP has been successfully applied for determining explicit solutions for linear MPC controllers, see, e.g., (Muñoz
de la Peña et al., 2004). In nonlinear MPC, the state and input

states are typically discretized to apply DP algorithms. The
main idea is to discretize the state space with a finite grid and
to express each state outside of the grid as an interpolation of
nearby grid elements. The same interpolation law is applied to
compute the cost of the current nongrid state as a function of the
costs of the nearby grid states. As such, by discretizing both the
state space for each time in the control horizon and the input
space, a transition diagram is obtained which approximates
the dynamics of the system in the continuous space. On this
discrete transition system, DP is applied to determine the path
with the smallest cost, which, for a given initial state xk|k ,
provides the control sequence U (k).
Solving MPC with DP via discretization suffers from the “curse
of dimensionality”, due to sampling of both state and input
spaces and the requirement for constructing the complete transition diagram, by evaluating the subsequent state and cost for
each sampled state and all the samples in the input space.
An alternative to approximate DP, namely sampling based
MPC, was developed within the area of Robotics, where typically optimization problems arising in control are non–convex,
due to either kinematic constraints or constraints posed by
obstacle avoidance. Sampling–based motion planning such as
rapidly-exploring random trees (RRTs) (LaValle, 1998) or randomized A* algorithms (Likhachev and Stentz, 2008), have
been extensively used to construct trees which connect an initial
state to a final state based on sampling states in the search
space and searching for feasible inputs to connect these states.
To approach issues related to the, possibly unfeasible, search
for an input after sampling only in the state space, a method
(SBMPC) was proposed in (Dunlap et al., 2010), which samples
the input space at each sampling period and creates trees that
contain feasible state trajectories. The optimal path to a goal in
the state–space is then searched for within the tree using goal–
directed search algorithms, such as LPA*. Such algorithms rely
on computing a heuristic measure of the distance from the
current sample to the goal. Selecting the heuristic is however,
not an obvious task for general nonlinear systems.
Therefore, a desirable feature of a sampling–based suboptimal
NMPC algorithm is a non–exponential growth in the tree generated through sampling of the state or input space at each
step in the horizon. Furthermore, it is also desirable to reduce
the dependency of the algorithm on the non–obvious selection
of a heuristic, which significantly impacts the performance of
the sampling–based strategy. To circumvent these issues, an
alternative suboptimal strategy for sampling the input space is
proposed in the next section, based on sequentially updating a
warm start feasible sequence of predicted inputs.
3. MAIN RESULT
In this section we present the proposed sampling–based algorithm for solving suboptimal NMPC problems and we provide
a detailed complexity analysis.
3.1 Prototype algorithm
Inspired by the suboptimal MPC and DP principles, we propose a sampling–based solution to Problem 2.1, which, by the
mechanism involved in the iterative improvement of the initial
feasible control sequence Uwarm (k), has a low increase with
the control horizon of the computational complexity.

j=0

j=1

j=2
uw
2|0

uw
0|0

u12|0
u22|0

uw
1|0

u32|0

j=1

j=0
u41|0

uw
0|0

u11|0

j=2

u31|0

uw
1|0

uw
2|0

=

u32|0

u21|0

j=1

j=0

u10|0

uw
1|0

j=2

= u41|0

uw
0|0 u2
0|0
u30|0

u40|0

3
uw
2|0 = u2|0

Fig. 1. Sampling principle.
Algorithm 1 Sampling–based suboptimal NMPC algorithm.
Input: N , {nj }j∈Z[0,N −1] , xk|k , X, U, XT , J(·, ·)
w
Uwarm (k) = {uw
k|k , . . . , uk+N −1|k }
Output: U (k), uk
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:

Jsub ← J(xk|k , Uwarm (k));
for all j = N − 1 : −1 : 0 do
Select nj samples uqj+k|k ∈ U, q ∈ Z[1,...,nj ] ;
for all q = 1 : 1 : nj do
if j ≥ 1 then
q
w
U (k)jq = {uw
k|k , . . . , uk+j−1|k , uk+j|k ,
w
uw
k+j+1|k , . . . , uk+N −1|k };
else
w
U (k)jq = {uqk|k , uw
k+1|k , . . . , uk+N −1|k };

if φ(xk|k , U (k)jq , i) ∈ X, ∀i ∈ Z[j+1,N −1] and
φ(xk|k , U (k)jq , N ) ∈ XT then
Jnew ← J(xk|k , U (k)jq );
if Jnew < Jsub then
Jsub ← Jnew ;
Uwarm (k) ← U (k)jq ;
U (k) = Uwarm (k), uk = uw
k|k ;

The principle behind the proposed sampling–based approach is
illustrated in Fig. 1 and formalized in Algorithm 1. In Fig. 1, the
iterative improvement of an initial cost provided by an initial
feasible control sequence Uwarm (k) is illustrated for the case
when N = 3. The algorithm keeps always Uwarm (k) as a
reference sequence, and it covers the horizon in a backward
fashion, in N iterations. Starting with j = N − 1, at each
iteration step, nj samples {uqk+j|k }q∈Z[1,nj ] are drawn from the
input constraint set U. For each sample, the reference sequence
Uwarm (k) is modified in the j th location, and a new sequence,
U (k)jq is obtained. If the following properties hold:

• U (k)jq is a feasible sequence, i.e., the constraints (3), (4)
and (5) hold,
• the new cost, Jnew = J(xk|k , U (k)jq ) decreases with
respect to J(xk|k , Uwarm (k)),
then Uwarm (k) is replaced by U (k)jq and the algorithm continues backwards with respect to the prediction time j, in a
similar manner. With this approach, at any point in time, if the
maximally allowed computational time is exceeded, a feasible,
improved control sequence exists and it can be utilized as a
suboptimal MPC solution.
By choosing to cover the control horizon backwards, we can
reuse at each step j, the states φ(xk|k , Uwarm (k), i) for all
i ∈ Z[1,j] , which, by the feasibility of Uwarm (k), already
satisfy the state and input constraints. This holds not only for
the original Uwarm (k), but for any subsequent improvement
of Uwarm (k). As such, also the stage costs up to the j th state
can be recovered from previous computations. The reusability
of previously computed costs and states by navigating the
control horizon in a backward manner resembles the working
principle of DP. This suggests intuitively that the proposed cost
improvement method could deliver good performance, which
is supported by results obtained in non–trivial case studies, see
Section 4. A formal analysis of convergence towards the DP
solution, as the discrete–time k increases, makes the object of
future research.
When k = 0, we can select an initial sequence Uwarm (0)
as the solution of the optimal MPC problem. In this case, we
can proceed with k = 1. Alternatively, we can select a feasible sequence Uwarm (0), an “oracle”, by randomly selecting
sequences of inputs in U until a feasible Uwarm (0) is found.
In this case, if it is feasible to afford the computational time,
we can proceed with Algorithm 1 in an attempt of obtaining an
improved sequence.
For k ∈ Z≥1 , to choose the input Uwarm (k) for Algorithm 1,
one can use the receding horizon principle of MPC. As such,
the input sequence Uwarm (k) = {uk|k−1 , . . . , uk+N −2|k−1 , u}
is a warm start at time k. If φ(xk−1|k−1 , U (k − 1), N ) ∈
XT and XT is positively invariant for the system xk+1 =
f (xk , kf (xk )), one can select u = kf (xk+N −1|k−1 ). In this
case, Uwarm (k) is a feasible solution, and therefore a candidate
warm start for every k ∈ Z≥1 . If there exists no terminal set XT
and no kf (·), then one can select u ∈ U such that Uwarm (k)
remains feasible, i.e., φ(xk−1|k−1 , Uwarm (k), N ) ∈ X. In these
circumstances, however, stability of the closed loop is not
guaranteed. Such an example is illustrated in Section 4.3.
Remark 3.1. Common sampling schemes are employed for
sampling of the input space U at each iteration, among which
we consider deterministic uniform sampling, which places each
sample at equal distance from each other, to cover uniformly the
space U. Alternatively, “true” random samples can be selected,
which are simpler to draw in higher dimensional spaces. Quasi–
random low–discrepancy sequences, see (Chakrabarty et al.,
2016), may be used as well, considering the fact that they
appear to be random for multiple purposes, such as Monte Carlo
simulations. Such sampling methods, e.g., Sobol or Halton sequences, have been shown, see, e.g., (de Dios Ortúzar et al.,
1994), to better cover the space than random sequences.
Remark 3.2. Assume that kf is a locally stabilizing control law
on XT , a sublevel set of Vf , which is positively invariant for the
system xk+1 = f (xk , kf (xk )) and Vf and L are, e.g., positive

definite quadratic functions. Considering that the sequence
U (k) provided by Algorithm 1 is a suboptimal solution solving
Problem 2.1, then, by the arguments in Section 2.2, Problem 2.1
is recursively feasible and it ensures stability of system (1) in
closed loop with uk = uw
k|k .
Remark 3.3. The working mechanism of Algorithm 1 resembles the working principle of the rollout algorithm described
in the survey paper (Bertsekas, 2005b). Therein, a rollout algorithm improves iteratively a base policy (here, the feasible
control sequence Uwarm (k)) to provide a suboptimal solution
to an optimal control problem via approximative DP and suboptimal control. The working principle of the rollout algorithm
with sampling based MPC, i.e., choosing at time k the iterated
solution of the previous time instance, k − 1, as a warm start,
and a sampling strategy, is not provided therein.
Remark 3.4. Due to the sampling procedure having a discontinuous behaviour, the inputs might be varying more than it is
safe for some applications. This problem might be alleviated
by either smoothening the inputs via interpolation of the obtained input sequence with the initial sequence, or by penalizing
∆uk = uk − uk−1 via constraints or via the cost function, such
that the input variability is limited to acceptable bounds.
3.2 Complexity analysis

In general, if we have p ∈ Z[2,n) processors, then the complexity of Algorithm 1 is


N (N + 1)
C = dn/pe c1
+ c2 N ,
(11)
2
where the term dn/pe appears due to the fact that a thread
can not engage in computations related to a subsequent time
horizon until all the threads have finalized the computations
related to the current time horizon j.
4. ILLUSTRATIVE EXAMPLES
The sampling–based suboptimal NMPC strategy proposed in
Section 3 is illustrated on three nonlinear systems, to highlight
various features of this method. All tests have been performed
on a system with the following specifications: Intel Core i73770 CPU 3.4GHz, 16GB RAM, 64-bit OS.
4.1 Cart–spring system

In order to analyze the complexity of Algorithm 1, the following assumptions are undertaken, for a given state x, input u and
input sequences U , U1 , U2 :
1) The cost of evaluating f (x, u) and performing a feasibility
test f (x, u) ∈ X is c1 ;
2) The cost of evaluating J(x, U ) is c2 ;
3) The cost of comparing J(x, U1 ) < J(x, U2 ) and changing
Jsub and Uwarm (k) if necessary, i.e., steps 11-13 in Algorithm 1, is negligible.
4) The current cost Jsub is instantaneously available for comparison with each of the nj samples according to step 11
in Algorithm 1 and each of the new sequences U (k)jq may
modify Jsub if the new cost Jnew is smaller than Jsub .
The complexity of Algorithm 1 for a given xk|k is the following:




N
−1
N
−1
X
X
C = c1 
(N − j)nj  + c2 
nj  .
(8)
j=0

N (N + 1)
+ c2 N.
(10)
2
The complexity bound given in (10), though quadratic in the
prediction horizon, yields a reasonable complexity, considering
that, in NMPC, a horizon N = 10 is considered a reasonably
large horizon.
C = c1

j=0

The method developed in this paper will first be applied to
a system incorporating an exponential nonlinearity, i.e., the
model of a cart with mass M , which is moving on a plane,
see (Raimondo et al., 2009). This cart is attached to a wall via
a spring with elastic constant k varying with the first state k =
k0 e−x1 , where x1 stands for the displacement of the carriage
from the equilibrium position. A damper acts as a resistor in
the system, with damping hd . The discretized nonlinear model
of the cart and spring system is the following:


x1 (k + 1)
x(k + 1) =
= f1 (x(k)) + F2 u(k),
(12)
x2 (k + 1)
where



x1 (k) + Ts x2 (k)
f1 (x(k)) =
,
ρ0 −x1
x2 (k) − Ts M
e
x1 (k) − Ts hMd x2 (k)
 T T
F2 = 0 Ms
,
(13)
where x2 is the velocity of the cart and u is an external force
which acts as an input to the system. The parameter values are
Ts = 0.4s, ρ0 = 0.33, M = 1, hd = 1.1.

If we assume nj ≤ n for all j ∈ Z[0,N −1] , then the complexity
(8) can be upper bounded by
N (N + 1)
C = nc1
+ c2 N n.
(9)
2
A possible reduction of the bound (9) might be attained, considering the fact that all the states subsequent to a non–feasible
state are no longer evaluated and checked for feasibility. This
means that, in step 9 of Algorithm 1, if φ(xk|k , U (k)jq , i) ∈
/X
for a specific i ∈ Z[j+1,N −1] , then φ(xk|k , U (k)jq , r) for all
r ∈ Z[i+1,N ] , are no longer evaluated, in which case steps 2)
and 3) are skipped all together.

The MPC controller has to steer the cart to the origin from
a non–zero initial state, while satisfying the input and state
constraints, which are
|u| ≤ 4.5, |x1 | ≤ 2.65,
(14)
and reducing the cost (2), where the stage cost and terminal
cost are quadratic functions, i.e., L(x, u) = xT Qx + uT Ru
and Vf (x) = xT P x. Choose the following parameters for the
MPC problem:


7.0814 3.3708
M = 4, P =
, Q = diag(1, 1), R = 1.
3.3708 4.2998

Consider now that many threads are available, from multiple
processors. Notice also that at each time in the horizon, all nj
computations can be performed separately. In these conditions,
assuming we have n threads available, then the complexity of
Algorithm 1 is upper bounded by

In (Raimondo et al., 2009) it is shown that the control law
u = kf (x) = − [ 0.8783 1.1204 ] f1 (x),
is locally stabilizing in the set
XT = {x|Vf (x) ≤ 4.7},

State trajectory

evaluation time[s]
0.16

6

terminal set
fmincon
nj=0

5

0.14
0.12

nj=5
nj=10

time [s]

4

nj=30

0.1
0.08

3
x

2

0.06

2

0.04

1

0.02

1

2

3

4

5

6

7

8

9

10

11

12

k

0

Fig. 4. Computation time.

−1
−2

−2

−1

0
x1

1

2

Fig. 2. State trajectory.
inputs
4
2

u

0
−2
−4
−6

2

4

6

8

10

12

14

16

18

20

k

Fig. 3. Inputs applied to the system.
which is a terminal set where the conditions for stability of
MPC mentioned in Section 2.2 are satisfied.

For the first experiment, fix N = 10. An initial Uwarm (0)
is provided by a random “oracle”. Consider nj = n, for all
j ∈ Z[0,N −1] , taking various values in the set {0, 5, 10, 30}.
For n = 0, Uwarm (0) is propagated though iterations without
any intervention or change from the sampling mechanism. This
serves as a reference, to notice the improvements brought in by
Algorithm 1. The results are illustrated in Fig. 2–5, where the
legend from Fig. 2 holds until Fig. 5. Iterations are considered
from k = 1 until k = 20. Though the different sampling
options proposed in Remark 3.1 provide in this example similar
outcomes, Halton points have been used here for illustration.
This choice is motivated by the practical feature that, adding
extra points only when they are required does not have impact
on the coverage of the set U.
In Fig. 2 and Fig. 3 the state trajectory and the inputs u(k)
applied to the system are illustrated. The constraint specifications (14) are satisfied for all cases, at all times. Notice the

Fig. 5. The cost Jsub .
Computational performance.
fmincon
sampling
bound (9)
bound (10)

20
15

Computational details for N<=20.

0.3

time [s]

25

time [s]

Algorithm 1 is applied for the MPC control of system (12). We
compare the results of this method with the results provided by
fmincon in Matlab, even though, as it will be seen later, the
optimization tool does not always provide the optimum, due
to local minima. The scalability of the algorithm is tested by
varying both the number of samples n and the control horizon
N . The tests illustrated in this paper are performed on a feasible
initial state x0|0 = [−2.5, 3]. The choice of the initial condition
does not influence greatly the results, which are similar for
other feasible initial states.

10

0.2

0.1

5
0

0
20

40

60

N

80

100

5

10

15

20

N

Fig. 6. Computational performance versus proposed bounds.
immediate smoothening of the trajectories and input sequence
even for a small n. In Fig. 4, the computational time, without
parallelization, is illustrated. At k = 1, the computational cost
of finding an oracle is included. Fig. 5 illustrates the values of
Jsub for each iteration. Notice, overall, that even for a small
number of samples, the performance of the closed loop system is significantly improved and the computational time is
promising, even for a non–parallel implementation. Also, as the
iteration k advances, the initially modest performance increases
significantly, due to the continuous improvement of Uwarm and
the receding horizon principle. Interestingly, even for n = 5,
at iteration k = 4, the cost Jsub (4) is smaller than the cost
computed via fmincon, which, due to local minima, provided
a feasible but not optimal solution.

Cost function
250
fmincon
nj=0

Jsub

200

nj=5

150

nj=20
nj=50

100
50
0

0

5

10

15

20

25
k

30

35

40

45

50

35

40

45

50

evaluation time[s]
1.4
1.2
1
t[s]

For the second experiment we aim to test the computational
complexity of Algorithm 1 in terms of the control horizon.
The implementation used here does not use parallelization. Fix
n = 10 and N takes values in the set {3, 10, 20, 50, 100}. The
results are illustrated in Fig. 6. With red, the computational
complexity for Algorithm 1 is depicted, and it is always smaller
than the bound (9), which is followed closely. It is expected that
on dedicated devices, the complexity of Algorithm 1 is smaller,
due to the processors not running in parallel threads related
to other system applications. With parallelization, further reduction in complexity is expected, as described in Section 3.2
and illustrated in Fig. 6. Notice that, for smaller horizons,
the complexity of Algorithm 1 is smaller than the complexity
of fmincon, even without parallelization. For N = 100,
fmincon provides solutions which are not feasible, while the
proposed Algorithm 1 still terminates in 21.6 seconds.

0.8
0.6
0.4

4.2 Buck–Boost power converter

0.2

Next, the bilinear model of a Buck–Boost power converter is
considered, as in (Spinu et al., 2011):


x(k)T C1
x(k + 1) = Ax(k) + Bu(k) +
u(k),
(15)
x(k)T C2
T

where x := [ vC iL ] ∈ X ⊂ R2 is the state vector consisting
of the voltage across the output capacitor and the current
T
through the filter inductor. The input u := [ d1 d2 ] ∈ U ⊂ R2
stands for the duty–cycle ratio of the control signal applied to
the switching node. The parameters are





− RH1 C 0
0 0
A = I2 + Ts
,
B
=
Ts ,
v
s
0
− RLL
L 0




0 0
0 − L1
C1 =
T
,
C
=
Ts ,
s
2
0 C1
0 0
with the values RL = 0.2Ω, C = 22µF , L = 220µH,
Ts = 10µs.
The aim of the control loop is to stabilize the system to the
T
T
equilibrium point xe := [ 20 0.5 ] , ue := [ 0.81 0.4 ] , under
2
the constraints iL ∈ R[0,3] , vC ∈ R[−0.1,22.5] , u ∈ R[0,1] . The
terminal controller



−0.0014 −0.3246
u = ue + K(x − xe ), K =
,
0.0001 −0.0055
stabilizing the system in the terminal set XT given in (Spinu
et al., 2011, set P in Section IV.C), and the quadratic cost with
the matrices


46.6617 42.8039
Q = diag(1, 2), R = diag(1, 1), P =
,
42.8039 69.4392
satisfies all the conditions for stability formulated in Section 2.2.
Similarly to the previous example, we apply Algorithm 1 for
the MPC control of system (15). Fix the initial state x0|0 =
T
[ 1 2 ] + xe , which is outside of the terminal set XT , and
N = 10. Uwarm (0) is given by an oracle and we use random
sampling of U. See in Fig. 7 the effect of varying n on the
cost function Jsub and the time necessary, per iteration, to
compute the corresponding control law. Notice the effect of
the unknown termination time on the evaluation time for the
optimization performed through fmincon and the relatively
equal computational time of Algorithm 1 over the iterations k.

0

0

5

10

15

20

25
k

30

Fig. 7. Performance of Algorithm 1 on the power converter.
All the constraints were satisfied for all the presented situations,
and it is expected that, through implementation on dedicated
multi–thread systems, the computational time for the control
law decreases to the extent of fitting the tight sampling period
of the power converter.
4.3 Wheeled mobile robot
The last example illustrates the methodology developed in this
paper for an obstacle avoidance task by a nonholonomic system
with trigonometric nonlinearities, due to kinematics, i.e., a
model of the wheeled mobile robot (WMR), as described in
(Kuhne et al., 2005):
"
#
x1 (k) + u1 (k) cos x3 (k)Ts
x(k + 1) = x2 (k) + u1 (k) sin x3 (k)Ts .
(16)
x3 (k) + u2 (k)Ts
In (16), the state x ∈ R3 describes the position and the orientation of the robot with respect to a global inertial frame
{O, X, Y }, and the input u ∈ R2 gives the linear and angular velocity, respectively. The parameter Ts = 0.1s is the
discretization period of system (16).
The MPC strategy aims at driving the WMR from an initial
T
state x0|0 = [ 0 6 0 ] to the origin of the inertial frame, i.e.,
T
xg = [ 0 0 0 ] , while satisfying the input constraints u1 ∈
R[−0.47,0.47] , u2 ∈ R[−3.77,3.77] . A quadratic cost function of
the form
N
−1
X
T
J(xk|k , U (k)) =xk+N |k P xk+N |k +
xTk+j|k Q(j)xk+j|k +
j=1

+

N
−1
X

uTk+j|k Ruk+j|k

j=0

is considered, with the parameters: Q(j) = 2j−1 Q, P =
50Q(N ), Q = diag(1, 1, 0.5), R = diag(0.1, 0.1), N = 5.
In this example we illustrate Algorithm 1 with the above parameters. Consider Uwarm (0) generated by an “oracle”, and
sampling of U based on random sequences. The result for
n = 30 is illustrated in Fig 8. fmincon could not be applied

State trajectory around obstacle.
6

Algorithm 1
SBMPC
obstacle

5

x2

4

improves in terms of cost an initially feasible control sequence.
This suboptimal NMPC strategy provides a promising computational complexity even for large control horizons, with good
perspectives for parallel implementation. Future work aims at
investigating the convergence of the proposed algorithm to the
optimal solution, the scalability in terms of system dimension
and the effect of parallelization on computational time for real–
life applications.

3

REFERENCES
2

1

0

-1.6

-1.4

-1.2

-1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

x1

Fig. 8. Trajectories of WMR.

total verification time[s]

0.8
0.7

time [s]

0.6
0.5
0.4
0.3
0.2
0.1
0
50

100

150

200

250

300

350

400

k

Fig. 9. Time performance of WMR.
due to feasibility issues related to non–convexity caused by the
presence of the obstacle. Notice the avoidance of the obstacle
of the WMR under the control law generated by Algorithm 1.
Due to the sampling mechanism, however, the inputs are not
smooth, which could be alleviated, e.g., by an interpolation
mechanism. For comparison we illustrate also the result using
SBMPC (Dunlap et al., 2010), with the same horizon N and
n = 30. Notice, in the state trajectory plot in Fig 8 the fact
that, after 400 iterations, the WMR did not reach xg yet, and
the input values are not yet 0, which means that the sampling
mechanism did not consider that the WMR is close to the goal.
In this case, it is recommended to sample more densely around
0 in the set U, rather than uniformly covering U with samples.
Notice in Fig 9 that, due to the building of a tree and the search
in a tree for the best path, SBMPC is more computationally
demanding. The complexity of Algorithm 1 with the given n
fits the sampling period Ts = 0.1s of the WMR, which makes
the method applicable for real–time control of this system.
5. CONCLUSION
In this paper, an algorithm based on sampling of the input space
at each time in the horizon was proposed, which iteratively

Bertsekas, D.P. (2005a). Dynamic programming and optimal
control, volume 1. Athena Scientific Belmont, MA.
Bertsekas, D.P. (2005b). Dynamic programming and suboptimal control: A survey from ADP to MPC. European Journal
of Control, 11(4), 310–334.
Chakrabarty, A., Dinh, V., Corless, M., Rundell, A., Zak,
S., and Buzzard, G. (2016).
Support vector machine informed explicit nonlinear model predictive control using low-discrepancy sequences. IEEE Transactions
on Automatic Control, preprint published on–line, DOI:
10.1109/TAC.2016.2539222, PP(99).
de Dios Ortúzar, J., Willumsen, L.G., et al. (1994). Modelling
transport. Wiley New Jersey.
Dunlap, D.D., Caldwell, C.V., and Collins, E.G. (2010). Nonlinear model predictive control using sampling and goaldirected optimization. In 2010 IEEE International Conference on Control Applications, 1349–1356. IEEE.
Grüne, L. and Pannek, J. (2011). Nonlinear model predictive
control. Theory and algorithms. London: Springer–Verlag.
Johansen, T.A. (2004). Approximate explicit receding horizon
control of constrained nonlinear systems. Automatica, 40(2),
293–300.
Kuhne, F., Lages, W.F., and Da Silva, J.G. (2005). Point stabilization of mobile robots with nonlinear model predictive
control. In IEEE International Conference Mechatronics and
Automation, 2005, volume 3, 1163–1168. IEEE.
LaValle, S.M. (1998). Rapidly-exploring random trees: A new
tool for path planning. TR 98-11, Computer Science Dept.,
Iowa State University.
Ławryńczuk, M. (2009). Computationally efficient nonlinear
predictive control based on state-space neural models. In International Conference on Parallel Processing and Applied
Mathematics, 350–359. Springer.
Lazar, M. and Heemels, W. (2009). Predictive control of hybrid systems: Input-to-state stability results for sub-optimal
solutions. Automatica, 45(1), 180–185.
Lee, J.M. and Lee, J.H. (2004). Approximate dynamic programming strategies and their applicability for process control: A review and future directions. International Journal of
Control Automation and Systems, 2, 263–278.
Likhachev, M. and Stentz, A. (2008). R* search. Lab Papers
(GRASP), 23.
Magni, L., Raimondo, D.M., and Allgöwer, F. (2009). Nonlinear model predictive control: Towards new challenging
applications, volume 384. Springer, Lecture notes in Control
and Information Sciences.
Mayne, D. and Rawlings, J. (2009). Model predictive control:
theory and design. Madison, WI: Nob Hill Publishing, LCC.
Muñoz de la Peña, D., Alamo, T., Bemporad, A., and Camacho,
E.F. (2004). A dynamic programming approach for determining the explicit solution of linear MPC controllers. In
IEEE Conference on Decision and Control.

Pannocchia, G., Rawlings, J.B., and Wright, S.J. (2011). Conditions under which suboptimal nonlinear MPC is inherently
robust. Systems & Control Letters, 60(9), 747–755.
Piovesan, J.L. and Tanner, H.G. (2009). Randomized model
predictive control for robot navigation. In IEEE International Conference on Robotics and Automation, 2009, 94–99.
Kobe, Japan.
Raimondo, D.M., Limon, D., Lazar, M., Magni, L., and Camacho, E.F. (2009). Min-max model predictive control of nonlinear systems: A unifying overview on stability. European
Journal of Control, 15(1), 5–21.
Scokaert, P.O., Mayne, D.Q., and Rawlings, J.B. (1999). Suboptimal model predictive control (feasibility implies stability). IEEE Transactions on Automatic Control, 44(3), 648–
654.
Spinu, V., Lazar, M., and van den Bosch, P. (2011). An explicit
state-feedback solution to constrained stabilization of dc-dc
power converters. In 2011 IEEE International Conference
on Control Applications (CCA), 1112–1118. IEEE.

