On the Approximation of Constrained Linear Quadratic
Regulator Problems and their Application to Model
Predictive Control - Supplementary Notes
Michael Muehlebach and Raffaello D’Andrea∗

arXiv:1803.05510v1 [cs.SY] 23 Feb 2018

Abstract
By parametrizing input and state trajectories with basis functions different approximations to the
constrained linear quadratic regulator problem are obtained. These notes present and discuss technical
results that are intended to supplement a corresponding journal article. The results can be applied in a
model predictive control context.

1

Outline

Sec. 2 summarizes the notation. The resulting approximations are given by convex finite-dimensional optimization problems that have a quadratic cost, linear equality constraints, and linear semi-infinite inequality
constraints. Sec. 3 suggests several strategies for dealing with these semi-infinite constraints. One of these
strategies is discussed in more details in Sec. 4, as it is found to be computationally efficient in a model
predictive control context. The remaining sections discuss several technical results. Sec. 5 and Sec. 6 are
related to an upper, respectively lower bound approximation of the cost of the underlying constrained linear
quadratic regulator problem. Sec. 7 presents a result that is related to the convergence of the different approximations. Sec. 8 shows recursive feasibility and closed-loop stability of the resulting model predictive
control algorithm. Sec. 9 and Sec. 10 highlight important properties of the basis functions parametrizing
input and state trajectories.

2

Notation

We are concerned with the approximation of the following optimal control problem

1
||x||2 + ||u||2 + ψ(xT )
2
s.t. ẋ = Ax + Bu, x(0) = x0 , x(T ) = xT ,
Cx x + Cu u ≤ b, xT ∈ X ,
x ∈ L2n , u ∈ L2m ,

J∞ := min

(1)

where the space of square integrable functions mapping from the interval I := (0, T ) to Rq is denoted by
L2q , where q is a positive integer; and the L2q -norm is defined as the map L2q → R,
Z
2
x → ||x||, ||x|| := xT x dt,
(2)
I
∗

Michael Muehlebach and Raffaello D’Andrea are with the Institute for Dynamic Systems and Control, ETH Zurich. The
contact author is Michael Muehlebach, michaemu@ethz.ch. This work was supported by ETH-Grant ETH-48 15-1.

1

with dt the Lebesgue measure. The function ψ : Rn → R is assumed to be positive definite and strongly
convex, A ∈ Rn×n , B ∈ Rn×m , Cx ∈ Rnc ×n , Cu ∈ Rnc ×m , and b ∈ Rnc are constant, and the set X is closed
and convex. The dynamics as well as the stage constraints are assumed to be fulfilled almost everywhere.
Thus we simply write
f = g, f ≤ g,
(3)
when we mean f (t) = g(t), respectively f (t) ≤ g(t) for all t ∈ I almost everywhere, with f, g ∈ L2q .
The weak derivative of x is denoted by ẋ.1 To simplify notation we abbreviate the domain of the objective
function by
X := L2n × L2m × Rn .
(4)
We assume throughout the article that the constraints in (1) are nonempty, i.e. there exists trajectories x
and u, fulfilling the dynamics, the initial condition, the constraints, and thus achieve a finite cost.
The main motivation for studying problem (1) comes from the fact that (1) often serves as a starting
point for model predictive control (MPC).
The problem (1) can be written in compact form as
J∞ = min ||x||2 + ||u||2 + ψ(xT )
s.t. (x, u, xT ) ∈ C ∩ D,

(5)

where the set C contains all trajectories x ∈ L2n , u ∈ L2m that fulfill the constraints and D the trajectories
that fulfill the dynamics and the initial condition, that is,
C := {(x, u, xT ) ∈ X | Cx x + Cu u ≤ b, xT ∈ X },
D := {(x, u, xT ) ∈ X|ẋ = Ax + Bu, x(0) = x0 , x(T ) = xT }.

(6)
(7)

It was shown that (5) can be approximated by two auxiliary optimization problems, where input and state
trajectories are parameterized with basis functions,
x̃(t) = (In ⊗ τ s (t))T ηx ,

ũ(t) = (Im ⊗ τ s (t))T ηu .

The parameter vectors are denoted by ηx ∈ Rns and ηu ∈ Rms , whereas τ s (t) := (τ1 (t), τ2 (t), . . . , τs (t)) ∈
Rs contains the first s basis functions, ⊗ denotes the Kronecker product, and Iq ∈ Rq×q refers to the
identity matrix for any integer q > 0. The superscript s refers to the number of basis functions used for
the approximation. For ease of notation the superscript s will be dropped whenever the number of basis
functions is clear from context. The basis functions are required to satisfy the following assumptions:
A1) The basis functions τi ∈ L21 , i = 1, 2, . . . , s are linearly independent and orthonormal with respect to
the standard L21 -scalar product.
A2) The basis functions fulfill τ̇ s (t) = Ms τ s (t) for all t ∈ I, for some Ms ∈ Rs×s .
The finite-dimensional subspace spanned by the first s basis functions will be denoted by X s ,
X s := {(x, u, xT ) ∈ X | ηx ∈ Rns , ηu ∈ Rms , x = (In ⊗ τ s )T ηx , u = (Im ⊗ τ s )T ηu }.

(8)

We can think of an element in X s not only as an element in X (i.e. a tuple of a finite-dimensional vector
and two square integrable functions), but also as a finite-dimensional vector given by the corresponding
parameter vectors ηx and ηu . To make this distinction explicit, we introduce the map π qs : L2q → Rqs ,
defined as
Z
x → (Iq ⊗ τ s )xdt,
(9)
I

The equations of motion imply that ẋ ∈ L2n , which can be used to conclude that x has a unique absolutely continuous
representative defined on the closure of I (a classical solution of the equations of motion). With x(t) we refer to the value this
unique absolutely continuous representative takes at time t ∈ [0, T ].
1

2

which maps an arbitrary element x ∈ L2n to its first s basis function coefficients. Similarly, we define
π s : X → Rns × Rms × Rn as
(x, u, xT ) → (π ns (x), π ms (u), xT ).
(10)
As a consequence, we write π s (x) for describing the finite dimensional representation of x ∈ X s , that is,
its representation in terms of the parameter vectors ηx and ηu . The adjoint map (π qs )∗ : Rqs → L2q is given
by
η → (Iq ⊗ τ )T η,
(11)
and is used to obtain the trajectory corresponding to the vector η ∈ Rqs , containing the first s basis function
coefficients. Similarly, we define (π s )∗ : Rns × Rms × Rn → X as
(ηx , ηu , xT ) → ((π ns )∗ (ηx ), (π ms )∗ (ηu ), xT ).

(12)

The composition (π s )∗ π s : X → X yields the projection of an element x ∈ X onto the subspace X s ⊂ X.
The sets C and D are approximated in two different ways,
CUs := C ∩ X s ,

(13)

CLs := {(x, u, xT ) ∈ X |

Z

δ p̃T (−Cx x − Cu u + b)dt ≥ 0

I

∀δ p̃ = (Inc ⊗ τ )T δηp : δ p̃ ≥ 0, δηp ∈ Rnc s ; xT ∈ X },
DUs
DLs

s

:= D ∩ X ,

:= {(x, u, xT ) ∈ X | δ p̃ = (In ⊗ τ )T δηp ,
Z
Z
T
˙
− δ p̃ xdt − δ p̃T (Ax + Bu)dt − δ p̃(0)T x0 + δ p̃(T )T xT = 0, ∀δηp ∈ Rns },
I

(14)
(15)

(16)

I

resulting in the following two approximations to (5)
Js := inf ||x̃||2 + ||ũ||2 + ψ(xT )
s.t. (x̃, ũ, xT ) ∈ CUs ∩ DUs ,

(17)

J˜s := min ||x̃||2 + ||ũ||2 + ψ(xT )
s.t. (x̃, ũ, xT ) ∈ CLs ∩ DLs .

(18)

and

By definition of the constraints CUs and DUs , the minimizer of (17) (for s ≥ s0 ) is required to be an element
of X s . Consequently, the problem (17) is equivalent to
Js = inf |ηx |2 + |ηu |2 + ψ(xT )
s.t. (ηx , ηu , xT ) ∈ π s (CUs ) ∩ π s (DUs ),

(19)

which corresponds to a convex finite-dimensional optimization problem. Likewise, it is found that the
convex finite-dimensional problem
J˜s = min |ηx |2 + |ηu |2 + ψ(xT )
s.t. (ηx , ηu , xT ) ∈ π s (CLs ) ∩ π s (DLs ),

(20)

is equivalent to (18) in the sense that its (unique) minimizer (ηx , ηu , x̄T ) is related to the minimizer (x, u, xT ) ∈
X of (18) by x = (π ns )∗ (ηx ), u = (π ms )∗ (ηu ), x̄T = xT , and achieves the same cost.
The proposed approximations can be applied in the context of MPC. We will show that by repeatedly
solving the infinite-horizon optimal control problem (19) (with I = (0, ∞), ψ = 0, xT = limt→∞ x(t) = 0),
3

recursive feasibility and closed-loop stability are inherent to the resulting MPC algorithm. For the sake of
completeness, (19) is written out as
min ηxT (Q ⊗ Is )ηx + ηuT (R ⊗ Is )ηu
s.t. (In ⊗

MsT

(21)

− A ⊗ Is )ηx − (B ⊗ Is )ηu = 0,

(In ⊗ τ (0))T ηx = x0 ,
(Cx ⊗ τ (t))T ηx + (Cu ⊗ τ (t))T ηu ≤ b, ∀t ∈ [0, ∞),
where the input and state costs are weighted with the matrices Q > 0 and R > 0, which is common in
MPC.
In the following an efficient numerical solution algorithm for (21) will be discussed.

3

Implementation of the semi-infinite constraint

The optimization (21) is a convex finite-dimensional optimization problem. However, it is not a quadratic
program, as it includes the semi-infinite constraint
(Cx ⊗ τ (t))T ηx + (Cu ⊗ τ (t))T ηu ≤ b, ∀t ∈ [0, ∞).

(22)

As a result, (21) cannot be solved by a standard quadratic programming solver. Three different approaches
to deal with the semi-infinite constraint are immediate:
1) global polyhedral approximation
2) sum-of-squares approximation
3) local polyhedral approximation (active-set approach)
The first is based on a fixed polyhedral approximation, leading to a quadratic program. The second is based
on exploiting polynomial basis functions for obtaining a characterization using linear matrix inequalities,
whereas the third is based on an iterative constraint sampling scheme, resulting in a local polyhedral approximation. In the following subsections we will focus on approach 1 and 2. We will describe an efficient
solution algorithm based on approach 3 in detail in Sec. 4.
It turns out that it is enough to check the constraint (22) over a compact time interval, instead of the
unbounded interval t ∈ [0, ∞). This is because the basis functions are assumed to be linearly independent
and exponentially decaying according to Assumption A1 and A2. A formal proof of this claim can be found
in Sec. 9. The compact time interval for which the constraint (22) has to be checked is denoted by [0, Tc ]
and depends on the choice of basis functions, on the order s, and in some cases also on the bound b (see
Sec. 9).

3.1

Global polyhedral approximation

In order to construct a global polyhedral approximation of the set CUs we assume that an upper bound on the
achievable cost Js , denoted by J¯s , is available. In order to simplify the discussion, we assume further that
for now Cx and Cu are row vectors, and that b is a scalar. The proposed scheme can be readily extended
to the case where Cx and Cu are matrices, and b is a vector. The approximation is based on constraint
sampling, where we tighten the constraint slightly to
Cx x̃(ti ) + Cu ũ(ti ) ≤ (1 − )b,

(23)

with  > 0, and where ti , denotes the constraint sampling instances, i = 1, 2, . . . , ni , which are to be
determined. The algorithm is based on the following two steps:
4

1) Compute
h(t) := max Cx x̃(t) + Cu ũ(t) − b
s.t. Cx x̃(ti ) + Cu ũ(ti ) ≤ (1 − )b, i = 1, . . . , ni ,
(x̃, ũ, lim x̃(te )) ∈ DUs ,
te →∞

||x̃|| + ||ũ||2 ≤ J¯s ,
2

for all times t ∈ Is , where Is contains a finite number of sampling instances (to be made precise
below).
2) Find the local peaks of h(t), denoted by t∗i . Add each t∗i to the constraint sampling points if h(t∗i ) >
−b/2. Repeat the procedure until h(t) ≤ −b/2 for all t ∈ Is .
Note that the function h(t) is again only evaluated at the discrete time points t ∈ Is . The index set Is has
to be chosen such that h(t) ≤ −b/2 for all t ∈ Is implies that h(t) ≤ 0 for all t ∈ [0, ∞). As remarked
earlier, due to the fact that the basis functions are exponentially decaying and linearly independent it is
enough to check h(t) ≤ 0 for all t ∈ [0, Tc ], for a fixed time Tc , as h(t) ≤ 0 for all t ∈ (Tc , ∞) will be
fulfilled automatically. Moreover, a Lipschitz constant of
Cx x̃(t) + Cu ũ(t) − b

(24)

can be found by using an upper bound on its time-derivative, that is, for example,
|τ (t)T M T (Cx ηx + Cu ηu )| ≤ |τ (t)||M T (Cx ηx + Cu ηu )|,

(25)

where the first term can be bounded for all t ∈ [0, ∞) due to the fact that the basis functions are exponentially decaying and the second term can be bounded using the fact that the cost Js is below J¯s . We therefore
choose the index set Is as
2L
Is = {tk < Tc | tk = k , k = 0, 1, 2, . . . },
(26)
b
where L denotes a Lipschitz constant of (24).
It is important to note that the optimization in step 1 imposes the dynamics and the upper bound J¯s on
the cost. Both constraints tend to reduce the number of constraint sampling points ti greatly. The initial
condition x0 enters the optimization as an optimization variable. The optimization in step 1 represents a
quadratically constrained linear program for each time instant t ∈ Is , and as such, it can be solved using
standard software packages. The whole procedure for determining the constraint sampling points is done
offline. Once these time instances are found, the optimization problem that is solved online reduces to a
quadratic program. The number of constraint sampling points ti is upper bounded by the cardinality of the
index set Is , and thus guaranteed to be finite. Due to the fact that the above procedure is greedy, it will not
necessarily lead to the smallest number of constraint sampling points.

3.2

Sum-of-squares approximation

In case exponentially decaying polynomials are used as basis functions, sum-of-squares techniques can be
applied. In particular, it is shown in [3] that the set
{η ∈ Rs | η T (1, t, . . . , ts−1 ) ≥ 0, ∀t ∈ [0, ∞)}

(27)

can be expressed using matrix inequalities that are linear in the coefficients η. In the case of exponentially
decaying polynomials it is therefore enough to approximate the exponential decay by a polynomial upper
bound (for example by appropriately truncating a Taylor series expansion at 0), in order to approximate
the constraint (22) in a slightly conservative manner. As a result, by applying the results from [3], the
optimization problem (21) is approximated by a semidefinite program that can be solved using standard
optimization routines.
5

4

A dedicated active-set approach

In the following section we present an efficient optimization routine for solving (19) (and likewise (21)).
The method is an extension of traditional active set methods and generalizes to optimization problems with
a linear quadratic cost function, linear equality constraints, and linear semi-infinite inequality constraints,
i.e.
ˆ c ) := min z T Hz
J(I
s.t. Aeq z = beq ,

(28)
(29)

lb ≤ (Inc ⊗ τ (t))T Cz z ≤ lu ,

∀t ∈ Ic ,

(30)

where Ic is any subset of the non-negative real line. In case of (21), the interval Ic is taken to be [0, ∞).
Note that an optimization problem, whose objective function has a linear part, can be brought to the form
(28) by completing the squares. It is assumed that the optimization problem (28) is feasible, that lb < 0
and lu > 0, and that the Hessian H is positive definite, which guarantees existence and uniqueness of the
corresponding minimizer2 .
The method is based on the observation that if the set Ic consists merely of a collection of time instants
(constraint sampling instances) ti , (28) reduces to a quadratic program that can be solved efficiently. Moreover, due to the fact that the basis functions fulfill Assumptions A1 and A2, a trajectory parametrized with
the basis functions has a finite number of maxima and minima, as is shown in Sec. 10. Consequently, (28)
has only a finite number of active constraints. The collection of the time instants corresponding to these
active constraints will be denoted by Ic∗ . If this finite collection of constraint sampling instants is known
ˆ c∗ ) = J(I
ˆ c ). In
ahead of time, one could simply solve (28) with respect to Ic∗ instead of Ic , resulting in J(I
k+1
k
k
k
ˆ
ˆ
addition, for any subset Ic of Ic it holds that J(Ic ) ≥ J(Ic ), and likewise, if Ic is a subset of Ic we have
ˆ c ) can be
ˆ k ). Hence, a monotonically increasing sequence J(I k ), bounded above by J(I
ˆ k+1 ) ≥ J(I
that J(I
c
c
c
k+1
k
k
constructed using any sequence of sets Ic that fulfill Ic ⊂ Ic ⊂ · · · ⊂ Ic . In particular, such sets are obtained by starting with an arbitrary initial guess Ic0 containing a finite number of constraint sampling points
(or even the empty set), and by adding at least one constraint violation point at each iteration. Moreover,
at each iteration the inactive constraints contained in the set Ick can be removed, as this will not alter the
ˆ k ). In that way, the number of constraint sampling instances contained
optimizer nor the optimal value J(I
c
k
in Ic remains finite. This motivates Alg. 1, which solves (28) up to a given tolerance, by constructing an
approximation to the set of active constraints Ic∗ .
ˆ k ) = J(I
ˆ ∗ ) = J(I
ˆ c ). In order to achieve constraint
Proposition 4.1 Alg. 1 converges, that is, limk→∞ J(I
c
c
violations smaller than  at most
ˆ c ) − J(I
ˆ 0 ))
4cτ (J(I
c
(31)
2
σ
steps are required, where cτ is defined as
cτ := sup |τ (t)|2 ,

(32)

t∈Ic

and σ denotes the smallest eigenvalue of the Hessian H.
ˆ k ) is monotonically increasing whenever
Proof From the above arguments it can be concluded that J(I
c
ˆ k ) < J(I
ˆ c ) and upper bounded by J(I
ˆ c ). Therefore the sequence J(I
ˆ k ) converges. It remains to show
J(I
c
c
ˆ ck ) converges to J(I
ˆ c ). The strong convexity of the objective function can be used to establish
that J(I
ˆ ck+1 ) − J(I
ˆ ck )),
|z k+1 − z k |2 ≤ 4σ −1 (J(I
2

(33)

This is due to the fact that the constraints describe a closed convex set and due to the strong convexity of the objective
function.

6

Algorithm 1 Iterative constraint sampling
Initialize: initial guess for the constraint sampling points: Ic0 = {t1 , t2 , . . . , tN }; maximum number of
iterations: MAXITER; constraint satisfaction tolerance: ;
1: k = 0
2: for k <MAXITER do
ˆ k)
3:
solve (28) for Ick → z k , J(I
c
4:
if infeasible then
5:
abort
6:
else if z k fulfills (30) for all t ∈ I (with tol. ) then
7:
algorithm converged
8:
return z k
9:
end if
10:
find at least one constraint violation instant → tc
11:
remove inactive time instants in Ick
12:
Ick+1 = Ick ∪ {tc }, k = k + 1
13: end for
ˆ k ) and J(I
ˆ k+1 ). As a result we can conclude
where z k+1 and z k are the minimizer corresponding to J(I
c
c
k
k
that z converges, and that limk→∞ z is well-defined and satisfies the constraint (30) (by a contradiction
ˆ c ), which,
ˆ k ) ≥ J(I
argument). It is therefore a feasible candidate for (28), implying that limk→∞ J(I
c
ˆ ck ) ≤ J(I
ˆ c ) for all integers k, leads to limk→∞ J(I
ˆ ck ) = J(I
ˆ c ).
combined with J(I
It remains to show that (31) is fulfilled. To that extent, let  > 0 denote the smallest constraint violation
that occurs within the first N steps. For all k ≤ N − 1, there exists the constraint violation point tkj , which
will be added to Ick , and therefore
 ≤ |(e(tkj ) ⊗ τ (tkj ))T (z k − z k+1 )|,

(34)

where e(tkj ) is a canonical unit vector. Combining the Cauchy-Schwarz inequality and the above bound on
|z k+1 − z k | results in
ˆ ck+1 ) − J(I
ˆ k )).
2 ≤ 4cτ σ −1 (J(I
(35)
c
By summing over the first N steps we arrive at
ˆ cN ) − J(I
ˆ c0 )) ≤ 4σ −1 cτ (J(I
ˆ c∗ ) − J(I
ˆ c0 )),
N 2 ≤ 4σ −1 cτ (J(I

(36)

ˆ ck ) is strictly increasing and bounded above by J(I
ˆ c∗ ). Dividing by 2 on both sides
since the sequence J(I
concludes the proof.

4.1

Implementation details

Alg. 1 can be naturally embedded in an active-set method. An introduction to active-set methods for solving
quadratic programs can be found for example in [1, Ch. 10]. Starting from an initial guess of the active conˆ c0 ) is solved:
straint sampling instants, which is denoted by Ic0 , the quadratic program with optimal cost J(I
0
This is done by initially assuming that all constraints in the set Ic are active. The resulting optimization
problem reduces to an equality constrained quadratic program, whose solution can be calculated by solving
a linear system of equations. The Lagrange multipliers corresponding to (30), which are denoted by µ(t),
t ∈ Ic0 , are evaluated subsequently. If all the constraints are indeed active, then the optimizer to the quadratic
ˆ c0 ) has been found. If, however, not all the constraints are found to be active, that is,
program with cost J(I
if there are some Lagrange multipliers that are zero, the standard active-set procedure, see [1, Ch. 10] is
used to find the subset of active constraints Ia ⊂ Ic0 . Provided that the active set Ia ⊂ Ic0 and the optimizer
ˆ 0 ) has been found, the constraint (30) is then checked for all t ∈ Ic .
to the quadratic program with cost J(I
c
7

If no constraint violations occur, the solution to (28) has been found. If constraint violations occur, the time
instant tc for which a violation occurs, is added to the set of active constraints Ia resulting in Ic1 = Ia ∪ {tc }.
The above procedure is then repeated until convergence.
Each iteration requires solving equality constrained quadratic programs of the type
min z T Hz

(37)

s.t. Aeq z = beq ,

(c(t) ⊗ τ (t))T Cz z = la (t),

∀t ∈ Ia ,

(38)

where c(t) ∈ Rnc , la (t) ∈ R, t ∈ Ia , and Ia ⊂ Ick describe the active constraints corresponding to (30). Due
to the fact that very few constraints are expected to be active, we use a range space approach, [1, p. 238].
To that extent, the equality constraint is eliminated and the optimizer z ∗ corresponding to (37) is rewritten
as
X
z ∗ = b̂ + ĤCzT
(c(t) ⊗ τ (t))T µ(t),
(39)
t∈Ia

where the dual variable µ(t) ∈ R, defined for t ∈ Ick , satisfies
X
(c(t) ⊗ τ (t))µ(t) = la (tj ) − (c(tj ) ⊗ τ (tj ))T Cz b̂,
(c(tj ) ⊗ τ (tj ))T Cz ĤCzT

(40)

t∈Ia

for all tj ∈ Ia ⊂ Ick , and µ(t) = 0 for all t ∈ Ick \ Ia , with
Ĥ := H −1 ATeq (Aeq H −1 ATeq )−1 Aeq H −1 − H −1 ,
b̂ := H −1 ATeq (Aeq H −1 ATeq )−1 beq .

(41)
(42)

The dual variable µ(t) is therefore obtained by solving (40), and the optimizer z ∗ is then determined via
(39). At each iteration, a single constraint is either added or removed. Therefore, in order to efficiently find
a solution to (40), a LDLT -decomposition of the matrix
{(c(tj ) ⊗ τ (tj ))T Cz ĤCzT (c(ti ) ⊗ τ (ti ))}(tj ,ti )∈Ia ×Ia

(43)

is computed and adapted at each step by performing rank-1 updates. The matrix Ĥ and the vector b̂ are
precomputed. For additional details regarding the regularity of (43), and issues related to cycling and
stalling we refer to [1, Ch. 10] and [4, p. 467].

4.2

Constraint check

It remains to explain how to efficiently check whether the constraint (30) is fulfilled for a given solution
candidate z. We assume that the interval Ic has the form Ic = [0, Tc ]. As it has been explained in Sec. 3,
the constraint check over the interval [0, ∞) reduces to the check over a compact interval, provided that the
basis functions fulfill Assumptions A1 and A2.
A straightforward approach would be to exploit the specific structure of the basis functions. For example, if the basis functions consist of exponentially decaying polynomials having a degree of at most 4,
determining the stationary points of (30) amounts to solving a quartic equation, which can be done analytically. As a result, it would be enough to check the constraints at these stationary points in order to determine
if the constraint is satisfied or not.
We propose a more general approach that is based on local Taylor approximations, and thus valid for
arbitrary basis functions compatible with Assumptions A1 and A2. In order to simplify the discussion, we
consider the special case of (30), where Cz is the identity and nc = 1. The resulting algorithm extends
naturally to the more general case. According to Taylor’s theorem we obtain the following identity
τ (t)T z = τ (0)T z + τ̇ (0)T zt + τ̈ (0)T z
8

t2
t3
+ τ (3) (t̄)T z ,
2
6

(44)

where t̄ ∈ [0, t]. As will become clear in the following, a third-order Taylor expansion represents a good
compromise between approximation quality and computational effort. The last term of the previous equation can be bounded by the Cauchy-Schwarz inequality leading to
|τ (3) (s)T z| ≤ sup |τ (t̄)||(M T )3 z| =: R(z),

∀s ∈ [0, t].

(45)

t̄∈Ic

As a consequence, the following upper and lower bounds are obtained
bl (t) ≤ τ (t)T z ≤ bu (t),

(46)

for all t ∈ Ic , with
t2
t3
− R(z) ,
2
6
2
t3
t
bu (t) := τ (0)T z + τ̇ (0)T zt + τ̈ (0)T z + R(z) .
2
6
bl (t) := τ (0)T z + τ̇ (0)T zt + τ̈ (0)T z

(47)
(48)

The situation is exemplarily depicted in Fig. 1 (right). Given that τ̇ (0)T z ≥ 0 the lower bound attains its
maximum at time
p
τ̈ (0)T z + (τ̈ (0)T z)2 + 2R(z)τ̇ (0)T z
tu :=
> 0,
(49)
R(z)
whereas if τ̇ (0)T z < 0 the upper bound attains its minimum at time
p
−τ̈ (0)T z + (τ̈ (0)T z)2 − 2R(z)τ̇ (0)T z
tl :=
> 0.
R(z)

(50)

Thus, if the lower bound exceeds lu or the upper bound drops below ll , that is, if
bl (tu ) > lu ,

or

bu (tl ) < ll ,

(51)

the constraint (30) is guaranteed to be violated at time tu , respectively at time tl . If this is not the case, we
are guaranteed that the constraint (30) is satisfied in the interval [0, ts ], with
ts := min{t1 , t2 | bu (t1 ) = lu , bl (t2 ) = ll }.

(52)

Thus, finding the value ts requires the solution of two cubic equations, which stems from the fact that a
third order Taylor approximation was used as a starting point. By shifting the parameter vector in time
by ts and repeating the above procedure, the constraint is either found to be satisfied for all t ∈ Ic or
a constraint violation is detected. Shifting the parameter vector by ts amounts in multiplying z with the
matrix exponential
T
eM ts z → z.
(53)
The procedure is illustrated by the flow chart depicted in Fig. 1 (left).
The efficiency of the proposed strategy can be improved via the following observation. A constraint
violation occurring close to t = 0, is often found within the first few iterations. However, if no constraint
violation occurs, the whole interval Ic needs to be traversed, which tends to increase computation. The
computational effort may be reduced by including additional conservative constraint satisfaction checks.
For example, upper and lower bounds can be tightened by a factor γ ∈ (0, 1) such that the satisfaction of
γll ≤ τ (ti )T z ≤ γlu

(54)

for certain time instances ti implies (30) (for all t ∈ Ic ). As a result, at each iteration, the above inequality
is checked. If it is found to be fulfilled, then constraint satisfaction can be guaranteed.
9

Ic = [0, Tc ], z
0→t

constraint
satisfaction

true

check (54)

ts

lu
compute tu , tl

bu (t)
τ (t)T z

tu

t
constraint
violation

true

check (51)

ll
compute ts

constraint
satisfaction

true

bl (t)
T

eM ts z → z
t + ts → t

t + ts > Tc

Figure 1: Left: Flow chart that illustrates the proposed constraint satisfaction check. Right: Illustration of
the upper and lower bounds bu (t) and bl (t) obtained from the Taylor expansion of τ (t)T z. In that case the
test (51) is indecisive, and constraint satisfaction over the interval [0, ts ] can be guaranteed.

5

Properties B1-B5

The sets C, CUs , and CLs have the following properties
B1) the sets CUs , CLs , and C are closed and convex.
B2) the sets π s (CUs ), π s (CLs ) are closed and convex.
B3) CUs ⊂ CUs+1 ⊂ C.
B4) CLs ⊃ CLs+1 ⊃ C.
B5) (π s )∗ π s (CUs ) ⊂ CUs , (π s )∗ π s (CLs ) ⊂ CLs .
We will sketch the proofs of Properties B1-B5 in the following. It follows from the linearity of the stage
constraint, and the fact that the terminal constraint X is convex, that the sets C, CUs , and CLs are likewise
convex. We will sketch the proof that the set CLs is closed. The argument can be translated to the set C by
using the following variational formulation
Z
C = {(x, u, xT ) ∈ X |
δpT (−Cx x − Cu u + b)dt ≥ 0
I
Z
2
∀δp ∈ Lnc : δp ≥ 0,
δpT bdt < ∞; xT ∈ X }. (55)
I

CLs

We will argue indirectly, i.e. that the complement of
is open. To that extent we choose (x, u, xT ) ∈
s
X \ CL . As a result, there exists a test function δ p̃, with δ p̃ ≥ 0, which is spanned by the first s basis
functions and is such that
Z
δ p̃T (−Cx x − Cu u + b)dt < 0.
(56)
I

10

For any x̂ ∈ L2n , û ∈ L2m , and x̂T ∈ Rn , with ||x̂ − x|| < , ||û − u|| < , and |x̂T − xT | < , it follows that
Z
Z
Z
T
T
δ p̃ (−Cx x̂ − Cu û + b)dt = δ p̃ (−Cx x − Cu u + b)dt + δ p̃T (Cx (x − x̂) + Cu (u − û))dt, (57)
I

I

I

where the last integral can be bounded by (using the Cauchy-Schwarz inequality)
||δ p̃||(|Cx |||x − x̂|| + |Cu |||u − û||) < ||δ p̃||(|Cx | + |Cu |).

(58)

We can infer from X being closed, that there exists an open ball centered at xT , which does not intersect
X . As a result, by choosing  small enough, it can be concluded that
Z
δ p̃T (−Cx x̂ − Cu û + b)dt < 0, x̂T 6∈ X ,
(59)
I

and therefore (x̂, û, xT ) ∈ X \ CLs , for all x̂ ∈ L2n , û ∈ L2m , x̂T ∈ Rn with ||x − x̂|| < , ||u − û|| < , and
|x̂T − xT | < . Hence, the complement of CLs is open, and therefore CLs is closed.
Note that the set CUs is given by the intersection of the set C with the linear subspace X s spanned by the
first s basis functions. Both of these sets are closed3 implying that CUs is closed as well.
The projection π s is linear, which asserts the convexity of the sets π s (CUs ) and π s (CLs ). Moreover, it is
surjective, and hence, by the open mapping theorem, it follows directly from CLs and CUs being closed that
π s (CLs ) and π s (CUs ) are closed as well.
The inclusion CUs ⊂ CUs+1 follows directly from CUs = C ∩ X s and the inclusion X s ⊂ X s+1 . In
other words, given (x, u, xT ) ∈ CUs , the parameter vectors ηx and ηu , corresponding to the state and input
trajectories x and u, can be extended with zeros resulting in trajectories x̂, û spanned by s+1 basis functions.
But x̂ = x and û = u and therefore (x, u, xT ) ∈ CUs+1 . The inclusion CLs ⊃ CLs+1 follows from the fact that
the dimension of the subspace to which the test functions δ p̃ are constrained increases with s.
The claim (π s )∗ π s (CUs ) ⊂ CUs follows from the fact that (π s )∗ π s is a projection from X onto X s ⊂ X
and that CUs ⊂ X s . The claim that (π s )∗ π s (CLs ) ⊂ CLs follows from the linearity of the stage constraints.
More precisely, it follows by noting that for any δ p̃ = (Inc ⊗ τ )T δηp and any x ∈ L2n ,
Z
δ p̃T Cx xdt = δηpT Cx π ns (x)
I
Z
= δ p̃T Cx (π ns )∗ π ns (x)dt
(60)
I

holds.

6

Properties C1-C5

The sets D, DUs , and DLs have the following properties
C1) the sets DUs , DLs , and D are closed and convex.
C2) the sets π s (DUs ) and π s (DLs ) are closed and convex.
C3) DUs ⊂ DUs+1 ⊂ D.
C4) DLs ⊃ DLs+1 ⊃ D.
C5) (π s )∗ π s (DUs ) ⊂ DUs , (π s )∗ π s (DLs ) ⊂ DLs .
3

X s is finite-dimensional, thus complete, and hence also closed.

11

We will sketch the proof of Properties C1-C5 below.
The convexity of the sets DUs , DLs , and D follows directly from the linearity of the dynamics.
The fact that the set DLs is closed can be seen by a similar argument used for showing closedness of C˜s
in Sec. 5, that is, showing that X \ DLs is open. The set D can be characterized using the variational equality,
Z
Z
T
− δ ṗ xdt − δpT (Ax + Bu)dt − δp(0)T x0 + δp(T )T xT = 0, ∀δp ∈ Hn ,
(61)
I

I

where Hn denotes the set of functions in L2n having a weak derivative in L2n . Thus, again a similar argument
can be applied to show that D is closed. It follows that DUs is closed, since DUs is defined as the intersection
of D with the closed set X s .
The linearity and surjectivity of the map π s implies that π s (DUs ) and π s (DLs ) are indeed closed (by the
open mapping theorem) and convex.
The inclusion DUs ⊂ DUs+1 ⊂ D for all s follows directly from the fact that X s ⊂ X s+1 ⊂ X and
s
DU = D ∩ X s . The inclusion DLs+1 ⊂ DLs for all s can be seen by noting that the variational equality
in (16) has to hold for variations spanned by more and more basis functions as s increases. The claim
that D is contained in DLs follows from the equivalence of (61) with the formulation in (7). The properties
(π s )∗ π s (DUs ) ⊂ DUs and (π s )∗ π s (DLs ) ⊂ DLs can be shown using the same arguments as in 5, where the latter
relies on the linearity of the dynamics.

7

Convergence results

We prove the following result for the case where I has infinite measure.
Lemma 7.1 Given that the basis functions form an algebra and that the basis functions are dense in the set
of smooth functions with compact support in I, it holds that
lim CUs = lim CLs .

s→∞

s→∞

Proof The idea of the proof is the following: We claim that lims→∞ CUs ⊃ lims→∞ CLs . We assume that the
claim is incorrect and show that this leads to a contradiction. Thus, we choose (x, u, xT ) ∈ lims→∞ CLs ,
such that there exists an open set U (bounded) and a k ∈ {1, 2, . . . , nc }, for which
Z
δv(−Cxk x − Cuk u + bk )dt < 0
(62)
I

holds for all smooth test functions δv : I → R, δv ≥ 0, with support in U , and δv(t0 ) > 0 for some t0 ∈ U .
Due to the smoothness of the test functions, δv(t0 ) > 0 readily implies that there is an open neighborhood
of t0 , denoted by N (t0 ), such that δv(t) > 0, ∀t ∈ N (t0 ). The above integral exists, since δv is bounded,
has compact support and x ∈ L2n , u ∈ L2m . We fix t0 ∈ U and pick one of these variations that is positive,
strictly positive at time t0 , and has support in U , which we name δp. Due to the fact that√the basis functions
are dense in the
√ set of smooth functions with compact support, there exists a sequence δ p̃i that converges
uniformly to δp. Due to the fact that the basis functions form an algebra, δ p̃i lies likewise in the span of
the basis functions. Moreover, for a given  > 0 (small enough) there exits an integer N > 0 such that
||δ p̃i − δp||∞ < C1 
holds for all integers i > N , where C1 > 0 is constant.
We claim that there is an integer p > 0 such that the basis function τp : I → R is nonzero for t0 .
This can be shown by a contradiction argument: If the claim was not true, then all basis functions τ p =
(τ1 , τ2 , . . . , τp ) would be zero at time t0 . The basis functions fulfill Assumption A2, and hence the first order
differential equation τ̇ p = Mp τ p , which is guaranteed to have unique solutions. From τ p (t0 ) = 0, τ̇ p (t0 ) =
12

0 we can infer that τ p (t) = 0 for all t ∈ (0, ∞) is the (unique) solution to τ̇ p = Mp τ p , contradicting
Assumption A1.
Thus, we can choose the basis function τp that is nonzero for t0 and hence also nonzero in a neighborhood around t0 , due to the smoothness of the basis functions. The same applies to the function τp2 , which
is likewise contained in the set of basis functions (the basis functions form an algebra that is closed under
multiplication). Moreover, due to the fact that the basis functions are orthonormal, it follows that
Z
τp2 dt = 1.
(63)
I

In the same way, the function δ p̃i (t)τp2 (t) is also contained in the set of basis functions, is non-negative for
all t ∈ I, and is bounded and integrable for all integers i > 0 (|τ p | is bounded, see 10). Thus, it follows that
the integral
Z
τp2 δ p̃i (−Cxk x − Cuk u + bk )dt

(64)

I

exists. By assumption (x, u, xT ) ∈ lims→∞ CLs , and therefore
Z
0 ≤ τp2 δ p̃i (−Cxk x − Cuk u + bk )dt
ZI
= τp2 δp(−Cxk x − Cuk u + bk )dt
ZI
+ τp2 (δ p̃i − δp)(−Cxk x − Cuk u + bk )dt,

(65)
(66)
(67)

I

where the last term can be bounded by
Z
C1

τp2 |Cxk x + Cuk u − bk |dt.

(68)

I

The above integral is bounded due to the fact that τp2 is bounded and integrable, x ∈ L2n , u ∈ L2m , and that
bk is bounded, and therefore (67) can be made arbitrarily small by sufficiently increasing i. However, this
leads to a contradiction, since (66) is strictly negative, according to (62). This proves that lims→∞ CUs ⊃
lims→∞ CLs . The desired result is then established due to the fact that CUs ⊂ CLs holds for all integers s > 0.

8

Recursive feasibility and closed-loop stability

Proposition 8.1 Provided that the optimization (21) is feasible at time t = 0, it remains feasible for all
times t > 0, and the resulting closed-loop system is guaranteed to be asymptotically stable.
Proof The proof is taken from [2] and included for completeness. The following notation is introduced:
The closed-loop state and input trajectories are denoted by x(t) and u(t). The predicted trajectories are
referred to as x̃(tp |t), ũ(tp |t), where tp > 0 denotes the prediction horizon. For tp = 0, the prediction
matches the true trajectory, that is x̃(0|t) = x(t), ũ(0|t) = u(t) for all t ∈ [0, ∞). The predictions x̃(tp |t),
ũ(tp |t) are obtained by solving (21) subject to the initial condition x0 = x(t), which yields the parameters
ηx and ηu defining x̃(tp |t) and ũ(tp |t) by
x̃(tp |t) = (In ⊗ τ (tp ))T ηx ,

ũ(tp |t) = (Im ⊗ τ (tp ))T ηu .

(69)

In order to highlight the dependence on the initial condition, the resulting optimal cost of (21) is denoted
by J MPC (x(t)).
13

By assumption, (21) is feasible at time t = 0. The resulting trajectories x̃(tp |0), ũ(tp |0) fulfill the
equations of motion, the initial condition x̃(0|0) = x(0), and the constraints and hence, the system evolves
according to x(t) = x̃(t|0), u(t) = ũ(t|0), ∀t ∈ [0, Td ). Due to the time-shift property of the basis functions
implied by Assumption A2, the feasible candidates
x̃(tp + Td |0) = (In ⊗ τ (tp )T )(In ⊗ exp(Ms Td )T )ηx ,
ũ(tp + Td |0) = (Im ⊗ τ (tp )T )(Im ⊗ exp(Ms Td )T )ηu

(70)

for the optimization at time Td can be constructed from the optimizer ηx and ηu at time 0. As a result,
recursive feasibility of (21) follows by induction.
We will show that the optimal cost J MPC acts as a Lyapunov function. The function J MPC is a valid
Lyapunov candidate since J MPC (x) > 0 for all x 6= 0 and J MPC (x) = 0 if and only if x = 0. Due to the fact
that the shifted trajectories x̃(tp +Td |0) and ũ(tp +Td |0) (as defined in (70)) are feasible for the optimization
at time Td , the following upper bound on J MPC (x(Td )) can be established
Z ∞
MPC
J
(x(Td )) ≤
x̃(tp |0)T Qx̃(tp |0) + ũ(tp |0)T Rũ(tp |0)dtp .
(71)
Td

The right-hand side can be rewritten as
Z Td
MPC
J
(x(0)) − x̃(tp |0)T Qx̃(tp |0) + ũ(tp |0)T Rũ(tp |0)dtp ,

(72)

0

resulting in
J MPC (x(Td )) − J MPC (x(0)) ≤
Z Td
x̃(tp |0)T Qx̃(tp |0) + ũ(tp |0)T Rũ(tp |0)dtp .
−

(73)

0

The right-hand side of (73) is guaranteed to be strictly negative, except for x(0) = 0, and thus, by induction,
J MPC (x(kTd )) is strictly decreasing, which concludes the proof.

9

Reduction of the semi-infinite constraint

The following section discusses the reduction of the semi-infinite constraint
al ≤ τ (t)T η ≤ au

(74)

over the interval t ∈ [0, ∞), with al , au ∈ R, al < 0, au > 0 to a compact interval. Thereby we consider
the symmetric case, where |al | = |au | first, before discussing the asymmetric case |al | 6= |au |. It is shown
that the length of this compact time interval depends only on the properties of the basis functions and on
the ratio between |au | and |al |. Both al and au are assumed to be finite.

9.1

The symmetric case

Proposition 9.1 Provided that the basis functions fulfill Assumptions A1 and A2 for all t ∈ [0, ∞) there
exists a positive real number Tc such that
sup |τ (t)T η| = max |τ (t)T η|
t∈[0,Tc ]

t∈[0,∞)

holds for all parameter vectors η ∈ Rs .
14

(75)

Proof Without loss of generality we restrict the parameter vectors to have unit magnitude, that is, |η| = 1.4
We prove the claim in 4 steps. We first derive an exponentially decaying upper bound on the Euclidean
norm of the basis function vector τ . We then use this bound to argue that the basis functions are linearly
independent over the interval [0, Ti ] (the scalar Ti will be determined). The third step consists of constructing
a lower bound on
max |τ (t)T η|
(76)
t∈[0,Ti ]

that holds for all parameter vectors η with |η| = 1. Linear independence of the basis functions on [0, Ti ] will
be used to argue that this lower bound is strictly positive. In the last step, we show that if t is sufficiently
large, |τ (t)T η| will be below this lower bound, which concludes the proof.
Step 1): The fact that Ms is asymptotically stable implies that there exists a quadratic Lyapunov function
that decays exponentially. This provides a means to establish the following bound
|τ (t)|2 ≤ C2 e−c2 t ,

∀t ∈ [0, ∞),

(77)

where C2 > 0, c2 > 0 are constant.
Step 2): We use the Gram matrix to argue that the basis functions are linearly independent over the
interval [0, Ti ]. According to [5, p. 2, Thm. 3] it holds that the basis functions are linearly independent in
the set [0, Ti ] if and only if the matrix
Z Ti
τ τ T dt
(78)
0

has full rank. This is the case if the bilinear form
Z
T
v

Ti

τ τ T dt v

(79)

0

is strictly positive for all v ∈ Rs with |v| = 1. Combining the fact that the basis functions are orthonormal
with the Cauchy-Schwarz inequality, leads to the following lower bound of the above bilinear form,
Z ∞
Z ∞
T 2
|τ |2 dt.
(80)
(v τ ) dt ≥ 1 −
1−
Ti

Ti

Using the upper bound (77) we therefore obtain
Z Ti
Z
T
T
v
τ τ dtv ≥ 1 −
0

∞

C2 e−c2 t dt = 1 −

Ti

C2 −c2 Ti
e
,
c2

(81)

for all v ∈ Rs with |v| = 1. Thus, we fix Ti > 0, such that
1>

C2 −c2 Ti
e
,
c2

(82)

implying that the matrix (78) is positive definite and has therefore full rank. Consequently, the basis functions are guaranteed to be linearly independent on the interval [0, Ti ].
Step 3): We claim that
c∗ := inf max |τ (t)T η|
(83)
|η|=1 t∈[0,Ti ]

is well-defined and strictly positive. To that extent we first prove that the function g : Rs → [0, ∞),
g(η) := max |τ (t)T η|
t∈[0,Ti ]
4

The claim holds trivially for η = 0; in case η 6= 0 we can always normalize η.

15

(84)

is continuous (in fact Lipschitz-continuous). Therefore we consider two parameter vectors η1 and η2 with
g(η1 ) ≤ g(η2 ) (without loss of generality). From the fact that g(η2 ) can be rewritten as g(η1 + (η2 − η1 ))
and by invoking the triangle inequality it can be inferred that
g(η1 + (η2 − η1 )) ≤ max |τ (t)T η1 | + |τ (t)T (η2 − η1 )|.
t∈[0,Ti ]

(85)

By noting that the maximum of the above sum is smaller than the sum of the summand’s maxima it can be
concluded that
g(η2 ) ≤ g(η1 ) + max |τ (t)T (η2 − η1 )|.
(86)
t∈[0,Ti ]

Moreover, by combining the Cauchy-Schwarz inequality and the bound (77) we obtain
p
|g(η2 ) − g(η1 )| ≤ C2 |η2 − η1 |,

(87)

showing that the function g is indeed (Lipschitz) continuous. As a result, the Bolzano-Weierstrass theorem
asserts that the infimum in (83), is attained and well-defined. It remains to argue that c∗ > 0. For the sake
of contradiction, we assume c∗ = 0. This implies the existence of the minimizer η ∗ , with |η ∗ | = 1, which
fulfills
max |τ (t)T η ∗ | = 0.
(88)
t∈[0,Ti ]

As a result, it follows that τ (t)T η ∗ is zero for all t ∈ [0, Ti ], which contradicts the fact that the basis functions
τ are linearly independent on [0, Ti ].
Step 4): From the upper bound (77) and the Cauchy-Schwarz inequality it follows that
p
c2 t
|τ (t)T η| ≤ C2 e− 2 , ∀t ∈ [0, ∞),
(89)
√
and for all η ∈ Rs with |η| = 1. Clearly, c∗ ≤ C2 and therefore we can choose the time Tc such that
c∗ =

p
c2 Tc
C2 e− 2 ,

(90)

implying
sup |τ (t)T η| < c∗ ≤ max |τ (t)T η|,
t∈[0,Ti ]

t∈(Tc ,∞)

(91)

for all η ∈ Rs with |η| = 1. This proves the claim.

9.2

The asymmetric case

Proposition 9.2 Provided that the basis functions fulfill Assumptions A1 and A2 for all t ∈ [0, ∞) there
exists a positive real number Tc such that
al ≤ τ (t)T η ≤ au ,

∀t ∈ [0, Tc ]

(92)

al ≤ τ (t)T η ≤ au ,

∀t ∈ [0, ∞),

(93)

implies
for any η ∈ Rs , where al , au ∈ R, al < 0, au > 0.
Proof We define f˜ := τ T η and establish upper and lower bounds on f˜. Without loss of generality we
assume η 6= 0. Combining Assumption A2 with the Caley-Hamilton theorem leads to
f˜(s) (t) + a1 f˜(s−1) (t) + · · · + as f˜(t) = η T τ (s) (t) + a1 η T τ (s−1) (t) + · · · + as η T τ (t)

= η T M s + a1 M s−1 + · · · + as M τ (t) = 0,
16

(94)
(95)

where a1 , a2 , . . . as are the coefficients of the characteristic polynomial of the matrix M . Thus, the trajectory
f˜ and its time derivatives fulfill the following set of differential equations




τ (0)T
0
1
0
... 0

 0

τ (0)T M T
0
1
... 0 




f,
f
(0)
=
(96)
f˙ =  ..

 η,

..
..
..
..


 .

.
.
.
.
−as −as−1 −as−2 . . . −a1
τ (0)T (M (s−1) )T
{z
}
{z
}
|
|
:=M̂

:=Ĥ

where f := (f˜, f˜(1) , . . . , f˜(s−1) ). The matrix M̂ is Hurwitz and therefore, due to the Lyapunov theorem,
there exists a symmetric matrix P > 0, P ∈ Rs×s that satisfies
P M̂ + M̂ T P + Q = 0

(97)

for any symmetric matrix Q > 0, Q ∈ Rs×s . We fix the positive definite matrix Q and consider the
quadratic Lyapunov function V (t) = f (t)T P f (t), where P satisfies (97). The time derivative of V can be
upper bounded by
λQ T
f Pf
λP
λQ
≤ − P V,
λ

V̇ = −f T Qf ≤ −λQ |f |2 ≤ −

(98)
(99)

where the minimum eigenvalue of Q is denoted by λQ and the maximum eigenvalue of P is denoted by λP .
As a result, this yields the upper bound
λQ

V (t) ≤ V (0)e− λP t ≤ λP λĤ

T Ĥ

λQ

|η|2 e− λP t ,

(100)

T

where λĤ Ĥ denotes the maximum eigenvalue of the matrix Ĥ T Ĥ. According to the proof of Prop. 9.1, we
may choose Ti such that
C2 −c2 Ti
1>
e
,
(101)
c2
implying that the basis functions are linearly independent on the interval [0, Ti ] (see (82) and Prop. 9.1 for
the definition of the constants c2 and C2 ). Linear independence can be used to establish the following lower
bound, c.f. Prop. 9.1:
Z
Z
Ti

Ti

|f˜|2 dt = η T

τ τ T dt η ≥ c3 |η|2 ,

(102)

where the constant c3 > 0 denotes the minimum eigenvalue of the matrix
Z Ti
τ τ T dt.

(103)

0

0

0

Without loss of generality it is assumed that |al | ≤ |au |. Choosing the real number Tc > Ti implies that the
constraint is imposed over the interval [0, Ti ] and therefore
Z Ti
|f˜|2 dt ≤ Ti |au |2 .
(104)
0

Combined with the above upper bound on |η|2 it follows that
|η|2 ≤

Ti |au |2
.
c3

17

(105)

This can be used to upper bound the squared magnitude of f˜(t), that is,
T

λP λĤ Ĥ |au |2 Ti − λQP t
1
2
˜
V (t) ≤
e λ ,
|f (t)| ≤
λP
λ P c3

(106)

where λP denotes the minimum eigenvalue of the matrix P . As a result, we may choose Tc such that the
above upper bound is below |al |2 . This may be achieved by choosing Tc > max{t̂, Ti }, where





|au |
c3 λ P
λP
2 ln
+ ln
t̂ := −
(107)
λQ
|al |
Ti λP λĤ T Ĥ
and ln refers to the natural logarithm.

10

Additional properties

In the following we will discuss some of the properties of the basis functions that fulfill Assumptions A1
and A2.

10.1

Conditions imposed on Ms

The fact that the basis functions are orthonormal on the interval I = (0, T ) implies
Z
d
T
T
τ (T )τ (T ) − τ (0)τ (0) =
(τ τ T )dt
dt
I
= Ms + MsT .

(108)
(109)

In case the interval I has infinite measure, that is, I = (0, ∞), the above formula reduces naturally to
−τ (0)τ (0)T = Ms + MsT .

10.2

(110)

Bounds on the Euclidean norm

In case the interval I has infinite measure, that is, I = (0, ∞), it holds that
|τ (t)| ≤ |τ (0)|,

∀t ∈ [0, ∞).

(111)

This results from
d
(τ (t)T τ (t)) = τ (t)T (Ms + MsT )τ (t)
dt
= −|τ (0)T τ (t)|2 ≤ 0,

(112)
(113)

for all t ∈ [0, ∞), which follows from (110).

10.3

Finite number of minima and maxima

In the following we will argue that a function that is not everywhere zero and spanned by basis functions
fulfilling Assumptions A1 and A2, has a finite number of minima and maxima. The function will be
denoted by f˜ := τ T η, η ∈ Rs , η 6= 0, where the function and the basis functions are defined over the
interval t ∈ I = (0, T ). If I has infinite measure, it follows from Sec. 9 that f˜ takes is maxima and minima
within a compact interval, and hence, without loss of generality, we assume in the following that I has finite
measure. Moreover, due to the Caley-Hamilton theorem it holds that
f˜(s) (t) + a1 f˜(s−1) (t) + · · · + as f˜(t) = 0,
18

(114)

for all t ∈ I, where the ak , k = 1, 2, . . . , s are the coefficients of the characteristic polynomial of the matrix
Ms . According to [6], the time derivative of f˜ has at most s − 1 zeros on any subinterval of length l, where
s
X
ak l k
k=1

k!

< 1.

(115)

This proves readily that f˜ has a finite number of minima and maxima in the interval I.

References
[1] R. Fletcher. Practical Methods of Optimization. John Wiley & Sons, second edition, 1987.
[2] M. Muehlebach and R. D’Andrea. Parametrized infinite-horizon model predictive control for linear
time-invariant systems with input and state constraints. Proceedings of the American Control Conference, pages 2669–2674, 2016.
[3] Y. Nesterov. Squared functional systems and optimization problems. In H. Frenk, K. Roos, T. Terlaky,
and S. Zhang, editors, High Performance Optimization, chapter 17, pages 405–440. Springer, 2000.
[4] J. Nocedal and S. J. Wright. Numerical Optimization. Springer, second edition, 2006.
[5] G. Sansone. Orthogonal Functions. Interscience Publishers, 1959.
[6] S. Yakovenko. On functions and curves defined by ordinary differential equations. Fields Institute
Communications, 24:497–525, 1999.

19

