Density control in ITER: an iterative learning control and robust control approach

Plasma density control for next generation tokamaks, such as ITER, is challenging because of multiple reasons. The response of the usual gas valve actuators in future, larger fusion devices, might be too slow for feedback control. Both pellet fuelling and the use of feedforward-based control may help to solve this problem. Also, tight density limits arise during ramp-up, due to operational limits related to divertor detachment and radiative collapses. As the number of shots available for controller tuning will be limited in ITER, in this paper, iterative learning control (ILC) is proposed to determine optimal feedforward actuator inputs based on tracking errors, obtained in previous shots. This control method can take the actuator and density limits into account and can deal with large actuator delays. However, a purely feedforward-based density control may not be sufficient due to the presence of disturbances and shot-to-shot differences. Therefore, robust control synthesis is used to construct a robustly stabilizing feedback controller. In simulations, it is shown that this combined controller strategy is able to achieve good tracking performance in the presence of shot-to-shot differences, tight constraints, and model mismatches.


PAPER • OPEN ACCESS
Density control in ITER: an iterative learning control and robust control approach

Introduction
The next generation experimental tokamak, ITER, is cur rently being built. This device is larger than its predecessors [1], which gives rise to multiple challenges, one of which is plasma density control. Density control is usually, in present day devices, achieved by feedback control of the gas valve openings based on the deviation between the requested and measured plasma density [2]. However, in ITER, the response of gas valves might be too slow for feedback control. Firstly, the valves are positioned far away from the vessel (i.e. approximately 20 m) delaying the response by the travel time through the pipe, being in the order of a second or more. Secondly, it is more difficult for neutral gas particles to penetrate the hot ITER edge plasma, having edge temperatures in the order of 1 keV, further reducing the effectiveness and response of gas fuelling [3,4]. Pellet injection is available as a second actuator which directly fires pellets of frozen fuel into the plasma, which penetrate further [5]. However, since these pellets are fast, they lead to localized density increase, and they may not yet be able to ablate in the early phase of the shot where the temperature and density are still low. Beside these two key actuators, other factors influence the L-mode density evo lution, including the vacuum vessel pumps and the presence of the plasma facing components. Optimal preparation for an ITER shot can only be achieved by accurate modelling of for example the ramp-up Plasma density control for next generation tokamaks, such as ITER, is challenging because of multiple reasons. The response of the usual gas valve actuators in future, larger fusion devices, might be too slow for feedback control. Both pellet fuelling and the use of feedforward-based control may help to solve this problem. Also, tight density limits arise during ramp-up, due to operational limits related to divertor detachment and radiative collapses. As the number of shots available for controller tuning will be limited in ITER, in this paper, iterative learning control (ILC) is proposed to determine optimal feedforward actuator inputs based on tracking errors, obtained in previous shots. This control method can take the actuator and density limits into account and can deal with large actuator delays. However, a purely feedforward-based density control may not be sufficient due to the presence of disturbances and shot-to-shot differences. Therefore, robust control synthesis is used to construct a robustly stabilizing feedback controller. In simulations, it is shown that this combined controller strategy is able to achieve good tracking performance in the presence of shot-to-shot differences, tight constraints, and model mismatches.
Keywords: plasma density control, ITER, robust control, feedforward control, pellet fuelling, gas fuelling, ramp-up (Some figures may appear in colour only in the online journal)

International Atomic Energy Agency
Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. phase. Nevertheless, such prep arations are only meaningful if the actual achieved density tracks the modelled one, avoiding too large density oscillations or deviations from the desired density evolution. As such, density control is especially relevant for larger devices like ITER.
The topic of this paper is the control of the particle density in ITER despite the issues mentioned above. This controller must be able to deal with the large time scale separation in actuators (pellets and gas), fast transitions and changes in dynamics of the plasma during the plasma ramp-up, and complications in the modelling of plasma fuelling. It will be shown that these complications make simple linear feedback control alone not suitable for ITER. Feedforward control is needed in addition to feedback control to resolve the aforementioned problems. A well-tuned feedforward actuator command is able to act effectively and is inherently non-causal, enabling corrections before they are needed, e.g. to inject gas in anticipation of a density change request. In present tokamaks, to achieve the needed control performance, the feedforward control signal is the result of meticulous tuning. For ITER, this is not advisable because of the much higher cost of a single shot compared to present day tokamaks. We therefore propose to use iterative learning control (ILC), a control method whereby the time trajectory of the actuator input signals is modified from preceding experiments (all having the same density reference) in such a way that the norm of the tracking error over the period of interest is reduced [6]. This can be achieved by using the result of one trial to design an improved feedforward signal for the next. ILC therefore particularly works well for repetitive control systems, which have the same reference and perform the same task over and over again. ILC has been successfully applied to for example industrial motion systems [7,8]. The use of ILC in tokamaks was proposed in [9], where its application to current density profile control is discussed. Furthermore, characteristics of tokamaks that make ILC particularly suited for tokamak plasma control are also highlighted. In particular, tokamaks are essentially repetitive systems as well, since they have to be operated in shots. Even though the shots differ due to different experiments, many parts of each shot (including the ramp-up) are often the same. Also, many quantities will be controlled in feedforward, in particular in the early operational stages of a new device due to lack of real-time measurements, or because feedback controllers have a longer development time. The approach can also be beneficial when restarting a device after a substantial upgrade. In ILC, actuator and operational constraints can be easily implemented, which is not the case for feedback control design. It is important to recognize, however, that even for two identically prepared tokamak shots, each shot will be slightly different due to for example different wall conditioning. Pure feedforward based control cannot deal with these non-repetitive disturbances. Feedback control is therefore necessary to deal with the possible loss in reference tracking performance and can work together with any type of feedforward, including ILC. To guarantee stability and feedback control performance during the whole rampup, the advanced robust H ∞ synthesis technique is used to synthesize the feedback controller [10]. With this technique, feedback controllers can be mathematically synthesized for systems with uncertainties. In figure 1 the proposed controller architecture is visualized. Here, the block indicated by Σ denotes the tokamak, r the reference signal of the average density, which is taken from a DINA simulation of the ITER ramp-up [11], and y 1 the average density that is achieved. The difference between this latter density and the reference is denoted by e 1 and is fed to the feedback controller. The actuator signal u 1 is the sum of the feedback control signal u fb and the feedforward control signal u ff . Clearly, when only feedforward control is used, u fb = 0.
In this paper we will apply the mentioned controller structure to simulations of the particle density evolution expected in ITER, and show that this control methodology is able to resolve the aforementioned problems. To test our controller, we use a control-oriented model of the plasma particle density evolution [12], which has been specially updated to include issues that govern the particle density evolution in ITER during ramp-up.
It is important to realize that the objective of this paper is not to make quantitative statements about the physics of the density evolution in ITER and the possibility for control. Much of the details of wall retention, scrape-off-layer (SOL) and fuelling physics for ITER are still uncertain at the time of writing. Instead, the purpose of this paper is to show that, based on a control-oriented model for the ITER density evolution, based on today's best guess of ITER density dynamics, the proposed control scheme can successfully control the ITER density during plasma ramp-up.
The remainder of this paper is structured as follows. Firstly, in section 2, a brief summary of the proposed ITER current ramp-up phase is given, describing its various phases, targets and constraints. This is followed by a brief summary of the controloriented model in section 3. The level of complexity in this model is aimed to be sufficient to capture all relevant fuelling dynamics needed to properly assess ITER density control during the rampup phase. Then the controller design is discussed and justified in sections 5 and 6, after which simulation results showing the performance of the controller are presented in section 7.

DINA simulations of plasma ramp-up phase
During the ramp-up phase, the plasma current, I p , and density, n, are increased, preparing the plasma to be heated to the temperatures needed to achieve thermo-nuclear fusion. Density control is critical to consistently achieve a reliable and stable ramp-up phase, since variations in the density cause variations in the temperature evolution, which in turn cause a different current density profile evolution, which may have important consequences for the confinement quality and MHD stability. Towards the end of the current ramp-up, a minimum density should reliably be achieved that allows for example the use of neutral beam injection (NBI) heating, and the triggering of the high-confinement mode (H-mode), required for achieving the ITER burn conditions. DINA simulations have been performed to model the plasma evolution during the ramp-up [11]. The results of those simulations will be used as starting point to tune a controloriented model. The time-evolution of key parameters of the simulation are shown in figure 2. The following features are relevant: at approximately 17s, the plasma configuration is diverted, and the end of the ramp up coincides with the L-H transition. At this point, having I p = 15 MA, a density of 4 · 10 19 electrons/m 3 is achieved. Note that the DINA code does not self-consistently solve the density profile evolution. In the remainder of this paper the density will be scaled such that n avg,20 = n avg /(1 × 10 20 (m −3 )) particles.
Although the DINA model of the ramp-up assumes a steady increase of the density (together with the plasma current), this assumption may not be easy to achieve. Firstly, as we will see in this paper, it will depend on the capabilities of the density control system. Secondly, a number of physics effects may cause perturbations. For example, right after breakdown the cold and low density plasma may be influenced strongly by the wall, which could act as a particle pump, yielding a drop in density [13]. Furthermore, the change from a limited to diverted configuration will alter the pumping of the plasma. Finally, as the (edge) temperature increases during the rampup, it may be more and more difficult for particles to be ionized deep enough into the plasma to be effective for fuelling. This effect, depending on the details of the transport of neutrals and ions in the plasma edge, may influence the fueling of ITER plasmas significantly. Therefore, in the remainder of this paper, we aim to include these complications in the controller design and testing.

Density limits
In this section, the plasma density limits during ramp-up will be specified. As this topic is still being investigated, the limits shown here are solely for the purpose of demonstration, and do not necessarily reflect the expectations for ITER. A lower limit on the density is set by the possibility of error field penetration. This low density limit is determined by the level of error fields due to the slight misalignment in the field coils, which may only be revealed after construction, but also depend on the success of the correction coils that are going to be used at ITER [14]. Hence, for ITER this limit is not precisely known yet. During limiter plasma operation, an upper limit for edge fuelling is set by a radiative collapse. When the plasma is diverted, an upper limit is set by the divertor detachment, which is to be avoided during initial operation. Both upper limits are closely linked to the density in the scrapeoff-layer (SOL) or far outer layer of the plasma [15]. During the ramp-up phase these upper limits are usually lower than the Greenwald density, n GW , the well-known upper limit for densities in tokamaks [16]. For simplicity the lower limit is chosen to be 0.2n GW . The Greenwald density limit is only a function of plasma current I p (MW), and is given by where a is the tokamak minor radius in meters.
The precise physics mechanisms of upper limits on the plasma density is still being investigated [17,18]. For now, in this paper, we assume the upper limit is a simple fraction of the Greenwald density, being n up = 0.45n GW throughout the ramp-up. It can be refined in the future as the physics behind this constraint are better understood. It is generally acknowledged that these physical density limits do not hold for pellet injection [5]. Therefore, in this work, we allow minor density violations due to pellets. In figure 3 the reference plasma density and both density limits are visualised.

Description of the control-oriented model
To assess the density control during the ramp-up a dynamic model of the plasma density evolution in response to various fueling actuators is needed. The model should capture the key aspects that will affect the evolution of the density. However, for a preliminary control study it is not necessary to model many of these physics processes in great detail, but it is more important for iterative controller design that the model is fast enough for rapid simulation (a few minutes per run at most). Furthermore, the model structure should be suitable for controller design. Here, a recently developed three inventory model approach [19] containing a plasma, wall, and vacuum inventory is used. This model was originally developed to improve the density control in TCV and ASDEX Upgrade [12]. For a detailed description of the model, including the model equations, we refer the reader to [12]. For the work presented in this paper, the model was updated to capture elements relevant specifically for ramp-up. For example, taking into account the fuelling of the plasma after breakdown, the possible impact of wall pumping and the creation of a diverted plasma during the ramp-up. The loss of fuelling efficiency with the increase of edge temperature, has also been included. The model is heuristic and has tunable parameters that can be adjusted to yield the density evolution one expects in ITER. However, since there is uncertainty in modeling the mentioned effects, it is difficult to decide which values to choose for these parameters. Therefore, many parameters can only be guessed based on experience, or extrapolated from data available for existing machines. Inevitably, there will be significant uncertainty in these parameters. The control methodology that will be introduced in section 5 and section 6 is therefore designed to guarantee performance in the presence of uncertainties, by taking them into account in the design of the controller.
In figure 4, schematic outline of the model with three inventories is given, including fluxes between each particle reservoir. There are three inventories that can contain particles: the plasma, the vacuum vessel, and the wall. The main sources and sinks to the entire model are: the neutral gas fuelling by means of gas valves, fuelling with pellets and the pumping of the vacuum.

Inventory description
The interaction between the three inventories can be described by a set of coupled differential equations. The inventories interact by exchanging particles via various physics processes, as shown in figure 4. Neutrals in the vacuum can be ionized and enter the plasma (ionization), ions in the plasma can recombine, neutralize and exit the plasma to the vacuum. Neutrals can move from the outer plasma, the scrape-off layer (SOL) either to the wall (SOL-wall) or the vacuum (SOLe flux) and from the wall they can recycle back into the vacuum (recycling). The vacuum and wall are regarded as zero dimensional, but the plasma is described by a 1D partial differential equation for the electron density profile. The details on how these processes are modelled can be found in [12]. A summary of the parameter values is given in the appendix. Two extensions relevant for ITER or the plasma ramp-up are discussed next.

Neutral screening effect
The decrease in effectiveness of gas fuelling with increasing edge temperature is modelled in an ad-hoc fashion. An ionization deposition function Λ iz (ρ) (as a function of normalized toroidal flux label ρ) is chosen to model the distribution of the neutral density decay in the plasma and the temperatureprofile dependent reaction rate. It changes in time as a function of the ionization depth λ iz . Here, it is assumed that λ iz scales with edge temperature variations T e,b − T e,b0 around a nominal edge electron temperature T e,b0 : For now, the proportionality k iz,T is chosen to be −1 · 10 −4 eV −1 . The mean ionization depth λ iz,0 is introduced in section 7.

Plasmaless model
Gas puffing prior to plasma breakdown is modelled by only taking into account the wall and vacuum inventories. Wall outgassing, gas puffing, and pumping can affect these inventories. These processes are modelled by removing the plasma dependent terms in the relevant inventory equations [12].

Model of the gas supply system
For ITER, the gas supply system consists of 21 m of pipe between the valve and the vacuum vessel, with mass flow controlled (MFC) valves. The gas supply is limited to 7.36 × 10 22 particles/s, corresponding to 200 Pa m 3 s −1 . The MFC system is modelled by a simple first order dynamical system with a constant time delay of 0.5 seconds and rise time of 0.4 seconds. Its transfer function P g (s) is given by:

Pellet fuelling model
The pellet deposition is modelled as depositing particles of fixed size on a fixed location in the plasma. For simulations, the pellet size is set to 6 × 10 21 , corresponding to the theoretical size [20]. Furthermore, the pellet system has a constant delay of 0.2 s. The pellet model has a maximum delivery frequency, but we assume a 100 percent reliability in the delivery of pellets by the launcher. The ITER pellet injection system cannot use the pellet velocity as control parameter to optimize the pellet ablation. The pellet size, however, can be changed by 20% on a time scale of roughly 3 s. Since this work presents a first control assessment, this is not included in the pellet model.

Model input and outputs and discretization
Galerkin projection yields an ordinary differential approximation of the partial differential equation system described in [12]. The vector [Γ valve Γ pellet ] T acts as inputs u to the model, having n u entries. The only controlled output of the system is the volume averaged density, which is a linear combination of the internal states [12]. As the underlying model is 1-dimensional, possible future outputs like the separatrix density (necessary for more elaborate density constraints) can be easily implemented. For digital implementation and computer simulation, the model is discretized in time on a time grid having n t grid points. In the following, discrete representations of inputs and time are denoted by subscript k.

Simple control applied to ITER density tracking
As anticipated in the introduction, simple linear (e.g. PID) controllers will not suffice for density control in ITER. In this section, we illustrate the shortcomings of these controllers in control-oriented simulations of ITER ramp-up using the model presented in the previous section. In figures 5 and 6 two simulations are shown with two different controllers: one PD controller and one PID controller with anti wind-up. Anti wind-up is a technique that prevents actuator saturation due to the integral action of the controller. Both controllers aim to control the density evolution using solely gas injection. The input of the controller is the tracking error, which was defined in figure 1. The controller needs to reduce this error, and by doing so the tracking performance is increased. Clearly, the time delay results in insufficient tracking performance a few seconds after breakdown for both controllers. The PID-controller performs the best early in the ramp-up, but has the largest error at the end even though anti windup is applied. Both controllers cause density oscillation in the early ramp-up. This performance reduction of a controller in the presence of delay is a common effect [10]. Reducing the controller gain helps to resolve these oscillations, but in this case, lower gains do not achieve sufficient performance, and violate the lower density limit. Additional pellet fuelling can only help to resolve this issue after a sufficiently high density and temperature is achieved, which is too late to prevent density limit violations in the early ramp-up.

Introduction to ILC
Iterative learning control provides a means to increase performance for complicated, repeated tasks [6]. In figure 7, the basic structure of an iterative learning control scheme is visualized. In this figure, a standard feedback interconnection of a plant Σ and controller K can be seen. The system Σ denotes the true system, i.e. the real tokamak, see figure 1. Upon carrying out the first experimental trial, this results in an error signal e 1 (t) in the time domain. Then, a new feedforward signal for the next trial (2) is designed off-line based on this error. This is Figure 5. Tracking evolution for a PID controller with anti-windup, and a PD without. Clearly, performance is insufficient for both controllers. The delay of the gas supply system limits the achievable performance, yielding either oscillations or a too slow response. The integral action of the PID controller results in a large error at the end of the ramp-up, even though anti-windup is applied. done using a quadratic programming (QP) algorithm, as will be explained later. In this paper, we will use so-called optimal ILC, where the change in each new feedforward signal is the optimizer of a (constrained) optimization problem [9]: Here, the vector u is the complete trajectory of the input introduced in section 3, whereas ∆ u denotes the deviation from the previous or predefined trajectory for u. The cost function is denoted by J, whereas A ineq and b ineq describe the inequality constraints on the actuator signals. Matrix A eq and vector b eq contain the equality constraints. In the coming section, all elements in this constrained optimization problem formulation will be derived.

Linear time varying model
To apply Optimal ILC, a reasonably accurate linear model is required. During plasma ramp-up the system dynamics change, for example due to increasing plasma current and decreased transport. Therefore, a simple linearization about one operational point would not be able to capture these physics sufficiently. Here, we will use the linearization of the previously derived density model along a suitable trajectory. The result will be a linear time varying model (LTV model) [21], which can be used to predict the variation in the output of the system to a change in input signal. For the nominal trajectory along which the nonlinear model will be linearized we choose the desired density evolution. Then, perturbations δ of the output y k , state x k , and input u k around this trajectory on time k can be used to write: Note that the system matrices also depend on k. Next, all inputs on all N time points can be stacked in a single vector, denoted by δ u : and in the same fashion for δ u , δ y, and δ x . Then, a state transition matrix T ∈ R Nny×Nnu , is used to write all outputs as a function of all inputs: where D ∈ R Nny×nx is another matrix, describing the effect of the initial condition x 0 . The matrices T and D can be derived by recursively writing the output of the system [9], and can also incorporate descriptions of the actuator delays. With this LTV description, we can write for the error of the next trial Substituting, we can derive an expression for the predicted error at the next trial as a function of the change in input ∆u = δ u j+1 − δ u j : Note that e j is defined as the error achieved using the true system output and not the model. Note the difference between 'trials' j, which are (simulations of) tokamak discharges, and time points k, which are elements in the time vector within a (simulated) discharge.

Cost function derivation
The next ILC trial should aim to reduce the tracking error, while simultainously avoiding large changes or oscillations in the input trajectories. Such objectives can be achieved by designing a suitable cost function. This cost function contains the output error, which is smoothed off-line using a fourth order Butterworth filter [22] to discard oscillations in the signal due to pellet injection. It is assumed that after a shot, the truly achieved n avg can be recovered from measurements. Besides this smoothed error, the cost function also penalizes the input itself. Furthermore, changes in the input trajectory are penalized, such that the final converged input will be close to the original one. Otherwise the linearization might be a poor approximation. Each cost function term is weighted by a weight factor, such that the relative importance of cost function terms can be tuned: The error weighting matrix W e is typically given by W e = ν e I + ν s I s , where ν e denotes the weight factor for the error. By structuring it this way, the error on each time step is penalized equally. A second weight factor and diagonal matrix, denoted by ν s and I s respectively, are introduced to only penalize the error on certain time instances, for example the L-H transition. The penalty on the control input effort is structured as: where ν ... denote the different penalties on use of gas and pellets. The operator ⊗ denotes the Kronecker product, and I ∈ R nuN×nuN . The very same structure is used for W ∆ u , but with weights ν d,gas and ν d,pel . The ILC weight factors are normalized on the error: ν e = 1. The tuning procedure is as follows: first all weights are set to zero except ν d,gas , ν d,pel , and ν e . As such, a trade-off between performance and numerical conditioning of the optimization problem can be obtained. Then, the other weights are tuned one by one to achieve their target. All weights are unchanged throughout the rest of the paper, and summarized in table 1. The error, e j+1 can be written in terms of δ u using (9). The cost function can be written in the form A∆ u − b by defining: By gathering the remainder of the cost function terms (the ones that do not explicitly depend on ∆ u) the following b can be derived: Next, the expression A∆ u − b 2 2 is expanded to arrive at: where H, f, and c denote the standard notation for quadratic optimization problems [23].

Constraints
The inequality constraints consist of actuator limits, the lower density limit based on the Greenwald density limit, and the upper limits, as discussed in section 2.1. The actuator limits can simply be written as: Where the actuator upper limits are denoted by u max . The density constraints use a prediction by the LTV model and actually achieved density trajectories that are recovered after a shot. For example, for the upper limits we can write: Note that as such, these constraints are also formulated in ∆ u. Combining all constraints, the following inequality can be derived:  Where I denote identity matrices of the appropriate sizes, and n g,l denotes the lower Greenwald limit. Next, we implement one equality constraint. This constraint prevents the use of pellets early in the ramp-up, which was discussed in section 1. This constraint simply consists of a matrix that select certain elements of the input vector, such that the following can be written: since if ∆ u is forced to be zero, and the nominal simulation does not use pellet inputs during the required time period, they will be constrained to zero for all consecutive shots.

Robust feedback control
The ILC scheme outlined above generates a feedforward signal, but can not cope with shot-to-shot differences nor with unexpected disturbances acting on the plasma. For this reason, it is beneficial to include a feedback controller to compensate for this [10]. A feedback controller is designed using robust  control (RC) theory [10], as a first implementation only for gas injection. Pellets are therefore only controlled in feedforward. The robust controller is designed simultaneously for four different snap-shots of the LTV system along the complete ramp-up using an LMI (linear matrix inequality) formulation of the H ∞ design problem [24]. It is assumed that if the controller robustly stabilizes these four systems, it will stabilize all systems in the LTV model. In figure 8 the Bode diagram of the open loop transfer KΣ for two LTV models in the ramp-up is shown: The beginning of the ramp-up for Σ 1 , and the very end for Σ 4 . Also, the effect of the delay on system 1 is shown as Σ 1,del . Clearly, the delay causes significant decrease of the phase margin. As discussed before, the large gas supply delay could result in destabilization of the feedback control system [10]. Therefore, the delay margin [25] is checked for all systems in the LTV model. The robust controller has a delay margin larger than 3 seconds for all systems in the LTV model. Therefore, stability of the closed loop during the whole ramp-up on the real tokamak system can be assumed.

Introduction and proof-of-principle simulation
In this section, simulation results are presented showing the performance of the previously derived controllers. For this purpose, we define the 'true model' as the model (section 3) with parameters that are assumed to corresponding to the real system. These parameters are not known for the control design. The control design uses instead a 'control model' representing our best knowledge of the true model. The true model has different pumping dynamics, and different values for the dimensionless mean ionization depth λ iz,0 , and the neutral screening proportionality k iz,T , as they are the main source of model uncertainty and significantly influence the dynamics. In table 2 the parameter variations are summarized. For the simulation, first, the control model is used to manually tune a feedforward control that results in sufficient reference tracking. In figure 9 it can be seen that the same feedforward, applied to the true model, yields a different average density evolution (in blue) that violates the constraints (in grey). Then, pure ILC, i.e. without feedback control, is applied for 10 consecutive ramp-ups, iteratively updating the gas valve and pellet time waveforms. The resulting trajectory is shown in figure 10. Indeed, the density evolution has converged to the desired reference trajectory. Note that pellets cause slight momentarily violations of the density constraints, which is deemed acceptable due to their rapid diffusion: The algorithm designs and assumes a constant inflow of particles due to pellets, which is only later divided into separate pellets in the simulation. The time traces of gas and pellet inputs that were necessary to achieve the improved density tracking are shown in figure 11. Note that before t = 10 s the pellet inputs are constrained to zero, as was discussed in section 5.3. Furthermore, gas is injected before t = 0 s to achieve sufficient tracking performance just after breakdown. The error two-norm evolution is plotted in figure 12. Because of the identical model mismatch for all ramp-ups, the convergence is smooth. The corresponding error in the time domain for all odd trials is plotted in figure 13. Around t = 18 s density control is difficult. Here,   pellets are just allowed, but are relatively large due to their fixed size. Furthermore, the density limits are particularly restrictive in this part of the ramp-up. Therefore, selecting the proper fuelling inputs is critical and ILC particularly shows its strength.

Combined control to handle shot-to-shot differences
Next, we simulate the effect of shot-to-shot differences.
Besides the previously applied constant model mismatch, each ramp-up now also has a different value of λ iz,0 and τ SOL in a range of two percent around the nominal value of the true system, which is summarized in table 2. In this case, the robust feedback controller is needed to achieve a consistent ramp-up tracking performance. The simulated shot-to-shot differences are significantly smaller than the constant control model mismatches, which is expected to be the case during real operation. The nominal mismatch is now constructed in such a way that the true tokamak system fuels easier than the model. Therefore, using the initial feedforward waveform, the achieved average density lies above the reference, and in fact above the constraint. This is indicated in figure 14. The first five ramp-ups are solely ILC based, in order to relax the effort of the feedback controller. The weight factors in the ILC cost function are unchanged with respect to the proof-of-principle simulation before. Then, 5 trials using both controllers are performed. The resulting density evolution is plotted in figure 15.
Despite the shot-to-shot differences, tracking performance is greatly improved. Again, pellets cause slight violation of the density limit.  Starting point of the simulations including shot-to-shot differences. Here, the model mismatch is such that the real system fuels more efficiently, resulting in a higher density for a given input compared to the model.  . Optimal inputs as determined by the ILC algorithm. The initial (dashed blue) and final gas input Γ valve (t) (solid blue), and the initial (dashed amber) and final pellet input Γ pellet (t)(solid amber) are plotted. Due to penalty on changes in input, the trajectories remain close to the original ones. Before t = 10 s, the pellets are constrained to zero. Gas is injected before t = 0 s to ensure sufficient tracking of the reference after breakdown.
For this simulation, convergence is not monotonic due to shot-to-shot differences. This is indicated in figure 17. The changes in inputs that had to be made to achieve the increase in performance are indicated in figure 16. Note that in this plot, only the gas feedfoward trace (i.e. without the feedback controller) is shown. Furthermore, weights are all relative to one another: As gas is becoming less efficient and the error is increasing, the algorithm will decide to use pellets over gas, since the weight on the error is higher than on the use of gas over pellets. The relative effort of feedback and feedforward control is summarized in figure 18. Here, the feedback controller clearly acts on the disturbances arising due to pellet injection, which is not seen by the ILC scheme due to the Butterworth filtering. Nevertheless, overall, ILC does most of the control.

Statistical analysis of shot-to-shot differences
Next, fifteen sets of ten ramp-ups are simulated to investigate convergence of the error 2-norm in the presence of shot-toshot differences. Again, τ SOL and λ iz,0 are varied in a range of two percent around their nominal value for the true system   Error convergence for fifteen sets of 10 ramp-ups in the case of a less effective fuelling for the true system. The shot-to-shot differences clearly have influence on the converged error, and the mean value is significantly larger than for the case without shot-toshot differences. Figure 15. Performance of ILC and RC combined after 10 simulated shots including simulated shot-to-shot differences in τ SOL and λ iz,0 . Tracking is greatly improved, and only the transient pellets cause density violation. Figure 16. Optimal inputs as determined by the ILC algorithm. The initial (dashed blue) and final gas input (solid blue), and the initial (dashed amber) and final pellet input (solid amber) are plotted. Note that for the gas input, this is solely the feedforward control action.
for each new ramp-up. First, the nominal values for the true system are chosen such that the fuelling is less effective. Then, the error 2-norm convergence of the error in figure 19 is achieved. For the case without shot-to-shot differences, the error-norm is slightly smaller, as is indicated in amber in the same figure. The shot-to-shot differences only have a small influence on the converged error in this case. The same holds for a positive model mismatch, where the true system has a more effective fuelling. This is indicated in figure 20.

Outlook
The work presented in this paper contains a systematic procedure for density control design using iterative learning control (ILC) for feedforward improvement and robust control (RC) for feedback control, that can handle the presently known constraints for actuators and physics. When the date of the first plasma in ITER approaches, new understanding of fuelling and particle transport mechanisms will help refine the model, while detailed knowledge of the actuator capabilities and other constraints will better define the boundary conditions for density control. With this new knowledge, the controller design based on the ILC+RC approach can be updated to yield predictions for the ITER density control.
Already today, more work can still be done on both the model and the control itself. The sensitivity of the controller to the different (heuristic) parameters in the model could be verified more extensively using parameter scans. Experimental data from JET or other tokamaks could be used to benchmark the model, the model can be included in the ITER plasma control system simulation Platform (PCSSP) [26] for systematic closed-loop testing in conjunction with other controllers active during the plasma ramp-up.
On the control side, one could tune ILC for multiple scenarios. Here, we only used a standard DINA simulation rampup scenario, but one could also use different cases. This work considered only the plasma ramp-up until the L-H transition. Future studies could focus on density control following the L-H transition and into the flat-top phase, where slow timescale feedback control is expected to play a more important role. The ILC approach lends itself well to other transient phases of the plasma, in particular the ramp-down, which could also be studied. Also, many more constraints can still be implemented. For example, the loss-of-fuelling effect could be severe, and one may need pellets in the feedback controller to deal with the destabilizing effects of the shot-to-shot differences. However, the neutral screening effect and the large gas transport delay complicate the choice between the two inputs for the controller. Also, note that in this paper the focus lies on n avg -control. However, in principle, it is also possible to choose a different parameter to control, such as the SOL or edge pedestal density. The controlled variable could also vary in time.
Finally, controller design is just one part of the control loop. In reality, the volume-averaged density on which the robust controller operates may not be known reliably at any point in time, due to diagnostics errors. Dynamic state reconstruction algorithms such as a Kalman filter could be used to estimate the average density, or even the density profile from measurements, as recently applied to TCV and AUG [12].

Conclusions
In this paper, a control-oriented transport model for the plasma density evolution in ITER is used to design a combined feedback and feedforward control solution for density control during the ITER ramp-up. The model from [12] was adapted for use in ITER, and it was shown that the model can reproduce the density evolution in the reference ramp-up scenario from DINA simulations. It has been shown that it is able to reproduce the expected changes in pumping, fuelling efficiency, and the severe gas delay during the ramp-up. It has been shown that simple proportional feedback control lacks performance for the ITER ramp-up scenarios. ITER density control needs advanced controllers mainly because of two reasons. First, the delay of the gas valve is so large that feedback alone inherently fails to achieve enough performance after breakdown. Second, feedforward is in theory able to achieve this performance, but is difficult to tune for a new device like ITER. Therefore, a self-learning algorithm based on iterative learning control (ILC) [6] has been implemented. This technique has been shown to resolve delay problems and achieve sufficient reference tracking for deterministic ramp-up scenarios. Unfortunately, on a real tokamak, not every ramp-up is the same. Therefore, a robustly stabilizing feedback controller has been developed that can deal with the disturbances due to small shot-to-shot differences. It has been shown that together with ILC, this robust controller can achieve convergence in the tracking error even though some model parameters are changed randomly, and is still able to track the desired reference trajectory after a small number of shots. It can therefore be concluded that, based on simulation results, iterative learning control in combination with robust control is a promising solution for ITER density ramp-up control. Error convergence for fifteen sets of 10 ramp-ups in the case of a more effective fuelling for the true system. In this case, the shot-to-shot differences have a slightly smaller influence on the converged error compared to the case in figure 19.

Acknowledgments
We would like to thank Dr. A. Polevoi for providing the DINA ITER ramp-up model. The first author is grateful for the many discussions on various plasma physics and tokamak operation topics during his stay at the ITER Organization that helped form this work. Also, the input of scientists from CEA Cadarache, especially Dr. Hugo Bufferand, on divertor detachment modelling which led to the implementation of the density constraint is kindly acknowledged. Furthermore, valuable remarks from Dr. Tom Oomen, Dr. Roland Toth, and Jurre Hanema from Eindhoven University of Technology on ILC and RC were gratefully received. This work was carried out with support from a Fusenet internship grant, and is supported in part by Dr. F. Felici's NWO (Dutch Science Foundation) VENI grant. ITER is the Nuclear Facility INB no. 174. The views and opinions expressed in this paper do not necessarily reflect those of the ITER Organization.

Appendix. Parameter values in the nominal model
In table A1 below the parameter values in the model are summarized. The mean ionization depth λ iz,0 , the pellet deposition location ρ pellet,dep , and the width of their deposition function ρ pellet,width are dimensionless like ρ, the normalized toroidal flux label. This label is defined via ρ tor = Φ/(πB 0 ) as ρ = ρ tor /ρ tor,B , where ρ tor,B is the toroidal flux label on the LCFS. More details on the equations where the parameters occur can be found in [12].