# A combination of algebraic, geometric and numerical methods in the contrast problem by saturation in magnetic resonance imaging

Bernard Bonnard, Mathieu Claeys, Olivier Cots, Alain Jacquemard, Pierre Martinon

<figure><img src="image_1.png" /></figure>

To cite this version:

Bernard Bonnard, Mathieu Claeys, Olivier Cots, Alain Jacquemard, Pierre Martinon. A combination of algebraic, geometric and numerical methods in the contrast problem by saturation in magnetic resonance imaging. 2014. hal-01001975

HAL Id: hal-01001975

https://inria.hal.science/hal-01001975

Submitted on 5 Jun 2014

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A COMBINATION OF ALGEBRAIC, GEOMETRIC AND NUMERICAL METHODS IN THE CONTRAST PROBLEM BY SATURATION IN MAGNETIC RESONANCE IMAGING

BERNARD BONNARD*, MATHIEU CLAEYS†, OLIVIER COTS‡, ALAIN JACQUEMARD§, AND PIERRE MARTINON¶

**Abstract.** In this article, the contrast imaging problem by saturation in nuclear magnetic resonance is modeled as a Mayer problem in optimal control. The optimal solution can be found as an extremal solution of the Maximum Principle and analyzed with the recent advanced techniques of geometric optimal control. This leads to a numerical investigation based on shooting and continuation methods implemented in the HamPath software. The results are compared with a direct approach to the optimization problem and implemented within the Bocop toolbox. In complement lmi techniques are used to estimate a global optimum. It is completed with the analysis of the saturation problem of an ensemble of spin particles to deal with magnetic fields inhomogeneities.

**Key words.** Geometric optimal control, Contrast imaging in NMR, Direct method, Shooting and continuation techniques, Moment optimization.

AMS subject classifications.

1. Introduction. The Bloch equations are a set of macroscopic equations which accurately describe the experimental model in nuclear magnetic resonance (NMR) and magnetic resonance imaging (MRI) [20]. The spin-1/2 particle is represented by a magnetization vector $M = (M_x, M_y, M_z)$ in the laboratory reference frame which evolves according to:

$$ (1.1) \qquad \frac{dM}{d\tau} = \gamma M \wedge B + R(M) $$

where $\gamma$ is the gyromagnetic ratio, $B = (B_x, B_y, B_z)$ is the applied magnetic field which decomposes into a strong polarizing magnetic field $B_z = B_0$ in the z-direction, while $B_x, B_y$ are the components of a Rf-magnetic field in the transverse direction and $R(M)$ is the dissipation of the form:

$$ \left( -\frac{M_x}{T_2}, -\frac{M_y}{T_2}, -\frac{(M_z - M_0)}{T_1} \right) $$

where $T_1, T_2$ are the longitudinal and transverse relaxation parameters characteristic of the species. The point $(0, 0, M_0)$ represents the equilibrium of the free motion and $M_0$ can be normalized to 1 and $N = (0, 0, 1)$ denotes the North pole of the Bloch ball: $|M| \le 1$.

Denoting $\omega_0 = -\gamma B_0$ the resonant frequency, and let $u(\tau) = -\gamma B_y, v(\tau) = -\gamma B_x$, the Bloch equations in the stationary frame of reference can be written in the matrix

<footnote>*Institut de Mathématiques de Bourgogne, Université de Bourgogne, 9 avenue Alain Savary, 21078 Dijon, France, bernard.bonnard@u-bourgogne.fr, on leave at Inria Sophia Antipolis Méditerranée, 06902 Sophia Antipolis, France.</footnote>

<footnote>†Department of Engineering, University of Cambridge, United Kingdom and LAAS, Université de Toulouse, France, mathieu.claeys@eng.cam.ac.uk</footnote>

<footnote>‡Inria Sophia Antipolis Méditerranée, 06902 Sophia Antipolis, France, olivier.cots@inria.fr</footnote>

<footnote>§Institut de Mathématiques de Bourgogne, Université de Bourgogne, 9 avenue Alain Savary, 21078 Dijon, France, alain.jacquemard@u-bourgogne.fr</footnote>

<footnote>¶Inria and Ecole Polytechnique, 91128 Palaiseau France, martinon@cmap.polytechnique.fr</footnote>
form:

$$ (1.2) \qquad \frac{d}{d\tau} \begin{bmatrix} M_x \\ M_y \\ M_z \end{bmatrix} = \begin{bmatrix} 0 & -\omega_0 & u(\tau) \\ \omega_0 & 0 & -v(\tau) \\ -u(\tau) & v(\tau) & 0 \end{bmatrix} \begin{bmatrix} M_x \\ M_y \\ M_z \end{bmatrix} - \begin{bmatrix} \frac{M_x}{T_2} \\ \frac{M_y}{T_2} \\ \frac{M_z}{T_1} - 1 \end{bmatrix} $$

The Bloch equations can be represented in a rotating frame of reference: $S(\tau) = \exp(\tau\omega\Omega_z)$, $M = S(\tau)q$, $q = (x, y, z)$,

$$ \Omega_z = \begin{bmatrix} 0 & -1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} $$

and introducing the Rf-field representation:

$$ u_1 = u \cos \omega\tau - v \sin \omega\tau $$

$$ v_1 = u \sin \omega\tau + v \cos \omega\tau, $$

one gets the Bloch equations in the moving frame:

$$ (\text{Bloch3}) \qquad \frac{d}{d\tau} \begin{bmatrix} x \\ y \\ z \end{bmatrix} = \begin{bmatrix} 0 & -\Delta\omega & u_2 \\ \Delta\omega & 0 & -u_1 \\ -u_2 & u_1 & 0 \end{bmatrix} \begin{bmatrix} x \\ y \\ z \end{bmatrix} - \begin{bmatrix} \frac{x}{T_2} \\ \frac{y}{T_2} \\ \frac{z-1}{T_1} \end{bmatrix} $$

where $\Delta\omega = \omega_0 - \omega$ is the resonance offset.

The control is bounded by $m, m = 2\pi \times 32.3\text{Hz}$ being the experimental intensity of the experiments. Assuming $\Delta\omega = 0$, and normalizing the time $t = \tau m$, denoting $\Gamma = 1/mT_2$, $\gamma = 1/mT_1$, $2\Gamma \ge \gamma \ge 0$, the system is normalized to:

$$ \frac{dx}{dt} = -\Gamma x + u_2 z $$

$$ (1.3) \qquad \frac{dy}{dt} = -\Gamma y - u_1 z $$

$$ \frac{dz}{dt} = \gamma(1-z) + u_1 y - u_2 x, $$

where $|u| \le 1$. This equation describes the evolution of the magnetization vector in NMR but in MRI, they are some distortion due to spatial position of the spin in the image and $B_0$, $B_1$ inhomogeneities, variation of $B_0$ producing a resonance offset and $B_1$-inhomogeneity implies a variation of the amplitude of the Rf-field.

In this article we shall restrict to consider the $B_1$-inhomogeneity, the main objective being to test the numerical methods. This leads to introduce in (1.3) a scaling factor $a_i \ge 0$ depending on the spatial position of the spin in the image and the equation transforms into

$$ \frac{dx}{dt} = -\Gamma x + a_i u_2 z $$

$$ (1.4) \qquad \frac{dy}{dt} = -\Gamma y - a_i u_1 z $$

$$ \frac{dz}{dt} = \gamma(1-z) + a_i(u_1 y - u_2 x), $$
where $|u| \le 1$, $q = (x, y, z)$.

The saturation problem in NMR consists to bring the magnetization vector $q$ from the North pole (0, 0, 1) to the center of the Bloch ball. A very important success of geometric control problem was an explicit solution to the saturation in minimum time [17] which allows to compute the physical limit to switch off the signal in NMR. This leads to introduce in MRI the saturation problem for an ensemble of N spins represented by $q_i$, $i = 1, \dots, N$ and the associated optimal control is formulated as the Mayer problem $\min_{v(\cdot)} \sum_{i=1}^{N} |q_i(t_f)|^2$ [14, 21], where $t_f$ is the transfer time and corresponds to compute a so-called robust control in the saturation problem, absorbing the magnetic fields inhomogeneities.

This study is an important issue in the research project combining theoretical analysis, numerical studies and experiments in the contrast problem by saturation in MRI [6, 9, 18] and improving preliminar results based on the grape algorithm [15]. The problem being to distinguish two species, e.g. oxygenated-desoxygenated blood, with different relaxation parameters ($\gamma_1, \Gamma_1$), ($\gamma_2, \Gamma_2$) by steering the first species to saturation: $q_1(t_f) = 0$, while maximizing the contrast $|q_2(tf)|$ of the second species.

Using the symmetry of revolution along the z-axis, we restrict our study to an ensemble of 2D-spin systems, controlled by the single control $u_1$:

$$ (1.5) \qquad \begin{aligned} \frac{dy}{dt} &= -\Gamma y - a_i u_1 z \\ \frac{dz}{dt} &= \gamma(1-z) + a_i u_1 y, \end{aligned} $$

where $|u_1| \le 1$, and the control associated to the ith-spin is scaled by $a_i > 0$.

In section 2 we present the theoretic contribution for the contrast problem. The mimizers are selected using the Maximum Principle and the optimal solutions are formed concatenating Bang solutions, where $|u_1| = 1$ and Singular ones. The problem is to find the BS-sequence. The analysis is similar in the multisaturation problem for an ensemble of two spins.

In section 3 the contrast problem is analyzed combining an indirect method where the solutions of the Maximum Principle are analyzed using a numerical and continuation method implemented in the HamPath<fnref n="1" /> software [12] and a direct method based on the Bocop<fnref n="2" /> toolbox [1]. The numerical solutions are only analyzed while the methods are detailed in [7] and the global optimum is estimated using lmi<fnref n="3" /> techniques [13, 19].

In section 4 we present algebraic-geometric tools to analyse the saturation problem of an ensemble of two spins, and numerical simulations using indirect and direct methods in relation with the problem of $B_1$-inhomogeneity.

## 2. The Maximum Principle and the geometric analysis of the extremals.

### 2.1. Geometric setting and Maximum Principle in a Mayer problem.

Consider a smooth control system of the form:

$$ \frac{dq}{dt} = F(q, u), \quad q \in \mathbf{Q}, \quad u \in \mathbf{U}. $$

<footnote n="1">¹http://cots.perso.math.cnrs.fr/hampath</footnote>

<footnote n="2">²http://bocop.org</footnote>

<footnote n="3">³http://homepages.laas.fr/henrion/software/gloptipoly3</footnote>
Given an initial point $q_0 \in \mathbb{Q}$, a transfer time $t_f$ and a terminal manifold $\mathbb{Q}_f$ to reach: $\mathbb{Q}_f = \{f(q) = 0\}$, a *Mayer problem* is an optimal problem:

$$c(q(t_f)) \rightarrow \min_{u(\cdot)}$$

In the geometric framework, one introduces the accessibility set at time $t$: $A(q_0, t) = \cup_{u(\cdot)} q(t, q_0, u)$ and the cost extended manifold $\mathbb{Q}_f^m = \{f = 0; c = m\}$, where $m$ is the minimal cost. In this representation, the minimal cost amounts to find a control $u^*$ with corresponding trajectory such that $q^* = q(t_f, q_0, u^*)$ belongs to the boundary of the accessibility set (see Fig. 1) and the Maximum Principle asserts the existence of an adjoint vector $p^*(t_f)$ perpendicular to $\Pi$, $\Pi$ being a hyperplane separating $K(q_0, t_f)$ from $\mathbb{Q}_f^m$, where $K(q_0, t_f)$ is a convex approximating cone of $A(q_0, t_f)$ at $q^*$. This leads to the following necessary optimality conditions [22].

<figure><img src="image_3.png" /><figcaption>FIG. 1. Optimal solution in a Mayer problem: $q^*$ lies on the boundary of the accessibility set $A(q_0, t_f)$ and terminal manifold $\mathbb{Q}_f^m$ with a minimal cost.</figcaption></figure>

**PROPOSITION 2.1.** If $u^*$ with corresponding trajectory $q^*$ on $[0, t_f]$ is optimal then there exists an adjoint vector $p^*(\cdot)$ such that denoting $H(q, p, u) = \langle p, F(q, u) \rangle$ the pseudo-Hamiltonian, we have for almost every $t \in [0, t_f]$:

$$(i) \dot{q}^* = \frac{\partial H}{\partial p}(q^*, p^*, u^*), \quad \dot{p}^* = -\frac{\partial H}{\partial q}(q^*, p^*, u^*)$$

$$(ii) H(q^*, p^*, u^*) = \max_{v \in U} H(q^*, p^*, v) \quad (\text{maximization condition})$$

and the following boundary conditions

$$(iii) f(q^*(t_f)) = 0$$

$$(iv) p^*(t_f) = p^0 \frac{\partial c}{\partial q}(q^*(t_f)) + \sum_{i=1}^{k} \sigma_i \frac{\partial f_i}{\partial q}(q^*(t_f)) \quad (\text{transversality condition}), \\ \sigma = (\sigma_1, \dots, \sigma_k) \in \mathbb{R}^k, p^0 \le 0$$

**DEFINITION 2.2.** The solutions of conditions (i) and (ii) of Prop. 2.1 are called extremals and BC-extremals if they satisfy the boundary conditions. Note that

$$M(q^*, p^*) = \max_{u \in U} H(q^*, p^*, u)$$

is constant and if $M(q^*, p^*) = 0$, an extremal is called exceptional.

*Application to MRI.* In the ideal contrast problem by saturation, since the Bloch ball is invariant for the dynamics, the state constraints are not active and the Maximum Principle can be applied. The state boundary conditions with $q = (q_1, q_2)$ and $f = 0$ is the set $q_1 = 0$. The cost function to minimize is $c(q(t_f)) = -|q_2(t_f)|^2$ and splitting the adjoint vector into $p = (p_1, p_2)$, the transversality condition is $p_2(t_f) = -2p^0 q_2(t_f)$, $p^0 \le 0$. If $p^0$ is nonzero, it can be normalized to $p^0 = -1/2$.
**2.2. Parameterization of the extremal curves.** The system is written as $\dot{q} = F_0(q) + u_1 F_1(q) + u_2 F_2(q)$, $|u| \le 1$. We denote $z = (q, p)$, $H_i(z) = \langle p, F_i(q) \rangle$, $i = 0, 1, 2$, the Hamiltonian lift of the system $\dot{z} = \vec{H}_0 + \sum_{i=1}^2 u_i \vec{H}_i$. If $(H_1, H_2) \ne 0$, the maximization condition in Prop. 2.1 leads to the following parameterization of the controls:

$$ (2.1) \qquad u_1 = \frac{H_1}{\sqrt{H_1^2 + H_2^2}}, \quad u_2 = \frac{H_2}{\sqrt{H_1^2 + H_2^2}}. $$

Define the switching surface:

$$ \Sigma: H_1 = H_2 = 0. $$

Plugging such a $u$ into the pseudo-Hamiltonian gives the true Hamiltonian $H_n = H_0 + \sqrt{H_1^2 + H_2^2}$. The corresponding extremals are called of order zero.

Besides those generic extremals, additional extremals are related to Lie algebraic properties of the system.

**2.3. Lie bracket computations.** The Lie bracket of two vector fields $F, G$, is computed with the convention:

$$ [F, G](q) = \frac{\partial F}{\partial q}(q)G(q) - \frac{\partial G}{\partial q}(q)F(q), $$

and if $H_F, H_G$, are the Hamiltonian lifts, the Poisson bracket is:

$$ \{H_F, H_G\}(z) = dH_F(\vec{H}_G)(z) = H_{[F,G]}(z). $$

To simplify the computations, each system is lifted on the semi-direct Lie product $GL(3, \mathbb{R}) \times_s \mathbb{R}^3$ acting on the $q$-space using the action $(A, a) = Aq + a$. The Lie bracket computation rule is $((A, a), (B, b)) = ([A, B], Ab - Ba)$ where $[A, B] = AB - BA$. We denote $F_0 = (A_0, a_0)$ with $A_0 = \text{diag}(-\Gamma_1, -\Gamma_1, -\gamma_1, -\Gamma_2, -\Gamma_2, -\gamma_2)$ and $a_0 = (0, 0, \gamma_1, 0, 0, \gamma_2)$, whereas the control fields $F_1, F_2$ are identified to $B_1 = \text{diag}(C_1, C_1)$ and $B_2 = \text{diag}(C_2, C_2)$, where $C_1$ and $C_2$ are the antisymmetric matrices $C_1 = E_{32} - E_{23}$, $C_2 = E_{13} - E_{31}$ with $E_{ij} = (\delta_{ij})$ (Kronecker symbol). See [6] for more details.

Next, we present in details the Lie brackets needed in our computations, each entry is formed by a couple $(v_1, v_2)$ and we use the notation omiting the indices. We set $\delta = \gamma - \Gamma$.

* Length 1:

$$ F_0 = (-\Gamma x, -\Gamma y, \gamma(1-z)) $$

$$ F_1 = (0, -z, y) $$

$$ F_2 = (z, 0, -x) $$

* Length 2:

$$ [F_0, F_1] = (0, \gamma - \delta z, -\delta y) $$

$$ [F_0, F_2] = (-\gamma + \delta z, 0, \delta x) $$

$$ [F_1, F_2] = (-y, x, 0) $$
• Length 3:

$$
\begin{align*}
& [[F_1, F_2], F_0] = 0 \\
& [[F_1, F_2], F_1] = F_2 \\
& [[F_1, F_2], F_2] = -F_1 \\
& [[F_0, F_1], F_1] = (0, -2\delta y, -\gamma + 2\delta z) \\
& [[F_0, F_1], F_2] = (\delta y, \delta x, 0) = [[F_0, F_2], F_1] \\
& [[F_0, F_2], F_2] = (-2\delta x, 0, 2\delta z - \gamma) \\
& [[F_0, F_1], F_0] = (0, -\gamma(\gamma - 2\Gamma) + \delta^2 z, -\delta^2 y) \\
& [[F_0, F_2], F_0] = (\gamma(\gamma - 2\Gamma) - \delta^2 z, 0, \delta^2 x)
\end{align*}
$$

**2.4. Stratification of the surface Σ : H₁ = H₂ = 0 and partial classification of the extremal flow near Σ.** Let z = (q, p) be a curve solution of
**H̃0 + u₁H̃1 + u₂H̃2.** Differentiating H₁ and H₂ along such a solution, one gets:

$$
\begin{equation} \tag{2.2}
\begin{aligned}
\dot{H}_1 &= \{H_0, H_1\} - u_2\{H_1, H_2\} \\
\dot{H}_2 &= \{H_0, H_2\} + u_1\{H_1, H_2\}
\end{aligned}
\end{equation}
$$

Hence we have:

PROPOSITION 2.3. Let $z_0 \in \Sigma_1 = \Sigma \setminus \{H_1, H_2\} = 0$ and define the control $u_s$ by:

$$
(2.3) \qquad u_s(z) = \frac{(-\{H_0, H_2\}(z), \{H_0, H_1\}(z))}{\{H_1, H_2\}(z)},
$$

and plugging such $u_s$ into $H$ defines the true Hamiltonian

$$
H_s(z) = H_0(z) + u_{s,1}(z)H_1(z) + u_{s,2}(z)H_2(z)
$$

which parameterize the singular solutions of the bi-input system contained in Σ₁. This gives the first stratum of the surface Σ. Moreover, the behaviours of the extremals of order zero near a point z₀ of Σ₁ can be easily analyzed using (2.2) and a nilpotent model where all Lie brackets at z₀ ∈ Σ₁ of length ≥ 3 are zero. See [3, 10] for similar computations. Denoting:

$$
\{H_0, H_1\}(z_0) = a_1, \quad \{H_0, H_2\}(z_0) = a_2, \quad \{H_1, H_2\}(z_0) = b
$$

and using polar coordinates $H_1 = r \cos \theta$, $H_2 = r \sin \theta$, then (2.2) becomes:

$$
\begin{equation} \tag{2.4}
\begin{aligned}
\dot{r} &= a_1 \cos \theta + a_2 \sin \theta \\
\dot{\theta} &= \frac{1}{r} (b - a_1 \sin \theta + a_2 \cos \theta).
\end{aligned}
\end{equation}
$$

To analyze this equation, we write:

$$
a_1 \sin \theta - a_2 \cos \theta = A \sin(\theta + \phi)
$$

with $A \tan \phi = -a_2/a_1$, $A = \sqrt{a_1^2 + a_2^2}$. Hence the equation $\dot{\theta} = 0$ leads to the relation

$$
A \sin(\theta + \phi) = b,
$$
which has two distinct solutions on $[0, 2\pi]$ denoted $\theta_0, \theta_1$ if and only if $A > |b|$, one solution if $A = |b|$ and zero solution if $A < |b|$. Moreover $\theta_1 - \theta_0 = \pi$ if and only if $b = 0$. Plugging $\theta_0, \theta_1$ in (2.4), one gets two solutions of (2.4). Hence we deduce:

LEMMA 2.4. If $\sqrt{a_1^2 + a_2^2} > |b|$ and $b \neq 0$, we have a broken extremal formed by concatenating two extremals of order zero at each point $z_0$ of $\Sigma_1$.

At such a point $z_0$ of $\Sigma_1$, the singular control (2.3) is such that

$$u_{s,1}^2 + u_{s,2}^2 = \frac{a_1^2 + a_2^2}{b^2} > 1$$

and hence is not admissible. Next we analyze more degenerated situations and one needs the following concept.

Goh condition. Higher order necessary optimality conditions along singular extremals in the bi-input case are related to finiteness of the index of the quadratic form associated to the intrinsic second order derivative [4]:

$$(2.5) \qquad \{H_1, H_2\} = 0 \text{ (Goh condition).}$$

Using $H_1 = H_2 = \{H_1, H_2\} = 0$ and (2.2), one gets the additional relations:

$$(2.6) \qquad \{H_1, H_2\} = \{H_0, H_1\} = \{H_0, H_2\} = 0.$$

Then differentiating again along a solution leads to the relations:

$$(2.7) \qquad \{\{H_1, H_2\}, H_0\} + u_1\{\{H_1, H_2\}, H_1\} + u_2\{\{H_1, H_2\}, H_2\} = 0$$

$$(2.8) \qquad \left\{ \begin{array}{l} \{\{H_0, H_1\}, H_0\} + u_1\{\{H_0, H_1\}, H_1\} + u_2\{\{H_0, H_1\}, H_2\} = 0 \\ \{\{H_0, H_2\}, H_0\} + u_1\{\{H_0, H_2\}, H_1\} + u_2\{\{H_0, H_2\}, H_2\} = 0 \end{array} \right.$$

This leads in general to three relations to compute two control components, and for a generic system such conditions are not satisfied [11], but in our case, according to Lie brackets computations, we have:

LEMMA 2.5. If $H_1 = H_2 = 0$, one has

$$\{\{H_1, H_2\}, H_0\} = \{\{H_1, H_2\}, H_1\} = \{\{H_1, H_2\}, H_2\} = 0$$

and (2.7) is satisfied.

To analyze (2.8), which can be written $\tilde{A} + \tilde{B}u = 0$ and if $\det \tilde{B} \neq 0$, the corresponding singular control is given by:

$$(2.9) \qquad u'_s(z) = -\tilde{B}^{-1}(z)\tilde{A}(z)$$

Using the relations:

$$H_1 = H_2 = \{H_1, H_2\} = \{H_0, H_1\} = \{H_0, H_2\} = 0,$$

the vector $p$ is orthogonal to $F_1, F_2, [F_1, F_2], [F_0, F_1], [F_0, F_2]$. Introducing:

$$A = \begin{pmatrix} A_1 \\ A_2 \end{pmatrix}, \quad B = \begin{pmatrix} B_1 & B_3 \\ B_2 & B_4 \end{pmatrix}, \quad C = (F_1, F_2, [F_1, F_2], [F_0, F_1], [F_0, F_2]),$$

with

$$A_1 = \det(C, [[F_0, F_1], F_0]), \quad A_2 = \det(C, [[F_0, F_2], F_0]),$$
and

$$B_1 = \det(C, [[F_0, F_1], F_1]), \quad B_2 = \det(C, [[F_0, F_2], F_1]),$$

$$B_3 = \det(C, [[F_0, F_1], F_2]), \quad B_4 = \det(C, [[F_0, F_2], F_2]),$$

the relation (2.8) leads to:

$$A + Bu = 0,$$

and if $\det B \neq 0$, one gets the singular control given by the feedback:

$$(2.10) \qquad u'_s(q) = -B^{-1}(q)A(q)$$

and the associated vector field:

$$Q'_s = F_0 + u'_{s,1}F_1 + u'_{s,2}F_2.$$

Moreover, the singular control has to be admissible: $|u'_s| \le 1$. We introduce the stratum:

$$\Sigma_2: H_1 = H_2 = \{H_1, H_2\} = \{H_0, H_1\} = \{H_0, H_2\} = 0 \setminus \det \tilde{B} = 0.$$

Hence:

LEMMA 2.6.

1. On the stratum $\Sigma_2$, there exists singular extremals satisfying Goh condition where the singular control is given by the feedback (2.9).

2. For the contrast problem:

$$(2.11) \qquad \det B = (x_1 y_2 - x_2 y_1)^4 (\delta_1 - \delta_2)(2\delta_1 z_1 - \gamma_1)(2\delta_2 z_2 - \gamma_2) \\ \qquad \qquad \qquad \qquad (2(\delta_1^2 \gamma_2 z_1 - \delta_2^2 \gamma_1 z_2) - \gamma_1 \gamma_2 (\delta_1 - \delta_2) - 2\delta_1 \delta_2 (\gamma_1 z_2 - \gamma_2 z_1)),$$

The behaviours of the extremals of order zero near a point $z_0 \in \Sigma_2$ is a complicated problem. Additional singular extremals can be contained in the surface:

$$\Sigma_3: H_1 = H_2 = \{H_1, H_2\} = \{H_0, H_1\} = \{H_0, H_2\} = \det \tilde{B} = 0,$$

and they can be computed easily since the corresponding controls has to force the surface $\det B = 0$ to be invariant. Some have an important physical meaning, due to the symmetry of revolution of the Bloch equations. They correspond to control the system, imposing $u_2 = 0$. In this case, one can restrict the system to

$$\mathbf{Q} = \{q = (q_1, q_2) \in \mathbb{R}^n : |q_1| \le 1, |q_2| \le 1, x_1 = x_2 = 0\}.$$

The computations of the corresponding extremals amount to replace in the relations; $H_2$ by $\varepsilon H_2$ and to impose $\varepsilon = 0$. The remaining relations are then:

$$H_1 = \{H_0, H_1\} = 0$$

and from (2.8) one gets the relations:

$$(2.12) \qquad \{\{H_0, H_1\}, H_0\} + u_{1,s}\{\{H_0, H_1\}, H_1\} = 0,$$

and thus, this defines the singular control:

$$(2.13) \qquad u_{1,s} = - \frac{\{\{H_0, H_1\}, H_0\}}{\{\{H_0, H_1\}, H_1\}}$$
and the associated Hamiltonian $H_{1,s} = H_0 + u_{1,s}H_1$. We have the following result:

**PROPOSITION 2.7.** *The extremals of the single-input case are extremals of the bi-input case with the additional condition: $x_1 = p_{x_1} = x_2 = p_{x_2} = 0$.*

In order to be optimal, they have to satisfy the generalized Lengendre-Clebsch condition:

$$ (2.14) \qquad \frac{\partial}{\partial u_1} \frac{d^2}{dt^2} \frac{\partial H}{\partial u_1} = \{H_1, \{H_1, H_0\}\}(z) \le 0 $$

Observe that if we impose $u_2 = 0$, the classification of the extremals near the switching surface, which reduces to $H_1 = 0$, is a standard problem [16].

Finally, another important property of the extremal flow, again a consequence of the symmetry of revolution is given next, in relation with Goh condition. It is a consequence of Noether integrability theorem.

**PROPOSITION 2.8.** *In the contrast problem, for the Hamiltonian vector field $\vec{H}_n$ which solutions are extremals of order zero, the Hamiltonian lift $H(z) = \{H_1, H_2\}(z) = (p_{x_1}y_1 - p_{y_1}x_1) + (p_{x_2}y_2 - p_{y_2}x_2)$ is a first integral.*

**2.5. The saturation problem for a single spin in minimum time.** The time minimal saturation problem for a single spin is presented in details in [9] and we recall briefly the results. Denoting $q_1 = (y_1, z_1)$ the coordinates of the first spin. The singular trajectories are the two lines: the vertical axis of revolution $y_1 = 0$, where the singular control is zero and the horizontal line $z_1^* = \gamma/2(\gamma - \Gamma)$ where the control is $u_1 = \gamma(2\Gamma - \gamma)/2(\gamma - \Gamma)y_1$ and the dynamics is given by:

$$ \dot{y}_1 = -\Gamma y_1 - \frac{\gamma^2(2\Gamma - \gamma)}{4(\gamma - \Gamma)^2 y_1}. $$

Consider the physical situation: $2\Gamma \ge 3\gamma$ and the horizontal line is such that $0 > z_1 > -1$. Observe that along this line the singular control $\to \infty$ when $y_1 \to 0$. The optimal policy is given by the following theorem and represented on Fig. 2, compared with the standard inversion sequence $\delta_+^h \delta_s^v$ used in practice in NMR.

**THEOREM 2.9.** *In the time minimal saturation problem for a single spin, if $2\Gamma \ge 3\gamma$, the optimal policy is of the form $\delta_+^h \delta_s^v$ ($\delta_+^h \delta_s^v$ (or the symmetric policy with respect to the z-axis of revolution) where $\delta_s^h, \delta_s^v$ are respectively the horizontal and vertical singular line.*

<figure><img src="image_4.png" /><figcaption>FIG. 2. Time minimal solution (left) compared with inversion sequence (right).</figcaption></figure>
The corresponding time minimal solution is denoted $T_{\min}$. An estimate of $T_{\min}$ is: $T_{\min} \approx t_1 + t_2$ where $t_1$ is the limit time to steer $y_1 = -1$ to $y_1 = 0$ along the singular horizontal line:

$$t_1 = \int_0^1 \frac{dy_1}{\Gamma y_1 + \frac{a}{y_1}} = \frac{1}{2\Gamma} \left( \ln \left( 1 + \frac{a}{\Gamma} \right) - \ln \left( \frac{a}{\Gamma} \right) \right), \quad a = \frac{\gamma^2(2\Gamma - \gamma)}{4(\gamma - \Gamma)^2},$$

and $t_2$ the time to steer $(0, z_1^*)$ to $(0, 0)$ along the vertical singular line is $t_2 = \ln(1 - z_1)/\gamma$. In this policy, an important role is played by the intermediate Bang arcs $\delta_+$ connecting the two singular arcs and this leads to the following definition.

**DEFINITION 2.10.** A bridge is a Bang arc $\delta_+$ or $\delta_-$ corresponding to an extremal $z_\pm$ such that both extremities belong to $\Sigma'_1: H_1 = \{H_1, H_0\} = 0$.

**3. Numerical simulations, the ideal contrast problem.** The ideal contrast problem by saturation in the single-input case, can be summarized this way:

$$(\text{ICPS}) \quad \left\{ \begin{array}{l} c(q(t_f)) = -|q_2(t_f)|^2 \rightarrow \min_{u(\cdot)}, \quad \text{fixed } t_f \\ \dot{q} = F_0(q) + u_1 F_1(q), \\ q(0) = q_0 \\ q_1(t_f) = 0 \end{array} \right.$$

where $q = (q_1, q_2)$, $q_i = (y_i, z_i) \in \mathbb{R}^2$, $|q_i| \le 1$, $i = 1, 2$. The initial condition for each spin is $q_i(0) = (0, 1)$. The vector fields $F_0$ and $F_1$ are given by:

$$F_0(q) = \sum_{i=1,2} (-\Gamma_i y_i) \frac{\partial}{\partial y_i} + (\gamma_i(1-z_i)) \frac{\partial}{\partial z_i},$$

$$F_1(q) = \sum_{i=1,2} -z_i \frac{\partial}{\partial y_i} + y_i \frac{\partial}{\partial z_i},$$

where $\Lambda_i = (\gamma_i, \Gamma_i)$ are the physical parameters representing each spin.

We present the simulations using the numerical methods (see [7] for a complete description of the algorithms).

The simulations correspond to the two following sets of experimental data, with the relaxation times in seconds and $T_{\min}$ the solution of the time minimal saturation problem for a single spin, from section 2.5.

**P$_1$: Fluid case.**

Spin 1: Cerebrospinal fluid: $T_1 = 2$, $T_2 = 0.3$;

Spin 2: Water: $T_1 = 2.5 = T_2$.

$T_{\min} = 26.17040$.

**P$_2$: Blood case.**

Spin 1: Deoxygenated blood: $T_1 = 1.35$, $T_2 = 0.05$;

Spin 2: Oxygenated blood: $T_1 = 1.35$, $T_2 = 0.2$.

$T_{\min} = 6.7981$.

Optimal solutions of the contrast problem are concatenations of bang and singular extremals. For the following sections, we introduce some notations. We note BS the sequence composed by one bang arc ($\delta_+$ or $\delta_-$) followed by one singular arc ($\delta_s$), and nBS, $n > 1$, the concatenation of n BS-sequences.

**3.1. First results with fixed final time.** The first difficulty comes from the discontinuities of the optimal control structure. We need to know the control structure (meaning the number of Bang-Singular sequences) before calling the multiple
shooting method. The indirect method also typically requires a reasonable estimate for the control switching times, as well as the states and costates values at the initial and switching times. We use the Bocop software based upon direct methods to obtain approximate optimal solutions in order to initialize the indirect shooting, within the HamPath code. We recall that the costate (or adjoint state) for Pontryagin's Principle corresponds to the Lagrange multipliers for the dynamics constraints in the discretized problem, and can therefore be extracted from the solution of the direct method.

The only a priori information is the value of the minimum time transfer $T_{min}$, used to set the final time $t_f$ in the $[T_{min}, 3T_{min}]$ range. We note $t_f = \lambda T_{min}$ with $\lambda$ in $[1, 3]$. The state variables are initialized as constant functions equal to the initial state, i.e. $y_1(\cdot) = y_2(\cdot) = 0$, $z_1(\cdot) = z_2(\cdot) = 1$. For the control variables we use the three constant initializations $u_1(\cdot) \in \{0.1, 0.25, 0.5\}$. The discretization method used is implicit midpoint (2nd order) with a number of time steps set to $\lambda \times 100$. In order to improve convergence, we add a small regularization term to the objective to be minimized, $\varepsilon_{reg} \int_0^{t_f} |u(t)|^2 dt$, with $\varepsilon_{reg} = 10^{-3}$.

We repeat the optimizations for $\lambda$ in $\{1.1, 1.5, 1.8, 2.0, 3.0\}$ with the three control initializations, see Table. 1. The solutions from Bocop are used to initialize the continuations in HamPath, and we discuss in the following sections the results obtained with the indirect method. Both methods confirm the existence of many local solutions, as illustrated on Fig. 3 for $\lambda = 1.5$, due in particular to symmetry reasons.

<table><caption>TABLE 1</caption><thead><tr><th>λ</th><th>1.1</th><th>1.5</th><th>1.8</th><th>2</th><th>3</th></tr></thead><tbody><tr><td>u<sub>init</sub> : 0.1</td><td>0.636 (++)</td><td>0.678 (+ - +)</td><td>0.688 (+ - +)</td><td>0.702 (- +)</td><td>0.683 (- + - +)</td></tr><tr><td>u<sub>init</sub> : 0.25</td><td>FAIL</td><td>0.661 (+ + - +)</td><td>0.673 (+ + - +)</td><td>0.691 (- + +)</td><td>0.694 (+ - +)</td></tr><tr><td>u<sub>init</sub> : 0.5</td><td>0.636 (++)</td><td>0.684 (++)</td><td>0.699 (- +)</td><td>0.697 (++)</td><td>0.698 (++)</td></tr></tbody></table>

**Fluid case: Batch optimizations (Direct method).** For each value of $\lambda$ we test the three initializations for the control $u$, and record the value of the objective (i.e. the contrast), as well as the control structure (i.e. the signs of bang arcs). CPU times for a single optimization are less than one minute on a Intel Xeon 3.2GHz.

<figure><img src="image_5.png" /><figcaption>FIG. 3. **Fluid case: Two local solutions for $\lambda = 2.0$.** Trajectories for spin 1 and 2 in the $(y,z)$-plane are portrayed in the first two subgraphs of each subplot. The corresponding control is drawn in the bottom subgraph. The two bang arcs have the same sign for the left solution, whereas for the right solution, the two bang arcs are of opposite sign.</figcaption></figure>
**3.2. Second order conditions.** According to proposition 3.2 from [8], the non-existence of conjugate points on each singular arc of a candidate solution is a necessary condition of local optimality. See [8] for details about conjugate points in the contrast problem. Here, we compute for each singular arc of all the solutions from subsection 3.1, the first conjugate point along the arc, applying the algorithm presented in subsection 4.3 from [8]. None of the solutions has a conjugate point on a singular arc. Hence all the solutions satisfy the second order necessary condition of local optimality. Fig. 4 represents the computations of the two conjugate points (since the structure is 2BS) of the best solution with $\lambda = 2.0$ from subsection 3.1.

<figure><img src="image_6.png" /><figcaption>FIG. 4. **Fluid case: second order conditions.** Second order necessary condition checked on the best solution with $\lambda = 2.0$ from subsection 3.1. The rank condition from the algorithm presented in subsection 4.3 from [8] is evaluated along the two singular arcs. See [2] for details on the concept of conjugate times. On the left subplot, for each singular arc, the curve is reparameterized so that the final time corresponds to the abscissa 1 (vertical blue dashed line); the determinant associated with the rank condition is plotted (top subgraph), so there is a conjugate time whenever it vanishes (vertical red dashed lines). One observes that conjugate times on each arc are located after the (normalized to 1) final time, satisfying necessary condition of local optimality of the trajectory. At the bottom, the smallest singular value of the matrix whose rank we test is plotted, extracting only the relevant information to detect the rank drops. On the right subplot is presented a zoom of top-left subgraph near the two conjugate times.</figcaption></figure>

**3.3. Influence of the final time.** Given that the initial point (the North pole) is a stationary point, the contrast is an increasing function of $t_f$ acting as a parameter. Indeed, applying a zero control at $t = 0$ leaves the system in its initial state so there is an inclusion of admissible controls between problems when the final time is increased (and the bigger the set of controls, the larger the maximum contrast). Having increasing bounded (by one, which is the maximum possible contrast given the final condition on spin no. 1) functions, it is natural to expect asymptotes on each branch.

In both cases $P_1$ and $P_2$, the contrast problem has many local solutions, possibly with different control structures. Besides, the structure of the best policy can change depending on the final time. The possible change of structure along a single path of zeros is emphasized in Fig. 5. In this figure, the branch made of 2BS solutions is represented in blue, whereas the 3BS branch is the dashed red line. We also show a crossing between two value functions of two different paths of zeros in Fig. 6.

Then for each solution of each branch the second order necessary condition is checked as in subsection 3.2: the first conjugate point of each singular extremal is computed. There is no failure in this test, hence all the solutions satisfy the necessary
second order condition of local optimality. Fig. 7 presents the second order conditions
along the extended path from Fig. 5.

<figure><img src="image_7.png" /><figcaption>FIG. 5. Fluid case: influence of the final time. On the left subgraph are shown the control laws of solutions at $\lambda = 2$ and $\lambda = 1.32$ from path from the right subplot. For $\lambda = 1.32$, we can see the saturating singular arc around the normalized time $\tau = 0.92$ (the time is normalized to be between 0 and 1 for each solution). The 2BS solution at $\lambda = 1.32$ is used to initialize a multiple shooting with a 3BS structure and then to perform a new homotopy from $\lambda = 1.32$ to $\lambda = 1$. On the right subgraph is portrayed the two homotopies: the first from $\lambda = 2$ to $\lambda = 1.32$ and the second to $\lambda = 1$, with one more BS sequence. The value function, the norm of the initial adjoint vector, the norm of the shooting function and the switching times along the path are given. The blue color represents 2BS solutions while the red color is for 3BS structures. The dashed red lines come from the extended path after the change of structure detected around $\lambda = 1.32$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIG. 6. Fluid case: influence of the final time. Crossing between two branches with 3BS solutions. The crossing is around $\lambda = 1.0484$, see top subgraph. Thus for $\lambda \le 1.0484$, the best solution, locally, has a 3BS structure of the form $\delta_+\delta_s\delta_+\delta_s\delta_-$ (bottom-left subgraph) while for $\lambda \in [1.0484, 1.351]$ the best solution is of the form $\delta_+\delta_s\delta_-\delta_s\delta_-$ (bottom-right subgraph). On the two bottom subgraphs, the trajectories for spin 1 and 2 in the $(y,z)$-plane are portrayed with the corresponding control, both for $\lambda = 1.0484$.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>FIG. 7. **Fluid case: influence of the final time.** Second order necessary condition checked along the extended path from Fig. 5. For all solutions from $\lambda = 1$ to $\lambda = 3$ are computed the first conjugate times along each singular arc. For $\lambda \in [1, 1.32]$, the structure is 3BS and there are 3 singular arcs. For $\lambda \in [1.32, 3]$, there are 2 singular arcs. Each singular interval is normalized in such a way the initial time is 0 and the final time is 1. The lower dashed horizontal line represents the final time 1. There is no conjugate time before the normalized final time 1 which means that all solutions satisfy the second order necessary condition of local optimality. Note that at a magenta cross, around (1.32, 1), the control of the first singular arc saturates the constraint $|u| = 1$, and so no conjugate time is computed after this time.</figcaption></figure>

**3.4. Sub-optimal syntheses in fluid and blood cases.** We give the syntheses of locally optimal solutions obtained in the blood and fluid cases. Note that in the special case $t_f = T_{min}$, for both cases the solution is 2BS and of the form $\delta_+\delta_s\delta_+\delta_s$.

For the fluid case, the left subplot of Fig. 8 represents all the different branches we obtained by homotopy on $\lambda$. The greatest two value functions intersect around $t_f = 1.048T_{min}$. The right subplot shows the sub-optimal synthesis. The best policy is:

$$ (3.1) \qquad \begin{aligned} & \delta_+\delta_s\delta_+\delta_s && \text{for } \lambda \in [1.000, 1.006], \\ & \delta_+\delta_s\delta_+\delta_s\delta_-\delta_s && \text{for } \lambda \in [1.006, 1.048], \\ & \delta_+\delta_s\delta_-\delta_s\delta_-\delta_s && \text{for } \lambda \in [1.048, 1.351], \\ & \delta_+\delta_s\delta_-\delta_s && \text{for } \lambda \in [1.351, 3.000]. \end{aligned} $$

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>FIG. 8. **Fluid case, sub-optimal synthesis.** Illustration on the left subplot, of local solutions (each branch corresponds to a control structure). The suboptimal synthesis is plotted on right sub-plot. The colors are blue for 2BS structure, red for 3BS and green for 4BS. The best policy is $\delta_+\delta_s\delta_+\delta_s\delta_-\delta_s$ for $\lambda \le 1.0484$, and $\delta_+\delta_s\delta_-\delta_s\delta_-\delta_s$ for $\lambda \in [1.0484, 1.351]$. Then, for $\lambda \in [1.351, 3]$, the best policy is 2BS and of the form $\delta_+\delta_s\delta_-\delta_s$.</figcaption></figure>
For the blood case, the results are excerpted from [12]. The left subplot of Fig. 9 shows the contrast for five different components of {$h = 0$}, for final times $t_f \in [1, 2]T_{\min}$. The three black branches are made only of BS solutions whereas the two others are made of 2BS and 3BS solutions. To maximise the contrast, the best policy, drawn as solid lines, is:

$$ (3.2) \qquad \begin{aligned} & \delta_+\delta_s\delta_+\delta_s && \text{for } \lambda \in [1.000, 1+\varepsilon], \varepsilon > 0 \text{ small} \\ & \delta_+\delta_s && \text{for } \lambda \in [1+\varepsilon, 1.294], \\ & \delta_+\delta_s\delta_-\delta_s\delta_-\delta_s && \text{for } \lambda \in [1.294, 2.000]. \end{aligned} $$

<figure><img src="image_12.png" /><figcaption>FIG. 9. Blood case, sub-optimal synthesis. Illustration on the left subplot, of local solutions (each branch corresponds to a control structure). Best policy as solid lines, local solutions as dashed lines. The suboptimal synthesis is plotted on right subplot. The colors are black for BS structure, blue for 2BS and red for 3BS. The best policy is BS for $t_f \in (1, 1.294]T_{\min}$ and 3BS of the form $\delta_+\delta_s\delta_-\delta_s\delta_-\delta_s$ for $t_f \in (1.294, 2]T_{\min}$. In the special case $t_f = T_{\min}$, the solution is 2BS and of the form $\delta_+\delta_s\delta_+\delta_s$</figcaption></figure>

**3.5. Sub-optimal syntheses compared to global results.** We now apply the lmi method to the contrast problem, described in [7], in order to obtain upper bounds on the true contrast. Comparing these bounds to the contrast of our solutions then gives an insight about their global optimality.

Table 2 shows the evolution of the upper bound on the contrast in function of lmi relaxation order, for the fluid case with $t_f = T_{\min}$. As expected, the method yields a monotonically non-increasing sequence of sharper bounds. Relaxations of orders 4 and 5 yield very similar bounds, but this should not be interpreted as a termination criterion for the lmi method.

<table><caption>TABLE 2</caption><thead><tr><th>r</th><th>√-J<sup>r</sup><sub>M</sub></th><th>N<sub>r</sub></th><th>t<sub>r</sub></th></tr></thead><tbody><tr><td>1</td><td>0.8474</td><td>63</td><td>0.7</td></tr><tr><td>2</td><td>0.7552</td><td>378</td><td>3</td></tr><tr><td>3</td><td>0.6226</td><td>1386</td><td>14</td></tr><tr><td>4</td><td>0.6069</td><td>3861</td><td>332</td></tr><tr><td>5</td><td>0.6040</td><td>9009</td><td>8400</td></tr></tbody></table>

Fluid case, $t_f = T_{\min}$: upper bounds on contrast $\sqrt{-J_M^r}$, numbers of moments $N_r$ and CPU times $t_r$ in function of relaxation order $r$.

Figs. 10 and 11 compare the tightest upper bounds found by the lmi method with the best candidate solutions found by Bocop and HamPath, in both the blood
and fluid cases. The figures also represent the relative gap between the methods defined as $(C_{LMI} - C_H)/C_H$, where $C_{LMI}$ is the lmi upper bound and $C_H$ is the contrast found with HamPath. As such, this measure characterizes the optimality gap between the methods. It does not, however, specify which of the method(s) could be further improved. At the fifth relaxation, the average gap is around 11% in the blood case, which, given the application, is satisfactory on the experimental level. For the fluid case, the average gap on the contrast is about 2% at the fifth relaxation, which strongly suggest that the solution is actually a global optimum. The gap is even below the 1% mark for $t_f \le 2T_{min}$.

<figure><img src="image_13.png" /><figcaption>FIG. 10. Fluid case. Best upper bounds (dashed line) by the lmi method compared with best solutions by HamPath (solid line), and relative gap between the two.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIG. 11. Blood case. Best upper bounds (dashed line) by the lmi method compared with best solutions by HamPath (solid line), and relative gap between the two.</figcaption></figure>

**4. The saturation of an ensemble of spin particles with $B_1$ inhomogeneity.** As for the saturation problem for a single spin, the main problem is to analyze the singular flow and in particular the behaviour of the extremals with respect to the surface $S: \{\{H_1, H_0\}, H_1\} = 0$ where the control saturates the constraints. The computation of this surface is a complicated task in the framework of algebraic geometry. Hence, we restrict our analysis for the case of two spins.
**4.1. Preliminaries.** The singular extremals are defined by:

$$H_1 = \{H_1, H_0\} = \{\{H_1, H_0\}, H_0\} + u_{1,s}\{\{H_1, H_0\}, H_1\} = 0.$$

Furthermore in the exceptional case, we have the additional constraint: $H_0 = 0$. Defining $H_s = H_0 + u_{1,s}H_1$, the singular extremals are solutions of the constrained Hamiltonian equation:

$$(4.1) \qquad \frac{dz}{dt} = \vec{H}_s(z), \quad \Sigma'_1: H_1 = \{H_1, H_0\} = 0,$$

those equations defining an Hamiltonian vector field on the surface $\Sigma'_1 \setminus S$, restricting the standard symplectic form $\omega = dp \wedge dq$.

Using the notation: $\mathcal{D} = \{\{H_1, H_0\}, H_1\}$ and $\mathcal{D}' = \{\{H_1, H_0\}, H_0\}$, the differential system (4.1) can be desingularized using the time reparameterization: $ds = dt/\mathcal{D}(z(t))$ and we get the equations:

$$\frac{dq}{ds} = \mathcal{D}F_0 - \mathcal{D}'F_1, \quad \frac{dp}{ds} = -p \left( \mathcal{D} \frac{\partial F_0}{\partial q} - \mathcal{D}' \frac{\partial F_1}{\partial q} \right),$$

restricted to the surface $\Sigma'_1 \setminus S$.

In the case of two spins, since the state-space is of dimension four, using the two constraints $H_1 = \{H_1, H_0\} = 0$, the equations (4.1) can be reduced to an equation of the form:

$$(4.2) \qquad \frac{dq}{dt} = F_0(q) - \frac{\mathcal{D}'(q, \lambda)}{\mathcal{D}(q, \lambda)} F_1(q)$$

where $\lambda$ is a one-dimensional parameter whose dynamics is deduced from the adjoint equation. Still, this equation is very intricate, and we shall focus on to the exceptional case with the additional constraint $H_0 = 0$. According to the Maximum Principle, this case corresponds to the situation where the transfer time $t_f$ is not fixed. With this additional constraint, the adjoint vector can be eliminated and we have:

**PROPOSITION 4.1.** In the exceptional case, the singular control is given by the feedback: $u_{1,s}^e = -D'(q)/D(q)$, where:

$$D = \det(F_0, F_1, [F_1, F_0], [[F_1, F_0], F_1]), \quad D' = \det(F_0, F_1, [F_1, F_0], [[F_1, F_0], F_0])$$

with the corresponding vector field on $\mathbf{Q}$: $X_e = F_0 - \frac{D'}{D}F_1$, which can be again desingularized using the reparameterization: $ds = dt/D(q(t))$ to produce the smooth vector field:

$$(4.3) \qquad X_e^r = DF_0 - D'F_1$$

**4.2. Analysis of the singular flow.**

**4.2.1. Algebraic computations.** The point $N = ((0,1), (0,1))$ is a singular point of $X_e^r$ and under a translation $N$ is taken as the origin of the coordinates. We have:

$$F_0 = (-\Gamma y_1, -\gamma z_1, -\Gamma y_2, -\gamma z_2),$$

$$F_1 = ((-z_1 + 1), y_1), (1 - \varepsilon)(-z_2 + 1), y_2)$$
where $(1 - \varepsilon)$ denotes the control rescaling of the second spin.

We have $D = (1 - \varepsilon)\tilde{D}$, where $\tilde{D}$ is a quadric which decomposes into $h_2 + h_3 + h_4$ where $h_i$ are the homogeneous part of degree $i$:

$$h_2 = (2\Gamma - \gamma)\bar{h}_2$$

$$\bar{h}_2 = \Gamma (2\Gamma - \gamma) ((\varepsilon - 1) y_1 + y_2)^2 + \gamma^2 (\varepsilon - 1)^2 z_1^2 - \gamma^2 (2 - 2\varepsilon + \varepsilon^2) z_2 z_1 + \gamma^2 z_2^2$$

$$h_3 = 2(\gamma - \Gamma)\bar{h}_3$$

$$\bar{h}_3 = (\gamma - 2\Gamma) (\gamma + 2\Gamma (\varepsilon - 1)^2) z_2 y_1^2 + (\gamma - 2\Gamma) (\gamma + 2\Gamma) (\varepsilon - 1) (y_2 z_1 + z_2 y_2) y_1 - \gamma^2 \varepsilon (\varepsilon - 2) z_2 z_1^2 + ((\gamma - 2\Gamma) (2\Gamma + (\varepsilon - 1)^2 \gamma) y_2^2 + \gamma^2 \varepsilon (\varepsilon - 2) z_2^2) z_1$$

$$h_4 = 4(\gamma - \Gamma)^2 \bar{h}_4$$

$$\bar{h}_4 = 2 (\gamma + (\varepsilon - 1)^2 \Gamma) z_2^2 y_1^2 + 4 (\varepsilon - 1) (\gamma + \Gamma) z_2 y_2 z_1 y_1 + 2 (\Gamma + (\varepsilon - 1)^2 \gamma) y_2^2 z_1^2$$

$$D' = 2\gamma^2(\Gamma - \gamma)(2\Gamma - \gamma)(1 - \varepsilon)(z_1 - z_2)((\varepsilon - 1) z_1 y_2 + z_2 y_1).$$

In particular we deduce (compare with [5] in the contrast problem):

**PROPOSITION 4.2.** *The quadric $D'$ reduces to a cubic form which is factorized into a linear and a quadratic (homogeneous) forms.*

**4.2.2. Singular analysis.** We assume $\gamma > 0$ and $2\Gamma > 3\gamma$. It implies $\gamma \neq \Gamma$ and $\gamma \neq 2\Gamma$. The main result is the following:

**THEOREM 4.3.** *Provided $\varepsilon \neq 1$ the equilibrium points of $X_e^r = DF_0 - D'F_1$ are all contained in $\{D = D' = 0\}$.*

*Proof.* Obviously, every point of $\{D = 0\} \cap \{D' = 0\}$ is a singularity of $X_e^r$.

Conversely, let us assume $\varepsilon \neq 1$. We first divide $X_e^r$ by $1 - \varepsilon$. We still assume that $\Gamma \neq 0$. We consider the equations $\{(X_e^r)_{y_1} = 0, (X_e^r)_{z_1} = 0, (X_e^r)_{y_2} = 0, (X_e^r)_{z_2} = 0\}$ and remark that the last third are divisible by $\gamma$. By homogeneity, changing $\gamma$ into $\gamma\Gamma$, we get rid of $\Gamma$. So we may assume $\Gamma = 1$. The resulting system is denoted $\Sigma_r$. We add the two polynomials $((\varepsilon - 1) z_1 y_2 + z_2 y_1) a_1 - 1$ and $(z_1 - z_2) a_2 - 1$, and the polynomials $\gamma g - 1$, $(\gamma - 1)g_1 - 1$, $(\gamma - 2)g_2 - 1$. We denote $\tilde{\Sigma}_r$ this new system, involving four new variables $g_1, g_2, a_1, a_2$. We compute a Gröbner basis with total degree with reverse lexicographic order on $(y_1, y_2, z_1, z_2, \varepsilon, g, g_1, g_2, a_1, a_2)$ and get $\{1\}$. Hence, provided $\gamma$ is different from 0, 1, 2, there is no singular point of $X_e^r$ outside of $\{D = 0\} \cap \{D' = 0\}$.

□

The remaining of the section is devoted to the singularity resolution. From the factorized form of $D'$ (Proposition 4.2) we get:

**PROPOSITION 4.4.** *$\{D = 0\} \cap \{D' = 0\}$ is an algebraic variety of algebraic dimension 2 whose components are located in the hyperplane $z_1 = z_2$ and in the hypersurface $(\varepsilon - 1) z_1 y_2 + z_2 y_1 = 0$.*

These components are studied in the following analysis, and explicitly expressed in Lemmas 4.5, 4.6, 4.7, 4.8.

* Case A: components of $\{D = 0\} \cap \{D' = 0\}$ in $z_1 = z_2$.

Under the constraint $z_1 = z_2$, we have a factorization $\tilde{D} = p_1 p_2$ with:

$$p_1 = 2(\gamma - \Gamma)z_1 + \gamma - 2\Gamma$$
and:

$$p_2 = \left( 2(\gamma - \Gamma)\left(\gamma + (\varepsilon - 1)^2\Gamma\right)z_1 + \Gamma(\varepsilon - 1)^2(\gamma - 2\Gamma)\right)y_1^2 + \\ (4(\gamma - \Gamma)(\gamma + \Gamma)(\varepsilon - 1)z_1 + 2\Gamma(\varepsilon - 1)(\gamma - 2\Gamma))y_2y_1 + \\ \left(2(\gamma - \Gamma)\left(\Gamma + (\varepsilon - 1)^2\gamma\right)z_1 + \Gamma(\gamma - 2\Gamma)\right)y_2^2$$

The first polynomial has one root $z_1 = z_{\gamma,\Gamma}$

$$z_{\gamma,\Gamma} = \frac{1}{2} \frac{2\Gamma - \gamma}{\gamma - \Gamma}$$

which corresponds to the plane-solution $\{(y_1, z_{\gamma,\Gamma}, y_2, z_{\gamma,\Gamma}), (y_1, y_2) \in \mathbb{R}^2\}$. We put:

$$d_2(y_1, y_2) = (\gamma + (\varepsilon - 1)^2 \Gamma) y_1^2 + 2(\varepsilon - 1)(\gamma + \Gamma) y_2 y_1 + (\Gamma + (\varepsilon - 1)^2 \gamma) y_2^2$$

The discriminant of $d_2$ with respect to $y_1$ is $-4(\varepsilon - 2)^2 \gamma \Gamma \varepsilon^2 y_2^2$ which is strictly negative provided $\varepsilon \neq 0$. So $d_2$ is non-zero outside $y_1 = y_2 = 0$.

So, provided $y_1^2 + y_2^2 \neq 0$, $d_2 \neq 0$, and $p_2 = 0$ is solved with respect to $z_1$. We get $z_1 = r_2(y_1, y_2)$ with

$$r_2(y_1, y_2) = \frac{\Gamma(2\Gamma - \gamma)((\varepsilon - 1)y_1 + y_2)^2}{2(\gamma - \Gamma)d_2(y_1, y_2)}$$

and $(y_1, r_2(y_1, y_2), y_2, r_2(y_1, y_2))$ (defined for $(y_1, y_2) \neq (0,0)$) vanishes both $D$ and $D'$.

Finally, if $y_1 = y_2 = 0$, we have the solution $(0, z, 0, z)$, $z \in \mathbb{R}$.

We summarize the case $z_1 = z_2$ in:

LEMMA 4.5. $\{D = 0\} \cap \{D' = 0\} \cap \{z_1 = z_2\}$ is the union of an affine plane $z_1 = z_2 = z_{\gamma,\Gamma}$, a rational surface $z_1 = z_2 = r_2(y_1, y_2)$ (defined for $(y_1, y_2) \neq (0,0)$), and the line $\{(0, z, 0, z), z \in \mathbb{R}\}$.

* Case B: components of $\{D = 0\} \cap \{D' = 0\}$ in $(\varepsilon - 1)z_1y_2 + z_2y_1 = 0$.

- Assume first that $y_1 = 0$ and $z_1 \neq z_2$. We have $z_1y_2 = 0$.

* If $y_1 = z_1 = 0$, then:

$$\tilde{D} = (\gamma - 2\Gamma)(\Gamma(\gamma - 2\Gamma)y_2^2 + \gamma^2 z_2^2)$$

Since $2\Gamma > \gamma$, $\{\tilde{D} = 0\} \cap \{y_1 = z_1 = 0\}$ corresponds to two lines intersecting at $N$.

* If $y_1 = y_2 = 0$, then let us put

$$d_1(z_1) = 2\varepsilon(\varepsilon - 2)(\gamma - \Gamma)z_1 + 2\Gamma - \gamma.$$

We have:

$$\tilde{D} = \gamma^2(z_2 - z - 1)(d_1(z_1)z_2 - (\varepsilon - 1)^2(2\Gamma - \gamma)z_1)$$

Observe that the polynomial $d_1$ vanishes if and only if $z_1$ equals $\tilde{z}_{\gamma,\Gamma}$ with

$$\tilde{z}_{\gamma,\Gamma} = \frac{1}{2} \frac{\gamma - 2\Gamma}{\varepsilon(\varepsilon - 2)(\gamma - \Gamma)}$$
and in this case, there is no solution such that $z_2 \neq z_1$.
Provided $d_1(z_1) \neq 0$, one gets $z_2 = r_1(z_1)$:

$$r_1(z_1) = \frac{(\varepsilon - 1)^2 (2\Gamma - \gamma) z_1}{d_1(z_1)}$$

which is a rational function of $z_1$. And the intersection with $\{D = 0\} \cap \{D' = 0\}$ is the curve $\{(0, z_1, 0, r_1(z_1)) \ z_1 \in \mathbb{R} \setminus \{\tilde{z}_{\gamma, \Gamma}\}\}$.

LEMMA 4.6. $\{D = 0\} \cap \{D' = 0\} \cap \{y_1 = 0\} \cap \{(z_1 - z_2) \neq 0\}$ is the union of two lines of $\{y_1 = z_1 = 0\}$ intersecting at $N$ and a rational curve $\{(0, z_1, 0, r_1(z_1)) \ z_1 \in \mathbb{R} \setminus \{\tilde{z}_{\gamma, \Gamma}\}\}$.

– Let us assume $y_1 \neq 0$.

We can eliminate $z_2$ using:

$$z_2 = \frac{z_1 y_2 (1 - \varepsilon)}{y_1}$$

and, substituting in $y_1^2 \tilde{D}$ we get the factorization $y_1^2 \tilde{D} = q_1 q_2$, with:

$$q_1 = \Gamma (\varepsilon - 1) (2\Gamma - \gamma) y_1^3 + \gamma^2 (\varepsilon - 1) z_1^2 y_1 + \gamma^2 (\varepsilon - 1)^2 z_1^2 y_2 \\ + (2\Gamma\varepsilon (\varepsilon - 2) (\gamma - \Gamma) z_1 - \Gamma (\gamma - 2\Gamma)) y_2 y_1^2$$

and:

$$q_2 = (\varepsilon - 1)(\gamma - 2\Gamma)y_1 + (2\varepsilon(2-\varepsilon)(\gamma - \Gamma)z_1 + \gamma - 2\Gamma)y_2 \\ = (\varepsilon - 1)(\gamma - 2\Gamma)y_1 + d_1(z_1)y_2$$

Provided $d_1 \neq 0$ (that is $z_1 \neq \tilde{z}_{\gamma, \Gamma}$), we solve $q_2 = 0$ with respect to $y_2$, and then we get the value of $(y_2, z_2)$:

$$\left( \frac{(\varepsilon - 1)(\gamma - 2\Gamma) y_1}{d_1(z_1)}, \frac{(\varepsilon - 1)^2 (2\Gamma - \gamma) z_1}{d_1(z_1)} \right)$$

LEMMA 4.7. $\{D = 0\} \cap \{D' = 0\} \cap \{(z_1 - z_2) y_1 d_1(z_1) \neq 0\}$ is a rational surface $(y_2 = \rho_2(y_1, z_1), z_2 = \rho_1(z_1)y_1 \neq 0)$ with $z_1 \neq \tilde{z}_{\gamma, \Gamma}$.
We put $d_3$

$$d_3 = (2\Gamma\varepsilon(\varepsilon - 2)(\gamma - \Gamma)z_1 - \Gamma(\gamma - 2\Gamma))y_1^2 + \gamma^2(\varepsilon - 1)^2 z_1^2$$

Its discriminant with respect to $y_1$ is:

$$-4 (2\Gamma - 4\gamma z_1\varepsilon + 2\gamma z_1\varepsilon^2 - \gamma + 4\Gamma z_1\varepsilon - 2\Gamma z_1\varepsilon^2)\Gamma\gamma^2 z_1^2 (\varepsilon - 1)^2 \\ -4 (2\Gamma - \gamma + 2\varepsilon(2-\varepsilon)(\Gamma - \gamma)z_1)\Gamma\gamma^2 z_1^2 (\varepsilon - 1)^2$$

and its sign changes when $z_1$ reaches $\tilde{z}_{\gamma, \Gamma}$.

Provided $d_3(y_1, z_1) \neq 0$, we solve $q_1$ with respect to $y_2$, and then we get the value of $(y_2, z_2)$:

$$\left( \frac{(\Gamma(2\Gamma - \gamma)y_1^2 + \gamma^2 z_1^2)(1 - \varepsilon)y_1}{d_3(y_1, z_1)}, \frac{(\Gamma(2\Gamma - \gamma)y_1^2 + \gamma^2 z_1^2)(\varepsilon - 1)^2 z_1}{d_3(y_1, z_1)} \right)$$

LEMMA 4.8. $\{D = 0\} \cap \{D' = 0\} \cap \{(z_1 - z_2) y_1 d_3(z_1) \neq 0\}$ is a rational surface with parameterization $(y_2 = \rho_3(y_1, z_1), z_2 = \rho_4(y_1, z_1))$.
• Analysis of the behaviours of the solutions of $X_e^r$ near O.

We set $\tilde{z}_i = 1 + z_i$ and we have the following approximations:

$$-D = (1 - \varepsilon)\tilde{D}, \quad \tilde{D} = h_1 + h_2,$$

$$h_1 = \gamma^2 \varepsilon (\varepsilon - 2) (\gamma - 2\Gamma) (\tilde{z}_1 - \tilde{z}_2)$$

$$h_2 = \Gamma (\varepsilon - 1)^2 (\gamma - 2\Gamma)^2 y_1^2 + 2\Gamma (\gamma - 2\Gamma)^2 (\varepsilon - 1) y_2 y_1 + \Gamma (\gamma - 2\Gamma)^2 y_2^2 - \gamma^2 (\varepsilon - 1)^2 (\gamma - 2\Gamma) \tilde{z}_1^2 - \gamma^2 (\gamma - 2\Gamma) \tilde{z}_2^2 + \gamma^2 (\varepsilon^2 + 2 - 2\varepsilon) (\gamma - 2\Gamma) \tilde{z}_1 \tilde{z}_2$$

$$-D' = 2\gamma^2(\Gamma - \gamma)(2\Gamma - \gamma)(1 - \varepsilon)(\tilde{z}_2 - \tilde{z}_1)[(-1 + \tilde{z}_1)y_2(\varepsilon - 1) + (-1 + \tilde{z}_2)y_1].$$

**Conclusion:** these computations allow to evaluate the equilibrium points and the behaviors of the solutions near such point, using linearization methods. A first step towards the global behavior is the following result.

**LEMMA 4.9.** *The surface $y_1 = y_2 = 0$ is foliated by lines solutions connecting O to the north pole N, the singular control being zero.*

**4.3. Simulations for $N = 2$ spins.** For a number of $N$ spins, the multisaturation problem in the single-input case can be summarized in this way:

$$(\text{MSS}) \qquad \left\{ \begin{array}{l} c(q(t_f)) = \frac{1}{N} \sum_{i=1,N} |q_i(t_f)|^2 \rightarrow \min_{u_1(\cdot)} t_f \text{ fixed} \\ \dot{q} = F_0(q) + u_1 F_1(q), \quad q(0) = q_0 \end{array} \right.$$

where $q = (q_1, \dots, q_N)$, $q_i = (y_i, z_i) \in \mathbb{R}^2$, $|q_i| \le 1$, $i = 1, \dots, N$. The initial condition for each spin is $q_i(0) = (0, 1)$. The vector field $F_0$ and $F_1$ are given by:

$$F_0(q) = \sum_{i=1,N} (-\Gamma y_i) \frac{\partial}{\partial y_i} + (\gamma(1-z_i)) \frac{\partial}{\partial z_i}, \quad F_1(q) = \sum_{i=1,N} a_i \left( -z_i \frac{\partial}{\partial y_i} + y_i \frac{\partial}{\partial z_i} \right),$$

where $\Lambda = (\gamma, \Gamma)$ are the physical parameters and $a = (a_i)_{i=1,N}$ is the set of positive scaling factors.

**Numerical results.** We present now the simulations made for $N = 2$ spins. As for the previous contrast problem, we first use the direct method from Bocop to obtain a rough solution that is refined by multiple shooting with HamPath. We take the blood case, with $(T_1, T_2) = (1.35, 0.05)$ for each spin. The final time is set to $t_f = \lambda T_{\min}$ where $T_{\min} = 6.7981$ is the minimum time to bring one spin to saturation. The scaling factors are taken uniformly in $[1 - \varepsilon_{\max}, 1]$, i.e. for $N = 2$ we have $(a_1, a_2) = (1, 1 - \varepsilon_{\max})$. The direct method uses an implicit midpoint 2nd order discretization with 250 steps, with all state variables initialized to 0 and all control variables initialized to 0.1. The indirect methods uses the data from the direct solution to initialize the multiple shooting. Fig. 12 and Table. 4.3 illustrate the solutions from both methods for $\lambda = 1.1$, $\varepsilon_{\max} = 0.3$. We observe that the control presents a 2BS structure with three switching times $t_1, t_2, t_3$, and that the solution from the direct method is quite close to the one refined by the indirect method.

**Homotopy on transfert time (parameter $\lambda$).** Starting from the first solutions obtained by the direct method, we now study with HamPath a homotopy for $\lambda \in [1, 2]$. The norm of the shooting function $S$ while increasing, stays below $10^{-10}$. Fig. 13 shows the evolution of the mean distance $(\frac{1}{N} \sum_{i=1}^N |q_i(t_f)|)$ of the spins to
<table><caption>TABLE 3

Blood case, $N = 2$, $\lambda = 1.1$, $\varepsilon_{\max} = 0.3$: comparison of the direct and indirect methods.</caption><thead><tr><th>Method</th><th>c(q(t<sub>f</sub>))</th><th>initial costate p(0)</th><th>t<sub>1</sub>, t<sub>2</sub>, t<sub>3</sub></th><th>CPU</th></tr></thead><tbody><tr><td>Bocop</td><td>2.56 &times; 10<sup>-4</sup></td><td>(2.57, -3.46, -2.96, -4.92) &times; 10<sup>-4</sup></td><td>(0.043, 0.338, 0.382)</td><td>0.66s</td></tr><tr><td>HamPath</td><td>2.56 &times; 10<sup>-4</sup></td><td>(1.82, -3.85, -3.47, -4.48) &times; 10<sup>-4</sup></td><td>(0.040, 0.347, 0.385)</td><td>0.43s</td></tr></tbody></table>

<figure><img src="image_15.png" /><figcaption>FIG. 12. Blood case with $N = 2$ spins, $\lambda = 1.1$, $\varepsilon_{\max} = 0.3$. Solution from direct (dots) and indirect (line) methods. Control structure is 2BS. The rough solution from the direct method is actually quite close to the refined solution from the indirect method.</figcaption></figure>

the origin for $\lambda \in [1, 2]$. We observe a linear decrease in log scale, which suggests a law of the form $D = e^{-C\lambda}$. We also display the evolution of the switching times, the control structure being 2BS over the whole homotopy. The decreasing duration of the first singular arc may indicate that this branch of solutions will stop or change at a certain limit value for $\lambda$. Fig. 14 shows the control for $\lambda = 1, 1.4989$ and 2.

<figure><img src="image_16.png" /><figcaption>FIG. 13. Homotopy on transfer time (parameter $\lambda$). The mean distance to the origin (left) decreases linearly in log scale. The switching times (right) indicate a decreasing duration $t_2 - t_1$ for the first singular arc.</figcaption></figure>

**Non-envelopping trajectories.** As we consider a set of $N$ spins with a $B_1$-inhomogeneity term $\varepsilon$ varying in $[0, \varepsilon_{\max}]$, we can wonder if the optimal solution is actually determined by the extremal cases alone, namely by taking only the two spins for $\{0, \varepsilon_{\max}\}$. To gain an insight on this point, we take the optimal control from the solution for $N = 2$ spins, and add a third spin in the middle with $\varepsilon = \varepsilon_{\max}/2$, see Fig. 15. We observe that the trajectory for the third spin is not envelopped by the two original spins, and that its final position is much farther from the origin. Furthermore, this third spin does not get closer to saturation as the transfer time increases, contrary to the two original spins. This indicates that a proper solution of the multisaturation problem is likely to require optimization with at least a significant number $N$ of spins.
<figure><img src="image_17.png" /><figcaption>FIG. 14. *Homotopy on transfer time (parameter λ): controls.* For λ = 1, 1.4989 and 2, the control structure remains 2BS, with the duration of the first singular arc decreasing for larger λ.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>FIG. 15. Case $N = 2 + 1$ spins, $\lambda = 1.1$ and homotopy on $\lambda$. Trajectories for all spins (left), zoom near the origin (center): the trajectory for the third spin added after optimization (dashed blue lines) is not envelopped by the two original spins (solid red lines), and its final position is much farther from the origin. Final distance to the origin (right) w.r.t. $\lambda$: the distance for the third spin does not decrease for longer transfer times.</figcaption></figure>

**Remark on conjugate times.** For every solution computed by the indirect shooting in HamPath, we check that there is no conjugate time over the singular arcs. Namely, all singular arcs satisfy this necessary condition for local optimality.

<figure><img src="image_19.png" /><figcaption>FIG. 16. Multisaturation with $N = 11$ spins: fluid case, $\epsilon_{max} = 0.3$. On the right, the trajectory of all spins for $\lambda = 1$ i.e. $t_f = T_{min}$. On the right, a closeup on the final positions of all spins for $\lambda = 1, 1.125$ and $1.25$. We observe that the spins tend to spread regularly around the origin, and get closer for larger transfer times.</figcaption></figure>

**5. Conclusion.** In this article, we have presented theoretical results to analyze the contrast problem and the multisaturation problem for an ensemble of two spins.
They are completed by numerical simulations using both direct and indirect methods and lmi techniques.

They are a preliminary step to a complete understanding of the contrast problem for a theoretical point of view and numerical simulations using the direct method to compensate the B₀ and B₁ inhomogeneities. See Fig. 16 for a preliminary result in the case of B₁ inhomogeneity.

REFERENCES

[1] F. J. Bonnans, P. Martinon & V. Grélard. *Bocop - A collection of examples*, Technical report, INRIA, 2012. RR-8053.

[2] B. Bonnard, J.-B. Caillau & E. Trélat, *Second order optimality conditions in the smooth case and applications in optimal control*, ESAIM Control Optim. Calc. Var., **13** (2007), no. 2, 207-236.

[3] B. Bonnard, J.-B. Caillau & E. Trélat, *Geometric optimal control of elliptic Keplerian orbits*, Discrete Contin. Dyn. Syst. Ser. B 5, (2005), no. 4, 929-956.

[4] B. Bonnard & M. Chyba, *Singular trajectories and their role in control theory*, vol. **40** of Mathematics & Applications, Springer-Verlag, Berlin (2003), xvi+357.

[5] B. Bonnard, M. Chyba, A. Jacquemard & J. Marriott, *algebraic geometric classification of the singular flow in the contrast imaging problem in nuclear magnetic resonance*, Math. Control Relat. Fields, **3** (2013), no. 4, 397-432.

[6] B. Bonnard, M. Chyba & J. Marriott, *Singular Trajectories and the Contrast Imaging Problem in Nuclear Magnetic resonance*, SIAM J. Control Optim., **51** (2013), no. 2, 1325-1349.

[7] B. Bonnard, M. Claeys, O. Cots & P. Martinon, *Geometric and numerical methods in the contrast imaging problem in nuclear magnetic resonance*, Acta Appl. Math., (2013), to appear.

[8] B. Bonnard & O. Cots, *Geometric numerical methods and results in the control imaging problem in nuclear magnetic resonance*, Math. Models Methods Appl. Sci., **24** (2014), no. 1, 187-212.

[9] B. Bonnard, O. Cots, S. Glaser, M. Lapert, D. Sugny & Y. Zhang, *Geometric optimal control of the contrast imaging problem in nuclear magnetic resonance*, IEEE Trans. Automat. Control, **57** (2012), no 8, 1957-1969.

[10] J.-B. Caillau & B. Daoud, *Minimum time control of the circular restricted three-body problem*, SIAM J. Control Optim., **50** (2011), no. 6, 3178-3202.

[11] Y. Chitour, F. Jean & E. Trélat, *Genericity results for singular curves*, Journal of Differential Geometry, **73** (2006), no. 1, 45-73.

[12] O. Cots, *Contrôle optimal géométrique : méthodes homotopiques et applications*. Phd thesis, Institut Mathématiques de Bourgogne, Dijon, France, 2012.

[13] D. Henrion, J. B. Lasserre & J. Löfberg, *GloptiPoly 3: Moments, Optimization and Semidefinite Programming*, Optim. Methods and Software, **24** (2009), no. 4-5, 761-779.

[14] J. S Li & N. Khaneja, *Control of inhomogeneous quantum ensembles*, Phys. Rev. A. 73, 030302 (2006).

[15] N. Khaneja, T. Reiss, C. Kehlet, T. Schulte-Herbruggen & S. J. Glaser, *Optimal control of coupled spin dynamics: design of nmr pulse sequences by gradient ascent algorithms*, J. Magn. Reson., **172** (2005), no 2, 296-305.

[16] I. Kupka, *Geometric theory of extremals in optimal control problems. i. the fold and Maxwell case*, Trans. Amer. Math. Soc., **299** (1987), no. 1, 225-243.

[17] M. Lapert, Y. Zhang, S. J. Glaser & D. Sugny, *Towards the time-optimal control of dissipative spin-1/2 particles in nuclear magnetic resonance*, J. Phys. B: At. Mol. Opt. Phys., **44** (2011), no. 15.

[18] M. Lapert, Y. Zhang, M. Janich, S. J. Glaser & D. Sugny, *Exploring the physical limits of saturation contrast in Magnetic Resonance Imaging Sci. Rep.* 2, 589 (2012).

[19] J. B. Lasserre, D. Henrion, C. Prieur & E. Trélat, *Nonlinear optimal control via occupation measures and LMI-relaxations*, SIAM J. Control Optim., **47**, (2008), no. 4, 1643-1666.

[20] M. H. Levitt, *Spin dynamics: basics of nuclear magnetic resonance*, John Wiley & Sons, 2001.

[21] J. S Li, *Control of inhomogeneous ensembles*, Phd dissertation, Harvard University, 2006.

[22] L. S. Pontryagin, V. G. Boltyanski, R. V. Gamkrelidze & E. F. Mishchenko, *Matematicheskaya teoriya optimalnykh protsessov*, "Nauka", Moscow, fourth ed. (1983).