Article Contents
Article Contents

# Time minimal saturation of a pair of spins and application in Magnetic Resonance Imaging

• * Corresponding author

The first author is partially supported by the ANR Project - DFG Explosys (Grant No. ANR-14-CE35-0013-01;GL203/9-1), the fourth author is supported by the Austrian FWF grant F5004, and the four authors are supported by the FMJH PGMO and from the support of EDF, Thales and Orange

• In this article, we analyze the time minimal control for the saturation of a pair of spins of the same species but with inhomogeneities of the applied RF-magnetic field, in relation with the contrast problem in Magnetic Resonance Imaging. We make a complete analysis based on geometric control to classify the optimal syntheses in the single spin case to pave the road to analyze the case of two spins. The ${\texttt {Bocop}}$ software is used to determine local minimizers for physical test cases and Linear Matrix Inequalities approach is applied to estimate the global optimal value and validate the previous computations. This is complemented by numerical computations combining shooting and continuation methods implemented in the ${\texttt {HamPath}}$ software to analyze the structure of the time minimal solution with respect to the set of parameters of the species. Symbolic computations techniques are used to handle the singularity analysis.

Mathematics Subject Classification: 49K15, 49M15, 81Q93, 14Q20.

 Citation:

• Figure 1.  In this example, $(\gamma, \Gamma) = (0.12, 0.5)$. (Left) Collinearity set, singular set and visualization of $N$, $O$, $O_1$, $S_1$, $S_1'$, $S_3$ and $S_3'$. (Right) Sign of $\alpha(q)$ in the domain $y\le 0$. The sign in the domain $y\ge 0$ is given by symmetry

Figure 2.  The points $S_1$, $S_2$ and $S_3$ with different orders. Left sub-graph: $(\gamma,\Gamma) = (0.1, 0.5)$ and right: $(\gamma,\Gamma) = (0.163,0.5)$

Figure 6.  Schematic time minimal synthesis to steer a single spin system from the North Pole N to any point of the Bloch ball in the reachable set, for parameters $(\Gamma,\gamma) \in C$. An arbitrary zoom has been used to construct the figure. The set of $\Sigma_i$ forms the switching surface $\Sigma$ dividing the $+1$ and $-1$ areas respectively in red and blue. The minimal time trajectory to steer the spin from N to O is $NS_1S_2S_2'O$, i.e. it is of the form $\sigma _+^N \sigma _s^h \sigma _+^b \sigma _s^v$ with horizontal $\sigma _s^h$ and vertical $\sigma _s^v$ singular arcs. The spin leaves the horizontal singular arc before the point $S_3$ (where the control saturates the constraint) producing a bridge $\sigma _+^b$ to reach the vertical singular line

Figure 7.  Schematic time minimal synthesis to steer a single spin system from the North Pole N to any point of the Bloch ball in the reachable set, for parameters $(\Gamma,\gamma) \in B$. An arbitrary zoom has been used to construct the figure. The set of $\Sigma_i$ forms the switching surface $\Sigma$ dividing the $+1$ and $-1$ areas respectively in red and blue. The minimal time trajectory to steer the spin from N to O is $NS_1'O$, i.e. it is of the form $\sigma _+^N \sigma _s^v$

Figure 8.  Schematic time minimal synthesis to steer a single spin system from the North Pole N to any point of the Bloch ball in the reachable set, for parameters $(\Gamma,\gamma) \in A_2$. For $(\Gamma,\gamma) \in A_1$, the synthesis is the same but there is no horizontal singular line inside the Bloch ball. An arbitrary zoom has been used to construct the figure. The minimal time trajectory to steer the spin from N to O is $NS_1'O$, i.e. it is of the form $\sigma _+^N \sigma _s^v$

Figure 3.  The domain of interest of the parameters in white with the sub-domains $A_1$, $A_2$, $B$ and $C$

Figure 4.  The sub-domains in $(\Gamma,\gamma/\Gamma)$ coordinates. One can notice that when $\Gamma$ tends to $0$, then the slopes of the curves $S_1 = S_2$ and $S_1 = S_3$, in $(\Gamma,\gamma)$ coordinates, are equal to $2/3$

Figure 5.  The slopes $T_2/T_1 = \gamma/\Gamma$ for the species from Table 1 with the particular case when $\omega_\mathrm{max} = 2 \pi\times 32.3$ Hz. This case is represented by a bullet while the slope is represented by a line

Figure 14.  Water case: $\rho = \bar{\rho}$, $\theta = 0.7854$ and $\varepsilon = 0.1$. Trajectories of spins 1 and 2, and control

Figure 9.  Deoxygenated blood case ($C_1$) with RF-inhomogeneity ($\varepsilon = 0.1$). Trajectories for spin 1 and 2 in the (y, z)-plane are portrayed in the first two subgraphs. The corresponding control is drawn in the right subgraph. The horizontal lines $z_1 = z_2 = z_s = \gamma/(2\delta)$, $\delta = \gamma-\Gamma$, is represented by dashed lines. Note that the last bang arc is not well captured by the direct solver because it is too short

Figure 10.  Water case ($C_4$) with RF-inhomogeneity ($\varepsilon = 0.1$). Trajectories for spin 1 and 2 in the (y, z)-plane are portrayed in the first two subgraphs

Figure 11.  Saturation problem of one spin and two spins for the cases $C_1$, $C_2$, $C_3$ and $C_4$. Relative error $\mathrm{err}(r) = (t_f-T_{LMI'}^r)/t_f$ where $r$ is the order of relaxation, $T_{LMI'}^r$ is the optimal value of (59) using the formulation (57) and $t_f$ is the final time computed with $\texttt{Bocop}$

Figure 12.  The relaxation parameters are given by $\rho = \bar{\rho}$ and (Left) $\theta = \theta_\mathrm{1a,1}^+$, (Right) $\theta = \theta_\mathrm{1a,1}^-$. Each subgraph represents the graph of the switching function $H_1$ along the extremal. For $\theta = \theta_\mathrm{1a,1}^-$, we observe that $H_1$ crosses $0$ twice. A singular arc must be added to continue the homotopy on $\theta$. A zoom in on the graph of $H_1$ is given and may be located thanks to the vertical red lines. Note that we add a singular arc, since in this case the singular extremal is time-minimizing for small time

Figure 13.  The homotopies H1a, H1b, H2 and H3 are presented respectively in blue, red, black and black. The homotopies H2 and H3 are presented only on the top-right subgraph. The plain and dashed lines distinguish the different structures. The top-left subgraph gives the cost (i.e. the final time $t_f$) with respect to $\theta$ for the homotopies H1a and H1b. One can see (it is more visible in the zoom) the intersection of the two paths of zeros in terms of cost. H1b is better for small value of $\theta$ and H1a is better for greater values. One can notice that this two paths of zeros are distinct from the bottom-left subgraph. This subgraph gives the norm of the initial adjoint vector with respect to the homotopic parameter. The norm of the shooting function along the paths H1a and H1b is given in the bottom-left subgraph.% One can see the very good accuracy. The strategies from H1a and H1b are compared with homotopies H2 and H3 on the top-right subgraph

Figure 15.  Fat case: $\rho = \bar{\rho}$, $\theta = 0.4636$ and $\varepsilon = 0.1$. Trajectories of spins 1 and 2, and control

Figure 16.  Fluid case: $\rho = \bar{\rho}$, $\theta = 0.1489$ and $\varepsilon = 0.1$. Trajectories of spins 1 and 2, and control

Figure 17.  H2: $\rho = \bar{\rho}$, $\theta = 0.2$ and $\varepsilon = 0.1$. Trajectories of spins 1 and 2, control and $H_1$

Table 1.  Matter name with relaxation times in seconds, ratio $T_2/T_1$ and value $\theta = {\rm atan}(\gamma/\Gamma)$

 Name $T_1$ $T_2$ $\frac{T_2}{T_1} = \frac{\gamma}{\Gamma}$ $\theta={\rm atan}(\frac{\gamma}{\Gamma})$ Water 2.5 2.5 1.0 0.7854 Fat 0.2 0.1 0.5 0.4636 Cerebrospinal Fluid 2.0 0.3 0.15 0.1489 Oxygenated blood 1.35 0.2 0.1481 0.1471 White cerebral matter 0.78 0.09 0.1154 0.1148 Gray cerebral matter 0.92 0.1 0.1087 0.1083 Brain 1.062 0.052 0.0490 0.0489 Deoxygenated blood 1.35 0.05 0.0370 0.0370 Parietal muscle 1.2 0.029 0.0242 0.0242

Table 2.  Cases treated numerically corresponding respectively to the Deoxygenated case ($C_1$), the Oxygenated case ($C_2$), the Cerebrospinal fluid case ($C_3$) and the Water case ($C_4$). The 5th (resp. 4th) column gives the final time found by ${\texttt {Bocop}}$ for the saturation of one spin (resp. two spins with $B_1$-inhomogeneity). The parameter $\varepsilon$ is fixed to $0.1$

 Case $\Gamma$ $\gamma$ $t_f$ (2 spins) $t_f$ (1 spin) $C_1$ 9.855$\times 10^{-2}$ 3.65 $\times 10^{-3}$ 44.769 42.685 $C_2$ 2.464$\times 10^{-2}$ 3.65 $\times 10^{-3}$ 113.86 110.44 $C_3$ 1.642$\times 10^{-2}$ 2.464$\times 10^{-3}$ 168.32 164.46 $C_4$ 9.855$\times 10^{-2}$ 9.855$\times 10^{-2}$ 15.0237 8.7445

Table 3.  Single spin saturation for the case $C_2$. $T_{LMI}^r$ is the optimal value of the hierarchy (59). Likewise, $T_{LMI'}^r$ is the optimal value of the hierarchy of SDP problems derived from (57). $t_f = 110.44$ is the best time found by the ${\texttt {Bocop}}$ software. The second and fourth columns are the relative errors between best solution found by the ${\texttt {Bocop}}$ software and the one found by moment/LMI techniques for each relaxation order $r$

 1st Formulation 2nd Formulation r $N_m$ $(t_f-T_{LMI}^r)/t_f$ $N_m$ $(t_f-T_{LMI'}^r)/t_f$ 1 25 0.8143 30 0.818 2 105 0.5164 105 0.5958 3 294 0.2611 252 0.4355 4 660 0.1491 495 0.1842 5 1287 0.0932 858 0.1284 6 2275 0.0643 1365 0.096 7 3740 0.0517 2040 0.0797 8 5814 0.0461 2907 0.0716

Table 4.  The homotopies are detailed on each line. The first column gives the name of the homotopy. The second colum "Init" gives the initial value of $\theta$, the associated structure with the reference to the figure which presents the trajectory with the control and the switching function. For this initial value, the solution is obtained from direct method and multiple shooting. The column "Transition" gives the same details but when a first change in the structure is detected during the homotopy thanks to the monitoring. The last column "End" gives again the same details at the end of the homotopy. For the homotopies H1a and H1b, the end occurs when a second change in the structure is detected while for H2 and H3, the end occurs when the final value of $\theta$ is reached. The red and green colors may help to understand the change in the structure, which is explained more in details in this section.

 Name Init Transition End H1a θmax ≈ 1.1071 θ1a, 1 ≈ 0.5069 θ1a, 2 ≈ 0.2722 > θmin = 0.02 σ-σsσ+σ0 σ-σsσ-σsσ+σ0 σ-σsσ-σ+σ0 H1b θmin θ1b, 1 ≈ 0.2752 θ1b, 2 ≈ 0.3018 < θmax σ-σsσ+σsσ+σ0 σ-σsσ+σ0 σ-σs -σ+σ0 H2 θ = 0.2 θ2, 1 ≈ 0.3618 θ = atan(1.5) ≈ 0.9828 σ+σsσ+σsσ+σ0 σ+σsσ+σ0 σ+σsσ+σ0 H3 θ = 0.25 θ3, 1 ≈ 0.4778 θ = θmax σ+σsσ+σsσ+σ0 σ+σsσ+σ0 σ+σsσ+σ0
•  [1] E. Allgower and K. Georg, Introduction to Numerical Continuation Methods, vol. 45 of Classics in Applied Mathematics, Soc. for Industrial and Applied Math., Philadelphia, PA, USA, 2003. doi: 10.1137/1.9780898719154. [2] J. T. Betts, Practical Methods for Optimal Control Using Nonlinear Programming, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, 2001. [3] F. J. Bonnans, P. Martinon and V. Grélard, Bocop - A Collection of Examples, Technical report, INRIA (2012) RR-8053. [4] B. Bonnard, J.-B. Caillau and E. Trélat, Second order optimality conditions in the smooth case and applications in optimal control, ESAIM: COCV, 13 (2007), 207-236.  doi: 10.1051/cocv:2007012. [5] B. Bonnard and M. Chyba, Singular Trajectories and Their Role in Control Theory, vol 40 of Mathematics and Applications, Springer-Verlag, Berlin, 2003. [6] B. Bonnard, M. Chyba, A. Jacquemard and J. Marriott, Algebraic geometric classification of the singular flow in the contrast imaging problem in nuclear magnetic resonance. Mathematical Control and Related Fields-AIMS, Special issue in the honor of Bernard Bonnard. Part II., 3 (2013), 397-432.  doi: 10.3934/mcrf.2013.3.397. [7] B. Bonnard, M. Chyba and J. Marriott, Singular trajectories and the contrast imaging problem in nuclear magnetic resonance, SIAM J. Control Optim., 51 (2013), 1325-1349.  doi: 10.1137/110833427. [8] B. Bonnard, M. Claeys, O. Cots and P. Martinon, Geometric and numerical methods in the contrast imaging problem in nuclear magnetic resonance, Acta Appl. Math., 135 (2015), 5-45.  doi: 10.1007/s10440-014-9947-3. [9] B. Bonnard and O. Cots, Geometric numerical methods and results in the control imaging problem in nuclear magnetic resonance, Math. Models Methods Appl. Sci., 24 (2014), 187-212.  doi: 10.1142/S0218202513500504. [10] B. Bonnard, O. Cots, S. Glaser, M. Lapert, D. Sugny and Y. Zhang, Geometric optimal control of the contrast imaging problem in nuclear magnetic resonance, IEEE Trans. Automat. Control, 57 (2012), 1957-1969.  doi: 10.1109/TAC.2012.2195859. [11] B. Bonnard, O. Cots, J. Rouot and T. Verron, Working Notes on the Time Minimal Saturation of a Pair of Spins and Application in Magnetic Resonance Imaging, http://hal.archives-ouvertes.fr/hal-01721845/ [12] B. Bonnard, J.-C. Faugére, A. Jacquemard, M. Safey El Din and T. Verron, Determinantal sets, singularities and application to optimal control in medical imagery, Proceedings of the 2016 ACM International Symposium on Symbolic and Algebraic Computation, 103–110, ACM, New York, 2016. [13] B. Bonnard and I. Kupka, Théorie des singularités de l'application entrée/sortie et optimalité des trajectoires singulières dans le problème du temps minimal, Forum Math., 5 (1993), 111-159.  doi: 10.1515/form.1993.5.111. [14] U. Boscain and B. Piccoli, Optimal Syntheses for Control Systems on 2-D Manifolds, Springer SMAI, 2004. [15] R. Bulirsch and J. Stoer, Introduction to Numerical Analysis, vol.12 of Texts in Applied Mathematics, Springer-Verlag, New York, 2$^{nd}$ edition, 1993. doi: 10.1007/978-1-4757-2272-7. [16] J.-B. Caillau, O. Cots and J. Gergaud, Differential continuation for regular optimal control problems, Optimization Methods and Software, 27 (2011), 177-196.  doi: 10.1080/10556788.2011.593625. [17] M. Claeys, J. Daafouz and D. Henrion, Modal occupation measures and {LMI} relaxations for nonlinear switched systems control, Automatica, 64 (2016), 143-154.  doi: 10.1016/j.automatica.2015.11.003. [18] S. Conolly, D. Nishimura and A. Macovski, Optimal control solutions to the magnetic resonance selective excitation problem, IEEE Trans. Med. Imaging, 5 (1986), 106-115. [19] M. Gerdts, Optimal Control of ODEs and DAEs, ed. De Gruyter, Berlin, 2012. doi: 10.1515/9783110249996. [20] D. Henrion, J. B. Lasserre and J. Löfberg, GloptiPoly 3: Moments, optimization and semidefinite programming, Optim. Methods and Software, 24 (2009), 761-779.  doi: 10.1080/10556780802699201. [21] A. J. Krener, The high order maximal principle and its application to singular extremals, SIAM J. Control Optim., 15 (1977), 256-293.  doi: 10.1137/0315019. [22] I. Kupka, Geometric theory of extremals in optimal control problems. Ⅰ. The fold and Maxwell case, Trans. Amer. Math. Soc., 299 (1987), 225-243.  doi: 10.2307/2000491. [23] M. Lapert, Y. Zhang, M. Braun, S. J. Glaser and D. Sugny, Singular extremals for the time-optimal control of dissipative spin 1/2 particles, Phys. Rev. Lett., 104 (2010), 083001. [24] M. Lapert, Y. Zhang, M. A. Janich, S. J. Glaser and D. Sugny, Exploring the physical limits of saturation contrast in magnetic resonance imaging, Sci. Rep., 589 (2012). [25] J.-B. Lasserre,  Moments, Positive Polynomials and Their Applications, Imperial College Press, London, 2010. [26] J.-B. Lasserre, D. Henrion, C. Prieur and E. Trélat, Nonlinear optimal control via occupation measures and LMI-relaxations, SIAM J. Control Optim., 47 (2008), 1643-1666.  doi: 10.1137/070685051. [27] M. H. Levitt, Spin Dynamics: Basics of Nuclear Magnetic Resonance, John Wiley and Sons, New York-London-Sydney, 2008. [28] L. Markus, Quadratic differential equations and non-associative algebras, Princeton Univ. Press, Princeton, N.J., 5 (1960), 185-213. [29] H. Maurer, Numerical solution of singular control problems using multiple shooting techniques, J. Optim. Theory Appl., 18 (1976), 235-257.  doi: 10.1007/BF00935706. [30] A. Mosek, The Mosek Optimization Toolbox for Matlab Manual, Version 7.1 Revision, 2015. [31] J. Nocedal and S. J. Wright, Numerical Optimization, Springer-Verlag, New York, 1999. doi: 10.1007/b98874. [32] L. S. Pontryagin, V. G. Boltyanskiĭ, R. V. Gamkrelidze and E. F. Mishchenko, The Mathematical Theory of Optimal Processes, Translated from the Russian by K. N. Trirogoff, edited by L. W. Neustadt, Interscience Publishers John Wiley and Sons, Inc., New York-London, 1962. [33] M. Putinar, Positive polynomials on compact semi-algebraic sets, Indiana Univ. Math. J., 42 (1993), 969-984.  doi: 10.1512/iumj.1993.42.42045. [34] H. Schättler, The local structure of time-optimal trajectories in dimension three under generic conditions, SIAM J. Contr. Opt., 26 (1988), 899-918.  doi: 10.1137/0326050. [35] H. J. Sussmann, Time-optimal control in the plane,, in Feedback Control of Linear and Nonlinear Systems, (Bielefeld/Rome, 1981), vol. 39 of Lecture Notes in Control and Inform. Sci., Springer, Berlin, (1982), 244–260. doi: 10.1007/BFb0006833. [36] H. J. Sussmann, Regular synthesis for time-optimal control of single-input real-analytic systems in the plane, SIAM J. Control and Opt., 25 (1987), 1145-1162.  doi: 10.1137/0325062. [37] E. Van Reeth, H. Ratiney, M. Tesch, D. Grenier, O. Beuf, S. J. Glaser and D. Sugny, Optimal control design of preparation pulses for contrast optimization in MRI, J. Magn. Reson., 279 (2017), 39-50.

Figures(17)

Tables(4)