# Difference between revisions of "Optimal control module"

(Added full example) |
(small updates to the optimal control toolbox - still needs more updating) |
||

Line 4: | Line 4: | ||

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

− | <math>H= | + | <center> |

+ | <math>\mathbf{H}=\mathbf{H}_0+\sum\limits_{k=1}^K c_k(t)\mathbf{H}_k</math> | ||

+ | </center> | ||

− | where <math> | + | where <math>\mathbf{H}_0</math> is the <tt>drift</tt> Hamiltonian, and <math>\mathbf{H}_k</math> are the control operators with time-dependent control amplitudes <math>c_k(t)</math>. 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 specification|spin system]] defined with Spinach: | An example of six control operators for an HCF [[Spin system specification|spin system]] defined with Spinach: | ||

Line 78: | Line 80: | ||

Spinach provides the GRAPE and Krotov methods of optimal control. | Spinach provides the GRAPE and Krotov methods of optimal control. | ||

− | |||

− | |||

− | |||

− | |||

− | |||

− | |||

− | |||

− | |||

− | |||

− | |||

− | |||

− | |||

− | |||

− | |||

− | |||

− | |||

− | |||

− | |||

− | |||

− | |||

− | |||

− | |||

− | |||

− | |||

− | |||

=Numerical optimisation= | =Numerical optimisation= | ||

Line 164: | Line 141: | ||

Note: this example uses the optim structure to collect all options into fminnewton - it is equivalent to using matlab's struct function. The cost_function is supplied as a nested function in this example, although a local function or simple function could also be used. | Note: this example uses the optim structure to collect all options into fminnewton - it is equivalent to using matlab's struct function. The cost_function is supplied as a nested function in this example, although a local function or simple function could also be used. | ||

+ | |||

+ | ==Optimal control toolbox function== | ||

+ | |||

+ | =Control functions= | ||

+ | |||

+ | [[control_sys.m]] | ||

+ | |||

+ | [[grape.m]] | ||

+ | |||

+ | [[krotov.m]] | ||

+ | |||

+ | [[penalty.m]] | ||

+ | |||

+ | =Optimisation functions= | ||

+ | |||

+ | [[fminkrotov.m]] | ||

+ | |||

+ | [[fminnewton.m]] | ||

+ | |||

+ | [[fminsimplex.m]] | ||

+ | |||

+ | [[hessprep.m]] | ||

+ | |||

+ | [[hessreg.m]] | ||

+ | |||

+ | [[lbfgs.m]] | ||

+ | |||

+ | [[linesearch.m]] | ||

+ | |||

+ | [[optim_report.m]] | ||

+ | |||

+ | [[quasinewton.m]] |

## Revision as of 19:02, 25 August 2016

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

## Contents

# Control operators

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

[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 amplitudes [math]c_k(t)[/math]. 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):

[math]J=\left\langle\texttt{target}\middle|\texttt{rho}\right\rangle[/math]

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

[math]1-J=1-\left\langle\texttt{target}\middle|\texttt{rho}\right\rangle[/math]

which would find a minimum at [math]1-J=0[/math]

# 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.power_level=2*pi*10e3; % Power level (rad/s) 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); function [cost,grad,hess] = cost_function(objective) [data,cost,grad,hess]=grape(spin_system,drift,optim.controls,waveform,... optim.time_step,optim.nsteps,optim.rho,optim.target,optim.power_level); [penalty_cost,penalty_grad,penalty_hess]=... penalty(waveform,'SNS',-ones(size(waveform)),ones(size(waveform))); cost=cost+penalty_cost; grad=grad+penalty_grad; hess=hess+penalty_hess; end

Note: this example uses the optim structure to collect all options into fminnewton - it is equivalent to using matlab's struct function. The cost_function is supplied as a nested function in this example, although a local function or simple function could also be used.