# RAPTOR tutorial: Control vector parametrization.

Here we explore how to parametrize the input (control) trajectory with a set of parameters. This is useful for optimization purposes, so you can parametrize a whole time-trajectory with just a few numbers.

## Contents

- Set up for time-varying H&CD actuators
- Define inputs via control vector parametrization:
- Generate input vector from parameters
- Set different time knots for each grid

The actuator trajectories are written as a sum of basis functions in time

where is the actuator index and is the basis function index. In the end, collect all together in a single vector.

We start by using use the same example as tutorial 3 but specify the inputs differently.

## Set up for time-varying H&CD actuators

This is identical to the set-up of tutorial 3

    close all hidden;
    clear; run('../RAPTOR_path.m');

    [config] = RAPTOR_config; % load default config
    config.grid.tgrid = [0:0.001:0.1]; % time grid
    config.grid.rhogrid = linspace(0,1,21); % spatial grid
    config.debug.iterplot = 0; %10; % plot profiles every N iterations
    config.debug.iterdisp = 5; % display progress every N iterations

    % heating & current drive parameters
    % set ECH/ECCD parameters
    config.echcd.params.active = true;
    config.echcd.params.rdep = [-1 -1]; % -1 to indicate it is commanded extrernally
    config.echcd.params.wdep = [0.3, 0.35]; % deposition width
    config.echcd.params.cd_eff = [0 1]; % CD efficiency
    config.echcd.params.uindices = [2 3 4 5]; % index in input vector
    % 2,3 are for the two powers, 4,5 is for the deposition location

    % set NBI parameters
    config.nbhcd.params.active = true;
    config.nbhcd.params.rdep = [0]; % -1 to indicate commanded extrernally
    config.nbhcd.params.wdep = [0.8]; % broad heating
    config.nbhcd.params.cd_eff = [0]; % no current drive
    config.nbhcd.params.uindices = [6]; % index in input vector

    % rerun the config file with the new config structure definition
    [model,params,init,g,v,U] = build_RAPTOR_model(config);

## Define inputs via control vector parametrization:

Define number of knots in time grid.

    params.cvp.tknots = 3; % 3 knots in time (start,middle,end)

Define a scaling factor for each input such that the numbers are "human".

    params.cvp.pscale = [1e6 1e6 1e6 1 1 1e6]; % [MA, MW, MW, [0,1], [0.1], MW] scales

Define the order for the polynomial bases of each vector.

    params.cvp.psporder = [1 0 0 1 1 0]; % piecew linear Ip, pw constant Paux, pwl rdep

## Generate input vector from parameters

specify parameter vector in order of appearance. Note the 'piecewise linear' (usporder=1) require 3 numbers, while 'piecewise constant' (usporder=0) require 2 numbers. In general = tknots+usporder-1

    p0 = [ [0.08,0.25,0.25], ... % Ip ramp to 250kA
        [0.0  0.5],... % Power step to 0.5MW
        [0.2  0.5],... % Power step to 0.5MW
        [0 0.4 0] ,... % Position sweep from 0 to 0.4 and back
        [0.2 0 0.2],...% Position sweep from 0.2 to 0 and back
        [0,1]]'; % 1MW NBI power.
    [U,dU_dp,dU_dt] = RAPTOR_cvp(p0,params);

    % inspect the results
    plot(params.tgrid,diag(1./params.cvp.pscale)*U);
    legend('Ip [MA]','P_{EC1} [MW]','P_{EC2} [MW]','r_{dep,EC1}','r_{dep,EC2}','P_{NBI} [MW]')

initial condition

    init.Ip0 = U(1,1);
    x0 = RAPTOR_initial_conditions(model,init,g(:,1),v(:,1));
    % Run RAPTOR
    simres = RAPTOR_predictive(x0,g,v,U,model,params);

         it  telaps newt     res   t[ms]  dt[ms]  Ip[kA] Icd[kA] Ibs[kA] Ioh[kA]      qe    qmin      q0   Vl[V] Te0[keV] Ti0[keV] ne0[e19]   f_ss
          1   0.077    4 2.9e-09       0       1      80   0.451    1.41    78.1    15.4    5.54    5.98 3.0e+00    0.40    0.40    1.00 1.7e+01
          6    0.22    2 5.9e-10       5       1      97     1.5     8.2    87.3    12.7    3.82    5.05 3.5e+00    0.79    0.79    1.00 1.5e+01
         11    0.34    2 3.0e-14      10       1     114    1.92    9.38     103    10.8    3.02    3.97 3.4e+00    0.98    0.98    1.00 1.5e+01
         16    0.46    2 1.5e-13      15       1     131    2.34    10.2     118    9.41     2.5    3.16 3.3e+00    1.17    1.17    1.00 1.5e+01
         21    0.58    2 1.1e-13      20       1     148    2.79    10.8     134    8.33    2.15     2.6 3.3e+00    1.36    1.36    1.00 1.4e+01
         26    0.69    2 7.1e-14      25       1     165    3.25    11.4     150    7.47    1.89    2.22 3.2e+00    1.56    1.56    1.00 1.4e+01
         31    0.81    2 5.1e-14      30       1     182    3.72    11.8     166    6.77    1.69    1.93 3.1e+00    1.76    1.76    1.00 1.4e+01
         36    0.94    2 5.6e-14      35       1     199    4.19    12.3     183    6.19    1.54    1.71 3.1e+00    1.96    1.96    1.00 1.4e+01
         41     1.1    2 2.6e-14      40       1     216    4.67    12.6     199    5.71    1.41    1.52 3.0e+00    2.15    2.15    1.00 1.4e+01
         46     1.2    2 4.0e-14      45       1     233    5.14    12.9     215    5.29    1.31    1.36 3.0e+00    2.35    2.35    1.00 1.4e+01
         51     1.3    5 3.9e-14      50       1     250    5.61    13.2     231    4.93    1.23    1.27 2.3e+00    6.36    6.36    1.00 1.2e+01
         56     1.5    3 4.3e-14      55       1     250    70.7    65.1     114    4.93     1.2    1.38 5.3e-01   13.70   13.70    1.00 3.3e+00
         61     1.6    2 4.5e-10      60       1     250      86    80.9      83    4.93    1.19    1.51 3.4e-01   15.92   15.92    1.00 2.2e+00
         66     1.8    2 5.3e-12      65       1     250    91.5    86.4      72    4.93    1.18    1.66 2.8e-01   16.86   16.86    1.00 1.9e+00
         71     1.9    2 1.7e-12      70       1     250    93.9    88.9    67.1    4.93    1.17    1.81 2.6e-01   17.45   17.45    1.00 1.7e+00
         76       2    2 9.1e-14      75       1     250    95.1    90.5    64.3    4.93    1.16    1.96 2.4e-01   17.93   17.93    1.00 1.7e+00
         81     2.1    2 5.0e-14      80       1     250    95.6    91.7    62.7    4.93    1.15     2.1 2.3e-01   18.38   18.38    1.00 1.6e+00
         86     2.3    2 5.9e-14      85       1     250    95.4    92.7    61.8    4.93    1.15     2.2 2.2e-01   18.80   18.80    1.00 1.6e+00
         91     2.4    2 7.4e-14      90       1     250    94.6    93.6    61.8    4.93    1.15    2.27 2.1e-01   19.20   19.20    1.00 1.5e+00
         96     2.5    2 4.8e-14      95       1     250    93.2    94.3    62.5    4.93    1.15    2.31 2.0e-01   19.55   19.55    1.00 1.5e+00
         it  telaps newt     res   t[ms]  dt[ms]  Ip[kA] Icd[kA] Ibs[kA] Ioh[kA]      qe    qmin      q0   Vl[V] Te0[keV] Ti0[keV] ne0[e19]   f_ss
        101     2.6    2 6.1e-14     100       1     250    91.3    94.9    63.8    4.93    1.14    2.32 1.9e-01   19.87   19.87    1.00 1.5e+00

plot results

    out = RAPTOR_out(simres,model,params);
    subplot(321); plot(out.time,out.U(1,:)); title('Ip');
    subplot(322); plot(out.time,out.U([2:3,6],:)); title('P_{aux}')
    subplot(323); plot(out.time,out.U([4:5],:)); title('\rho_{dep}')
    subplot(324); plot(out.time,out.f_ss); title('f_{ss} (steady stateness)')
    subplot(325); [cs,h] = contour(out.time,out.rho,out.q,[1 2 3 4],'color','k'); clabel(cs,h,'labelspacing',72); title('rational q locations');
    xlabel('t [s]'); ylabel('\rho');
    subplot(326); [cs,h] = contour(out.time,out.rho,out.te/1e3,[0.5 1 2 4 8],'color','k'); clabel(cs,h,'labelspacing',72); title('T_e contours [keV]')
    xlabel('t [s]'); ylabel('\rho');

## Set different time knots for each grid

    params.cvp.psporder = [1 1 1 1 1 1]; % all piecewise linear

    tknots1 = linspace(params.tgrid(1),params.tgrid(end),3)';
    tknots2 = linspace(params.tgrid(1),params.tgrid(end),4)';
    tknots3 = params.tgrid(1) + [0,0.45,0.5,1]'*params.tgrid(end);
    params.cvp.tknots = {tknots1,tknots2,tknots2,tknots2,tknots2,tknots3}; % define as cell array, one cell for each actuator

    p0 = [ [0.08,0.25,0.25], ... % Ip ramp to 250kA
        [0.0  0.2  0.2 0.5],... % Power step to 0.5MW
        [0.2  0.2  0.5 0.5],... % Power step to 0.5MW
        [0 0.4 0.4 0] ,... % Position sweep from 0 to 0.4 and back
        [0.2 0 0 0.2],...% Position sweep from 0.2 to 0 and back
        [0,0,1,1]]'; % 1MW NBI power.
    [U,dU_dp,dU_dt] = RAPTOR_cvp(p0,params);

    % *plot results*
    figure
    subplot(411); plot(params.tgrid,U(1,:));     title('Ip');
    subplot(412); plot(params.tgrid,U([2:3],:)); title('P_{aux}')
    subplot(413); plot(params.tgrid,U(6,:));     title('P_{aux}')
    subplot(414); plot(params.tgrid,U([4:5],:)); title('\rho_{dep}')
    xlabel('t[s]');

    % end of tutorial