#### 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 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] 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 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 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] 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 ,(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 . A regular extremal is formed by bang arcs defined by , while the singular extremals satisfy . The optimal control problem reduces to find the sequence , where is a bang arc with and a singular arc. For linear systems , under the mild Kalman controllability condition rank , 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 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 where the state variable belongs to the Bloch ball which is invariant for the dynamics since the parameters belong to the set . 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: 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 , the goal is to reach in a given transfer time the final state (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 where is a mapping. The class of admissible control is the set 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 . We fix a terminal manifold defined by 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 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 where is the usual scalar product in , and is the adjoint vector. is called the pseudo-Hamiltonian, and we introduce the maximized Hamiltonian .

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 such that the following conditions are satisfied:
with the maximality condition
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 ,
(ii)*Mayer problem*: at the final time ,

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 . 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 . The following results are well known.

*Definition 2.1. *Relaxing , a singular trajectory of the system on is a control trajectory pair such that the input-output mapping
is singular.

Proposition 2.2. *One then deduces the following results.*(1)*If is not singular on , then the input-output mapping is open at for the -norm topology.*(2)*If is singular on , then there exists such that satisfies a.e. on the weak maximum principle*

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 . Let us denote by the Hamiltonian lifts, then the pseudo-Hamiltonian takes the form . The maximality condition (2.4) leads to the following parameterization of the extremal controls: outside the switching surface . Plugging such into the Hamiltonian gives the true Hamiltonian . The smooth solutions of the corresponding Hamiltonian field are called extremals of order zero.

The relation leads to the following interpretation.

Proposition 2.3. *Extremals of order zero correspond to the singularity of the input-output mapping
*

##### 2.4. Single-Input Affine Case

The system is of the form and . Applying the maximality condition, one gets two types of extremals.(i)Regular extremals: the control is given by , , and . If , 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 , identically. Differentiating twice with respect to time, we get the relations and from this second condition, one derives when the denominator is not vanishing the corresponding singular control

Plugging such into the pseudo-Hamiltonian defines the singular flow solution of the Hamiltonian vector field denoted as and starting at from the two constraints and . They have the following interpretation.

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

##### 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 and if it is strict, it is called the strict Legendre-Clebsch condition. In the singular case, one has . 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.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 , and one can consider the extension of the control system by setting , , 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 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 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 of

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 . The variational equation
is called the Jacobi equation. A Jacobi field is a nontrivial solution , and it is said to be vertical at time if where is the canonical projection.

The following geometric result is crucial [13].

Proposition 3.2. *Let be the fiber , and let 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 . Moreover, the rank of the restriction to of is at most .*

This leads to the following definition.

*Definition 3.3. *We define the exponential mapping
as the projection on the state space of the integral curve of with initial condition , where can be restricted by linearity to the sphere . If is the reference extremal, a time is said to be geometrically conjugate to 0 if the mapping is not of rank at , and the associated point is said to be geometrically conjugate to . We denote by 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 generated by the Jacobi fields , vertical at and such that is orthogonal to . At a conjugate time , one has

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 , associated to is by definition a singularity of the input-output mapping . One assumes that on each subinterval , , the singularity of on is of codimension one.(A3) We are in the nonexceptional case where along the reference extremal.

Theorem 3.4. *Under our assumptions, the first geometric conjugate time is the first time such that if , the input-output mapping is not open at , and if , it becomes open at for the -norm topology.*

Corollary 3.5. *Consider the control system , , and let be a reference extremal of order zero solution of such that assumptions (A2), (A3) are satisfied, then the first time 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 , . 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 , 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 , 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 , then we can embed locally the reference extremal into a central field formed by all extremal trajectories starting from . One restricts the adjoint vector with the normalization and the assumption (A3). One introduces where is the train of Lagrangian manifolds generated by the Hamiltonian flow , with 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
**
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 to the center of the Bloch ball. Due to the symmetry of revolution, one can restrict the problem to the 2D-single input system , 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 to any point of the Bloch ball. By symmetry, one can only consider the domain . 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 and the Lie brackets up to order three are where , The singular trajectories are contained in the set defined by which leads to . Hence, the singular locus is the union of the two singular lines: the line corresponding to the axis of revolution and the horizontal axis . Eliminating , the singular control is given by where and . Computing, one gets(i)for , this simplifies into and . Hence, the singular control is zero, and the dynamics is governed by the equation ,(ii)for , we get and , which gives where . The flow on the horizontal direction is given by

##### 4.2. Parameter Conditions

The interesting case is when and . This leads to In this case, when . Moreover, we have the following estimate. Setting , we get near 0 that , and the trajectory reaching zero in is of order , the singular control being of order . In particular, one deduces the following.

Lemma 4.1. *The singular control in the neighborhood of the point , along the horizontal direction is of order 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 from the maximum principle, and from the Legendre-Clebsch condition (2.15), we deduce that . More precisely, if we introduce , one arrives at(i)hyperbolic case , the singular direction is small time minimizing;(ii)elliptic case , the singular trajectory is small time maximizing.

Applying this test when , one has(i)the horizontal singular line is a fast direction;(ii)the vertical singular line is fast provided .

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 with a bang pulse, a zero control is then applied to reach the target along the horizontal axis, using in particular the 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 . From a straightforward computation, one deduces that it will form an oval curve defined by , which shrinks into a point if . Outside this set, is defined by and , the sign of being given by . This form allows to compare the time to travel different trajectories not crossing by using the Stokes theorem. We have also 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 , 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 , and we denote by and bang arcs for the system , . 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 and where is the switching function evaluated along an extremal arc . According to the maximum principle, near , the extremal is of the form (resp., ) if (resp., ).

*Fold case:* denoting the second-order derivative along a bang arc , the fold case is the situation where and . We have three cases.(i)Hyperbolic Case. at the contact point with , and . The connection with a singular extremal with 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 , and . 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 or if .

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 The singular control is , and it is not admissible for . 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 and the control bound is large enough such that the bang arc starting from the north pole intersects the singular arc , 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 . We must have We denote by the solution of the variational equation such that , where the equation is integrated backwards from time to 0. By construction (since and the derivative is zero) and hence at time , is orthogonal to and to . Therefore, and are collinear. We denote by the angle between and measured counterclockwise. One deduces that switching occurs when . Moreover, on the singular set .

We have by definition and in the analytic case, the ad-formula gives, for small , 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 identified with the set of matrices of : acting on the subspace of vector in . Lie brackets of two affine vector fields , being defined by the computation of the exponential is a standard exercise in linear algebra which amounts to compute a Jordan normal form of , 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 and . Using the symmetry of revolution with respect to the -axis, one can restrict the study to the domain . 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.

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.

#### 5. The Contrast Problem

For simplicity, we will limit our analysis to the single-input case. The system can be written shortly as where , representing the state of each spin 1/2 particle. From Bloch equation, we get 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 while the singular control is defined by the Lie brackets being computed in Section 2.

In order to be optimal, the singular extremals have to satisfy the generalized Legendre-Clebsch condition and the transversality condition at the final time. Decomposing the adjoint vector in the form where is associated to the spin , one has and if , it can be normalized by homogeneity to . The case 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 , which gives in the singular case. With this restriction, the adjoint vector can be eliminated, and the singular control in this exceptional case is given by where with the corresponding vector fields defined by Observe that in the general case, a similar computation shows that the singular extremals are solutions of an equation of the form where is a one-dimensional time-dependent parameter whose dynamics is deduced from the adjoint equation.

##### 5.3. Desingularization

Meromorphic differential equations of the form (5.8) can be desingularized as using the time reparameterization . According to the formula (5.4), the singular control can explode near , 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 , provided .

##### 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
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 , and the equation is where . We define the action of on by the action of the symplectic change of coordinates
a feedback acting trivially.

We have the following theorem [23].

Theorem 5.4. *The following diagram is commutative:
**
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 .

##### 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 restricted to the surface 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: and the differential equation (5.14) can be desingularized using the time reparameterization . One gets the equations which restricted to are a semicovariant.

Geometric invariants are related to the surface 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 -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 , , and let be a smooth singular control such that with corresponding trajectory . We have the following.

Theorem 6.1. *Under generic assumptions, the first conjugate time is the time along the singular arc such that a policy is time extremizing in a -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 or (the maximum bound being normalized to 1), and by symmetry of revolution, one can assume that . At the final time , one must have , and moreover . 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 , and the singular case called cheap control is the case where , and some control components are not penalized. A regularization process consists of adding quadratic penalties in the cost 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 , 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: , . For the regularized problem, one considers only normal extremals associated to the pseudo-Hamiltonian The control is computed in the normal case by solving , . Saturation occurs if , and the control is given by .

Another regularization process is to use the refined cost . 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 and the cost by a cost function of the form where is a weight parameter. Another advantage of the homotopy method is, once the structure of the extremal trajectory is determined by the continuation method, the structure being a sequence of saturating controls and normal extremals, the true sequence 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 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 . We use the normalized coordinates , which entails that the initial state is by definition the north pole of the Bloch sphere. The normalized control is defined as , leading to , while the normalized time is given by . In these new coordinates, the system takes the form with and . From an experimental point of view, 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 ms and ms, while for the second spin 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 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)**

**(a)**

**(b)**

**(a)**

**(b)**

**(a)**

**(b)**

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 The objective of this continuation is to determine the structure of the extremal control at the limit . 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 where is the minimum time to transfer the first spin to 0. Figure 8 displays the control and the state variable when . 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 , 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 .

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

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

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