Time-dependent SOLPS-ITER simulations of the tokamak plasma boundary for model predictive control using SINDy

Time-dependent SOLPS-ITER simulations have been used to identify reduced models with the sparse identification of nonlinear dynamics (SINDy) method and develop model-predictive control of the boundary plasma state using main ion gas puff actuation. A series of gas actuation sequences are input into SOLPS-ITER to produce a dynamic response in upstream and divertor plasma quantities. The SINDy method is applied to identify reduced linear and nonlinear models for the electron density at the outboard midplane ne,sepOMP and the electron temperature at the outer divertor Te,sepdiv . Note that Te,sepdiv is not necessarily the peak value of T e along the divertor. The identified reduced models are interpretable by construction (i.e. not black box), and have the form of coupled ordinary differential equations. Despite significant noise in Te,sepdiv , the reduced models can be used to predict the response over a range of actuation levels to a maximum deviation of 0.5% in ne,sepOMP and 5%–10% in Te,sepdiv for the cases considered. Model retraining using time history data triggered by a preset error threshold is also demonstrated. A model predictive control strategy for nonlinear models is developed and used to perform feedback control of a SOLPS-ITER simulation to produce a setpoint trajectory in ne,sepOMP using the integrated plasma simulator framework. The developed techniques are general and can be applied to time-dependent data from other boundary simulations or experimental data. Ongoing work is extending the approach to model identification and control for divertor detachment, which will present transient nonlinear behavior from impurity seeding, including realistic latency and synthetic diagnostic signals derived from the full SOLPS-ITER output.


Introduction
Tokamak fusion reactors require active control to ensure safe and reliable operation. Robust controllers are already part of standard tokamak operation, e.g. feedback control using a set of external magnetic coils to stabilize the vertical position. In the area of shape, current, and position control, a range of control strategies have been applied, including classic proportional integral derivative (PID), linear-quadratic regulator (LQR) [1,2], H-infinity [3], and model predictive control (MPC) techniques [4]. Feed-forward or PID methods are regularly used to control core parameters such as the line-averaged density and plasma pressure using heating systems and gas puff actuators [5]. Active control of the divertor particle and energy fluxes will be required in power-producing tokamak reactors, as the plasma-facing components cannot withstand the unmitigated loads. Recently, significant advances have been made in developing active control of the boundary plasma and conditions at the divertor using fast actuation of extrinsic impurity gas valves [6][7][8].
These efforts have demonstrated successful control of the divertor fluxes, in some cases while maintaining acceptable core confinement. Core-edge integration must be considered as reduction of the divertor power loads is often associated with a deleterious effect on core confinement due to radiation in the confined plasma region [9][10][11][12][13]. To date, divertor heat and particle flux control generally uses PID methods with empirically tuned gains. This tuning requires many pulses under relatively forgiving conditions, which may not be possible in future burning plasma reactors, although it may be achievable to estimate the gains and transfer them to another experiment [14].
MPC offers potential advantages over PID control in that nonlinear models can be applied, and the model can be retrained during operation. The response model used in the controller can be developed using data from both simulations and experiments. Developing a reliable reduced model of a high-dimensional nonlinear phenomena, such as tokamak boundary detachment, is a significant challenge. Evaluation of the nonlinear model and minimization of a cost function for optimal actuation also make MPC more computationally expensive than PID approaches. Here we present methods for developing physics-informed reduced models of the boundary plasma response to gas puff actuation, which can be rapidly evaluated in a model-based control loop.
Presently, boundary plasma transport models span a wide range of complexity and computational cost. Kinetic transport simulations can capture plasma turbulence, but often do not contain the detailed plasma-neutral and plasma-surface interactions that are important in the open field line scrape-off-layer (SOL) and divertor regions. In addition, they are typically too computationally expensive to evolve the kinetic profiles over timescales relevant to cross-field transport. Analytic descriptions such as the two-point-model [15] aid in the understanding of complex SOL physics, but require parameters that describe the volumetric momentum and pressure loss along a field line that must come from experimental or 2D simulation data and do not generally contain actuator variables such as the gas puff or account for time-dependence [16]. Fluid plasma boundary transport simulations (e.g. SOLPS [17], UEDGE [18], etc), on the other hand, can contain detailed plasmaneutral and plasma-surface models while using ad-hoc models to represent cross-field turbulent transport. In this work the coupled 2D fluid plasma and kinetic neutral transport code SOLPS-ITER is applied [19].
While SOLPS-ITER can simulate the boundary plasma response to gas puff actuation, the simulations are still too computationally expensive to include in tokamak plasma control loop. To overcome this, a data-driven methodology is used that yields interpretable (not black box) models given a dynamic response to actuator inputs from SOLPS-ITER. The methods described here are not specific to a given boundary transport code, and would also apply to experimental datasets given sufficient temporal diagnostic resolution. To develop the general approach, the complex and computationally intensive problem of divertor detachment control using extrinsic impurity seeding is reduced to the simpler and less nonlinear problem of control of upstream (midplane at the core-SOL interface) and downstream (divertor) plasma quantities using main ion gas puff actuation.
System identification is performed using the method of Sparse Identification of Nonlinear Dynamics (SINDy) [20]. SINDy has been applied to complex systems such as plasma convection [21], magnetohydrodynamics [22], fluid dynamics [23][24][25][26][27][28][29][30][31][32][33][34][35], and MPC [36]. The method yields a system of differential equations that describes the dynamical evolution of a set of simulation or experimental data. While SINDy has often been applied to identify the 'true' dynamical equations from data, here we are interested in finding reduced models that describe the evolution and actuation response to a small number of variables selected from the full SOLPS output, similar to the approach taken in [37]. The success of parameterizing SOLPS [38,39], in terms of core conditions and external actuators, suggests an opportunity for reducing the problem from 2D partial differential equations to a small number of ODEs that can be rapidly evaluated in a controller.

SOLPS-ITER simulation setup
The simulations of the plasma boundary are performed using version 3.0.8 of the SOLPS-ITER code. SOLPS-ITER is a coupling of the B2.5 code, which solves a set of 2D fluid plasma conservation equations, and the EIRENE kinetic neutral transport code which describes plasma-neutral and plasmamaterial interactions [19,40]. The numerical approach for the plasma solution is a finite volume method on a field aligned structured grid, while the EIRENE code uses a Monte-Carlo approach on an unstructured mesh, as illustrated in figure 1. The plasma fluid equations have the form of convectiondiffusion equations, which are each solved implicitly, with Picard iteration to handle the nonlinear coupling. A detailed description of the numerical methods is presented elsewhere, e.g. [41]. In a simplified description, the evolution equations for the 2D axisymmetric distribution of the state variables of plasma density, parallel momentum, and internal energy x = [n, mnV ∥ , 3 2 nT] are given as ∂n ∂t Here n, T, and V are the plasma density, temperature, and velocity, respectively, with subscripts e and i used to denote electrons and ions. The particle and parallel momentum flux densities are Γ n and Γ m , and q + 3 2 pV is the internal energy flux density. The static pressure is p = nT, and Q represents collisional energy exchange between the plasma species. The source terms S include plasma-neutral interactions calculated by EIRENE as well as friction and viscous terms. Additional auxiliary equations are used to describe charge continuity and the cross-field components of the velocity. In general the source terms as well as the diffusivities are nonlinear functions of the plasma parameters.
Although SOLPS-ITER is widely used for steady-state analysis, time-dependent simulations are relatively rare. Previous work on transient response includes simulation of edge localized modes [42,43], transition to detachment [44,45], and gas injection [46,47]. In many cases only the steady-state solution is of interest, and significant development effort has gone into accelerating the time advance to a steady-state solution, with non-physical dynamics in the transient phase [48,49].
To produce meaningful dynamics, the acceleration schemes described in [48] are deactivated, and the time-step is reduced and number of nonlinear relaxation iterations between the equations increased until the transient results converge. For the cases shown here, dt = 1 × 10 −5 s is used, with 16 internal iterations (the 'nstg2' input variable in the SOLPS code). Presently the timescale of the neutral particles leaving the gas puff (and recycling surfaces) to be ionized in the plasma is assumed to be much shorter than the timescales of the plasma evolution. This means that EIRENE can be run in a quasi steady-state manner, with the source terms in the plasma equations updated every B2.5 timestep. Studies with timedependent neutral evolution will be presented in a future publication. Recent advances have investigated alternate timeadvance strategies for SOLPS [50]. Here we focus on general methods for developing reduced models and control strategies given a time-series dataset, which can be applied to other numerical implementations or experimental data. The magnetic topology is taken from an EFIT [51] reconstruction of shot 174 310 of the DIII-D tokamak at 3500 ms. The plasma-facing component geometry, magnetic topology, and puff and pump locations are shown in figure 1. The simulation setup was not chosen to match specific experimental conditions, but instead to use values which yield approximate typical conditions relevant to fusion plasmas while simplifying the numerical model. The input power crossing the core boundary is 6 MW, equally split between ions and electrons. Only deuterium ions are considered, with a core boundary condition of 7.5 × 10 20 ions s −1 (corresponding to neutral beams fuelling), and 3 cm decay length boundary conditions are used at the outer radial surfaces. A gas puff of molecular deuterium is introduced, Γ D2 , with an initial magnitude of 2.5 × 10 21 atoms s −1 . Neutral particles are removed from the pumping surface indicated in figure 1, which has a recycling coefficient of 0.5, all other surfaces are fully recycling. Crossfield diffusivities of D ⊥ = 0.3 m 2 s −1 and χ ⊥ = 1.0 m 2 s −1 are used, corresponding to the typical values used for ITER simulations using SOLPS-ITER [52]. The resulting power decay length is ∼4 mm. Cross-field plasma drifts are not activated.
The simulations take approximately 5-20 s per timestep, using 16 core MPI parallelization, with the lower range corresponding to low density high temperature plasmas where the neutral mean free path is reduced. The number of Monte Carlo particles in EIRENE is adjusted such that approximately 50 percent of the time is spent in EIRENE (30 000 particle for divertor recycling and gas puff sources and 3000 for outer boundary sources). The resulting steady-state profiles of electron density and temperature at the outboard midplane, and the outer divertor are shown in figure 2.

Generation of system dynamics
To develop reduced models, a subset of state variables are selected from the full 2D time series data of the SOLPS state. To illustrate a relevant set of variables for a reduced model capturing dynamic behavior of upstream and downstream conditions, we select the separatrix electron density at the outboard midplane (n OMP e,sep ) and the separatrix electron temperature at the outer divertor (T div e,sep ). The former describes an important quantity for tokamak operation that is actuated by a main ion gas puff, and is directly related to the divertor conditions [15]. The latter is of interest for survival of the divertor, and has been shown to be proportional to a large number of quantities of interest to describe the divertor fluxes and other quantities from experimental and simulation data [38,39]. The approach is then to use SOLPS-ITER to simulate the response of the two reduced model state variables caused by actuation of the main ion gas puff, then use these time series data to fit a reduced dynamic model.
From the initial conditions shown in figure 2, the gas puff actuator is varied to produce a plasma response. Figure 3 shows results from a set of simulations where the gas puff magnitude was slowly varied between zero and 1 × 10 23 atoms s −1 . The blue lines indicate the system-state trajectories, while the black crosses are results from steady-state periods. The state variables have a nonlinear relationship to the input variable, particularly at low values of the gas puff. A baseline actuation level of 2.5 × 10 21 atoms s −1 , indicated by the vertical dashed line, is selected as an initial state from where the solution is converged to steady-state. Dynamic actuation near this value is useful for identifying approximate linear and nonlinear models. The actuation range is limited to 0.05 × 10 22 atoms s −1 of the steady-state level. The n OMP e,sep series has a noise level of approximately 0.5%, while the T div e,sep series has a larger noise of approximately 10%. This noise is largely due to Monte-Carlo noise from the EIRENE simulation, which can be reduced at the cost of computational time, but is useful to assess the robustness of the system identification methods.
The phase space diagram for this sequence is shown in figure 5 as the blue curve. It can be seen that the variables are tightly coupled, with a relationship n OMP e,sep = f(T div e,sep ) ∝ (T div e,sep ) −1/2 , as indicated by the dashed black line. This scaling is expected from the 2PM relationship between these variables, which can be written as The upstream parallel heat flux density q || ∝ P SOL /χ 1/2 ⊥ is largely determined by the power crossing the separatrix P SOL and the choice of cross-field diffusivity χ ⊥ , which are held constant across each simulation. The terms (1 − f cool ) and (1 − f mom ) are ad-hoc variables describing the loss of energy and momentum in the flux tube, which decrease from unity when there is significant radiation, and as the plasma temperature drops and plasma-neutral interactions become more significant [38]. With main ion puffing, it is expected that the loss terms do not have a large impact until T div e,sep becomes very small [53], although deviations from the exact 2PM model scaling are observed in the dataset. Generally the loss terms are found to be strongly correlated with T div e,sep , from both experimental and simulation data [38,39,52]. A dynamic form of the two point model, using 0D particle inventories to add gas puff actuator variables, has been derived in [16]. The relatively simple algebraic relationships between the selected state variables suggest that identifying ODE systems describing the dynamics might be possible, given proper selection of the form of the right-hand side of the ODE, as discussed in section 3.2. Note that the tight functional relationship between the two variables means that a functional evaluation could be used in the simplest case, fitting only the dynamics of a single variable. However, this approach will not work for more complex scenarios.
In the present work, we focus on developing and demonstrating methods for training linear and nonlinear reduced models using data from main ion gas actuation while holding other simulation parameters constant. It is worthwhile, however, to consider future applications of the methods to more complex systems. The addition of nitrogen impurity seeding in the simulation results in significant deviation from the (T div e,sep ) −1/2 scaling. This is shown in figure 5 as green and red lines for increasing and decreasing nitrogen puff, respectively. The impurity seeded trajectories in phase space are highly nonlinear and exhibit hysteresis. Identifying a model (or multiple Figure 5. Phase space diagram for several sets of SOLPS-ITER simulations using dynamic gas puff actuation sequences. Blue, orange, and yellow curves use main ion actuation and correspond to different cross-field diffusivities in the simulation. Fits to n OMP e,sep ∝ (T div e,sep ) −1/2 are shown as dashed black lines. Green and red curves are for simulations using nitrogen impurity seeding, for increasing and decreasing levels, respectively. models) capturing the impurity response for a reactor detachment controller will be challenging. Phase space curves for additional simulations with χ ⊥ = 0.5 m 2 s −1 and 2.0 m 2 s −1 are shown as orange and yellow traces, respectively. The distinct curves, still approximately maintaining the 2PM scaling, indicate that parameterization of the reduced model coefficients by transport quantities might be possible. Identification of impurity seeding actuation models and parameterized models is beyond the scope of this paper.

System identification using SINDy
System identification is performed using the method of Sparse Identification of Nonlinear Dynamics (SINDy) [20]. SINDy is used to identify a system of sparse ODEs that describe the dynamical evolution of a set of simulation or experimental data. The original SINDy approach has been extended to include actuator response and model predictive control application [36,54], and enforcing known analytic constraints such as energy conservation [22,23], as well as other advancements as described in [55]. In the SINDy approach, the state variable dynamics with actuation are approximated aṡ where x and u denote respectively the state and actuator variables and Θ(x, u) represents a library of candidate functions over the output state x and input actuator u for the right-hand side of the ODE, with (potentially) sparse coefficients Ξ. The choice of these candidate functions can be informed by the form of the underlying equations, e.g. in [20] trigonometric functions are shown to result in a sparse set of equations when the underlying ODE has the form of sin(x). Here we select a polynomial basis, motivated by the form of the 2D and 1D fluid plasma equations. Here p n (x, u) denotes the collection of degree n multivariate monomials in the elements in x and u, e.g. p 0 (x, u) = 1, , etc, for a system with two state variables and a single actuator.
The derivativeẋ can be approximated numerically, in this case by a first order finite difference, from the time series data x. The coefficients Ξ are then selected to accurately describeẋ by solving a sparsity-regularized least-squares problem, i.e.
where a regularization term R(·) is included to promote sparsity in Ξ so that spurious high-order terms that marginally improve the fit are discarded. This optimization problem can be solved using methods such as sequentially thresholded least squares (STLS) [56], sparse relaxed regularized regression (SR3) [57], or LASSO [58], and the exact form of R(·) depends on the choice of the optimization solver. Here we apply STLS method, which promotes the sparsity of Ξ by solving the least-squares problem without regularization, removing the candidate functions associated to coefficients of magnitude smaller than a prescribed threshold λ > 0, and repeating this process until convergence. The choice of λ affects the sparsity and accuracy of the fit, with larger values leading to a sparser Ξ with a potentially higher residual. The 'best' choice of λ depends on the application, and can be evaluated by calculating the Pearson correlation coefficient for a range of λ, as discussed below.
By promoting sparsity, the SINDy method generally results in interpretable equations which can be easily integrated and evaluated for stability and convergence, as opposed to blackbox methods for generating reduced models such as neural networks. In practice, numerical evaluations ofẋ from noisy data x can be challenging, and methods such as polynomial filtering derivatives (e.g, Savitzky-Golay) or total variational derivatives [59] are used. Recent developments in kernelbased approaches [60,61], and weak form ODEs [62] are also being investigated. In this work we simply apply a windowed smoothing to the time-series data.
An example of the procedure applied to the input dataset described in section 3.1 is shown in figure 4. A single model is identified using the four SOLPS-ITER datasets and integrated from the initial condition of each series using the actuation shown in figure 4(c). The result is shown in figures 4(a) and (b) as black dashed lines. The resulting L 2 error in the model is shown in figure 4(d).
To produce this model, the PySINDy package, version 1.7.1, was used [55,63]. The input actuator u = Γ D2 and state variables x 1 = n OMP e,sep , x 2 = 1/T div e,sep were normalized by 1 × 10 18 , 1 × 10 16 , and 1 × 10 −4 , respectively. As discussed in [20], the fit can be improved by selecting a candidate function library informed by the underlying form of the dynamical equations. The choice of x 2 = 1/T div e,sep is motivated by the inverse relationship between the state variables, and was found to improve the accuracy of the model results with the data as compared to x 2 = T div e,sep . The numerical derivative was calculated using a Savitzky-Golay smoothed finite difference with 51 time points on a uniform grid with interval t = 10 −4 sec, with a threshold parameter λ = 0.1. The window was selected to reduce the noise in the evaluation of the derivative of the temperature data, with the smoothed dataset shown as a dark orange curve in figure 4(b) (time averaged using the same window as the finite difference). The resulting model has the forṁ where the subscript of ξ indicates the power of the the state variable and actuator terms [x 1 , x 2 , u]. The threshold λ = 0.1 was selected to yield a linear system. A nonlinear model would have higher order terms, e.g. ξ 120 x 1 x 2 2 . The impact of the sparsity parameter λ can be assessed by calculating the Pearson correlation coefficient ρ between the numerical evaluation ofẋ and the identified model Θ(x, u)Ξ for a range of λ. The correlation and the number of non-zero terms in the model are shown in figure 6(a). For this small actuation range, the model can be made sparse without significant loss of accuracy. Using λ = 0.1, the resulting linear model captures the features of the simulation dynamics with L 2 error in n OMP e,sep of < 0.5 %, and < 10 % and 5 % in the raw and smoothed T div e,sep series, respectively, as shown in figure 4(d). The linear model works well to describe actuation magnitudes within a small range ∆u of the initial value. For actuation sequences with larger actuation range, sampling a larger portion of the nonlinearity shown in figure 3, a nonlinear model is required. This is illustrated in figure 7. A set of three SOLPS-ITER simulation datasets are concatenated, corresponding to actuations of a multi-frequency sinusoid, a triangle wave, and a series with u ∝ sin(t) + At, as shown in figure 7. The ∆u of each series is 2 × 10 21 , 5 × 10 21 , and 5 × 10 22 atoms s −1 , respectively. Applying the linear model above, results in the dashed black lines for n OMP e,sep in figures 7(a) and (d), and the green lines for T div e,sep in figures 7(b) and (d). The model error increases with ∆u, with poor performance over the final series. A nonlinear model can be trained over the same sequences, using λ = 10 −3 . The results are shown as dashed red lines for n OMP e,sep and purple lines for T div e,sep . The nonlinear model better describes the dynamics over the high ∆u sequence with raw error less than 10%. Figure 8 shows the identification of a nonlinear model using λ = 10 −4 for a larger dataset with a large actuation range. The Pearson correlation coefficient and number of nonzero model terms for this series as functions of λ are shown in figure 6(b). In this case, the model correlation (accuracy) clearly decreases with the number of terms. It is found that to integrate the ODE system over all time series with acceptable L 2 error ( figure 8(d)), the number of terms can be reduced by only ≈ 10 from the full, non-sparse system. This indicates a possible limitation of the method for this simulation data, in that it may be difficult to identify a sparse system of equations that has high accuracy when integrated over a long time range. The lack of sparsity alone is not a problem, but higher order terms can make the system more unstable. However, as long as the prediction window with sufficient accuracy exceeds the control horizon with feedback then the model may be acceptable. A different selection of candidate functions for the ODE right-hand side (e.g. fractional powers or using x 2 = (T div e,sep ) −1/2 as motivated by the form of the 2D and 2PM equations,) other sparse regularization solvers (e.g. LASSO, SR3), or more noise robust approaches may also yield sparser and lower order models. Adaptive retraining of the model may also be an effective approach.

Model retraining
The model training using SINDy can be applied based on an error threshold to trigger adaptive retraining over an increasing or windowed dataset. The SINDy procedure is not computationally expensive for the dimension of the datasets shown here, taking ≈5 ms to identify a model over 1 second of simulation data in a coupled MATLAB and pySINDY implementation that has not yet been optimized for performance. Figure 9 shows an example of the retraining procedure. A model is identified on the first 100 ms of data using SINDy, after which the model is re-initialized every 150 ms (corresponding to feedback) and simulated forward in time. The simulated error signals are indicated by the colored traces in figures 9(d) and (e). If the error in the simulated signal exceeds a preset threshold, 5 % in n OMP e,sep and 10 % in T div e,sep , a retraining is triggered. Each retraining is indicated by a vertical red line in figure 9. The simulated behavior past the error threshold is drawn as grey lines to show the model behavior. In this example, the retraining uses the dataset from t = 0 sec, but a shorter window could also be used. Initially a linear model is identified, which exceeds the error threshold when ∆u grows too large (t ≈ 0.52 s). A series of retrainings are triggered as the actuation range increases until the final model (identified at t ≈ 1.1 s) does not exceed the error thresholds for the remainder of the sequence.

Model predictive control
The reduced models can be used to design actuation sequences to control the state variables to match a given setpoint trajectory. If the reduced model is linear, traditional techniques such as LQR and PID are suitable. More generally, an MPC strategy, which is applicable for nonlinear models, can recover the optimal linear actuation, albeit with greater computational cost. Here we apply an MPC scheme using a rolling time horizon to yield a series of piece-wise constant actuation levels to best match a setpoint trajectory, with or without feedback. This approach also allows for straightforward inclusion of constraints in the optimization process, for example, minimum and maximum actuation levels and rates of change. For simplicity, we assume a fully observable system without latency in the actuator or measurement. Figure 10 schematically illustrates the MPC strategy. Given a setpoint trajectory x * (t), a series of piecewise constant actuation levels u j of duration ∆t act are used to drive the model over ∆t predict . The optimal actuation levels are those that minimize a cost function J(x * , u j ). The optimal actuation is then applied across ∆t control , and the process is repeated. The general form of the cost function used is given as The first term is the L 2 norm of the difference in the predicted x and setpoint trajectories The other terms represent penalties on the actuation level J u = ∑ Q u j |u j | 2 and J ∆u = ∑ Q ∆u j |∆u j | 2 . Constraints, e.g. minimum and maximum actuation levels, can be implemented via the last term J con. . Consideration of core-edge integration can be introduced using the cost function, by including an appropriate state variable (e.g. core radiation or position of peak radiation along the separatrix leg) which has constraint terms that inform the optimal actuation while the setpoint trajectory is defined in the target variable (e.g. divertor heat flux). This will be explored in the future. Given the tight correlation between the state variables n OMP e,sep and T div e,sep when using main ion actuation, it may not be possible to reach an arbitrary state such as low T div e,sep and high n OMP e,sep in a reasonable time and actuation level. An example application of the MPC approach is shown in figure 11. A 0.4 sec series of SOLPS-ITER data is used to train a linear model. The MPC strategy is used to determine an actuator sequence of piecewise constant segments of ∆t act = 25 ms to produce a density setpoint trajectory n OMP e,sep ∝ sin(t), shown as the thick black line in figure 11(a). With a linear model, the MPC actuation ( figure 11(c), blue line) considering only the first term in the cost function (equation (7)) yields a piecewise approximation of the LQR actuation (red line). In both cases the n OMP e,sep prediction well matches the setpoint trajectory. The setpoint trajectory was chosen to result in negative actuation when evaluated using a linear model. The green curves show the MPC strategy applied with an additional constraint, u > 0, avoids this nonphysical actuation while resulting in increased error during this time.
MPC has been used to perform feedback control of a SOLPS-ITER simulation, as shown in figure 12. The MPC strategy is implemented as a module of the Integrated Plasma Simulator (IPS) [64]. The IPS manages the overall time advance and program execution, with the main loop launching a SOLPS-ITER simulation, evaluating the cost function and optimal actuation, and updating the input files appropriately. The setpoint trajectory (black), predicted optimal actuation (red), and resulting SOLPS-ITER n OMP e,sep (blue) and T div e,sep (orange), are shown in figures 12(a) and (b). The actuation input for each interval is shown in figure 12(c). The setpoint trajectory is well reproduced by the simulation, with reasonable agreement between the prediction from a nonlinear model and the actual state evolution.

Summary and conclusions
The SOLPS-ITER transport code has been used to simulate the dynamic response of the plasma boundary to gas puff actuation. The high dimensional dataset is used to identify reduced models for a limited set of state variables of interest, here n OMP e,sep and T div e,sep . The SINDy method is applied to identify a set of linear and nonlinear ODEs that best match the SOLPS-ITER dynamics, resulting in interpretable models that can be evaluated using standard analysis and integration methods.
The choice of state variables are found to have a tight functional relationship over the range simulated, which is expected from the simplified system (i.e. main ion fuelling vs impurity seeding) and the variables connected along a flux tube. The nonlinear system can be well matched over a small actuation range using a linear model, or a nonlinear model can be fit over a wider actuation range. Given the low computational cost of model identification, it is expected that retraining can be performed periodically in real time with the trained model is no longer accurate. The linear and nonlinear models are used in an MPC strategy to match given setpoint trajectories, with or without constraints. The controller has been implemented as a module in the IPS, and used with feedback on a SOLPS-ITER simulation to achieve dynamic density control.
As a next step, the time-dependent SOLPS-ITER dynamics should be validated against experimental data. Given the challenges of validating even steady-state SOLPS-ITER simulations against experimental data, it is not expected that close quantitative agreement will be obtained. Capturing additional physics such as time-dependent neutral dynamics or crossfield drifts may make accumulation of sufficient simulation data intractable. However, if the simulations capture the correct trend in the plasma response to actuation, then control using such model with feedback may be achievable. Identifying models from experimental data, or supplementing the simulations with experimental data may also result in improved models.
Further necessary development includes implementation of synthetic diagnostics, which can be directly compared to experimental data. This is regularly performed on the full 2D state output from SOLPS-ITER, however saving the necessary 2D state for many timesteps is not practical. The synthetic diagnostics should be directly evaluated in the main program loop, or by performing multiple code evaluations as is already done in the IPS feedback control implementation described above. Experimentally realistic actuator models including latency must also be developed and applied.
The sensitivity of the model to hyperparameters such as the input normalizations, smoothing windows for evaluatinġ x, and the threshold parameter λ must be assessed for general application to other devices and experimental conditions. In general, the hyperparameters should be scanned to find Paretooptimal models, which is facilitated by the relatively low computational cost of SINDy. Methods shown to be more robust to noise such as weak-form SINDy and kernel based model identifications are being evaluated and will be presented in the future. Spurious identifications of high order terms can prevent the integration of the ODEs to within a given error threshold over an acceptable prediction windows. Improved selection of state variables and the SINDy function library may improve the model sparsity. Another potential approach is to develop a 'trusted' nonlinear model that has trained on a large dataset including long time history of simulation and experimental data. Adaptive retraining would then only fit the nonzero terms in this model, possibly also using linear models and dynamic mode decomposition [65] to span ranges with poor performance of the nonlinear model.