Difference between revisions of "Optimal control module"

From Spinach Documentation Wiki
Jump to: navigation, search
(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=H_0+\sum\limits_k c_kH_k(t)</math>
+
<center>
 +
<math>\mathbf{H}=\mathbf{H}_0+\sum\limits_{k=1}^K c_k(t)\mathbf{H}_k</math>
 +
</center>
  
where <math>H_0</math> is the <tt>drift</tt> Hamiltonian, and <math>H_k(t)</math> are the control Hamiltonians with control amplitudes <math>c_k</math>. The optimal control finds the optimal set of control amplitudes with respect to the overlap of an initial state and a target state. The controls are passed to either <tt>[[grape.m]]</tt> or <tt>[[krotov.m]]</tt> within a cell array structure. The control amplitudes are passed in the form a a multi-row vector <tt>waveform</tt>, rows corresponding to control channels and columns corresponding to time increments.
+
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.
 
==GRAPE==
 
Gradient Ascent Pulse Engineering (GRAPE) objective function and gradient. Propagates the system through a user‐supplied shaped pulse from a given initial state and projects the result onto the given final state. The real part of the projection is returned, along with its gradient with respect to amplitudes of all operators in every step of the shaped pulse.
 
 
    [diag_data,total_objective,total_grad,total_hess]=...
 
          [[grape]](spin_system,drift,controls,waveform,time_step,nsteps,...
 
          initial_state_array,target_state_array,power_level)
 
 
Four derivative calculation methods are available: Hausdorff series, second order central finite difference, fourth order central finite difference and Sophie Schirmer's expm algorithm. The default (Sophie's algorithm) is fast and accurate to machine precision. Sophie Schirmer’s expm algorithm is the only option for the Newton‐Raphson method.
 
 
==Krotov==
 
Krotov and Zhu‐Rabitz algorithms for phase‐sensitive state‐to‐state transfer problems. Tannor version as described in D. Tannor, V. Kazakov, V. Orlov, "Time‐dependent Quantum Molecular Dynamics", Plenum, 1992; Zhu‐Rabitz version as described in W. Zhu, H. Rabitz, J. Chem. Phys. 1998, 109, 385.
 
 
    [waveform,diag_data]=[[krotov]](spin_system,drift,controls,waveform,time_step,...
 
          nsteps,initial_state,target_state,lambda,power_level,options);
 
 
The function returns the optimized waveform and a similar diagnostics data structure as grape.m.
 
 
==Penalty terms==
 
penalty terms for pulse waveforms. Returns the penalty function and its gradient for the waveform, which should be supplied as row vector or a horizontal stack thereof.
 
 
    [penalty_term,penalty_grad]=...
 
          [[penalty]](waveform,type,floor_bound,ceiling_bound);
 
 
If a stack of waveforms is supplied, individual rows of the stack are treated as independent and the penalties are summed together.
 
  
 
=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.

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.

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