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

Authors: Bernard Bonnard; Steffen J. Glaser; Dominique Sugny
Journal: Advances in Mathematical Physics
 (2012)

---

## Abstract

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.


---

## Body

## 1. Introduction

The dynamics of a spin 1/2 particle in nuclear magnetic resonance (NMR) is described by the Bloch equation(1.1)dM⃗dt=M⃗∧ω⃗-R(M⃗-M⃗0),
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](1.2)dxdt=(Ax+a)+∑i=1mui(Bix+bi),
where x∈ℝn, (A,B1,…,Bm) are real constant n×n matrices, and a,b1,…,bm∈ℝ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:minu(·)C(x(tf)) where tf 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(1.3)ẋ=∂H∂p,ṗ=-∂H∂x,
where H(x,p,u)=〈p,F(x,u)〉, ẋ=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](1.4)H(x,p,u)=max⁡v∈UH(x,p,v).
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 systemsẋ=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 ẋ=Ax+Bu, under the mild Kalman controllability condition rank [B,…,An-1B]=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 input-output 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(1.5)Ex0,tf:u(⋅)⟼x(tf,x0,u),
where x(·) is the response at time tf of the system to the control u(·), and x0 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(1.6)dxdt=-Γx+u2z,dydt=-Γy-u1z,dzdt=γ(1-z)+u1y-u2x,
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=u1+iu2=|u|eiφ 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:
(1.7)dq1dt=F1(q1,Λ1,u),dq2dt=F2(q2,Λ2,u),
written shortly as dx/dt=F(x,u) where x=(q1,q2) 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 x0=((0,0,1),(0,0,1)), the goal is to reach in a given transfer time tf the final state q1(tf)=0 (corresponding to zero magnetization of the first spin) while maximizing |q2(tf)|2 [11].


## 2. Necessary Optimality Conditions: Maximum Principle

We consider a control system of the form(2.1)dx(t)dt=F(x(t),u(t)),
where F:ℝn×ℝm↦ℝn is a C∞ mapping. The class of admissible control 𝒰 is the set {u:[0,tf]→U} of bounded and measurable maps where U⊂ℝm is the control domain. Given an initial condition x0∈ℝn and u(·)∈𝒰, we denote by x(t,x0,u) the solution of the control system initiating from x0 and defined on a subinterval of [0,tf]. We fix a C∞ terminal manifold ℳ defined by g(x)=0 where g:ℝn→ℝk. We will consider the two following optimal control problems.(i)
Time Minimum Control Problem. Reach in minimum time fromx0 the terminal manifold ℳ.(ii)
Mayer Problem. Steerx0 to ℳ while minu(·)∈𝒰C(x(tf)) where tf is a fixed transfer time, and the cost function C is a C∞ regular function from ℝn into ℝ.We denote the accessibility set at timetf by A(x0,tf)=∪u(·)∈𝒰x(tf,x0,u).Let us consider the Hamiltonian function(2.2)H(x,p,u)=〈p,F(x,u)〉,
where 〈·,·〉 is the usual scalar product in ℝn, and p∈ℝn∖{0} is the adjoint vector. H is called the pseudo-Hamiltonian, and we introduce the maximized Hamiltonian M(x,p)=maxv∈UH(x,p,v).It is well known that the Pontryagin maximum principle is a set of necessary conditions for our optimal control problem [3, 4].


### 2.1. The Pontryagin Maximum Principle

Letu*(·) be an admissible control whose corresponding trajectory x*(·)=x(·,u*,x0) is optimal, then there exists a vector function p*(·)∈ℝn∖{0} such that the following conditions are satisfied:(2.3)ẋ*=∂H∂p(x*,p*,u*),ṗ*=-∂H∂x(x*,p*,u*),
with the maximality condition(2.4)H(x*(t),p*(t),u*(t))=M(x*(t),p*(t)),
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 tf,
(2.5)g(x*(tf))=0,p*(tf)⊥∂g∂x(x*(tf)),(ii)
Mayer problem: at the final time tf,
(2.6)g(x*(tf))=0,p*(tf)=p0∂C∂x(x*(tf))+∑i=1kδi∂gi∂x(x*(tf)),δ=(δ1,…,δk)∈Rk,p0≤0.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(tf) to be in the boundary of the accessibility set A(x0,tf).The respective boundary conditions define the so-called BC-extremals:(i)
time minimum case: p*(tf)⊥ℳ,(ii)
Mayer problem: let us fix m∈ℝ and introduce the manifold ℳm={x,g(x)=0,C(x)=m}. One must have for the optimal solution p*(tf)⊥ℳm, where m is such that the cost is minimum.


### 2.2. 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.
Relaxingu∈U, a singular trajectory of the system ẋ=F(x,u) on [0,tf] is a control trajectory pair (x,u) such that the input-output mapping
(2.7)Ex0,tf:u(⋅)∈L∞[0,tf]⟼x(tf,x0,u)
is singular.Proposition 2.2.
One then deduces the following results.(1)
If(x,u) is not singular on [0,tf], then the input-output mapping is open at u(·) for the L∞[0,tf]-norm topology.(2)
If(x,u) is singular on [0,tf], then there exists p(·)∈ℝn∖{0} such that (x,p,u) satisfies a.e. on [0,tf] the weak maximum principle(2.8)ẋ=∂H∂p,ṗ=-∂H∂x,∂H∂u=0.Next, we make computations of extremals which will be used in the sequel.


### 2.3. Bi-Input Affine Case

A system such that the control enters linearly is called affine. We consider a situation for whichF(x,u)=F0(x)+u1F1(x)+u2F2(x), where the Fi are C∞-vector fields and u=(u1,u2), with the control domain U:|u|≤1. Let us denote by Hi=〈p,Fi(x)〉 the Hamiltonian lifts, then the pseudo-Hamiltonian H(x,p,u) takes the form H0+∑i=12uiHi. The maximality condition (2.4) leads to the following parameterization of the extremal controls:(2.9)u1=H1H12+H22,u2=H2H12+H22,
outside the switching surface Σ:H1=H2=0. Plugging such u into the Hamiltonian gives the true Hamiltonian Hn=H0+(H12+H22)1/2. The smooth solutions of the corresponding Hamiltonian field are called extremals of order zero.The relationu12+u22=1 leads to the following interpretation.Proposition 2.3.
Extremals of order zero correspond to the singularity of the input-output mapping(2.10)Ex0,tf:u∈L∞[0,tf]∩{|u|=1}⟼x(tf,x0,u).


### 2.4. Single-Input Affine Case

The system is of the formF(x,u)=F(x)+uG(x) and |u|≤1. Applying the maximality condition, one gets two types of extremals.(i)
Regular extremals: the control is given byu(t)=sign[HG(z(t))], z=(x,p), and HG=〈p,G〉. If u=±1, it is called a bang arc and if the number of switchings is finite, it is called bang-bang.(ii)
Singular extremals: since the system is linear inu, the condition (2.4) leads in the singular case to the condition HG=0, identically. Differentiating twice with respect to time, we get the relations
(2.11){HG,HF}={{HG,HF},HF}+u{{HG,HF},HG}=0,
and from this second condition, one derives when the denominator is not vanishing the corresponding singular control(2.12)uS=-{{HG,HF},HF}{{HG,HF},HG}.Plugging suchuS into the pseudo-Hamiltonian defines the singular flow solution of the Hamiltonian vector field denoted as HS and starting at t=0 from the two constraints HG={HG,HF}=0 and |uS|≤1. They have the following interpretation.Proposition 2.4.
Singular extremals correspond to the singularities of the input-output mapping(2.13)Ex0,tf:u∈L∞[0,tf]⟼x(tf,x0,u).


### 2.5. 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(2.14)∂2H∂u2≤0,
and if it is strict, it is called the strict Legendre-Clebsch condition. In the singular case, one has ∂2H/∂u2=0. A refined necessary condition has to be satisfied for optimality which is deduced from the higher-order maximum principle [12]. It is called the generalized Legendre-Clebsch condition(2.15)∂∂ud2dt2∂H∂u={HG,{HG,HF}}≤0.


### 2.6. 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 u12+u22=1, and one can consider the extension of the control system by setting u1=sinα, u2=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 thatG is nonzero, there exists a coordinate system (x1,…,xn) on an open set V such that G=∂/∂xn. Denoting x=(x′,xn), the system splits into(2.16)ẋ′=F(x′,xn),ẋn=F0′(x′)+u,
where the system F′ defined on an open set V′ with xn, 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 HamiltonianH′(x′,p′,xn)=〈p′,F′(x′,xn)〉, one has(2.17)ddt∂H∂u={HG,HF}=-∂H′∂xn,∂∂ud2dt2∂H∂u={{HG,HF},HG}=-∂2H′∂xn2.
These formulas establish the relation between the Legendre-Clebsch and the generalized Legendre-Clebsch of the corresponding systems.


## 2.1. The Pontryagin Maximum Principle

Letu*(·) be an admissible control whose corresponding trajectory x*(·)=x(·,u*,x0) is optimal, then there exists a vector function p*(·)∈ℝn∖{0} such that the following conditions are satisfied:(2.3)ẋ*=∂H∂p(x*,p*,u*),ṗ*=-∂H∂x(x*,p*,u*),
with the maximality condition(2.4)H(x*(t),p*(t),u*(t))=M(x*(t),p*(t)),
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 tf,
(2.5)g(x*(tf))=0,p*(tf)⊥∂g∂x(x*(tf)),(ii)
Mayer problem: at the final time tf,
(2.6)g(x*(tf))=0,p*(tf)=p0∂C∂x(x*(tf))+∑i=1kδi∂gi∂x(x*(tf)),δ=(δ1,…,δk)∈Rk,p0≤0.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(tf) to be in the boundary of the accessibility set A(x0,tf).The respective boundary conditions define the so-called BC-extremals:(i)
time minimum case: p*(tf)⊥ℳ,(ii)
Mayer problem: let us fix m∈ℝ and introduce the manifold ℳm={x,g(x)=0,C(x)=m}. One must have for the optimal solution p*(tf)⊥ℳm, where m is such that the cost is minimum.


## 2.2. 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.
Relaxingu∈U, a singular trajectory of the system ẋ=F(x,u) on [0,tf] is a control trajectory pair (x,u) such that the input-output mapping
(2.7)Ex0,tf:u(⋅)∈L∞[0,tf]⟼x(tf,x0,u)
is singular.Proposition 2.2.
One then deduces the following results.(1)
If(x,u) is not singular on [0,tf], then the input-output mapping is open at u(·) for the L∞[0,tf]-norm topology.(2)
If(x,u) is singular on [0,tf], then there exists p(·)∈ℝn∖{0} such that (x,p,u) satisfies a.e. on [0,tf] the weak maximum principle(2.8)ẋ=∂H∂p,ṗ=-∂H∂x,∂H∂u=0.Next, we make computations of extremals which will be used in the sequel.


## 2.3. Bi-Input Affine Case

A system such that the control enters linearly is called affine. We consider a situation for whichF(x,u)=F0(x)+u1F1(x)+u2F2(x), where the Fi are C∞-vector fields and u=(u1,u2), with the control domain U:|u|≤1. Let us denote by Hi=〈p,Fi(x)〉 the Hamiltonian lifts, then the pseudo-Hamiltonian H(x,p,u) takes the form H0+∑i=12uiHi. The maximality condition (2.4) leads to the following parameterization of the extremal controls:(2.9)u1=H1H12+H22,u2=H2H12+H22,
outside the switching surface Σ:H1=H2=0. Plugging such u into the Hamiltonian gives the true Hamiltonian Hn=H0+(H12+H22)1/2. The smooth solutions of the corresponding Hamiltonian field are called extremals of order zero.The relationu12+u22=1 leads to the following interpretation.Proposition 2.3.
Extremals of order zero correspond to the singularity of the input-output mapping(2.10)Ex0,tf:u∈L∞[0,tf]∩{|u|=1}⟼x(tf,x0,u).


## 2.4. Single-Input Affine Case

The system is of the formF(x,u)=F(x)+uG(x) and |u|≤1. Applying the maximality condition, one gets two types of extremals.(i)
Regular extremals: the control is given byu(t)=sign[HG(z(t))], z=(x,p), and HG=〈p,G〉. If u=±1, it is called a bang arc and if the number of switchings is finite, it is called bang-bang.(ii)
Singular extremals: since the system is linear inu, the condition (2.4) leads in the singular case to the condition HG=0, identically. Differentiating twice with respect to time, we get the relations
(2.11){HG,HF}={{HG,HF},HF}+u{{HG,HF},HG}=0,
and from this second condition, one derives when the denominator is not vanishing the corresponding singular control(2.12)uS=-{{HG,HF},HF}{{HG,HF},HG}.Plugging suchuS into the pseudo-Hamiltonian defines the singular flow solution of the Hamiltonian vector field denoted as HS and starting at t=0 from the two constraints HG={HG,HF}=0 and |uS|≤1. They have the following interpretation.Proposition 2.4.
Singular extremals correspond to the singularities of the input-output mapping(2.13)Ex0,tf:u∈L∞[0,tf]⟼x(tf,x0,u).


## 2.5. 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(2.14)∂2H∂u2≤0,
and if it is strict, it is called the strict Legendre-Clebsch condition. In the singular case, one has ∂2H/∂u2=0. A refined necessary condition has to be satisfied for optimality which is deduced from the higher-order maximum principle [12]. It is called the generalized Legendre-Clebsch condition(2.15)∂∂ud2dt2∂H∂u={HG,{HG,HF}}≤0.


## 2.6. 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 u12+u22=1, and one can consider the extension of the control system by setting u1=sinα, u2=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 thatG is nonzero, there exists a coordinate system (x1,…,xn) on an open set V such that G=∂/∂xn. Denoting x=(x′,xn), the system splits into(2.16)ẋ′=F(x′,xn),ẋn=F0′(x′)+u,
where the system F′ defined on an open set V′ with xn, 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 HamiltonianH′(x′,p′,xn)=〈p′,F′(x′,xn)〉, one has(2.17)ddt∂H∂u={HG,HF}=-∂H′∂xn,∂∂ud2dt2∂H∂u={{HG,HF},HG}=-∂2H′∂xn2.
These formulas establish the relation between the Legendre-Clebsch and the generalized Legendre-Clebsch of the corresponding systems.


## 3. 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 isẋ=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,tf] of(3.1)ẋ=∂H∂p(x,p,u),ṗ=-∂H∂x(x,p,u),∂H∂u=0.Our first assumption is the strong-Legendre condition:(A1)
The quadratic form∂2H/∂u2 is negative definite along the reference solution.Using the implicit function theorem, the extremal control is then locally defined as a smooth functionu(z), and plugging this function into H defines a smooth true Hamiltonian, still denoted as H, and the extremal is solution of ż=H⃗(z) with initial condition z0=(x0,p0).


### 3.1. The Geometric Concept of Conjugate Point

Definition 3.1.
Letz=(x,p) be the reference extremal defined on [0,tf]. The variational equation
(3.2)δż=dH⃗(z(t))δz
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.
LetL0 be the fiber Tx0*X, and let Lt=exp[tH⃗(L0)] be its image by the one-parameter subgroup generated by H, then Lt 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 Lt of Π is at most n-1.This leads to the following definition.Definition 3.3.
We define the exponential mapping(3.3)exp⁡x0,t(p0)=Π(z(t,x0,p0))
as the projection on the state space of the integral curve of H⃗ with initial condition z0=(x0,p0), where p0 can be restricted by linearity to the sphere |p0|=1. If z=(x,p) is the reference extremal, a time tc>0 is said to be geometrically conjugate to 0 if the mapping p0↦expx0,t is not of rank (n-1) at t=tc, and the associated point x(tc) is said to be geometrically conjugate to x0. We denote by t1c the first conjugate time and C(x0) is the conjugate locus formed by the set of first conjugate points.An algorithm can be deduced.Testing Conjugacy
Letz(t)=(x(t),p(t)) be the reference extremal, and consider the vector space of dimension n-1 generated by the Jacobi fields δzi=(δxi,δpi), i=1,…,n-1 vertical at t=0 and such that δpi(0) is orthogonal to p0. At a conjugate time tc, one has
(3.4)rank⁡[δx1(tc),…,δxn-1(tc)]<n-1.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 trajectoryz(t), t∈[0,tf] associated to u(·) is by definition a singularity of the input-output mapping Ex0,tf. One assumes that on each subinterval [t0,t1], 0<t0<t1≤tf, the singularity of E on u|[t0,t1] is of codimension one.(A3)
We are in the nonexceptional case whereH≠0 along the reference extremal.Theorem 3.4.
Under our assumptions, the first geometric conjugate timet1c is the first time t such that if t<t1c, the input-output mapping Ex0,t is not open at u|[0,t], and if t>t1c,  it becomes open at u|[0,t] for the L∞-norm topology.Corollary 3.5.
Consider the control systemF(x,u)=F0(x)+u1F1(x)+u2F2(x), |u|≤1, and let (x,p,u) be a reference extremal of order zero solution of H0+(H12+H22)1/2 such that assumptions (A2), (A3) are satisfied, then the first time t1c such that the input-output mapping becomes open along x(·) which is the first geometric conjugate time.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 systemẋ=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.


### 3.2. Focal Type Conditions

The concept of conjugate point is associated to optimal control problems with fixed end-points conditions. In the contrast problem, it has to be adapted to the problem with fixed initial conditionx(0)=x0, but with variable final condition x(tf)∈ℳm. From the maximum principle, the transversality condition defined a Lagrangian manifold z(tf)=(x(tf),p(tf))∈ℳm⊥ which again defines along a singular extremal a train of Lagrangian manifolds Lt, integrating backwards t≤tf. Optimality can be deduced from the geometric concept of focal points corresponding to the singularities of (Lt,Π).


### 3.3. Extremal Fields and Hamilton-Jacobi-Bellman Equation

For simplicity, we consider here the time-minimal control problem. We pick up a reference extremal trajectoryt↦x(t), t∈[0,tf] 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,tf], then we can embed locally the reference extremal into a central field ℱ formed by all extremal trajectories starting from x0=x(0). One restricts the adjoint vector with the normalization H=1 and the assumption (A3). One introduces(3.5)L={(H=1)∩(∪t≥0Lt)},
where Lt is the train of Lagrangian manifolds generated by the Hamiltonian flow H⃗, with L0=Tx0*X being the initial condition. By construction, ℱ is the projection of L on the state space X.Proposition 3.6.
Excludingx0, there exists an open neighborhood W of the reference trajectory and two smooth mappings V:W↦ℝ, û:W→U such that for each (x,u)∈W×U, the maximization condition
(3.6)H(x,dV(x),û(x))≥H(x,dV(x),u).
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.


## 3.1. The Geometric Concept of Conjugate Point

Definition 3.1.
Letz=(x,p) be the reference extremal defined on [0,tf]. The variational equation
(3.2)δż=dH⃗(z(t))δz
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.
LetL0 be the fiber Tx0*X, and let Lt=exp[tH⃗(L0)] be its image by the one-parameter subgroup generated by H, then Lt 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 Lt of Π is at most n-1.This leads to the following definition.Definition 3.3.
We define the exponential mapping(3.3)exp⁡x0,t(p0)=Π(z(t,x0,p0))
as the projection on the state space of the integral curve of H⃗ with initial condition z0=(x0,p0), where p0 can be restricted by linearity to the sphere |p0|=1. If z=(x,p) is the reference extremal, a time tc>0 is said to be geometrically conjugate to 0 if the mapping p0↦expx0,t is not of rank (n-1) at t=tc, and the associated point x(tc) is said to be geometrically conjugate to x0. We denote by t1c the first conjugate time and C(x0) is the conjugate locus formed by the set of first conjugate points.An algorithm can be deduced.Testing Conjugacy
Letz(t)=(x(t),p(t)) be the reference extremal, and consider the vector space of dimension n-1 generated by the Jacobi fields δzi=(δxi,δpi), i=1,…,n-1 vertical at t=0 and such that δpi(0) is orthogonal to p0. At a conjugate time tc, one has
(3.4)rank⁡[δx1(tc),…,δxn-1(tc)]<n-1.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 trajectoryz(t), t∈[0,tf] associated to u(·) is by definition a singularity of the input-output mapping Ex0,tf. One assumes that on each subinterval [t0,t1], 0<t0<t1≤tf, the singularity of E on u|[t0,t1] is of codimension one.(A3)
We are in the nonexceptional case whereH≠0 along the reference extremal.Theorem 3.4.
Under our assumptions, the first geometric conjugate timet1c is the first time t such that if t<t1c, the input-output mapping Ex0,t is not open at u|[0,t], and if t>t1c,  it becomes open at u|[0,t] for the L∞-norm topology.Corollary 3.5.
Consider the control systemF(x,u)=F0(x)+u1F1(x)+u2F2(x), |u|≤1, and let (x,p,u) be a reference extremal of order zero solution of H0+(H12+H22)1/2 such that assumptions (A2), (A3) are satisfied, then the first time t1c such that the input-output mapping becomes open along x(·) which is the first geometric conjugate time.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 systemẋ=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.


## 3.2. Focal Type Conditions

The concept of conjugate point is associated to optimal control problems with fixed end-points conditions. In the contrast problem, it has to be adapted to the problem with fixed initial conditionx(0)=x0, but with variable final condition x(tf)∈ℳm. From the maximum principle, the transversality condition defined a Lagrangian manifold z(tf)=(x(tf),p(tf))∈ℳm⊥ which again defines along a singular extremal a train of Lagrangian manifolds Lt, integrating backwards t≤tf. Optimality can be deduced from the geometric concept of focal points corresponding to the singularities of (Lt,Π).


## 3.3. Extremal Fields and Hamilton-Jacobi-Bellman Equation

For simplicity, we consider here the time-minimal control problem. We pick up a reference extremal trajectoryt↦x(t), t∈[0,tf] 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,tf], then we can embed locally the reference extremal into a central field ℱ formed by all extremal trajectories starting from x0=x(0). One restricts the adjoint vector with the normalization H=1 and the assumption (A3). One introduces(3.5)L={(H=1)∩(∪t≥0Lt)},
where Lt is the train of Lagrangian manifolds generated by the Hamiltonian flow H⃗, with L0=Tx0*X being the initial condition. By construction, ℱ is the projection of L on the state space X.Proposition 3.6.
Excludingx0, there exists an open neighborhood W of the reference trajectory and two smooth mappings V:W↦ℝ, û:W→U such that for each (x,u)∈W×U, the maximization condition
(3.6)H(x,dV(x),û(x))≥H(x,dV(x),u).
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.


## 4. 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,(4.1)dydt=-Γy-uz,dzdt=γ(1-z)+uy,|u|≤m.
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 pointq0=(y,z), find the time minimal solution in a small neighborhood V of q0, according to the Lie algebraic structure of the system at q0;(2)
glue together the local solutions to get the global time optimal solution. More precisely, we will describe the global synthesis to steerN=(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.


### 4.1. Computations of Lie Brackets and Singular Trajectories Properties

One has(4.2)F=-Γy∂∂y+γ(1-z)∂∂z,G=-z∂∂y+y∂∂z
and the Lie brackets up to order three are(4.3)[G,F]=(-γ+δz)∂∂y+δy∂∂z,
where δ=γ-Γ,(4.4)[[G,F],F]=(γ(γ-2Γ)-δ2z)∂∂y+δ2y∂∂z,[[G,F],G]=2δy∂∂y+(γ-2δz)∂∂z.
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 z0=γ/(2δ). Eliminating p, the singular control is given by D′+usD=0 where D′=det(G,[[G,F],F]) and D=det(G,[[G,F],G]). Computing, one gets(i)
fory=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 ż=γ(1-z),(ii)
forz0=γ/(2δ), we get D′=-yγ(2δ-γ) and D=-2δy2, which gives
(4.5)uS=γ(γ-2δ)2δy,
where γ-2δ=2Γ-γ≥0. The flow on the horizontal direction is given by(4.6)ẏ=-Γy-γ22Γ-γ4δ2y.


### 4.2. Parameter Conditions

The interesting case is whenδ<0 and |γ/2δ|<1. This leads to(4.7)Γ>3γ2.
In this case, uS→+∞ when y→0-. Moreover, we have the following estimate. Setting α=γ2(2Γ-γ)/(4γ2), we get near 0 that ẏ≃-α/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 pointy=0, along the horizontal direction is of order O(1/t) and is in the L1 but not L2 category.


### 4.3. 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 haveH≥0 from the maximum principle, and from the Legendre-Clebsch condition (2.15), we deduce that {HG,{HG,HF}}≤0. More precisely, if we introduce D′′=det(G,F)=γz(z-1)+Γy2, one arrives at(i)
hyperbolic caseDD′′>0, the singular direction is small time minimizing;(ii)
elliptic caseDD′′<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 providedz0<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<z0 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)+Γy2=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.
Assumem=+∞, 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].


### 4.4. Classification of Regular Extremals NearΣ′

The singular extremals are contained in the subsetHG={HG,HF}=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)={HG,HF}(w(t))≠0 where Φ(t)=HG(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., δ-δ+) if Φ̇(t)>0 (resp., Φ̇(t)<0).Fold case: denoting Φ̈±(t)={{HG,HF},HF}±{{HG,HF},HG}(z(t)) 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 |uS|<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 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: δ+δ-δ+ if Φ̇+(t)>0 or δ-δ+δ- if Φ̇-(t)<0.This classification is far from being sufficient to analyze locally our problems. In particular, we have the following situations.


### 4.5. Saturating Singular Case

It is a transition between the hyperbolic and the parabolic cases. The singular controluS 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(4.8)ẋ=1-y2,ẏ=-(x+1)+u,|u|≤1.
The singular control is uS=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.


### 4.6. 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 z0=γ/(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 pointB, 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].


### 4.7. 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(4.9)〈p(0),G(q(0))〉=〈p(t),G(q(t))〉=0.
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(4.10)V(t)=e-tad(F+εG)G(q(t)),
and in the analytic case, the ad-formula gives, for small t,(4.11)V(0)=∑n≥0(-t)nn!adn(F+εG)G(q(t)).
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,ℝ)×Sℝ2 identified with the set of matrices of GL(3,ℝ):(4.12){(10gw),g∈GL(2,R),w∈R2},
acting on the subspace of vector in ℝ3:(1q). Lie brackets of two affine vector fields (A′,a′), (B′,b′) being defined by(4.13)((A′,a′),B′,b′)=([A′,B′],A′b′-B′a′),
the computation of the exponential exp[-tad(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.


### 4.8. 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 Figure1. 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.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.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.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 T1=740ms, T2=60ms, and ωmax/(2π)=32.3Hz. 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.


## 4.1. Computations of Lie Brackets and Singular Trajectories Properties

One has(4.2)F=-Γy∂∂y+γ(1-z)∂∂z,G=-z∂∂y+y∂∂z
and the Lie brackets up to order three are(4.3)[G,F]=(-γ+δz)∂∂y+δy∂∂z,
where δ=γ-Γ,(4.4)[[G,F],F]=(γ(γ-2Γ)-δ2z)∂∂y+δ2y∂∂z,[[G,F],G]=2δy∂∂y+(γ-2δz)∂∂z.
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 z0=γ/(2δ). Eliminating p, the singular control is given by D′+usD=0 where D′=det(G,[[G,F],F]) and D=det(G,[[G,F],G]). Computing, one gets(i)
fory=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 ż=γ(1-z),(ii)
forz0=γ/(2δ), we get D′=-yγ(2δ-γ) and D=-2δy2, which gives
(4.5)uS=γ(γ-2δ)2δy,
where γ-2δ=2Γ-γ≥0. The flow on the horizontal direction is given by(4.6)ẏ=-Γy-γ22Γ-γ4δ2y.


## 4.2. Parameter Conditions

The interesting case is whenδ<0 and |γ/2δ|<1. This leads to(4.7)Γ>3γ2.
In this case, uS→+∞ when y→0-. Moreover, we have the following estimate. Setting α=γ2(2Γ-γ)/(4γ2), we get near 0 that ẏ≃-α/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 pointy=0, along the horizontal direction is of order O(1/t) and is in the L1 but not L2 category.


## 4.3. 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 haveH≥0 from the maximum principle, and from the Legendre-Clebsch condition (2.15), we deduce that {HG,{HG,HF}}≤0. More precisely, if we introduce D′′=det(G,F)=γz(z-1)+Γy2, one arrives at(i)
hyperbolic caseDD′′>0, the singular direction is small time minimizing;(ii)
elliptic caseDD′′<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 providedz0<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<z0 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)+Γy2=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.
Assumem=+∞, 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].


## 4.4. Classification of Regular Extremals NearΣ′

The singular extremals are contained in the subsetHG={HG,HF}=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)={HG,HF}(w(t))≠0 where Φ(t)=HG(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., δ-δ+) if Φ̇(t)>0 (resp., Φ̇(t)<0).Fold case: denoting Φ̈±(t)={{HG,HF},HF}±{{HG,HF},HG}(z(t)) 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 |uS|<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 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: δ+δ-δ+ if Φ̇+(t)>0 or δ-δ+δ- if Φ̇-(t)<0.This classification is far from being sufficient to analyze locally our problems. In particular, we have the following situations.


## 4.5. Saturating Singular Case

It is a transition between the hyperbolic and the parabolic cases. The singular controluS 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(4.8)ẋ=1-y2,ẏ=-(x+1)+u,|u|≤1.
The singular control is uS=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.


## 4.6. 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 z0=γ/(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 pointB, 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].


## 4.7. 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(4.9)〈p(0),G(q(0))〉=〈p(t),G(q(t))〉=0.
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(4.10)V(t)=e-tad(F+εG)G(q(t)),
and in the analytic case, the ad-formula gives, for small t,(4.11)V(0)=∑n≥0(-t)nn!adn(F+εG)G(q(t)).
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,ℝ)×Sℝ2 identified with the set of matrices of GL(3,ℝ):(4.12){(10gw),g∈GL(2,R),w∈R2},
acting on the subspace of vector in ℝ3:(1q). Lie brackets of two affine vector fields (A′,a′), (B′,b′) being defined by(4.13)((A′,a′),B′,b′)=([A′,B′],A′b′-B′a′),
the computation of the exponential exp[-tad(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.


## 4.8. 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 Figure1. 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.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.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.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 T1=740ms, T2=60ms, and ωmax/(2π)=32.3Hz. 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.


## 5. The Contrast Problem

For simplicity, we will limit our analysis to the single-input case. The system can be written shortly as(5.1)dxdt=F(x)+uG(x),
where x=(q1,q2), qi=(yi,zi) representing the state of each spin 1/2 particle. From Bloch equation, we get(5.2)F=∑i=12(-yiΓi∂∂yi+γi(1-zi)∂∂zi),G=∑i=12(-zi∂∂yi+yi∂∂zi),
where Λi=(Γi,γi) are the relaxation parameters characterizing each spin.


### 5.1. BC Singular Extremals

According to the general computations of Section2, a singular extremal is contained in the constrained set(5.3)Σ′′:〈p,G(x)〉=〈p,[G,F](x)〉=0,
while the singular control is defined by(5.4)uS=-〈p,[[G,F],F](x)〉〈p,[[G,F],G](x)〉,
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=(p1,p2) where pi is associated to the spin i, one has(5.5)p2(tf)=-2p0q2(tf),p0≤0,
and if p0≠0, it can be normalized by homogeneity to p0=-1/2. The case p0=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 timetf is fixed to the optimal time.


### 5.2. Exceptional Singular Extremals

If the transfer time is not fixed, then according to the maximum principle, this leads to the additional constraintM=maxu(·)[HF+uHG]=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(5.6)uSe=-D′(x)D(x),
where(5.7)D=det⁡(F,G,[G,F],[[G,F],G]),D′=det⁡(F,G,[G,F],[[G,F],F]),
with the corresponding vector fields defined by(5.8)dxdt=F(x)-D′(x)D(x)G(x).
Observe that in the general case, a similar computation shows that the singular extremals are solutions of an equation of the form(5.9)dxdt=F(x)-D′(x,λ)D(x,λ)G(x),
where λ is a one-dimensional time-dependent parameter whose dynamics is deduced from the adjoint equation.


### 5.3. Desingularization

Meromorphic differential equations of the form (5.8) can be desingularized as(5.10)dxds=DF(x)-D′G(x),
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.


### 5.4. Feedback Classification and the Contrast Problem

An important object to analyze our problem is the action of the feedback groupGf on the set of systems. One will restrict our presentation to the single-input affine case, and we denote by 𝒞=(F,G) the set of all such (smooth) systems on a state space X≃ℝn.Definition 5.2.
Let(F,G), (F′,G′) be two elements of 𝒞. They are called feedback equivalent if there exists a smooth diffeomorphism φ of ℝn on a feedback u=α(x)+uβ(x) where α, β are smooth, β is invertible such that
(5.11)(i)F′=φ*F+φ*(G⋅α),(ii)G′=φ*(G⋅β),
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 𝒢.Definition 5.3.
Let(F,G)∈𝒞 and let λ be the map which associates to (F,G) the constrained differential equation whose solutions are singular extremals. The constraint is given by Σ′′:HG={HG,HF}=0, and the equation is FS=F+uSG where uS=-{{HG,HF},HF}/{{HG,HF},HG}. We define the action of (φ,α,β) on (Σ,FS) by the action of the symplectic change of coordinates
(5.12)φ⃗:x=φ(X),p=P∂φ-1∂X,
a feedback acting trivially.We have the following theorem [23].Theorem 5.4.
The following diagram is commutative:(5.13)C⟶λ(C)Gf↓↓GfC⟶λ(C).
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.


### 5.5. Classification of the DistributionG

A feedbacku=α(x)+β(x)v acts on the control direction G as G→βG, which corresponds to the classification of the one-dimensional distribution x→SpanG(x).


### 5.6. Collinear Set

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


### 5.7. The Singular Trajectories

It is defined by the constrained Hamiltonian equations given by(5.14)ẋ=F(x)+uS(x,p)G(x),ṗ=-p(∂F∂x+uS(x,p)∂G∂x(x)),
restricted to the surface(5.15)Σ′′:HG={HG,HF}=0,
where uS 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:(5.16)D={{HG,HF},HG},D′={{HG,HF},HF},
and the differential equation (5.14) can be desingularized using the time reparameterization ds=dt/(𝒟(x(t),p(t))). One gets the equations(5.17)dxds=DF-D′G,dpds=-p(D∂F∂x-D′∂G∂x),
which restricted to Σ′′ are a semicovariant.Geometric invariants are related to the surface{𝒟=𝒟′=0}∩Σ′′ and to the singularities of (5.17) in this surface. The geometric framework to analyze such nonisolated singularities is presented in [24, 25].


## 5.1. BC Singular Extremals

According to the general computations of Section2, a singular extremal is contained in the constrained set(5.3)Σ′′:〈p,G(x)〉=〈p,[G,F](x)〉=0,
while the singular control is defined by(5.4)uS=-〈p,[[G,F],F](x)〉〈p,[[G,F],G](x)〉,
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=(p1,p2) where pi is associated to the spin i, one has(5.5)p2(tf)=-2p0q2(tf),p0≤0,
and if p0≠0, it can be normalized by homogeneity to p0=-1/2. The case p0=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 timetf is fixed to the optimal time.


## 5.2. Exceptional Singular Extremals

If the transfer time is not fixed, then according to the maximum principle, this leads to the additional constraintM=maxu(·)[HF+uHG]=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(5.6)uSe=-D′(x)D(x),
where(5.7)D=det⁡(F,G,[G,F],[[G,F],G]),D′=det⁡(F,G,[G,F],[[G,F],F]),
with the corresponding vector fields defined by(5.8)dxdt=F(x)-D′(x)D(x)G(x).
Observe that in the general case, a similar computation shows that the singular extremals are solutions of an equation of the form(5.9)dxdt=F(x)-D′(x,λ)D(x,λ)G(x),
where λ is a one-dimensional time-dependent parameter whose dynamics is deduced from the adjoint equation.


## 5.3. Desingularization

Meromorphic differential equations of the form (5.8) can be desingularized as(5.10)dxds=DF(x)-D′G(x),
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.


## 5.4. Feedback Classification and the Contrast Problem

An important object to analyze our problem is the action of the feedback groupGf on the set of systems. One will restrict our presentation to the single-input affine case, and we denote by 𝒞=(F,G) the set of all such (smooth) systems on a state space X≃ℝn.Definition 5.2.
Let(F,G), (F′,G′) be two elements of 𝒞. They are called feedback equivalent if there exists a smooth diffeomorphism φ of ℝn on a feedback u=α(x)+uβ(x) where α, β are smooth, β is invertible such that
(5.11)(i)F′=φ*F+φ*(G⋅α),(ii)G′=φ*(G⋅β),
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 𝒢.Definition 5.3.
Let(F,G)∈𝒞 and let λ be the map which associates to (F,G) the constrained differential equation whose solutions are singular extremals. The constraint is given by Σ′′:HG={HG,HF}=0, and the equation is FS=F+uSG where uS=-{{HG,HF},HF}/{{HG,HF},HG}. We define the action of (φ,α,β) on (Σ,FS) by the action of the symplectic change of coordinates
(5.12)φ⃗:x=φ(X),p=P∂φ-1∂X,
a feedback acting trivially.We have the following theorem [23].Theorem 5.4.
The following diagram is commutative:(5.13)C⟶λ(C)Gf↓↓GfC⟶λ(C).
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.


## 5.5. Classification of the DistributionG

A feedbacku=α(x)+β(x)v acts on the control direction G as G→βG, which corresponds to the classification of the one-dimensional distribution x→SpanG(x).


## 5.6. Collinear Set

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


## 5.7. The Singular Trajectories

It is defined by the constrained Hamiltonian equations given by(5.14)ẋ=F(x)+uS(x,p)G(x),ṗ=-p(∂F∂x+uS(x,p)∂G∂x(x)),
restricted to the surface(5.15)Σ′′:HG={HG,HF}=0,
where uS 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:(5.16)D={{HG,HF},HG},D′={{HG,HF},HF},
and the differential equation (5.14) can be desingularized using the time reparameterization ds=dt/(𝒟(x(t),p(t))). One gets the equations(5.17)dxds=DF-D′G,dpds=-p(D∂F∂x-D′∂G∂x),
which restricted to Σ′′ are a semicovariant.Geometric invariants are related to the surface{𝒟=𝒟′=0}∩Σ′′ and to the singularities of (5.17) in this surface. The geometric framework to analyze such nonisolated singularities is presented in [24, 25].


## 6. 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 aC0-synthesis in the limit case m≫Γ,γ. This result is based on [9, 26].


### 6.1. Synthesis for a Reference Singular Arc

We consider the control systemdx/dt=F(x)+uG(x), |u|≤1, and let uS be a smooth singular control such that |uS|<1 with corresponding trajectory x(t). We have the following.Theorem 6.1.
Under generic assumptions, the first conjugate timet1c is the time along the singular arc such that a policy δ±δSδ± is time extremizing in a C0-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 isu=+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 tf, one must have p2(tf)=q2(tf), and moreover HG(tf)=0. The final point concentrates a switching. Hence, one gets the following.Lemma 6.2.
The simplest possible sequence is of the formBS,  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.


### 6.2. 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 ẋ=Ax+Bu, and the cost can be written as ∫0tftuRudt. The regular case corresponds to R>0, and the singular case called cheap control is the case where R≥0, and some control components (u1,u2,…,up) are not penalized. A regularization process consists of adding quadratic penalties in the cost ε(∑i=1pui2) 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:minu(·)(C(x(tf))+(1-λ)∫0tfu(t)2dt), λ∈[0,1]. For the regularized problem, one considers only normal extremals associated to the pseudo-Hamiltonian(6.1)Hλ=HF+uHG-12(1-λ)u2.
The control is computed in the normal case by solving ∂Hλ/∂u=0, un=HG/(1-λ). Saturation occurs if |un|>1, and the control is given by (SignHG)/(1-λ).Another regularization process is to use the refined costminu(·)(C(x(tf))+(1-λ)∫0tf|u|2-λdt). 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 q1(tf)=0 and the cost -|q22(tf)| by a cost function of the form(6.2)Cα=α|q12(tf)|-|q22(tf)|,
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.


### 6.3. 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(6.3)dMxdt=-MxT2+ωyMz,dMydt=-MyT2-ωxMz,dMzdt=(M0-Mz)T1+ωxMy-ωyMx,
where the coordinates of the magnetization vector are (Mx,My,Mz), T1 and T2 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)=(Mx,My,Mz)/M0, 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=(ux,uy)=(2π/ωmax)(ωx,ωy), leading to |u|≤2π, while the normalized time is given by τ=ωmaxt/(2π). In these new coordinates, the system takes the form(6.4)ẋ=-Γx+uyz,ẏ=-Γy-uxz,ż=γ(1-z)+uxy-uyx,
with Γ=2π/(ωmaxT2) and γ=2π/(ωmaxT1). 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 T1=2000 ms and T2=200 ms, while for the second spin T1=T2=2500 ms.We next present some numerical computations. We first analyze the structure of the singular flow. For that purpose, we plot in Figure3 the projection of this flow onto the planes (y1,z1) and (y2,z2). 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 ((-1-z02,z0),(-1-z02,z0)) where z=z0 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.Projection of the singular flow onto the planes(y1,z1) and (y2,z2) 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 105). 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.
(a)(b)It is the same as Figure3, but the singular flow is followed for longer times to show the attracting property of the north pole for the singular flow.
(a)(b)It is the same as Figure3, but the conjugate points are computed along the singular extremals. Their positions are represented by red crosses.
(a)(b)Zoom of the results of Figure5.
(a)(b)In order to improve the contrast, a locally optimal BS-sequence solution of the maximum principle is computed as follows. The transfer timetf is fixed, and the Mayer problem is regularized by using a cost of the form(6.5)C(x(tf))+(1-λ)∫0tf|u|2-λdt.
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 λ; the transfer time is fixed to 2Tmin where Tmin 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 tf=1.1Tmin, 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 2Tmin.Evolution of the control fieldu for different values of the homotopy parameter λ, which is equal to 0, 0.6, 0.8, 0.9 from (a) to (d). The transfer time normalized to 1 is equal to 2Tmin.
(a)(b)(c)(d)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).
(a)(b)(c)It is the same as Figure8, but for a control field with a BSBS structure. Note the similarity of this result with the one of Figure 8.
(a)(b)(c)It is the same as Figure9 but for a control field of the form BSBSBS with tf=1.1Tmin.
(a)(b)(c)Figure 11
Evolution of the final contrast as a function of the transfer time.We conclude this paper by illustrating our numerical results with the first results of a contrast experiment. We consider two surfaces as displayed in Figure12 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.Experimental results on the contrast problems. The inner disk mimics the spin 1, while the outside ring mimics the spin 2. The two surfaces are separated by a thin black circle. (a) is a reference image before the application of the control field when the two spins are at the north pole of the Bloch sphere. (b) represent the contrast after a real experiment. The panel (b) corresponds to|My|, while the panel (c) to |Mz|, that is, the remaining y and z magnetization components after the optimized pulse. A color has been associated to each value of the contrast between 0 and 1, 0 and 1 corresponding, respectively, to the colors black and white. In the experiment, relaxation parameters are taken to be T1=1.3 ms, T2=50 ms for the inner sample and T1=1.3 ms, T2=200 ms for the outside sample.
(a)(b)(c)


## 6.1. Synthesis for a Reference Singular Arc

We consider the control systemdx/dt=F(x)+uG(x), |u|≤1, and let uS be a smooth singular control such that |uS|<1 with corresponding trajectory x(t). We have the following.Theorem 6.1.
Under generic assumptions, the first conjugate timet1c is the time along the singular arc such that a policy δ±δSδ± is time extremizing in a C0-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 isu=+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 tf, one must have p2(tf)=q2(tf), and moreover HG(tf)=0. The final point concentrates a switching. Hence, one gets the following.Lemma 6.2.
The simplest possible sequence is of the formBS,  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.


## 6.2. 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 ẋ=Ax+Bu, and the cost can be written as ∫0tftuRudt. The regular case corresponds to R>0, and the singular case called cheap control is the case where R≥0, and some control components (u1,u2,…,up) are not penalized. A regularization process consists of adding quadratic penalties in the cost ε(∑i=1pui2) 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:minu(·)(C(x(tf))+(1-λ)∫0tfu(t)2dt), λ∈[0,1]. For the regularized problem, one considers only normal extremals associated to the pseudo-Hamiltonian(6.1)Hλ=HF+uHG-12(1-λ)u2.
The control is computed in the normal case by solving ∂Hλ/∂u=0, un=HG/(1-λ). Saturation occurs if |un|>1, and the control is given by (SignHG)/(1-λ).Another regularization process is to use the refined costminu(·)(C(x(tf))+(1-λ)∫0tf|u|2-λdt). 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 q1(tf)=0 and the cost -|q22(tf)| by a cost function of the form(6.2)Cα=α|q12(tf)|-|q22(tf)|,
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.


## 6.3. 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(6.3)dMxdt=-MxT2+ωyMz,dMydt=-MyT2-ωxMz,dMzdt=(M0-Mz)T1+ωxMy-ωyMx,
where the coordinates of the magnetization vector are (Mx,My,Mz), T1 and T2 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)=(Mx,My,Mz)/M0, 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=(ux,uy)=(2π/ωmax)(ωx,ωy), leading to |u|≤2π, while the normalized time is given by τ=ωmaxt/(2π). In these new coordinates, the system takes the form(6.4)ẋ=-Γx+uyz,ẏ=-Γy-uxz,ż=γ(1-z)+uxy-uyx,
with Γ=2π/(ωmaxT2) and γ=2π/(ωmaxT1). 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 T1=2000 ms and T2=200 ms, while for the second spin T1=T2=2500 ms.We next present some numerical computations. We first analyze the structure of the singular flow. For that purpose, we plot in Figure3 the projection of this flow onto the planes (y1,z1) and (y2,z2). 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 ((-1-z02,z0),(-1-z02,z0)) where z=z0 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.Projection of the singular flow onto the planes(y1,z1) and (y2,z2) 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 105). 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.
(a)(b)It is the same as Figure3, but the singular flow is followed for longer times to show the attracting property of the north pole for the singular flow.
(a)(b)It is the same as Figure3, but the conjugate points are computed along the singular extremals. Their positions are represented by red crosses.
(a)(b)Zoom of the results of Figure5.
(a)(b)In order to improve the contrast, a locally optimal BS-sequence solution of the maximum principle is computed as follows. The transfer timetf is fixed, and the Mayer problem is regularized by using a cost of the form(6.5)C(x(tf))+(1-λ)∫0tf|u|2-λdt.
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 λ; the transfer time is fixed to 2Tmin where Tmin 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 tf=1.1Tmin, 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 2Tmin.Evolution of the control fieldu for different values of the homotopy parameter λ, which is equal to 0, 0.6, 0.8, 0.9 from (a) to (d). The transfer time normalized to 1 is equal to 2Tmin.
(a)(b)(c)(d)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).
(a)(b)(c)It is the same as Figure8, but for a control field with a BSBS structure. Note the similarity of this result with the one of Figure 8.
(a)(b)(c)It is the same as Figure9 but for a control field of the form BSBSBS with tf=1.1Tmin.
(a)(b)(c)Figure 11
Evolution of the final contrast as a function of the transfer time.We conclude this paper by illustrating our numerical results with the first results of a contrast experiment. We consider two surfaces as displayed in Figure12 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.Experimental results on the contrast problems. The inner disk mimics the spin 1, while the outside ring mimics the spin 2. The two surfaces are separated by a thin black circle. (a) is a reference image before the application of the control field when the two spins are at the north pole of the Bloch sphere. (b) represent the contrast after a real experiment. The panel (b) corresponds to|My|, while the panel (c) to |Mz|, that is, the remaining y and z magnetization components after the optimized pulse. A color has been associated to each value of the contrast between 0 and 1, 0 and 1 corresponding, respectively, to the colors black and white. In the experiment, relaxation parameters are taken to be T1=1.3 ms, T2=50 ms for the inner sample and T1=1.3 ms, T2=200 ms for the outside sample.
(a)(b)(c)



---

*Source: 857493-2012-01-16.xml*