A Review of Geometric Optimal Control for Quantum Systems in Nuclear Magnetic Resonance

We present a geometric framework to analyze optimal control problems of uncoupled spin 1/2 particles occurring in nuclear magnetic resonance. According to the Pontryagin’s maximum principle, the optimal trajectories are solutions of a pseudo-Hamiltonian system. This computation is completed by sufficient optimality conditions based on the concept of conjugate points related to Lagrangian singularities. This approach is applied to analyze two relevant optimal control issues in NMR: the saturation control problem, that is, the problem of steering in minimum time a single spin 1/2 particle from the equilibrium point to the zero magnetization vector, and the contrast imaging problem. The analysis is completed by numerical computations and experimental results.


Introduction
The dynamics of a spin 1/2 particle in nuclear magnetic resonance NMR is described by the Bloch equation where M is the magnetization vector of the particle, ω is the radio-frequency control magnetic field, and R is the relaxation matrix 1, 2 .Such a system is an example of bilinear system in control theory 3 where H x, p, u p, F x, u , ẋ F x, u being the control system and p the adjoint vector.The function H depending upon u is called the pseudo-Hamiltonian and the optimal control has to satisfy for almost every time the maximization condition 4 The solutions of 1.3 and 1.4 are called extremals, and they can be classified into two types: i singular extremals if they satisfy the relation ∂H/∂u 0, ii regular extremals if the control takes its values in the boundary of the control domain.
General extremals are the concatenation of singular and regular arcs.An important example in our study is the case of single-input control systems ẋ F x uG x , where the control domain is |u| ≤ 1.A regular extremal is formed by bang arcs defined by u t sign p t , G x t , while the singular extremals satisfy p t , G x t 0. The optimal control problem reduces to find the sequence BSBS . .., where B is a bang arc with u ±1 and S a singular arc.For linear systems ẋ Ax Bu, under the mild Kalman controllability condition rank B, . . ., A n−1 B n, there exists no singular arc, and optimal solutions satisfy a bang-bang principle 5 .The situation is quite different in a bilinear system, which corresponds to a kind of universal nonlinear model, since in the analytic case, the inputoutput of any system can be approximated by the one of a bilinear system, see 6 .One objective of this paper is to show the ubiquity of the singular arcs in the optimal control of spin 1/2 particles.This point can already be detected in the example of a single spin 7 .The corresponding time-optimal control problem will be discussed in details in this paper since the discussion is surprisingly very intricate.Also it will serve as an introduction to the contrast problem in NMR imaging.
Except in some cases, for example, time minimal control for linear systems, the maximum principle is only a necessary optimality condition, and sufficient conditions are related to the concept of extremal field and Hamilton-Jacobi-Bellman equation.This leads to the introduction in optimal control of the difficult notion of conjugate points.Based on recent works 8, 9 , a review of this concept is presented in this paper for singular extremals.The important part consists in clarifying the relation between the geometric concept associated to Lagrangian singularities of the projection of the extremal flow on the state space and the problem of optimality which is characterized by the openness properties in a suitable topology of the input-output mapping where x • is the response at time t f of the system to the control u • , and x 0 is the initial state.Similar results exist in the bang-bang case, but they will not be used in the present paper, see 10 .In two final sections, the different concepts are used to analyze the two following problems.
i Time Minimum Problem for a Single Spin.Under suitable coordinates, the system takes the form where the state variable q x, y, z belongs to the Bloch ball |q| ≤ 1 which is invariant for the dynamics since the parameters belong to the set Λ : 2Γ ≥ γ ≥ 0. The control field is u u 1 iu 2 |u|e iϕ where the amplitude is bounded by a given m.A complete derivation of 1.6 will be done in Section 5.
ii Contrast Imaging Problem.In this case, a simplified model is formed by coupling two systems 1.6 with different relaxation parameters: written shortly as dx/dt F x, u where x q 1 , q 2 belongs to the product of two Bloch balls.Both spins interact through the same magnetic field represented by u • .The associated optimal control problem is the following.Starting from the equilibrium point x 0 0, 0, 1 , 0, 0, 1 , the goal is to reach in a given transfer time t f the final state q 1 t f 0 corresponding to zero magnetization of the first spin while maximizing |q 2 t f | 2 11 .

Necessary Optimality Conditions: Maximum Principle
We consider a control system of the form where The class of admissible control U is the set {u : 0, t f → U} of bounded and measurable maps where U ⊂ R m is the control domain.Given an initial condition x 0 ∈ R n and u • ∈ U, we denote by x t, x 0 , u the solution of the control system initiating from x 0 and defined on a subinterval of 0, t f .We fix a C ∞ terminal manifold M defined by g x 0 where g : R n → R k .We will consider the two following optimal control problems.It is well known that the Pontryagin maximum principle is a set of necessary conditions for our optimal control problem 3, 4 .

The Pontryagin Maximum Principle
Let u * • be an admissible control whose corresponding trajectory x * • x •, u * , x 0 is optimal, then there exists a vector function p * • ∈ R n \ {0} such that the following conditions are satisfied: with the maximality condition where M is a constant which can be chosen positive in the time minimal case.The following boundary conditions are satisfied for the two respective problems: i time minimum case: at the final time t f , ii Mayer problem: at the final time t f ,

2.6
The boundary conditions on the adjoint vector are called transversality conditions.

Geometric Interpretation
Both optimal problems satisfy the same Hamiltonian dynamics defined by 2.3 , 2.4 and the triples x • , p • , u • solutions of this system are called extremals.Geometrically they correspond to necessary conditions for the terminal point x t f to be in the boundary of the accessibility set A x 0 , t f .The respective boundary conditions define the so-called BC-extremals: i time minimum case: ii Mayer problem: let us fix m ∈ R and introduce the manifold M m {x, g x 0, C x m}.One must have for the optimal solution p * t f ⊥ M m , where m is such that the cost is minimum.

Singular Extremals and the Weak Maximum Principle
The weak maximum principle consists in replacing the maximality condition 2.4 by ∂H/∂u x, p, u 0. The following results are well known.
Definition 2.1.Relaxing u ∈ U, a singular trajectory of the system ẋ F x, u on 0, t f is a control trajectory pair x, u such that the input-output mapping is singular.
Proposition 2.2.One then deduces the following results.
1 If x, u is not singular on 0, t f , then the input-output mapping is open at u • for the L ∞ 0, t f -norm topology.Next, we make computations of extremals which will be used in the sequel.

Bi-Input Affine Case
A system such that the control enters linearly is called affine.We consider a situation for which F x, u

Single-Input Affine Case
The ii Singular extremals: since the system is linear in u, the condition 2.4 leads in the singular case to the condition H G 0, identically.Differentiating twice with respect to time, we get the relations and from this second condition, one derives when the denominator is not vanishing the corresponding singular control Plugging such u S into the pseudo-Hamiltonian defines the singular flow solution of the Hamiltonian vector field denoted as H S and starting at t 0 from the two constraints H G {H G , H F } 0 and |u S | ≤ 1.They have the following interpretation.

Proposition 2.4. Singular extremals correspond to the singularities of the input-output mapping
2.13

Higher-Order Maximum Principle and the Generalized Legendre-Clebsch Condition
From the maximality condition, one has the necessary optimality condition called the Legendre-Clebsch condition

The Relation between the Bi-Input and the Single Input Cases
An important relation in our analysis comes from the work of 9 .For an extremal of order zero, one has u 2 1 u 2 2 1, and one can consider the extension of the control system by setting u 1 sin α, u 2 cos α, and α v where x, α is the extended state space and v is the new control.
For the single-input case, one can define a reduction of the system, using the so-called Goh transformation.Assuming that G is nonzero, there exists a coordinate system x 1 , . . ., x n on an open set V such that G ∂/∂x n .Denoting x x , x n , the system splits into ẋ F x , x n , ẋn F 0 x u, 2.16 where the system F defined on an open set V with x n , taken as the control variable, is called the reduced system.For both cases, the extremals are in correspondence, due to the intrinsic interpretation of the singular trajectories as singularities of the input-output mapping.Moreover, introducing the reduced Hamiltonian H x , p , x n p , F x , x n , one has

2.17
These formulas establish the relation between the Legendre-Clebsch and the generalized Legendre-Clebsch of the corresponding systems.

Second-Order Optimality Condition
We will present the results needed in our analysis to get sufficient second-order optimality conditions under generic assumptions in the singular case.The smooth system is ẋ F x, u , u ∈ U, and x ∈ X, and the pseudo-Hamiltonian can be written as H z, u p, F x, u , where z x, p .One can consider a reference singular extremal z, u solution on 0, Our first assumption is the strong-Legendre condition: A1 The quadratic form ∂ 2 H/∂u 2 is negative definite along the reference solution.

Advances in Mathematical Physics
Using the implicit function theorem, the extremal control is then locally defined as a smooth function u z , and plugging this function into H defines a smooth true Hamiltonian, still denoted as H, and the extremal is solution of ż H z with initial condition z 0 x 0 , p 0 .

The Geometric Concept of Conjugate Point
Definition 3.1.Let z x, p be the reference extremal defined on 0, t f .The variational equation is called the Jacobi equation.A Jacobi field is a nontrivial solution δz δx, δp , and it is said to be vertical at time t if δx t dΠ z t δz t 0 where Π : x, p → x is the canonical projection.
The following geometric result is crucial 13 .Proposition 3.2.Let L 0 be the fiber T * x 0 X, and let L t exp t H L 0 be its image by the oneparameter subgroup generated by H, then L t is a Lagrangian manifold whose tangent space at z t is spanned by the Jacobi fields vertical at t 0.Moreover, the rank of the restriction to This leads to the following definition.Definition 3.3.We define the exponential mapping exp x 0 ,t p 0 Π z t, x 0 , p 0 3.3 as the projection on the state space of the integral curve of H with initial condition z 0 x 0 , p 0 , where p 0 can be restricted by linearity to the sphere |p 0 | 1.If z x, p is the reference extremal, a time t c > 0 is said to be geometrically conjugate to 0 if the mapping p 0 → exp x 0 ,t is not of rank n − 1 at t t c , and the associated point x t c is said to be geometrically conjugate to x 0 .We denote by t 1c the first conjugate time and C x 0 is the conjugate locus formed by the set of first conjugate points.
An algorithm can be deduced.

Testing Conjugacy
Let z t x t , p t be the reference extremal, and consider the vector space of dimension n − 1 generated by the Jacobi fields δz i δx i , δp i , i 1, . . ., n − 1 vertical at t 0 and such that δp i 0 is orthogonal to p 0 .At a conjugate time t c , one has rank δx 1 t c , . . ., δx n−1 t c < n − 1.

3.4
Seminal works in optimal control are to relate this concept to the optimality properties 8, 13 .
In relation with Mayer problem, this question is translated into openness of the input-output mapping.One needs to introduce the ad hoc generic assumptions completing A1 .
A2 The extremal trajectory z t , t ∈ 0, t f associated to u • is by definition a singularity of the input-output mapping E x 0 ,t f .One assumes that on each subinterval t 0 , t 1 , 0 < t 0 < t 1 ≤ t f , the singularity of E on u| t 0 ,t 1 is of codimension one.
A3 We are in the nonexceptional case where H / 0 along the reference extremal.Combined with the CotCot code 14 , this gives the practical point of view used to test optimality in the time minimal case and the contrast problem along an extremal of order zero, since nonopenness of the input-output mapping along any such arc is a necessary optimality condition.
A more complicate and subtle question is to consider the same problem along singular arcs for the single-input affine control system ẋ F x uG x , |u| ≤ 1.The corresponding theoretical results are presented in 9 based on the relation between the bi and single input cases using Goh transformation explained above.The conclusion is that, again in this case, the conjugate point analysis can be applied to test optimality in such situations.

Focal Type Conditions
The concept of conjugate point is associated to optimal control problems with fixed endpoints conditions.In the contrast problem, it has to be adapted to the problem with fixed initial condition x 0 x 0 , but with variable final condition x t f ∈ M m .From the maximum principle, the transversality condition defined a Lagrangian manifold z t f x t f , p t f ∈ M ⊥ m which again defines along a singular extremal a train of Lagrangian manifolds L t , integrating backwards t ≤ t f .Optimality can be deduced from the geometric concept of focal points corresponding to the singularities of L t , Π .

Extremal Fields and Hamilton-Jacobi-Bellman Equation
For simplicity, we consider here the time-minimal control problem.We pick up a reference extremal trajectory t → x t , t ∈ 0, t f such that A1 , A2 , and A3 are satisfied.Moreover, we assume that the reference trajectory is one-to-one and that there exists no conjugate point on 0, t f , then we can embed locally the reference extremal into a central field F formed by all extremal trajectories starting from x 0 x 0 .One restricts the adjoint vector with the normalization H 1 and the assumption A3 .One introduces where L t is the train of Lagrangian manifolds generated by the Hamiltonian flow H, with L 0 T * x 0 X being the initial condition.By construction, F is the projection of L on the state space X.Proposition 3.6.Excluding x 0 , there exists an open neighborhood W of the reference trajectory and two smooth mappings V : W → R, u : W → U such that for each x, u ∈ W × U, the maximization condition H x, dV x , u x ≥ H x, dV x , u .

3.6
The reference trajectory is optimal with respect to all smooth trajectories of the system, with the same extremities and contained in W. The Lagrangian manifold L is a graph, and V is the generating mapping {x ∈ W, p ∂V/∂x}.
In the next sections, our geometric tools will be applied to our two case studies.

The Single Spin 1/2 Case
One considers the time minimal control of steering the state 0, 0, 1 to the center of the Bloch ball.Due to the symmetry of revolution, one can restrict the problem to the 2D-single input system q F uG, Despite the apparent simplicity of the equations, the problem is very rich, and most of the geometric tools of the 2D-optimal control theory are needed to analyze the problem.One will use 13, 15 as general references for the concepts and techniques.To handle the problem, one introduces two steps: 1 given a point q 0 y, z , find the time minimal solution in a small neighborhood V of q 0 , according to the Lie algebraic structure of the system at q 0 ; 2 glue together the local solutions to get the global time optimal solution.More precisely, we will describe the global synthesis to steer N 0, 0, 1 to any point of the Bloch ball.By symmetry, one can only consider the domain y ≤ 0. This amounts to compute the switching locus Σ formed by points where the optimal solutions switch and the separating locus L where there exist two different optimal solutions which intersect.

Computations of Lie Brackets and Singular Trajectories Properties
One has and the Lie brackets up to order three are

4.4
The singular trajectories are contained in the set S defined by det G, G, F 0 which leads to y −2δz γ 0. Hence, the singular locus is the union of the two singular lines: the line y 0 corresponding to the axis of revolution and the horizontal axis z 0 γ/ 2δ .Eliminating p, the singular control is given by D

Parameter Conditions
The interesting case is when δ < 0 and |γ/2δ| < 1.This leads to In this case, u S → ∞ when y → 0 − .Moreover, we have the following estimate.Setting α γ 2 2Γ − γ / 4γ 2 , we get near 0 that ẏ −α/y, and the trajectory reaching zero in y < 0 is of order O √ t , the singular control being of order O 1/ √ t .In particular, one deduces the following.
Lemma 4.1.The singular control in the neighborhood of the point y 0, along the horizontal direction is of order O 1/ √ t and is in the L 1 but not L 2 category.

Optimality Problem
For small time, the optimality status of the singular arc can be deduced from the generalized Legendre-Clebsch conditions.In order to be optimal, we have H ≥ 0 from the maximum principle, and from the Legendre-Clebsch condition 2.15 , we deduce that {H G , {H G , H F }} ≤ 0.More precisely, if we introduce D det G, F γz z − 1 Γy 2 , one arrives at i hyperbolic case DD > 0, the singular direction is small time minimizing; ii elliptic case DD < 0, the singular trajectory is small time maximizing.
Applying this test when Γ > 3γ/2, one has i the horizontal singular line is a fast direction; ii the vertical singular line is fast provided z 0 < z < 1.
In particular, this test excludes the standard policy in NMR called the inversion recovery sequence 16, 17 .In this case, the spin is first steered close to the south pole 0, −1 with a bang pulse, a zero control is then applied to reach the target along the horizontal axis, using in particular the −1 ≤ z < z 0 time maximizing singular arc.
In order to complete the optimality analysis, one introduces for 2D-system the following clock-form ω.The collinear set C is defined by det F, G 0. From a straightforward computation, one deduces that it will form an oval curve defined by γz z − 1 Γy 2 0, which shrinks into a point if γ 0. Outside this set, ω pdx is defined by ω F 1 and ω G 0, the sign of dω being given by y γ − 2δz .This form allows to compare the time to travel different trajectories not crossing C by using the Stokes theorem.We have also dω 0 on the singular set S. Using the form ω, one deduces the following lemma.

Lemma 4.2. Assume m
∞, then the small broken arc formed by the concatenation of small pieces of horizontal and vertical singular lines is time optimal.
The analysis is completed by the optimality properties of bang-bang extremals.The two main tools are the classification of the regular extremals near the switching surface Σ {w q, p / p, G q 0}, see, for instance, 18 , and the concept of conjugate points for bang-bang extremals which is described in 19, 20 .

Classification of Regular Extremals Near Σ
The singular extremals are contained in the subset H G {H G , H F } 0, and we denote by δ and δ − bang arcs u ±1 for the system dx/dt F uG, |u| ≤ 1.We denote by δ S a singular arc and by δ 1 δ 2 an arc δ 1 followed by an arc δ 2 .We have the following.
Ordinary switching arc: it is a time t such that a bang arc switches with the condition Φ t 0 and Φ t {H G , H F } w t / 0 where Φ t H G w t is the switching function evaluated along an extremal arc w • .According to the maximum principle, near Σ , the extremal is of the form δ δ − resp., the second-order derivative along a bang arc ±1, the fold case is the situation where Φ ± t Φ± t 0 and Φ± / 0. We have three cases.i Hyperbolic Case. at the contact point with Σ , Φ t > 0 and Φ− t < 0. The connection with a singular extremal with |u S | < 1 and satisfying the strong generalized condition is possible.The extremals near such a point are of the form δ ± δ S δ ± .
ii Elliptic Case. at the contact point with Σ , Φ t < 0 and Φ− t > 0. The connection with a singular extremal is not possible, and locally every extremal is bang-bang but Advances in Mathematical Physics 13 with no uniform bound on the number of switchings.From the theory of conjugate points in the 2D-bang-bang case, every time optimal solution is locally of the form δ δ − or δ − δ .
iii Parabolic Case.Here, Φ t and Φ− t have the same sign at the contact point.One can check that the singular extremal is not admissible, and every extremal curve is locally bang-bang with at most two switchings: This classification is far from being sufficient to analyze locally our problems.In particular, we have the following situations.

Saturating Singular Case
It is a transition between the hyperbolic and the parabolic cases.The singular control u S saturates at a point B. This situation was analyzed in the literature, see, for instance, 15 .At the point B, there is a birth of a switching curve which can be estimated using the techniques of local models 13 .B is identified to 0, the singular arc is normalized to the x-axis, and the model can be written as The singular control is u S 1 x, and it is not admissible for x > 0. Near 0, we have optimal policies of the form δ δ S δ δ − .Even more complicated situation can occur, and the most interesting is the one related to the interaction between the two singular lines, which is described next.

The SiSi Singularity
Assume that Γ > 3γ/2 and the control bound m is large enough such that the bang arc u m starting from the north pole intersects the singular arc z 0 γ/ 2δ , Then for the problem of reaching zero magnetization, the horizontal singular arc is saturating at a point denoted as B and from the analysis of 7 ; they are local optimal policies of the form δ m δ Sh δ m δ Sv where δ m is an arc associated to u m, and δ Sh and δ Sv are, respectively, singular horizontal and vertical arcs.Note in particular that the horizontal singular arc is not followed up to the saturating point B. Hence, we introduce the following definition.Definition 4.3.One calls a bridge between the horizontal and the vertical singular arcs the bang arc such that the concatenation singular-bang-singular arc is optimal.
In our study, the existence of a bridge is related to the following phenomenon.At the saturating point B, there is a birth of a switching locus which meets the horizontal singular arc.
Next, we present a geometric method to evaluate switching points for 2D-systems.This method is effective if the system is bilinear.The geometric process is described in details in 15 .

Effective Evaluation of Switching Points
Instead of using the adjoint equation to determine the switching sequences, we introduce the following coordinate invariant point of view.Assume 0, t to be two consecutive switching times on an arc δ or δ − where the control is u ε ±1.We must have p 0 , G q 0 p t , G q t 0. 4.9 We denote by V • the solution of the variational equation V ∂F/∂q ε ∂G/∂q V such that V t G q t , where the equation is integrated backwards from time t to 0. By construction p 0 , V 0 0 since p t , V t 0 and the derivative is zero and hence at time t 0, p 0 is orthogonal to G q 0 and to V 0 .Therefore, V 0 and G q 0 are collinear.We denote by θ t the angle between G q 0 and V t measured counterclockwise.One deduces that switching occurs when θ t 0 mod π .Moreover, θ 0 on the singular set det G, G, F 0. We have by definition V t e −tad F εG G q t , 4.10 and in the analytic case, the ad-formula gives, for small t,

4.11
The computation can be made explicit in the bilinear case F q Aq a, G q Bq.The system can be lifted into an invariant system on the semidirect product Lie group GL 2, R × S R 2 identified with the set of matrices of GL 3, R : acting on the subspace of vector in R 3 : 1 q .Lie brackets of two affine vector fields A , a , B , b being defined by the computation of the exponential exp −t ad F εG is a standard exercise in linear algebra which amounts to compute a Jordan normal form of ad F εG , see 21 for the details of the computations.

Global Synthesis
In order to complete the global time minimal synthesis with initial point at the north pole, we must glue together the local syntheses, and the final result is represented in Figure 1.We have assumed that Γ > 3γ/2 and m Γ.Using the symmetry of revolution with respect to the z-axis, one can restrict the study to the domain y ≤ 0. The north pole being a singular point for the optimal control, the first applied control is taken as u m and will form a boundary arc of the accessibility domain.The collinear set is a parabola connecting O to the north pole.
The main feature of the synthesis is the SiSi singularity.The switching locus Σ is composed of the arc denoted Σ 1 starting from the north pole and reaching the horizontal singular arc at A, the horizontal line segment Σ 2 between A and B, the singular control saturating at B, the switching locus Σ 3 due to the saturating phenomenon ending on the vertical axis at D, and the vertical singular line segment denoted as Σ 4 between D and O.The bang arc u −m starting from A separates the synthesis in two domains, one with a bang-bang policy δ m δ −m and the other where the optimal policy contains a nontrivial horizontal singular line.Due to the symmetry of revolution, the separating locus L reduces to the z-axis.
Note that different works have shown that such optimal control field can be implemented with a very good accuracy in NMR experiments, the standard error being of the order of few percents 7, 22 .This point is illustrated in Figure 2 for the saturation control problem where a reasonable match between theory and experiments has been obtained.This result confirms that the optimized pulse sequence can really be implemented with modern NMR spectrometers.

The Contrast Problem
For simplicity, we will limit our analysis to the single-input case.The system can be written shortly as dx dt F x uG x , 5.1  In the upper panel, the small insert represents a zoom of the optimal trajectory near the origin.The dotted line is the switching curve originating from the horizontal singular line.The solid curve is the optimal trajectory near the origin.The lower panel displays the corresponding control laws.τ is the reduced time defined by τ ω max / 2π t.
where x q 1 , q 2 , q i y i , z i representing the state of each spin 1/2 particle.From Bloch equation, we get where Λ i Γ i , γ i are the relaxation parameters characterizing each spin.

BC Singular Extremals
According to the general computations of Section 2, a singular extremal is contained in the constrained set while the singular control is defined by the Lie brackets being computed in Section 2.
In order to be optimal, the singular extremals have to satisfy the generalized Legendre-Clebsch condition p, G, G, F ≤ 0 and the transversality condition at the final time.Decomposing the adjoint vector in the form p p 1 , p 2 where p i is associated to the spin i, one has and if p 0 / 0, it can be normalized by homogeneity to p 0 −1/2.The case p 0 0 has the following straightforward interpretation.
Lemma 5.1.The time-optimal solution of the first spin system is solution of the contrast problem provided the transfer time t f is fixed to the optimal time.

Exceptional Singular Extremals
If the transfer time is not fixed, then according to the maximum principle, this leads to the additional constraint M max u • H F uH G 0, which gives p, F x 0 in the singular case.With this restriction, the adjoint vector can be eliminated, and the singular control in this exceptional case is given by where with the corresponding vector fields defined by Observe that in the general case, a similar computation shows that the singular extremals are solutions of an equation of the form where λ is a one-dimensional time-dependent parameter whose dynamics is deduced from the adjoint equation.

Desingularization
Meromorphic differential equations of the form 5.8 can be desingularized as using the time reparameterization ds dt/D x t .According to the formula 5.4 , the singular control can explode near D 0, in particular saturating the constraint |u| ≤ m.This is connected to the bridge phenomenon described in Section 3. Some of the singular extremals can cross the singular surface D 0, provided D 0.

Feedback Classification and the Contrast Problem
An important object to analyze our problem is the action of the feedback group G f on the set of systems.One will restrict our presentation to the single-input affine case, and we denote by C F, G the set of all such smooth systems on a state space X R n .
Definition 5.2.Let F, G , F , G be two elements of C. They are called feedback equivalent if there exists a smooth diffeomorphism ϕ of R n on a feedback u α x uβ x where α, β are smooth, where ϕ * Z is the vector field image given by ϕ * Z ∂ϕ −1 /∂x Z • ϕ .This action defines a group structure on the set of triplets ϕ, α, β ; this group is denoted as G.
Definition 5.3.Let F, G ∈ C and let λ be the map which associates to F, G the constrained differential equation whose solutions are singular extremals.The constraint is given by Σ : H G {H G , H F } 0, and the equation is We define the action of ϕ, α, β on Σ, F S by the action of the symplectic change of coordinates ϕ : x ϕ X , p P ∂ϕ −1 ∂X , 5.12 a feedback acting trivially.
We have the following theorem 23 .

5.13
In other words, λ is a covariant.
Moreover, under mild assumption, λ is a complete covariant.In particular, an important geometric problem is to relate the set of parameters Γ 1 , γ 1 , Γ 2 , γ 2 to geometric invariants of the singular flow and optimality properties.The research program is the following.

Classification of the Distribution G
A feedback u α x β x v acts on the control direction G as G → βG, which corresponds to the classification of the one-dimensional distribution x → Span G x .

Collinear Set
The collinear set C is the set of points where F and G are collinear and is feedback invariant.Geometrically, it will form a one-dimensional curve.

The Singular Trajectories
It is defined by the constrained Hamiltonian equations given by

restricted to the surface
where u S is given by the formulae 2.14 .This set of equations defines a Hamiltonian vector field on Σ , restricting the standard symplectic form ω dx ∧ dp.
We introduce the two Poisson brackets: and the differential equation 5.14 can be desingularized using the time reparameterization ds dt/ D x t , p t .One gets the equations 5.17 which restricted to Σ are a semicovariant.

Advances in Mathematical Physics
Geometric invariants are related to the surface {D D 0}∩Σ and to the singularities of 5.17 in this surface.The geometric framework to analyze such nonisolated singularities is presented in 24, 25 .

Geometric and Numerical Methods
The object of this section is to combine geometric and numerical methods to analyze the contrast problem.The first step is to construct along a reference singular trajectory a C 0synthesis in the limit case m Γ, γ.This result is based on 9, 26 .

Synthesis for a Reference Singular Arc
We consider the control system dx/dt F x uG x , |u| ≤ 1, and let u S be a smooth singular control such that |u S | < 1 with corresponding trajectory x t .We have the following.Theorem 6.1.Under generic assumptions, the first conjugate time t 1c is the time along the singular arc such that a policy δ ± δ S δ ± is time extremizing in a C 0 -neighborhood of the reference singular arc.
If we apply this result to the contrast problem, one can deduce the simplest result about an extremal policy which provides a local optimal solution.The initial point is the north pole of the Bloch ball, and such a point is singular for the singular flow.Hence, the first control to be applied is u 1 or u −1 the maximum bound m being normalized to 1 , and by symmetry of revolution, one can assume that u 1.At the final time t f , one must have p 2 t f q 2 t f , and moreover H G t f 0. The final point concentrates a switching.Hence, one gets the following.Lemma 6.2.The simplest possible sequence is of the form BS, where B is a bang arc and S is a singular arc.
More complicated sequences can occur, and from our geometric analysis, this can be related to two phenomena which can be easily computed: i existence of a conjugate time, ii saturation of the control and birth of a bridge, which leads to a BSBS policy.

Numerical Continuation Method
A standard regularization process is the one used in the LQ-control which is an application of Tychonov theorem to control system 27 .The dynamics is linear of the form ẋ Ax Bu, and the cost can be written as t f 0 t uRu dt.The regular case corresponds to R > 0, and the singular case called cheap control is the case where R ≥ 0, and some control components u 1 , u 2 , . . ., u p are not penalized.A regularization process consists of adding quadratic penalties in the cost ε p i 1 u 2 i to get a regular problem for which an optimal solution is computed.This can be done by solving the shooting equation on the set of extremals solution of the maximum principle.The cheap case is obtained by taking ε → 0, and the solutions converge thanks to Tychonov theorem.In our nonlinear situation, the convergence is more intricate, but the contrast problem which is a Mayer problem can be interpreted as a cheap control problem.The regularization amounts to the standard homotopy in the cost: min u • C x t f 1 − λ t f 0 u t 2 dt , λ ∈ 0, 1 .For the regularized problem, one considers only normal extremals associated to the pseudo-Hamiltonian The control is computed in the normal case by solving ∂H λ /∂u 0, u n H G / 1 − λ .Saturation occurs if |u n | > 1, and the control is given by Sign H G / 1 − λ .
Another regularization process is to use the refined cost min u The analysis of the convergence is related to the comparison of the limit between the normal extremals of the regularized system and the original singular extremals.This point goes however beyond the scope of this paper.This regularizing process can replace the so-called GRAPE algorithm which is a standard gradient method commonly used in NMR 28 .To use the GRAPE algorithm in this situation, one replaces the boundary condition q 1 t f 0 and the cost −|q 2 2 t f | by a cost function of the form where α is a weight parameter.Another advantage of the homotopy method is, once the structure of the extremal trajectory is determined by the continuation method, the structure being a sequence of saturating controls and normal extremals, the true sequence BSBS . . .for  the contrast problem is computed, replacing each normal arc by a singular one.This leads to a true BC-extremal sequence satisfying the maximum principle and which is computed using a multiple shooting method.Such an algorithm is implemented in the HAMPATH code 29 .Moreover, the conjugate point test can be used for every singular arc of the sequence.

Numerical Results
We apply in this section the different tools presented above on a realistic NMR contrast problem corresponding to the cerebrospinal fluid/water case 11 .As stated below, we consider a simple model of two uncoupled spin 1/2 particles with different relaxation parameters.Each spin 1/2 particle is governed by the Bloch equation where the coordinates of the magnetization vector are M x , M y , M z , T 1 and T 2 are the relaxation times, and ω x , ω y are the two components of the magnetic field ω which is bounded by |ω| ≤ ω max .We use the normalized coordinates x, y, z M x , M y , M z /M 0 , which entails that the initial state is by definition the north pole 0, 0, 1 of the Bloch sphere.The normalized control is defined as u u x , u y 2π/ω max ω x , ω y , leading to |u| ≤ 2π, while the normalized time is given by τ ω max t/ 2π .In these new coordinates, the system takes the form with Γ 2π/ ω max T 2 and γ 2π/ ω max T 1 .From an experimental point of view, ω max / 2π can be chosen up to 15 000 Hz, but the value 32.3 Hz will be considered in this work.In the case of the cerebrospinal fluid/water sample, the relaxation parameters for the first spin describing the fluid are T 1 2000 ms and T 2 200 ms, while for the second spin T 1 T 2 2500 ms.We next present some numerical computations.We first analyze the structure of the singular flow.For that purpose, we plot in Figure 3 the projection of this flow onto the planes y 1 , z 1 and y 2 , z 2 .In this example, we assume that a bang pulse of large amplitude has been first applied to the system.The initial point of the singular flow is of coordinates , z 0 where z z 0 is the horizontal singular line of the first spin.Note that this first bang is necessary so that the singular trajectory of the first spin reaches the center of the Bloch ball, as expected in the contrast problem.One clearly sees in Figure 3 the excellent contrast that can be obtained with this very simple control field of the form bang-singular.Note that some singular control fields diverge as displayed in Figure 3.The attracting property of the north pole for the singular flow is represented on Figure 4.This singular flow is followed for longer times in Figure 5, which shows that this flow is attracted toward the north pole of the Bloch ball.We have also computed the position of the conjugate points for each singular extremal.Figures 5 and 6 show that the structure bang-singular is not optimal since the first conjugate point occurs before the time where the magnetization of the first spin vanishes.A more complicated pulse sequence composed of different bang and singular pulses has therefore to be used.In order to improve the contrast, a locally optimal BS-sequence solution of the maximum principle is computed as follows.The transfer time t f is fixed, and the Mayer problem is regularized by using a cost of the form The objective of this continuation is to determine the structure of the extremal control at the limit λ 1.We use the Hampath code 29 in the numerical simulations.Once this structure obtained, an extremal solution is computed for the initial optimal control problem in which the exact switching times are computed from a multiple shooting method.The results which are displayed in Figures 7-11 depend on the chosen transfer time.In Figure 7, we represent the evolution of the control field for different values of the homotopy parameter Advances in Mathematical Physics λ; the transfer time is fixed to 2T min where T min is the minimum time to transfer the first spin to 0. Figure 8 displays the control and the state variable when λ 0.93.This gives the structure BSBS of the extremal trajectory, as shown by the dots in this figure.The exact control is then computed from this solution in the limit Mayer case as can be seen in Figure 9 with the corresponding state trajectories.In Figure 10, we represent the same results but with t f 1.1T min , and the extremal solution has a more complicate sequence of the form BSBSBS.In Figure 11, we represent the evolution of the final contrast as a function of the transfer time.This contrast varies between 0.69 and 0.78 when the time increases from 1.1 to 2T min .
We conclude this paper by illustrating our numerical results with the first results of a contrast experiment.We consider two surfaces as displayed in Figure 12 filled in with spins 1 or 2 in a homogeneous manner.We apply the optimal control field, and we associate a color to the final modulus of the magnetization vector of the spin 2. This color is white if the modulus is equal to 1, black if it is zero, and a grey variant is in between.One clearly sees in Figure 12 the excellent contrasts that can be obtained in this example.
u s D 0 where D det G, G, F , F and D det G, G, F , G .Computing, one gets i for y 0, this simplifies into D 0 and D −z γ −2δz .Hence, the singular control is zero, and the dynamics is governed by the equation ż γ 1 − z , ii for z 0 γ/ 2δ , we get D −yγ 2δ − γ and D −2δy 2 , which gives γ − 2δ 2Γ − γ ≥ 0. The flow on the horizontal direction is given by ẏ −Γy − γ 2 2Γ − γ 4δ 2 y .4.6

Figure 1 :
Figure 1: Global synthesis from the north pole.The switching locus Σ is formed by the union of the surfaces Σ i .Due to the symmetry of revolution, the cut locus corresponds to the z-axis.

Figure 2 :
Figure 2: Plot of the optimal trajectories left and of the inversion recovery sequence right in the plane y, z .Experimental parameters are taken to be T 1 740 ms, T 2 60 ms, and ω max / 2π 32.3 Hz.The experimentally measured trajectories are represented by filled squares and open diamonds.In the upper panel, the small insert represents a zoom of the optimal trajectory near the origin.The dotted line is the switching curve originating from the horizontal singular line.The solid curve is the optimal trajectory near the origin.The lower panel displays the corresponding control laws.τ is the reduced time defined by τ ω max / 2π t.

Figure 3 :
Figure 3: Projection of the singular flow onto the planes y 1 , z 1 and y 2 , z 2 in the cerebrospinal fluid/water case.The trajectories are plotted in black solid line and in red dashed line .The control fields of the dashed extremals diverge.The trajectories have been plotted up to the explosion of the field the absolute value is then larger than 10 5 .The horizontal and vertical dashed lines represent the y-and z-axes of the Bloch ball.Note the contrast that can be obtained between the two spins, which illustrates the role of the singular flow in this problem.

Figure 4 :
Figure 4: It is the same as Figure 3, but the singular flow is followed for longer times to show the attracting property of the north pole for the singular flow.

Figure 5 :Figure 6 :Figure 7 :
Figure 5: It is the same as Figure 3, but the conjugate points are computed along the singular extremals.Their positions are represented by red crosses.

Figure 8 :
Figure 8: Trajectories of the magnetization vector for the first and second spins for the control field represented in the bottom panel.This field is the solution of the homotopy problem for λ 0.93.The red points indicate the possible structure of the solution in terms of bang and singular extremals see the text .

Figure 9 :
Figure 9: It is the same as Figure 8, but for a control field with a BSBS structure.Note the similarity of this result with the one of Figure 8.

Figure 10 :Figure 11 :Figure 12 :
Figure 10: It is the same as Figure 9 but for a control field of the form BSBSBS with t f 1.1T min .

2
Advances in Mathematical Physicswhere x ∈ R n , A, B 1 , . . ., B m are real constant n × n matrices, and a, b 1 , . . ., b m ∈ R n .In this paper, in view of applications, we will concentrate our analysis to two different cases, a time optimal control and a Mayer control problem which consist of minimizing a cost function C : min u • C x t f where t f is the transfer time.In both cases, according to the Maximum Principle 4 , if the control u • is constrained to a control domain U, minimizers are solutions of the Hamiltonian equations Advances in Mathematical Physicsi Time Minimum Control Problem.Reach in minimum time from x 0 the terminal manifold M.ii Mayer Problem.Steer x 0 to M while min u • ∈U C x t f where t f is a fixed transfer time, and the cost function C is a C ∞ regular function from R n into R. is the usual scalar product in R n , and p ∈ R n \ {0} is the adjoint vector.H is called the pseudo-Hamiltonian, and we introduce the maximized Hamiltonian M x, p max v∈U H x, p, v .
Theorem 3.4.Under our assumptions, the first geometric conjugate time t 1c is the first time t such that if t < t 1c , the input-output mapping E x 0 ,t is not open at u| 0,t , and if t > t 1c , it becomes open at u| 0,t for the L ∞ -norm topology.
(A2), (A3) are satisfied, then the first time t 1c such that the input-output mapping becomes open along x • which is the first geometric conjugate time.