# Optimal Control

This section provides a brief introduction to the package's interface for sequence optimization. We use the Cramer-Rao bound (CRB) to assess a sequence's performance and optimize the amplitudes ($ω_1$) and durations ($T_\text{RF}$) of RF-pulses to reduce the CRB, assuming a Balanced Hybrid-State Free Precession Pulse Sequence. For computational efficiency, the derivatives of the CRB wrt. $ω_1$ and $T_\text{RF}$ are calculated with the adjoint state method common in the optimal control literature.

For this tutorial, we use the following packages:

```
using MRIgeneralizedBloch
using LinearAlgebra
BLAS.set_num_threads(1) # single threaded is faster in this case
using Optim # provides the optimization algorithm
using Plots
```

Here, we optimize the pulse sequence for a predefined set of parameters:

```
m0s = 0.15
R1f = 0.5 # 1/s
R2f = 15 # 1/s
Rx = 30 # 1/s
R1s = 3 # 1/s
T2s = 10e-6 # s
ω0 = 0 # rad/s
B1 = 1; # in units of B1_nominal
```

and we optimize

`Npulses = 200;`

pulses, spaced

`TR = 3.5e-3; # s`

apart. The cyle duration of

`Npulses * TR`

`0.7000000000000001`

seconds is shorter than the optimal duration, which is in the range of 4-10s. We here use a small `Npulses`

to speed up the computations. The Linear Approximation of the generalized Bloch model is precomputed with

`R2slT = precompute_R2sl();`

In the calculation of the CRB, we account for following gradients:

`grad_list = (grad_m0s(), grad_R1f(), grad_R2f(), grad_Rx(), grad_R1s(), grad_T2s(), grad_ω0(), grad_B1());`

and we sum up the CRB of all parameters, weighted by the following vector:

`weights = transpose([0, 1, 0, 0, 0, 0, 0, 0, 0]);`

Note that the vector `weights`

has one more entry compared to the `grad_list`

vector, as the first derivative is always wrt. $M_0$, regardless of `grad_list`

. Here, we only optimize for the CRB of $m_0^s$, while accounting for a fit of all 9 model parameters.

We take some initial guess for the pulse train:

`α = abs.(sin.((1:Npulses) * 2π/Npulses));`

initialize with a constant `TRF = 300μs`

:

`TRF = 300e-6 .* one.(α);`

and define the first RF pulse as a 500μs inversion pulse by modifying vectors accordingly and by creating a bit vector that indicates the position of the inversion pulse:

```
α[1] = π
TRF[1] = 500e-6
isInversionPulse = [true; falses(length(α)-1)];
```

We note that inversion pulses are not optimized by this toolbox. We calculate the initial $ω_1$

`ω1 = α ./ TRF;`

and plot the initial control:

```
p1 = plot(TR*(1:Npulses), α ./ π, ylabel="α/π")
p2 = plot(TR*(1:Npulses), 1e6TRF, ylim=(0, 1e3), xlabel="t (s)", ylabel="TRF (μs)")
p = plot(p1, p2, layout=(2, 1), legend=:none)
```

With above defined weights, the function `MRIgeneralizedBloch.CRB_gradient_OCT`

returns the CRB

```
(CRBm0s, grad_ω1, grad_TRF) = MRIgeneralizedBloch.CRB_gradient_OCT(ω1, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, T2s, R2slT, grad_list, weights, isInversionPulse=isInversionPulse)
CRBm0s
```

`6341.950600296484`

along with the gradients:

```
p1 = plot(TR*(1:Npulses), grad_ω1 .* ((-1) .^ (1:Npulses)), ylabel="∂CRB(m0s) / ∂ω1 (s/rad)")
p2 = plot(TR*(1:Npulses), grad_TRF .* ((-1) .^ (1:Npulses)), ylabel="∂CRB(m0s) / ∂TRF (1/s)", xlabel="t (s)")
p = plot(p1, p2, layout=(2, 1), legend=:none)
```

Note that we remove the oscillating nature of the gradient for the display.

In this example, we limit the control to the following bounds

```
ω1_min = 0 # rad/s
ω1_max = 2e3π # rad/s
TRF_min = 100e-6 # s
TRF_max = 500e-6; # s
```

and the function `MRIgeneralizedBloch.bound_ω1_TRF!`

modifies `ω1`

and `TRF`

to comply with these bounds and returns a single vector in the range `[-Inf, Inf]`

that relates to the bounded control by a `tanh`

transformation:

`x0 = MRIgeneralizedBloch.bound_ω1_TRF!(ω1, TRF; ω1_min = ω1_min, ω1_max = ω1_max, TRF_min = TRF_min, TRF_max = TRF_max)`

```
400-element Vector{Float64}:
Inf
-1.683988287115545
-1.4729988693565415
-1.3209392186669002
-1.2011921567736858
-1.1019049796994895
-1.0167581958717071
-0.9419826140711725
-0.8751499388454517
-0.814605017748036
⋮
-1.1102230246251565e-16
-1.1102230246251565e-16
-1.1102230246251565e-16
-1.1102230246251565e-16
-1.1102230246251565e-16
-1.1102230246251565e-16
-1.1102230246251565e-16
-1.1102230246251565e-16
-1.1102230246251565e-16
```

Further, we initialize a gradient of the same length:

`G = similar(x0);`

and define the cost function:

```
function fg!(F, G, x)
ω1, TRF = MRIgeneralizedBloch.get_bounded_ω1_TRF(x; ω1_min = ω1_min, ω1_max = ω1_max, TRF_min = TRF_min, TRF_max = TRF_max)
(F, grad_ω1, grad_TRF) = MRIgeneralizedBloch.CRB_gradient_OCT(ω1, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, T2s, R2slT, grad_list, weights, isInversionPulse=isInversionPulse)
F = abs(F)
F += MRIgeneralizedBloch.second_order_α!(grad_ω1, grad_TRF, ω1, TRF; idx = .!isInversionPulse, λ=1e4)
F += MRIgeneralizedBloch.RF_power!(grad_ω1, grad_TRF, ω1, TRF; idx = .!isInversionPulse, λ=1e-3, Pmax=3e6, TR=TR)
F += MRIgeneralizedBloch.TRF_TV!(grad_TRF, ω1, TRF; idx = .!isInversionPulse, λ=1e3)
MRIgeneralizedBloch.apply_bounds_to_grad!(G, x, grad_ω1, grad_TRF; ω1_min = ω1_min, ω1_max = ω1_max, TRF_min = TRF_min, TRF_max = TRF_max)
return F
end;
```

We perform the optimization with the package Optim.jl, which requires the cost function `fg!(F, G, x)`

to take the cost, the gradient, and the control as input variables and to over-write the gradient in place. The cost function calculates the gradient of the CRB with above described optimal control code and we, further, add some regularization terms: `MRIgeneralizedBloch.second_order_α!`

penalizes the curvature of α, which smoothes the flip angle train and helps ensuring the hybrid state conditions. The penalty `MRIgeneralizedBloch.RF_power!`

penalizes the power deposition of the RF-pulse train if $\Sigma_i(ω_1^2[i] ⋅ T_{\text{RF}}[i]) / T_{\text{cycle}} ≥ P_{\max}$ and helps with compliance to safety limits. Assuming a reasonable `λ`

, the optimization will converge to an average RF-power deposition equal to or less than `Pmax`

in units of (rad/s)². Heuristically, the value `Pmax=3e6`

(rad/s)² proofed to be a reasonable choice for 3T systems. The penalty `MRIgeneralizedBloch.TRF_TV!`

penalizes fast fluctuations of $T_\text{RF}$. This penalty is justified by the knowledge that fluctuations of the control have negligible effect if they are fast compared to the biophysical time constants. We note, however, that this penalty is not required and rather ensure *beauty* of the result and speeds up convergence.

With all this in place, we can start the actual optimization

```
result = optimize(Optim.only_fg!(fg!), # cost function
x0, # initialization
BFGS(), # algorithm
Optim.Options(
iterations=10_000, # larger number as we use a time limit
time_limit=(15*60) # in seconds
)
)
```

```
* Status: success (objective increased between iterations)
* Candidate solution
Final objective value: 9.077900e+01
* Found with
Algorithm: BFGS
* Convergence measures
|x - x'| = 3.01e+02 ≰ 0.0e+00
|x - x'|/|x'| = 0.00e+00 ≤ 0.0e+00
|f(x) - f(x')| = 2.78e-10 ≰ 0.0e+00
|f(x) - f(x')|/|f(x')| = 3.07e-12 ≰ 0.0e+00
|g(x)| = 7.03e-09 ≤ 1.0e-08
* Work counters
Seconds run: 445 (vs limit 900)
Iterations: 1194
f(x) calls: 3447
∇f(x) calls: 3447
```

After transforming the optimized control back into the space of bounded $ω_1$ and $T_\text{RF}$ values

```
ω1, TRF = MRIgeneralizedBloch.get_bounded_ω1_TRF(result.minimizer; ω1_min = ω1_min, ω1_max = ω1_max, TRF_min = TRF_min, TRF_max = TRF_max)
α = ω1 .* TRF;
```

we analyze the CRB(m0s):

```
(CRBm0s, grad_ω1, grad_TRF) = MRIgeneralizedBloch.CRB_gradient_OCT(ω1, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, T2s, R2slT, grad_list, weights, isInversionPulse=isInversionPulse)
CRBm0s
```

`65.55456648506446`

and observe a substantial reduction. Further, we plot the optimized control:

```
p1 = plot(TR*(1:Npulses), α ./ π, ylabel="α/π")
p2 = plot(TR*(1:Npulses), 1e6TRF, ylim=(0, 1e3), xlabel="t (s)", ylabel="TRF (μs)")
p = plot(p1, p2, layout=(2, 1), legend=:none)
```

To further analyze the results, we can calculate and plot all magnetization components:

```
m = calculatesignal_linearapprox(α, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, T2s, R2slT; output=:realmagnetization)
m = vec(m)
xf = [m[i][1] for i=1:length(m)]
yf = [m[i][2] for i=1:length(m)]
zf = [m[i][3] for i=1:length(m)]
xs = [m[i][4] for i=1:length(m)]
zs = [m[i][5] for i=1:length(m)]
p = plot(xlabel="t (s)", ylabel="m (normalized)")
plot!(p, TR*(1:Npulses), xf ./(1-m0s), label="xᶠ")
plot!(p, TR*(1:Npulses), yf ./(1-m0s), label="yᶠ")
plot!(p, TR*(1:Npulses), zf ./(1-m0s), label="zᶠ")
plot!(p, TR*(1:Npulses), xs ./ m0s , label="xˢ")
plot!(p, TR*(1:Npulses), zs ./ m0s , label="zˢ")
```

And we can also plot the dynamics of the free spin pool on the Bloch sphere:

`p = plot(xf, zf, xlabel="xf", ylabel="zf", framestyle = :zerolines, legend=:none)`

As yᶠ is close to zero in this particular case, we neglect it in this 2D plot.

*This page was generated using Literate.jl.*