Using the detachment controller

N.b. currently this code has been developed for and tested only in 1D and with a single processor.

Overview

The DetachmentController is designed to adjust plasma parameters dynamically to maintain the detachment front at a desired position. It uses a proportional-integral-derivative (PID) controller approach, customized with specific logic for detachment control in plasma simulations. The controller modifies the source terms for energy or particle density based on the difference between the current and desired detachment front locations.

Key Components

• Constructor: Initializes the controller with simulation options, including setpoints, control parameters (gain, integral time, derivative time), actuator type (power or particles), and other relevant settings.

• transform(Options& state): Main method called at each simulation timestep to calculate the control action. It performs several steps:

  1. Computes the current location of the detachment front based on plasma densities.

  2. Calculates the error between the desired and current detachment front locations.

  3. Updates the control signal based on the PID algorithm.

  4. Applies the control signal to adjust the source terms accordingly.

• outputVars(Options& state): Outputs diagnostic variables related to the control process, such as the detachment front location, control signal, and PID terms.

• restartVars(Options& state): Manages state variables for restarting the simulation from saved states, ensuring the controller can smoothly continue from a previous point, without perturbations due to reinitialising internal variables.

Usage

To use the DetachmentController in a BOUT++ simulation, you must provide a setpoint for the detachment front location, choose the actuator type, and configure the PID controller parameters (gain, integral time, derivative time). The input for a minimum working example would be something like

[hermes]
components = (d+, d, e, detachment_controller, ...)

[detachment_controller]
; Desired location of the detachment front, in metres from the target
detachment_front_setpoint = 1.0
; What is the main neutral species?
neutral_species = "d+"
; What should we adjust to control the detachment front?
actuator = "power"
; What species should we look at for `source_shape`
species_for_source_shape = e
; What species should we control, and how much control should apply
; to each species?
species_list = e, d+
scaling_factors_list = 1.0, 1.0

; What value to use for Kp? See the rest of the documentation
controller_gain = 1.0

[Pe]
function = 1.0
source_shape = (power_flux*2/3 / (mesh:length_xpt))*H(mesh:y_xpt - y)

[Pd+]
function = Pe:function
; No Pd+ source or source_shape

Top-level overview of the algorithm

• Detachment Front Location Calculation: Determines the point where the neutral density exceeds the electron density, which is taken as the position of the detachment front.

• PID Controller Implementation: Adjusts the control signal based on the error between the current and target detachment front locations. Supports both velocity and position form PID control.

• Actuator Type: Defines whether the control action adjusts power or particle sources in the simulation.

• Control Signal Application: Modifies the source terms for specified plasma species based on the control signal, influencing the plasma dynamics to achieve the desired detachment front location.

Fine details

Finding the detachment front

We define the detachment front as the point where the neutral density (N_neutral_species) becomes larger than the electron density (Ne). There’s lots of other (better) definitions, but this one is the easiest to calculate. For each timestep, the controller starts at the divertor target and checks whether a cell has \(n_N > n_e\). The first time it finds a cell fulfilling that condition, it calculates to first order the position where the cross-over would occur, by comparing the last point where \(n_N > n_e\) (point 1) with the first point where \(n_e > n_N\) (point 2);

\[x_D = \frac{(x_1 \times n_{e,2} - n_{e,1} \times x_2) - (x_1 \times n_{n,2} - n_{n,1}*x_2)}{(n_{n,1} - n_{n,2}) - (n_{e,1} - n_{e,2})}\]

Note that, since the \(y\) values aren’t stored directly, we compute it from \(dy\) while iterating along cells.

Calculating the feedback response via PID control

Once we have the detachment position, we need to control its position. Firstly, we define the error \(e = x_S - x_D\) where \(x_S\) is our setpoint and \(x_D\) is the actual position of the detachment front. We also fetch the current time \(t\).

If starting a new simulation, we don’t actually want to start controlling the detachment front straight away, but instead wait for settling_time (in seconds). Until then, we just set control = initial_control.

Once we’re through the settling_time, we only apply control if the time has changed by some min_time_for_change and if the absolute change of the error is at least min_error_for_change. Usually we set min_time_for_change to be a small, but non-zero number, which prevents us from dividing by zero when evaluating the derivative of the error. Often the min_error_for_change is zero, although you can experiment with this and see if it gives you improved performance of the controller. Once these conditions are met, we update the control value.

There are two options for the control algorithm to use — the ‘position’ and ‘velocity’ forms of the PID. These are, respectively,

\[C = C_0 + s\times K_c\times \left[e + \frac{1}{\tau_I} \int e\cdot dt +\tau_D \frac{d e}{d t} \right]\]

and

\[C = C_{prev} + s\times K_c\times \left[\Delta e + \frac{\Delta t}{\tau_I} e +\tau_D \Delta \left(\frac{d e}{d t}\right) \right]\]

where * \(C\) is the control value applied to the source * \(s\) is the response sign (\(+1\) for particles or \(-1\) for power) * \(K_c\) is the controller gain * \(\tau_I\) is the integral time * \(\tau_D\) is the derivative time * \(\int e\cdot dt\) is the error integral * \(\frac{d e}{d t}\) is the error derivative * \(C_0\) is the control offset (usually zero, unless you’re running in proportional-only mode and need a finite control when the error is zero) * \(C_{prev}\) is the previous value of \(C\)

In the position form, we directly update the control value, while in the velocity form we calculate the change of the control value and add it to the previous control value. These two approaches should be largely equivalent, although the velocity form has the advantage of avoiding integral windup when the control is at either its min or max value. Conversely, with the position form we can apply other anti-windup schemes which explicitly change the error integral. A particularly useful anti-windup method (only available in position form) is reset_integral_on_first_crossing, where the error integral is reset when the detachment front first reaches the setpoint.

In either scheme, the control value is bounded by minval_for_source_multiplier and maxval_for_source_multiplier, which can be used either to ensure that the system stays within physical, engineering or numerical-stability constraints.

Understanding PID control and tuning the coefficients

See en.wikipedia.org/wiki/Proportional-integral-derivative_controller. Most relevant sections are

• Controller theory

  • Proportional term

  • Integral term

  • Derivative term

• Loop tuning

• Ziegler-Nichols method

Input Parameters

• ``detachment_front_setpoint``: The desired position of the detachment front from the divertor target, measured in meters from the divertor target. It represents the target location where the control system aims to maintain the front.

• ``velocity_form``: A boolean flag indicating whether to use the velocity form (if true) or the position form (if false) of the PID controller.

• ``min_time_for_change``: The minimum time interval, in seconds, before the control signal can be updated. This parameter prevents divide-by-zero errors when evaluating the error derivative.

• ``min_error_for_change``: The minimum change in error required before the control signal is updated. This term is mostly left for experimenting.

• ``minval_for_source_multiplier``: The minimum value that the control signal (source multiplier) can take.

• ``maxval_for_source_multiplier``: The maximum value for the control signal (source multiplier).

• ``species_for_source_shape``: Specifies the plasma species from which to select the source shape. The feedback source is the product of the control signal (source multiplier) and the source shape. If actuator='power' this will set source_shape=Ps::source_shape where s=species_for_source_shape, and if actuator='particles' this will set source_shape=Ns::source_shape.

• ``neutral_species``: Indicates the main neutral species in the plasma. This species is used to determine the location of the detachment front.

• ``actuator``: Defines the actuator to be adjusted to control the detachment front position. Options include ‘power’ for energy sources or ‘particles’ for particle density sources.

• ``initial_control``: The initial value for the source multiplier (control signal) at the start of the simulation.

• ``control_offset``: The expected control value when the error equals zero. It serves as a baseline around which the control signal is modulated. Only used in position form.

• ``settling_time``: The time allowed for the system to settle before activating certain control terms, measured in seconds. It delays the start of control actions to ensure initial transients have subsided.

• ``ignore_restart``: A flag to ignore the restart file, mainly useful for development purposes. It forces the controller to initialize from the provided settings rather than loading a previous state.

• ``reset_integral_on_first_crossing``: Resets the error integral to zero when the detachment front first reaches the desired position. This feature can help prevent integral wind-up and improve control stability.

• ``controller_gain`` (Kc): The proportional gain of the PID controller. It determines the strength of the response to the error.

• ``integral_time``: The integral time of the PID controller, which influences the rate at which the integral term accumulates error over time.

• ``derivative_time``: The derivative time of the PID controller, affecting how strongly the controller reacts to the rate of change of the error.

• ``buffer_size``: The number of points to store for calculating derivatives. It determines the size of the window over which the derivative of the error is computed. Usually a value of around 3 to 5 seems to filter out the worst of the noise due to small time-steps. Very large values will potentially introduce a destabilising lag in the derivative response (and also will eventually start to impact computational performance).

• ``species_list``: A comma-separated list of species to which the PI-controlled source will be applied.

• ``scaling_factors_list``: A comma-separated list of scaling factors corresponding to each species listed in species_list. These factors adjust the magnitude of the control action applied to each species.

• ``diagnose``: Enables the output of additional diagnostic variables related to the control process if set to true.

• ``debug``: Controls the level of debugging information printed to the screen. A value of 0 disables debugging output, 1 enables basic output, and 2 provides extensive debugging information.