# RAPTOR tutorial: Feedback controllers

In this tutorial we define feedback controllers for some plasma quantities. This means some inputs are not known a-priori, but are determined on-line by the controllers. Instead, the inputs represent references to the system.

In this example we will control the central heating power to obtain a reference . The controller implemented in `beta_controller.m` for this purpose is a standard PI controller. Inspect this file for more details, and use it as a template for your own controller.

Since some other actuators (for example Ip) are NOT controlled in feedback, we also need to define an "open-loop" controller, that will pass the inputs directly to the actuator command.

## Contents

- Define controllers
- Plot results

Load and modify default parameters

    close all hidden;
    clear; run ../RAPTOR_path.m;

    [config] = RAPTOR_config;
    config.grid.tgrid = [0:0.001:0.1];
    config.grid.rhogrid = linspace(0,1,21);
    config.debug.iterdisp = 10; % plot every 10;
    % actuators
    config.echcd = echcd('echcd_gaussian'); % use default with 2 EC sources
    config.echcd.params.active = true;
    config.nbhcd = nbhcd('none'); % no NBH

recall that this gives two EC sources, one central ECH and one off-axis ECCD.

    disp(config.echcd.params)

           active: 1
             rdep: [0 0.4000]
             wdep: [0.3000 0.3500]
           cd_eff: [0 1]
        eta_eccd0: 7.0000e+14
          cd_trap: 0
         uindices: [2 3]
            check: 0

note that the EC actuator indices are [2,3] respectively (1=Ip by default).

## Define controllers

We need two controllers: Controllers are defined very similarly to heating&current drive actuators. First, call RAPTOR_contollers with controller function handles as inputs:

    config.controllers = RAPTOR_controllers({@openloop_controller,@beta_controller});

This loads the default parameters for two controllers. We now modify the controller parameters.

    config.controllers.beta_controller.params.Ki = 2e9; % i gain
    config.controllers.beta_controller.params.Kp = 2e6; % p gain
    config.controllers.beta_controller.params.t_fb = 0.05;  % switch-on time of feedback controller
    config.controllers.beta_controller.params.iterplot = 10;
    config.controllers.beta_controller.params.fname = ... % file to store internal controller data.
        fullfile(config.envopts.RAPTOR_path,'personal','beta_controller_data');

Now define mapping between various controllers and actuators For each controller, specify a 2x1 cell which contains the indices of commanded actuators respectively inputs (u) received by that controller. For example:

Open-loop controller: command actuators [1,3] (Ip and ECCD power) in open-loop from inputs [1,2]

    map_openloop = {[1 3],[1 2]}; %

Beta feedback controller command actuator 2 (ECH power) by input 3 (the reference beta)

    map_betafb = {2,3};

Combine them both into the `controllers.mapping` parameter

    config.controllers.mapping = {map_openloop,map_betafb};

Now inspect the `controllers` structure

    disp(config.controllers)

        openloop_controller: [1x1 struct]
            beta_controller: [1x1 struct]
                    mapping: {{1x2 cell}  {1x2 cell}}

Prepare the model

    [model,params,init,g,v,U] = build_RAPTOR_model(config);

    % Define inputs
    % The inputs are now defined, following the above controller mapping, as
    % |u=[Ip,P_{ECCD},|$\beta_{ref}$|]|

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

    % Off-axis power (feedforward)
    U(2,:) = 1e6*ones(size(params.tgrid)); % off-axis power
    U(2,params.tgrid<0.05) = linspace(0,1e6,sum(params.tgrid<0.05)); % ramp

    % Beta reference: 1% (feedback)
    U(3,:) = 0.01*ones(size(params.tgrid)); % beta reference (%)

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

    % Run simulation
    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.17    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     1.3    3 8.3e-15      10       1     169    1.48    7.65     160    7.27    3.59     3.7 5.3e+00    0.65    0.65    1.00 2.6e+01
         21     1.8    3 5.2e-12      20       1     250    7.96    18.1     224    4.93    2.39     3.3 3.4e+00    1.71    1.71    1.00 2.0e+01
         31     2.3    2 2.8e-10      30       1     250    20.5    31.5     198    4.93    1.86    3.39 1.6e+00    3.37    3.37    1.00 1.0e+01
         41     2.8    2 8.0e-12      40       1     250    36.6    43.1     170    4.93    1.62    3.68 9.5e-01    4.95    4.95    1.00 6.3e+00
         51     3.4    4 3.4e-14      50       1     250    56.2    52.5     141    4.93     1.5     4.3 6.0e-01   13.09   13.09    1.00 4.5e+00
         61     3.9    2 5.6e-09      60       1     250    96.6     115    38.4    4.93    1.42    5.98 3.3e-01   12.84   12.84    1.00 2.7e+00
         71     4.5    2 2.2e-10      70       1     250    98.7     131    19.8    4.93    1.36    7.67 2.9e-01   15.21   15.21    1.00 2.5e+00
         81       5    2 2.7e-12      80       1     250    96.1     130      24    4.93    1.31    8.48 2.7e-01   15.38   15.38    1.00 2.3e+00
         91     5.5    2 6.5e-12      90       1     250    95.4     129    25.1    4.93    1.27     8.8 2.5e-01   15.36   15.36    1.00 2.1e+00
    saved controller data to /afs/ipp-garching.mpg.de/home/f/ffelici/matlab/RAPTOR/personal/beta_controller_data
        101       6    2 7.0e-14     100       1     250    95.2     130    24.8    4.93    1.23    9.78 2.2e-01   15.40   15.40    1.00 2.0e+00

## Plot results

    out = RAPTOR_out(simres,model,params);
    clf;
    subplot(321); plot(out.time,out.U(1,:)); title('Ip');
    subplot(322); plot(out.time,out.U([2:3],:)); title('P_{aux}')
    subplot(323); plot(out.time,[out.beta;U(3,:)]); title('\beta and \beta_{ref}')
    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');

end of tutorial