# RAPTOR tutorial on trajectory optimization

This tutorial will explore more advanced features of RAPTOR and its use for trajectory optimization.

## Contents

- Simulate nominal ramp-up
- Prepare optimization problem
- Set up linear constraints and parameter bounds
- Cost function/constraint parameters
- Cost function definition
- Constraint definition
- Defining the SQP solver options
- Run SQP optimization algorithm
- Compare nominal and optimized trajectories.
- compare loop voltages

## Simulate nominal ramp-up

First define the environment, load the params structure and modify some default values to suit our needs.

    % Load and modify default parameters
    close all hidden;
    clear; run('../RAPTOR_path.m');

    [config] = RAPTOR_config;

    % modify some parameters
    config.grid.rhogrid = linspace(0,1,11); % spatial grid
    config.grid.tgrid = [0:0.001:0.1]; % time grid
    config.debug.iterplot = 0; % no plot, no disp this time (so we can see the optimizer output)
    config.debug.iterdisp = 0; %
    config.numerics.restol = 1e-8; % relax solution tolerance to be a bit faster

    % set ECH/ECCD parameters
    config.echcd.params.active = true;
    config.echcd.params.rdep = [0.4]; % off axis
    config.echcd.params.wdep = [0.4];
    config.echcd.params.cd_eff = [2]; % off-axis current drive
    config.echcd.params.uindices = [2]; % indices in input vectors

    config.nbhcd = nbhcd('none'); % no nbhcd

    % generate model, params, init, geometry g, kinetic profiles v
    [model,params,init,g,v,U] = build_RAPTOR_model(config);

Define a "nominal" ramp-up trajectory with constant ramp-rate and sudden switch-on of current drive actuators and simulate this nominal case.

    % Ip
    U0(1,:) = 250e3*ones(size(params.tgrid)); % input Ip trace: ramp from 80kA to 250kA
    U0(1,params.tgrid<=0.025) = linspace(80e3,250e3,sum(params.tgrid<=0.025));

    % Off-axis power
    U0(2,:) = 2e6*ones(size(params.tgrid)); % off-axis power
    U0(2,params.tgrid<0.05) = 1;

Define initial conditions and run simulation

    % Initial conditions
    init.Ip0 = U0(1,1);
    x0 = RAPTOR_initial_conditions(model,init,g(:,1),v(:,1));

    % Run simulation of nominal case
    [simres0] = RAPTOR_predictive(x0,g,v,U0,model,params);
    out0 = RAPTOR_out(simres0,model,params);

## Prepare optimization problem

We need to first define the discrete set of parameters which determine the actuator trajectories, via control vector parametrization.

    % tknots: 5 linearly spaced knot points.
    % must be a vector.
    tknots = linspace(0,1,5)'*params.tgrid(end);
    params.cvp.tknots = tknots;
    params.cvp.pscale = [1e6 1e6]; % MA, MW for human scale.
    params.cvp.psporder = [1 0]; % piecewise linear Ip, piecewise constant Paux
    params.cvp.uind = [1 1];
    params.cvp.transp_matrix = [];

    % initial condition for parameters
    p0 = [zeros(1,numel(tknots)),zeros(1,numel(tknots)-1)]';

    % generate basis functions in time for these settings.
    [~,dU_dp,~,nnp] = RAPTOR_cvp(p0,params);
    % nnp is the number of parameters for each input trajectory
    disp(nnp)
    % so U(1,:) is parametrized by 5 numbers, and U(2,;) by 4.

         5     4

    % Now fit parameters p to initial U0 by least squares i.e.
    % find p such that U0=dU_dp*p
    UUP = [squeeze(dU_dp(1,:,:)),squeeze(dU_dp(2,:,:))]'; % fit matrix
    p0 = UUP \ reshape(U0',numel(U0),1);
    p0 = max(p0,0); % (to avoid roundoff problems)

    % checking the result:
    p0
    % generate U with these p0;
    U(1,:) = squeeze(dU_dp(1,:,:))'*p0;
    U(2,:) = squeeze(dU_dp(2,:,:))'*p0;

    % check whether they match
    figure(1);
    clf;subplot(211); plot(params.tgrid,U(1,:),'b-',params.tgrid,U0(1,:),'r--'); %compare
        title('I_p'); xlabel('time [s]'); ylabel('[A]')
        subplot(212); plot(params.tgrid,U(2,:),'b-',params.tgrid,U0(2,:),'r--'); %compare
        title('P_{aux}'); xlabel('time [s]'); ylabel('W')

    p0 =

        0.0800
        0.2500
        0.2500
        0.2500
        0.2500
        0.0000
        0.0000
        2.0000
        2.0000

## Set up linear constraints and parameter bounds

    % flag to evolve forward sensitivity equation, required for optimization
    % problem.
    params.numerics.calc_sens = 1;
    % set constraints on parameter vector p
    % Linear equality constraints Aeq*p = beq
    % Linear inequality constraints Aineq*p <= bineq
    Aineq = []; bineq = []; Aeq = []; beq = [];

    % lower and upper bounds lbou <= p <= ubou (element-wise)
    lbou = zeros(numel(p0),1); ubou = zeros(numel(p0),1); % initialize
    lbou(1:nnp(1)) = 0.05; % 50kA minimum current
    lbou(nnp(1)+1:end) = 1e-6; % 0 power minimum
    ubou(1:nnp(1)) = 0.5; % 500kA maximum limit
    ubou(nnp(1)+1:end) = 2; % 2MW maximum Paux

Some parameters are fixed: ensure this by specifying equal upper and lower bounds. Could also do this by equality constraints directly.

    pfixmask = false(sum(nnp),1); % mask which is "true" for fixed parameters.
    pfixmask(1) = true; % fix init Ip
    pfixmask(nnp(1)) = true; % fix final Ip
    pfixmask(nnp(1)+1) = true; % fix init power
    pfixmask(sum(nnp)) = true; % fix final power

    % impose fixed elements of p
    lbou(pfixmask) = p0(pfixmask); % impose fixed p0 elements through bounds
    ubou(pfixmask) = p0(pfixmask); % impose fixed p0 elements through bounds

## Cost function/constraint parameters

set up optimization function parameters by calling the function `OL_cost_constraint` without input arguments

    OL_params = OL_cost_constraint;
    disp(OL_params)
    % This parameter structure contains "cost" and "constraint" fields with
    % separate parameters to choose cost function and nonlinear constraints.

                 cost: [1x1 struct]
           constraint: [1x1 struct]
         RAPTOR_model: []
        RAPTOR_params: []
                   x0: []
                   g0: []
                   v0: []
                check: 0

## Cost function definition

The cost function is configured through the `OL_params.cost` field The general form reads

where

- (proximity to target profile)
- Stationarity factor.
- Ohmic flux swing
- Distance from target loop voltage profile.

and each norm is weighted by a matrix .

    disp(OL_params.cost)

               nu_iota: 0
                nu_ss1: 1
                nu_ss2: 0
              nu_psioh: 0
                nu_upl: 0
              nu_shear: 0
           iota_target: 0
            upl_target: 0
                W_iota: 1
                 W_ss1: 1
                 W_ss2: 1
               W_shear: 1
                 W_upl: 1
        iterpointsplot: 0
               it_cost: Inf

In this case we only desire a stationary profile at the end of the simulation, i.e. flat loop voltage profile.

    % cost function parameters
    OL_params.cost.nu_iota  = 0;
    OL_params.cost.nu_ss1   = 1; % target zero final loop voltage (stationary state)
    OL_params.cost.nu_psioh = 0;
    OL_params.cost.nu_upl   = 0;

## Constraint definition

The implemented constraints are on the rotational transform profile or on the (edge) loop voltage. They are simply imposed through `OL_params.constraint`

    disp(OL_params.constraint);

              iotalimit: []
             Uedgelimit: []
                 Li_max: 1.2000
                 Li_min: 0.6500
           epsilon_iota: 1.0000e-06
          epsilon_Uedge: 1.0000e-06
        epsilon_Li3edge: 1.0000e-06
             constrplot: 0
             constrdisp: 0

We choose to only require q>1 at all times constraint function parameters

    OL_params.constraint.iotalimit = 1;
    OL_params.constraint.Uedgelimit = [];
    OL_params.constraint.constrdisp = 0;

We also pass the RAPTOR model, parameter, and initial conditions to the OL_params structure, so that the optimization routines know what RAPTOR model to simulate.

    OL_params.RAPTOR_model  = model;
    OL_params.RAPTOR_params = params;

    OL_params.x0 = x0;
    OL_params.g0 = g;
    OL_params.v0 = v;

    OL_params.check = 0;

## Defining the SQP solver options

As a final step, we need to configure the optimization solver as well

    % SQP solver options

    % A custom plot function handle
    h_plotfun = @(p,optimValues,state) OL_plot_function(p,optimValues,state,model,params,OL_params);

    % These are the plot functions to use during the optimization run.
    PlotFcns= {h_plotfun,@optimplotfval,@optimplotconstrviolation};

    % SQP solver options: pass the plot functions, and some other parameters.
    opt = optimset('PlotFcns',PlotFcns,...
                    'disp','iter',...
                    'GradObj','on',...
    		'GradConstr','on',...
                    'algorithm','sqp',...
                    'maxiter',30,...
                    'Tolfun',1e-4,...
                    'TolX',1e-4);

## Run SQP optimization algorithm

    % Define the cost function handle, i.e. calling OL_cost_constraint with
    % 'cost' as second argument.
    h_costfun = @(p) OL_cost_constraint(p,'cost',OL_params);

    % Define the constraint function handle, i.e. calling OL_cost_constraints
    % with 'constraint' as second argument.
    h_constraints = @(p) OL_cost_constraint(p,'constraint',OL_params);
    % h_constraints = []; % no constraints - uncomment this to have
    % unconstrained optimization.

    % Launch optimization problem
    disp('please wait, optimizing...');
    [poptimal,fval] = fmincon(h_costfun,p0,Aineq,bineq,Aeq,beq,lbou,ubou,h_constraints,opt);

    % compute trajectories corresponding to optimal p
    [U,dU_dp,ddU_dt_dp] = RAPTOR_cvp(poptimal,params);
    simres_opt = RAPTOR_predictive(x0,g,v,U,model,params);
    out_opt = RAPTOR_out(simres_opt,model,params);

    please wait, optimizing...
                                                              Norm of First-order
     Iter F-count            f(x) Feasibility  Steplength        step  optimality
        0       1    2.118623e-01   2.662e-05                           4.142e+00
        1      13    2.090542e-01   2.609e-05   1.977e-02   7.227e-03   2.559e+00
        2      15    1.913071e-01   5.860e-04   7.000e-01   1.972e-01   2.385e+00
        3      16    1.237641e-01   1.116e-04   1.000e+00   1.068e-01   2.682e-01
        4      17    1.284364e-01   2.675e-05   1.000e+00   8.713e-02   1.107e+00
        5      18    1.299853e-01   8.142e-06   1.000e+00   8.512e-02   4.427e-01
        6      19    1.281024e-01   2.365e-06   1.000e+00   2.720e-01   2.438e-01
        7      20    1.295596e-01   5.770e-07   1.000e+00   3.956e-02   1.757e-01
        8      21    1.291161e-01   7.156e-08   1.000e+00   5.107e-02   1.368e-01
        9      22    1.232456e-01   9.290e-09   1.000e+00   2.582e-01   1.654e-01
       10      23    9.579152e-02   1.446e-05   1.000e+00   1.034e+00   9.405e-01
       11      25    1.012049e-01   6.545e-06   7.000e-01   3.234e-01   9.041e-01
       12      29    8.796228e-02   5.868e-06   3.430e-01   3.001e-01   7.654e-01
       13      30    8.807077e-02   2.761e-06   1.000e+00   7.193e-02   4.101e-01
       14      31    8.906067e-02   1.254e-06   1.000e+00   4.188e-02   4.298e-01
       15      32    9.196925e-02   2.353e-07   1.000e+00   7.056e-02   1.028e-01
       16      33    9.286400e-02   1.744e-08   1.000e+00   2.182e-02   3.607e-02
       17      34    9.293442e-02   7.429e-10   1.000e+00   2.540e-03   1.033e-02
       18      35    9.293617e-02   1.111e-10   1.000e+00   1.927e-04   5.011e-03
       19      36    9.293617e-02   1.111e-10   7.000e-01   1.042e-04   5.011e-03

    Local minimum possible. Constraints satisfied.

    fmincon stopped because the size of the current step is less than
    the selected value of the step size tolerance and constraints are
    satisfied to within the default value of the constraint tolerance.

## Compare nominal and optimized trajectories.

    disp('optimized case')
    RAPTOR_plot_GUI({out0,out_opt});

    optimized case

## compare loop voltages

    clf;
    subplot(221);
    hp(1) = plot(out0.time,out0.U(1,:)/1e3,'k'); hold on;
    hp(2) = plot(out_opt.time,out_opt.U(1,:)/1e3,'r');
    legend(hp,{'nominal','optimal'},'location','southeast');
    title('Ip [kA]'); xlabel('time [s]')

    subplot(222);
    plot(out0.time,out0.U(2,:)/1e6,'k');hold on;
    plot(out_opt.time,out_opt.U(2,:)/1e6,'r');
    title('P_{ECCD} [MW]'); xlabel('time [s]')

    subplot(223);
    plot(out0.rho,out0.upl(:,end),'k');hold on;
    plot(out_opt.rho,out_opt.upl(:,end),'r');
    ylabel('[V]');xlabel('\rho')
    title('Final loop voltage profile')

    subplot(224);
    plot(out0.rho,out0.q(:,end),'k');hold on;
    plot(out_opt.rho,out_opt.q(:,end),'r');
    ylabel('q');xlabel('\rho')
    title('Final q profile')