# Difference between revisions of "Optimal control module"

This section contains the description of the functions in the Optimal Control module of Spinach.

# Control operators

A bilinear control system should be split into a drift Hamiltonian (that which is beyond control of the experiment), and control Hamiltonians.

$\mathbf{H}=\mathbf{H}_0+\sum\limits_{k=1}^K c_k(t)\mathbf{H}_k$

where $\mathbf{H}_0$ is the drift Hamiltonian, and $\mathbf{H}_k$ are the control operators with time-dependent control amplitudes $c_k(t)$. The optimal control finds the optimal set of control amplitudes with respect to the overlap of an initial state and a target state.

An example of six control operators for an HCF spin system defined with Spinach:

    % 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;


The control operators should be defined as a cell array being a field of a structure optim

    optim.controls={LxH,LyH,LxC,LyC,LxF,LyF};


The structure optim should have optimisation and control field names added, to be an input of fminnewton.m, then to be parsed by optim_tols.m.

# Initial and target states

In the context of state-to-state transfer, an initial and target state must be defined. Example syntax for a three spin system stating with all magnetisation on the first spin with the target state of all magnetisation on the third spin:

    % Initial normalised state
rho=state(spin_system,{'Lz'},{1});
optim.rho=rho/norm(rho);

% Desired normalised target state
target=state(spin_system,{'Lz'},{3});
optim.target=target/norm(target);


The metric to measure the overlap of the current state from the optimised waveform with the target state is the fidelity (or the cost):

$J=\left\langle\texttt{target}\middle|\texttt{rho}\right\rangle$

It is useful to define a function to minimise, rather than maximise (as above). For this reason, the optimal control module uses a functional metric of overlap as the infidelity, defined as

$1-J=1-\left\langle\texttt{target}\middle|\texttt{rho}\right\rangle$

which would find a minimum at $1-J=0$

# Initial guess

The GRAPE optimal control relies on the piecewise constant approximation, where the total time allowed for the optimisation is split into nsteps time points. Although the Krotov optimal control does not rely on the piecewise constant approximation, time should also be split into (relatively small) time points, which numerically approximate continuous time.

An initial guess needs to be provided to any optimisation. In Spinach, the initial guess should be a matrix of size kxn, where k is the number of control operators, and n is the number of time points.

Syntax example for the above HCF spin system:

    % Power level
optim.power_level=2*pi*10e3;

% Pulse duration (seconds)
optim.pulse_duration=0.100;

% Number of time steps
optim.nsteps=50;

% Calculate the time step
optim.time_step=pulse_duration/nsteps;

% Take a random guess for the pulse
guess=randn(numel(optim.controls),optim.nsteps)/10;


Note that it is useful to scale the initial guess, in this case by dividing by 10.

# Objective function

The fundamental basis of the optimal control module is to define a cost function (also named an objective function), which returns the cost, gradient, and Hessian.

    function [cost,grad,hess] = cost_function(objective)


Spinach provides the GRAPE and Krotov methods of optimal control.

# Numerical optimisation

fminnewton.m contains an implementation of steepest descent, quasi‐Newton (SR1, BFGS and LBFGS), and Newton‐Raphson minimization algorithms. Works in much the same way as Matlab's built‐in fmincon, but allows for much greater flexibility, such as passing the initial Hessian guess. An example syntax could be:

    options = struct('method','newton','terminal_fidelity',1-eps,'max_iterations',Inf,'tol_x',eps,'tol_gfx',eps);
fminnewton(@cost_function,guess,options);


The options structure is very similar to Matlab's built‐in optimset and is created with Spinach's optim_tols.m.

## Example

The following is a simple code to optimise state transfer between to species in a chain of three spins (assuming spin_system is already defined):

    % Drift Hamiltonian
drift=hamiltonian(assume(spin_system,'nmr'));

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;

rho=state(spin_system,{'Lz'},{1});
target=state(spin_system,{'Lz'},{3});

optim.controls={LxH,LyH,LxC,LyC,LxF,LyF};     % Control operators
optim.rho=rho/norm(rho);                      % Initial normalised state
optim.target=target/norm(target);             % Desired normalised target state
optim.pulse_duration=0.100;                   % Pulse duration (seconds)
optim.nsteps=50;                              % Number of time steps
optim.time_step=pulse_duration/nsteps;        % Calculate the time step

% Take a random guess for the pulse
guess=randn(numel(optim.controls),optim.nsteps)/10;

optim.method=newton;
optim.terminal_fidelity=1-eps;
optim.max_iterations=Inf;
optim.tol_x=eps;
optim.tol_gfx=eps;

fminnewton(@cost_function,guess,optim);

optim.time_step,optim.nsteps,optim.rho,optim.target,optim.power_level);

penalty(waveform,'SNS',-ones(size(waveform)),ones(size(waveform)));

cost=cost+penalty_cost;