Symmetry reduction for dynamic programming ?
John Maidens a , Axel Barrau b , Silvère Bonnabel c , and Murat Arcak a
a

Department of Electrical Engineering & Computer Sciences, University of California, Berkeley, 253 Cory Hall, Berkeley,
CA, 94720 USA

arXiv:1801.03237v1 [cs.SY] 10 Jan 2018

b
c

SAFRAN TECH, Groupe Safran, Rue des Jeunes Bois - Châteaufort, 78772 Magny Les Hameaux, France

MINES ParisTech, PSL Research University, Centre for Robotics, 60 bd Saint-Michel, 75006 Paris, France

Abstract
We present a method of exploiting symmetries of discrete-time optimal control problems to reduce the dimensionality of
dynamic programming iterations. The results are derived for systems with continuous state variables, and can be applied to
systems with continuous or discrete symmetry groups. We prove that symmetries of the state update equation and stage costs
induce corresponding symmetries of the optimal cost function and the optimal policies. We then provide a general framework
for computing the optimal cost function based on gridding a space of lower dimension than the original state space. This
method does not require algebraic manipulation of the state update equations; it only requires knowledge of the symmetries
that the state update equations possess. Since the method can be performed without any knowledge of the state update map
beyond being able to evaluate it and verify its symmetries, this enables the method to be applied in a wide range of application
problems. We illustrate these results on two six-dimensional optimal control problems that are computationally difficult to
solve by dynamic programming without symmetry reduction.
Key words: Optimal control; Dynamic programming; Invariant systems; Model reduction; Path planning; Medical applications

1

Introduction

The dynamic programming algorithm for computing optimal control policies has, since its development, been
known to suffer from the “curse of dimensionality” (Bellman, 1957). Its applicability in practice is typically limited to systems with four or five continuous state variables because the number of points required to grid a
space of n continuous state variables increases exponentially with the state dimension n. This complexity has
led to a collection of algorithms for approximate dynamic
programming, which scale to systems with larger state
dimension but lack the guarantees of global optimality
of the solution associated with the original dynamic programming algorithm (Bellman and Dreyfus, 1959; Bertsekas, 2012; Powell, 2007, 2016).
? A preliminary version of this work was presented at the
2017 American Control Conference. Corresponding author
J. Maidens.
Email addresses: maidens@eecs.berkeley.edu (John
Maidens), axel.barrau@mines-paristech.fr (Axel
Barrau), silvere.bonnabel@mines-paristech.fr (Silvère
Bonnabel), arcak@eecs.berkeley.edu (Murat Arcak).

Preprint submitted to Automatica

In practice, many real-world systems exhibit symmetries
that can be exploited to reduce the complexity of system
models. Symmetry reduction has found applications in
fields ranging from differential equations (Clarksonz and
Mansfield, 1994; Bluman and Kumei, 2013) to model
checking (Emerson and Sistla, 1996; Kwiatkowska et al.,
2006). In control engineering, symmetries have been exploited to improve control of mechanical systems (Marsden et al., 1990; Bloch et al., 1996; Bullo and Murray,
1999), develop more reliable state estimators (Barrau
and Bonnabel, 2014), study the controllability of multiagent systems (Rahmani et al., 2009) and to reduce the
complexity of stability and performance certification for
interconnected systems (Arcak et al., 2016; Rufino Ferreira et al., 2017). Symmetry reduction has also been applied to the computation of optimal control policies for
continuous-time systems in (Grizzle and Marcus, 1984;
Ohsawa, 2013) and Markov decision processes (MDPs)
in (Zinkevich and Balch, 2001; Narayanamurthy and
Ravindran, 2007).
In this paper, we present a theory of symmetry reduction for the optimal control of discrete-time, stochastic
nonlinear systems with continuous state variables. This
reduction allows dynamic programming to be performed

11 January 2018

crete time systems. Then we introduce our problem and
provide the reader with a primer on the classical theory
of symmetries. We also introduce the notion of invariant
control systems with invariant costs.

in a lower-dimensional state space. Since the computational complexity of a dynamic programming iteration
increases exponentially with state dimension, this reduction significantly decreases computational burden. Further, our proposed method does not rely on an explicit
transformation of the state update equations, making
the method applicable in situations where a such a transformation is difficult or impossible to find analytically.

2.1

We begin by introducing dynamic programming for finite horizon optimal control following the notation of
(Bertsekas, 2005). We consider a discrete-time dynamical system

We present two theorems that summarize our method
of symmetry reduction. Theorem 4 describes how symmetries of the system dynamics imply symmetries of the
optimal cost and optimal policy functions. Theorem 6
then describes a method of computing the cost function
based on reduced coordinate system that depends on
fewer state variables.

k = 0, 1, . . . , N − 1

xk+1 = fk (xk , uk , wk ),

(1)

where xk ∈ X ⊆ Rn is the system state, uk ∈ U ⊆ Rm is
the control variable to be chosen at time k, wk ∈ W ⊆
R` are independent continuous random variables each
with density pk , and N ∈ Z+ is a finite control horizon.
Associated with this system is an additive cost function

This paper builds on the work we presented in the conference paper (Maidens et al., 2017). The most substantial
improvement is the additional theoretical results presented in Sections 2.3 and 3.2. The conference version
presented an ad hoc symmetry reduction for a magnetic
resonance imaging (MRI) application, but did not provide a general methodology for computing the coordinate reduction. This paper addresses this shortcoming
by presenting a general method based on the moving
frame formalism, which leads to the general symmetry
reduction formula presented in Theorem 6. Additionally, the MRI example has been reworked to match this
new formalism, and the numerical implementation and
graphs of the numerical solution have been improved. We
have also included two new extensions of this formalism
to the case of equivariant costs in Section 3.3 and to the
synchronization problem of stochastic dynamic systems
on matrix groups in Section 3.4, along with examples to
illustrate the algorithm in these contexts.

gN (xN ) +

N
−1
X

gk (xk , uk , wk )

k=0

that we wish to minimize through our choice of
uk . We define a control system to be a tuple S =
QN −1
(X , U, W, p, f, g, N ) where p =
k=0 pk is the joint
density of the random variables wk .
We consider a class of control policies π = {µ0 , . . . , µN −1 }
where µk : X → U maps observed states to admissible
control inputs. Given an initial state x0 and a control
policy π, we define the expected cost under this policy as
"
Jπ (x0 ) = E gN (xN ) +

This paper is organized as follows: in Section 2 we introduce notation and provide background information both
on dynamic programming for optimal control, and on
the mathematical theory of symmetries. In Section 3, we
derive our main theoretical results, that is, we prove that
control system symmetries induce symmetries of the optimal cost function and optimal control policy, and then
leverage the result to present a general method of performing dynamic programming in reduced coordinates.
In Section 4 we apply the algorithm to a cooperative control problem for two Dubins vehicles using a Lie group
formulation. In Section 5 we apply symmetry reduction
to compute the solution of an optimal control problem
arising in dynamic MRI acquisition. Code to reproduce
the computational results in this paper is available at
https://github.com/maidens/Automatica-2017.
2

Dynamic programming for optimal control of
stochastic systems

N
−1
X

#
gk (xk , µk (xk ), wk ) .

k=0

An optimal policy π ∗ is defined as one that minimizes
the expected cost:
Jπ∗ (x0 ) = min Jπ (x0 )
π∈Π

where Π denotes the set of all admissible control policies.
The optimal cost function, denoted J ∗ (x0 ), is defined to
be the expected cost corresponding to an optimal policy.
As in (Bertsekas, 2005), we use min to denote the infimum value regardless of whether there is a policy π ∈ Π
that achieves a minimum. In the example problems presented in Sections 4 and 5, the existence of a minimum
is guaranteed by compactness or finiteness arguments
respectively. In the general case, the optimal cost function J ∗ can be computed using the dynamic programming algorithm regardless of the existence of minimizers, but the existence of an optimal policy π ∗ requires
that a minimum be achieved for each xk ∈ X .

Dynamic Programming and Symmetries

In this section, we first recall the main features of dynamic programming for optimal control of stochastic dis-

2

wk ∈ W we have:

We quote the following result due to Bellman from (Bertsekas, 2005):

φ−1 ◦ fk (φ(xk ), χ(uk ), ψ(wk )) = fk (xk , uk , wk ),
k = 0, 1, . . . , N − 1
gk (φ(xk ), χ(uk ), ψ(wk )) = gk (xk , uk , wk ),
k = 0, 1, . . . , N − 1,
gN (φ(xN )) = gN (xN ), and
pk (ψ(wk ))| det Dψ(wk )| = pk (wk )
k = 0, 1, . . . , N − 1

Proposition 1 (Dynamic Programming) For every initial state x0 , the optimal cost J ∗ (x0 ) is equal to
J0 (x0 ), given by the last step of the following algorithm,
which proceeds backward in time from period N − 1 to
period 0:
JN (xN ) = gN (xN )
"
Jk (xk ) = min E gk (xk , uk , wk ) + Jk+1
uk ∈U




fk (xk , uk , wk )

#

where Dψ denotes the Jacobian of ψ.

(2)

The rationale is simple: For any fixed α ∈ G, consider
the change of variables Xk = φα (xk ), Uk = χα (uk ),
Wk = ψα (wk ). Then, we have

k = 0, 1, . . . , N − 1,

where the expectation is taken with respect to the probability distribution of wk defined by the density p. Furthermore, if there exists u∗k minimizing the right hand
side of (2) for each xk and k, then the policy π ∗ =
{µ∗0 , . . . , µ∗N −1 } where µ∗k (xk ) = u∗k is optimal.

Xk+1 = fk (Xk , Uk , Wk ),

and for k = 0, 1, . . . , N −1 we have also gk (Xk , Uk , Wk ) =
gk (xk , uk , wk ). As a result, if u1 , . . . , uN −1 is a series
of controls that minimize J(x0 ), then one can expect
U1 , . . . , UN −1 to minimize J(X0 ), under some assumptions on the noise. As a result, the optimal control
problem needs only be solved once for all initial conditions belonging to the set {φα (x0 )|α ∈ G}, reducing the
initial n dimensional problem to a n − r dimensional
problem. The present paper derives a theory for such
symmetry reduction in dynamic programming, and
provides various examples of engineering interest.

The intermediate functions Jk (xk ) for k > 0 computed
in this manner represent the optimal cost of the tail
subproblem beginning at xk . The optimal cost of the
entire problem is given by the function J ∗ (x0 ) = J0 (x0 )
obtained when this recursion terminates.

2.2

k = 0, 1, . . . , N − 1,

Invariant system with invariant costs

2.3

We first recall the definition of a transformation
group for a control system, as in (Martin et al., 2004;
Jakubczyk, 1998; Respondek and Tall, 2002). See Olver
(1999) for the more general theory.

Cartan’s moving frame method

To find a reduced coordinate system in which to perform dynamic programming, we will use the moving
frame method of Cartan (Cartan, 1937). In general, this
method only results in a local coordinate transformation
as it relies on the implicit function theorem. In this paper we will focus on a transformation within a single coordinate chart. For many practical problems, including
both examples in this paper, the transformation computed using this method extends to all of X with the exception of a lower-dimensional submanifold of the state
space. In such cases, only a single coordinate chart is required for the purpose of gridding the entire state space
for dynamic programming.

Definition 2 (Transformation group) A
transformation group on X × U × W is set of tuples
hα = (φα , χα , ψα ) parametrized by elements α of a Lie
group G having dimension r, such that the functions
φα : X → X , χα : U → U and ψα : W → W are all C 1
diffeomorhpisms and satisfy:
• φe (x) = x, χe (u) = u, ψe (w) = w when e is the
identity of the group G and
• φa∗b (x) = φa ◦ φb (x), χa∗b (u) = χa ◦ χb (u), ψa∗b (x) =
ψa ◦ ψb (x) for all a, b ∈ G where ∗ denotes the group
operation and ◦ denotes function composition.

We briefly introduce the moving frame method following
the presentation in (Bonnabel et al., 2008). Consider
an r-dimensional transformation group (with r ≤ n)
acting on X via the diffeomorphisms (φα )α∈G . Assume
we can split φα as (φaα , φbα ) with r and n − r components
respectively so that φaα is an invertible map. Then, for
some c in the range of φa , we define a coordinate cross
section to the orbits C = {x : φae (x) = c}. This cross
section is an n−r-dimensional submanifold of X . Assume
moreover that for any point x ∈ X , there is a unique
group element α ∈ G such that φα (x) ∈ C. Such α will

To simplify notation we will sometimes suppress the subscripts α. In the present paper, we will consider the following class of systems and cost functions.
Definition 3 (Invariant control system with invariant costs) A control system S is G-invariant with
G-invariant costs if for all α ∈ G, xk ∈ X , uk ∈ U and

3

any xk ∈ X, and uk ∈ U we have

be denoted γ(x), and the map γ : X → G will be called
moving frame.

"

#

E gk (xk , uk , wk ) + Jk+1 (fk (xk , uk , wk ))

A moving frame can be computed by solving the normalization equation:

"

Z
=

#
gk (xk , uk , wk ) + Jk+1 (fk (xk , uk , wk )) pk (wk )dwk

W

"

Z

φaγ(x) (x) = c.

#
gk (φ(xk ), χ(uk ), ψ(wk )) + Jk+1 (φ−1 ◦ fk (φ(xk ), χ(uk ), ψ(wk ))) pk (wk )dwk

=
W

"

Z
=

#
gk (φ(xk ), χ(uk ), ψ(wk )) + Jk+1 (fk (φ(xk ), χ(uk ), ψ(wk ))) pk (wk )dwk

W

Define the following map ρ : X → Rn−r as

"

Z
=

#
gk (φ(xk ), χ(uk ), ψ(wk )) + Jk+1 (fk (φ(xk ), χ(uk ), ψ(wk ))) pk (ψ(wk ))| det Dψ(xk )|dwk

W

ρ(x) = φbγ(x) (x).

#

"

Z
=

gk (φ(xk ), χ(uk ), w̃k ) + Jk+1 (fk (φ(xk ), χ(uk ), w̃k )) pk (w̃k )dw̃k
ψ(W)

Note that, for all α ∈ G we have ρ(φα (x)) = ρ(x), that
is, the components of ρ are invariant to the group action
on the state space. Further, due to our assumptions, the
restriction of ρ to C is injective. We denote this restricted
function ρ̄, and it will serve as a reduced coordinate system to solve the invariant optimal control problem.

3

"

Z
=

#
gk (φ(xk ), χ(uk ), wk ) + Jk+1 (fk (φ(xk ), χ(uk ), wk )) pk (wk )dwk

W

"

#

= E gk (φ(xk ), χ(uk ), wk ) + Jk+1 (fk (φ(xk ), χ(uk ), wk ))

where the change of variables w̃k is defined via w̃k =
ψ(wk ) and the tildes are subsequently dropped. Therefore, from the one-step dynamic programming principle
we have

Main Results
"

#

Jk (xk ) = min E gk (xk , uk , wk ) + Jk+1 (fk (xk , uk , wk ))
uk ∈U

In order to combat the “curse of dimensionality” associated with performing dynamic programming in highdimensional systems, we describe a method to reduce
the system’s dimension by exploiting symmetries in the
dynamics and stage costs.

"

#

= min E gk (φ(xk ), χ(uk ), wk ) + Jk+1 (fk (φ(xk ), χ(uk ), wk ))
uk ∈U

"
=

#

min E gk (φ(xk ), ũk , wk ) + Jk+1 (fk (φ(xk ), ũk , wk ))

ũk ∈χ(U )

= Jk (φ(xk )).

3.1

Symmetries imply equivalence classes of optimal
policies

Thus J ∗ = J ∗ ◦ φ. Now, if π ∗ = {µ∗0 , . . . , µ∗N −1 } is an
optimal policy and we denote x̃k = φ(xk ) then for any
k ∈ {0, . . . , N − 1} we have

Theorem 4 (Symmetries of the optimal cost and
policy) Let G be a group and let S be a G-invariant
control system with G-invariant costs. Then the optimal
cost functions Jk (x0 ) satisfy the symmetry relations

Jk (x̃k ) = Jk (xk )
"

#

= E g(xk , µ∗k (xk ), wk ) + Jk+1 (fk (xk , µ∗k (xk ), wk ))
"

#

= E gk (φ(xk ), χ(µ∗k (xk )), wk ) + Jk+1 (fk (φ(xk ), χ(µ∗k (xk )), wk ))
"

Jk = Jk ◦ φα

#

= E gk (φ(xk ), χ ◦ µ∗k ◦ φ−1 (φ(xk )), wk ) + Jk+1 (fk (φ(xk ), χ ◦ µ∗k ◦ φ−1 (φ(xk )), wk ))
"

#

= E gk (x̃k , χ ◦ µ∗k ◦ φ−1 (x̃k ), wk ) + Jk+1 (fk (x̃k , χ ◦ µ∗k ◦ φ−1 (x̃k ), wk )) .

for any k = 0, . . . , N and any α ∈ G. Furthermore, if
π ∗ = {µ∗0 , . . . , µ∗N −1 } is an optimal policy then so is π̃ ∗ :=
∗
−1
{χα ◦ µ∗0 ◦ φ−1
α , . . . , χα ◦ µN −1 ◦ φα } for any α ∈ G.

Thus π̃ ∗ := {χ ◦ µ∗0 ◦ φ−1 , . . . , χ ◦ µ∗N −1 ◦ φ−1 } is an
optimal policy. 2
3.2

PROOF. We prove this by induction on k proceeding
backward from the base case k = N . First, note that

Dynamic programming can be performed using reduced coordinates

Theorem 4 readily implies the problem can be reduced,
as all states along an orbit of G are equivalent in terms
of cost, and that there are equivalence classes of optimal
policies. So it suffices to only consider the cost corresponding to a single representative of each equivalence

JN (xN ) = gN (xN ) = gN (φ(xN )) = JN (φ(xN )).
Now, suppose that for some k ∈ {0, . . . , N − 1} we have
Jk+1 (xk+1 ) = Jk+1 (φ(xk+1 )) for all xk+1 ∈ X . Then for

4

class, and to find a single representative of the optimal
policy within each class. This can now easily be done
using the injective map ρ̄ : C → Rn−r .

want the cost function J to be equivariant too, we will
need ϕα (·) to be linear. Thus we will simply assume that
ϕα is of the form ϕα (J) = l(α)J, that is, it is a scaling
of the cost, where l : R>0 → R>0 . For simplicity’s sake,
we consider here the problem to be noise free. Along the
lines of the preceding sections it is easily proved that

For x̄ ∈ ρ̄(C) ⊂ Rn−r , let z ∈ C be such that x̄ = ρ̄(z),
and define
J¯k (x̄) = Jk (z).
The following result shows that the functions Jk on the
n-dimensional space X ⊆ Rn are completely determined
by the values of J¯k on the subset ρ̄(C) of Rn−r .

Jk (φ(x)) = ϕ ◦ Jk (x)
as already noticed in (Alvarez and Stokey, 1998) for
the case of homogeneous costs. Symmetry reduction can
then be applied. We now give two tutorial examples.

Corollary 5 For any x ∈ X and k = 0, . . . , N , the cost
function Jk for the full problem can be computed in terms
of the lower-dimensional cost function J¯k as

Example 1 Consider the linear system
xx+1 = Axk + Buk

Jk (x) = Jk (φγ(x) (x)) = J¯k (x̄),

with quadratic costs gk = xTk Qxk + uTk Ruk . The system is invariant to scalings, φα (x) = αx, χα (u) = αu,
and the cost is equivariant letting ϕα (J) = αJ, where
α ∈ G = R>0 . The unit sphere is a cross section to
the orbits, and the normalization equation yields γ(x) =
1/||x||. Applying the results above, we see that the controls
that minimize J(x0 ) are ||x0 ||u∗1 , · · · , ||x0 ||u∗N −1 , where
u∗1 , · · · , u∗N −1 are those minimizing J( ||xx00 || ). This agrees
with the well known fact that the optimal controller for
the problem above is the linear quadratic controller, and
is indeed of the linear form uk = −Fk xk .

where x̄ := ρ̄(φγ(x) (x)) is well defined as φγ(x) (x) ∈ C.
It is thus sufficient to have evaluated J¯ at all points of
ρ̄(C) ⊂ Rn−r to be able to instantly evaluate J at any
−1
point of X . An optimal policy π̄ ∗ = {µ̄∗k : ρ̄(C) → U}N
k=0
can also be lifted to an optimal policy on the original
state space via this method:
µ∗k (x) = µ̄∗k (ρ̄(φγ(x) (x))).
This allows optimal trajectories of the system to be computed in the original coordinates X via the lifted policy
π ∗ = {µ∗0 , . . . , µ∗N −1 }.

Example 2 Consider the following system and costs
xx+1 = Axk + Buk ,

Theorem 6 (Dynamic programming in reduced
coordinates) The reduced coordinates are in one to
one correspondance with the cross-section C. For any
x̄ ∈ ρ̄(C), let z ∈ C satistisfy ρ̄(z) = x̄. Then in the reduced coordinates, the sequence J¯k can be computed recursively via

where ||uk ||1 denotes the L1 norm of uk and h is a map
satisfying h(ax) = ah(x) for a > 0. Such costs may
arise when one tries to force some controls to zero to
create sparsity, a method known as L1 regularization.
This problem is challenging, particularly for nonconvex
h. But according to the theory above, it is sufficient to
solve it numerically for initial conditions lying on the
unit sphere of the state space.



J¯k (x̄) = minuk ∈U E gk (z, uk , wk ) + J¯k+1 (ρ(fk (z, uk , wk ))) .

3.4

PROOF. We have

uk ∈U

uk ∈U

3.3

Optimal formation control on Lie groups

We now apply the theory presented in Section 2.2 where
the state space X ⊆ Rn is the Cartesian product of
matrix Lie groups. Note that, straightforward modifications arise along the way as the state space and noise
space are not vector spaces as in the theory above. The
methodology is then applied to the synchronization of
two non-holonomic cars in the presence of uncertainties.

"
#


J¯k (x̄) = Jk (z) = min E gk (z, uk , wk ) + Jk+1 fk (z, uk , wk )
"
#


= min E gk (z, uk , wk ) + J¯k+1 ρ ◦ fk (z, uk , wk ) .

gk = h(xk ) + ||uk ||1

2

Case of equivariant costs
We model the system as a collection of K agents, where
the state of each agent evolves on a r-dimensional matrix
Lie group G. We assume that the evolution of the state
of agent j proceeds according to the equation

So far, we have considered the costs to be invariant under transformation. We now briefly discuss how this can
be generalized. The cost gk is said to be equivariant if
there exists a family of diffeomorphisms ϕα such that
gk (φα (xk ), χα (uk ), ψα (wk )) = ϕα ◦gk (xk , uk , wk ). As we

j
Xk+1
= Xkj M (ujk )Wkj

5

(3)

These dynamics can be embedded in the threedimensional special Euclidean matrix Lie group
G = SE(2), by defining the state

where Xk , M (uk ), Wk are all square matrices belonging to G, uk is a control that lives in some finite dimensional vector space, and Wk is the noise. The control objective is to reach a desired configuration, that
is, a desired value for the relative configurations of the
agents (X 1 )−1 X 2 , . . . , (X K−1 )−1 X K , see e.g., Sarlette
et al. (2010) for more information.



cos θk − sin θk zk



Xk = 
 sin θk cos θk yk  ,
0
0
1

Systems of this form are naturally invariant to left multiplication of all X j by some matrix A ∈ G:


input matrix



AX 1


 . 
φA (X) =  .. 


AX K



cos( L1 vk tan sk ) − sin( L1 vk tan sk ) vk


1
1

M (vk , sk ) = 
 sin( L vk tan sk ) cos( L vk tan sk ) 0  ,
0
0
1

where X = (X 1 , . . . , X k ) ∈ G K . Letting χ(u1 , . . . , uK ) ≡
(u1 , . . . , uK ), ψ(W 1 , . . . , W K ) ≡ (W 1 , . . . , W K ), and
the costs be of the form g̃((X 1 )−1 X 2 , . . . , (X K−1 )−1 X K )+
h(u1 , . . . , uK ), we get an invariant system with invariant
costs.
One can define a cross section to the orbits by letting the
first agent coordinates be equal to the identity matrix,
that is, C = {X ∈ G K : X 1 = I}. The normalization
equation is given by I = φaγ(X) (X) = γ(X)X 1 , hence the
moving frame is given by γ(X) = (X 1 )−1 . The invariants
are computed as

0




0
,
1

We wish to compute a control policy for a two-vehicle
system, with states X 1 and X 2 , where the controls can
only take a finite number of values, and with terminal
cost
"
#

1
2
1 −1 2
J(X0 , X0 ) = E gN (XN ) XN


(X 1 )−1 X 2




..
=
.
.


(X 1 )−1 X K

p
2 + X 2 − 1|, that
where gN (X) = arccos(X11 )2 + | X13
23
is, we want the vehicles to have the same heading, and
follow each other at unit distance. Thanks to the theory
developed above, the stochastic control problem is reduced from problem with a six dimensional state space
to a problem with a three dimensional state space only.

The optimal stochastic control problem can then be
solved in the reduced coordinate system defined by ρ,
reducing the state space from dimension Kr to (K −1)r.
4


cos wk − sin wk


Wk =  sin wk cos wk
0
0

with state update equation of the form (3).



ρ(X) = φbγ(X)

and noise matrix

For numerical simulations, the cost functions J¯k were
computed on a fixed grid of dimension 51 × 51 × 65
using turning radius parameter L = 1, input sets vk ∈
{−0.1, 0, −0.1} and sk ∈ {−1, 0, −1} Globally optimal
input and state trajectory sequences corresponding to
iT
h
are
the initial condition x0 = 0.1 0 12 π −0.1 0 32 π
shown in Figures 1 and 2. These are compared against
a deterministic version of the model with wk = 0 in
Figures 3 and 4.

Application I: Cooperative formation control
for two stochastic Dubins vehicles

We consider two identical Dubins vehicles each with dynamics
zk+1 = zk + vk cos θk
yk+1 = yk + vk sin θk
1
θk+1 = θk + vk tan sk + wk
L

5

where yk and zk denote the two-dimensional position of
the vehicle, θk denotes the heading of the vehicle, vk is
a velocity input, sk is a steering angle input, and wk
is independent, identically-distributed zero-mean Gaussian noise with variance σ 2 , and L is a parameter that
determines the vehicle’s steering radius.

Application II: MRI Fingerprinting

Magnetic resonance imaging (MRI) has traditionally focused on acquisition sequences that are static, in the
sense that sequences typically wait for magnetization
to return to equilibrium between acquisitions. Recently,

6

���

���

���

���

�����

�
�
�������
�������

�����

�����

�����
���

����

�
�
�������
�������

���

����

����

����
�

�

��

��

��

��

�

�

��

����

��

��

��

����

Fig. 1. Optimal input sequence for cooperative stochastic
Dubins vehicle model with σ = 0.3.

Fig. 3. Optimal input sequence for cooperative deterministic
Dubins vehicle model.

���

���
���

���

�

���

�

���

����

����
����
����

���

���

���

���

�
����
����

Fig. 2. Optimal state sequence for the cooperative stochastic
Dubins vehicle model with σ = 0.3.

����

����

���

���

�

researchers have demonstrated promising results based
on dynamic acquisition sequences, in which spins are
continuously excited by a sequence of random input
pulses, without allowing the system to return to equilibrium between pulses. Model parameters corresponding
to T1 and T2 relaxation, off-resonance and spin density
are then estimated from the sequence of acquired data.
This technique, termed magnetic resonance fingerprinting (MRF), has been shown to increase the sensitivity,
specificity and speed of magnetic resonance studies (Ma
et al., 2013; Davies et al., 2014).

Fig. 4. Optimal state sequence for the cooperative deterministic Dubins vehicle model.

of T1 and T2 relaxation rates and the sequence of radiofrequency (RF) input pulses, used to excite the spins.

The following model was introduced in the conference
paper (Maidens et al., 2017). In this paper, an optimal
control was computed via dynamic programming on a
very sparse six-dimensional grid. Now using our symmetry reduction technique, we exploit symmetry reduction
to provide a much more accurate optimal input sequence
computed on a finer five-dimensional grid.

This technique could be further improved by replacing
randomized input pulse sequences with sequences that
have been optimized for informativeness about model
parameters. To this end, we present a model of MR spin
dynamics that describes the measured data as a function

7

in the joint distribution of Y = (y0 , . . . , yN ). The Fisher
information I can be expressed as a quadratic function
of the sensitivities of xk with respect to θ1 :

We model the spin dynamics via the equations

xk+1





0
θ2 0 0







= Uk 
 0 θ2 0  xk +  0 
1 − θ1
0 0 θ1

(4)


1/γ 0 0
 ∂
∂ T
I=
xk 
0 1/γ 0
 ∂θ1 xk

∂θ1
k=0
0 0 0
N
X

where the states x1,k and x2,k describe the transverse
magnetization (orthogonal to the applied magnetic field)
and x3,k describes the longitudinal magnetization (parallel to the applied magnetic field). To simplify the presentation, off-resonance is neglected in this model. Control inputs Uk ∈ SO(3) describe flip angles corresponding to RF excitation pulses that rotate the state about
the origin. Between acquisitions, transverse magnetization decays according to the parameter θ2 = e−∆t/T2 and
the longitudinal magnetization recovers to equilibrium
(normalized such that the equilibrium is x0 = [0 0 1]T )
according to the parameter θ1 = e−∆t/T1 where ∆t is
the sampling interval.

where the sensitivities
tivity equations:

∂
∂θ1 xk+1

We assume that data are acquired immediately following
the RF pulse, allowing us to make a noisy measurement
of the transverse magnetization. We also assume that the
measured data are described by a multivariate Gaussian
random variable
"
#
100
yk =
x k + vk
010

∂
∂θ1 xk

satisfy the following sensi-

 




0
θ2 0 0
0 00
 
 ∂



 



= Uk 
 0 θ2 0  ∂θ1 xk + Uk 0 0 0 xk +  0  .
−1
0 0 θ1
0 01

It should be noted that for system (4), the objective
function I has many local optima as a function of the input sequence Uk . Thus, in contrast with (Maidens et al.,
2016a) which consider optimal experiment design for hyperpolarized MRI problems, for this model, local search
methods provide little insight into what acquisition sequences are good. In contrast with the MRI model presented in (Maidens and Arcak, 2016), where global optimal experiment design heuristics are developed for linear
dynamical systems, in this model the decision variables
Uk multiply the state vector xk , making the output yk
a nonlinear function of the sequence U = (U0 , . . . Uk−1 ).
Thus we must use dynamic programming to find a solution.

where
# vk is a zero-mean Gaussian noise with covariance
"
γ 0
. This model results from a time discretization
0 γ
of the Bloch equations (Bloch, 1946; Nishimura, 2010)
under a time scale separation assumption that specifies
that the RF excitation pulses act on a much faster time
scale than the relaxation time constants T1 and T2 . A
simplified two-state version of this model was considered
in (Maidens et al., 2016b), where the transverse magnetization was modelled using a single state describing the
magnitude of [x1,k , x2,k ]T .

5.1

We see from the model (4) that magnetization in the
transverse direction decays while magnetization in the
longitudinal direction grows. However only the transverse component of the magnetization can be measured.
Thus there is a trade-off between making measurements
(which leads to loss of magnetization) and magnetization recovery. This is the trade-off that we hope to manage through the optimal design of an input sequence Uk .

Model

To present this problem in the formalism we have introduced, we define an augmented state vector

"
xk =

We wish to quantify the informativeness of an acquisition sequence based on the information about the T1 relaxation parameter θ1 that is contained in the resulting
data set. More formally, we wish to choose Uk ∈ SO(3)
to maximize the Fisher information about θ1 contained

xk
∂
∂θ1 xk

#
∈ R6 .

We can write the dynamics of the augmented state as
a control system with f and g defined in Equation (5).
This system has a one-dimensional group of symmetries

8


0 0 0 θ2 0


 Uk 0 0 0  0 θ2




0 0

0
0
0


fk (xk , Uk , wk ) = 

0 0
0 0 0


0 0 0 U   0 0
k 

0 0
0 0 0


00 0 0 0 0


0 0 0 0 0 0






0
0
0
0
0
0
T 
 xk
gk (xk , Uk , wk ) = −xk 

1
0 0 0 γ 0 0


0 0 0 0 1 0


γ
00 0 0 0 0


defined by

cos(α)

 sin(α)


 0
φα (xk ) = 

 0

 0

0

cos(α)

χα (Uk ) = 
 sin(α)
0

0

0

0






 0 
0 0 0 0






1 − θ1 
θ1 0 0 0 

 xk + 



0 θ2 0 0 
 0 



 0 
0 0 θ2 0 



−1
1 0 0 θ1

(5)

Next, we compute the invariants ρ(x) using
ρ(x) = φbγ(x)

sin(atan2(x1 , x2 )) cos(atan2(x1 , x2 ))


0
0


=
0
0


0
0


ψα (wk ) = wk

for any α ∈ R/2πZ.

5.2





where atan2 denotes the multi-valued inverse tangent
function

y
if x > 0,

arctan( x )

y


arctan(
)
+
π
if x < 0 and y ≥ 0,

x


arctan( xy ) − π if x < 0 and y < 0,
atan2(y, x) =
+ π2
if x = 0 and y > 0,



π

if x = 0 and y < 0,
−


 2
undefined
if x = 0 and y = 0.


0

cos(α) 0
0
0
0


0
1
0
0
0
 xk

0
0 cos(α) − sin(α) 0

0
0 sin(α) cos(α) 0

0
0
0
0
1



− sin(α) 0
cos(α) sin(α) 0





cos(α) 0
 Uk − sin(α) cos(α) 0
0
1
0
0
1
− sin(α) 0

0 0 0 0

0


p

x21 + x22

0

0
1
0
0
0


0

0
0
0


cos(atan2(x1 , x2 )) − sin(atan2(x1 , x2 )) 0 x

sin(atan2(x1 , x2 )) cos(atan2(x1 , x2 )) 0

0
0
1
0

0







x3


 1

√

=  x2 +x2 (x2 x4 − x1 x5 ) .
 1 2

 √ 1 (x x + x x )
 x2 +x2 1 4
2 5 
 1 2

x6

Dynamic programming in reduced coordinates

To perform dynamic programming in a reduced coordinate system, we begin by defining the cross-section
C = {x : x1 = 0, x2 > 0}, and computing the moving
frame γ(x). To do so, we solve

Further, ρ restricted to the cross-section C is injective
with inverse ρ̄−1 : R+ × R4 → C given by ρ̄−1 (x̄) =
h
iT
0 x̄1 x̄2 x̄3 x̄4 x̄5 . The theory above tells us we can
thus solve the optimal stochastic control problem in a 5
dimensional state space, reducing the original 6 dimensional problem of 1 dimension.

0 = φaγ(x) (x) = x1 cos γ(x) − x2 sin γ(x).
Isolating γ yields

Note that this reduction can be computed using only the
state transformation group φ, without reference to the

γ(x) = atan2(x1 , x2 )

9

state update equation f . Thus the same reduction can
be applied to any system with the same symmetries.
5.3

Results

To implement this algorithm, we discretize the reduced
five-dimensional state space and two-dimensional input
space via grids of size 6 × 10 × 15 × 15 × 15 and 16 ×
8 respectively. The code was written in the Julia language and parallelized to allow evaluation of Jk in parallel across grid points (Maidens et al., 2016b). The implementation is publicly available at https://github.
com/maidens/Automatica-2017.
Optimal input and state trajectories for the model corresponding to the initial condition at the equilibrium
x0 = [0 0 1 0 0 0]T are plotted in Figures 5 and 6.

(a) Magnetizations

(b) Sensitivities
Fig. 6. Optimal state sequence for the MR fingerprinting
model. Here we have plotted the longitudinal and transverse
components of both the magnetization (states x1 , x2 , and
x3 ) and the sensitivities (states x4 , x5 , and x6 ) where the
transverse component is computed as the Euclidean norm of
the vectors (x1 , x2 ) and (x4 , x5 ) respectively.

Fig. 5. Optimal input sequence for the MR fingerprinting
model. The angles α, β and δ represent rotations about the
z, y and x axes respectively, resulting in an control input
Uk = Rz (αk )Ry (βk )Rx (γk ).

In contrast with the results from (Maidens et al., 2016b)
where we considered a simplified version of the model,
for this full model we no longer find that the optimal flip
angle sequence converges to a cyclic pattern, rather it
appears irregular. However, state sequence of longitudinal magnetizations and transverse magnetization magnitudes appears to converge to a constant sequence. This
is likely because in this work we assumed Gaussian noise
in the inputs in contrast with the Rician noise assumed
in the previous work, therefore it is no longer necessary
to conserve magnetization across multiple time steps before generating a reliable measurement.
6

moderate state dimension. We have applied this technique to compute globally optimal trajectories to a sixdimensional original MRI model with a one-dimensional
group of symmetries and for a six-dimensional stochastic
Dubins vehicle model with a three-dimensional group of
symmetries by reducing the dimension of the state space
to five and three dimensions respectively. Since computation time for dynamic programming depends exponentially on the state space dimension, this technique
enables the computation of optimal control policies for
systems in which it was previously infeasible.

Conclusion

We have presented a method of reducing the complexity
of dynamic programming for systems in which the state
dynamics, stage costs and transition probabilities are invariant under a group of symmetries. This allows us to
compute globally optimal control policies for systems of

Acknowledgements
Research supported in part by the National Science
Foundation under grant ECCS-1405413.

10

volume 25, pages 183–198, 1998.
M. Kwiatkowska, G. Norman, and D. Parker. Symmetry
reduction for probabilistic model checking. In T. Ball
and R. B. Jones, editors, Computer Aided Verification: 18th International Conference, CAV 2006, Seattle, WA, USA, pages 234–248. Springer, 2006.
D. Ma, V. Gulani, N. Seiberlich, K. Liu, J. L. Sunshine,
J. L. Duerk, and M. A. Griswold. Magnetic resonance
fingerprinting. Nature, 495(7440):187–192, 03 2013.
J. Maidens and M. Arcak. Semidefinite relaxations in optimal experiment design with application to substrate
injection for hyperpolarized MRI. In Proceedings of
the American Control Conference (ACC), pages 2023
– 2028, 2016.
J. Maidens, J. W. Gordon, M. Arcak, and P. E. Z. Larson. Optimizing flip angles for metabolic rate estimation in hyperpolarized carbon-13 MRI. IEEE Transactions on Medical Imaging, 35:2403 – 2412, 2016a.
J. Maidens, A. Packard, and M. Arcak. Parallel dynamic programming for optimal experiment design in
nonlinear systems. In Proceedings of the Conference
on Decision and Control (CDC), pages 2894 – 2899,
2016b.
J. Maidens, A. Barrau, S. Bonnabel, and M. Arcak. Symmetry reduction for dynamic programming and application to MRI. In Proceedings of the American Control Conference (ACC), 2017. To appear.
J. Marsden, R. Montgomery, and T. Ratiu. Reduction,
Symmetry, and Phases in Mechanics. Number 436
in Memoirs of the American Mathematical Society.
American Mathematical Soc., 1990.
P. Martin, P. Rouchon, and J. Rudolph. Invariant tracking. ESAIM: Control, Optimisation and Calculus of
Variations, 10(1):1–13, 2004.
S. M. Narayanamurthy and B. Ravindran. Efficiently
exploiting symmetries in real time dynamic programming. In Proceedings of the 20th International Joint
Conference on Aritifical Intelligence, IJCAI’07, pages
2556–2561. Morgan Kaufmann, 2007.
D. G. Nishimura. Principles of Magnetic Resonance
Imaging. Lulu, 2010.
T. Ohsawa. Symmetry reduction of optimal control systems and principal connections. SIAM Journal on
Control and Optimization, 51(1):96–120, 2013.
P. J. Olver. Classical invariant theory, volume 44. Cambridge University Press, 1999.
W. B. Powell. Approximate Dynamic Programming:
Solving the Curses of Dimensionality.
WileyInterscience, 2007.
W. B. Powell. Perspectives of approximate dynamic programming. Annals of Operations Research, 241(1):
319–356, 2016.
A. Rahmani, M. Ji, M. Mesbahi, and M. Egerstedt.
Controllability of multi-agent systems from a graphtheoretic perspective. SIAM Journal on Control and
Optimization, 48(1):162–186, 2009.
W. Respondek and I. A. Tall. Nonlinearizable singleinput control systems do not admit stationary symmetries. Systems & Control Letters, 46(1):1–16, 2002.

References
F. Alvarez and N. L. Stokey. Dynamic programming
with homogeneous functions. Journal of economic theory, 82(1):167–189, 1998.
M. Arcak, C. Meissen, and A. Packard. Symmetry reduction. In Networks of Dissipative Systems: Compositional Certification of Stability, Performance, and
Safety, pages 55–62. Springer, 2016.
A. Barrau and S. Bonnabel. The invariant extended
Kalman filter as a stable observer. Arxiv preprint,
page arXiv:1410.1465, 2014.
R. Bellman. Dynamic Programming. Princeton University Press, 1957.
R. Bellman and S. Dreyfus. Functional approximations
and dynamic programming. Mathematical Tables and
Other Aids to Computation, 13(68):247–251, 1959.
D. P. Bertsekas. Dynamic Programming and Optimal
Control, Volume I. Athena Scientific, 3rd edition,
2005.
D. P. Bertsekas. Dynamic Programming and Optimal
Control, Volume II: Approximate Dynamic Programming. Athena Scientific, 4th edition, 2012.
A. M. Bloch, P. S. Krishnaprasad, J. E. Marsden, and
R. M. Murray. Nonholonomic mechanical systems
with symmetry. Arch. Rational Mech. Anal., 136:21–
99, 1996.
F. Bloch. Nuclear induction. Phys. Rev., 70:460–474,
1946.
G. Bluman and S. Kumei. Symmetries and Differential
Equations. Applied Mathematical Sciences. Springer,
2013.
S. Bonnabel, P. Martin, and P. Rouchon. Symmetrypreserving observers. IEEE Transactions on Automatic Control, 53(11):2514–2526, 2008.
F. Bullo and R. M. Murray. Tracking for fully actuated
mechanical systems: a geometric framework. Automatica, 35(1):17–34, 1999.
É. Cartan. La théorie des groupes finis et continus et
la géométrie différentielle: traitées par la méthode du
repère mobile. Cahiers scientifiques. Gauthier-Villars,
1937.
P. A. Clarksonz and E. L. Mansfield. Symmetry reductions and exact solutions of a class of nonlinear heat
equations. Physica D: Nonlinear Phenomena, 70(3):
250 – 288, 1994.
M. Davies, G. Puy, P. Vandergheynst, and Y. Wiaux. A
compressed sensing framework for magnetic resonance
fingerprinting. SIAM Journal on Imaging Sciences, 7
(4):2623–2656, 2014.
E. A. Emerson and A. P. Sistla. Symmetry and model
checking. Formal Methods in System Design, 9(1):
105–131, 1996.
J. Grizzle and S. Marcus. Optimal control of systems
possessing symmetries. IEEE Transactions on Automatic Control, 29(11):1037–1040, 1984.
B. Jakubczyk. Symmetries of nonlinear control systems
and their symbols. In Canadian Math. Conf. Proceed,

11

A. Rufino Ferreira, C. Meissen, M. Arcak, and
A. Packard. Symmetry reduction for performance certification of interconnected systems. IEEE Transactions on Control of Networked Systems, 2017.
A. Sarlette, S. Bonnabel, and R. Sepulchre. Coordinated
motion design on lie groups. IEEE Transactions on
Automatic Control, 55(5):1047–1058, 2010.
M. Zinkevich and T. Balch. Symmetry in Markov decision processes and its implications for single agent
and multi agent learning. In In Proceedings of the 18th
International Conference on Machine Learning, pages
632–640. Morgan Kaufmann, 2001.

12

