Nonlinear Stabilization via Control Contraction Metrics: a
Pseudospectral Approach for Computing Geodesics

arXiv:1607.04340v2 [] 4 Nov 2017

Karen Leung

Abstract— Real-time nonlinear stabilization techniques are
often limited by inefficient or intractable online and/or offline
computations, or a lack guarantee for global stability. In this
paper, we explore the use of Control Contraction Metrics
(CCM) for nonlinear stabilization because it offers tractable
offline computations that give formal guarantees for global
stability. We provide a method to solve the associated online
computation for a CCM controller - a pseudospectral method to
find a geodesic. Through a case study of a stiff nonlinear system,
we highlight two key benefits: (i) using CCM for nonlinear
stabilization and (ii) rapid online computations amenable to
real-time implementation. We compare the performance of a
CCM controller with other popular feedback control techniques, namely the Linear Quadratic Regulator (LQR) and
Nonlinear Model Predictive Control (NMPC). We show that
a CCM controller using a pseudospectral approach for online
computations is a middle ground between the simplicity of LQR
and stability guarantees for NMPC.

I. I NTRODUCTION
Control design for a general nonlinear dynamical system
continues to be a challenging problem. In general, there
are three primary benchmarks when constructing an ‘ideal’
nonlinear controller: (i) a globally stabilizing controller, (ii)
tractable offline computations and (iii) fast and tractable
online computations amenable for real-time implementation.
A common solution is to linearize the dynamics about a
particular operating point or desired trajectory and apply a
linear control synthesis techniques, e.g. the Linear Quadratic
Regulator (LQR). The offline and online computations are
tractable however, the system will in general only be locally
stable and may become unstable far from the point of linearization. Thus in general, global stability cannot be guaranteed for LQR. At the other end of the spectrum is Nonlinear
Model Predictive Control (NMPC) which solves a finitehorizon optimal control problem (OCP) at each iteration [1].
Although NMPC generally provides good performance, the
associated online optimization problems are generally nonconvex due to dynamic constraints, and can be difficult to
solve in real time.
An elegant theory of necessary and sufficient conditions
for global stabilizability can be given in terms of control
Lyapunov functions (CLF) [2]. However, the search for
a CLF is generally non-convex, and existing constructive
methods such as backstepping rely on the system having a
particular triangular structure [3]. In general, unless there is
some structure to exploit, the search for a CLF is essentially
intractable.
This work was supported by the Australian Research Council. The authors
are with the Australian Centre for Field Robotics (ACFR), University of
Sydney, NSW 2006, Australia kleu6124@uni.sydney.edu.au,

i.manchester@acfr.usyd.edu.au

Ian R. Manchester

Contraction theory is an attractive alternative because it
combines the simplicity and tractability of linear analysis
with formal guarantees of global stability [4]. The concept of
a Control Contraction Metric (CCM) was introduced in [5],
[6] to generalize contraction analysis to constructive control
design, and gives sufficient conditions for every trajectory
of a nonlinear system to be stabilizable, a property called
universal stabilizability. The conditions are general enough
to be necessary and sufficient for linear systems and feedback
linearizable systems. The offline search for a CCM is a
convex optimization problem [6] and for polynomial systems
can be solved using sum-of-squares (SOS) programming [7].
Given a CCM, a state feedback controller can be constructed via an integration along a shortest path (a geodesic)
between the system’s current state and the target state [5],
[6]. Therefore the online computation of this scheme requires
solving an optimization problem to find a geodesic with
respect to the CCM. This is analogous to NMPC; solving
an optimization problem at each time step. However note
that the CCM approach is generally simpler due to the lack
of dynamic constraints. Indeed, for many nonlinear systems
a state-independent (a.k.a. “flat”) CCM can be found, for
which all geodesics are straight lines [8], but nontrivial metrics can also be found which can offer better performance.
This paper considers the computation of geodesics in the
context of CCM and nonlinear stabilization. In particular,
the application of CCM to a stiff nonlinear system which
are often difficult to control. Indirect methods to finding a
geodesic involves deriving an explicit differential equation
(the “geodesic equation”) via the Euler-Lagrange equation,
however for non-trivial metrics this will be difficult to solve
[9]. Other methods that have been explored recently include
the “Phase Flow Method” [10], fast marching [11], and graph
cuts [12]. Although indirect methods would offer greater
accuracy and confidence that first-order optimality conditions
are met, disadvantages include small radii of convergence
and the need to analytically derive the particular necessary
conditions for each instance of a problem.
Direct methods involve discretizing a shortest path or
optimal control problem, and then solving the resulting
nonlinear problem (NLP) using generic methods and this is
the approach taken in this paper. Popular methods include
multiple shooting (e.g. [13]), direct collocation (e.g. [14]),
and global pseudospectral (e.g. [15], [16], [17]).
In this paper, we follow a well studied Chebyshev pseudospectral approach [18], [19], [20], [16], [21] to solve the
geodesic problem but with modifications to cater for our
specific framework. We offer a criteria that describes the
accuracy of the geodesic. By forming an efficient and rapid

solution to the geodesic problem for a CCM, we enable a
fast construction of a CCM controller amenable to real-time
applications. We will compare the performance of the CCM
controller with the LQR and NMPC framework which are
often used in real-time applications.
The structure of the paper is as follows: in Section II we
define the problem and recall relevant facts about metrics,
geodesics, and contraction analysis; in Section III we detail
the pseudospectral method applied to the geodesic problem.
An example stiff nonlinear system is introduced in Section
IV. In Section V we present the benefits of using a CCM
controller over LQR and NMPC and in Section VI we showcase the efficiency and accuracy of using a pseudospectral
method in approximating a geodesic over multiple shooting.
Finally we offer some brief conclusions.
II. P RELIMINARIES AND P ROBLEM F ORMULATION
In this paper we consider nonlinear control affine systems,
possibly time dependent, of the form:
ẋ = f (x,t) + B(x,t)u

(1)

where x(t) ∈ Rn is the state, u(t) ∈ Rm is the control, and
t ∈ R+ := [0, ∞). The function f : Rn × R+ → Rn is assumed
to be smooth, and B : Rn × R+ → Rn×m has columns bi (x,t),
i = 1, 2, ..., m.
The objective is to design a tractable state-feedback control
policy that can globally exponentially stabilize any feasible
trajectory of the system (1). A system is called universally
exponentially stabilizable if this is possible [6].
A. Metrics and Geodesics
A Riemannian metric equips a smooth manifold – which
for this paper remains Rn but can be more general – with
an inner product on each tangent space, i.e δxT M(x)δx [22],
[23]. This gives a local definition of length and orthogonality.
Consider a smooth curve γ : [0, 1] → Rn with ∂∂ γs 6= 0 ∀s ∈
[0, 1]. Define e(s, γ) = γsT (s)M(γ(s))γs (s) where γs = ∂∂ γs .
Then the length of the curve with end points at γ(0) and
γ(1) is defined as:
`(γ) =

Z 1p

e(s, γ) ds

Z 1

arg min
γ∈C[0, 1]

e(s, γ) ds

and this can be extended to piecewise-smooth curves by
summing over smooth pieces.
A Riemannian metric defines a distance between points –
the Riemannian distance – given by the infimum of lengths
of curves joining the end points. If the resulting metric
space is complete, then the Hopf-Rinow theorem states that a
minimizing path exists, which is a geodesic [23]. Geodesics
are analogues of straight lines in curved space. A special
property of a geodesic is that it has constant speed, i.e. e(s, γ)
is independent of s ∈ [0, 1]. If γ is a geodesic, we have
Z 1

e(s, γ) ds =: E(γ) ,
0

where E(γ) is referred to as the energy of γ. So paths of minimum energy also have minimum length, and for numerical
optimization the energy is preferred due to smoothness.

s.t γ(0) = x∗ (t),

γ(1) = x(t) . (3)

0

In this paper, we will use a pseudospectral method to
discretize this problem, and then apply a quasi-Newton
algorithm to the resulting nonlinear program.
B. Control Contraction Metrics
Contraction analysis, as presented in [4], is based upon
the study of the differential (linearized) dynamics of (1).
Roughly speaking, if all solutions of a system are locally
stable, then all solutions are globally stable [4]. This can be
established by way of a contraction metric, i.e. a Riemannian
metric
M(x,t) for which the associated differential lengths
p
δxT M(x,t)δx shrink exponentially with time.
A CCM is a Riemannian metric for which differential
lengths can be made to shrink by control action and more
details about contraction theory can be found in [6]. A
central result in [6] is that the search for a CCM, M(x,t),
is equivalent to finding a dual CCM W (x,t) = M −1 (x,t).
The search for W (x,t) (and ρ(x,t)) is done by solving the
following Linear Matrix Inequality (LMI):
−Ẇ +WAT + AW − ρBBT < −2λW
where A is the Jacobian matrix and B is the control matrix.
Note that the LMI is jointly convex in W (x,t) and ρ(x,t).
If (1) has polynomial dynamics, then sum-of-squares (SOS)
programming [24] is a computationally tractable method
to find W (x,t) and ρ(x,t). Thus the offline computation
required for a CCM controller is a convex problem and is
tractable through SOS programming.
Using W (x,t) and ρ(x,t) satisfying the LMI, a control
signal u(t) that stabilizes the system to the desired trajectory
x? (t) is given by (4) where u? (t) is the nominal control input.

(2)

0

`(γ)2 =

Classically, a minimal geodesic would be found by using
the Euler-Lagrange equation to obtain necessary conditions
in the form of an ordinary differential equation [9], but
solving the associated two-point boundary value problem is
generally non-trivial. Alternatively, the search for a geodesic
can be expressed as a direct optimization over smooth paths:

u(t) = u? (t) −

1
2

Z 1

ρ(γ(s),t)B(t)T M(γ(s),t)

0

∂γ
ds
∂s

(4)

The integral is computed over γ(s), a geodesic with respect
to the metric M(x,t) connecting x∗ (t) = γ(0) to x(t) = γ(1).
Alternatively, one can utilize the fact that the Riemannian
energy to the target state forms a control Lyapunov function,
and use any control technique requiring a CLF (e.g. pointwise min-norm control) [6, Sec. IV.A]. Either way, the main
computational task to solve in real-time is computation of a
non-trivial minimal geodesic associated with the CCM, i.e. a
solution to Problem (3). It is precisely this computation that
we address in the following sections.
III. P SEUDOSPECTRAL M ETHOD
Since (3) is an infinite dimensional problem over all
smooth curves, the problem needs to be discretized to make
it amenable to numerical solutions. Since geodesics are

smooth paths on smooth manifolds, polynomials are a natural
class of basis functions. In the pseudospectral context, the
discretization points are called collocation points or nodes.
A good choice of polynomial basis, collocation points and
integration scheme are essential for successful application.
In parametrizing the geodesic γ(s), we will represent each
state {x1 , x2 , x3 , ..., xn } along the geodesic with (5). ci j is
the jth coefficient of γi , D is the maximum degree of the
polynomial and φ j (s) is the jth polynomial basis function.
D

γ(s) = [γ1 (s), γ2 (s), ..., γn (s)]T , γi (s) =

∑ ci j φ j (s).

(5)

j=0

The Chebyshev Pseudospectral method is a popular technique in solving optimal control problems [20]. The name
comes from using Chebyshev polynomials to define the
nodes. In particular, the Chebyshev-Gauss-Lobatto (CGL)
nodes are represented by a simple analytic formula and can
be mapped to any interval through an affine transformation.
Further, the Clenshaw-Curtis quadrature (CCQ) scheme is
based on CGL nodes and have weights that are given by
a simple analytic formula [16], [25]. The polynomial basis
associated with the Chebyshev Pseudospectral method are
Lagrange interpolating polynomials. In this paper however,
we analyze a variation of this Chebyshev Pseudospectral
method - we show that for our case study of a stiff nonlinear
system, a Chebyshev polynomial basis is more effective at
solving the geodesic problem than Lagrange polynomials.
A. Discretizaton of the Problem
We follow the pseudospectral method outlined in [20] but
our problem is simpler due to the lack of dynamic constraints
and Chebyshev polynomials are used as the basis, rather than
Lagrange polynomials. Here, we will describe the specifics
of our formulation.
Following on from (5), we can write (6) where Tk (s) is the
kth degree Chebyshev polynomial of the first kind defined on
the [0, 1] interval. c is the decision vector containing all the
coefficients:
c = (c10 , c11 , ...c1D , c20 , c21 , ...c2D , ...cn0 , cn1 , ...cnD ) .
Since the coefficients appear linearly, the differentiation
needed for gradient descent can be derived analytically.
D

γi (c; s) =

∑ ci j Tj (s),

j=0

∂ γi
:= γsi (c; s) =
∂s

D

∑ ci j

j=0

dTk (s)
(6)
ds

Recall D is the degree of the polynomial and n is the
dimension of the state space in (1). So the number of
unknowns for the problem will be (D+1)×n. We will define
N as the stopping index for our CGL nodes {s0 , s1 , ..., sN }.
Integration of a function along a curve is approximated using
the CCQ quadrature scheme evaluated at the CGL nodes. We
refer the reader to [20] for the formulation. Essentially, we
can evaluate (3) for our CCM using this scheme where wk
are the CCQ weights.
Z 1
0

N

γsT M(γ,t)γs ds =

∑ γs (sk )T M(γ(sk ),t)γs (sk )wk
k=0

(7)

The CCQ representation shown in (7) is exact when N = D
and if the integrand is polynomial. However, we cannot
guarantee that M(x,t) is polynomial, even though W (x,t)
is solved through SOS programming. As such, we require
N > D to improve the accuracy of the quadrature scheme.
This is the primary motivation behind selecting a Chebyshev
polynomial basis over Lagrange polynomials. Lagrange polynomials are defined by the number of nodes, meaning that
we always have D = N. Hence we cannot guarantee accuracy
for our quadrature scheme and as we increase the number
of nodes, the size of our problem grows. While Chebyshev
polynomials are defined independent of N. An exploration
for a sufficiently large N for our particular stiff system is
given in Section VI. Chebyshev polynomials also give rise
to a simple endpoint constraint since Tk (0) = 0 if k is even
or Tk (0) = −1 if k is odd, and Tk (1) = 1.
We can now substitute (6) into our optimization problem
(3) and see that our constraints (8) are now linear. Here,
P = [T0 (0), T1 (0), T2 (0), ..., TD (0)] is a row vector and 1 a
row vector of ones of the same length. BDiagn (F) is a block
diagonal matrix with the matrix F down the diagonal n times.


 ∗ 
BDiagn (P)
x (t)
C (c) :=
=0
(8)
c−
BDiagn (1)
x(t)
|
{z
}
| {z }
Ac

bc

As such, the optimization problem in (3) can be expressed
by (9). This is a smooth nonlinear optimization problem with
linear constraints and c is the optimization variable. We will
refer to objective function as E(c) in future sections.
N

argmin
c∈Rn(D+1)

s.t

∑ γs (c, sk )T M(γ(c, sk ),t)γs (c, sk ) wk
k=0

(9)

C (c) = 0

B. Quasi-Newton Optimization
This section will describe important aspects of the optimization procedure used to solve (9), in particular, give an
analytic representation of the gradient. To improve computation costs we employ the standard BFGS quasi-Newton
strategy to approximate the Hessian using the gradient. To
ensure first-order optimality conditions are met, the KKT matrix is solved at each iteration of the optimization algorithm
[26]. Every component in the KKT matrix is readily known
or can be computed with ease; Ac and bc are described in
(8), H is approximated using BFGS and g can be calculated
analytically. We use the relation Ṁ = −MẆ M (the dot
represents a differentiation in any argument). Thus each
component of g is given by (10). (The arguments are omitted
for a concise representation.)
N
∂E
dTk
∂W
= ∑ 2γsT M
wk − γsT M
Mγs T j wk
∂ ci j k=1
ds
∂ xi

c,sk

(10)

The initialization for the optimization variable c is simply
a straight line connecting x∗ (t) to x(t).
Hyperparameters of the quasi-Newton descent algorithm
include the tolerance β = 10−10 to terminate the descent,
initial step size of the line search α0 = 1, termination
condition for the line search c̄ = 0.1 and the rescaling factor

τ = 0.1 for the backtracking line search. However, these
parameters are not too expensive, and are set at typical values
to ensure reliable convergence.
The expensive parameters are the degree of polynomial
used and the number of nodes. The minimal degree and
number of nodes depend on the CCM and the start and end
points. If the CCM is constant, then geodesics are trivially
straight lines. Alternatively, if the CCM is state dependent
then geodesic will tend to be nonlinear and hence a higher
degree and more nodes would be needed.
C. Adaptive Degree Selection
In this section, we discuss a method in selecting the minimal D and N to obtain a geodesic with sufficient accuracy.
Recall that a geodesic has constant energy along the path.
Since this constraint is not enforced during the optimization,
the deviation of e(s, γ) from the constant energy can be used
as an independent validation that the solution found is a
geodesic. The constant energy constraint cannot be enforced
beforehand because the energy is not known a priori. When
our optimization algorithm converges to the solution c∗ , let
E(c∗ ) := E ∗ be the optimal energy.
Then to measure the difference between the energy along
the geodesics and E ∗ we use (11) as a measure of accuracy.
This can be interpreted as relative root-squared-error because
we want to take into account the increases in nonlinearity as
the metric is more nonlinear, or the points are further apart
(i.e., E ∗ is larger).
1
E = ∗
E

is over a convex set.
  
  
−x1 + x3
0
ẋ1
ẋ2  = x2 − x2 − 2x1 x3 + x3  + 0 u
1
1
ẋ3
−x2

(12)

Sum-of-squares programming was carried out on Yalmip
[28] and the semidefinite programming solver MOSEK. The
metric is enforced to match LQR locally about zero [6] by
having W (0,t) = P−1 where P is the solution to the algebraic
Riccati equation (with Q = R = I). A CCM was found and
it is of the form W (x,t) = W0 + W1 x1 + W2 x12 where Wi are
constant matrices and ρ(x,t) = ρ0 + ρ1 x1 + ρ2 x12 . Note that
W and ρ are quadratic in one variable, which is also the
quadratic term in (12). It was shown in [6] that metrics of
this form are complete, and thus a minimal geodesic exists
between every pair of points.
V. C ONTROLLER D ESIGN
In this section, we compare the performance of using
a CCM controller to two popular methods for feedback
control: LQR and NMPC. A CCM controller is constructed
at each time step by solving the geodesic problem using
the proposed pseudospectral method and integrating the
differential controller along this path. The LQR controller
was linearized at the origin while a multiple-shooting NMPC
regime was implemented on ACADO [29], an open-source
C++ software that solves optimal control problems. ACADO
approximates u as piecewise constant, and the state trajectories as piecewise linear.
A. Starting Far From the Origin

s

Z 1

(e(s, γ) − E ∗ )2 ds

(11)

0

For a true geodesic E = 0. Given the facts: (i) set of all
polynomials on [0, 1] spans the set of smooth curves, (ii) by
increasing the maximum degree we enlarge the set of curves
that can be represented, and (iii) the curve of minimal energy
has constant speed, we propose to simply increase the degree
until E is reduced below a certain threshold. It was found
that E < 10−6 was a sufficient tolerance.
As mentioned before, we require D < N in our optimization. Since we also want to keep N small to reduce
computation costs, we consider values for N in the form
of N = D + a. This offers smallest possible N values without
sacrificing too much on computation time. Section VI-A will
explore a suitable value for a.
IV. I LLUSTRATIVE E XAMPLE
To illustrate how a geodesic is computed and used to
construct a CCM controller, the system given in (12) is used
as the subject of this case study. This system was investigated
in [6] and [27]; it is an interesting system because it is not
feedback linearizable, local controllers fail globally and the
dynamics are stiff and highly unstable. In [27], the problem
of uniting a locally optimal and globally stabilizing controller
is explored. In the CCM framework, this problem can be
easily formulated because the search for W (x,t) and ρ(x,t)

First we look at cases where the initial condition is relatively far away from LQR’s point of linearization to investigate the region of stability each controller exhibits. When
starting at [4, 4, 6]T and the target state x∗ (t) = [0, 0, 0]T , the
LQR controller fails to stabilize the system, even though
it is successful at other nearby points, e.g., [5, 5, 5]T . This
highlights the weakness in global stability of LQR. On the
other hand, a CCM controller was successful in stabilizing
the system and Figure 1 illustrates this. Another starting
point further out at [9, 9, 9]T (CCM 9 in Figure 1) was also
tested and the CCM was still able to stabilize the system.
The stiffness of the system is evident, indicated by the very
rapid and slow dynamics. ACADO was not able to compute
a solution from this initial condition, which we believe is
due to the stiffness of the ODE: the slow dynamics require
a long time horizon, and the fast dynamics require a short
sampling time, leading to a very large number of required
nodes over a long time horizon.
B. Starting Close to the Origin
On the other hand, considering initial states close to the
origin where the stiffness is not as prevalent, we find that
all CCM, LQR and NMPC controllers produce very similar
results. For an initial state of x(0) = [1, 1, 1]T and target
state of x∗ (t) = [0, 0, 0]T , the state trajectories are almost
indistinguishable since they are all locally optimal in an
LQ sense. Clearly LQR is superior in terms of computation
time since it only uses elementary matrix operations and

TABLE I
AVERAGE COMPUTATION TIME TO COMPUTE A GEODESIC WITH
DIFFERENT STARTING POINTS AND x∗ (t) = 0.
x(t)
[1, 1, 1]T
[3, 3, 3]T
[5, 5, 5]T
[7, 7, 7]T
[9, 9, 9]T

Fig. 1. States starting at [4, 4, 6]T using LQR and CCM controllers (green
and red). States starting at [9, 9, 9]T (blue) is stabilized using a CCM
controller.

Fig. 2. Comparison of average computation time of a LQR, CCM and
NMPC controller at each iteration during a simulation with states starting
at [1, 1, 1]T .

it is often the technique used when operating close to
the point of linearization. However, Figure 2 compares the
computation time between LQR, CCM and NMPC. We see
that the computation of CCM is reasonably low and with
the advantage that is can also stabilize far from the origin.
As such, CCM poses as a nice middle ground in terms of
computational costs and global stability. To summarize: for
this system, LQR was effective only for regions local to the
point of linearization, while multiple shooting NMPC using
ACADO required many discretization points and did not
succeed for larger initial conditions. CCM gave stabilizing
results with starting points far from and near the origin.
The following section will discuss ways the geodesic can be
computed in order to rapidly construct the CCM controller.
VI. F INDING A G EODESIC
We have seen that using a CCM controller was effective
in stabilizing (12) while LQR and NMPC was unsuccessful.
This section will two explore ways to compute the geodesics
which allows the construction of a CCM controller possible.
As such, we shall see that the proposed method, as outlined
in Section III-A, is very fast and accurate, making it possible
for real-time applications.
A. Pseudospectral Method
First, we numerically validate our choice of Chebyshev
polynomials over Lagrange interpolating polynomials. For
relatively close end points that do not require high degrees,
Chebyshev polynomials yield faster computation times than
using Lagrange polynomials. For example, computation time

D
4
4
5
6
7

Avg. comp. time
0.00811
0.00848
0.01487
0.02069
0.02542

for Lagrange polynomials was on the order of 10−2 while
it was on the order of 10−3 for Chebyshev polynomials
(with endpoints, [0, 0, 0]T and [1, 1, 1]T ). This implies that
Lagrange polynomials will perform even worse for higher
degrees/nodes As such, since a Chebyshev polynomial basis
offer more rapid convergence, this is our choice of basis.
First, a value for a must be fixed. With starting states
at [9, 9, 9]T , the effect on E was tested with different
values of a and this is given in Figure 3 (left). a = 4
was chosen because it gave results that were below the
tolerance (10−6 ) for low degrees without a huge cost in
computation time (Figure 3 right). The results were implemented in Julia, running on a 2012 MacBook Pro with a
2.90GHz Intel Core i7 processor. The code is available at
https://github.com/karenl7/GeodesicCCM.
With a = 4 fixed, the degree, E and the average computation time for various initialization points are given in
Table I. It can be seen that only a degree 7 polynomial was
needed for a far away starting point at [9, 9, 9]T . As a point of
comparison for the same end points N = D = 70 was needed
for Lagrange polynomials to achieve acceptable accuracy,
at an average computation time of 19.7 seconds. Hence we
choose Chebyshev polynomials over Lagrange polynomials.
Since computation time depends on the degree and number
of nodes used, possible offline computations to speed up the
online computation could involve pre-computing optimal D
and a values for regions in the state space. These values can
be accessed, rather than found adaptively during real-time
implementation. However, this could get very complex for
higher dimensional problems and when x∗ (t) is not constant.
Nonetheless, the proposed Chebyshev pseudospectral
method with Chebyshev polynomials as basis functions offer
a very quick and efficient way of approximating a geodesic
with a sufficient accuracy.
B. Multiple Shooting using ACADO
The geodesic problem can also be expressed as a simple
OCP by introducing trivialR dynamical constraints and this
is given by ẋ = u, J = 01 uT M(x)u dt. The computation
time to find the geodesic using ACADO was significantly
greater than the Chebyshev pseudospectral method. Even
converging to a solution for a path connecting [1, 1, 1]T to
[0, 0, 0]T took multiple orders of magnitude greater than the
pseudospectral method and with significantly less accuracy.
For example, with 100 segments, it took roughly 8.8 seconds
and E > 10−2 . Increasing the number of segments would
certainly improve the accuracy, but this would severely
increase computation time, and would be impractical for real-

Fig. 3. Left: E versus D with various values of a. End points are x(t) = [9, 9, 9]T and x∗ (t) = [0, 0, 0]T . Right: Average computation time versus D for
same end points.

time implementation. This is perhaps due to the need to invert
W (x,t) to obtain M(x,t) and take derivatives in the ACADO
toolkit which in general is expensive. However, with the
pseudospectral framework, we could utilize the formula for
the gradient. As such, compared to the multiple shooting
method in ACADO, the proposed pseudospectral method is
significantly faster and more accurate even with less nodes.
VII. CONCLUSIONS
We examined the online computation required for a CCM
controller. The benefits associated with a CCM controller
include tractable offline computations, guarantees on global
stability and, as presented in this paper, rapid and tractable
online computations. It was shown that for a particular stiff
system, a CCM controller was able to stabilize the system
beyond LQR’s region of stability and that the stiffness made
multiple-shooting for NMPC very inefficient. Thus CCM has
the potential to compute nonlinear stabilizing controllers in
real-time for a larger range of initialization points. A pseudospectral method using Chebyshev polynomials to solve
the geodesic problem was found to be fast and tractable,
amenable for real-time implementation. We conclude that the
CCM/pseudospectral method offers a viable alternative for
certain difficult nonlinear stabilization problems, since it encompasses the simplicity of LQR and the global performance
of NMPC.
R EFERENCES
[1] F. Allgöwer and A. Zheng, Nonlinear model predictive control.
Birkhauser, 2012, vol. 26.
[2] E. D. Sontag, “A universal construction of Artstein’s theorem on
nonlinear stabilization,” Systems & Control Letters, vol. 13, no. 2,
pp. 117–123, Aug. 1989.
[3] M. Krstic, I. Kanellakopoulos, and P. Kokotovic, Nonlinear and
adaptive control design. Wiley, 1995, vol. 222.
[4] W. Lohmiller and J.-J. E. Slotine, “On Contraction Analysis for Nonlinear Systems,” Automatica, vol. 34, no. 6, pp. 683–696, June 1998.
[5] I. R. Manchester and J.-J. E. Slotine, “Control Contraction Metrics and
Universal Stabilizability,” in Proceedings of the IFAC World Congress,
Cape Town, South Africa, 2014.
[6] ——, “Control contraction metrics: Convex and intrinsic criteria for
nonlinear feedback design,” IEEE Transactions on Automatic Control,
In Press 2017.
[7] E. M. Aylward, P. A. Parrilo, and J. E. Slotine, “Stability and
robustness analysis of nonlinear systems via contraction metrics and
SOS programming,” Automatica, vol. 44, no. 8, pp. 2163–2170, 2008.
[8] I. R. Manchester, J. Z. Tang, and J. E. Slotine, “Unifying classical and
optimization-based methods for robot tracking control with control
contraction metrics,” International Symposium on Robotics Research
(ISRR), pp. 1 – 16, 2015.
[9] V. I. Arnold, Mathematical methods of classical mechanics. Springer
Science & Business Media, 1989, vol. 60.

[10] L. Ying and E. J. Candes, “Fast geodesics computation with the phase
flow method,” Journal of computational physics, vol. 220, no. 1, pp.
6–18, 2006.
[11] R. Kimmel and J. A. Sethian, “Computing geodesic paths on manifolds,” Proceedings of the National Academy of Sciences, vol. 95,
no. 15, pp. 8431–8435, 1998.
[12] Y. Boykov and V. Kolmogorov, “Computing geodesics and minimal
surfaces via graph cuts,” in Proc. IEEE International Conference on
Computer Vision, 2003.
[13] M. Diehl, H. J. Ferreau, and N. Haverbeke, “Efficient numerical
methods for nonlinear MPC and moving horizon estimation,” in
Nonlinear Model Predictive Control. Springer, 2009, pp. 391–417.
[14] D. A. Benson, G. T. Huntington, T. P. Thorvaldsen, and A. V.
Rao, “Direct trajectory optimization and costate estimation via an
orthogonal collocation method,” Journal of Guidance, Control, and
Dynamics, vol. 29, no. 6, pp. 1435–1440, 2006.
[15] G. N. Elnagar and M. A. Kazemi, “Pseudospectral Chebyshev optimal
control of constrained nonlinear dynamical systems,” Computational
Optimization and Applications, vol. 11, no. 2, pp. 195–217, 1998.
[16] F. Fahroo and I. M. Ross, “Direct trajectory optimization by a
Chebyshev pseudospectral method,” Journal of Guidance, Control, and
Dynamics, vol. 25, no. 1, pp. 160–166, 2002.
[17] D. Garg, M. Patterson, W. W. Hager, A. V. Rao, D. A. Benson, and
G. T. Huntington, “A unified framework for the numerical solution of
optimal control problems using pseudospectral methods,” Automatica,
vol. 46, no. 11, pp. 1843–1851, 2010.
[18] P. Williams, “Application of pseudospectral methods for receding
horizon control,” Journal of Guidance, Control, and Dynamics, vol. 27,
no. 2, pp. 310–314, 2004.
[19] F. Fahroo and I. M. Ross, “Pseudospectral methods for infinite-horizon
nonlinear optimal control problems,” Journal of Guidance, Control,
and Dynamics, vol. 31, no. 4, pp. 927–936, 2008.
[20] Q. Gong, I. M. Ross, and F. Fahroo, “A chebyshev pseudospectral
method for nonlinear constrained optimal control problems,” Joint 48th
IEEE Conference on Decision and Control and 28th Chinese Control
Conference, 2009.
[21] I. M. Ross and M. Karpenko, “A review of pseudospectral optimal
control: from theory to flight,” Annual Reviews in Control, vol. 36,
no. 2, pp. 182–197, 2012.
[22] J. M. Lee, Riemannian manifolds: an introduction to curvature.
Springer Science & Business Media, 2006, vol. 176.
[23] W. M. Boothby, An introduction to differentiable manifolds and
Riemannian geometry. Academic press, 1986.
[24] P. A. Parrilo, “Semidefinite programming relaxations for semialgebraic
problems,” Mathematical programming, vol. 96, no. 2, pp. 293–320,
2003.
[25] L. N. Trefethen, “Is Gauss quadrature better than Clenshaw-Curtis?”
SIAM review, vol. 50, no. 1, pp. 67–87, 2008.
[26] J. Nocedal and S. Wright, Numerical optimization. Springer Science
& Business Media, 2006.
[27] V. Andrieu and C. Prieur, “Uniting two control lyapunov functions
for affine systems,” IEEE Transactions on Automatic Control, vol. 55,
no. 8, pp. 1923–1927, 2010.
[28] J. Löfberg, “Yalmip: A toolbox for modeling and optimization in
matlab,” in Computer Aided Control Systems Design, 2004 IEEE
International Symposium on. IEEE, 2004, pp. 284–289.
[29] B. Houska, H. J. Ferreau, and M. Diehl, “ACADO toolkit-An opensource framework for automatic control and dynamic optimization,”
Optimal Control Applications and Methods, vol. 32, no. 3, pp. 298–
312, May 2011.

