Bernard Bonnard
email: bernard.bonnard@u-bourgogne.fr.

Mathieu Claeys
email: mclaeys@laas.fr.

Olivier Cots
email: olivier.cots@inria.fr.

Pierre Martinon
email: marti-non@cmap.polytechnique.fr.

Geometric and numerical methods in the contrast imaging problem in nuclear magnetic resonance

Keywords: Geometric optimal control, Contrast imaging in NMR, Direct method, Shooting and continuation techniques, Moment optimization

In this article, the contrast imaging problem 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 techniques of geometric control. This leads to a numerical investigation based on so-called indirect methods using the HamPath software. The results are then compared with a direct method implemented within the Bocop toolbox. Finally lmi techniques are used to estimate a global optimum.

Introduction

The control in Nuclear Magnetic Resonance (NMR) of a spin-1/2 particle is described by the Bloch equation [START_REF] Levitt | Spin dynamics : basics of nuclear magnetic resonance[END_REF] dM x dτ = -M x /T 2 + ω y M z -∆ ωM y

dM y dτ = -M y /T 2 -ω x M z + ∆ ωM x dM z dτ = (M 0 -M z )/T 1 + ω x M y -ω y M x (1)
where the state variable is the magnetization vector M = (M x , M y , M z ), the control is the magnetic field ω = (ω x , ω y , 0) bounded by |ω| ≤ ω max , ω max = 2π × 32.3Hz being the experimental intensity of the experiments [START_REF] Lapert | Exploring the physical limits of saturation contrast in Magnetic Resonance Imaging[END_REF], ∆ ω is the resonance offset and τ is the time. In order to normalize our system, we introduce: q = (x, y, z) = (M x , M y , M z )/M 0 and q belongs to the Bloch ball |q| ≤ 1. The time is normalized to t = τ ω max and we introduce the control u = ω/ω max , |u| ≤ 1. In this paper, we assume ∆ ω = 0, leading to the following normalized system:

dx dt = -Γ x + u 2 z dy dt = -Γ y -u 1 z dz dt = γ(1 -z) + u 1 y -u 2 x, (2)
where the parameters Γ = 1/(ω max T 2 ) and γ = 1/(ω max T 1 ) satisfy 2Γ ≥ γ ≥ 0.

In the contrast problem, we consider two uncoupled spin-1/2 systems corresponding to different particles, each of them solution of the Bloch equation ( 2) with respective coefficients (γ 1 ,Γ 1 ), (γ 2 ,Γ 2 ), and controlled by the same magnetic field. By denoting each system by qi = f (q i ,Λ i , u), Λ i = (γ i ,Γ i ) and q i = (x i , y i , z i ), this leads to consider the system written shortly as q = F(q, u), where q = (q 1 , q 2 ) is the state variable and where q means the derivative of q with respect to the time t. In the numerical simulations, we shall consider two cases encountered in the experiments, the physical parameters being the relaxation times given in seconds.

Definition 1

The contrast problem by saturation is the following optimal control problem (OCP): starting from the equilibrium point q 0 = ((0, 0, 1) , (0, 0, 1)), reach in a given transfer time t f the fixed state q 1 (t f ) = 0 (corresponding to the zero magnetization of the first spin, called saturation) while maximising |q 2 (t f )| 2 , the contrast being |q 2 (t f )|.

We are interested in this article in solving the contrast problem by saturation for a range of values of t f . The contrast problem can be stated as a Mayer problem given by the following smooth conditions:

1. A system q = F(q, u), q ∈ Q ⊂ R n , with

Q = {q = (q 1 , q 2 ) ∈ R n : |q 1 | ≤ 1, |q 2 | ≤ 1},
and fixed initial state q(0) = q 0 . The control belongs to the control domain U = {u ∈ R : |u| ≤ 1}.

2. A terminal manifold to reach,

Q f = {q = (q 1 , q 2 ) ∈ R n : q 1 = 0, |q 2 | ≤ 1} ⊂ Q,
defined by f (q) = 0, where f : R n → R k . 3. A cost to be minimized of the form min u(•) c(q(t f )) where c : R n → R is a regular mapping and t f is the final time.

The contrast problem by saturation can be summarized this way:

(OCP)

       c(q(t f )) = -|q 2 (t f )| 2 -→ min u(•) , t f fixed q = F(q, u), q(0) = q 0 q 1 (t f ) = 0 (OCP)
In practical experiments we consider two cases:

a) The bi-input case where q = (q 1 , q 2 ) ∈ R 6 ∩ Q and |u| = u 2 1 + u 2 2 ≤ 1. b) The single-input case where the system is restricted to x 1 = x 2 = 0, the control field is restricted to the real field, i.e., u 2 = 0, and each spin is restricted to the plane q i = (y i , z i ).

The contrast imaging problem was analyzed in a series of recent articles using geometric optimal control [START_REF] Bonnard | Singular trajectories and their role in control theory[END_REF][START_REF] Bonnard | algebraic geometric classification of the singular flow in the contrast imaging problem in nuclear magnetic resonance[END_REF] or numerical indirect methods [START_REF] Bonnard | Geometric numerical methods and results in the control imaging problem in nuclear magnetic resonance[END_REF], and leads to experimental work described in [START_REF] Bonnard | Geometric optimal control of the contrast imaging problem in nuclear magnetic resonance[END_REF][START_REF] Lapert | Exploring the physical limits of saturation contrast in Magnetic Resonance Imaging[END_REF], but restricted mainly to the single-input case. In such studies, the contrast problem is essentially reduced to the analysis of the socalled singular trajectories of the system modeling the problem. The objective of this article is to complete the geometric and the numerical analysis of the problem. Concerning the geometric analysis, we have mainly consider in our previous works the case where the phase of the RF magnetic field is zero, and thus u 2 = 0. One objective of this article will be to fill the gap between the bi-input and single-input studies. From the numerical point of view, the indirect method using the HamPath software will be complemented with the Bocop toolbox [5] based on the direct approach.

A feature of the contrast problem is the existence of many locally optimal solutions which can be computed by the previously mentioned direct or indirect methods. An important question is then to assert global optimality. Hence another objective of this article is to use a moment/linear matrix inequality (LMI) technique [START_REF] Henrion | GloptiPoly 3: Moments, Optimization and Semidefinite Programming[END_REF][START_REF] Lasserre | Positive polynomials and their applications[END_REF] to compute such an estimate using the GloptiPoly toolbox [START_REF] Henrion | GloptiPoly 3: Moments, Optimization and Semidefinite Programming[END_REF].

This article is organized in four sections. The first section recalls and completes the geometric framework. The second section presents in details the three numerical methods. The numerical results are described and compared in the two final sections.

Necessary optimality conditions, geometric framework and results

The aim of this section is to recall the theoretical foundations. We use [START_REF] Bonnard | Singular trajectories and their role in control theory[END_REF] as a reference in geometric control, whereas the application to the contrast problem is excerpted from [START_REF] Bonnard | Singular Trajectories and the Contrast Imaging Problem in Nuclear Magnetic resonance[END_REF].

Maximum principle

Proposition 1 If u * with corresponding trajectory q * on [0,t f ] is optimal then the following necessary optimality conditions are satisfied. Define the pseudo-Hamiltonian H(q, p, u) = p, F(q, u) , p being the adjoint vector, there exists p * (•) such that for almost every t ∈ [0,t f ], we have

(i) q * = ∂ H ∂ p (q * , p * , u * ), ṗ * = - ∂ H ∂ q (q * , p * , u * )
(ii) H(q * , p * , u * ) = max v∈U H(q * , p * , v) (maximization condition)

and the following boundary conditions

(iii) f (q * (t f )) = 0 (iv) p * (t f ) = p 0 ∂ c ∂ q (q * (t f )) + k ∑ i=1 σ i ∂ f i ∂ q (q * (t f )) (transversality condition), σ = (σ 1 , . . . , σ k ) ∈ R k , p 0 ≤ 0
Definition 2 The solutions of conditions (i) and (ii) of Prop. 1 are called extremals and BC-extremals if they satisfy the boundary conditions. Note that M(q * , p * ) = max v∈U H(q * , p * , v) is constant and if M(q * , p * ) = 0, an extremal is called exceptional.

1.2 Application to the contrast problem

State space

Since the Bloch ball is invariant for the dynamics of each spin particle, the state constraints are not active and the maximum principle can be applied.

Boundary conditions

In the contrast problem, q = (q 1 , q 2 ), f = 0 is the set q 1 = 0, and the cost to minimize is c(q(t f )) = -|q 2 (t f )| 2 . Hence, splitting the adjoint vector into p = (p 1 , p 2 ), we deduce from the transversality condition that p 2 (t f ) = -2p 0 q 2 (t f ), p 0 ≤ 0. If p 0 is nonzero, it can be normalized to p 0 = -1/2.

Parameterization of the extremal curves

Bi-input case

The system is written as q = F 0 (q) + u 1 F 1 (q) + u 2 F 2 (q), |u| ≤ 1. The maximization condition in Prop. 1 leads then to the following parameterization of the extremal controls:

u 1 = H 1 H 2 1 + H 2 2 , u 2 = H 2 H 2 1 + H 2 2 , (3)
where H i = p, F i (q) are Hamiltonian lifts. Define the switching surface:

Σ : H 1 = H 2 = 0.
Plugging such a u into the pseudo-Hamiltonian gives the true Hamiltonian H n =

H 0 + H 2 1 + H 2 2 .
Outside Σ , the smooth solutions of the corresponding vector field #-H n are called extremals of order zero. Let (q, p, u) be a normal order zero extremal, using polar coordinates u 1 = cos α, u 2 = sin α. If α ≡ constant then the extremal has a stationary phase control. In the other case we say a non-stationary phase control.

Single-input case

Consider the case where u 2 = 0 and the system is written q = F 0 (q) + u 1 F 1 (q), where q belongs to a 4-dimensional space Q and |u 1 | ≤ 1. The application of the maximization condition yields are two types of extremals:

-Regular extremals: The control is given by u 1 (t) = sgn (H 1 (z(t))), where z = (q, p) and t ∈ [0,t f ]. If the number of switchings is finite, it is called bang-bang. -Singular extremals: The singular case is the situation where H 1 (z(t)) = 0. The singular control can be computed as follows. If F, G are two vector fields, the Lie bracket is:

[F, G](q) = ∂ F ∂ q (q)G(q) - ∂ G ∂ q (q)F(q). If #- H F ,
#-H G are two Hamiltonian vector fields the Poisson bracket is:

{H F , H G }(z) = dH F ( #- H G )(z). If H F = p, F(q) , H G = p, G(q) then we have: {H F , H G }(z) = p, [F, G](q) .
Applied to the contrast problem and differentiating twice H 1 (z(t)) = 0 with respect to time leads to:

{H 1 , H 0 } = {{H 1 , H 0 }, H 0 } + u 1 {{H 1 , H 0 }, H 1 } = 0.
From the second condition, we derive when the denominator is not vanishing the corresponding singular control

u 1,s = - {{H 1 , H 0 }, H 0 } {{H 1 , H 0 }, H 1 } . ( 4
)
We denote Σ 1 :

H 1 = 0 and Σ ′ 1 : H 1 = {H 1 , H 0 } = 0.
Plugging such u 1,s into the pseudo-Hamiltonian defines the Hamiltonian vector field #-H s . Singular extremals such that {{H 1 , H 0 }, H 1 } = 0 are called of order two and they are the solutions of #-H s starting at t = 0 from the surface Σ ′ 1 and satisfying |u 1,s | ≤ 1. Optimal solutions of the contrast problem are concatenations of bang and singular extremals. For the following sections, we introduce some notations. We note by BS the sequence composed by one bang arc (δ + or δ -) followed by one singular arc (δ s ), both non empty, and nBS, n > 1, the concatenation of n BS-sequences, all non empty. The complexity of a sequence is measured by the number of concatenated arcs. By noticing that in our cases one has Γ i , γ i ≪ 1, i = 1, 2, this leads to the following lemma: Lemma 1 The simplest BC-extremal in the contrast problem is of the form δ + δ s , and so has a BS structure.

Generalized Legendre-Clebsch condition

In the singular case the higher-order maximum principle [START_REF] Krener | The high order maximal principle and its application to singular extremals[END_REF] is used to complement Prop. 1 and we have: Proposition 2 In the singular case, the generalized Legendre-Clebsch condition

∂ ∂ u 1 d 2 dt 2 ∂ H ∂ u 1 = {H 1 , {H 1 , H 0 }} ≤ 0 (5)
is a necessary condition for optimality in the Mayer problem.

Saturation problem

The saturation problem consists in bringing the first spin to zero and is analyzed in details in [START_REF] Bonnard | Singular Trajectories and the Contrast Imaging Problem in Nuclear Magnetic resonance[END_REF]25]. We recall briefly the results, denoting (y 1 , z 1 ) the coordinates of the first spin. The singular trajectories are the two lines the z 1 -axis of revolution y 1 = 0 and the horizontal line

z 1 = γ 1 /2(γ 1 -Γ 1 ). Along this horizontal direction, the flow is ẏ1 = -Γ 1 y 1 - γ 2 1 (2Γ 1 -γ 1 ) 4(γ 1 -Γ 1 ) 2 y 1
. Consider the physical situation where 2Γ 1 > 3γ 1 and the horizontal line is such that 0 > z 1 > -1. The optimal policy to steer the north pole to zero is of the form δ + δ h s δ + δ v s , starting from an arc δ + (where u = +1) whereas the second arc δ + connects the horizontal arc δ h s to δ v s . The corresponding minimal time is denoted T min , and in the contrast by saturation one must have t f ≥ T min . An estimate of T min can be obtained as follows. Assuming ω max unbounded, the time to steer y 1 = -1 to y 1 = 0 along the singular horizontal line is according to the dynamics

τ 1 = 1 0 dy 1 Γ 1 y 1 + γ 2 1 (2Γ 1 -γ 1 ) 4(γ 1 -Γ 1 ) 2 y 1 ,
which can be easily computed. If τ 2 is the time to steer (0, z 1 ) to (0, 0) along the vertical singular line, then T min ≃ τ 1 + τ 2 .

In the limit case t f = T min in the contrast problem by saturation, the optimal policy is a 2BS structure. We have the following proposition.

Proposition 3

The saturation solution of the first spin system can be embedded as an extremal solution of the contrast problem by saturation with p 0 = 0 in the transversality condition of Prop. 1.

A geometric property

The following property is crucial: Property 1 Consider the rotation of angle θ with respect to the z i -axis written as exp(θ R) and defined by

X i = x i cos θ + y i sin θ , Y i = -x i sin θ + y i cos θ , Z i = z i ( 6
)
where i = 1, 2. Then, the system is invariant up to the rotation of the control

v 1 = u 1 cos θ + u 2 sin θ , v 2 = -u 1 sin θ + u 2 cos θ . (7)
Clearly, the boundary conditions in the contrast problem are invariant for the rotation, and we have the following proposition (the first assertion being due to Noether theorem in Hamiltonian form)

Proposition 4
We have:

i) The Hamiltonian lift H R (z) = p, Rq is a first integral along the optimal solution. ii) Starting from the point N = ((0, 0, 1), (0, 0, 1)), there is a one-parameter family of optimal solutions, each deduced by the θ -rotation (6), ( 7) respectively on the state and on the control.

Classification of the extremal curves

The geometric control methods are efficient to stratify the set of extremals, and they are complemented by the geometric analysis of the extremal curves. This is explained next based on preliminary works [7,16] for the first point and [START_REF] Bonnard | algebraic geometric classification of the singular flow in the contrast imaging problem in nuclear magnetic resonance[END_REF] for the second point.

The framework

The system is written as

q = F 0 (q) + u 1 F 1 (q) + u 2 F 2 (q)
and an extremal control of order zero is given by:

u = (H 1 , H 2 ) H 2 1 + H 2 2 .
This can be directly plugged into the Hamiltonian lift H 0 + u 1 H 1 + u 2 H 2 and leads to the Hamiltonian:

H n = H 0 + H 2 1 + H 2 2
, whose smooth solutions are the extremals of order zero. The associated differential equation can be regularized using the reparameterization ds = dt

H 2 1 + H 2 2 .
This leads to analyze a dynamical system with non isolated singularities located in the surface Σ : H 1 = H 2 = 0. This approach allows to investigate globally the extremal flow of order zero and the construction of non smooth solutions formed by connecting extremals of order zero at point of Σ , provided the Weierstrass-Erdmann conditions are satisfied:

p(t + ) = p(t -), H n (t + ) = H n (t -),
where t is the time of contact with the switching surface. See [START_REF] Bonnard | algebraic geometric classification of the singular flow in the contrast imaging problem in nuclear magnetic resonance[END_REF] for the related analysis in the single-input case.

Spherical coordinates

Introducing for each spin system the spherical coordinates:

x i = ρ i sin ϕ i cos θ i , y i = ρ i sin ϕ i sin θ i , z i = ρ i cos ϕ i ,
the system takes the form

dρ i dt = γ i cos ϕ i -ρ i (δ i cos 2 ϕ i + Γ i ) dϕ i dt = -γ i sin ϕ i ρ i + δ i sin ϕ i cos ϕ i + (-u 1 sin θ i + u 2 cos θ i ) dθ i dt = -cot ϕ i (u 1 cos θ i + u 2 sin θ i ),
where

δ 1 = γ 1 -Γ 1 , i = 1, 2. Setting θ ′ 1 = θ 1 , θ ′ 2 = θ 1 -θ 2 , one gets: dρ 1 dt = γ 1 cos ϕ 1 -ρ 1 (δ 1 cos 2 ϕ 1 + Γ 1 ) dϕ 1 dt = -γ 1 sin ϕ 1 ρ 1 + δ 1 sin ϕ 1 cos ϕ 1 + (-u 1 sin θ ′ 1 + u 2 cos θ ′ 2 ) dθ ′ 1 dt = -cot ϕ 1 (u 1 cos θ ′ 1 + u 2 sin θ ′ 1 ) dρ 2 dt = γ 2 cos ϕ 2 -ρ 2 (δ 2 cos 2 ϕ 2 + Γ 2 ) dϕ 2 dt = -γ 2 sin ϕ 2 ρ 2 + δ 2 sin ϕ 2 cos ϕ 2 -u 1 sin(θ ′ 1 -θ ′ 2 ) + u 2 cos(θ ′ 1 -θ ′ 2 ) dθ ′ 2 dt = -cot ϕ 1 (u 1 cos θ ′ 1 + u 2 sin θ ′ 1 ) + cot ϕ 2 (u 1 cos(θ ′ 1 -θ ′ 2 ) + u 2 sin(θ ′ 1 -θ ′ 2 )).
Setting

v 1 = u 1 cos θ ′ 1 +u 2 sin θ ′ 1 , v 2 = -u 1 sin θ ′ 1 +u 2 cos θ ′ 2 the system takes the form dρ 1 dt = γ 1 cos ϕ 1 -ρ 1 (δ 1 cos 2 ϕ 1 + Γ 1 ) dϕ 1 dt = -γ 1 sin ϕ 1 ρ 1 + δ 1 sin ϕ 1 cos ϕ 1 + v 2 dθ ′ 1 dt = -cot ϕ 1 v 1 dρ 2 dt = γ 2 cos ϕ 2 -ρ 2 (δ 2 cos 2 ϕ 2 + Γ 2 ) dϕ 2 dt = -γ 2 sin ϕ 2 ρ 2 + δ 2 sin ϕ 2 cos ϕ 2 + v 1 sin θ ′ 2 + v 2 cos θ ′ 2 dθ ′ 2 dt = -cot ϕ 1 v 1 + cot ϕ 2 (v 1 cos θ ′ 2 -v 2 sin θ ′ 2 ).
Observe that the Hamiltonian H n is invariant for a θ -rotation on the set of controls. A diffeomorphism on the state Q = f (q) induces a Mathieu symplectic transformation:

Q = f (q), p = t ∂ f ∂ q (q)P.
In particular, the transformation:

θ ′ 1 = θ 1 , θ ′ 2 = θ 1 -θ 2 in
duces on the adjoint variable the transformation:

p θ 1 = p θ ′ 1 + p θ ′ 2 , p θ 2 = -p θ ′ 2
In the spherical coordinates the Hamiltonian H n is computed with:

H 2 1 + H 2 2 =p 2 ϕ 1 + p 2 ϕ 2 + 2(p ϕ 1 p ϕ 2 + p θ 1 p θ 2 cot ϕ 1 cot ϕ 2 ) cos(θ 1 -θ 2 )+ p 2 θ 1 cot 2 ϕ 1 + p 2 θ 2 cot 2 ϕ 2 + 2(p ϕ 1 p θ 2 cot ϕ 2 -p ϕ 2 p θ 1 cot ϕ 1 ) sin(θ 1 -θ 2 ).
One deduces the following:

Proposition 5 θ ′ 1 is a cyclic variable and p θ ′ 1 = p θ 1 + p θ 2 is a first integral.
Computed explicitely, this gives:

p θ 1 + p θ 2 = (-p x 1 y 1 + p y 1 x 1 ) + (-p x 2 y 2 + p y 2 x 2 ) = -p, Rq ,
where R is the matrix introduced in section 1.6. Using Lie brackets, this has the following important interpretation:

p, Rq = p, [F 1 , F 2 ](q) = {H 1 , H 2 }(q, p).

Lie bracket computations

To compute Lie brackets, one observes that the system can be lifted onto the Lie product G × G, where G is the semidirect Lie product GL(3, R) × s 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). One uses the notation F 0 = (A 0 , a 0 ) with A 0 = diag (-Γ 1 , -Γ 1 , -γ 1 , -Γ 2 , -Γ 2 , -γ 2 )
and a 0 = (0, 0, γ 1 , 0, 0, γ 2 ), whereas the control fields

F 1 , F 2 are identified to B 1 = diag (C 1 ,C 1 ) and B 2 = 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 i j = (δ i j ) (Kronecker symbol). Note that C 1 and C 2 generates the Lie algebra so(3) since [C 1 ,C 2 ] = C 3 = E 21 -E 12 .
According to the Lie bracket computation rule on the semidirect Lie product, we can reduce the computations to matrix Lie brackets. We have the following results, see [START_REF] Bonnard | algebraic geometric classification of the singular flow in the contrast imaging problem in nuclear magnetic resonance[END_REF] for the details. Consider first the case of a single spin and let A = diag (-Γ , -Γ , -γ) and denote δ = γ -Γ . We have:

Lemma 2 There are two cases i) If δ = 0, the Lie algebra generated by {A,C 1 ,C 2 } is gl(3, R). ii) If δ = 0, the Lie algebra generated by {A,C 1 ,C 2 } is RA ⊕ so(3).
The computations for the coupled system give the following result.

Lemma 3 Denoting δ i = γ i -Γ i ,
0 , B 1 , B 1 } is RA 0 ⊕ sl(3, R) ⊕ sl(3, R).
Next, we present in details the Lie brackets needed in our computations. Each entry is formed by a couple of 3 dimensional vectors and we use the notation omiting the indices.

-Length 1: [START_REF] Betts | Practical methods for optimal control using nonlinear programming[END_REF] Stratification of the surface Σ : H 1 = H 2 = 0 and partial classification of the extremal flow near Σ Preliminary results about this classification are coming from [7,16]. We make a self contained presentation building on [7, 16] with additional results.

F 0 = (-Γ x, -Γ y, γ(1 -z)) F 1 = (0, -z, y) F 2 = (z, 0, -x) -Length 2: [F 0 , F 1 ] = (0, γ -δ z, -δ y) [F 0 , F 2 ] = (-γ + δ z, 0, δ x) [F 1 , F 2 ] = (-y, x, 0) -Length 3: [[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δ y, -γ + 2δ z) [[F 0 , F 1 ], F 2 ] = (δ y, δ x, 0) = [[F 0 , F 2 ], F 1 ] [[F 0 , F 2 ], F 2 ] = (-2δ x, 0, 2δ z -γ) [[F 0 , F 1 ], F 0 ] = (0, -γ(γ -2Γ ) + δ 2 z, -δ 2 y) [[F 0 , F 2 ], F 0 ] = (γ(γ -2Γ ) -δ 2 z, 0, δ 2 x) 1.7.
Let z = (q, p) be a curve solution of #-

H 0 + u 1 #- H 1 + u 2 #- H 2 .
Differentiating H 1 and H 2 along such a solution, one gets:

Ḣ1 = {H 0 , H 1 } -u 2 {H 1 , H 2 } Ḣ2 = {H 0 , H 2 } + u 1 {H 1 , H 2 } (8)
and extremals of order zero are such that the control (u 1 , u 2 ) is defined by

u i = H i √ H 2 1 +H 2 2 , i = 1, 2.
Let z 0 ∈ Σ 2 = Σ \ {H 1 , H 2 } = 0 and define the control u s by:

u s (z) = (-{H 0 , H 2 }(z), {H 0 , H 1 }(z)) {H 1 , H 2 }(z) . (9)
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).
The behaviours of the extremals of order zero near Σ can be analyzed using a nilpotent model where all Lie brackets at z 0 ∈ Σ 2 of length ≥ 3 are zero. Denoting:

{H 0 , H 1 }(z 0 ) = a 1 , {H 0 , H 2 }(z 0 ) = a 2 , {H 1 , H 2 }(z 0 ) = b
and using polar coordinates

H 1 = r cos θ , H 2 = r sin θ , we deduce from (8): ṙ = a 1 cos θ + a 2 sin θ θ = 1 r (b -a 1 sin θ + a 2 cos θ ). (10)
To analyze this equation, we write:

a 1 sin θ -a 2 cos θ = A sin(θ + ϕ) with A tan ϕ = -a 2 /a 1 , A = a 2 1 + a 2 2 .
Hence the equation θ = 0 leads to the relation A sin(θ + ϕ) = b, which has two distinct solutions on [0, 2π[ denoted θ 0 , θ 1 if and only if A > |b|, one solution if A = |b| and zero solution if A < |b|. Moreover θ 1 -θ 0 = π if and only if b = 0. Plugging θ 0 , θ 1 in (10), we have:

Lemma 4 If a 2 1 + a 2 2 >
|b| and b = 0, we have a broken extremal formed by concatenating two extremals of order zero at each point z 0 of Σ 2 .

At such a point z 0 of Σ 2 , the singular control given by ( 9) is such that

u 2 s,1 + u 2 s,2 = a 2 1 +a 2 2
b 2 > 1 and hence is not admissible. This describes the behaviours of extremals of order zero near a generic point of z 0 . Next we analyze more degenerated situations. One needs the following concepts.

Goh condition. Higher order necessary optimality conditions along singular extremals are related to finitness of the index of the quadratic form associated to the intrinsic second order derivative, see [START_REF] Bonnard | Singular trajectories and their role in control theory[END_REF]. First we have:

{H 1 , H 2 } = 0 (Goh condition). (11)
Using [START_REF] Bonnard | Singular trajectories and their role in control theory[END_REF], one gets the relations:

H 1 = H 2 = {H 1 , H 2 } = 0 and
{H 1 , H 2 } = {H 0 , H 1 } = {H 0 , H 2 } = 0. ( 12
)
Then differentiating, we obtain:

{{H 1 , H 2 }, H 0 } + u 1 {{H 1 , H 2 }, H 1 } + u 2 {{H 1 , H 2 }, H 2 } = 0 (13) {{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 ( 14
)
This leads in general to three relations to compute two control components, see [START_REF] Chitour | Genericity results for singular curves[END_REF] for a discussion of the non genericity of such relations. However, according to Lie brackets computations, we have in our case:

Lemma 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 ( 13) is satisfied.

To analyze [START_REF] Bulirsch | Introduction to numerical analysis[END_REF], which can be written à + Bu = 0, we use:

H 1 = H 2 = {H 1 , H 2 } = {H 0 , H 1 } = {H 0 , H 2 } = 0. Hence p is orthogonal to F 1 , F 2 , [F 1 , F 2 ], [F 0 , F 1 ], [F 0 , F 2 ]. Introducing: A = A 1 A 2 , B = B 1 B 3 B 2 B 4 , 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 ]), A 2 = det(C, [[F 0 , F 2 ], F 0 ]), and
B 1 = det(C, [[F 0 , F 1 ], F 1 ]), B 2 = det(C, [[F 0 , F 2 ], F 1 ]), B 3 = det(C, [[F 0 , F 1 ], F 2 ]), B 4 = det(C, [[F 0 , F 2 ], F 2 ]
), the relation ( 14) leads to:

A + Bu = 0, and if det B = 0, one gets the singular control given by the feedback:

u ′ s = -B -1 A (15)
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 | ≤ 1. One denotes: Σ ′ 2 : H 1 = H 2 = {H 1 , H 2 } = {H 0 , H 1 } = {H 0 , H 2 } = 0 \ det B = 0 and singular extremals solutions of #- H ′ s = #- H 0 + u ′ s,1 #- H 1 + u ′ s,2 #- H 2 project onto solutions of Q ′ s .
Computations of det B. We give det B in the general case:

det B =(x 1 y 2 -x 2 y 1 ) 4 (δ 1 -δ 2 )(2δ 1 z 1 -γ 1 )(2δ 2 z 2 -γ 2 ) 2(δ 2 1 γ 2 z 1 -δ 2 2 γ 1 z 2 ) -γ 1 γ 2 (δ 1 -δ 2 ) -2δ 1 δ 2 (γ 1 z 2 -γ 2 z 1 ) , (16)
in the fluid case (P 1 ):

det B = -(180625(x 1 y 2 -x 2 y 1 ) 4 (34z 1 + 3) 2 )/12573377307869184
and in the blood case (P 2 ):

det B = -244140625(x 1 y 2 -x 2 y 1 ) 4 (23z 2 + 2)(52z 1 + 1) (13208z 1 -2921z 2 + 162) /2900177174423310336.
Additional singular extremals can be contained in the surface:

Σ ′′ 2 : H 1 = H 2 = {H 1 , H 2 } = {H 0 , H 1 } = {H 0 , H 2 } = det B = 0.
In particular there are those associated to extremals such that u 2 = 0. Indeed, this corresponds to replace H 2 by εH 2 and to impose ε = 0. The remaining relations are then:

H 1 = {H 0 , H 1 } = 0
and from ( 14) one gets the relations:

{{H 0 , H 1 }, H 0 } + u 1 {{H 0 , H 1 }, H 1 } = 0. ( 17
)
Therefore this defines the singular control:

u ′′ 1,s = - {{H 0 , H 1 }, H 0 } {{H 0 , H 1 }, H 1 } (18)
and the associated Hamiltonian

H ′′ 1,s = H 0 + u ′′ 1,s H 1 , with u ′′ 1,s = u 1,s
, where u 1,s is defined section 1.3.2. In particular we have the following result: Proposition 6 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.
The classification of regular extremals near the switching surface in the singleinput case is coming from [24], see also [START_REF] Bonnard | Singular trajectories and their role in control theory[END_REF].

To resume our analysis, we have:

Theorem 1
The singular extremals in the contrast problem are classified by the following stratification of vector fields:

1. In Σ 2 : H 1 = H 2 = 0 \ {H 1 , H 2 } = 0
, the singular control is given by (9):

u s (z) = (-{H 0 , H 2 }(z), {H 0 , H 1 }(z)) {H 1 , H 2 }(z) . 2. In Σ ′ 2 : H 1 = H 2 = {H 1 , H 2 } = {H 0 , H 1 } = {H 0 , H 2 } = 0 \ det B = 0
, the singular control is given by the feedback (15):

u ′ s = -B -1 A. 3. In Σ ′′ 2 : H 1 = H 2 = {H 1 , H 2 } = {H 0 , H 1 } = {H 0 , H 2 } = det B = 0,
there exists singular controls such that u 2 = 0 and u 1 is given by (18):

u 1,s = u ′′ 1,s = - {{H 0 , H 1 }, H 0 } {{H 0 , H 1 }, H 1 } .
Remark 1 Observe that additional singular extremals can exist in Σ ′′ 2 . The respective codimension are: two for Σ 2 , five for Σ ′ 2 and six for Σ ′′ 2 .

The numerical methods

In this section we present the three numerical methods used for the resolution of the contrast problem.

Bocop

The so-called direct approach transforms the infinite dimensional optimal control problem (OCP) into a finite dimensional optimization problem (NLP). This is done by a discretization in time applied to the state and control variables, as well as the dynamics equation. These methods are usually less precise than indirect methods based on Pontryagin's Maximum Principle, but more robust with respect to the initialization. Also, they are more straightforward to apply, hence their wide use in industrial applications. We refer the reader to for instance [START_REF] Betts | Practical methods for optimal control using nonlinear programming[END_REF] and [START_REF] Nocedal | Numerical optimization[END_REF] for more details on direct transcription methods and NLP algorithms.

Summary of the time discretization:

t ∈ [0,t f ] → {t 0 = 0, . . . ,t N = t f } z(•), u(•) → X = {z 0 , . . . , z N , u 0 , . . . , u N-1 ,t f } Criterion → min c(z N ) Dynamics → (ex : Euler) z i+i = z i + h f (z i , u i ) Admissible Controls → u i ∈ U Boundary Conditions → Φ(z 0 , z N ) = 0
We therefore obtain a nonlinear programming problem on the discretized state and control variables

(NLP) min F(X) = c(z N ) LB ≤ C(X) ≤ UB
All tests were run using the Bocop software [5]. The discretized nonlinear optimization problem is solved by the well-known IPOPT solver [START_REF] Wächter | On the implementation of a primal-dual interior point filter line search algorithm for large-scale nonlinear programming[END_REF] with MUMPS [START_REF] Amestoy | A fully asynchronous multifrontal solver using distributed dynamic scheduling[END_REF], while the derivatives are computed by sparse automatic differentiation with ADOL-C [38] and COLPACK [START_REF] Gebremedhin | Exploiting sparsity in jacobian computation via coloring and automatic differentiation: a case study in a simulated moving bed process[END_REF].

In the numerical experiments, we used a midpoint (implicit 2nd order) discretization with 1000 to 5000 time steps. Execution times on a Xeon 3.2GHz CPU were a few minutes. 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.

HamPath in the single-input case

In this subsection, we recall the method used to solve the contrast problem in the single-input case, based on the HamPath code, see [START_REF] Bonnard | Geometric numerical methods and results in the control imaging problem in nuclear magnetic resonance[END_REF]18]. The HamPath software is based upon indirect methods: simple and multiple shooting, differential continuation (or homotopy) methods and computation of the solutions of the variational equations, needed to check second order conditions of local optimality.

Multiple shooting method for the contrast problem in the single-input case

We introduce the λ variable such that t f = λ T min . For a fixed λ ≥ 1, we must solve the contrast problem by multiple shooting method since the optimal control is discontinuous (see 1.3.2). Shooting methods consist in finding a zero of a specific function and use Newton like algorithms. While simple shooting leads to solve a two points boundary value problem, multiple shooting takes into account intermediate conditions. This leads to the first numerical difficulty.

Difficulty 1 Before calling the multiple shooting method, we need to know the structure of the optimal solution and we need a good approximation of the solution to make the method converge.

Fix t f = λ T min with λ > 1, assume to simplify that the optimal solution z(t) = (q(t), p(t)), t ∈ [0,t f ], is made of a single BS-sequence and denote t 1 , 0 < t 1 < t f , the switching time between the regular arc and the singular arc. The generalization to the case nBS is straightforward. By construction, z(•) is a concatenation of the integral curve of #-

H joining z(0) to z 1 = z(t 1 ) with u 1 (t) = sgn (H 1 (z(t))), t ∈ [0,t 1 )
, and of the extremal (of order two) solution of #-H s starting at z 1 , where z 1 must satisfy (see subsection 1.3.2): [START_REF] Gebremedhin | Exploiting sparsity in jacobian computation via coloring and automatic differentiation: a case study in a simulated moving bed process[END_REF] and then

H 1 (z 1 ) = {H 1 , H 0 }(z 1 ) = {{H 1 , H 0 }, H 0 }(z 1 ) + u 1 {{H 1 , H 0 }, H 1 }(z 1 ) = 0,
z 1 ∈ Σ ′ 1 .
Remark 2 We use the same notation z 1 for the third component of the first spin state q 1 = (x 1 , y 1 , z 1 ) and for the first switching point between bang and singular arcs. In the first case, z 1 is scalar and in the second case, z 1 = z(t 1 ) = (q(t 1 ), p(t 1 )) lives in the cotangent space T * Q, thus no confusion is possible.

Let fix the notations for both bang and singular extremals. We write z(•,t 0 , z 0 ) or z(•,t 0 , q 0 , p 0 ) the solution of ż(•) = #-H(z(•)), z(t 0 ) = z 0 = (q 0 , p 0 ), with q 0 = (q 1,0 , q 2,0 ), p 0 = (p 1,0 , p 2,0 ) and t 0 = 0. Under these notations each intial spin state is q i,0 = (y i,0 , z i,0 ) = (0, 1), i = 1, 2, and the adjoint vector is given by p i,0 = (p y i,0 , p z i,0 ).

The singular extremal is denoted z s (•,t 1 , z 1 ) solution of żs (•) = #- H s (z s (•)), z s (t 1 ) = z 1 ,
where z s = (q s , p s ) with q s = (q 1,s , q 2,s ) and p s = (p 1,s , p 2,s ). Under these notations the boundary conditions from subsection 1.2.2 become:

q 1,s (t f ,t 1 , z 1 ) = (0, 0) q 2,s (t f ,t 1 , z 1 ) = p 2,s (t f ,t 1 , z 1 ). (20)
Hence, the bang and singular arcs z(•,t 0 , z 0 ), z 0 = (q 0 , p 0 ) and z s (•,t 1 , z 1 ) are entirely defined by (t 0 , z 0 ,t 1 , z 1 ,t f ).

Remark 3 In the contrast problem, t 0 , t f = λ T min and q 0 are fixed and known. For z 1 , we can either consider it as a dependant variable if we replace z 1 by the value z(t 1 ,t 0 , z 0 ), or as an independant variable and add the matching conditions:

z(t 1 ,t 0 , z 0 ) = z 1 , (21)
expressing the continuity of the extremal at the switching time.

We define now the multiple shooting function Sλ (p 0 ,t 1 , z 1 ) which maps (p 0 ,t 1 , z 1 ) to the conditions ( 19), ( 20) and [START_REF] Henrion | Optimal switching control design for polynomial systems: an LMI approach[END_REF]. See [START_REF] Bulirsch | Introduction to numerical analysis[END_REF] for details about multiple shooting methods with application in optimal control and [31] for applications in the Bang-Singular-Bang case. The function Sλ is given by:

Sλ : R 3n+1 → R 3n+2   p 0 t 1 z 1   → Sλ (p 0 ,t 1 , z 1 ) =       H 1 (z 1 ) {H 1 , H 0 }(z 1 ) q 1,s (t f ,t 1 , z 1 ) q 2,s (t f ,t 1 , z 1 ) -p 2,s (t f ,t 1 , z 1 ) z(t 1 ,t 0 , z 0 ) -z 1       , ( 22
)
Remark 4 The shooting system [START_REF] Henrion | GloptiPoly 3: Moments, Optimization and Semidefinite Programming[END_REF] has an extra condition which is redundant and can be eliminated. See [3] for the general case where the shooting system has more equations than unknowns and where the Gauss-Newton method is used to compute a zero of the shooting function.

Lemma 6 In the shooting system [START_REF] Henrion | GloptiPoly 3: Moments, Optimization and Semidefinite Programming[END_REF], the transversality condition:

q 2,s (t f ,t 1 , z 1 ) -p 2,s (t f ,t 1 , z 1 ) = 0 R 2 reduces to p 2,s (t f ,t 1 , z 1 ), q 2,s (t f ,t 1 , z 1 ) -|q 2,s (t f ,t 1 , z 1 )| 2 = 0 ( 23
)
and the system is square.

Proof Denote the singular extremal z s = (q s , p s ), with q s = (q 1,s , q 2,s ), p s = (p 1,s , p 2,s ) and q i,s = (y i,s , z i,s ), p i,s = (p y i,s , p z i,s ), i = 1, 2. The condition H 1 (z 1 ) = 0 implies that

H 1 (z s [t f ]) = 0 where [t f ] stands for (t f ,t 1 , z 1 ), since z s (•) is contained in Σ ′ 1 . Besides, the terminal condition q 1,s [t f ] = 0 R 2 implies that H 1 (z s [t f ]) = -p y 2,s [t f ]z 2,s [t f ] + p z 2,s [t f ]y 2,s [t f ] = det(q 2,s [t f ], p 2,s [t f ]),
i.e. q 2,s and p 2,s are colinear at t f . Hence we have:

     H 1 (z 1 ) = 0 q 1,s [t f ] = 0 R 2 q 2,s [t f ] -p 2,s [t f ] = 0 R 2 ⇔      H 1 (z 1 ) = 0 q 1,s [t f ] = 0 R 2 p 2,s [t f ], q 2,s [t f ] -|q 2,s [t f ]| 2 = 0 ⊓ ⊔
The square system replaces the system [START_REF] Henrion | GloptiPoly 3: Moments, Optimization and Semidefinite Programming[END_REF] and we define the multiple shooting function S λ by:

S λ : R 3n+1 → R 3n+1   p 0 t 1 z 1   → S λ (p 0 ,t 1 , z 1 ) =       H 1 (z 1 ) {H 1 , H 0 }(z 1 ) q 1,s [t f ] p 2,s [t f ], q 2,s [t f ] -|q 2,s [t f ]| 2 z(t 1 ,t 0 , z 0 ) -z 1       . ( 24
)
The multiple shooting method for the contrast problem consists in finding a zero of the multiple shooting function S λ , i.e. in solving S λ (p 0 ,t 1 , z 1 ) = 0.

Algorithmic of the HamPath code The fortran hybrid Newton method hybrj (from the minpack library [32]) is used to solve the nonlinear system S λ (p 0 ,t 1 , z 1 ) = 0. The hybrj code implements the Powell hybrid algorithm [START_REF] Powell | A hybrid method for nonlinear equations[END_REF] based upon Newton-Raphson method combined with Levenberg-Marquardt algorithm. The Jacobian of S λ must be provided to hybrj and must be invertible to get convergence. In [3], the authors analyze the convergence of the shooting algorithm and the invertibility of the shooting function for optimal control problems with singular arcs. Define the function

E(t 0 ,t 1 ,t f , z 0 , z 1 , ez 0 , ez 1 ) : R 3+4×2n → R 3n+1 by S λ (p 0 ,t 1 , z 1 ) = E(t 0 ,t 1 ,t f , (q 0 , p 0 ), z 1 , z(t 1 ,t 0 , z 0 ), z s (t f ,t 1 , z 1 )),
where

z 0 = (q 0 , p 0 ), ez 0 = z(t 1 ,t 0 , z 0 ) and ez 1 = z s (t f ,t 1 , z 1 ).
Remark 5 The user must only provide the function E (S λ is not needed) and the two ordered true-Hamiltonians H and H s (the regular and singular controls must be computed), i.e. H for the first arc and H s for the second arc. Then, the code assembles automatically the multiple shooting function and its Jacobian which are given to hybrj.

We detail how S ′ λ (p 0 ,t 1 , z 1 ) is computed. The Jacobian is given by:

∂ S λ ∂ p 0 = ∂ E ∂ z 0 + ∂ E ∂ ez 0 • ∂ z ∂ z 0 (t 1 ,t 0 , z 0 ) 0 I n ∂ S λ ∂t 1 = ∂ E ∂t 1 + ∂ E ∂ ez 0 • ∂ z ∂t f (t 1 ,t 0 , z 0 ) + ∂ E ∂ ez 1 • ∂ z s ∂t 0 (t f ,t 1 , z 1 ) ∂ S λ ∂ z 1 = ∂ E ∂ z 1 + ∂ E ∂ ez 1 • ∂ z s ∂ z 0 (t f ,t 1 , z 1 ), (25)
where

∂ ∂t 0 , ∂ ∂t f and ∂ ∂ z 0
mean respectively the partial derivative with respect to the initial time, the final time and the initial condition. All the partial derivatives of E are computed by Automatic Differentiation (AD) with the tapenade software [START_REF] Hascoët | The Tapenade Automatic Differentiation tool: principles, model, and specification[END_REF]. The derivative ∂ z 0 z(t 1 ,t 0 , z 0 ) is the solution at t 1 of the variational equation:

δ ż(t) = d #- H(z(t))δ z(t), δ z(t 0 ) = I 2n , z(t 0 ) = z 0 , whereas ∂ z 0 z s (t f ,t 1 , z 1 )
is the solution at t f of the variational equation:

δ ż(t) = d #- H s (z s (t))δ z(t), δ z(t 1 ) = I 2n , z s (t 1 ) = z 1 . Then ∂ t f z(t 1 ,t 0 , z 0 ) is simply equal to #- H(z(t 1 ,t 0 , z 0 )) while ∂ t 0 z s (t f ,t 1 , z 1 ) is the solu- tion at t f of: δ ż(t) = d #- H s (z s (t))δ z(t), δ z(t 1 ) = - #- H s (z 1 ), z s (t 1 ) = z 1 .
HamPath uses again AD to compute

#- H, d #- H and d #- H s .

Homotopy method on the final time

Continuation techniques or homotopic methods are well known and widely used: see for example [START_REF] Allgower | Introduction to numerical continuation methods[END_REF] for theoretical and numerical details. We just point out some facets encountered in optimal control, especially for the contrast problem.

In the contrast problem, natural parameters are the relaxation times and the final time t f . In this article, our goal is to give the optimal policy with respect to t f . Thus we define the homotopic function h : R N × R → R by h(x, λ ) = S λ (x), λ ≥ 1, where S λ is the multiple shooting function (24) and x = (p 0 ,t 1 , z 1 ) in the BS case.

The classical difficulties about homotopic methods consist in assuring that a curve in h -1 (0) exists, is sufficiently smooth and will intersect a fixed target homotopic level in a finite length. Suppose h is continuously differentiable and that we know x 0 such that h(x 0 , λ 0 ) = 0 and

rank ∂ h ∂ x (x 0 , λ 0 ) = N.
Suppose also that 0 is a regular value of h. Then a continuously differentiable curve starting from (x 0 , λ 0 ) exists and is either diffeomorphic to a circle or the real line. The curves in h -1 (0) are disjoints, and we call each branch of h -1 (0) a path of zeros.

Remark 6 For a 2BS structure, we must introduce two others arcs and the variable x is x = (p 0 ,t 1 , z 1 ,t 2 ,t 2 , z 3 ). Hence, in the BS and 2BS cases, the paths of zeros we compare are not defined in the same spaces.

Besides, it is possible that a 2BS solution is continuously deformed in a way that a BS or a 3BS structure appeared. The continuation process would stop at the change of structure. A BS-sequence must be removed if a bang arc tends to disappear, or added if a singular arc tends to saturate the constraint |u| = 1. Finally, the following difficulty is crucial.

Difficulty 2 An important issue in the contrast problem is to deal with the many local solutions.

For a given value λ , we must compare the cost associated to each component of h -1 (0) ∩ {λ = λ }. This global aspect is responsible for a possible loss of regularity on the value function λ → c(q(t f )), t f = λ T min and on the optimal path of zeros.

Algorithmic of the HamPath code Unlike well-known prediction-correction methods, see [START_REF] Allgower | Introduction to numerical continuation methods[END_REF], we merely follow the path of zeros by integrating the differential system with a high order Runge-Kutta scheme, without any correction. The Jacobian of the homotopic function is computed using variational equations as in subsection 2.2.1. See [START_REF] Bonnard | Geometric numerical methods and results in the control imaging problem in nuclear magnetic resonance[END_REF][START_REF] Caillau | Differential continuation for regular optimal control problems[END_REF] for details about the algorithmic.

lmi

The direct and indirect methods provide local optimal solutions. By comparing the different paths of zeros, one obtains a strong candidate solution whose global optimality must be asserted. This can be done by the moment/lmi techniques described below.

The moment approach is a global optimization technique that transforms a nonlinear, possibly infinite-dimensional optimization problem into convex and finitedimensional relaxations in the form of Linear Matrix Inequalities (lmi). We follow [START_REF] Lasserre | Nonlinear optimal control via occupation measures and LMI-relaxations[END_REF] for the specific case of optimal control with bounded controls and [START_REF] Lasserre | Positive polynomials and their applications[END_REF] for the main steps of the method.

An overview of the method is given in Fig. 1. The first step consists in linearizing the problem by casting it as a Linear Program (LP) on a measure space, a problem often referred to as a generalized moment problem. This can be performed by the use of so-called occupation measures, encoding admissible trajectories. The second step is to exploit the problem's structure, here given by its polynomial data, to manipulate the measures by their moment sequences. This leads to a Semi-Definite Program (SDP) with countably many decision variables, one for each moment. The third and last step is to truncate this last problem to a finite set of those moments, leading to a relaxation in the form of Linear Matrix Inequalities (lmi). Those relaxations can be solved by off-the-shelf software.

In the sequel, we show J ≥ J LP ≥ J SDP ≥ J r LMI , which is the relation of interest if the method is used to yield upper bounds on the contrast. We also point to

J = inf . . .

Optimisation problem

J LP = inf µ . . .

Linear Program on measures

J SDP = inf y . . .

Semi-Definite Program on moments

J r LMI = inf (y k ) k≤2r . . .

lmi relaxation of order r

Step 1

Step 2

Step 3

Fig. 1 Overview of the moment approach.

Step 1: the initial optimization problem is cast as a Linear Program on a measure space.

Step 2: the latter is itself recast as a Semi-Definite Program on moment sequences.

Step 3: truncating the sequences to moments of degree smaller than 2r leads to the Linear Matrix Inequality relaxation of order r.

relevant theorems in the literature showing that for the problem at hand, the cost of the lmi relaxations asymptotically converge to that of the original problem, i.e. lim r→∞ J r LMI → J.

Step 1: formulation as a measure LP

The first step is to embed problem (OCP) into a linear program on measures by the use of occupation measures. We first set up the notations. For z ∈ Z ⊂ R n , R[z] denotes the ring of polynomials in the variables z. M + (Z) is the set of finite, positive Borel measures supported on compact set Z and Z v(z) dµ(z) denotes the integration of a continuous function v ∈ C(Z) with respect to µ ∈ M + (Z). We also write simply v dµ should no ambiguity arise. For concise notations, define compact set K = [0,t f ] × Q × U. Consider then the following linear program on measures, an instance of the generalized moment problem:

J LP = inf µ,ϕ -|q 2 | 2 dϕ s.t. v(t f , •) dϕ -v(0, 0) = ∂ v ∂t + ∂ v ∂ q F(q, u) dµ, ∀v ∈ R[t, q], µ ∈ M + (K), ϕ ∈ M + (Q f ).
(LP)

Proposition 7 (Measure embedding) To every pair (u(t), q(t)) admissible for (OCP) corresponds a pair (µ, ϕ) admissible for (LP) achieving the same cost, hence J ≥ J LP .

Proof For each admissible pair (u(t), q(t)), define its corresponding time occupation measure µ ∈ M + (K) by:

µ(A, B, C) := [0,t f ]∩A δ q(t) (B) δ u(t) (C) dt.
Here, δ x is the Dirac measure supported at x, and A, B and C are Borel subsets of resp. [0,t f ], Q and U. That is, µ measures the time "spent" by the admissible triplet (t, u(t), x(t)) on Borel subsets of K. Similarly, we define the final state occupation measure ϕ ∈ M + (Q f ) for the same admissible pair simply as ϕ := δ q(t f ) . Then, evaluating each polynomial test function v along the admissible trajectory yields, by the chain rule:

v(t f , q(t f )) -v(0, q(0)) = t f 0 dv(t, q(t)) = t f 0 ∂ v ∂t (t, q(t)) + ∂ v ∂ q (t, q(t)) F(q(t), u(t)) dt.
For the constructed occupation measures, this last temporal integration becomes the weak dynamics of (LP). The proposition follows.

⊓ ⊔

Note that the set of admissible vector fields for dynamics F is convex for any q ∈ Q, such that theorem 3.6 (ii) of [START_REF] Lasserre | Nonlinear optimal control via occupation measures and LMI-relaxations[END_REF] holds and J = J LP .

In [START_REF] Henrion | Optimal switching control design for polynomial systems: an LMI approach[END_REF], an alternative formulation as a measure LP is given for switched systems. For these systems, controls are modeled by measures, instead of the above construction defining occupation measures supported on the control space. The alternative formulation is believed to be more computationally effective for the moment approach, as it involves measures defined on Euclidean spaces of lower dimension.

Observe that in the single input case, problem (OCP) is bilinear, with the control taking its values inside the polytope defined by U = conv{-1, +1}. The control problem can then be cast as a switched system, the first mode being driving the system with u(t) = +1, and the second one by u(t) = -1. Following [START_REF] Henrion | Optimal switching control design for polynomial systems: an LMI approach[END_REF], this leads to:

J LP ′ = inf µ 1 ,µ 2 ,ϕ -|q 2 | 2 dϕ s.t. ∀v ∈ R[t, q] : v(t f , •) dµ f -v(0, 0) = ∂ v ∂t + ∂ v ∂ q (F 0 + F 1 ) dµ 1 + ∂ v ∂t + ∂ v ∂ q (F 0 -F 1 ) dµ 2 , µ 1 , µ 2 ∈ M + ([0,t f ] × Q), ϕ ∈ M + (Q f ).
(LP')

In the following, both measure formulations will be tested on the contrast problem.

Step 2: recasting as a moment SDP

Unfortunately, there is no generic tractable method to solve problems (LP) or (LP'), and additional structure on problem data is required. For optimal control problem (OCP), this structure is provided by the polynomial cost and dynamics, as well as the basic semi-algebraic characterization of the compact sets K and Q f . It is then possible to manipulate measures by their moments in a given polynomial basis, which yields a semi-definite program on countably many moments, with cost J SDP = J LP .

In the following, we concentrate on (LP) in the single-input case, as the discussion is straightforward to adapt for the bi-input case or for alternative formulation (LP').

We first set up the notations. Define the moment of order k ∈ N n of a measure µ supported on Z ⊂ R n by using the multi-index notation

y µ k = Z z k 1 1 • • • z k n n dµ(z) := Z z k dµ(z). ( 26
)
Conversely, µ is said to be a representing measure for an arbitrary sequence of reals y µ k if µ satisfies ( 26) for all k ∈ N n . The degree of monomial z k is denoted by deg

z k = |k| 1 = ∑ n i=1 k i . N n d denotes the set {k ∈ N n : |k| 1 ≤ n}. Then, with a multi-indexed sequence of reals y = (y k ) with k ∈ N n , let L y : R[z] → R be the (Riesz) linear functional of f = ∑ k f k z k defined by f → L y ( f ) = ∑ k f k y k , f ∈ R[z].
Note that for ease of presentation, we represent the array of coefficients of polynomial f (z) by the same letter f . When confusion may arise, we write the polynomial explicitly by f (z). Also observe that when the sequence y is generated by a measure such as in [START_REF] Lapert | Exploring the physical limits of saturation contrast in Magnetic Resonance Imaging[END_REF], L y ( f ) = f dµ by linearity. Define the moment matrix of order d ∈ N associated with an arbitrary sequence of reals y as the real symmetric matrix M d (y) whose (multi-indexed) (i, j)th entry reads

M d (y)[i, j] = L y z i+ j = y i+ j , ∀i, j ∈ N n d . (27)
Similarly, define the localizing matrix of order d associated with y and g ∈ R[z] as the real symmetric matrix M d (g * y) whose (i, j)th entry reads

M d (g(z) * y)[i, j] = L y g(z) z i+ j = ∑ k g k y i+ j+k , ∀i, j ∈ N n d . (28)
These constructions are central for the following well-known fact:

Proposition 8 Let Z be a basic, semi-algebraic set defined by Z := {z ⊂ R n : g i (z) ≥ 0, i = 1 . . . n Z } . Then, a necessary condition for a sequence of reals y µ to have a representing measure µ ∈ M + (Z) is:

M d (y µ ) 0, M d (g i * y µ ) 0, ∀d ∈ N.
Proof Take an arbitrary polynomial p(z) of degree d. Then for a positive measure µ supported on Z, we have necessarily 0 ≤ p 2 (z) dµ = p ′ M d (y)p and also 0 ≤ p 2 (z) g i (z) dµ = p ′ M d (g i (z) * y)p. Since this holds for any p(x), the result follows.

⊓ ⊔

Returning back to (LP), notice that K can be represented by the following basic, semi-algebraic representation

K =            t ∈ R, q ∈ R 4 , u ∈ R : g 1 (t, q, u) := t (t f -t) ≥ 0 g 2 (t, q, u) := 1 -q 2 11 -q 2 12 ≥ 0 g 3 (t, q, u) := 1 -q 2 21 -q 2 22 ≥ 0 g 4 (t, q, u) := 1 -u 2 ≥ 0            , (29)
and similarly, for Q f :

Q f =        q ∈ R 4 : g f 1 (q) := q 11 = 0 g f 2 (q) := q 12 = 0 g f 3 (q) := 1 -q 2 21 -q 2 22 ≥ 0        . ( 30
)
Then, consider the following problem posed on real sequences y µ and y ϕ :

J SDP = inf y µ , y ϕ L y ϕ (-|q 2 | 2 ) s.t. L y ϕ (v(t f , •)) -v(0, 0) = L y µ ∂ v ∂t + ∂ v ∂ q F(q, u) , ∀v ∈ R[t, q], M d (y µ ) 0, M d (g i * y µ ) 0, ∀d ∈ N, i = 1 . . . n K , M d (y ϕ ) 0, M d (g f i * y ϕ ) 0, ∀d ∈ N, i = 1 . . . n Q f .
(SDP) As a simple consequence of Prop. 8, we have Proposition 9 Problem (SDP) is a relaxation of (LP), i.e. J LP ≥ J SDP .

Note that sets K and Q f satisfy the assumptions of [START_REF] Lasserre | Positive polynomials and their applications[END_REF]Th. 3.8], such that J LP = J SDP .

Step 3: truncation of moment sequences

The last step is certainly the most obvious one. Problem (SDP) has countably many decision variables, and truncating those to a finite set leads directly to a convex and finite-dimensional program in the form of linear matrix inequalities.

Given that all polynomial data are quadratic, consider the following lmi problems for r ≥ 1:

J r LMI = inf (y µ k ) |k|≤2r , (y ϕ k ) |k|≤2r L y ϕ (-|q 2 | 2 ) s.t. ∀v ∈ R[t, q] : L y ϕ (v(t f , •)) -v(0, 0) = L y µ ∂ v ∂t + ∂ v ∂ q F(q, u) , M r (y µ ) 0, M r-s i (g i * y µ ) 0, i = 1 . . . n K , M r (y ϕ ) 0, M r-s f i (g f i * y ϕ ) 0, i = 1 . . . n Q f . (

Summary of the LMI method

The moment/LMI method approach for optimization consist in reformulating an optimization problem as a linear program on measures. When the data is polynomial, a hierarchy of LMI relaxations can be constructed, whose cost converge to that of the original problem.

The strong feature of the method is that those LMIs generate lower bounds on the true cost, and can therefore be used as certificates of global optimality. On the other hand, the weak points of the method are its poor algorithmic complexity for unstructured problem, as well as for the special case of optimal control, the unavailability of a generic method to recover controls.

Note that the passage to a given LMI relaxation starting from measure problem (LP) or (LP') can be fully automated with high-level commands using the GLOP-TIPOLY toolbox [START_REF] Henrion | GloptiPoly 3: Moments, Optimization and Semidefinite Programming[END_REF].

Numerical simulations, single-input case

We present here the results about the single-input case (see subsection 1.3.2) of the problem (OCP). From the experimental point of view we are interested in the following cases, the parameters being the relaxation times given in seconds and the minimal time to steer the first spin from the north pole to the center of the Bloch ball. Recalling difficulty 1, we need to know the control structure (meaning the number of Bang-Singular sequences) before calling the multiple shooting method. A possible strategy is to use continuation (or homotopy) techniques introducing for example Tychonov type regularization, in order to detect the nBS sequence and find a smooth approximation of the non-smooth solution. In [START_REF] Bonnard | Geometric numerical methods and results in the control imaging problem in nuclear magnetic resonance[END_REF]18] this kind of strategy is applied. Here, we use the Bocop software based upon direct methods (detailed in subsection 2.1) to obtain approximate optimal solutions used to initialize indirect shooting. In both cases P 1 and P 2 , the contrast problem has many local solutions, see difficulty 2, possibly with different control structures. Besides, the structure of the best policy can change depending on the final time, as detailed below in 3.5 and 3.6, see also [18] for full details.

In this section, we first detail in four steps how we can obtain a potential local optimal solution for a fixed t f chosen in the [T min , 3T min ] range. In subsection 3.1, optimization by the direct method provides first approximations, and in 3.2 the control structure is given to initialize indirect shooting in 3.3 and to check second order necessary conditions in 3.4. Then, we are interested in the influence of the final time on the best policy in subsection 3.5. We have to monitor crossing between value functions, change of structure along one same path and second order necessary conditions. The simulations from subsection 3.1 to subsection 3.5 are detailed only in the fluid case. In subsection 3.6, we give for the fluid and blood cases, the best sub-optimal syntheses we have and finally we compare these sub-optimal syntheses with estimates of global optima obtained from lmi techniques in 3.7. These comparisons are used to assert the global optimality of the solutions obtained from direct and indirect approaches.

Before presenting the results, we explain how we introduce the λ variable. To do that, we change the parameterization of our Mayer problem introducing the time τ such that λ T min τ = t. Using the notations from the introduction, the new problem becomes: min u(•) c(q(τ f )), q = λ T min F(q, u), q(τ 0 ) = q 0 , f (q(τ f )) = 0, τ 0 = 0 and τ f = 1. Then in the BS case, h(x, λ ) = S λ (x), x = (p 0 , τ 1 , z 1 ) where

S λ : R 3n+1 → R 3n+1   p 0 τ 1 z 1   → S λ (p 0 , τ 1 , z 1 ) =       H 1 (z 1 ) {H 1 , H 0 }(z 1 ) q 1,s ([τ f ]) p 2,s ([τ f ]), q 2,s ([τ f ]) -|q 2,s ([τ f ])| 2 z(τ 1 ,t 0 , z 0 , λ ) -z 1       , (31)
where [τ f ] stands for (τ f , τ 1 , z 1 , λ ).

Step 1: direct method

We present here the results for the direct approach. The only a priori information is the value of the minimum time transfer T min from [18], used to set the final time t f in the [T min , 3T min ] range. We note t f = λ T min with λ in [START_REF] Allgower | Introduction to numerical continuation methods[END_REF]3]. The state variables are initialized as constant functions, with the values

y 1 (•) = 0, z 1 (•) = 0.5, y 2 (•) = 0, z 2 (•) = 1.
For the control variables we also use a constant function, with four differents initializations u 1 (•) ∈ {0.05, 0.1, 0.25, 0.5}. The discretization method used is implicit midpoint (2nd order) with 500 to 1000 time steps depending on the problem. In order to improve convergence, we add a small regularization term to the objective to be minimized, ε

t f 0 |u(t)| 2 dt, with ε = 10 -5 .
We repeat the four optimizations (with the different control initializations) for λ in {1.1, 1.5, 1.8, 2.0, 3.0}. Bocop converges for all initializations except when λ = 3, where only one optimization out of four converges. Table 1 groups together the results from all the optimizations. For each solution, we give the contrast and the corresponding rank between each different initializations for a same value of λ . We also indicate the norm of the approximate initial adjoint vector p(0), taken from the multipliers for the discretized dynamics in the (NLP) problem. For both λ = 1.8 and λ = 2.0, three different initializations lead to the same solution while for λ = 1.1 and λ = 1.5, the four initializations converge to four different solutions. Note that solutions 11b and 11c seem to be identical, and similarly solutions 20a and 20b do not differ a lot, as shown on Figs. 2 and3. This illustrate Difficulty 2 about the many local solutions, due in particular to symmetry reasons.

Table 1 Fluid case, step 1: Direct method. Results of optimizations from direct methods. The final time t f = λ T min is taken such as λ ∈ {1.1, 1.5, 1.8, 2.0, 3.0}. The control is initialized by a constant in {0.05, 0.1, 0.25, 0.5}. For each solution, the contrast and the norm of the adjoint vector (coming from the multipliers for the discretized dynamics in the (NLP) problem) are given. For a fixed λ we present the rank of the solutions based on the contrast level. Note that if two ranks are equal, the two optimizations have converged to the same solution. Finally a name is given for each different solution. 1 obtained with Bocop. 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. Note that the first singular arcs, for each control subgraph, do not end at the same times.

λ u 1 (•) init. Contrast Rank |p(0)| Name 1.
-1 0 1 -1 0 1 y 1 z 1 -1 0 1 -1 0 1 y 2 z 2 0 0.2 0.4 0.6 0.8 1 -1 0 1 τ u 1 -1 0 1 -1 0 1 y 1 z 1 -1 0 1 -1 0 1 y 2 z 2 0 0.2 0.4 0.6 0.8 1 -1 0 1 τ u 1
Fig. 3 Fluid case, step 1: Direct method. Solutions 20a (left subplot) and 20b (right subplot) from Table 1 obtained with Bocop. 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 solution 20a, whereas for solution 20b, the two bang arcs are of opposite sign.

Step 2: determination of the structure

The structure is obtained by examination of the control plots. This gives the normalized switching times τ i and the state-costate

z i = z(τ i ) associated, i = 1, • • • , m -1, or similarly t i since t i = λ T min τ i ,
where m is the number of bang and singular arcs (a 2BS structure has m = 4 arcs). The variables τ i and z(τ i ), i = 1, • • • , m -1, with the initial adjoint vector p 0 = p(0) are required to initialize the multiple shooting method, see the definitions of the shooting function in the BS case at eqs. (24) (with time t) and (31) (with normalized time τ). Table 2 sums up the different structures for all the distinct solutions found in Table 1. We can already notice that the best structure found by the direct approach is 3BS for λ = 1.1 and then 2BS for greater values of λ .

Table 2 Fluid case, step 2: determination of the structure. Same as Table 1 but with only the distinct solutions and reorganized. For each solution, a name and the contrast are given. For a fixed λ we present the rank of the solutions based on the contrast level. The structure is denoted by the signs of the bang arcs, and we normalize the first bang arc to be positive.

λ

Step 3: multiple shooting method

Thanks to steps 1 and 2, we have good approximations of the solutions in order to initialize the multiple shooting method described at 2.2.1, which solve the following shooting equation (see eq. ( 31)):

S λ (p 0 , τ 1 , z 1 , • • • , τ m-1 , z m-1 ) = 0,
where m is the number of arcs, τ 0 = 0 and τ f = 1, and the structure is m 2 BS. However, since indirect shooting methods are very sensitive to the initial point, we perform for each of the 13 different solutions from Table 2 not one but many shoots. Note that one call to the method is about a few seconds. To do so we add a random small perturbation ε i ∈ [-0.1, 0.1] to τ i , i = 1, • • • , m -1, before each call to the method. The approximate new state-costate z ε i = z(τ ε i ), τ ε i = τ i + ε i is obtained thanks to the discretized solution from step 1. We finally solve the shooting equation

S λ (p 0 , τ ε 1 , z ε 1 , • • • , τ ε m-1 , z ε m-1 ) = 0 (32)
about a hundred times, with different random perturbations. Notice that with this algorithm, one initialization from Table 2 can provide several solutions. Table 3 groups together the distinct solutions obtained after convergence of the multiple shooting method. Four initializations (11b, 11c, 11d and 15b) from Table 2 never made the method converge, whereas for instance, the two initializations 20a and 20b gave 12 different solutions. Note that for a given λ , the best direct optimization from Table 2 always provides the best solution after shooting.

Remark 7 We keep every solution from Table 3 because of difficulty 2.

Table 3 Fluid case, step 3: multiple shooting. Results after many calls to the multiple shooting method. The shooting equation solved is (32). This table gives, the value of λ , the initial name, the contrast and rank from Table 2 obtained by direct optimizations. The new contrast and rank after shooting is also provided with the structure of the solutions and the norm of the shooting function after convergence. Note that for a given λ , the best direct optimization from Table 2 Fig. 4 groups together solutions 1/12 and 6/12 for λ = 2.0 from Table 3. The switching times (vertical dashed lines) from direct optimization 20a are added on the control subplot to compare the initialization 20a with the solutions obtained after convergence of the multiple shooting method.

Step 4: second order conditions

According to proposition 3.2 from [START_REF] Bonnard | Geometric numerical methods and results in the control imaging problem in nuclear magnetic resonance[END_REF], the non-existence of conjugate points on each singular arc of a candidate solution is a necessary condition of local optimality. See [START_REF] Bonnard | Geometric numerical methods and results in the control imaging problem in nuclear magnetic resonance[END_REF] for details about conjugate points in the contrast problem. Here, we compute for each singular arc of all the solutions from Table 3, the first conjugate point along

-1 0 1 -1 0 1 y 1 z 1 -1 0 1 -1 0 1 y 2 z 2 0 0.2 0.4 0.6 0.8 1 -1 0 1 τ u 1 -1 0 1 -1 0 1 y 1 z 1 -1 0 1 -1 0 1 y 2 z 2 0 0.2 0.4 0.6 0.8 1 -1 0 1 τ u 1
Fig. 4 Fluid case, step 3: multiple shooting. Solutions 1/12 (left subplot) and 6/12 (right subplot) for λ = 2.0 from Table 3, with the switching times (vertical dashed lines) from direct optimization 20a. 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 arc, applying the algorithm presented in subsection 4.3 from [START_REF] Bonnard | Geometric numerical methods and results in the control imaging problem in nuclear magnetic resonance[END_REF]. 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. 5 represents the computations of the two conjugate points (since the structure is 2BS) of solution 1/12 with λ = 2.0 from Table 3. 3. The rank condition from the algorithm presented in subsection 4.3 from [START_REF] Bonnard | Geometric numerical methods and results in the control imaging problem in nuclear magnetic resonance[END_REF] is evaluated along the two singular arcs. See [START_REF] Bonnard | Second order optimality conditions in the smooth case and applications in optimal control[END_REF] 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.

Step 5: influence of the final time (homotopy method)

Given that the initial point (the North pole) is a stationary point, the constrast 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.

For all solutions from Table 3, provided by multiple shooting, a first homotopy on λ is made in the range [START_REF] Allgower | Introduction to numerical continuation methods[END_REF]3]. Fig. 6 shows the path of zeros in the left subgraph with the value function and the norm of the shooting function, for solution 2/12, λ = 2, from Table 3. The right subgraph presents the value functions of solutions 5/12, 9/12 and 12/12 after the homotopies, showing the branches of different local solutions. 3. The blue color represents 2BS solutions. 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. A first homotopy from λ = 2 to λ = 3 is made. A second from λ = 2 to λ = 1 stops around λ = 1.32 since the norm of the shooting function becomes greater than a threshold fixed to 10 -3 . This threshold is an option given to the homotopic method. The right subgraph displays the value functions after homotopies on solutions 5/12, 9/12 and 12/12 for λ = 2.0 from Table 3. All four homotopies stop before reaching λ = 1. These four path of zeros of local solutions have a distinctive feature. For a fixed λ , the time interval of the second bang arc of the best solution among the four path, is included in the time interval of the second bang arc of the second best solution, and so on.

The possible change of structure along a path of zeros is emphasized in Fig. 7. In the right subgraph of this figure, the path of zeros from Fig. 6 is extended, after the detection of a saturating singular arc (see left subgraph). Note that there is no solution from Table 3 on this extended path. Hence a new branch has been found by this detection. Note also that the 2BS solution are represented in blue, the 3BS in red and the 4BS in green. Then for each solution of each branch the second order necessary condition is checked as in subsection 3.4: 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. 8 presents the simulations for solution 2/12, λ = 2, from Table 3. Finally, we detect some crossing between value functions of different path of zeros as we can see in Fig. 9. 3 is portrayed. Compare to Fig. 6. 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 λ = 1.32. Note that there is no solution from Table 3 on the 3BS part of the extended path. Hence a new branch has been found. Fig. 8 Fluid case, step 5: influence of the final time. Second order necessary condition checked along the extended path from Fig. 7. For all solutions from λ = 1 to λ = 3 are computed the first conjugate times along each singular arc. For λ ∈ [1, 1.32], the structure is 3BS and there are 3 singular arcs. For λ ∈ [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. This is true for all extended path of zeros coming from solutions from Table 3. Note that the magenta crosses, plotted just before λ = 1.32, give the time when a singular arc saturates the constraint |u| = 1. 3. Note that the 3BS solution 1/2 with λ = 1.1 and all the 2BS solutions 1/7 with λ = 1.5, 1/6 with λ = 1.8 and 1/12 with λ = 2 give the same extended path after a change of structure. The crossing is around λ = 1.0484. Thus for λ ≤ 1.0484, the best solution has a 3BS structure of the form δ + δ s δ + δ s δ -δ s (bottom- left subgraph) while for λ ∈ [1.0484, 1.351] the best solution is of the form δ + δ s δ -δ s δ -δ s (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 λ = 1.0484.

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 δ + δ s δ + δ s .

For the fluid case, the left subplot of Fig. 10 represents all the different branches obtained in step 5 by homotopy on λ . The greatest two value functions intersect around t f = 1.048T min . The right subplot shows the sub-optimal synthesis. The best policy is:

δ + δ s δ + δ s for λ ∈ [1.000, 1.006], δ + δ s δ + δ s δ -δ s for λ ∈ [1.006, 1.048], δ + δ s δ -δ s δ -δ s for λ ∈ [1.048, 1.351], δ + δ s δ -δ s for λ ∈ [1.351, 3.000]. (33)
For the blood case, the results are excerpted from [18]. The left subplot of Fig. 11 shows the contrast for five different components of {h = 0}, for final times t f ∈ [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: Fig. 11 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 ∈ (1, 1.294)T min and 3BS of the form δ + δ s δ -δ s δ -δ s for t f ∈ (1.294, 2]T min . In the special case t f = T min , the solution is 2BS and of the form δ + δ s δ + δ s 3.7 Step 6: sub-optimal syntheses compared to global results (lmi method)

δ + δ s δ + δ s for λ = 1.000, δ + δ s for λ ∈ (1.000, 1.294], δ + δ s δ -δ s δ -δ s for λ ∈ [1.
We now apply the lmi method to the contrast problem, 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 4 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 . Both formulations (LP) and (LP') are compared. 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. Formulation (LP') allows for faster computation times at high relaxation orders, which was the principal motivation for its introduction. One can also notice that (LP') leads to sharper bounds than (LP). This is explained by the higher number of linear constraints admitted by the former at a given relaxation order. Indeed, F(q, u) in (LP) is bilinear in its arguments whereas F 0 (q) ± F 1 (q) in (LP') is simply linear. This last formulation therefore admits test functions v of one extra degree at any given relaxation order. Fig. 12 compares the evolution of the upper bounds for different values of t f ∈ [T min , 3T min ], for the blood and fluid cases respectively. For the whole end-time interval, (LP') does produce sharper bounds than (LP), although this is not guaranteed.

Formulation (LP) Formulation (LP') r -J r M N r t r -J r M N r t
Finally, Fig. 13 compares the tightest upper bounds found by the lmi method against 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 ≤ 2 T min .

Application to the homogeneity problem in Magnetic Resonance Imaging

In our study, we have restricted the analysis to the maximum achievable contrast by saturation for a single voxel of the image but in MRI, the quality of the contrast is affected by B 0 and B 1 (Rf-field) inhomogeneity so that the spins in the sample in different spatial positions experience different magnetic fields. This is clear from Bloch equation written in the fixed laboratory frame as:

dM dt = γM ∧ B -R(M)
where M is the magnetization vector, γ is the gyromagnetic ratio, B is the magnetic field which decomposes into B 0 + B 1 (t) where B 0 is the strong magnetic polarizing field in the z-direction and B 1 (t) = (B x (t), B y (t)) is the applied Rf-field in the orthogonal direction, while R(M) is the dissipation of the form:

M x T 2 , M y T 2 , (M z -M 0 ) T 1 ,
due to the relaxation parameters. denoting ω 0 = -γB 0 , u(t) = -γB y (t), v(t) = -γB x (t), and normalizing M to q = (x, y, z) = (M x , M y , M z )/M 0 , the Bloch equation takes the form: dq dt = Aq -R(q), with:

A =   0 -ω 0 u ω 0 0 -v -u v 0   , R(q) =   x/T 2 y/T 2 (z -1)/T 1   .
Expressing (x, y, z) in a moving reference frame:

x = x ′ cos ω 0 ty ′ sin ω 0 t, y = x ′ sin ω 0 t, +y ′ cos ω 0 t, z = z ′ , while the control Rf-field is written in the resonant form: u = u 1 cos ω 0 t + v 1 sin ω 0 t, v = -u 1 sin ω 0 t + v 1 cos ω 0 t, the Bloch equation becomes:

dq ′ dt = A ′ q ′ -R(q ′ ), A ′ =   0 0 u 1 0 0 -v 1 -u 1 v 1 0   .
In general, taking into account B 0 and B 1 inhomogeneity, one obtains in MRI, an ensemble of spin systems [START_REF] Li | Control of inhomogeneous ensembles[END_REF]:

dq ′ (t, s) dt = A s q ′ (t, s) -R(q ′ (t, s)), A s =   0 -ω s ε s u 1 ω s 0 -ε s v 1 -ε s u 1 ε s v 1 0   ,
where ω s ∈ [-B, B] represents the resonance offset due to B 0 inhomogeneity and ε s represents the Rf-inhomogeneity and both are determined experimentally according to spatial position of the spin in the sample.

In particular, in the contrast problem, each vortex corresponds to one of the two spins with different relaxation parameters Λ i = (1/T i 1 , 1/T i 2 ), i = 1, 2. For a simple vortex, one has ω s = 0 while ε s can be normalized to 1, and the maximum achievable contrast by saturation in a given transfer time t f can be numerically computed by the methods presented in this article, direct method (Bocop) vs indirect method (HamPath) or both.

A second step in the analysis is to use a combination of both techniques to improve the homogeneity quality of the image versus the maximum contrast. This leads to introduce an additional optimization problem concerning the ensemble of spin systems: min |q ′ (t f , s)q ′ d (t f , s)| where | • | represents an adequate norm distance between the response of the system and the desired contrast. For instance, considering the single-input case and two voxels, one can take q ′ d (t, s) = (0, c(t f )) where c(t f ) is the minimum of the two contrast by saturation of the voxels.

Conclusion

In this paper, we investigate several methods in optimal control to tackle successfully the contrast problem by saturation in nuclear magnetic resonance.

First of all, the application of geometric control techniques allows for an in-depth classification of the structures of the many local optima. As a result, the solutions in the single input are identified to be concatenations of bang-singular arcs. Those solutions are also valid for the bi-entry case, although additional structures, which we partly characterize, also satisfy first order conditions. The fine characterization of the candidate solutions is the basis for the first numerical method considered, an indirect approach by multiple shooting implemented through the HamPath software. When combined with homotopy techniques, this allows finding a family of solutions in function of the problems relaxation time.

The second numerical method used in the paper is a direct approach implemented by the Bocop software. Several initializations are proposed, as they often lead to different locally optimal solutions. The method often finds quality solutions, in a more straightforward way that the indirect method. In addition, those solutions provide an interesting alternative to regularization techniques for initializing the indirect multiple shootings, to obtain numerical solutions with high accuracy. The candidate trajectories are also continued via homotopy by the indirect method, to identify whole new branches of potential solutions.

The last method considered is a moment optimization technique implemented by the GloptiPoly toolbox, leading to convex semi-definite relaxations of the optimal control problem in the form of Linear Matrix Inequalities. The lower bounds provided by the method are compared against the upper bound given by the local synthesis. We report close agreement between those bounds in the single input case, which validates the candidate solutions as globally optimal.

We also present some preliminary simulations on the bi-input case. They seem to indicate little difference with respect to the single input case, both in terms of control structure and objective value. Complete investigation of the bi-input case remains an open problem.

P 1 : 2 :

 12 Fluid case. Spin 1: Cerebrospinal fluid: T 1 = 2, T 2 = 0.3; Spin 2: Water: T 1 = 2.5 = T 2 . P 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.

P 1 : 2 :

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

1 Fig. 2

 12 Fig.2Fluid case, step 1: Direct method. Solutions 11b (left subplot) and 11c (right subplot) from Table1obtained with Bocop. 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. Note that the first singular arcs, for each control subgraph, do not end at the same times.

Fig. 5

 5 Fig.5Fluid case, step 4: second order conditions. Second order necessary condition check on solution 1/12 with λ = 2.0 from Table3. The rank condition from the algorithm presented in subsection 4.3 from[START_REF] Bonnard | Geometric numerical methods and results in the control imaging problem in nuclear magnetic resonance[END_REF] is evaluated along the two singular arcs. See[START_REF] Bonnard | Second order optimality conditions in the smooth case and applications in optimal control[END_REF] 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.

Fig. 6

 6 Fig. 6 Fluid case, step 5: influence of the final time. On the left subgraph: homotopy on solution 2/12, λ = 2.0 from Table3. The blue color represents 2BS solutions. 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. A first homotopy from λ = 2 to λ = 3 is made. A second from λ = 2 to λ = 1 stops around λ = 1.32 since the norm of the shooting function becomes greater than a threshold fixed to 10 -3 . This threshold is an option given to the homotopic method. The right subgraph displays the value functions after homotopies on solutions 5/12, 9/12 and 12/12 for λ = 2.0 from Table3. All four homotopies stop before reaching λ = 1. These four path of zeros of local solutions have a distinctive feature. For a fixed λ , the time interval of the second bang arc of the best solution among the four path, is included in the time interval of the second bang arc of the second best solution, and so on.

Fig. 7

 7 Fig.7Fluid case, step 5: influence of the final time. On the left subgraph are shown the control laws of solutions at λ = 2 and λ = 1.32 from path from Fig.6. For λ = 1.32, we can see the saturating singular arc around τ = 0.92. The 2BS solution at λ = 1.32 is used to initialize a multiple shooting with a 3BS structure and then to perform a new homotopy from λ = 1.32 to λ = 1. On the right subgraph: complete homotopy on solution 2/12, λ = 2.0 from Table 3 is portrayed. Compare to Fig.6. 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 λ = 1.32. Note that there is no solution from Table3on the 3BS part of the extended path. Hence a new branch has been found.

1 Fig. 9

 19 Fig. 9 Fluid case, step 5: influence of the final time. Crossing bewteen branches coming from 3BS solution 1/2 with λ = 1.1 and 2BS solution 2/12 with λ = 2, after change of structure, from Table3. Note that the 3BS solution 1/2 with λ = 1.1 and all the 2BS solutions 1/7 with λ = 1.5, 1/6 with λ = 1.8 and 1/12 with λ = 2 give the same extended path after a change of structure. The crossing is around λ = 1.0484. Thus for λ ≤ 1.0484, the best solution has a 3BS structure of the form δ + δ s δ + δ s δ -δ s (bottom- left subgraph) while for λ ∈ [1.0484, 1.351] the best solution is of the form δ + δ s δ -δ s δ -δ s (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 λ = 1.0484.

Fig. 12

 12 Fig.12Step 6: LMI. Upper bounds on contrast, for fluid case on left subgraph and blood case on right subgraph, in function of final time for formulations (LP) (straight lines) and (LP') (dashed) for relaxation orders 3 to 5.

Fig. 13 Fig. 14 Fig. 15 Fig. 16

 13141516 Fig.13Step 6: LMI. Best upper bounds (dotted line) by the lmi method compared with best solutions by HamPath (straight), and relative gap between the two, for fluid case on top subgraph and blood case on bottom subgraph.

  we have: i) The Lie algebra generated by {B 1 , B 2 } is a Lie algebra isomorphic to so(3). ii) Assume |δ 1 | and |δ 2 | are nonzero and distinct. Then the Lie algebra generated by {A

  always provides the best solution after shooting.

	λ 1.1	Init. Name Init. Contrast Init. rank 11a 0.6460 1/4	New Contrast New Rank 0.6452 1/2	Structure 1 -1 -1	|S λ | 1.12e-13
	1.1	11a	0.6460	1/4	0.6361	2/2	1 -1 -1	6.02e-14
	1.5	15c	0.6930	1/4	0.6919	1/7	1 -1	5.08e-14
	1.5	15a	0.6899	2/4	0.6888	2/7	1 1	2.03e-14
	1.5	15d	0.6877	3/4	0.6860	3/7	1 -1 -1	2.73e-13
	1.5	15c	0.6930	1/4	0.6837	4/7	1 -1	1.01e-12
	1.5	15c	0.6930	1/4	0.6759	5/7	1 -1	2.76e-14
	1.5	15a	0.6899	2/4	0.6741	6/7	1 1	3.16e-14
	1.5	15a	0.6899	2/4	0.6664	7/7	1 1	3.43e-12
	1.8	18b	0.7025	1/2	0.7014	1/6	1 -1	1.64e-13
	1.8	18a	0.6915	2/2	0.6980	2/6	1 -1 1	1.21e-10
	1.8	18b	0.7025	1/2	0.6927	3/6	1 -1	3.50e-13
	1.8	18a	0.6915	2/2	0.6894	4/6	1 -1 1	2.36e-12
	1.8	18b	0.7025	1/2	0.6844	5/6	1 -1	1.38e-13
	1.8	18b	0.7025	1/2	0.6764	6/6	1 -1	3.02e-13
	2.0	20a	0.7050	1/2	0.7040	1/12	1 -1	2.23e-14
	2.0	20b	0.7017	2/2	0.7005	2/12	1 1	2.56e-14
	2.0	20a	0.7050	1/2	0.6977	3/12	1 -1	5.25e-14
	2.0	20a	0.7050	1/2	0.6951	4/12	1 -1	3.40e-13
	2.0	20b	0.7017	2/2	0.6918	5/12	1 1	1.72e-13
	2.0	20a	0.7050	1/2	0.6912	6/12	1 -1	7.88-07
	2.0	20b	0.7017	2/2	0.6879	7/12	1 1	8.75e-07
	2.0	20a	0.7050	1/2	0.6865	8/12	1 -1	2.47e-13
	2.0	20b	0.7017	2/2	0.6833	9/12	1 1	2.43e-11
	2.0	20a	0.7050	1/2	0.6829	10/12	1 -1	1.69e-13
	2.0	20a	0.7050	1/2	0.6783	11/12	1 -1	2.13e-13
	2.0	20b	0.7017	2/2	0.6751	12/12	1 1	1.10e-12
	3.0	30a	0.7079	1/1	0.7013	1/2	1 -1	7.18e-14
	3.0	30a	0.7079	1/1	0.6923	2/2	1 -1	2.23e-13

  294, 2.000]. 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 subplot. The colors are blue for 2BS structure, red for 3BS and green for 4BS. The best policy is δ + δ s δ + δ s δ -δ s for λ ≤ 1.0484, and δ + δ s δ -δ s δ -δ s for λ ∈ [1.0484, 1.351]. Then, for λ ∈ [1.351, 3], the best policy is 2BS and of the form δ + δ s δ -δ s .

													(34)
		0.7								0.7
		0.68								0.68
		0.66								0.66
	contrast	0.64							contrast	0.64
		0.62								0.62
		0.6								0.6
		0.58								0.58
		1	1.5	2		2.5		3		1	1.5	2	2.5	3
				λ								λ
	0.4 0.42 0.48 0.5 Fig. 10 1 0.44 0.46 Contrast	1.2	1.4	1.6	1.8	2	Contrast	1 0.38 0.4 0.42 0.46 0.48 0.5 0.44	1.5	2		2.5	3
				λ							λ

Table 4

 4 Comparison of formulations (LP) and (LP') for the fluid case with t f = T

	r

min : upper bounds on contrast -J r M , numbers of moments N r and CPU times t r in function of relaxation order r.

where s i = deg(g i )/2 if deg(g i ) is even and (deg(g i ) + 1)/2 otherwise. Similarly,

is even and (deg(g i ) + 1)/2 otherwise. Note that since M r+1 0 implies M r 0, only the highest order moment matrices and localizing matrices are necessary. In addition, as only additional constraints are added at each higher relaxation order, we have the following result: Proposition 10 (Lasserre's relaxation hierarchy)

It is also shown in [START_REF] Lasserre | Nonlinear optimal control via occupation measures and LMI-relaxations[END_REF] that the costs of the lmi relaxations actually converge to that of the cost of (SDP) when r → ∞.

Example: construction of the first relaxation

We conclude the presentation of the method by explicitly constructing the first relaxation of (OCP) for the single input case, under measure formulation (LP). Define the indexing of moments in the monomial basis by y µ i jkℓmn := t i q j 11 q k 12 q ℓ 21 q m 22 u n dµ, y ϕ jkℓm := q j 11 q k 12 q ℓ 21 q m 22 dϕ.

Then, the cost of (LMI) is explicited by

where the decision variables are members of the finite sets {y i jkℓmn : i

For the weak dynamics of (LMI), all polynomial test functions of degree 1 or lower will lead to linear constraints involving moments of order two or less. Plugging v = 1,t, q 11 , . . . , q 22 into (LMI) leads to

The semi-definite constraints for the moment matrices are, applying [START_REF] Lasserre | Positive polynomials and their applications[END_REF]:

All inequalities defining K and Q f in ( 29) and ( 30) are quadratic, hence their corresponding localizing moment matrix given by ( 28) is a scalar at the first relaxation, leading to

Finally, the equalities g 1 (q) = 0 and g 2 (q) = 0 defining Q f can be restated as the following linear equalities: We present here some preliminary results on the bi-input problem in the fluid case, with the direct method Bocop. We use the same constant initialization for the state variables, and for the control variables u x and u y we test the constant initializations -0.5, -0.1, 0.1 and 0.5. We recall that the control must satisfy the constraint |u| ≤ 1 at all times. We use the same regularisation as in the mono-input case, namely adding a term ε Unlike the many local solutions we found for the mono-input problem (cf Table. 1), here the first batch of optimizations seems to reveal only one branch of solutions, with a 2BS structure. We present in more details the solution for t f = 2T min , which gives a contrast of 0.7053. The control and the trajectories of the two spins are drawn on Fig. 14. We observe that during the first BS sequence, both spins stay in the same plane, and the phase of the control stays constant as well. Then during the second bang arc, the phase of the control changes and the spins move into separate planes, before rejoining again at the end of the bang arc. Finally, the second singular arc has again a constant phase for the control, with both spins staying in a plane until the final time. Figs. [START_REF] Caillau | Differential continuation for regular optimal control problems[END_REF] and 16 show the phase of the control and spins, with the norm of the control plotted to indicate the 2BS structure.

Classification of the singular arcs. On this solution, the singular arcs correspond to the Σ ′′ 2 case from Th.1. We numerically check that x 1 y 2y 1 x 2 = 0 over the two singular arcs, meaning that both spins remain in the same plane. We also check that the phase of the control is constant at these times, and is orthogonal to the plane of the spins, which indeed leads to θ = 0 when expressing the dynamics in spherical coordinates (see 1.7.2). The singular arcs in this bi-input solution are actually similar to the singular arcs in the mono-input case. Indeed, if we observe closely the spin trajectories, we find that this bi-input solution is similar to the mono-input solution 20a from Fig. 3. The contrast values confirm this similarity, with 0.7050 for the monoinput and 0.7053 for the bi-input.

When applied to the bi-entry problem, only formulation (LP) is available for the moment/LMI approach, as the quadratic constraint on the controls prevents the use of formulation (LP'). In addition to the increase in the state and control space sizes, this makes only lower relaxations order available for the bi-entry case. In practice, the method is limited to the third relaxation order, for which the upper bounds are essentially identical to those of the mono-entry problem, in both the blood and fluid case. This strongly limits the possible improvements of controlling the problem with two inputs, if any at all.