The B1 peak amplitude of the RF pulse and the specific absorption rate (SAR) can be reduced by a variation of the amplitude of the slice selective gradient (Gs) using the VERSE principle^1,2. Performing a simulation study, this idea was used to formulate an optimal control problem (gVERSE)³, but the results show a strong dependency on the chosen software and are restricted to a very coarse spatial resolution. Alternatively a trust region Newton method⁴ in an adjoint-based matrix-free variant allows for an unconstrained optimization of both RF pulse and gradient shape on a high spatial and temporal resolution⁵. However, the optimized results could exceed physical constraints and therefore may not be implementable on a real MR scanner. To overcome this limitations, we extend the framework to a trust-region semi-smooth Newton method to include physical constrains such as the peak B1 amplitude of the RF pulse and the slew rate of the gradient system.

The optimal control approach minimizes the discrepancy between the numerical Bloch simulation of the magnetizaion pattern $$$M(T;z)$$$ at refocusing time $$$T$$$ and each spatial position $$$z$$$ and the desired magnetization pattern $$$M_d(z)$$$ with three additional cost terms modeling the power of the RF pulse $$$B_{1,y}$$$, the slew rate $$$\dot{G}_s$$$ and the final amplitude of $$$G_s$$$, respectively:

$$\small\,J(M,u)=\frac{1}{2}\int_{-z}^z|M(T;z)-M_d(z)|_2^2\,dz+\frac\alpha2\int_0^{T}|B_{1,y}(t)|_2^2\,dt+\frac\beta2\int_0^{T}|\dot{G}_s(t)|_2^2\,dt+\frac\zeta2 |G_s(T)|_2^2.$$

Hard constraints on the peak B1 of the RF and the peak slew rate of the gradient are added using pointwise constraints: $$|B_{1}|_\infty\leq\,\!B_{1,\max}\qquad\text{and}\qquad|\dot{G}_s|_\infty\leq\dot{G}_{s,\max}.$$

Due to these inequalities, the first-order necessary conditions are no longer a system of smooth equations that can be solved by Newton's method. However, they can be reformulated as a non-smooth equation and solved by a generalized (semi-smooth) Newton method⁶, which is proven to converge superlinearly. This yields fast local convergence, whereas gradient-type methods only show linear convergence and are thus much slower to reach the optimum. For efficiency, the method is implemented in a matrix-free fashion, allowing a fine temporal resolution. Exact derivative information of the gradient and Hessian facilitate the efficient computation of optimal pulses. For convergence of the controls, the method is embedded into a trust-region framework⁴.

The described method was used to calculate and implement low SAR RF pulses with and without physical constraints (peak B1 of 13$$$\mu$$$T, peak slew rate of 175Tm-1s-1). The desired magnetization is defined on a spatial grid of 1001 points ($$$\pm$$$0.1m) and the shapes are computed with a temporal resolution of 0.01ms for a duration $$$T$$$ of 1.44ms. The initial guess for the controls are a standard Hamming windowed sinc with a time-bandwidth-product of 2 and a standard trapezoidal $$$G_s$$$-shape (shown in Figure 1 top left). To suppress excitation artefacts during refocusing, the RF pulse is fixed to zero after 1ms while the gradient shape is optimized until the final time $$$T$$$. To balance slice profile accuracy and SAR, we choose $$$\alpha=5\cdot 10^{-4}$$$ (SAR penalty), $$$\beta=3\cdot10^{-4}$$$ (slew rate penalty) and $$$\zeta=1.25$$$ (final $$$G_s$$$ amplitude penalty) for all computed results. The optimized shapes are implemented on a 3T MR scanner (Magnetom Skyra, Siemens Healthcare, Erlangen, Germany) using a modified GRE sequence (TE=3ms, TR=1000ms, FOV=200x200mm, matrix=320x320) to measure the slice profile of a water filled phantom.Figure 1 compares a conventional sinc-based RF pulse and a trapezoidal slice selective gradient with unconstrained and constrained joint optimized RF and gradient shapes. While the Bloch simulations after refocusing show a well-defined slice profile for all three examples, only the constrained shapes do not exceed the prescribed constraints and can be implemented without further modifications (e.g. stretching, scaling or clipping) on the above mentioned MR scanner. Figure 2 shows magnitude and phase of the reconstructed experimental 3T phantom measurements validating the Bloch simulations in Figure 1. Table 1 compares the B1 peak amplitude, peak slew rate and overall B1 power together with the root-mean-square error (RMSE) and mean absolute error (MAE) of the Bloch simulation compared to the ideal magnetization given in Figure 1. The optimized shapes fulfill the prescribed hard constraints while having the same low RMSE and MAE of the unconstrained solution. The presented approach demonstrates the joint optimization of RF and slice selective gradient shape with prescribed physical and technical constraints on a very fine grid. Compared to unconstrained methods, this guarantees practical implementability of optimized shapes without compromising the slice profile accuracy at reduced peak B1 and B1 power, especially for short RF pulses with a sharp slice profile. The computed pulses and shapes are therefore applicable in a wide range of imaging situations in MRI.