Table Of ContentConstrained Minimum-Energy Optimal Control of the
Dissipative Bloch Equations
Dionisis Stefanatosa, Jr-Shin Lib
0
1 aPrefecture of Kefalonia, Argostoli, Kefalonia 28100, Greece
0 bWashington University,St. Louis, MO 63130, USA
2
n
a
J
5 Abstract
2
In this letter, we apply optimal control theory to design minimum-energy π/2
] andπ pulsesfortheBlochsysteminthepresenceofrelaxationwithconstrained
C
control amplitude. We consider a commonly encountered case in which the
O transverse relaxation rate is much larger than the longitudinal one so that the
. latter can be neglected. Using the Pontryagin’s Maximum Principle, we de-
h
rive optimal feedback laws which are characterized by the number of switches,
t
a depending on the control bound and the coordinates of the desired final state.
m
Keywords: Maximum Principle, Bloch Equations
[
1
v
1. Introduction
1
6
4 Optimal controltheory [1] has been extensively used recently for the design
4 ofpulsesthatoptimizetheperformanceofvariousNuclearMagneticResonance
1. (NMR)andquantumsystemslimitedbythepresenceofrelaxation[2,3,4,5,6,
0 7, 8, 9, 10, 11, 12, 13], the dissipation due to random interactions between the
0 systemanditsenvironment. Inthisletter,weemploytoolsfromoptimalcontrol
1
toderiveminimum-energyπ/2andπ pulsesforasimpleNMRsystemdescribed
:
v by the Bloch equations. In particular, we study the case where transverse
i relaxationdominates the dissipationofthe systemandthe controlamplitude is
X
bounded.
r
a The problem of optimal control of the Bloch equations and its closely re-
latedcorrespondingproblemforatwo-levelquantumsystemhavereceivedcon-
siderable attention. D’ Alessandro and Dahleh [14] considered the problem of
minimum-energy optimal control for a two-level quantum system without dis-
sipation. Boscain and Mason [15] examined the time minimal problem for a
spin-1/2 particle in a magnetic field neglecting dissipation. Sugny, Kontz and
Jauslin [8], Bonnard and Sugny [9] and Bonnard, Chyba and Sugny [10] stud-
Email addresses: [email protected] (DionisisStefanatos),
[email protected] (Jr-ShinLi)
Preprint submitted toElsevier January 25, 2010
ied extensively the problem of time-optimal control for a dissipative two-level
quantum system.
In our recent work, we studied the problem of designing minimum-energy
π/2 and π pulses for the Bloch system dominated by the transverse relaxation
with unlimited control amplitude [12]. This dissipative system is of great prac-
tical importance as it is a very good approximation for many applications of
interest. Inthisarticle,weextendthispreviousworktoconsiderthecasewhere
the control amplitude is limited, which accounts for realistic limitations of the
experimentalsetupandalsomakesthe problemmoreinteresting fromacontrol
theoretic perspective.
In the next section, we formulate the related optimal control problems of
such pulse designs. The solutions of these problems are presented in Section 3,
which is the main contribution of this article. Then in Section 4, we present
some examples to demonstrate our analytical results.
2. Optimal Control of Dissipative Bloch Systems
In a resonant rotating frame, the Bloch equations with the longitudinal re-
laxation neglected are of the form [16]
z˙ = u x u y (1)
y x
−
x˙ = Rx u z (2)
y
− −
y˙ = Ry+u z, (3)
x
−
where r = (x,y,z) is the magnetization vector, u ,u are the transverse com-
x y
ponents of the magnetic field and R> 0 is the transverse relaxation rate. The
above equations constitute a dissipative bilinear control system. By the follow-
ing change of variables (see Fig. 1(a)) and time rescaling
a = lnr=ln( x2+y2+z2)
tanθ = x2+y2p/z
tanφ = py/x
t = Rt ,
new old
we arrive at a new system
a˙ = sin2θ (4)
−
θ˙ = u sinθcosθ (5)
−
φ˙ = vcotθ (6)
where u = (u /R)sinφ (u /R)cosφ,v = (u /R)cosφ+(u /R)sinφ are the
x y x y
−
normalized components of transverse magnetic field perpendicular and parallel
to r⊥ =(x,y), respectively.
Note thatv does notaffect the angleθ ofthe pulse. It just rotatesr around
z-axis,resulting in a waste of energy. Thus, optimality requires that v =0 and
2
y
z
r r
θ
u φ y
θ
r
⊥
O z
x
(a) u r⊥ (b) u x-axis
⊥ k
Figure 1: The optimal transverse magnetic field u is perpendicular to r⊥ and its phase is
constant (panel a). Without loss of generality, the experimental setup can be arrangedsuch
that u x-axis. In this case, φ= π/2 and r rotates in yz-plane (panel b). For convenience
k
we display z in the horizontal axis and y in the vertical, so that (r,θ) have the common
configurationofpolarcoordinates ontheplane.
hence φ = constant. If for example we choose u to be parallel to the x-axis,
then φ = π/2 and r rotates in yz-plane, see Fig. 1(b). This is the case that
we consider in this letter. Equations (4) and (5) are sufficient to describe this
rotation.
The optimal control problem that we would like to pursue is formulated as
follows. Considerthedynamicalsystemasin(4),(5),startingfrom(r(0),θ(0))=
(1,0) (corresponding to (a(0),θ(0)) = (0,0)) and for a specified final value
a(τ) < a(0) = 0 (equivalent to r(τ) < r(0) = 1), what is the optimal bounded
control u(t) with u(t) m R+, 0 t τ, that accomplishes the transfer
| | ≤ ∈ ≤ ≤
between the above starting point and the target point (a(τ),θ(τ) = π/2orπ),
whileminimizing theenergyE = τ u2(t)/2dt? Notethatforthe transfersthat
0
we study here, it must be m > 1/2, otherwise Eq. (5) reaches an equilibrium
R
point θ0 <π/2. Also, the final time τ is unspecified.
3. Derivation of the Optimal Control
The control Hamiltonian for the addressed problem is
H = u2/2+λ (u sinθcosθ) λ sin2θ (7)
θ a
− − −
whereλ ,λ aretheLagrangemultipliers. AccordingtoPontryagin’smaximum
θ a
principle[1],thenecessaryconditionsforoptimalityof(u(t),θ(t),a(t),λ (t),λ (t))
θ a
are
λ˙ = ∂H/∂θ=λ cos2θ+λ sin2θ (8)
θ θ a
−
3
F
2
F’
2
y I
y
F1 F’1 O z A O z A
(a) Caseofπ pulse (b) Caseofπ/2pulse
Figure 2: The final point F1′ can be reached by either following the minimum-energy path
AF1′, or traveling along AF1 up to point I and then leaving the system relax to F1′ (panel
a). Analogously, F2′ can be reached by either following the minimum-energy path AF2′, or
travelingalongAF2 andthenleavingthesystemrelaxtoF2′ (panel b). .
λ˙ = ∂H/∂a=0 (9)
a
−
u = argmaxH(u,θ,a,λ ,λ ) (10)
θ a
u
From(9),weimmediatelyknowthatλ isaconstant. Additionally,theoptimal
a
(u,θ,a,λ ,λ ) satisfies [1]
θ a
H(u,θ,a,λ ,λ )=0, 0 t τ. (11)
θ a
≤ ≤
Let E = min τ u2(t)/2dt be the minimum cost corresponding to the op-
u 0
timal solution. Using calculus of variations we find that small changes in the
R
a-coordinate of the final point δa , and small changes in the final time δτ,
τ
produce the following change in the minimum cost
δE =λ (τ)δa H(τ)δτ. (12)
a τ
−
Therefore
λ (τ)=∂E/∂a =∂E/∂r dr /da =∂E/∂r r , (13)
a τ τ τ τ τ τ
· ·
where r =eaτ is the radius of the final point.
τ
For the transfers that we examine here it is
∂E/∂r 0 (14)
τ
≥
i.e. thelargerthe finalr, themoreenergyisneededtoquicklyrotatethevector
′
before it dissipates. To see this, we refer to Fig. 2. Let E,E be the minimum
′ ′
energiesnecessarytoreachthefinalpointsF1(rτ,π),F1(rτ,π),respectively,with
′ ′
rτ ≤rτ,seeFig. 2(a). Thecorresponding′ minimum-energypathsareAF1,AF1.
An alternative way to reach the point F1 is the following: travel along AF1 up
4
′
to the point I, with zI = zF1′ = −rτ, then set u =′0′ and wait until dissipation
eliminates the y-coordinate, see (3). The energy E spent for this travel is the
′′ ′
portion of E necessary to reachI, so E E. By the definition of E it is also
′ ′′ ′ ≤ ′
E E , and thus E E. Analogously, let E,E be the minimum energies
≤ ≤ ′ ′
necessary to reach the final points F2(rτ,π/2),F2(rτ,π/2), respectively, with
′
r r again, see Fig. 2(b). The corresponding minimum-energy paths are
τ ≤ τ ′ ′
AF2,AF2. An alternative way to reach the point F2 is the following: travel
along AF2 up to the point F2, then set u= 0 and wait until dissipation brings
′
the system at the point F , see (3). The energy spent for this travel is the
2
′ ′
necessary energy to reach F2, i.e. E. Then, by definition of E , it is E E.
≤
Thus (14) is true and from (13) we have that λ (τ) 0. But λ =constant, so
a a
≥
we can set
λ =κ2/2, κ 0. (15)
a
≥
Having determined the sign of λ , we first examine the case of unbounded
a
control and then we use the developed intuition to study the general case of
bounded control.
3.1. Unbounded Control
Whenthecontroluisunbounded,thenfrom(10)weconcludethat∂H/∂u=
0. This condition and Eq. (7) yield
u=λ . (16)
θ
Using (16) and (15), the condition (11) becomes
λ2 2λ sinθcosθ κ2sin2θ =0. (17)
θ− θ −
The optimal u is then given by the following feedback law
u(θ)=λ =sinθ(cosθ+ cos2θ+κ2). (18)
θ
p
Note that only the positive solution of the quadratic equation has physical
meaning (corresponds to increasing θ) for the transfers that we study here.
Using (16) and (18), the validity of (8) can be easily verified. Inserting (18) in
(5) we obtain the differential equation for the optimal trajectory
θ˙ =sinθ cos2θ+κ2 (19)
p
Eliminating time between (4) and (19) we obtain
da sinθ
= . (20)
dθ −√cos2θ+κ2
Integratingtheaboveequationfromthestartingpoint(0,0)tothepoint(lnr,θ),
we find the optimal trajectory
cosθ+√cos2θ+κ2
r(θ)= (21)
1+√1+κ2
5
2.5 0.8
0.7
2
0.6
0.5
1.5
R)
u ( y0.4
1 0.3
0.2
0.5
0.1
00 0.5 θ (rad) 1 1.5 00 0.1 0.2 0.3 0.4 0z.5 0.6 0.7 0.8 0.9 1
(a) Optimalπ/2pulse (b) Optimaltrajectory
4
3.5
1
3
2.5 0.8
R)
u ( 2 y0.6
1.5
0.4
1
0.2
0.5
00 0.5 1 θ (1.r5ad) 2 2.5 3 −00.6 −0.4 −0.2 0 0z.2 0.4 0.6 0.8 1
(c) Optimalπ pulse (d) Optimaltrajectory
Figure3: minimum-energyπ/2(panela)andπ(panelc)pulsesforrτ =0.6. Thecorrespond-
ingtrajectoriesarealsoshown(panelsb,d).
Setting θ = π/2,π in the above equation, we find the optimal κ for the π/2
τ
and π pulses, as a function of the radius r of the final point
τ
2rτ 2√rτ
κ = , κ = (22)
π/2 1 r2 π 1 r
− τ − τ
The energy of the optimal pulses is calculated from
τ u2(t) θτ u2(θ)
E = dt= dθ (23)
0 2 0 2θ˙(θ)
Z Z
using (18) and (19). The results for θ =π/2,π are
τ
1 1+r
τ
E = , E = (24)
π/2 1 r2 π 1 r
− τ − τ
Using (24) in (13), it is easy to verify the validity of (15) with κ given by (22).
In Fig. 3 we plot the optimal π/2 and π pulses for r = 0.6, as well as the
τ
corresponding trajectories. Observe that optimal u(θ) is small close to θ =0,π
6
(z-axis), directions that are protected against relaxation, while it is large close
to θ =π/2(y-axis),where dissipationis maximizedandthus r mustbe rotated
faster.
3.2. Bounded Control
We now move on to the case where the control is bounded, i.e. u m,
| | ≤
withm>1/2aspointedoutbefore. ThecontrolHamiltonian(7)isaquadratic
formwithrespecttouthattakesitsmaximumvalueatu=λ ,if λ m,and
θ θ
| |≤
at the boundary point u = m if λ > m. The other boundary point, u = m,
θ
−
corresponds to decreasing θ and has no physical meaning for the transfers that
we examine here. Initially, the situation is as in the previous case where the
optimality condition u = λ holds and the optimal control is given by (18).
θ
Angle θ increases following (19) and u,λ change accordingly. Now suppose
θ
that at some point the control reaches the maximum allowable value m. From
(18) we see that this happens at the angles that satisfy the equation
sinθ(cosθ+ cos2θ+κ2)=m, (25)
p
which is equivalent to the following quadratic equation for cotθ
2 κ2
cot2θ cotθ+1 =0. (26)
− m − m2
If θ1,θ2 are the solutions of (26) in [0,π], then it is easy to show that for
θ (θ1,θ2) the following inequality holds
∈
sinθ(cosθ+ cos2θ+κ2)>m. (27)
p
In the interval θ (θ1,θ2) where the above inequality is true, the relation
∈
u = λ gives u(θ) > m which is not permissible. Thus, the optimal control in
θ
this interval is u(θ)=m. From (11) we can find the Lagrange multiplier λ for
θ
the same interval, which is
m2+κ2sin2θ
λ (θ)= (28)
θ
2(m sinθcosθ)
−
and λθ(θ1) = m from (16). This λθ satisfies the optimality condition (8) with
θ evolving in time according to (5) and u(θ)=m. It is not hard to verify using
(27), (28) that λθ > m for θ (θ1,θ2), thus u(θ) = m is indeed the optimal
∈
control in this interval. We call the solutions θ1,θ2 of (26), where a change in
the optimal control occurs, the switching angles. After the second switching,
the optimal control is given again by (18), with the same κ as in the initial
phase, since λ =κ2/2 is constant along the optimal trajectory.
a
As explained above, the switching angles determine the optimal feedback
law. For m 1, Eq. (26) has real solutions for κ √m2 1. For 1/2<m<1
≥ ≥ −
it has real solutions for every κ 0. We examine separately these two cases
≥
7
D
1
D B
2
D
3
y
C C O A
1 2 z
Figure 4: Switching curves AB and BC1 for the case m 1. The outermost curve AC1
≥
defines the reachable set from the starting point A(1,0) when u m. Curve BC2 is the
optimaltrajectoryforκ=√m2 1andθ [θB,π]. Theperpendic≤ularaxisθ=π/2crosses
− ∈
thecurvesAC1,BC1,BC2 atthepointsD1,D2,D3,respectively.
3.2.1. m 1
≥
For κ>√m2 1 there are two switching angles, given by
−
1 √κ2 m2+1
θ1,2 =cot−1 ± − (29)
m
!
In this letter the range of the function cot−1 is considered to be [0,π], so
cot−1(x)=π cot−1( x), x<0. (30)
− −
In the case κ = √m2 1, the two angles obtain the common value θ =
B
cot−1(1/m). The angle−of the first switching point is in the range θ1 [0,θB],
∈
whilethatofthe secondθ2 [θB,π]. Wefindtheequationofthe firstswitching
∈
curve, i.e. the curve composed by the points (r1,θ1). Before the switching the
optimal control is given by (18), so each of the points (r1,θ1) belongs to an
optimal curve of the form (21), i.e.
cosθ1+√cos2θ1+κ2
r1(θ1)= . (31)
1+√1+κ2
But this κ is related to the switching angle θ1 through (29), which we re-write
as
κ2 =(mcotθ1 1)2+m2 1 (32)
− −
sinceonlyκ2 appearsin(31). Thesetwoequationsdeterminethefirstswitching
curve, with the angle θ1 in the interval θ1 [0,θB].
∈
After the first switching, the optimal control takes the value u(θ) = m. It
maintains this value until the second switching, at angle θ2. Between the two
switchings the evolution is given by
da sin2θ
= , (33)
dθ −m sinθcosθ
−
8
as we derive from (4), (5) with u=m. Integrating the above equation from θ1
to θ2 we find the radius of the second switching point
2m sin2θ1 f(θ1,θ2)
r2(θ2)=r1(θ1) − exp (34)
r2m−sin2θ2 (cid:20)−√4m2−1(cid:21)
where
f(θ1,θ2)=cot−1 2mcotθ2−1 cot−1 2mcotθ1−1 (35)
√4m2 1 − √4m2 1
(cid:18) − (cid:19) (cid:18) − (cid:19)
and
θ2 =cot−1(2/m cotθ1). (36)
−
Therefore,to everyfirstswitching point(r1,θ1) correspondsa secondswitching
point (r2,θ2) with angle given by (36) and radius given by (34), (35). These
points compose the second switching curve. The two switching curves are plot-
ted in Fig. 4, curves AB and BC1. The joint point is B(rB,θB), where
√m2+1
r = . (37)
B
m+1
The outermost curve AC1 corresponds to the trajectory traveledfor u(θ)=
m,θ [0,π]. Theequationforthis trajectorycanbe foundbysetting(r1,θ1)=
∈
(1,0) (starting point A) in (34), (35). It is
r3(θ3)= 2m exp 1 cot−1 2mcotθ3−1 (38)
r2m−sin2θ3 (cid:20)−√4m2−1 (cid:18) √4m2−1 (cid:19)(cid:21)
whereweused(r3,θ3)todenoteapointonthecurveandθ3 [0,π]. Thepoints
∈
between this curve and the horizontal axis define the reachable set from A for
a specific control bound m. This curve meets the axes θ =π,π/2 at the points
C1(rC1,π),D1(rD1,π/2), where
π
r = exp (39)
C1 −√4m2 1
(cid:18) − (cid:19)
1 1
r = exp π cot−1 (40)
D1 −√4m2 1 − √4m2 1
(cid:20) − (cid:18) − (cid:19)(cid:21)
Thesearethepointswiththelargestradiusalongtheseaxes,thatcanbereached
from A(1,0). Using (31), (32), (34), (35) and (36) we find that the second
switching curve BC1 crosses the axis θ =π/2 at the point D2(rD2,π/2) where
√m2+2 1 m2 1
r = exp cot−1 − (41)
D2 1+√m2+1 −√4m2 1 √4m2 1
(cid:20) − (cid:18) − (cid:19)(cid:21)
As we mentioned above, switching takes place only for κ > √m2 1. For
−
κ √m2 1 there is no switching and the optimal trajectory is given by (21).
≤ −
9
Forthesevaluesofκ,theoptimaltrajectorycrossestheaxisθ =π,π/2atpoints
with radius
√κ2+1 1 κ
rπ = √κ2+1+−1, rπ/2 = √κ2+1+1 (42)
respectively, both increasing functions of κ 0. In Fig. 4 we plot the optimal
≥
trajectory without switching for the largest permissible value κ=√m2 1
−
cosθ0+√cos2θ0+m2 1
r0(θ0)= − (43)
m+1
andforθ0 ∈[θB,π]. Itcrossestheaxesθ =π,π/2atthepointsC2(rC2,π),D3(rD3,π/2)
where
m 1 √m2 1
r = − , r = − (44)
C2 m+1 D3 m+1
Thesearethelargestradiuspointsalongtheseaxes,thatcanbereachedwithout
switching.
Usingthe constructionshowninFig. 4wecanfindthe switchingpointsand
the optimal control for any final point of the form F1(rτ,π) or F2(rτ,π/2). If
F1 C1C2 then the optimal trajectory after the second switching is
∈
cosθ+√cos2θ+κ2
r(θ)=r2(θ2) (45)
cosθ2+√cos2θ2+κ2
where
κ2 =(1 mcotθ2)2+m2 1 (46)
− −
andS2(r2,θ2) is the secondswitching point. Eq. (45) is found after integrating
(20) from S2 to F1, while (46) holds because S2 BC1, see (29). The final
∈
point F1(rτ,π) belongs to this curve, so we find
cosθ2+√cos2θ2+κ2
r2(θ2)=rτ 1+√1+κ2 (47)
−
Plotting this curve for θ2 [θB,π], it crosses the second switching curve at the
∈
second switching point S2. The first switching angle can be found from
θ1 =cot−1(2/m cotθ2). (48)
−
Having determined the two switching angles and the optimal κ, the optimal
feedbackcontrolis alsodetermined. IfF1 C2O thennoswitchingis necessary
∈
and the optimal control is given by (18) with κ=κ given by (22).
π
ForF2 D1D2 thereisonlyoneswitching. Theoptimaltrajectoryafterthe
∈
switching is
2m sin2θ1 f(θ1,θ)
r(θ)=r1(θ1) − exp (49)
2m sin2θ −√4m2 1
r − (cid:20) − (cid:21)
10
