# RAPTOR tutorial : Local linearization

In this tutorial we first simulate the time-evolution of the plasma ramp-up using a RAPTOR predictive run. From this, we extract a local linear model describing the (linear) response of the state to a small input perturbation, and compare this response to the nonlinear response of the original model.

## Contents

- Nominal (non-perturbed) ramp-up simulation
- Generate discrete-time, linear-time-varying state-space model
- Check response of lin vs. nonlin system: perturb ECCD power

## Nominal (non-perturbed) ramp-up simulation

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

    [config] = RAPTOR_config; % load default config

    config.echcd.params.active = true;
    config.echcd.params.rdep = [0.4]; % one ECCD source at rho=0.4;
    config.echcd.params.wdep = [0.35];  % wdep =0.35
    config.echcd.params.cd_eff = [2];  % double current drive efficiency to see effect better.
    config.echcd.params.uindices = [2]; % index of power for each actuator in input vector

    config.nbhcd.params.uindices = 3; % for correct indexing of inputs

    % Other modifications
    config.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 = 10; % display progress every N iterations

    % Create RAPTOR model with new parameters
    [model,params,init,g,v,U] = build_RAPTOR_model(config);

    % inputs
    % plasma current
    U(1,:)= 200e3*ones(size(params.tgrid)); % input Ip trace: ramp from 80kA to 200kA
    U(1,params.tgrid<0.02) = linspace(80e3,200e3,sum(params.tgrid<0.02));
    % EC actuator power: off-axis current drive, 1MW start at 40ms
    U(2,:) = zeros(size(params.tgrid)); U(2,params.tgrid>0.04) = 1e6;

    % initial conditions
    % Plasma current
    init.Ip0 = U(1,1);
    % Initial condition for state
    x0 = RAPTOR_initial_conditions(model,init,g(:,1),v(:,1));
    % run RAPTOR
    simres = RAPTOR_predictive(x0,g,v,U,model,params);
    % get outputs
    out = RAPTOR_out(simres,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.063    4 5.1e-10       0       1      80       0    1.41    78.6    15.4    5.59    5.59 3.4e+00    0.21    0.21    1.00 1.2e+01
         11    0.32    2 2.6e-10      10       1     143       0    4.06     139    8.61    3.19    3.19 5.3e+00    0.40    0.40    1.00 2.0e+01
         21    0.56    3 9.0e-13      20       1     200       0    6.24     194    6.16    2.32    2.32 4.0e+00    0.64    0.64    1.00 1.6e+01
         31    0.79    2 1.8e-11      30       1     200       0    6.75     193    6.16    1.79    1.79 2.7e+00    0.83    0.83    1.00 1.1e+01
         41       1    2 9.2e-13      40       1     200       0    6.67     193    6.16    1.38    1.38 2.4e+00    0.95    0.95    1.00 8.4e+00
         51     1.3    2 2.7e-10      50       1     200    72.3    25.2     102    6.16     1.3     1.3 5.0e-01    3.20    3.20    1.00 3.2e+00

## Generate discrete-time, linear-time-varying state-space model

To generate a linear-time varying state space model, we linearize the nonlinear state equation around the nominal trajectory

or

with

and

    nt = numel(params.tgrid);
    Ad = cell(1,nt); Bd = Ad; % allocate empty cells
    for ik = 1:nt
        % apply scalings for numerical conditioning of inverse (same scalings
        % as in RAPTOR_predictive_step.m
       iA = out.df_dxk1{ik};
       Ad{ik} = -iA \ out.df_dxk{ik};
       Bd{ik} = -iA \ out.df_du{ik};
    end
    Cd = eye(size(Ad{1})); % all states as outputs
    Dd = 0; % no direct feedthrough

## Check response of lin vs. nonlin system: perturb ECCD power

    % Perturb input
    Upert = U; % copy old U
    Upert(1,:) = U(1,:)+1e4*randn(size(U(2,:))); % perturb Ip with noise
    Upert(2,:) = U(2,:)+1e5*rand(size(U(2,:)));  % perturb P_EC with noise

    % Plot perturbed and nominal inputs
    clf; subplot(211)
    plot(params.tgrid,U(1,:)/1e3,'b',params.tgrid,Upert(1,:)/1e3,'k--');
    title('I_p'); ylabel('[kA]'); xlabel('time [s]')

    subplot(212);
    plot(params.tgrid,U(2,:)/1e6,'b',params.tgrid,Upert(2,:)/1e6,'k--');
    title('ECCD at \rho=0.4'); ylabel('[MW]'); xlabel('time [s]')
    legend('nominal','perturbed','location','southeast');

Perturbed Nonlinear Simulation

    % new nonlinear simulation
    simrespert = RAPTOR_predictive(x0,g,v,Upert,model,params);
    outpert=RAPTOR_out(simrespert,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.061    4 1.3e-09       0       1    80.4   0.199    1.41    78.8    15.3     5.7     5.7 3.3e+00    0.22    0.22    1.00 1.2e+01
         11    0.38    4 3.0e-11      10       1     141   0.311    4.97     136    8.74    3.38    3.38 1.8e+00    0.46    0.46    1.00 2.9e+01
         21    0.75    4 2.5e-11      20       1     185   0.768    7.03     178    6.64    2.48    2.48 8.6e-01    0.70    0.70    1.00 3.3e+01
         31     1.1    3 3.7e-14      30       1     213    1.22    7.55     204    5.78     1.9     1.9 4.3e+00    0.91    0.91    1.00 1.9e+01
         41     1.4    3 1.8e-11      40       1     197    0.24    7.05     189    6.27    1.46    1.46 2.0e+00    0.99    0.99    1.00 8.7e+00
         51     1.7    3 2.2e-11      50       1     200    75.6    26.7    97.4    6.17    1.38    1.38 1.2e+00    3.37    3.37    1.00 1.0e+01

Linear simulation

This takes the perturbation and returns the state perturbation

    % linear model
    dU = Upert-U; % delta U
    dX = zeros(size(Ad{1},2),nt); % allocate
    for ik = 1:(nt-1) % run simulation
        dX(:,ik+1) = Ad{ik}*dX(:,ik) + Bd{ik}*dU(:,ik);
    end
    Xpertlin = Cd*dX + simres.X; % linear prediction of state (X+\delta X)

    % use internal RAPTOR evaluation functions to get the correponding te and iota
    iotalin = eval_iota(Xpertlin,simres.G,simres.V,model,false);
    telin   = eval_te(Xpertlin,simres.G,simres.V,model,false);

Plots comparing linear and nonlinear models

    clf; set(gcf,'position',[50 50 700 600])
    subplot(211)
    contour(out.time,out.rho,out.te/1e3,[0.2,0.5 1,2,3,4],'b'); hold on;
    contour(outpert.time,outpert.rho,outpert.te/1e3,[0.2,0.5,1,2,3,4],'k','linewidth',2);
    hold on;
    [cs,h] = contour(out.time,out.rho,telin/1e3,[0.2,0.5,1,2,3,4],'r:','linewidth',2);
    clabel(cs,h);
    title('T_e profile contours'); xlabel('t'); ylabel('\rho');

    hleg=legend('nominal','perturbed nonlinear','peturbed linear');
    set(hleg,'location','Southeast','orient','hor')

    subplot(212)
    contour(out.time,out.rho,1./out.iota,[1,6/5,4/3,3/2,2,3,4],'b'); hold on;
    contour(outpert.time,outpert.rho,1./outpert.iota,[1,6/5,4/3,3/2,2,3,4],'k','linewidth',2);
    [cs,h] = contour(out.time,out.rho,1./iotalin,[1,6/5,4/3,3/2,2,3,4],'r-','linewidth',1);
    clabel(cs,h);
    title('q profile contours'); xlabel('t'); ylabel('\rho');

    % end of tutorial