Optimal control module
This section contains instructions on setting up optimal control calculations in Spinach. The package contains an implementation of the gradient ascent pulse engineering (GRAPE) method, as well as its LBFGS-GRAPE and Newton-GRAPE enhancements. Ensemble optimal control is supported. Multiple examples are available in examples/optimal_control directory of the standard distribution.
Contents
Stage 1: specify the spin system
The usual spin system specification should be done. For example, if we are interested in performing magnetisation transfer from one spin to another in a peptide chain in a liquid, the NMR calculation for that chain needs to be set up:
% Magnetic field sys.magnet=9.4; % Spin system sys.isotopes={'15N','1H','13C','13C','13C','15N'}; % Chemical shifts, ppm inter.zeeman.scalar={0.0, 0.0, 55.0, 30.0, 175.0, 0.0}; % Scalar couplings, Hz (literature values) inter.coupling.scalar=cell(6); inter.coupling.scalar{1,3}=-11; inter.coupling.scalar{2,3}=140; inter.coupling.scalar{3,4}=35; inter.coupling.scalar{3,5}=55; inter.coupling.scalar{3,6}=7; inter.coupling.scalar{5,6}=-15; % Basis set bas.formalism='sphten-liouv'; bas.approximation='IK-0'; bas.level=4; % Spinach housekeeping spin_system=create(sys,inter); spin_system=basis(spin_system,bas);
Stage 2: get initial and target states
You need to request your initial state and your desired target states from Spinach. The states must then be normalised. For example, if we are looking to move all magnetisation from Lz on spin 2 to Lz on spin 5, the following needs to be written:
% Set up and normalise the initial state rho_init=state(spin_system,{'Lz'},{2}); rho_init=rho_init/norm(full(rho_init),2);
% Set up and normalise the target state rho_targ=state(spin_system,{'Lz'},{5}); rho_targ=rho_targ/norm(full(rho_targ),2);
The first line of each block calls the usual state.m function, the the second line divides the resulting state vector by its 2-norm. If you have multiple source states that must be transferred to multiple destination states, put the corresponding vectors into cell arrays.
Stage 3: get drift and control operators
In an optimal control problem, the Hamiltonian is split into the part that the instrument cannot change (called drift Hamiltonian) and the operators whose coefficients the instrument can control (called control operators):
[math]\mathbf{H}=\mathbf{H}_0+\sum\limits_{k=1}^K c_k(t)\mathbf{H}_k[/math]
where [math]\mathbf{H}_0[/math] is the drift
Hamiltonian, and [math]\mathbf{H}_k[/math] are the control operators with time-dependent control coefficients [math]c_k(t)[/math]. In a typical NMR experiment, the drift Hamiltonian contains chemical shifts and J-couplings, and control operators are Lx and Ly spin operators. All of these should be requested from Spinach as usual: the Hamiltonian is obtained from hamiltonian.m (if relaxaton is present, you should also call relaxation.m), and spin operators are obtained by calling operator.m function. In our peptide system, we would like Lx and Ly operators on carbon, nitrogen and protons, therefore:
% Drift Hamiltonian drift=hamiltonian(assume(spin_system,'nmr')); % Control operators LpH=operator(spin_system,'L+','1H'); LpC=operator(spin_system,'L+','13C'); LpF=operator(spin_system,'L+','19F'); LxH=(LpH+LpH')/2; LyH=(LpH-LpH')/2i; LxC=(LpC+LpC')/2; LyC=(LpC-LpC')/2i; LxF=(LpF+LpF')/2; LyF=(LpF-LpF')/2i;
If you have an ensemble with a different drift Hamiltonian in each system, put those Hamiltonians into a cell array.
Stage 4: put together the control structure
Spinach requires the user to specify all optimal control related options in a single structure called control. The following fields are available:
Optimal control parameter structure | ||
---|---|---|
Field name | Mandatory/Optional | Content |
control.fidelity | O | A character string specifying the type of fidelity functional. Specify 'real' for the real part of the overlap between the final state and the target state, 'imag' for the imaginary part, and 'square' for the absolute square of the overlap. The default is 'real'. |
control.operators | M | A cell array of matrices of control operators. If phase-amplitude optimisation or plotting is intended, these must be in the order {L1x,L1y,L2x,L2y,...}. |
control.drifts | M | A cell array containing drift Hamiltonian matrices. If a single system is treated, there should be only one matrix; if there is an ensemble (with, for example, different J-couplings in each member), the cell array should contain multiple matrices. |
control.rho_init | M | A cell array of initial states. For example, if state vector A should be moved into state vector X, and state vector B into state vector Y, your specification for the initial states should be {A,B}. |
control.rho_targ | M | A cell array of target states. For example, if state vector A should be moved into state vector X, and state vector B into state vector Y, your specification for the target states should be {X,Y}. |
control.pulse_dur | M | Control pulse duration in seconds. |
control.pulse_nsteps | M | Number of time slices in the control pulse, an even integer. |
control.pwr_levels | M | A vector of control power levels across the ensemble, in rad/s. These correspond to things like non-uniform B1 field: the amplifier outputs some nominal power, but it is felt differently in different parts of the ensemble. |
control.offsets | O | A vector of generalised offsets across the ensemble, in Hz. These correspond to things like resonance offsets that may be different for different signals in the spectrum. |
control.offset_operator | O | Offset operator, usually something like Lz on the spins of interest. For each member of the ensemble, this operator will be multiplied by the corresponding element of control.offsets array, converted to rad/s and added to the drift Hamiltonian of that ensemble member. |
control.basis | O | It is occasionally advantageous to run waveform optimisations not point-by-point, but with respect to coefficients in front of of some smooth functions, for example sinusoids. Those functions should be specified as rows of this array, which should have dimensions [nfunctions x control.pulse_npoints]. The functions should be orthonormal. |
control.penalties | O | A cell array of character strings specifying the penalty functionals to be added to the fidelity functional. These are needed to make sure the waveform does not fly out into infinity. See penalty.m for details. The allowed values are 'NS' (norm square), 'SNS' (spillout norm square) and 'DNS' (derivative norm square). Default is 'SNS'. |
control.p_weights | O | A vector of weight multipliers for the penalty functinals specified in control.penalties - default is 'SNS' penalty with weight 100. |
control.u_bound | O | Upper bound for the 'SNS' penalty, as a fraction of the power transmitted by the amplifier. Default value is +1. |
control.l_bound | O | Lower bound for the 'SNS' penalty, as a fraction of the power transmitted by the amplifier. Default value is -1. |
control.plotting | O | A cell array of character strings specifying the plotting options for the diagnostic output. Possible values are: 'xy_controls' (coefficient of each control as a function of time), 'phi_controls' (phase of each X,Y control pair as a function of time),'amp_controls' (amplitude of each X,Y control pair as a function of time),'correlation_order' (correlation order dynamics computed by trajan.m),'coherence_order' (coherence order dynamics computed by trajan.m,'local_each_spin' (probability density at each spin computed by trajan.m,'robustness' (fidelity histogram for the ensemble). |
Note that there are four independent ensembles: the first one for source-target pairs, the second one for drifts, the third one for control powers and the fourth one for offsets. These ensembles are kroneckered together. So, if you have two state-target pairs, ten different drift Hamiltonians, five power levels and five offsets, your ensemble would have 2x10x5x5=500 systems in it!
For the spin system specified above, the control structure would looke something like:
% Define control parameters control.drifts={H}; % Drift control.operators={LxH,LyH,LxC,LyC,LxN,LyN}; % Controls control.rho_init={rho_init}; % Starting state control.rho_targ={rho_targ}; % Destination state control.pwr_levels=2*pi*10e3; % Pulse power control.pulse_dur=20e-3; % Pulse duration control.pulse_nsteps=200; % Time points control.penalties={'SNS'}; % Penalty control.p_weights=100; % Penalty weight % Control trajectory analysis plots control.plotting={'correlation_order','local_each_spin',... 'xy_controls'};
Finally, call the optimcon.m function
% Spinach housekeeping spin_system=optimcon(spin_system,control);
to update the spin_system data structure.
Stage 5: specify the initial guess
Random numbers are usually fine. If the optimisation is carried out point-by-point (i.e. control.basis option is not specified), the dimension of the initial guess should be [n_controls x pulse_nsteps], for example:
% Take a random guess for the waveform guess=randn(numel(control.operators),control.pulse_nsteps)/3;
If the basis is specified, the dimension of the initial guess should be [n_controls x n_functions]. The initial guess is specified as a fraction of the amplifier power output, i.e. it should normally be within [-1,+1] range.
When an optimisation is to be carried out with respect to just the pae of the control sequence, the number of rows in the guess waveform should be equal to half the number of control operators, e.g.
% Take a random guess for the phases guess=randn(numel(control.operators)/2,control.pulse_nsteps)/3;
Stage 6: run the optimization
The following functions return the fidelity, its gradient and its Hessian:
- grape_xy.m
- Optimisation of waveforms expressed as Lx, Ly coefficients.
- grape_phase.m
- Optimisation of waveforms expressed in amplitude-phase coordinates with respect to phase.
In principle, the derivatives these functions return can be supplied to any optimisation function, including those supplied with Matlab. Spinach does, however, include optimisation tools that are designed specifically for optimal control settings:
- fminnewton.m
- Gradient descent, quasi-Newton and Newton-Raphson methods.
- fminsimplex.m
- Nelder-Mead simplex and multidirectional search.
For the purposes of optimal control optimisations, the following sets of settings provide a quick way to get started. The best balance between calculation speed and convergence rate is struck by LBFGS-GRAPE. To run Cartesian coordinate optimisation, use:
% Optimisation parameters optim.method='lbfgs'; % Optimisation method optim.extremum='maximum'; % Extremum type % Run the optimization fminnewton(@grape_xy,guess,optim,spin_system);
To run phase optimisation, declare an amplitude profile (same size as your phase guess) and call grape_phase.m instead:
% Optimisation parameters optim.method='lbfgs'; % Optimisation method optim.extremum='maximum'; % Extremum type % Use constant amplitude profile amp_profile=ones(size(guess)); % Run the optimization fminnewton(@grape_phase,guess,optim,amp_profile,spin_system);
To use the Newton-Raphson method, specify 'newton' instead of 'lbfgs'. Using gradient ascent ('grad_ascent') is almost never a good idea. The long list of possible optimization parameters and settings is provided in the documentation to fminnewton.m function.
See also
Control functions
grape_phase.m, grape_xy.m, ensemble.m, optimcon.m, grape.m, penalty.m
Optimisation functions
fminnewton.m, fminsimplex.m, hessprep.m, hessreg.m, lbfgs.m, linesearch.m, objeval.m, optim_report.m, optim_tols.m, quasinewton.m
Other relevant functions
hess_reorder.m, dirdiff.m, step.m
Version 2.0, authors: Ilya Kuprov, David Goodwin