John Maidens

Axel Barrau

Silvère Bonnabel

Murat Arcak

Symmetry reduction for dynamic programming

Keywords: Optimal control, Dynamic programming, Invariant systems, Model reduction, Path planning, Medical applications

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 rely on explicitly transforming the state update equations and can therefore be applied in circumstances where this is difficult. We illustrate these results on two six-dimensional optimal control problems that are computationally difficult to solve by dynamic programming without symmetry reduction.

Introduction

The dynamic programming algorithm for computing optimal control policies has, since its development, been known to suffer from the "curse of dimensionality" [START_REF] Bellman | Dynamic Programming[END_REF]. 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 [START_REF] Bellman | Functional approximations and dynamic programming[END_REF][START_REF] Bertsekas | Dynamic Programming and Optimal Control, Volume II: Approximate Dynamic Programming[END_REF][START_REF] Powell | Approximate Dynamic Programming: Solving the Curses of Dimensionality[END_REF][START_REF] Powell | Perspectives of approximate dynamic programming[END_REF].

In practice, many real-world systems exhibit symmetries that can be exploited to reduce the complexity of 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 (Silvère Bonnabel), arcak@eecs.berkeley.edu (Murat Arcak). system models. Symmetry reduction has found applications in fields ranging from differential equations [START_REF] Clarksonz | Symmetry reductions and exact solutions of a class of nonlinear heat equations[END_REF][START_REF] Bluman | Symmetries and Differential Equations[END_REF] to model checking [START_REF] Emerson | Symmetry and model checking[END_REF][START_REF] Kwiatkowska | Symmetry reduction for probabilistic model checking[END_REF]. In control engineering, symmetries have been exploited to improve control of mechanical systems [START_REF] Bloch | Nonholonomic mechanical systems with symmetry[END_REF][START_REF] Bullo | Tracking for fully actuated mechanical systems: a geometric framework[END_REF], develop more reliable state estimators [START_REF] Barrau | The invariant extended Kalman filter as a stable observer[END_REF], study the controllability of multiagent systems [START_REF] Rahmani | Controllability of multi-agent systems from a graphtheoretic perspective[END_REF] and to reduce the complexity of stability and performance certification for interconnected systems [START_REF] Arcak | Symmetry reduction[END_REF][START_REF] Ferreira | Symmetry reduction for performance certification of interconnected systems[END_REF]. Symmetry reduction has also been applied to the computation of optimal control policies for continuoustime systems in [START_REF] Grizzle | Optimal control of systems possessing symmetries[END_REF][START_REF] Ohsawa | Symmetry reduction of optimal control systems and principal connections[END_REF] 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 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. Fur-ther, 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.

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. This paper builds on the work we presented in the conference paper [START_REF] Maidens | Symmetry reduction for dynamic programming and application to MRI[END_REF]. 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 5, along with examples to illustrate the algorithm in these contexts. 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 symmetry reduction to compute the solution of an optimal control problem arising in dynamic MRI acquisition. In Section 5 we apply the algorithm to the problem of coordination on Lie groups. Simulations are performed on a cooperative control problem for two Dubins vehicles. Code to reproduce the computational results in this paper is available at https://github.com/maidens/Automatica-2017.

Dynamic Programming and Symmetries

In this section, we first recall the main features of dynamic programming for optimal control of stochastic discrete 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.

Dynamic programming for optimal control of stochastic systems

We consider a discrete-time dynamical system

x k+1 = f k (x k , u k , w k ), k = 0, 1, . . . , N -1 (1)
where x k ∈ X ⊆ R n is the system state, u k ∈ U ⊆ R m is the control variable to be chosen at time k, w k ∈ W ⊆ R are independent continuous random variables each with density p k , and N ∈ Z + is a finite control horizon. Associated with this system is an additive cost function

g N (x N ) + N -1 k=0 g k (x k , u k , w k )
that we wish to minimize through our choice of u k . We define a control system to be a tuple S = (X , U, W, p, f, g, N ) where p = N -1 k=0 p k is the joint density of the random variables w k .

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 x 0 and a control policy π, we define the expected cost under this policy as

J π (x 0 ) = E g N (x N ) + N -1 k=0 g k (x k , µ k (x k ), w k ) .
An optimal policy π * is defined as one that minimizes the expected cost:

J π * (x 0 ) = min π∈Π J π (x 0 )
where Π denotes the set of all admissible control policies. The optimal cost function, denoted J * (x 0 ), is defined to be the expected cost corresponding to an optimal policy. An optimal policy π * and the optimal cost function J * can be computed using the dynamic programming algorithm. We quote the following result due to Bellman from [START_REF] Bertsekas | Dynamic Programming and Optimal Control[END_REF]:

Proposition 1 (Dynamic Programming) For every initial state x 0 , the optimal cost J * (x 0 ) of the basic problem is equal to J 0 (x 0 ), given by the last step of the following algorithm, which proceeds backward in time from period N -1 to period 0:

J N (x N ) = g N (x N ) J k (x k ) = min uk∈U E g k (x k , u k , w k ) + J k+1 f k (x k , u k , w k ) k = 0, 1, . . . , N -1, (2)
where the expectation is taken with respect to the probability distribution of w k . Furthermore, if u * k = µ * k (x k ) minimizes the right hand side of (2) for each x k and k, then the policy π * = {µ * 0 , . . . , µ * N -1 } is optimal.

Invariant system with invariant costs

We first recall the definition of a transformation group for a control system, as in [START_REF] Martin | Invariant tracking[END_REF][START_REF] Jakubczyk | Symmetries of nonlinear control systems and their symbols[END_REF][START_REF] Respondek | Nonlinearizable singleinput control systems do not admit stationary symmetries[END_REF]. See [START_REF] Olver | Classical invariant theory[END_REF] for the more general theory.

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.

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, x k ∈ X , u k ∈ U and w k ∈ W we have:

φ -1 • f k (φ(x k ), χ(u k ), ψ(w k )) = f k (x k , u k , w k ), k = 0, 1, . . . , N -1 g k (φ(x k ), χ(u k ), ψ(w k )) = g k (x k , u k , w k ), k = 0, 1, . . . , N -1, g N (φ(x N )) = g N (x N ), and p k (ψ(w k ))| det Dψ(w k )| = p k (w k ) k = 0, 1, . . . , N -1.
The rationale is simple: For any fixed α ∈ G, consider the change of variables

X k = φ α (x k ), U k = χ α (u k ), W k = ψ α (w k )
. Then, we have

X k+1 = f k (X k , U k , W k ), k = 0, 1, . . . , N -1,
and for k = 0, 1, . . . , N -1 we have also

g k (X k , U k , W k ) = g k (x k , u k , w k ).
As a result, if u 1 , . . . , u N -1 is a series of controls that minimize J(x 0 ), then one can expect U 1 , . . . , U N -1 to minimize J(X 0 ), 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 {φ α (x 0 )|α ∈ G}, reducing the initial n dimensional problem to a n -r dimensional problem. The present paper derives a proper theory for such symmetry reduction in dynamic programming, and provides various examples of engineering interest.

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 [START_REF] Cartan | La théorie des groupes finis et continus et la géométrie différentielle: traitées par la méthode du repère mobile[END_REF]. In general, this method only results in a local coordinate transformation as it relies on the implicit function theorem. However, for many practical problems the transformation computed using this method extends globally. To simplify the exposure we will present a "global" version of the method.

We briefly introduce the moving frame method following the presentation in [START_REF] Bonnabel | Symmetrypreserving observers[END_REF]. 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 : φ a e (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 be denoted γ(x), and the map γ : X → G will be called moving frame.

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

φ a γ(x) (x) = c.
Define the following map ρ : X → R n-r as

ρ(x) = φ b γ(x) (x).
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.

Main Results

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.

Symmetries imply equivalence classes of optimal policies

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 J k (x 0 ) satisfy the symmetry relations

J k = J k • φ α
for any k = 0, . . . , N and any α ∈ G. Furthermore, if

π * = {µ * 0 , . . . , µ * N -1 } is an optimal policy then so is π * := {χ α • µ * 0 • φ -1 α , . . . , χ α • µ * N -1 • φ -1 α } for any α ∈ G.
PROOF. First, note that

J N (x N ) = g N (x N ) = g N (φ(x N )) = J N (φ(x N )).
Now, suppose that for some k ∈ {0, . . . , N -1} we have

J k+1 (x k+1 ) = J k+1 (φ(x k+1
)) for all x k+1 ∈ X . Then for any x k ∈ X, and u k ∈ U we have

g k (φ(x k ), χ(u k ), wk ) + J k+1 (f k (φ(x k ), χ(u k ), wk )) p k ( wk )d wk = W g k (φ(x k ), χ(u k ), w k ) + J k+1 (f k (φ(x k ), χ(u k ), w k )) p k (w k )dw k = E g k (φ(x k ), χ(u k ), w k ) + J k+1 (f k (φ(x k ), χ(u k ), w k ))
Therefore,

J k (x k ) = min uk∈U E g k (x k , u k , w k ) + J k+1 (f k (x k , u k , w k )) = min uk∈U E g k (φ(x k ), χ(u k ), w k ) + J k+1 (f k (φ(x k ), χ(u k ), w k )) = min ũk∈χ(U ) E g k (φ(x k ), ũk , w k ) + J k+1 (f k (φ(x k ), ũk , w k )) = J k (φ(x k )).
Thus

J * = J * • φ. Now, if π * = {µ * 0 , . . . , µ * N -1
} is an optimal policy and we denote xk = φ(x k ) then for any k ∈ {0, . . . , N -1} we have

Jk(xk) = 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)) = E gk(φ(xk), χ • µ * k • φ -1 (φ(xk)), wk) + Jk+1(fk(φ(xk), χ • µ * k • φ -1 (φ(xk)), wk)) = E gk(xk, χ • µ * k • φ -1 (xk), wk) + Jk+1(fk(xk, χ • µ * k • φ -1 (xk), wk)) Thus π * := {χ • µ * 0 • φ -1 , . . . , χ • µ * N -1 • φ -1
} is an optimal policy. 2 3.2 Dynamic programming can be preformed 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 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 → R n-r .

For x ∈ ρ(C) ⊂ R n-r , let z ∈ C be such that x = ρ(z), and define Jk (x) = J k (z).

The following result shows that the functions J k on the n-dimensional space X ⊆ R n are completely determined by the values of Jk on the subset ρ(C) of R n-r .

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

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

where x := ρ(φ γ(x) (x)) is well defined as φ γ(x) (x) ∈ C.
It is thus sufficient to have evaluated J at all points of ρ(C) ⊂ R n-r to be able to instantly evaluate J at any point of X .

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

Jk (x) = min uk∈U E g k (z, u k , w k ) + Jk+1 (ρ(f k (z, u k , w k ))) .
PROOF. We have

Jk (x) = J k (z) = min uk∈U E g k (z, u k , w k ) + J k+1 f k (z, u k , w k ) = min uk∈U E g k (z, u k , w k ) + Jk+1 ρ • f k (z, u k , w k ) . 2

Case of equivariant costs

So far, we have considered the costs to be invariant. However, equivariance can obviously also be considered. Indeed, the cost g k is said to be equivariant if there exists a family of diffeomorphisms ϕ α such that

g k (φ α (x k ), χ α (u k ), ψ α (w k )) = ϕ α •g k (x k , u k , w k ).
As we 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

J k (φ(x)) = ϕ • J k (x)
as already noticed in [START_REF] Alvarez | Dynamic programming with homogeneous functions[END_REF] for the case of homogeneous costs. Symmetry reduction can then be applied. We now give two tutorial examples.

Example 1 Consider the linear system

x x+1 = Ax k + Bu k with quadratic costs g k = x T k Qx k + u T k Ru k .
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 equations yield γ(x) = 1/||x||. Applying the results above, we see that the controls that minimize

J(x 0 ) are ||x 0 ||u * 1 , • • • , ||x 0 ||u * N -1 , where u * 1 , • • • , u * N -1 are those minimizing J( x0 ||x0||
). 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

u k = -F k x k .
Example 2 Consider the following system and costs

x x+1 = Ax k + Bu k , g k = h(x k ) + ||u k || 1
where ||u k || 1 denotes the L 1 norm of u k 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 L 1 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.

Application to 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, 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 T 1 and T 2 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 [START_REF] Ma | Magnetic resonance fingerprinting[END_REF][START_REF] Davies | A compressed sensing framework for magnetic resonance fingerprinting[END_REF]. 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 of T 1 and T 2 relaxation rates and the sequence of radiofrequency (RF) input pulses, used to excite the spins.

The following model was introduced in the conference paper [START_REF] Maidens | Symmetry reduction for dynamic programming and application to MRI[END_REF]. 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 we exploit symmetry reduction to provide a much more accurate optimal input sequence computed on a finer five-dimensional grid.

We model the spin dynamics via the equations

x k+1 = U k     θ 2 0 0 0 θ 2 0 0 0 θ 1     x k +     0 0 1 -θ 1     (3)
where the states x 1,k and x 2,k describe the transverse magnetization (orthogonal to the applied magnetic field) and x 3,k describes the longitudinal magnetization (parallel to the applied magnetic field). To simplify the presentation, off-resonance is neglected in this model. Control inputs U k ∈ 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 x 0 = [0 0 1] T ) according to the parameter θ 1 = e -∆t/T1 where ∆t is the sampling interval.

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

f k (x k , U k , w k ) =             0 0 0 U k 0 0 0 0 0 0 0 0 0 0 0 0 U k 0 0 0                         θ 2 0 0 0 0 0 0 θ 2 0 0 0 0 0 0 θ 1 0 0 0 0 0 0 θ 2 0 0 0 0 0 0 θ 2 0 0 0 1 0 0 θ 1             x k +             0 0 1 -θ 1 0 0 -1             g k (x k , U k , w k ) = -x T k             0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 γ 0 0 0 0 0 0 1 γ 0 0 0 0 0 0 0             x k (4)
measured data are described by a multivariate Gaussian random variable

y k = 1 0 0 0 1 0 x k + v k
where v k is a zero-mean Gaussian noise with covariance γ 0 0 γ . This model results from a time discretization of the Bloch equations [START_REF] Bloch | Nuclear induction[END_REF][START_REF] Nishimura | Principles of Magnetic Resonance Imaging[END_REF] 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 T 1 and T 2 . 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 [x 1,k , x 2,k ] T .

We see from the model (3) 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 U k .

We wish to quantify the informativeness of an acquisition sequence based on the information about the T 1 relaxation parameter θ 1 that is contained in the resulting data set. More formally, we wish to choose U k ∈ SO(3) to maximize the Fisher information about θ 1 contained in the joint distribution of Y = (y 0 , . . . , y N ). The Fisher information I can be expressed as a quadratic function of the sensitivities of x k with respect to θ 1 :

I = N k=0 ∂ ∂θ 1 x T k     1/γ 0 0 0 1/γ 0 0 0 0     ∂ ∂θ 1 x k
where the sensitivities ∂ ∂θ1 x k satisfy the following sensitivity equations:

∂ ∂θ1 x k+1 = U k     θ 2 0 0 0 θ 2 0 0 0 θ 1     ∂ ∂θ1 x k + U k     0 0 0 0 0 0 0 0 1     x k +     0 0 -1     .
It should be noted that for system (3), the objective function I has many local optima as a function of the input sequence U k . 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 U k multiply the state vector x k , making the output y k a nonlinear function of the sequence U = (U 0 , . . . U k-1 ). Thus we must use dynamic programming to find a solution.

Model

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

x k = x k ∂ ∂θ1 x k ∈ R 6 .
We can write the dynamics of the augmented state as a control system with f and g defined in Equation (4).

This system has a one-dimensional group of symmetries defined by

φ α (x k ) =             cos(α) -sin(α) 0 0 0 0 sin(α) cos(α) 0 0 0 0 0 0 1 0 0 0 0 0 0 cos(α) -sin(α) 0 0 0 0 sin(α) cos(α) 0 0 0 0 0 0 1             x k χ α (U k ) =     cos(α) -sin(α) 0 sin(α) cos(α) 0 0 0 1     U k     cos(α) sin(α) 0 -sin(α) cos(α) 0 0 0 1     ψ α (w k ) = w k
for any α ∈ R/2πZ.

Dynamic programming in reduced coordinates

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

0 = φ a γ(x) (x) = x 1 cos γ(x) -x 2 sin γ(x).
Isolating γ yields

γ(x) = atan2(x 1 , x 2 ).
Next, we compute the invariants ρ(x) using

ρ(x) = φ b γ(x) =          sin(atan2(x1, x2)) cos(atan2(x1, x2)) 0 0 0 0 0 0 1 0 0 0 0 0 0 cos(atan2(x1, x2)) -sin(atan2(x1, x2)) 0 0 0 0 sin(atan2(x1, x2)) cos(atan2(x1, x2)) 0 0 0 0 0 0 1          x =            x 2 1 + x 2 2 x3 1 √ x 2 1 +x 2 2 (x2x4 -x1x5) 1 √ x 2 1 +x 2 2 (x1x4 + x2x5) x6           
Further, ρ restricted to the cross-section C is injective with inverse ρ-1 : R + × R 4 → C given by ρ-1 (x) = 0 x1 x2 x3 x4 x5

T . 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.

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 J k 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 x 0 = [0 0 1 0 0 0] T are plotted in Figures 1 and2. 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.

Optimal Formation Control on Lie Groups

We now apply the theory presented in Section 2.2 where the state space X ⊆ R n 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.

General problem

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

X j k+1 = X j k M (u j k )W j k (5)
where X k , M (u k ), W k are all square matrices belonging to G, u k is a control that lives in some finite dimensional vector space, and W k 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., [START_REF] Sarlette | Coordinated motion design on lie groups[END_REF] for more information.

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

φ A (X) =      AX 1 . . . AX K     
where X = (X 1 , . . . , X k ) ∈ G K . Letting χ(u 1 , . . . , u K ) ≡ (u 1 , . . . , u K ), ψ(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(u 1 , . . . , u K ), 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 are 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

ρ(X) = φ b γ(X) =      (X 1 ) -1 X 2 . . . (X 1 ) -1 X K     
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.

Application: cooperative formation control for two stochastic Dubins vehicles

We consider two identical Dubins vehicles each with dynamics

z k+1 = z k + v k cos θ k y k+1 = y k + v k sin θ k θ k+1 = θ k + 1 L v k tan s k + w k
where y k and z k denote the two-dimensional position of the vehicle, θ k denotes the heading of the vehicle, v k is a velocity input, s k is a steering angle input, and w k is independent, identically-distributed zero-mean Gaussian noise with variance σ 2 , and L is a parameter that determines the vehicle's steering radius.

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

X k =     cos θ k -sin θ k z k sin θ k cos θ k y k 0 0 1     , input matrix M (v k , s k ) =     cos( 1 L v k tan s k ) -sin( 1 L v k tan s k ) v k sin( 1 L v k tan s k ) cos( 1 L v k tan s k ) 0 0 0 1     ,
and noise matrix

W k =     cos w k -sin w k 0 sin w k cos w k 0 0 0 1     ,
with state update equation of the form (5).

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 J(X 1 0 , X 2 0 ) = E g N (X 1 N ) -1 X 2 N where g N (X) = arccos(X 11 ) 2 + | X 2 13 + X 2 23 -1|, that 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.

For numerical simulations, the cost functions Jk were computed on a fixed grid of dimension 51 × 51 × 65 using turning radius parameter L = 1, input sets v k ∈ {-0.1, 0, -0.1} and s k ∈ {-1, 0, -1} Globally optimal input and state trajectory sequences corresponding to the initial condition x 0 = 0.1 0 1 2 π -0.1 0 3 2 π T are shown in Figures 3 and4. These are compared against a deterministic version of the model with w k = 0 in Figures 5 and6.

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 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.

E

  gk(xk, uk, wk) + Jk+1(fk(xk, uk, wk)) = W gk(xk, uk, wk) + Jk+1(fk(xk, uk, wk)) pk(wk)dwk = W gk(φ(xk), χ(uk), ψ(wk)) + Jk+1(φ -1 • fk(φ(xk), χ(uk), ψ(wk))) pk(wk)dwk = W gk(φ(xk), χ(uk), ψ(wk)) + Jk+1(fk(φ(xk), χ(uk), ψ(wk))) pk(wk)dwk = W gk(φ(xk), χ(uk), ψ(wk)) + Jk+1(fk(φ(xk), χ(uk), ψ(wk))) pk(ψ(wk))| det Dψ(xk)|dwk = ψ(W)

Fig. 1 .

 1 Fig. 1. 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 U k = Rz(α k )Ry(β k )Rx(γ k ).

  Fig. 2. 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. 4 .

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

Fig. 6 .

 6 Fig. 5. Optimal input sequence for cooperative deterministic Dubins vehicle model.

Acknowledgements

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