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

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.

1. Introduction

The dynamics of a spin 1/2 particle in nuclear magnetic resonance (NMR) is described by the Bloch equation𝑑𝑀⃗=⃗⃗⃗𝑀⃗𝑑𝑡𝑀∧𝜔−𝑅𝑀−₀, (1.1) where 𝑀⃗ is the magnetization vector of the particle, 𝜔⃗ is the radio-frequency control magnetic field, and 𝑅 is the relaxation matrix [1, 2]. Such a system is an example of bilinear system in control theory [3]𝑑𝑥=()+𝑑𝑡𝐴𝑥+𝑎^𝑚_𝑖=1𝑢_𝑖𝐵_𝑖𝑥+𝑏_𝑖, (1.2) where 𝑥∈ℝ^𝑛, (𝐴,𝐵₁,…,𝐵_𝑚) are real constant 𝑛×𝑛 matrices, and 𝑎,𝑏₁,…,𝑏_𝑚∈ℝ^𝑛. 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 𝐶∶min_𝑢(⋅)𝐶(𝑥(𝑡_𝑓)) where 𝑡_𝑓 is the transfer time. In both cases, according to the Maximum Principle [4], if the control 𝑢(⋅) is constrained to a control domain 𝑈, minimizers are solutions of the Hamiltonian equations𝑥̇=𝜕𝐻,̇𝜕𝑝𝑝=−𝜕𝐻𝜕𝑥, (1.3) where 𝐻(𝑥,𝑝,𝑢)=⟨𝑝,𝐹(𝑥,𝑢)⟩, 𝑥̇=𝐹(𝑥,𝑢) being the control system and 𝑝 the adjoint vector. The function 𝐻 depending upon 𝑢 is called the pseudo-Hamiltonian and the optimal control has to satisfy for almost every time the maximization condition [4]𝐻𝑥,𝑝,𝑢=max_𝑣∈𝑈𝐻𝑥,𝑝,𝑣. (1.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 𝜕𝐻/𝜕𝑢=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 𝑥̇=𝐹(𝑥)+𝑢𝐺(𝑥), where the control domain is |𝑢|≤1. A regular extremal is formed by bang arcs defined by 𝑢(𝑡)=sign[⟨𝑝(𝑡),𝐺(𝑥(𝑡))⟩], while the singular extremals satisfy ⟨𝑝(𝑡),𝐺(𝑥(𝑡))⟩=0. The optimal control problem reduces to find the sequence 𝐵𝑆𝐵𝑆…, where 𝐵 is a bang arc with 𝑢=±1 and 𝑆 a singular arc. For linear systems 𝑥̇=𝐴𝑥+𝐵𝑢, under the mild Kalman controllability condition rank [𝐵,…,𝐴^𝑛−1𝐵]=𝑛, 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𝐸^{𝑥0,𝑡𝑓}∶𝑢(⋅)⟼𝑥𝑡_𝑓,𝑥₀,𝑢, (1.5) where 𝑥(⋅) is the response at time 𝑡_𝑓 of the system to the control 𝑢(⋅), and 𝑥₀ 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−𝑧+𝑢₁𝑦−𝑢₂𝑥, (1.6) where the state variable 𝑞=(𝑥,𝑦,𝑧) belongs to the Bloch ball |𝑞|≤1 which is invariant for the dynamics since the parameters belong to the set Λ∶2Γ≥𝛾≥0. The control field is 𝑢=𝑢₁+𝑖𝑢₂=|𝑢|𝑒^𝑖𝜑 where the amplitude is bounded by a given 𝑚. 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) written shortly as 𝑑𝑥/𝑑𝑡=𝐹(𝑥,𝑢) where 𝑥=(𝑞₁,𝑞₂) belongs to the product of two Bloch balls. Both spins interact through the same magnetic field represented by 𝑢(⋅). The associated optimal control problem is the following. Starting from the equilibrium point 𝑥₀=((0,0,1),(0,0,1)), the goal is to reach in a given transfer time 𝑡_𝑓 the final state 𝑞₁(𝑡_𝑓)=0 (corresponding to zero magnetization of the first spin) while maximizing |𝑞₂(𝑡_𝑓)|² [11].

2. Necessary Optimality Conditions: Maximum Principle

We consider a control system of the form𝑑𝑥(𝑡)(𝑥(𝑡)(𝑡𝑑𝑡=𝐹,𝑢)), (2.1) where 𝐹∶ℝ^𝑛×ℝ^𝑚↦ℝ^𝑛 is a 𝐶^∞ mapping. The class of admissible control 𝒰 is the set {𝑢∶[0,𝑡_𝑓]→𝑈} of bounded and measurable maps where 𝑈⊂ℝ^𝑚 is the control domain. Given an initial condition 𝑥₀∈ℝ^𝑛 and 𝑢(⋅)∈𝒰, we denote by 𝑥(𝑡,𝑥₀,𝑢) the solution of the control system initiating from 𝑥₀ and defined on a subinterval of [0,𝑡_𝑓]. We fix a 𝐶^∞ terminal manifold ℳ defined by 𝑔(𝑥)=0 where 𝑔∶ℝ^𝑛→ℝ^𝑘. We will consider the two following optimal control problems.(i)Time Minimum Control Problem. Reach in minimum time from 𝑥₀ the terminal manifold ℳ.(ii)Mayer Problem. Steer 𝑥₀ to ℳ while min_{𝑢(⋅)∈𝒰}𝐶(𝑥(𝑡_𝑓)) where 𝑡_𝑓 is a fixed transfer time, and the cost function 𝐶 is a 𝐶^∞ regular function from ℝ^𝑛 into ℝ.

We denote the accessibility set at time 𝑡_𝑓 by 𝐴(𝑥₀,𝑡_𝑓)=∪_{𝑢(⋅)∈𝒰}𝑥(𝑡_𝑓,𝑥₀,𝑢).

Let us consider the Hamiltonian function𝐻𝑥,𝑝,𝑢=𝑝,𝐹(𝑥,𝑢), (2.2) where ⟨⋅,⋅⟩ is the usual scalar product in ℝ^𝑛, and 𝑝∈ℝ^𝑛⧵{0} is the adjoint vector. 𝐻 is called the pseudo-Hamiltonian, and we introduce the maximized Hamiltonian 𝑀(𝑥,𝑝)=max_𝑣∈𝑈𝐻(𝑥,𝑝,𝑣).

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

Let 𝑢^∗(⋅) be an admissible control whose corresponding trajectory 𝑥^∗(⋅)=𝑥(⋅,𝑢^∗,𝑥₀) is optimal, then there exists a vector function 𝑝^∗(⋅)∈ℝ^𝑛⧵{0} such that the following conditions are satisfied:𝑥̇^∗=𝜕𝐻𝑥𝜕𝑝^∗,𝑝^∗,𝑢^∗,𝑝^∗̇=−𝜕𝐻𝑥𝜕𝑥^∗,𝑝^∗,𝑢^∗, (2.3) with the maximality condition𝐻𝑥^∗(𝑡),𝑝^∗(𝑡),𝑢^∗(𝑡)=𝑀𝑥^∗(𝑡),𝑝^∗(𝑡), (2.4) where 𝑀 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 𝑡_𝑓, 𝑔𝑥^∗𝑡_𝑓=0,𝑝^∗𝑡_𝑓⟂𝜕𝑔𝑥𝜕𝑥^∗𝑡_𝑓, (2.5)(ii)Mayer problem: at the final time 𝑡_𝑓, 𝑔𝑥^∗𝑡_𝑓𝑝=0,_∗𝑡_𝑓=𝑝₀𝜕𝐶𝑥𝜕𝑥^∗𝑡_𝑓+^𝑘_𝑖=1𝛿_𝑖𝜕𝑔_𝑖𝑥𝜕𝑥^∗𝑡_𝑓,𝛿=𝛿₁,…,𝛿_𝑘∈ℝ^𝑘,𝑝₀≤0. (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 (𝑥(⋅),𝑝(⋅),𝑢(⋅)) solutions of this system are called extremals. Geometrically they correspond to necessary conditions for the terminal point 𝑥(𝑡_𝑓) to be in the boundary of the accessibility set 𝐴(𝑥₀,𝑡_𝑓).

The respective boundary conditions define the so-called BC-extremals:(i)time minimum case: 𝑝^∗(𝑡_𝑓)⟂ℳ,(ii)Mayer problem: let us fix 𝑚∈ℝ and introduce the manifold ℳ_𝑚={𝑥,𝑔(𝑥)=0,𝐶(𝑥)=𝑚}. One must have for the optimal solution 𝑝^∗(𝑡_𝑓)⟂ℳ_𝑚, where 𝑚 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 (𝜕𝐻/𝜕𝑢)(𝑥,𝑝,𝑢)=0. The following results are well known.

Definition 2.1. Relaxing 𝑢∈𝑈, a singular trajectory of the system 𝑥̇=𝐹(𝑥,𝑢) on [0,𝑡_𝑓] is a control trajectory pair (𝑥,𝑢) such that the input-output mapping 𝐸^{𝑥0,𝑡𝑓}∶𝑢(⋅)∈𝐿^∞0,𝑡_𝑓⟼𝑥𝑡_𝑓,𝑥₀,𝑢 (2.7) is singular.

Proposition 2.2. One then deduces the following results.(1)If (𝑥,𝑢) is not singular on [0,𝑡_𝑓], then the input-output mapping is open at 𝑢(⋅) for the 𝐿^∞[0,𝑡_𝑓]-norm topology.(2)If (𝑥,𝑢) is singular on [0,𝑡_𝑓], then there exists 𝑝(⋅)∈ℝ^𝑛⧵{0} such that (𝑥,𝑝,𝑢) satisfies a.e. on [0,𝑡_𝑓] the weak maximum principle𝑥̇=𝜕𝐻,̇𝜕𝑝𝑝=−𝜕𝐻,𝜕𝑥𝜕𝐻𝜕𝑢=0. (2.8)

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 which 𝐹(𝑥,𝑢)=𝐹₀(𝑥)+𝑢₁𝐹₁(𝑥)+𝑢₂𝐹₂(𝑥), where the 𝐹_𝑖 are 𝐶^∞-vector fields and 𝑢=(𝑢₁,𝑢₂), with the control domain 𝑈∶|𝑢|≤1. Let us denote by 𝐻_𝑖=⟨𝑝,𝐹_𝑖(𝑥)⟩ the Hamiltonian lifts, then the pseudo-Hamiltonian 𝐻(𝑥,𝑝,𝑢) takes the form 𝐻₀+∑²_𝑖=1𝑢_𝑖𝐻_𝑖. The maximality condition (2.4) leads to the following parameterization of the extremal controls:𝑢₁=𝐻₁𝐻²₁+𝐻²₂,𝑢₂=𝐻₂𝐻²₁+𝐻²₂, (2.9) outside the switching surface Σ∶𝐻₁=𝐻₂=0. Plugging such 𝑢 into the Hamiltonian gives the true Hamiltonian 𝐻_𝑛=𝐻₀+(𝐻²₁+𝐻²₂)^1/2. The smooth solutions of the corresponding Hamiltonian field are called extremals of order zero.

The relation 𝑢²₁+𝑢²₂=1 leads to the following interpretation.

Proposition 2.3. Extremals of order zero correspond to the singularity of the input-output mapping 𝐸^{𝑥0,𝑡𝑓}∶𝑢∈𝐿^∞0,𝑡_𝑓∩{|𝑢|=1}⟼𝑥𝑡_𝑓,𝑥₀,𝑢. (2.10)

2.4. Single-Input Affine Case

The system is of the form 𝐹(𝑥,𝑢)=𝐹(𝑥)+𝑢𝐺(𝑥) and |𝑢|≤1. Applying the maximality condition, one gets two types of extremals.(i)Regular extremals: the control is given by 𝑢(𝑡)=sign[𝐻_𝐺(𝑧(𝑡))], 𝑧=(𝑥,𝑝), and 𝐻_𝐺=⟨𝑝,𝐺⟩. If 𝑢=±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 in 𝑢, the condition (2.4) leads in the singular case to the condition 𝐻_𝐺=0, identically. Differentiating twice with respect to time, we get the relations 𝐻_𝐺,𝐻_𝐹=𝐻_𝐺,𝐻_𝐹,𝐻_𝐹+𝑢𝐻_𝐺,𝐻_𝐹,𝐻_𝐺=0, (2.11) and from this second condition, one derives when the denominator is not vanishing the corresponding singular control𝑢_𝑆𝐻=−^𝐺,𝐻_𝐹,𝐻_𝐹𝐻_𝐺,𝐻_𝐹,𝐻_𝐺. (2.12)

Plugging such 𝑢_𝑆 into the pseudo-Hamiltonian defines the singular flow solution of the Hamiltonian vector field denoted as 𝐻_𝑆 and starting at 𝑡=0 from the two constraints 𝐻_𝐺={𝐻_𝐺,𝐻_𝐹}=0 and |𝑢_𝑆|≤1. They have the following interpretation.

Proposition 2.4. Singular extremals correspond to the singularities of the input-output mapping 𝐸^{𝑥0,𝑡𝑓}∶𝑢∈𝐿^∞0,𝑡_𝑓⟼𝑥𝑡_𝑓,𝑥₀,𝑢. (2.13)

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𝜕²𝐻𝜕𝑢²≤0, (2.14) and if it is strict, it is called the strict Legendre-Clebsch condition. In the singular case, one has 𝜕²𝐻/𝜕𝑢²=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𝜕𝑑𝜕𝑢𝑑𝑡²²𝜕𝐻=𝐻𝜕𝑢^𝐺,𝐻_𝐺,𝐻_𝐹≤0. (2.15)

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 𝑢²₁+𝑢²₂=1, and one can consider the extension of the control system by setting 𝑢₁=sin𝛼, 𝑢₂=cos𝛼, and 𝛼̇=𝑣 where (𝑥,𝛼) is the extended state space and 𝑣 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 𝐺 is nonzero, there exists a coordinate system (𝑥₁,…,𝑥_𝑛) on an open set 𝑉 such that 𝐺=𝜕/𝜕𝑥_𝑛. Denoting 𝑥=(𝑥′,𝑥_𝑛), the system splits into𝑥̇′=𝐹𝑥^′,𝑥_𝑛,𝑥̇_𝑛=𝐹^′₀𝑥^′+𝑢, (2.16) where the system 𝐹′ defined on an open set 𝑉′ with 𝑥_𝑛, 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 𝐻′(𝑥′,𝑝′,𝑥_𝑛)=⟨𝑝′,𝐹′(𝑥′,𝑥_𝑛)⟩, one has𝑑𝑑𝑡𝜕𝐻=𝐻𝜕𝑢^𝐺,𝐻_𝐹=−𝜕𝐻′𝜕𝑥_𝑛,𝜕𝑑𝜕𝑢𝑑𝑡²²𝜕𝐻=𝐻𝜕𝑢^𝐺,𝐻_𝐹,𝐻_𝐺𝜕=−²𝐻′𝜕𝑥²_𝑛. (2.17) 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 is 𝑥̇=𝐹(𝑥,𝑢), 𝑢∈𝑈, and 𝑥∈𝑋, and the pseudo-Hamiltonian can be written as 𝐻(𝑧,𝑢)=⟨𝑝,𝐹(𝑥,𝑢)⟩, where 𝑧=(𝑥,𝑝). One can consider a reference singular extremal (𝑧,𝑢) solution on [0,𝑡_𝑓] of𝑥̇=𝜕𝐻,̇𝜕𝑝𝑥,𝑝,𝑢𝑝=−𝜕𝐻,𝜕𝑥𝑥,𝑝,𝑢𝜕𝐻𝜕𝑢=0. (3.1)

Our first assumption is the strong-Legendre condition:(A1) The quadratic form 𝜕²𝐻/𝜕𝑢² is negative definite along the reference solution.

Using the implicit function theorem, the extremal control is then locally defined as a smooth function 𝑢(𝑧), and plugging this function into 𝐻 defines a smooth true Hamiltonian, still denoted as 𝐻, and the extremal is solution of 𝑧̇=𝐻⃗(𝑧) with initial condition 𝑧₀=(𝑥₀,𝑝₀).

3.1. The Geometric Concept of Conjugate Point

Definition 3.1. Let 𝑧=(𝑥,𝑝) be the reference extremal defined on [0,𝑡_𝑓]. The variational equation 𝛿𝑧̇=𝑑𝐻⃗(𝑧(𝑡))𝛿𝑧 (3.2) is called the Jacobi equation. A Jacobi field is a nontrivial solution 𝛿𝑧=(𝛿𝑥,𝛿𝑝), and it is said to be vertical at time 𝑡 if 𝛿𝑥(𝑡)=𝑑Π_𝑧(𝑡)𝛿𝑧(𝑡)=0 where Π∶(𝑥,𝑝)↦𝑥 is the canonical projection.

The following geometric result is crucial [13].

Proposition 3.2. Let 𝐿₀ be the fiber 𝑇^∗_𝑥0𝑋, and let 𝐿_𝑡=exp[𝑡𝐻⃗(𝐿₀)] be its image by the one-parameter subgroup generated by 𝐻, then 𝐿_𝑡 is a Lagrangian manifold whose tangent space at 𝑧(𝑡) is spanned by the Jacobi fields vertical at 𝑡=0. Moreover, the rank of the restriction to 𝐿_𝑡 of Π is at most 𝑛−1.

This leads to the following definition.

Definition 3.3. We define the exponential mapping exp_𝑥0,𝑡𝑝₀=Π𝑧𝑡,𝑥₀,𝑝₀ (3.3) as the projection on the state space of the integral curve of 𝐻⃗ with initial condition 𝑧₀=(𝑥₀,𝑝₀), where 𝑝₀ can be restricted by linearity to the sphere |𝑝₀|=1. If 𝑧=(𝑥,𝑝) is the reference extremal, a time 𝑡_𝑐>0 is said to be geometrically conjugate to 0 if the mapping 𝑝₀↦exp_𝑥0,𝑡 is not of rank (𝑛−1) at 𝑡=𝑡_𝑐, and the associated point 𝑥(𝑡_𝑐) is said to be geometrically conjugate to 𝑥₀. We denote by 𝑡_1𝑐 the first conjugate time and 𝐶(𝑥₀) is the conjugate locus formed by the set of first conjugate points.

An algorithm can be deduced.

Testing Conjugacy
Let 𝑧(𝑡)=(𝑥(𝑡),𝑝(𝑡)) be the reference extremal, and consider the vector space of dimension 𝑛−1 generated by the Jacobi fields 𝛿𝑧_𝑖=(𝛿𝑥_𝑖,𝛿𝑝_𝑖), 𝑖=1,…,𝑛−1 vertical at 𝑡=0 and such that 𝛿𝑝_𝑖(0) is orthogonal to 𝑝₀. At a conjugate time 𝑡_𝑐, one has rank𝛿𝑥₁𝑡_𝑐,…,𝛿𝑥_𝑛−1𝑡_𝑐<𝑛−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 𝑧(𝑡), 𝑡∈[0,𝑡_𝑓] associated to 𝑢(⋅) is by definition a singularity of the input-output mapping 𝐸^{𝑥0,𝑡𝑓}. One assumes that on each subinterval [𝑡₀,𝑡₁], 0<𝑡₀<𝑡₁≤𝑡_𝑓, the singularity of 𝐸 on 𝑢|_{[𝑡0,𝑡1]} is of codimension one.(A3) We are in the nonexceptional case where 𝐻≠0 along the reference extremal.

Theorem 3.4. Under our assumptions, the first geometric conjugate time 𝑡_1𝑐 is the first time t such that if 𝑡<𝑡_1𝑐, the input-output mapping 𝐸^𝑥0,𝑡 is not open at 𝑢|_[0,𝑡], and if 𝑡>𝑡_1𝑐, it becomes open at 𝑢|_[0,𝑡] for the 𝐿^∞-norm topology.

Corollary 3.5. Consider the control system 𝐹(𝑥,𝑢)=𝐹₀(𝑥)+𝑢₁𝐹₁(𝑥)+𝑢₂𝐹₂(𝑥), |𝑢|≤1, and let (𝑥,𝑝,𝑢) be a reference extremal of order zero solution of 𝐻₀+(𝐻²₁+𝐻²₂)^1/2 such that assumptions (A2), (A3) are satisfied, then the first time 𝑡_1𝑐 such that the input-output mapping becomes open along 𝑥(⋅) 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 system 𝑥̇=𝐹(𝑥)+𝑢𝐺(𝑥), |𝑢|≤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 condition 𝑥(0)=𝑥₀, but with variable final condition 𝑥(𝑡_𝑓)∈ℳ_𝑚. From the maximum principle, the transversality condition defined a Lagrangian manifold 𝑧(𝑡_𝑓)=(𝑥(𝑡_𝑓),𝑝(𝑡_𝑓))∈ℳ^⟂_𝑚 which again defines along a singular extremal a train of Lagrangian manifolds 𝐿_𝑡, integrating backwards 𝑡≤𝑡_𝑓. Optimality can be deduced from the geometric concept of focal points corresponding to the singularities of (𝐿_𝑡,Π).

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 trajectory 𝑡↦𝑥(𝑡), 𝑡∈[0,𝑡_𝑓] 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,𝑡_𝑓], then we can embed locally the reference extremal into a central field ℱ formed by all extremal trajectories starting from 𝑥₀=𝑥(0). One restricts the adjoint vector with the normalization 𝐻=1 and the assumption (A3). One introduces𝐿=(𝐻=1)∩∪_𝑡≥0𝐿_𝑡, (3.5) where 𝐿_𝑡 is the train of Lagrangian manifolds generated by the Hamiltonian flow 𝐻⃗, with 𝐿₀=𝑇^∗_𝑥0𝑋 being the initial condition. By construction, ℱ is the projection of 𝐿 on the state space 𝑋.

Proposition 3.6. Excluding 𝑥₀, there exists an open neighborhood 𝑊 of the reference trajectory and two smooth mappings 𝑉∶𝑊↦ℝ, 𝑢∶𝑊→𝑈 such that for each (𝑥,𝑢)∈𝑊×𝑈, the maximization condition 𝐻(𝑥,𝑑𝑉(𝑥),𝑢(𝑥))≥𝐻(𝑥,𝑑𝑉(𝑥),𝑢). (3.6) The reference trajectory is optimal with respect to all smooth trajectories of the system, with the same extremities and contained in 𝑊. The Lagrangian manifold 𝐿 is a graph, and 𝑉 is the generating mapping {𝑥∈𝑊,𝑝=𝜕𝑉/𝜕𝑥}.

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 𝑞̇=𝐹+𝑢𝐺,𝑑𝑦𝑑𝑡=−Γ𝑦−𝑢𝑧,𝑑𝑧()|𝑢|𝑑𝑡=𝛾1−𝑧+𝑢𝑦,≤𝑚. (4.1) 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 𝑞₀=(𝑦,𝑧), find the time minimal solution in a small neighborhood 𝑉 of 𝑞₀, according to the Lie algebraic structure of the system at 𝑞₀;(2)glue together the local solutions to get the global time optimal solution. More precisely, we will describe the global synthesis to steer 𝑁=(0,0,1) to any point of the Bloch ball. By symmetry, one can only consider the domain 𝑦≤0. This amounts to compute the switching locus Σ formed by points where the optimal solutions switch and the separating locus 𝐿 where there exist two different optimal solutions which intersect.

4.1. Computations of Lie Brackets and Singular Trajectories Properties

One has𝜕𝐹=−Γ𝑦()𝜕𝜕𝑦+𝛾1−𝑧,𝜕𝜕𝑧𝐺=−𝑧𝜕𝜕𝑦+𝑦𝜕𝑧 (4.2) and the Lie brackets up to order three are[]=𝜕𝐺,𝐹−𝛾+𝛿𝑧𝜕𝜕𝑦+𝛿𝑦𝜕𝑧, (4.3) where 𝛿=𝛾−Γ,[[𝐺,𝐹],𝐹]=𝛾𝛾−2Γ−𝛿²𝑧𝜕𝜕𝑦+𝛿²𝑦𝜕,]]𝜕𝜕𝑧[[𝐺,𝐹,𝐺=2𝛿𝑦+𝜕𝜕𝑦𝛾−2𝛿𝑧.𝜕𝑧 (4.4) The singular trajectories are contained in the set 𝑆 defined by det(𝐺,[𝐺,𝐹])=0 which leads to 𝑦(−2𝛿𝑧+𝛾)=0. Hence, the singular locus is the union of the two singular lines: the line 𝑦=0 corresponding to the axis of revolution and the horizontal axis 𝑧₀=𝛾/(2𝛿). Eliminating 𝑝, the singular control is given by 𝐷′+𝑢_𝑠𝐷=0 where 𝐷′=det(𝐺,[[𝐺,𝐹],𝐹]) and 𝐷=det(𝐺,[[𝐺,𝐹],𝐺]). Computing, one gets(i)for 𝑦=0, this simplifies into 𝐷′=0 and 𝐷=−𝑧(𝛾−2𝛿𝑧). Hence, the singular control is zero, and the dynamics is governed by the equation 𝑧̇=𝛾(1−𝑧),(ii)for 𝑧₀=𝛾/(2𝛿), we get 𝐷′=−𝑦𝛾(2𝛿−𝛾) and 𝐷=−2𝛿𝑦², which gives 𝑢_𝑆=𝛾𝛾−2𝛿2𝛿𝑦, (4.5) where 𝛾−2𝛿=2Γ−𝛾≥0. The flow on the horizontal direction is given by𝑦̇=−Γ𝑦−𝛾²2Γ−𝛾4𝛿²𝑦. (4.6)

4.2. Parameter Conditions

The interesting case is when 𝛿<0 and |𝛾/2𝛿|<1. This leads toΓ>3𝛾2. (4.7) In this case, 𝑢_𝑆→+∞ when 𝑦→0⁻. Moreover, we have the following estimate. Setting 𝛼=𝛾²(2Γ−𝛾)/(4𝛾²), we get near 0 that 𝑦̇≃−𝛼/𝑦, and the trajectory reaching zero in 𝑦<0 is of order 𝑂(√𝑡), the singular control being of order 𝑂(1/√𝑡). In particular, one deduces the following.

Lemma 4.1. The singular control in the neighborhood of the point 𝑦=0, along the horizontal direction is of order 𝑂(1/√𝑡) and is in the 𝐿¹ but not 𝐿² 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 have 𝐻≥0 from the maximum principle, and from the Legendre-Clebsch condition (2.15), we deduce that {𝐻_𝐺,{𝐻_𝐺,𝐻_𝐹}}≤0. More precisely, if we introduce 𝐷^′′=det(𝐺,𝐹)=𝛾𝑧(𝑧−1)+Γ𝑦², one arrives at(i)hyperbolic case 𝐷𝐷^′′>0, the singular direction is small time minimizing;(ii)elliptic case 𝐷𝐷^′′<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 𝑧₀<𝑧<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≤𝑧<𝑧₀ time maximizing singular arc.

In order to complete the optimality analysis, one introduces for 2D-system the following clock-form 𝜔. The collinear set 𝐶 is defined by det(𝐹,𝐺)=0. From a straightforward computation, one deduces that it will form an oval curve defined by 𝛾𝑧(𝑧−1)+Γ𝑦²=0, which shrinks into a point if 𝛾=0. Outside this set, 𝜔=𝑝𝑑𝑥 is defined by 𝜔(𝐹)=1 and 𝜔(𝐺)=0, the sign of 𝑑𝜔 being given by 𝑦(𝛾−2𝛿𝑧). This form allows to compare the time to travel different trajectories not crossing 𝐶 by using the Stokes theorem. We have also 𝑑𝜔=0 on the singular set 𝑆. Using the form 𝜔, one deduces the following lemma.

Lemma 4.2. Assume 𝑚=+∞, 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 Σ′={𝑤=(𝑞,𝑝)/⟨𝑝,𝐺(𝑞)⟩=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 subset 𝐻_𝐺={𝐻_𝐺,𝐻_𝐹}=0, and we denote by 𝛿₊ and 𝛿₋ bang arcs 𝑢=±1 for the system 𝑑𝑥/𝑑𝑡=𝐹+𝑢𝐺, |𝑢|≤1. We denote by 𝛿_𝑆 a singular arc and by 𝛿₁𝛿₂ an arc 𝛿₁ followed by an arc 𝛿₂. We have the following.

Ordinary switching arc: it is a time 𝑡 such that a bang arc switches with the condition Φ(𝑡)=0 and Φ̇(𝑡)={𝐻_𝐺,𝐻_𝐹}(𝑤(𝑡))≠0 where Φ(𝑡)=𝐻_𝐺(𝑤(𝑡)) is the switching function evaluated along an extremal arc 𝑤(⋅). According to the maximum principle, near Σ′, the extremal is of the form 𝛿₊𝛿₋ (resp., 𝛿₋𝛿₊) if Φ̇(𝑡)>0 (resp., Φ̇(𝑡)<0).

Fold case: denoting Φ̈_±(𝑡)={{𝐻_𝐺,𝐻_𝐹},𝐻_𝐹}±{{𝐻_𝐺,𝐻_𝐹},𝐻_𝐺}(𝑧(𝑡)) the second-order derivative along a bang arc ±1, the fold case is the situation where Φ_±(𝑡)=Φ̇_±(𝑡)=0 and Φ̈_±≠0. We have three cases.(i)Hyperbolic Case. at the contact point with Σ^′, Φ̇₊(𝑡)>0 and Φ̇₋(𝑡)<0. The connection with a singular extremal with |𝑢_𝑆|<1 and satisfying the strong generalized condition is possible. The extremals near such a point are of the form 𝛿_±𝛿_𝑆𝛿_±.(ii)Elliptic Case. at the contact point with Σ′, Φ̇₊(𝑡)<0 and Φ̇₋(𝑡)>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, Φ̇₊(𝑡) and Φ̇₋(𝑡) 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 Φ̇₊(𝑡)>0 or 𝛿₋𝛿₊𝛿₋ if Φ̇₋(𝑡)<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 control 𝑢_𝑆 saturates at a point 𝐵. This situation was analyzed in the literature, see, for instance, [15]. At the point 𝐵, there is a birth of a switching curve which can be estimated using the techniques of local models [13]. 𝐵 is identified to 0, the singular arc is normalized to the 𝑥-axis, and the model can be written as𝑥̇=1−𝑦²,𝑦̇=−(𝑥+1)+𝑢,|𝑢|≤1. (4.8) The singular control is 𝑢_𝑆=1+𝑥, and it is not admissible for 𝑥>0. Near 0, we have optimal policies of the form 𝛿₊𝛿_𝑆𝛿₊𝛿₋.

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 𝑚 is large enough such that the bang arc 𝑢=𝑚 starting from the north pole intersects the singular arc 𝑧₀=𝛾/(2𝛿), Then for the problem of reaching zero magnetization, the horizontal singular arc is saturating at a point denoted as 𝐵 and from the analysis of [7]; they are local optimal policies of the form 𝛿_𝑚𝛿_𝑆ℎ𝛿_𝑚𝛿_𝑆𝑣 where 𝛿_𝑚 is an arc associated to 𝑢=𝑚, and 𝛿_𝑆ℎ and 𝛿_𝑆𝑣 are, respectively, singular horizontal and vertical arcs. Note in particular that the horizontal singular arc is not followed up to the saturating point 𝐵. 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 𝐵, 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, 𝑡 to be two consecutive switching times on an arc 𝛿₊ or 𝛿₋ where the control is 𝑢=𝜀=±1. We must have𝑝(0),𝐺𝑞(0)=𝑝(𝑡),𝐺𝑞(𝑡)=0. (4.9) We denote by 𝑉(⋅) the solution of the variational equation 𝑉̇=((𝜕𝐹/𝜕𝑞)+𝜀(𝜕𝐺/𝜕𝑞))𝑉 such that 𝑉(𝑡)=𝐺(𝑞(𝑡)), where the equation is integrated backwards from time 𝑡 to 0. By construction ⟨𝑝(0),𝑉(0)⟩=0 (since ⟨𝑝(𝑡),𝑉(𝑡)⟩=0 and the derivative is zero) and hence at time 𝑡=0, 𝑝(0) is orthogonal to 𝐺(𝑞(0)) and to 𝑉(0). Therefore, 𝑉(0) and 𝐺(𝑞(0)) are collinear. We denote by 𝜃(𝑡) the angle between 𝐺(𝑞(0)) and 𝑉(𝑡) measured counterclockwise. One deduces that switching occurs when 𝜃(𝑡)=0[mod(𝜋)]. Moreover, 𝜃̇=0 on the singular set det(𝐺,[𝐺,𝐹])=0.

We have by definition𝑉(𝑡)=𝑒^{−𝑡ad(𝐹+𝜀𝐺)}𝐺𝑞(𝑡), (4.10) and in the analytic case, the ad-formula gives, for small 𝑡,𝑉(0)=_𝑛≥0(−𝑡)^𝑛𝑛!ad^𝑛(𝐹+𝜀𝐺)𝐺𝑞(𝑡). (4.11) The computation can be made explicit in the bilinear case 𝐹(𝑞)=𝐴𝑞+𝑎, 𝐺(𝑞)=𝐵𝑞. The system can be lifted into an invariant system on the semidirect product Lie group GL(2,ℝ)×_𝑆ℝ² identified with the set of matrices of GL(3,ℝ):()10𝑔𝑤,𝑔∈GL2,ℝ,𝑤∈ℝ², (4.12) acting on the subspace of vector in ℝ³∶¹_𝑞. Lie brackets of two affine vector fields (𝐴′,𝑎′), (𝐵′,𝑏′) being defined by𝐴^′,𝑎^′,𝐵^′,𝑏^′=𝐴^′,𝐵^′,𝐴^′𝑏^′−𝐵^′𝑎^′, (4.13) the computation of the exponential exp[−𝑡ad(𝐹+𝜀𝐺)] is a standard exercise in linear algebra which amounts to compute a Jordan normal form of ad(𝐹+𝜀𝐺), 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 Figure 1. We have assumed that Γ>3𝛾/2 and 𝑚≫Γ. Using the symmetry of revolution with respect to the 𝑧-axis, one can restrict the study to the domain 𝑦≤0. The north pole being a singular point for the optimal control, the first applied control is taken as 𝑢=𝑚 and will form a boundary arc of the accessibility domain. The collinear set is a parabola connecting 𝑂 to the north pole. The main feature of the synthesis is the SiSi singularity. The switching locus Σ is composed of the arc (denoted Σ₁) starting from the north pole and reaching the horizontal singular arc at 𝐴, the horizontal line segment Σ₂ between 𝐴 and 𝐵, the singular control saturating at 𝐵, the switching locus Σ₃ due to the saturating phenomenon ending on the vertical axis at 𝐷, and the vertical singular line segment denoted as Σ₄ between 𝐷 and 𝑂. The bang arc 𝑢=−𝑚 starting from 𝐴 separates the synthesis in two domains, one with a bang-bang policy 𝛿_𝑚𝛿_−𝑚 and the other where the optimal policy contains a nontrivial horizontal singular line. Due to the symmetry of revolution, the separating locus 𝐿 reduces to the 𝑧-axis.

Figure 1 

Global synthesis from the north pole. The switching locus Σ is formed by the union of the surfaces Σ𝑖. Due to the symmetry of revolution, the cut locus corresponds to the 𝑧-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 (𝑦,𝑧). Experimental parameters are taken to be 𝑇1=740ms, 𝑇2=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𝜋))𝑡.

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) where 𝑥=(𝑞₁,𝑞₂), 𝑞_𝑖=(𝑦_𝑖,𝑧_𝑖) representing the state of each spin 1/2 particle. From Bloch equation, we get𝐹=²_𝑖=1−𝑦_𝑖Γ_𝑖𝜕𝜕𝑦_𝑖+𝛾_𝑖1−𝑧_𝑖𝜕𝜕𝑧_𝑖,𝐺=²_𝑖=1−𝑧_𝑖𝜕𝜕𝑦_𝑖+𝑦_𝑖𝜕𝜕𝑧_𝑖, (5.2) where Λ_𝑖=(Γ_𝑖,𝛾_𝑖) are the relaxation parameters characterizing each spin.

5.1. BC Singular Extremals

According to the general computations of Section 2, a singular extremal is contained in the constrained setΣ^′′∶𝑝,𝐺(𝑥)=𝑝,[𝐺,𝐹](𝑥)=0, (5.3) while the singular control is defined by𝑢_𝑆]𝑥=−𝑝,[[𝐺,𝐹,𝐹]()𝑝,[[𝐺,𝐹],𝐺](𝑥), (5.4) the Lie brackets being computed in Section 2.

In order to be optimal, the singular extremals have to satisfy the generalized Legendre-Clebsch condition ⟨𝑝,[𝐺,[𝐺,𝐹]]⟩≤0 and the transversality condition at the final time. Decomposing the adjoint vector in the form 𝑝=(𝑝₁,𝑝₂) where 𝑝_𝑖 is associated to the spin 𝑖, one has𝑝₂𝑡_𝑓=−2𝑝₀𝑞₂𝑡_𝑓,𝑝₀≤0, (5.5) and if 𝑝₀≠0, it can be normalized by homogeneity to 𝑝₀=−1/2. The case 𝑝₀=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 𝑡_𝑓 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 constraint 𝑀=max_𝑢(⋅)[𝐻_𝐹+𝑢𝐻_𝐺]=0, which gives ⟨𝑝,𝐹(𝑥)⟩=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) where([],],([],],𝐷=det𝐹,𝐺,𝐺,𝐹[[𝐺,𝐹,𝐺])𝐷′=det𝐹,𝐺,𝐺,𝐹[[𝐺,𝐹,𝐹]) (5.7) with the corresponding vector fields defined by𝑑𝑥(𝑥)−𝐷𝑑𝑡=𝐹^′(𝑥)𝐷(𝑥)𝐺(𝑥). (5.8) Observe that in the general case, a similar computation shows that the singular extremals are solutions of an equation of the form𝑑𝑥(𝑥)−𝐷𝑑𝑡=𝐹^′(𝑥,𝜆)𝐷()𝐺(𝑥)𝑥,𝜆, (5.9) 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) using the time reparameterization 𝑑𝑠=𝑑𝑡/𝐷(𝑥(𝑡)). According to the formula (5.4), the singular control can explode near 𝐷=0, in particular saturating the constraint |𝑢|≤𝑚. This is connected to the bridge phenomenon described in Section 3. Some of the singular extremals can cross the singular surface 𝐷=0, provided 𝐷′=0.

5.4. Feedback Classification and the Contrast Problem

An important object to analyze our problem is the action of the feedback group 𝐺_𝑓 on the set of systems. One will restrict our presentation to the single-input affine case, and we denote by 𝒞=(𝐹,𝐺) the set of all such (smooth) systems on a state space 𝑋≃ℝ^𝑛.

Definition 5.2. Let (𝐹,𝐺), (𝐹^′,𝐺^′) be two elements of 𝒞. They are called feedback equivalent if there exists a smooth diffeomorphism 𝜑 of ℝ^𝑛 on a feedback 𝑢=𝛼(𝑥)+𝑢𝛽(𝑥) where 𝛼, 𝛽 are smooth, 𝛽 is invertible such that (i)(),(),𝐹′=𝜑∗𝐹+𝜑∗𝐺⋅𝛼ii𝐺′=𝜑∗𝐺⋅𝛽 (5.11) where 𝜑∗𝑍 is the vector field image given by 𝜑∗𝑍=(𝜕𝜑⁻¹/𝜕𝑥)(𝑍∘𝜑). This action defines a group structure on the set of triplets (𝜑,𝛼,𝛽); this group is denoted as 𝒢.

Definition 5.3. Let (𝐹,𝐺)∈𝒞 and let 𝜆 be the map which associates to (𝐹,𝐺) the constrained differential equation whose solutions are singular extremals. The constraint is given by Σ^′′∶𝐻_𝐺={𝐻_𝐺,𝐻_𝐹}=0, and the equation is 𝐹_𝑆=𝐹+𝑢_𝑆𝐺 where 𝑢_𝑆=−{{𝐻_𝐺,𝐻_𝐹},𝐻_𝐹}/{{𝐻_𝐺,𝐻_𝐹},𝐻_𝐺}. We define the action of (𝜑,𝛼,𝛽) on (Σ,𝐹_𝑆) by the action of the symplectic change of coordinates 𝜑⃗∶𝑥=𝜑(𝑋),𝑝=𝑃𝜕𝜑⁻¹𝜕𝑋, (5.12) a feedback acting trivially.

We have the following theorem [23].

Theorem 5.4. The following diagram is commutative: (𝒞)𝐺𝒞⟶𝜆_𝑓↓↓𝐺_𝑓𝒞⟶𝜆(𝒞). (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 (Γ₁,𝛾₁,Γ₂,𝛾₂) to geometric invariants of the singular flow and optimality properties. The research program is the following.

5.5. Classification of the Distribution 𝐺

A feedback 𝑢=𝛼(𝑥)+𝛽(𝑥)𝑣 acts on the control direction 𝐺 as 𝐺→𝛽𝐺, which corresponds to the classification of the one-dimensional distribution 𝑥→Span𝐺(𝑥).

5.6. Collinear Set

The collinear set 𝐶 is the set of points where 𝐹 and 𝐺 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) restricted to the surfaceΣ^′′∶𝐻_𝐺=𝐻_𝐺,𝐻_𝐹=0, (5.15) where 𝑢_𝑆 is given by the formulae (2.14).

This set of equations defines a Hamiltonian vector field on Σ^′′, restricting the standard symplectic form 𝜔=𝑑𝑥∧𝑑𝑝.

We introduce the two Poisson brackets:𝒟=𝐻_𝐺,𝐻_𝐹,𝐻_𝐺,𝒟^′=𝐻_𝐺,𝐻_𝐹,𝐻_𝐹, (5.16) and the differential equation (5.14) can be desingularized using the time reparameterization 𝑑𝑠=𝑑𝑡/(𝒟(𝑥(𝑡),𝑝(𝑡))). One gets the equations𝑑𝑥𝑑𝑠=𝒟𝐹−𝒟^′𝐺,𝑑𝑝𝒟𝑑𝑠=−𝑝𝜕𝐹𝜕𝑥−𝒟^′𝜕𝐺,𝜕𝑥 (5.17) 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 a C⁰-synthesis in the limit case 𝑚≫Γ,𝛾. This result is based on [9, 26].

6.1. Synthesis for a Reference Singular Arc

We consider the control system 𝑑𝑥/𝑑𝑡=𝐹(𝑥)+𝑢𝐺(𝑥), |𝑢|≤1, and let 𝑢_𝑆 be a smooth singular control such that |𝑢_𝑆|<1 with corresponding trajectory 𝑥(𝑡). We have the following.

Theorem 6.1. Under generic assumptions, the first conjugate time 𝑡_1𝑐 is the time along the singular arc such that a policy 𝛿_±𝛿_𝑆𝛿_± is time extremizing in a C⁰-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 𝑢=+1 or 𝑢=−1 (the maximum bound 𝑚 being normalized to 1), and by symmetry of revolution, one can assume that 𝑢=+1. At the final time 𝑡_𝑓, one must have 𝑝₂(𝑡_𝑓)=𝑞₂(𝑡_𝑓), and moreover 𝐻_𝐺(𝑡_𝑓)=0. The final point concentrates a switching. Hence, one gets the following.

Lemma 6.2. The simplest possible sequence is of the form 𝐵𝑆, where 𝐵 is a bang arc and 𝑆 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 𝑥̇=𝐴𝑥+𝐵𝑢, and the cost can be written as ∫^𝑡𝑓₀^𝑡𝑢𝑅𝑢𝑑𝑡. The regular case corresponds to 𝑅>0, and the singular case called cheap control is the case where 𝑅≥0, and some control components (𝑢₁,𝑢₂,…,𝑢_𝑝) are not penalized. A regularization process consists of adding quadratic penalties in the cost 𝜀(∑^𝑝_𝑖=1𝑢²_𝑖) 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_𝑢(⋅)(𝐶(𝑥(𝑡_𝑓))+(1−𝜆)∫^𝑡𝑓₀𝑢(𝑡)²𝑑𝑡), 𝜆∈[0,1]. For the regularized problem, one considers only normal extremals associated to the pseudo-Hamiltonian𝐻_𝜆=𝐻_𝐹+𝑢𝐻_𝐺−12(1−𝜆)𝑢². (6.1) The control is computed in the normal case by solving 𝜕𝐻_𝜆/𝜕𝑢=0, 𝑢_𝑛=𝐻_𝐺/(1−𝜆). Saturation occurs if |𝑢_𝑛|>1, and the control is given by (Sign𝐻_𝐺)/(1−𝜆).

Another regularization process is to use the refined cost min_𝑢(⋅)(𝐶(𝑥(𝑡_𝑓))+(1−𝜆)∫^𝑡𝑓₀|𝑢|^2−𝜆𝑑𝑡). 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 𝑞₁(𝑡_𝑓)=0 and the cost −|𝑞²₂(𝑡_𝑓)| by a cost function of the form𝐶_𝛼=𝛼𝑞²₁𝑡_𝑓−𝑞²₂𝑡_𝑓, (6.2) 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 𝐵𝑆𝐵𝑆… 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) where the coordinates of the magnetization vector are (𝑀_𝑥,𝑀_𝑦,𝑀_𝑧), 𝑇₁ and 𝑇₂ are the relaxation times, and (𝜔_𝑥,𝜔_𝑦) are the two components of the magnetic field 𝜔 which is bounded by |𝜔|≤𝜔_max. We use the normalized coordinates (𝑥,𝑦,𝑧)=(𝑀_𝑥,𝑀_𝑦,𝑀_𝑧)/𝑀₀, 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 𝑢=(𝑢_𝑥,𝑢_𝑦)=(2𝜋/𝜔_max)(𝜔_𝑥,𝜔_𝑦), leading to |𝑢|≤2𝜋, while the normalized time is given by 𝜏=𝜔_max𝑡/(2𝜋). In these new coordinates, the system takes the form𝑥̇=−Γ𝑥+𝑢_𝑦𝑧,𝑦̇=−Γ𝑦−𝑢_𝑥̇()𝑧,𝑧=𝛾1−𝑧+𝑢_𝑥𝑦−𝑢_𝑦𝑥, (6.4) with Γ=2𝜋/(𝜔_max𝑇₂) and 𝛾=2𝜋/(𝜔_max𝑇₁). 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 𝑇₁=2000 ms and 𝑇₂=200 ms, while for the second spin 𝑇₁=𝑇₂=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 (𝑦₁,𝑧₁) and (𝑦₂,𝑧₂). 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−𝑧²₀,𝑧₀),(−1−𝑧²₀,𝑧₀)) where 𝑧=𝑧₀ 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.

(a)

(b)

Figure 3 

Projection of the singular flow onto the planes (𝑦1,𝑧1) and (𝑦2,𝑧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 105). The horizontal and vertical dashed lines represent the 𝑦- and 𝑧-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)

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.

(a)

(b)

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.

(a)

(b)

Figure 6 

Zoom of the results of Figure 5.

In order to improve the contrast, a locally optimal BS-sequence solution of the maximum principle is computed as follows. The transfer time 𝑡_𝑓 is fixed, and the Mayer problem is regularized by using a cost of the form𝐶𝑥𝑡_𝑓+(1−𝜆)^𝑡𝑓₀|𝑢|^2−𝜆𝑑𝑡. (6.5) 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 2𝑇_min where 𝑇_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 𝑡_𝑓=1.1𝑇_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 2𝑇_min.

(a)

(b)

(c)

(d)

Figure 7 

Evolution of the control field 𝑢 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 2𝑇min.

(a)

(b)

(c)

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

(a)

(b)

(c)

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.

(a)

(b)

(c)

Figure 10 

It is the same as Figure 9 but for a control field of the form BSBSBS with 𝑡𝑓=1.1𝑇min.

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

(a)

(b)

(c)

Figure 12 

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 |𝑀𝑦|, while the panel (c) to |𝑀𝑧|, that is, the remaining 𝑦 and 𝑧 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 𝑇1=1.3 ms, 𝑇2=50 ms for the inner sample and 𝑇1=1.3 ms, 𝑇2=200 ms for the outside sample.

Acknowledgments

The authors acknowledge O. Cots and Y. Zhang for useful discussions and for providing them with some numerical and experimental results of this paper. B. Bonnard and D. Sugny acknowledge support from the PEPS INSIS optimal control of spin dynamics in nuclear magnetic resonance imaging.