Inter-discharge optimization for fast, reliable access to ASDEX Upgrade advanced tokamak scenario

Rapid inter-discharge simulation and optimization using the RAPTOR code have allowed the development of a reliable and reproducible early heating strategy for an advanced tokamak (AT) scenario on ASDEX Upgrade. Solving for electron heat and current diffusion in RAPTOR with ad-hoc formulas for heat transport and electron cyclotron current drive (ECCD) efficiency is found to robustly recover the coupled dynamics of Te and q, while maintaining model parameters fixed for all discharges. The pedestal top boundary condition in pre-shot simulations is set by a newly derived scaling law for the electron pressure at ρ = 0.8, using a data set of previous AT discharges. RAPTOR simulations have allowed to develop an understanding of the onset of 3/2 tearing modes, which were observed to have a detrimental impact on confinement when low magnetic shear conditions are present at the rational surface during the high-β phase. Delaying the NBI heating, by a specific time interval found via simulations, has led to avoiding these modes. A non-linear optimization scheme has been applied to optimize the ECCD deposition radii to reach a stationary state with qmin>1 at the beginning of the flat-top phase, while ensuring a non-zero magnetic shear at q = 1.5 throughout the high-β phase, and has been successfully tested in experiment. However, further experiments, aiming for qmin>1.5 , have highlighted limitations of the present feedforward control approach in the presence of shot-to-shot variations that are not included in the applied model. Application of real-time model-based control is proposed to overcome model-reality mismatches in future work.


Introduction
Advanced tokamak (AT) scenarios aim for an elevated q profile, maximizing the bootstrap current fraction ( f bs = I bs /I p , note: j bs ∼ q∇p), to explore the potential of a steady state nuclear fusion reactor.Furthermore, these scenarios often achieve improved confinement relative to the IPB98(y,2) scaling law [1], further enhancing the bootstrap current fraction.On ASDEX Upgrade (AUG), an advanced scenario has been developed that applies counter-ECCD near the magnetic axis (where a high current drive efficiency can be achieved), to form an off-axis peak in the total current density.This strategy allows to study an advanced scenario with plasma parameters approaching reactor-relevant conditions.
During the plasma current ramp-up, the inductive current in the plasma is increased by using the central solenoid to apply a positive loop voltage at the edge of the plasma.The radial variation of the loop voltage profile U pl (ρ) throughout the plasma then provides a driving force for penetration of the ohmic current density.Heating the plasma during the plasma current ramp-up phase allows to reduce the resistivity of the plasma at an early stage, hence slowing down the current diffusion towards the center of the plasma.This early heating strategy slows down the tendency of an ohmic plasma to evolve towards a peaked current density profile (in stationary state, a flat loop voltage profile would correspond to a current density profile with a radial dependence self-similar to T 3/2 e ), prolonging the broad current density profile that corresponds to a more elevated q profile.This approach is of particular interest for advanced scenarios that aim to achieve a stationary, elevated q profile.Once the plasma current reaches the flat-top value, the elevated q profile can be maintained in a stationary state, if the relaxed ohmic current density profile is complemented with off-axis auxiliary current drive sources and a broad bootstrap current density profile.In H-mode plasmas, the bootstrap current driven in the pedestal region provides an important contribution.
On present-day tokamaks, accessing a stationary, elevated q profile early in the discharge allows for a longer time window to study the physics of the advanced scenario, lasting over multiple current diffusion time scales.Studying the stationary phase is important as stationarity is required for a reactorrelevant operating point.Early heating scenarios are of even greater importance for advanced scenarios on future machines like ITER and DEMO, where the current diffusion time scale is much larger due to the bigger size and the higher temperatures of the plasma.A late heating strategy, allowing the plasma to evolve to an inductive q profile before auxiliary current drive sources are applied to elevate the q profile, would result in a very long waiting time before the plasma reaches the stationary state conditions envisioned for the burning plasma state, hence compromising the commercial prospects of a fusion reactor based on an advanced scenario.
Disadvantages of the early heating strategy are a reduced robustness with respect to the onset of tearing modes, the impact of poorly known initial conditions and the significant sensitivity to the onset timing of heating and fueling actuators.As a consequence, finding the correct timing of the various actuators during ramp-up and the early flat-top phase gives rise to a delicate balancing act.For the JET hybrid scenario, the dependence of discharge stability on small scenario modifications has been reported in [2,3].Stationary states with different stationary current density profiles and differing confinement quality can be accessed depending on the ramp-up scenario (even for identical actuator inputs during the stationary state) [4].The different confinement quality is thought to be caused by the presence of NTMs and related differences in the q profile.Analysis of the reconstructed q profiles in [4] suggests that the tearing modes appear as the q profile drops through the corresponding rational surface, with the low local value of the magnetic shear facilitating the NTM onset.In [5], a scenario was developed which includes an ohmic L-mode phase before the high-performance phase, to let the current profile relax and to get rid of the impact of initial conditions.Afterwards, the q profile is modified in a more controllable way, as the dynamics are less sensitive to actuator onset timings when starting from a stationary rather than a transient state.Note that the importance of optimized actuator timings is not limited to scenarios with heating during the ramp-up phase.On DIII-D, the mostly inductive ITER baseline scenario is heated only after the flat-top phase is reached.The onset of a disrupting 2/1 tearing mode has been correlated to a specific feature of the current density profile, namely the formation of a steep well in the j par profile around the q = 2 surface [6].The stability of the disruptive 2/1 modes was found to be significantly dependent on the early current diffusion evolution.By manipulating plasma current, heating and density traces, a recipe avoiding the modes could be developed [7].
While the plasma dynamics are highly sensitive to actuator time traces during transient phases, profile simulators allow to evaluate the deterministic impact of actuator modifications.RAPTOR is a lightweight, modular, 1-dimensional transport solver, allowing to solve for the time evolution of poloidal flux ψ, electron and/or ion temperature T e,i and density n e,i [8,9].A set of reduced physics models and an efficient numerical implementation enable the high computational speed required for applications in real-time control and non-linear scenario optimization.
The application of RAPTOR for ramp-up optimization was first proposed in [10], and applied to optimize access to the hybrid ITER scenario in [11].Optimization of the stationary operating point of the ITER hybrid scenario has been performed in [12].A framework for ramp-down optimizations was first described in [13], and later applied to study optimal termination scenarios for AUG [14] and DEMO [15].In [16], RAPTOR, coupled to the QLKNN-hyper-10D transport model [17] and ADAS cooling factor data to evaluate impurity line radiation [18], was used to optimize the rampup phase for WEST.By successfully modeling how increased Nitrogen injection can peak the current density and enhance the core temperature, a ramp-up with improved MHD stability and increased margin to tungsten contamination has been proposed.This example illustrates the strength of a modular transport code, where one can select the models required to describe the physics at play.
The aim of the model-based scenario development and optimization in this paper is to make the counter-ECCD AUG early heating scenario more robust with respect to tearing modes, while reaching the desired, elevated q profile as early as possible in the discharge, at the beginning of the flattop phase.Reaching a reactor-relevant operating point in a time-efficient way, this work contributes to the investigation of the feasibility of the advanced scenario tokamak concept.Furthermore, the capability of RAPTOR to predict relatively complex features of the profile dynamics is validated to experimental data, increasing confidence regarding its application for simulation of ITER and DEMO.As will be shown, the use of trajectory optimization helps to achieve the desired q profile evolution, reaching q min < 1.5 early in the discharge, while maintaining q > 1 and ensuring non-zero magnetic shear near rational surfaces like 3/2 and 2/1.We will also quantify the limitations of feedforward trajectory optimization in the presence of shot-to-shot variations not captured by the model.For future experiments, real-time, model-based control is proposed to make the scenario more robust to model-reality mismatches.
The remainder of this paper is structured as follows.Section 2 discusses the RAPTOR set-up used for pre-and postshot simulations of AUG discharges, as applied for an extensive set of discharges performed within the advanced scenario program and the ITER baseline program.Furthermore, the origin of the experimental data used for benchmarking is presented.In section 3, we describe how RAPTOR simulations and optimizations have been leveraged to improve our understanding of the dynamics and performance of the AUG 1 MA counter-ECCD advanced scenario, in particular regarding the onset of tearing modes.In section 3.1 we discuss the characteristics of the desired stationary operating point.We then show how inter-discharge simulations and optimizations in RAPTOR have allowed to develop a reliable and reproducible early heating scenario, reaching a stationary state early in the flat-top phase.Section 3.2 presents how post-shot simulations led to a hypothesis on the origin of a 3/2 tearing mode observed in a previous early heating attempt.Based on this hypothesis, characteristics are defined, constraining the time evolution of a stable approach strategy towards the desired stationary operating point.Predictive, pre-shot simulations have been used to design a ramp-up strategy avoiding the onset of the deleterious 3/2 mode, as discussed in section 3.3.We show how this strategy has been successfully tested in experiment.Optimization of the EC deposition radii is discussed in section 4, aiming for relaxation of the q profile immediately after the flat-top phase has been reached.The optimization problem is formulated mathematically in section 4.1.In section 4.2, the aim is a stationary q profile with 1.35 < q min < 1.5, applying the strategy developed in section 3.3 to avoid the 3/2 mode.Section 4.3 attempts to avoid the creation of a q = 1.5 surface by pursuing q min > 1.5.In the section 5 we illustrate how RAPTOR simulations can be leveraged to gain insight in the sensitivity of the optimum and the physics characteristics of the operating point.Finally, an outlook towards real-time solutions and overall conclusions are presented in sections 6 and 7 respectively.

RAPTOR workflow for inter-discharge simulations for AUG
For all the AUG simulations presented in this paper, RAPTOR solves for both electron heat and current density transport, while T i , n e and Z eff are imposed from experimental measurements.The fact that the ion temperature has a less significant impact on the q profile dynamics compared to the electron temperature justifies the omission of T i predictions.While n e and Z eff do have an important impact on the current diffusion evolution, through current drive efficiency, bootstrap current and neoclassical conductivity, discharge-to-discharge variations are limited, so that predictive evaluation could be avoided for our purposes.The spatial coordinate ρ is defined as the square root of the enclosed toroidal magnetic flux, normalized with respect to the square root of the total toroidal magnetic flux enclosed by the last closed flux surface.A relatively coarse time step of 50 ms was found to be sufficient to simulate the salient features of the time evolution of T e and q, allowing to simulate a full AUG discharge (∼10 s) within 2 min on a single CPU.

Experimental reconstruction of equilibrium and radial plasma profiles
To compare simulations with experiments, this paper relies on the Integrated Data Analysis (IDA), a tool based on Bayesian probability theory developed at IPP [19,20].The radial profiles of T e and n e are inferred on a normalized poloidal flux coordinate ρ pol with a temporal resolution of 1 ms, combining information from Thomson scattering, the interferometer, the lithium beam and the ECE.Local error bars reflect at which radii data is available and whether inconsistencies arise between the data from different diagnostics.
The Integrated Data Analysis Equilibrium (IDE) aims to produce improved post-discharge equilibrium reconstructions and current profile estimates (equivalently: safety factor q estimates), by coupling an inverse kinetic Grad-Shafranov equilibrium solver to a predictive current diffusion solver.Current density source profiles are evaluated with TORBEAM [21] for EC driven current and with RABBIT [22] for NBI driven current.

Equilibrium geometry.
For all discharges discussed in this paper, the kinetic equilibrium reconstruction code IDE [23] has been run.The time-varying equilibrium geometry obtained by IDE is interfaced to the format required by RAPTOR, by running the CHEASE fixed-boundary equilibrium solver [24], based on EQDSK files of the IDE equilibrium at various times along the ramp-up, flat-top and ramp-down phases of the discharges.The geometric data required by RAPTOR is extracted from the CHEASE output files, as described in [25].The time step between equilibria in RAPTOR is increased from 100 ms to 200 ms to 1 s from the early ramp-up (where the geometry changes rapidly) to the flat-top (where the geometry stays mostly unchanged).As the equilibrium data processing is performed in an automated way, the user can easily increase or decrease the number of IDE equilibria to be loaded and processed.Typically, the ramp-up simulation starts at t = 0.2 s where I p has a value around 350 kA.Starting the simulation early in a relatively cold plasma (with fast current diffusion) reduces the impact of the initial conditions on the plasma state in the simulation.
In these simulations, the time derivative Φb (toroidal flux enclosed by the last closed flux surface) has been neglected (see equations ( 1) and ( 2) in [9]).For the present study, this assumption has only a small impact and has been verified a posteriori in a few cases.However, for faster ramp-ups or faster plasma shape changes, the impact of Φb can be significant and the corresponding term should be included in the diffusion equations.

Heating and current drive sources.
To evaluate auxiliary current drive, as required for the current diffusion equation evolved within the IDE, the TORBEAM [21] and RABBIT [22] codes have been run.The time dependent heating and current density deposition profiles for respectively electron cyclotron and neutral beam injection are used as inputs to the RAPTOR simulations: • The NBI profiles for electron heating p nb,e and current density j nb are directly provided as sources for the RAPTOR electron heat and poloidal flux equations.• A post-shot processing algorithm is run based on the TORBEAM heat deposition profiles for the individual EC sources.The workflow is illustrated in figure 1: at each time step and for each gyrotron, a Gaussian is fit through the heating deposition profile (black dashed curves on figure 1).From these Gaussians one can extract the time traces for deposition radius ρ dep and width w dep that are shown in figure 1.A smoothed ρ dep (t) trace and a time-averaged w dep are provided as inputs to the RAPTOR Gaussian electron cyclotron module (described in [10]), resulting in Gaussian heating deposition profiles for the EC sources corresponding to the individual gyrotrons.Each Gaussian deposits the full gyrotron power P gyro (t) injected in the plasma, assuming full absorption of the microwaves.The current drive deposition profiles are calculated based on the ad-hoc formula equation (1), introduced in [10] For all EC sources that drive counter-ECCD, the tuning parameter c cd is set to −11. Figure 2 shows how this value of c cd allows for a reasonable match of the total driven I ec and the TORBEAM prediction, indicating that the formula successfully captures the main trends of current drive efficiency degradation when depositing further off axis (the maximum deposition radius of the outermost EC source is given for reference in each panel).Note that the I ec offset between RAPTOR and TORBEAM for discharge 36087 (2-4 s) and 40398 (3-5 s) is due to an overprediction of the electron temperature in RAPTOR with respect to IDA, due to a confinement evolution observed in experiment and not captured by the reduced transport model.The significant counter-ECCD near the magnetic axis can lead to the formation of a current hole, reducing the total current density in the center to small or negative values.A zero or negative integrated value I p (ρ) = ´jpar dA causes divergence of q(ρ), which prevents the RAPTOR implicit solver from converging and ends the simulation.Furthermore, even when negative central current would be expected from conventional current diffusion, negative current densities can be avoided by a resistive kink MHD instability [26], consistent with observations in advanced scenarios in JET [27] and JT-60U [28].
The formation of negative central current in our simulations can be avoided by imposing a lower limit on the deposition radius of the innermost EC sources.For the simulation corresponding to figure 1, a lower limit ρ dep > 0.08 was imposed on the EC deposition radii.This obviously impacts the inner q profile for ρ ≲ 0.15.However, the impact on the simulated q profile further outward is expected to be small (since q depends mainly on the local j par (ρ) and the integrated value I pl (ρ)).

Electron heat diffusivity.
Experimental observation and theoretical understanding allow to construct simple formulas that contain key dependencies of the electron heat diffusivity on plasma parameters.For example, [29] found the ratio of χ e to q 2 to be constant in the core region for a sequence of TCV L-mode plasmas.For the RAPTOR code, a simple adhoc electron heat diffusivity model was introduced in [10] and extended in [30], proposing the analytical formula equation (2) to evaluate the radial distribution of χ e based on local profile quantities Formula equation (2) consists of three terms: (1) a small neoclassical contribution; (2) an anomalous contribution; (3) a Gaussian diffusion term in the center to reproduce the experimental observation of profile flattening.In [30], an additional factor T e0 [keV] cT e was added to the anomalous diffusion term, to capture the degradation of confinement for increased input power.This exponent can be loosely related to the power degradation exponent in confinement scaling laws 3 .
3 Assume a stationary plasma where the electron heat flux equals the total input power P: Qe = neχe ∂Te ∂ρ = P.A power degradation exponent αP for electron thermal energy confinement time τ th e = W th e /P ∼ P αP implies   Table 1 summarizes the values of the model parameters applied for the simulations presented in this paper.Note that no additional transport was added in the core, indicating the experimental T e profiles are relatively peaked in the center (most simulated and reconstructed q profiles remain above unity, so no sawtooth instabilities are expected).Furthermore, we did not activate the shear-dependence in the anomalous transport term (as the shear-dependent factor was set to unity, we did not include it in equation ( 2)).

Te boundary condition. The boundary condition for
the T e equation is imposed at ρ = 0.8.In the region ρ = [0.81], a linear temperature pedestal is imposed.Note however that the current diffusion equation is solved on the full radial domain ρ = [0 1], so that the bootstrap current driven by the pedestal is properly taken into account.The time trace for T e | ρ=0.8 is taken from the IDA inference of T e measurements [19].Within the simulation, the H-mode pedestal builds up as a direct consequence of the increase in T e boundary condition at ρ = 0.8.The LH transition timing is hence not provided explicitly.
Figure 3 illustrates the quality of the performance of the RAPTOR T e (ρ, t) predictions.The simple transport formula provides a robust prediction of core electron temperatures, when T e | ρ=0.8 is provided as an input.

Routine post-shot full-discharge simulations for model validation.
By routinely performing post-shot simulations and comparing the modeled T e (ρ, t) and q(ρ, t) data with the experimentally inferred profile dynamics by IDA and IDE, the applied set of reduced physics models is validated.A fulldischarge simulation is presented in figure 4: the main dynamics of T e and q are successfully captured, from ramp-up to ramp-down.The under-prediction of the outer q profile with respect to the IDE reconstruction during the flat-top phase is W th e ∼ P 1+αP .As W th e = ´3 2 neTedV, one can loosely relate T e0 ∼ P 1+αP (assuming constant density).Introducing χe ∼ T cT e e0 in the power balance yields P ∼ neT cT e e0 Te0 a ∼ (P 1+αP ) cT e +1 , hence (1 + αP)(cT e + 1) = 1 or cT e = 1 1+αP − 1.Note that typical values for the power degradation exponent αP = −0.5, −0.6 and −0.7 hence correspond to the respective Te exponents cT e = 1.0, 1.5 and 2.3, in line with cT e = 1.2 applied within this paper [30].
due to unexpected inconsistencies between the current diffusion evolution and the q profile resulting from the IDE reconstruction, where magnetic measurements and pressure constraints are taken into account.These measurements are not directly considered in the RAPTOR modeling.The nature of these inconsistencies is still under investigation.
In appendix B of [31], a more extensive set of comparisons between RAPTOR post-shot simulations and IDE reconstructions is shown, over a range of discharges from the AUG advanced scenario campaign with counter-ECCD (simulations are shown for most of the discharges mentioned in this paper).
For the T e and q profile estimates in figure 4, the respective local error bars provided by IDA and IDE are represented by a blue shaded area.These error bars should be interpreted with care.The T e error bounds indicate where the metric χ 2 (sum of the squares of the data residuals) increases by one, as explained in [20].Note that for the IDE q profile estimates presented in figure 4, no internal current measurements were available, leading to large error bars.Furthermore, the error bar is calculated without considering the current diffusion constraint (only the magnetics and the pressure constraints determine this error bar).Further details on the IDE error bar can be found in [20].
The successful model validation gives confidence that the proposed set-up can be applied to interpret discharges, guiding further scenario development.In section 2.3, we discuss how predictive, pre-shot simulations can be performed, enabling us to simulate the impact of adjustments to pre-programmed actuator traces and to find optimal actuator traces by solving a dynamic optimization problem.Since the total heating power is one of the input variables to the proposed scaling law, a pre-shot estimate of the ohmic heating power is required.In appendix A, we present an ad-hoc formula that can be used to evaluate the expected ohmic power, based on the pre-programmed plasma current and auxiliary power.A more consistent approach, using the ohmic power evaluated inside RAPTOR could be implemented.However, the present implementation was found to perform well for the discharges studied in this paper.Note that even though the ohmic heating power is usually negligible with respect to the total auxiliary heating power, the ohmic contribution can be important in a relatively cold plasma with little auxiliary heating.Furthermore, early in the ramp-up simulation, before the onset of the auxiliary heating sources, ohmic heating is the only heat source.
To derive a simple scaling law for n e T e | ρ=0.8 , we test the assumption of a power law expression with input variables n el , P aux + P oh and I p , i.e.
Equation ( 4) can be rewritten, taking the natural logarithm on both sides with y the response variable of the scaling law.Linear regression allows to obtain the power law coefficients α 0 , α n , α P and α I , based on the available data set {n e T e | ρ=0.8 ; I p ; P aux + P oh ; n el }.The available data consists of all IDA time points of the discharges performed earlier to develop the scenario.No discrimination is made between data points during ramp-up, flat-top and ramp-down phases, nor between H-and L-mode phases (the early heating discharges presented in this paper typically only have a short L-mode time window).Since the RAPTOR simulations and optimizations reported here have been performed in parallel to experimental developments, only discharges executed before a given simulation was performed could be used, hence the scaling law evolved with time.After a discharge is performed, the power law coefficients are updated, making use of the additional data from the latest discharge.The implementation allows the user to select any set of discharges to derive the scaling law.An example is shown in figure 5, where based on a set of 8 AUG advanced scenario discharges, the following scaling law is derived:  for IDA (blue solid) and RAPTOR (red dashed); (c) q at ρ = [0.40.6 0.8] and q 95 , for IDE (blue solid) and RAPTOR (red dashed), q 95 of CHEASE equilibria (black dots); (d) Te profile at various times, for IDA (blue solid) and RAPTOR (red dashed), the IDA error bar is indicated with a blue shaded area; (d) q profile at various times, for IDE (blue solid) and RAPTOR (red dashed), the IDE error bar is indicated with a blue shaded area.
To further illustrate the performance of the fit, the n e T e | ρ=0.8 predictions for six of the advanced scenario discharges are compared to the IDA time traces in figure 6.Note that even without discriminating H-and L-mode data points, the full time evolution during the ramp-up phase is well captured, encouraging further application of the scaling law approach.

Comparison pre-and post-shot simulation.
Predictive pre-shot simulation and optimization rely on the fact that distinct features can be reliably predicted.For the purposes of this paper, the time evolution of the q profile will be the main focus of our interest.In figure 7, pre-and post-shot simulations are compared, for an advanced scenario discharge (40398) executed on AUG with counter-ECCD.For both the evolution of T e (ρ, t) and q(ρ, t), pre-and post-shot simulations are in close agreement.For the predictive simulation, the T e boundary condition is imposed based on the scaling law approach introduced earlier in this section, deriving scaling law exponents based on data available from previous discharges.The evaluation of n e T e | ρ=0.8 is used to calculate the T e | ρ=0.8 boundary condition, assuming no shot-to-shot variation for n e | ρ=0.8 if fueling remains unchanged, as justified   The density evolution (including n e | ρ=0.8 ) for the pre-shot simulation is imposed from the IDA data available from an earlier discharge (39342), accounting for a later NBI and fueling onset timing by delaying the rising edge of the n e (ρ, t) time evolution.For these discharges, the neutral flux density is feedback controlled by means of deuterium fueling, as a proxy for the electron density at the separatrix.In the present work we make the assumption that the density evolution across the Hmode barrier is mainly determined by this separatrix boundary condition, and that changes of the central density profile due to local changes of the heating mix are small.The flat-top density profile is very similar for both discharges 39342 and 40398 (see n e time traces at ρ = 0 and ρ = 0.8 in figure 7), justifying the applied procedure.
Note however that there is some discrepancy between the q profile predicted by RAPTOR and the q profile of the IDE reconstruction (blue).The impact of a different initial state vanishes around 1 s, but later on the RAPTOR outer q profile is slightly lower compared to the IDE reconstruction.This discrepancy is similar to the q profile mismatch shown in figure 4 and requires further investigation.Note however that it does not modify significantly the dependence of the core q profile on actuators, which is the focus of the following sections.

Relation to other modeling frameworks
Note that the simulations presented in this paper rely on a loose coupling between the transport solver RAPTOR, the equilibrium solver CHEASE and the RABBIT module to evaluate heating and current drive from NBI that is evaluated as part of the IDE equilibrium reconstruction tool-chain.
As mentioned before, the reconstructed IDE equilibrium is interfaced to RAPTOR through the fixed-boundary equilibrium solver CHEASE.Note that when pre-shot simulations are performed, optimizing actuator time traces while maintaining IDE equilibria from a previous discharge, an inconstancy can arise between the kinetic profiles simulated by RAPTOR and the underlying equilibrium.However, iterative application of the RAPTOR transport solver and the CHEASE equilibrium solver allows to obtain a consistent solution of the kinetic profile evolution and the equilibrium geometry evolution.Usually the impact of iteration to a consistent solution on the resulting kinetic profile evolution is modest, if the equilibrium changes are not too large (illustrated e.g. in [32]).
While the electron cyclotron current drive efficiency is based on formula equation (1), the NBI current deposition, as well as the respective NBI heat source profiles to ions and electrons, are evaluated based on a RABBIT run performed by IDE for a previous discharge with similar plasma density.The omission of the direct coupling to separate heating and current drive codes within RAPTOR allows to accelerate scenario optimization.To verify the consistency of the obtained scenario, the simulation can be repeated in a transport solver that features direct coupling to (reduced) heating and current drive models.The ASTRA-based simulator proposed in [33] includes such direct coupling to the TORBEAM [21] and RABBIT [22] codes, ensuring consistency between the evolution of the kinetic profiles and the evaluation of the deposition profiles of respectively EC and NBI.Both the RAPTORbased simulator presented in this paper and the ASTRA-based simulator in [33] rely on experimental data from previous discharges to impose the density profile evolution (which is justified when the line-averaged density is feedback controlled and the density profile shape does not vary significantly from shot to shot).
Note that this type of loose coupling is what is usually used in real-time, where each module (transport, equilibrium, heat source models, . ..) contribute to the global solution of the next time step.A demonstration example is the kinetic equilibrium reconstruction performed in real-time with RAPTOR and LIUQE on TCV [34].Leveraging the large current drive efficiency near the plasma center, near-axis counter-ECCD allows to obtain an elevated q profile in a regime approaching dimensionless parameters close to what is envisioned for DEMO (36087, I p = 1 MA: β N ∼ 2.6, H 98y,2 ∼ 1.15, q 95 ∼ 3.9).Pursuing reactor-relevant operating regimes allows to quantitatively validate the available theoretical models, increasing confidence for their application on future machines, while experimentally determining the operational limits, e.g.β N limits for the onset of (resistive) MHD modes.To match the total imposed plasma current I p , in the presence of a large, negative EC current drive contribution, the ohmic current density is raised over the entire plasma radius.As the EC current density is centered at radii ρ < ∼0.3, an off-axis peak in the total parallel current density results, consistent with the desired elevated q profile [35].While a counter-ECCD scenario is interesting to pursue a reactorrelevant advanced scenario operating regime on present-day machines, the large negative auxiliary current is at odds with the aim of approaching non-inductive conditions.As a consequence, the elevated q profile needs to be sustained in a different way in a fusion reactor.The application of off-axis co-ECCD to maintain q > 1 for the ITER hybrid scenario has been modeled and optimized in RAPTOR in [12].

Model-based scenario development for fast and reliable access to advanced scenario
Let us list some of the properties that make this scenario attractive for a fusion reactor, and formulate some open physics questions that the AUG advanced scenario program aims to answer: • Maintaining q > 1, sawtooth instabilities can be avoided.
The absence of sawtooth seed islands allows to raise β N into the regime where NTMs are metastable [36,37].In [5,38], an upper limit has been found β N ∼ 2.7, above which ideal 2/1 modes are triggered.Note that low magnetic shear at rational surfaces can destabilize ideal MHD infernal modes, as investigated for TCV in [39].Furthermore, low magnetic shear at rational surfaces can increase the classical tearing mode stability parameter ∆ ′ , thus making the current density profile more prone to the onset of NTMs, as described in [40] for DIII-D q 0 > 1 scenarios.Similarly, [4] observed the onset of NTMs when q min drops through the corresponding rational surface, for AUG early heating scenarios.Even in the absence of sawteeth, so-called triggerless NTMs can originate [41], while ELM crashes can provide an alternative seed mechanism [42], especially for 2/1 NTMs, as the q = 2 is in close proximity to the plasma edge during operation at high plasma current.
The proximity to the edge furthermore increases the likelihood for mode locking and subsequent disruption [43].
• An elevated q profile allows to increase β pol , maximizing the bootstrap fraction ( f bs = I bs Ip ∼ β pol ∼ β N q/ϵ, with inverse aspect ratio ε).
• A high pedestal pressure can be obtained (a comparison for the electron pedestal pressure to the AUG ITER baseline scenario is shown in figure 3.5 of [31]), as advanced scenarios operate at high β pol , increasing the critical threshold of peeling-ballooning modes through an increased Shafranov shift [44].
• In [45], the impact of the divertor neutral density n 0,div on global confinement quality has been identified, highlighting the importance of edge and scrape-off layer physics.Data analysis and pedestal stability modeling indicate that an increased n 0,div causes an outward shift of the location of the maximum pressure gradient in the pedestal, where the maximum pressure coincides with a larger value of the safety factor q.This eventually leads to a larger destabilizing drive for ballooning modes and a reduced pedestal top value (while increased n 0,div is correlated with a measured increase of n e,ped , a more significant reduction of T e,ped and T i,ped is incurred [45]).To obtain good and reproducible performance in AUG AT discharges, the neutral flux density (measured by ionization gauge [46]) is feedback controlled, using deuterium fueling as actuator [47].As discussed in section 2.3.2,we assume that the separatrix density, controlled by proxy of the neutral flux, is the main physics driver determining the core density profile.While some minor density profile changes have been observed after ECCD modifications, see e.g.section 4.3.2, the assumption that the density profile evolution is similar from shot to shot, turned out to be robust.• A central physics aim of the AUG advanced scenario program is the identification of drivers for enhanced ion transport, studying the respective impact of q profile, E × Bshear, ion heating and β pol [35,48].It is well known that counter-ECCD at low density can lead to the formation of an ITB in the electron heat channel (e-ITB) [49].By operating at a high plasma current of I p = 1 MA, the present scenario reaches a higher density and a more significant coupling between ion and electron temperatures compared to advanced scenarios at lower I p [35], avoiding the formation of e-ITBs and the central accumulation of heavy ions.Note that the simple transport formula equation ( 2), is successfully applied over a wide range of discharges, as shown in figure 3, while no shear dependence of confinement is needed to match the experimental T e profiles (a sheardependence would be expected in the presence of e-ITBs).
The need for further optimization of this rather challenging scenario has been discussed in [35].The 1 MA discharge discussed in that paper (shot 36087) applies external heating only after the ohmic current density reaches a relaxed state at full current.Now we will discuss a first attempt to reach this scenario with an early heating strategy.for 36087 and 39342.Note that for the late heating discharge 36087 the off-axis minimum of q is formed while locally q < 1.5, while for the late heating discharge 39342 the off-axis minimum of q is formed while locally q > 1.5, leading to a time point after 1.8 s when q = 3/2 and the local magnetic shear is zero.

Understanding NTM-triggering in an early heating scenario attempt
While the plasma dynamics are highly sensitive to actuator time traces during transient phases, profile simulators allow to evaluate the deterministic impact of actuator modifications.Figure 8 presents RAPTOR post-shot simulations of the AUG discharges 36087 (late heating) and 39342 (early heating attempt), providing insight in the trade-off between late and early heating.Figure 10 shows the experimental time traces for β pol and H 98y,2 .Note that during the flat-top phase, these discharges control the value of β pol to a pre-programmed time trace, by modulation of neutral beam power injection.
• For the late heating discharge 36087, the q profile descends towards q min ∼ 1, before slowly rising towards the desired elevated state.According to the RAPTOR simulation q min ∼ 1 around 2 s, and the elevated q min is achieved around 3.5 s (see panel (c) in figure 8).• By applying the neutral beam heating earlier, during the plasma current ramp-up, the q profile of the early heating discharge decreases monotonically in time towards the final state, reaching the desired, elevated q profile earlier in time.
The early heating approach is however more prone to the onset of tearing modes.As visible in the magnetic spectrum for shot 39342 in figure 9, a 3/2 mode is triggered around 1.7 s.After the onset time of the tearing mode, confinement degrades: the central electron temperature, measured by the ECE diagnostic, drops from T e ∼ 8.0 keV to T e ∼ 5.5 keV and the experimental H 98y,2 factor reduces from ∼1.1 to ∼1.0 (see figure 10).Since the onset of tearing modes during the early flat-top phase had been anticipated when setting the reference for β feedback control, the plasma is maintained at β pol ∼ 1.2 for several seconds, allowing the tearing mode to disappear before raising β pol (and β N ) further around 3.5 s.
The β pol ramp leads to a confinement increase in the final phase of the discharge, with H 98y,2 ∼ 1.15 (see panel (b) in figure 10).
As illustrated in figure 8, the RAPTOR simulation hints that a 3/2 NTM is triggered due to the formation of an offaxis q min ∼ 1.5 at large radius (ρ ∼ 0.4), roughly 1 s after the flat-top plasma current has been reached.As the auxiliary heating and current drive start shaping the q profile before the profile has decreased and equilibrated, an off-axis minimum is created with q min > 1.5.The gradual, monotonic decrease of the q profile continues during the flat-top phase (while the plasma has already a relatively high β N ∼ 2.1) and brings q min towards the rational value q = 3/2.Low magnetic shear around the rational surface could explain the formation of a 3/2 NTM, similar to the observations in [4,40].

Predict NBI onset timing to avoid deleterious 3/2 mode
In this section, a first practical example of inter-discharge optimization is presented.Based on predictive RAPTOR simulations, quantitative changes to the discharge program have been tested in experiment.While in the present section, this procedure relies on running different simulations with manual adjustments to an actuator trace, an automated optimization scheme will be introduced in section 4.

Physics interpretation based on post-shot simulation.
The starting point is the set of interpretative simulations of previous discharges, namely the post-shot simulation of the late heating discharge 36087 and the early heating discharge 39342.As argued in the previous section, these simulations allow us to formulate a hypothesis for the cause of the onset of the tearing mode in the early heating attempt: having a minimum in the q profile drop through the rational surface q = 3/2 at relatively large radius (ρ > 0.2), with β N > 2, is undesirable, as it might cause the onset of a confinementdegrading NTM.

Definition physics aim for next discharge.
Based on the interpretation of the previous discharge, a qualitative change to the discharge program for the next shot has been proposed, namely a delay of the NBI, with the following rationale: delaying the NBI injection extends the period during which the plasma is maintained at relatively low temperature, in L-mode, allowing the q profile to decrease more rapidly in time.
The aim of the predictive step, relying on pre-shot RAPTOR simulations, is to find the optimum time delay of the NBI heating, allowing for the formation of the off-axis q min at a time when locally the q profile is already below q = 1.5.By avoiding an off-axis minimum of q with q min = 3/2, we aim to make the scenario more robust with respect to the onset of a tearing mode.
As the q profile dynamics are extremely sensitive to the timing of heating during the ramp-up phase, finding the optimal NBI time delay would be difficult in a model-free approach and might have required several discharge iterations, attempting different shots with varying NBI time delay.In the present work, pre-shot simulations have been used to obtain a quantitative estimate of the change in the q profile dynamics, as explained in the following section.The basic idea is to reduce the number of discharges required for scenario development by testing actuator updates first in simulations.

Predictive, pre-shot simulation.
The interpretative step is followed by a predictive step, resulting in quantitative changes to the pre-programmed actuator time traces for the following discharge.
Predictive RAPTOR simulations have been executed, delaying the onset timing of the NBI sources that are activated during the ramp-up.After attempting simulations with various time delays, ∆t = 200 ms has been selected 5 .The predictive RAPTOR simulation implementing a ∆t = 200 ms delay of the NBI heating during the ramp-up phase is shown in red dash-dotted lines in figure 11.The time traces for T i (ρ, t), n e (ρ, t) and boundary condition T e (ρ = 0.8, t), are constructed based on measurements from the previous attempt to establish an early heating scenario (39342).Starting from discharge 39342, the traces are kept constant between 0.5 s and 0.7 s, reflecting the delayed impact of the NBI onset.After 1.27 s, when β pol control is initiated, the original traces are used again, without shift in time.
At 1.2 s, the RAPTOR simulation indicates that the q profile has decreased further compared to the previous early heating attempt, as anticipated due to the longer L-mode phase.The q profile drops through q = 1.5 before β pol is raised and before auxiliary heating and current drive leads to the formation of an off-axis q min .Furthermore, q min = 1.5 occurs at smaller radius.If an NTM is triggered at this earlier time and lower β pol value, the impact on confinement is expected to be smaller [53].At later times, non-zero magnetic shear is maintained at q = 1.5.
The eventual output of the predictive step is a proposed actuator trace for the following discharge.In this case, the NBI trace proposed by the RAPTOR simulation has been included in the discharge program for the subsequent shots.In the following section, we will assess whether the proposed actuator trace modification has been successful.

Experimental validation of the optimized NBI actuator
trace.The present section gathers relevant experiential data to illustrate that the proposed modification to the NBI onset time has successfully increased the robustness of the early heating scenario for the onset of confinement-degrading tearing modes.Note that for some of these scenarios also the pre-programmed radial deposition radii of ECCD have been optimized in RAPTOR, as will be discussed in section 4.
All the discharges presented in this section include the 200 ms delay that has been proposed by RAPTOR in the discharge program.In the magnetic spectra of discharges 40029 and 40030 (the spectrogram for discharge 40030 is shown in figure 12, a more extensive overview of spectrogram data can be found in appendix B of [31]), no n = 2 signature is visible during the early flat-top phase.For later shots, a n = 2 mode reappears in the spectrogram, usually initiating around 0.8-1.1 s (40187: n = 2 around 0.8 s; 40188: n = 2 around 1.0 s; 40192: n = 2 around 1.1 s in figure 13; 40398: n = 2 around 1.1 s; 40825: n = 2 around 1.1 s).However, as the mode is formed at a smaller radius and at an earlier time, the impact of the mode is reduced.No clear impact of these early 3/2 modes on the central electron temperature is visible in these experiments.While a non-zero shear is maintained at q = 1.5 during the early flat-top phase, the n = 2 activity disappears, usually around 1.8 s.Excellent confinement can be achieved, even while the early 3/2 mode is present.In discharge 40188, H 98y,2 ∼ 1.15 as β pol control overshoots with β pol ∼ 1.35 around 1.6 s (reducing to β pol ∼ 1.25 afterwards, with confinement degrading as the steep ∇T i around ρ ∼ 0.55 relaxes; confinement recovers during the late β pol ramp after 3 s6 ).For further details on the achieved parameters in the discharges discussed here we refer the reader to table C1, figure C2 (time traces H 98y,2 , β pol and β N ) and figure C1 (time traces T i at various radii) in appendix C.
In the absence of the deleterious NTM, the β pol reference is raised during the early flat-top phase, to β pol ∼ 1.35 − 1.40 for discharges 40192 and 40398.Discharge 40192 features excellent confinement (H 98y,2 ∼ 1.10) and high ion temperatures (T i0 ∼ 8 keV), until a 3/2 NTMs abruptly reduces core electron and ion temperatures around 3.1 s (likely triggered by a growing n = 3 mode, see spectrogram figure 13).The mode disappears at 4 s and reappears at 4.8 s.The continued sporadic occurrence of (short) time windows with 3/2 modes during the later flat-top phase indicates that the scenario continues to be metastable regarding NTMs.Discharge 40398 reaches similar core T i and T e values, with a 3/2 mode appearing around 3.8 s.
Finally, we can compare the pre-shot RAPTOR simulation that has been used to propose the 200 ms time delay for the NBI onset with a post-shot simulation of the discharge 40192, for which the proposed delay has been executed.The post-shot simulation of discharge 40192 (green dotted lines in figure 11) recovers quite closely the pre-shot predicted q profile evolution (in red dash-dotted lines).While this is an indication for the validity of our approach, it should be stressed that the present feedforward approach is incapable of handling shotto-shot variations that affect the relevant dynamics and that are not captured by the applied reduced physics models.

Model-based optimization of stationary, elevated q min scenarios
After improving the robustness of early flat-top phase against the onset of tearing modes, the RAPTOR code is used to optimize the ECCD deposition profile, to propose the number of EC sources and the deposition radii to be applied in the experiment.A distinctive feature of the RAPTOR code is the availability of automated optimization routines, as introduced in [10].A fast run time and analytical gradients make a full time dependent non-linear optimization problem between shots computationally tractable.For further development of the 1 MA counter-ECCD scenario, the optimizer is used to assess the following questions: (i) For a given amount of total EC power, what is the maximum q min that can be maintained stable?(ii) What is the optimal EC deposition profile, allowing to reach the final elevated q profile as early as possible in the discharge?

Mathematical formulation of the optimization problem
For the optimization problems presented below, the cost function is a measure for the stationarity of the plasma: Minimizing the radially integrated squared loop voltage gradient, integrated over a time window (t ∈ [t 1 t 2 ]), we minimize the driving force for current diffusion, as explained in [10,54], maximizing the stationarity of the plasma state.
Two state constraints have been used, imposing minimum values of the safety factor q min and the magnetic shear s min for the final state in the simulation t f .The mathematical formulation of these constraints as integrals has been introduced in [10] ) With these cost and constraint functions, the optimal control problem can be defined, as explained in more detail in [10] C q>q min ⩽ 0 and C s>s min ⩽ 0 (state constraints).(10d) The RAPTOR state evolution equation f defines the time evolution of the plasma state x(t) (effectively the radial distribution of T e and ψ), over the simulation time interval t ∈ [t 0 t f ].The optimization variables p contain a subset of the deposition radii of the EC sources, as explained in the next section.This optimal control problem is solved by application of the sequential quadratic programming (SQP) algorithm, implemented in the Matlab function fmincon, as discussed in [10].
4.2.Optimization for early stationary q profile with q min < 1.5

Set-up of the predictive simulation used within the optimization routine.
The automated optimization scheme is used to find the optimum EC deposition profile and the adequate number of EC sources to achieve the desired q profile at an early time.The optimized NBI onset time t = 0.7 s, found in section 3.3, is included to avoid the onset of tearing modes during the early flat-top phase.
For the RAPTOR simulations performed within the optimization routine, the following set-up is used, starting from the interpretative post-shot simulation of discharge 39342.
• Density n e (ρ, t), ion temperature T i (ρ, t) and neutral beam deposition p nbe (ρ, t), p nbi (ρ, t), j nb (ρ, t) are set identical to the corresponding profile evolution in the post-shot simulation of the earlier discharge 39342.The profiles are maintained time-invariant between t = 0.5 s and t = 0.7 s, reflecting the delay in NBI power, as discussed in section 3.3.• The pedestal temperature T e | ρ=0.8 is calculated with a scaling law, as discussed in section 2.3.Obviously, the scaling law is based on the data from the set of discharges performed before the discharge that is being optimized.• A fixed number of EC heating and current drive sources are used, each providing counter-ECCD, delivering a total power P ec .The optimizer adjusts the EC deposition profile for a given total EC power.Various optimizations have been performed, with different numbers of EC sources activated (corresponding to different amounts of total power P ec ).
The EC sources start delivering power simultaneously in the early flat-top, around 1.2 s, right before β pol control by NBI modulation is activated.If the optimizer is not able to find a feasible solution for a given number of activated EC sources, a novel optimization attempt is launched where an additional EC source is activated in the RAPTOR simulation.According to these RAPTOR simulations, four active EC sources are needed to satisfy the constraints q > q min and s > s min (equation (10d)), depositing a total power P ec = 2.8 MW.
The radii of two EC sources are kept fixed at ρ dep = 0.08, to force some EC heating close to the magnetic axis, avoiding tungsten accumulation [55] 7 .The radii of the two other EC sources are used as optimization variables p = [ρ dep 1 ρ dep 2 ] T .By varying the off-axis distribution of EC current density, the optimizer can tailor the q profile.Both ρ dep 1 and ρ dep 2 are constant in time.

Cost and constraint function settings for the optimization routine.
To aim for a q profile at 3.5 s that is elevated, but has non-zero shear through the q = 3/2 surface, constraints impose q min = 1.35 for equation (8) and s min = 0 for equation (9) (weight functions W q min and W s min are used to activate the constraint over ρ = [0.1 1] and ρ = [0.25 1] respectively).The second constraint enforces the desired q profile to increase monotonically from the minimum outward, avoiding an oscillatory behavior where a local extremum can potentially cause zero magnetic shear at the rational surface q = 3/2.The cost function equation (7) integrates a measure of stationarity over the time window 1.5 s to 3.5 s, with the aim of reaching the final q profile as early as possible.Before 1.5 s, the delayed NBI recipe developed in section 3.3 allows for a reliable access into the high β pol phase.The present optimization makes sure that the final state is consistent with the stability constraints to avoid tearing modes (non-zero shear through q = 1.5 8 ) and the goal of an elevated q profile (q min > 1.35), while reaching this state as early as possible once the current ramp is completed and the reference high β pol value is reached (around 1.5 s).

4.2.3.
Pre-shot discharge optimization.Figure 14 presents the initial condition provided to the optimizer in blue (p = [0.40.4] T ), while the obtained optimum is shown in red dashed lines (p = [0.220.30] T ).The optimum q profile at 3.5 s indeed satisfies the constraints, while the q min value is maintained relatively stationary, within the interval q ∈ [1.35 1.45] after 1.5 s, indicating that a q profile close to the final state is achieved immediately from the onset of the EC power.The resulting pre-shot simulation hence combines the stable pathway avoiding or minimizing the impact of tearing modes with an elevated q min , right from the start of the high β pol phase and will be tested experimentally in section 4.2.4.
Due to the fast run time of the RAPTOR code, various sensitivity studies can be run at low computational expense.An example is shown in figure 15: the optimized simulation is repeated with a T e | ρ=0.8 20% higher and lower compared to the trace predicted by the scaling law.The corresponding q profile at 3.5 s and q min time trace is significantly affected: an increase in temperature pedestal leads to a larger amount of bootstrap current driven in the pedestal, leading to a more elevated q profile with q min above 1.5.Reversely, the degraded pedestal would lead to a q min around 1.2.The complex interaction between the current in the pedestal and the core will be discussed in more detail in section 5.3.At this point, we can conclude that feedforward optimization of the q profile relies on good reproducibility of the pedestal top conditions.

Experimental validation of the optimized discharge program.
The pre-programmed actuator settings for discharge 40398 combine the optimum NBI onset timing proposed in section 3.3 and the optimized EC deposition radii found in section 4.2.3.Figure 15 compares an interpretative, post-shot simulation (green solid lines) to the pre-shot predicted optimum state evolution (red dashed lines).Even though the outer EC sources have deposited further inward than anticipated, the time evolution of the post-shot reconstruction recovers well the general features desired from the optimum: the post-shot simulation shows a q profile with q min between q = 1.35 and q = 1.5 and with little time dependence after 1.5 s.The scaling law evaluation of T e | ρ=0.8 , evaluated before the discharge and used in the pre-shot RAPTOR simulation, matches well the experimental IDA reconstruction.
Let us now compare the performance of the discharge that executed a discharge program optimized with RAPTOR with earlier AT scenario attempts.In figure 16, post-shot RAPTOR simulations are compared for:    Discharge 40398 successfully combines the features of the access scenario described in section 3.3 with the early achievement of an elevated stationary q profile: q min decreases below q = 1.5 early during the flat-top phase, aiming for improved stability margin with respect to a potential 3/2 tearing mode and remains relatively constant after t = 1.5 s (also illustrated by the magnetic shear at the 3/2 surface in panel (d) of figure 16).As has been mentioned in section 3.3, a 3/2 mode is triggered around 1.1 s in discharge 40398.However, the mode disappears around 1.8 s and has no clear impact on confinement (see figure 17), consistent with our guiding hypothesis that the impact of a mode triggered at lower radius and lower β pol is more benign.
The access to an early stationary state in discharge 40398 is confirmed when comparing the RAPTOR post-shot simulation with a post-shot simulation in ASTRA [52] and with the kinetic equilibrium reconstruction of IDE, as shown in figure 18: after 1.5 s, the q profile values at various radii show a relatively constant time evolution, indicating the successful achievement of an early elevated q min .For comparison, the slow time evolution of the q profile values for the late heating discharge 36087, in RAPTOR and IDE, is also shown in figure 18.By finding the optimal allocation (timing, power and radial distribution) of the available actuators, we were able to achieve the high performance phase of the advanced scenario (with elevated q min and elevated β pol ) early during the flat-top phase, while avoiding the onset of tearing modes, as confirmed by the H-factor and β traces in figure 17. Figure C1 in appendix C shows that a high central ion temperature T i0 ∼ 8 keV is reached early, around 1.5 s, and maintained constant.

Towards q min >1.5
To avoid or minimize the impact of tearing modes, the previous section explored a ramp-up strategy for which at high β pol and for ρ > 0.25, a non-zero magnetic shear is maintained through q = 3/2.While the proposed strategy was successfully tested in experiment, allowing for fast and reliable access to a stationary elevated q profile, the stationary flat-top operating point is metastable with respect the onset of 3/2 NTMs (e.g.40192, with short-lived 3/2 NTMs leading to abrupt reductions of the core T e and T i ).Evidently, the onset of 3/2 modes could be avoided by increasing the q profile entirely above q = 1.5.In the present section, we solve an optimization problem to assess the feasibility of this approach with the available heating and current drive resources.Next, the optimized actuator settings are tested in experiment.

Pre-shot discharge optimization.
We use a forward simulation set-up identical to the previous section.One more EC heating and current drive source is added, which, according to our simulations, allows to maintain a stationary safety factor with q min > 1.5.This leads to an increase of the EC power to P ec = 3.6 MW (which impacts the pedestal temperature predicted by the scaling law).Within the present optimization problem, only one EC source is kept fixed at ρ dep = 0.08.The radii of the four other EC sources are used as optimization variables p = [ρ dep 1 ρ dep 2 ρ dep 3 ρ dep 4 ] T .The radii ρ dep i are constant in time.
The constraint equation ( 8) is used to impose q min = 1.55 (with weight functions W q min activated over ρ = [0.1;1]).We choose q min = 1.55 to maintain some margin between q min and q = 1.5 (avoiding low-shear conditions close to the 3/2 rational surface).The cost function equation ( 7) is identical to the one applied in the previous section, integrating a measure of stationarity over the time window 1.5 s-3.5 s.The optimizer effectively looks for the deposition radii that result in a final q profile fulfilling the constraint on q min , while maintaining the q profile as stationary as possible from 1.5 s onwards.
Figure 19 presents the initial condition provided to the optimizer in blue (p = [0.30.3 0.4 0.4] T ), while the obtained optimum is shown in red (p = [0.160.23 0.29 0.38] T ).By optimizing the deposition of the counter-ECCD, the optimizer manages to maintain the q profile above q = 1.55.Note that a broad ECCD deposition profile is proposed, resulting in a large region with q ∼ 1.55 and low magnetic shear, extending up to ρ ∼ 0.4.Between 1 and 2 s the q min trace drops below q = 1.5, which means that a 3/2 mode could still be triggered  Results of forcing fast access to stationary q min > 1.55, by automatically optimizing EC deposition radii.For the initial (blue solid) and optimized (red dashed) state evolution in RAPTOR: (a) q profile at t = 3.5 s (the grey shaded area highlights q < 1.55); (b) time evolution of the minimum of q for ρ > 0.25 (the grey shaded area highlights t > 1.5 s, where the cost function is active); (c) EC power deposition profile; (d) EC power and Te| ρ=0.8 boundary condition (from scaling law).
during this time interval.However, such a mode would disappear when the q profile is increased and the q = 3/2 surface ceases to exist.
In figure 20, the optimized simulation is shown in the red dashed lines, while the black dash-dotted and dotted lines provide a sensitivity study, presenting the simulation resulting for respectively a 20% decrease or increase of T e | ρ=0.8 .The pedestal temperature clearly has a significant impact on the q profile, with q min around 1.3 for the degraded pedestal and q min above 1.7 for the enhanced pedestal.

Experimental validation of the optimized discharge program.
The pre-programmed actuator settings for discharge 41102 combine the optimum NBI onset timing proposed in section 3.3 and the optimized EC deposition radii found in section 4.3.1.Figure 20 compares an interpretative, post-shot simulation to the pre-shot predicted optimum state evolution.The post-shot simulation includes the temperature pedestal and density evolution as measured in experiment (after post-discharge analysis in the IDA), as well as the neutral beam deposition profiles predicted by RABBIT based on the plasma profiles observed in experiment.Even though the ECCD profile provides a broad profile as proposed by the optimizer, the post-shot simulation indicates a reduction of q min significantly below q = 1.5, reaching a stationary value around q min ∼ 1.0.Utilizing the flexibility of the RAPTOR code, further simulations have been performed to understand this discrepancy between pre-shot expectation and post-shot interpretation.
To disentangle the impact of various physics drivers on the resulting q profile, a set of additional simulations has been performed, adjusting incrementally the inputs to the post-shot simulation.Relevant inputs of the post-shot simulation have been swapped with the corresponding input used in the optimized pre-shot simulation, based on expectations before the discharge was launched.
• The measured T e | ρ=0.8 is below the temperature pedestal height predicted by the scaling law, leading to a lower current in the pedestal region (by the end of the simulation window, T e | ρ=0.8 is about 15% below the scaling law prediction) 9 .The green dash-dotted profile in figure 21 indicates that increasing the pedestal boundary condition to the higher pedestal temperature expected before the discharge leads to a significant increase of the q profile, both for the outer profile and in the center.• Post-shot analysis has identified two additional mechanisms for increased current drive in the plasma center: (i) A more significant peaking of the n e (ρ, t) profile in the center leads to a local peak in the bootstrap current, driving additional current for ρ < 0.25 (see panels (c) and (e) in figure 21).(ii) The neutral beam current density profile is more peaked compared to the pre-shot assumption (the central electron and ion temperatures are respectively higher and lower compared to the previous discharges 40192 and 40398, leading to more NBI driven current in the core according to the RABBIT post-shot evaluation).The black dotted profile in figure 21 represents a post-shot simulation that ingests the density profile evolution, the neutral beam deposition profile evolution and the temperature pedestal evolution assumed before the experiment.We observe that central peaking of density and NBI deposition profiles has a significant impact on the q profile in the plasma center, where the total current density is small due to the negative ECCD applied in these experiments.The simulation indicates that if the pre-shot assumptions for pedestal temperature, density and NBI deposition had been achieved, the value of q min would have remained above q = 1.5.
This study highlights the sensitivity of post-shot simulations to the IDA reconstruction (T e | ρ=0.8 and n e ) and the RABBIT neutral beam profile calculations, including local features of the profiles.A recent effort aims to extend the IDA framework with physics-based priors, based on ASTRA simulations, which would avoid spurious oscillations in the reconstructed profile gradients in poorly diagnosed regions [56].

Experiments at reduced plasma current.
Several shots have attempted to apply the increased amount of EC power, to achieve q min > 1.5.Various of these attempts ended prematurely due to the occurrence of a disruptive 2/1 NTM.A salient example is discharge 40825.In the spectrogram of this discharge (not shown here), a sequence of 2/1 modes, likely triggered by ELMs, is visible, before the mode finally locks and leads to a disruption around 2.1 s.While this behavior can be expected due to the metastable nature of the q profile in this scenario, it is unclear whether the evolution of the current density profile decreases the stability.Alternatively, the scenario adjustments could be responsible for changing the ELM size and thus the NTM seed.The RAPTOR post-shot simulation indicates that the outer q profile of discharge 40825 is further radially outward compared to earlier discharges like 40398.Closer radial proximity to the plasma edge could lead to more significant coupling to the ELMs, even if the q profile is not more unstable.Note that the verification of this hypothesis requires further analysis, similar to the work on sawtooth-triggered NTMs in [57].
Based on this hypothesis, a reduction of the plasma current to I p = 900 kA is proposed.Predictive RAPTOR modeling indicates that the increase of q 95 should lead to a decrease of the q = 2 radius, reducing the potential impact of ELM seeds, while increasing the local magnetic shear tends to decrease the classical tearing mode stability parameter ∆ ′ [40], increasing the required seed island size.Maintaining the heating and current drive actuator time traces unchanged, while halting the plasma current rise when I p = 900 kA is reached, a set of experiments have been performed where 2/1 modes could be successfully avoided.For discharge 41400, a stable flat-top phase without any confinement-degrading modes, with β pol ∼ 1.45, β N ∼ 2.3, could be maintained, with an H 98y,2 factor around unity.Both IDE and RAPTOR indicate that q min > 1.5 was not maintained, while the pre-shot RAPTOR run indicated q min > 1.5 can be maintained.Similar to the illustration for discharge 41102 in figure 21, the low q min in the post-shot simulation compared to the pre-shot run is due to a reduced pedestal temperature and a significant peaking of the central density.While unpredicted shot-to-shot variations like these highlight some limitations of model-based shot-to-shot optimizations, we will show in the following section how sensitivity studies on various input variables can be easily run, increasing insight in the robustness and the characteristics of the operating point of a given scenario.

Sensitivity study on the q min > 1.5 scenario
In the present section we illustrate how a set of RAPTOR simulations allows to gain insight in the sensitivity of the optimum.The nominal simulation is the optimum obtained in section 4.3, aiming for q > 1.55 (figure 19).

Sensitivity to initial q profile
The standard initial conditions used for all simulations in this paper are taken from the IDA T e profile and IDE q profile for discharge 39342 at 0.2 s.We choose to keep the same initial condition at t = 0.2 s for all shots, to allow comparison of the different shots without differences induced by a different initial state (which is poorly diagnosed).For the sensitivity study in figure 22, an alternative initial q profile is constructed (more elevated, with q min ∼ 4, indicating a broader current density distribution in the core).Since we start simulating early in the ramp-up, when the plasma is still relatively cold, the difference in q profile has almost entirely disappeared when the NBI is started around 0.7 s and the plasma enters into H-mode.We conclude that the initial state has only a minor impact on the q profile evolution after 0.7 s.

Sensitivity to EC deposition radii
The sensitivity of the optimized q profile to variations in the EC deposition profile is analyzed.For each of a series of 50 simulations, the optimization vector p = [0.160.23 0.29 0.38] T is perturbed, adding a random number to each of the optimization variables, encoding the EC deposition For the RAPTOR optimum state evolution forcing q min > 1.55 (magenta dashed line) obtained in section 4.3, the sensitivity of the q profile with respect to variations in the EC deposition radii is studied.The profiles corresponding to 50 perturbed simulations (adding a random number to each of the optimization variables) are plotted in light grey, together with a mean, an upper and a lower standard deviation for the resulting q profile.(a) q profile; (b) EC driven current density; (c) integrated enclosed EC driven current.
radii (pseudorandom values drawn from the standard normal distribution, with standard deviation 0.02).Figure 23 shows the optimum q profile at 3.5 s in magenta, as well as the corresponding EC deposition profile.The profiles corresponding to the 50 perturbed simulations are plotted in light grey, together with the mean, the upper and the lower standard deviation.We observe relatively significant variations of the q profile at 3.5 s: minor variations in the EC deposition radii can cause significant differences in the q profile: the profiles of the upper and lower standard deviation enclose a band between q = 1.4 and 1.7 around ρ ∼ 0.2.The q profile is quite sensitive to changes in the EC deposition profile since the total current density is quite small with respect to the large (negative) values of the EC current density.

Sensitivity to pedestal temperature: coupling pedestal-central q profile
The impact of the pedestal temperature height on the core q profile of the optimized scenario at 3.5 s (introduced in figure 19) is analyzed in some more detail in this section.The 1 MA counter-ECCD scenario features a combination of a high pedestal with low net current in the core.Figure 24 shows time traces of the integrated plasma current density in the central region (ρ < 0.4, i.e. the region where ECCD is deposited) and in the pedestal region (ρ > 0.8).It is interesting to note that the current (ohmic and bootstrap) driven in the pedestal (∼190 kA) is about two thirds of the net current in the central region (∼290 kA).This is a particular feature of the counter-ECCD scenario: a broad current density profile is created by adding a large negative ECCD in the core to the large inductively driven current (for ρ < 0.4: NBI, bootstrap and ohmic current drive about ∼530 kA, of which about 350 kA is ohmic current, while the ECCD amounts to about ∼−240 kA).The individual current density contributions (and enclosed plasma currents) of the various inductive and non-inductive current sources are illustrated in figure 25.
As shown in figure 24, changing the pedestal temperature boundary condition T e | ρ=0.8 by adding or subtracting 20%, impacts the electron pressure gradient in the pedestal, hence modifying the amount of bootstrap current.The increase/reduction of T e | ρ=0.8 and I bs ped is compensated with a reduction/increase of the plasma loop voltage and I oh in the core.The applied transport model predicts a similar T e0 (transport is enhanced for increasing T e and increasing q, yielding a similar T e0 despite a higher T e | ρ=0.8 ).We can see how the increase/reduction of pedestal bootstrap current is counteracted by a respective reduction/increase of plasma current in the center, through a reduction/increase of the ohmic current (similar T e0 : similar σ neo , similar EC current drive efficiency, similar j bs ).In the outer core region (ρ ∈ [0.4 0.8]) the total current density seems mostly unchanged (figure 24 bottom right), as the reduction/increase of the loop voltage is compensated by increase/reduction of neoclassical conductivity σ neo due to the temperature increase/decrease.
To summarize: the excess/lack of bootstrap current driven in the pedestal ∆I p ped due to an increase/decrease of T e | ρ=0.8 leads to a respective decrease/increase of the plasma current inside ρ ∼ 0.4.For a given geometry and on-axis toroidal magnetic field, the q profile is directly proportional to the inverse of the enclosed net plasma current, more specifically, q ∼ 1 Ip(ρ) with I p (ρ) = ´jpar dA.When, like for the present scenario, the net I p (ρ) in the center is relatively small, the redistribution of ∆I p ped due to a changed T e | ρ=0.8 is expected to have a large impact on the central q profile.This is indeed what we observe in figure 24: as T e | ρ=0.8 is increased/decreased by 20%, q min is significantly affected, achieving respective values q min ∼ 1.75 and q min ∼ 1.30.

Outlook: from inter-discharge optimization to real-time control
In the present paper we have presented how the fast RAPTOR simulator, combining an ad-hoc transport model and a pedestal pressure scaling law, both tuned based on experimental data from previous discharges, has been applied to optimize preprogrammed actuator settings like the NBI onset timing and the ECCD deposition radii.
In the presence of shot-to-shot variations that are not included in the applied model, real-time control provides an interesting avenue to complement the feedforward optimization techniques presented here, through adequate control actions that compensate model-reality mismatches.For the development of real-time control schemes, fast models of the plasma response can be used, either offline, for the development of feedback controllers [58], or online, for model predictive control (MPC) [59,60].
In [58], an interpretative METIS [61] simulation of a DIII-D discharge has been used to derive a linear plasma response model that can be used for the design of feedback controllers for q(ρ) and β N .Furthermore, [58] has illustrated in simulations the robustness of such feedback controllers in the presence of disturbances to plasma parameters like confinement quality, plasma density and the effective charge Z eff .
In MPC schemes, the feedforward actuator signals are continuously updated in real-time, by repeatedly solving a quadratic programming optimization problem, making use of a linear plasma response model and including the relevant actuator constraints.In RAPTOR, a set of linearized models at various times along the reference trajectory can be calculated at a low computational expense, due to the availability of analytical Jacobians [10].MPC based on linearized RAPTOR models has been tested in simulations [62] and successfully applied in profile control experiments on TCV [60].Similar experiments on DIII-D have reported improved q profile tracking by combining feedforward optimization and on-line linear MPC [59,63].
As feedback controllers are aware of the measured plasma state, acting to reduce the offset between the observed state and the reference, the fidelity requirements for models used for real-time control are less stringent compared to pre-shot optimization.Application of feedback control strategies that rely on linear plasma models might be more challenging during the ramp-up phase, where the dynamics are more transient and non-linear.However, even if feedforward control is maintained during ramp-up, feedback control could reduce a remaining offset to the desired reference during flat-top.
An important challenge for real-time control of the q profile is to obtain a reliable estimate of the q profile.Improved q profile estimates can be obtained by coupling real-time equilibrium reconstruction with a dynamic state observer for the internal plasma profiles that combines transport modeling with the measurements available in real-time, as demonstrated for TCV in [34].A real-time RAPTOR dynamic state observer including evaluation of heating and current drive sources and real-time polarimetry measurements, as proposed for AUG in [64], would enable application of the real-time control techniques described above for the further development of AUG advanced scenarios.
If the pedestal is well-diagnosed, a linearized model around the nominal non-linear RAPTOR simulation could be leveraged to rapidly calculate the change in EC power and deposition radii required to maintain q min > 1.5, accounting for the difference between the expected pedestal temperature T e ped and the measured trace.As an intermediate step towards realtime control, an iterative learning control (ILC) approach [65,66] would be well-suited to iteratively approach the reference q min .The ILC algorithm is similar to MPC and allows to update the actuator trajectories before launching the discharge, using a post-shot interpretative simulation of the previous discharge as proxy for the expected state evolution in the absence of actuator trajectory changes.Note that, like MPC, ILC applies quadratic programming and a linear plasma model, rather than solving a non-linear optimization problem.However, shot-to-shot variations and unexpected events limit the effectiveness of ILC.Eventually, one could control the ECCD deposition based on the real-time estimate of the effective pedestal height and the q profile, both consistently evolved within a dynamic state observer.

Conclusion
Fast inter-discharge optimization has allowed to develop a stable and reproducible AT early heating scenario, reaching a relaxed elevated q profile and stationary high-performance confinement conditions early during flat-top, as experimentally tested in AUG experiments.
To efficiently maintains an elevated q profile, even at large plasma current I p = 1 MA, an AT scenario has been developed with central counter-ECCD.While the counter-ECCD scheme will not be applied in a fusion reactor, it allows to obtain an elevated q profile in a regime approaching reactor-relevant dimensionless parameters.The aim of this scenario includes the identifications of the physics drivers for improved confinement and the determination of operational limits, e.g.β N limits for the onset of (resistive) MHD modes.
Early heating of the plasma allows to slow down current diffusion, hence prolonging the broad internal current density profile characteristic of the early ramp-up phase.As a result, an elevated q profile can be achieved early during the flattop phase, allowing to study the physics of the AT scenario over multiple resistive times τ R on present-day devices.Early heating scenarios are of even greater importance for advanced scenarios on ITER and DEMO, where the much larger τ R is prohibitive for a late heating scenario.However, the development of early heating scenarios is often hampered by the appearance of tearing modes.
A set of reduced physics models in RAPTOR was successfully applied to model and optimize the ramp-up phase of the 1 MA counter-ECCD AT scenario on AUG.Solving electron heat and current diffusion in RAPTOR with ad-hoc formulas for the turbulent heat transport and the EC current drive efficiency was found to robustly recover the coupled dynamics of T e and q, while maintaining the small number of tuning parameters fixed over an extensive set of AT discharges.The set-up relies on the IDE tool-chain to obtain equilibrium geometry data, neutral beam heating and current drive deposition profiles (RABBIT [22]) and electron cyclotron deposition radii and widths for the individual EC sources (TORBEAM [21]).To allow for a pre-shot evaluation of the T e | ρ=0.8 , used as a boundary condition in RAPTOR, a new scaling law was derived for the pedestal top electron pressure, with line average density, heating power and I p as input variables, leading to a good scaling law fit n e T e | ρ=0.8 = 0.51n 0.82 el 10 19 m −3 (P oh MW + P aux MW ) 0.53 I 1.71 p MA with R 2 = 0.96 (for the linear fit on the logarithmic scale), when using a data set of discharges operating within the AT operating regime.Assuming the pedestal density remains unchanged with respect to a previous discharge, when fueling and NBI changes are minor (with control of n e sep by proxy of the neutral particle flux), the T e boundary condition at ρ = 0.8 can be readily calculated.While the predictive simulator should be sufficiently reliable to predict the main effects of the various actuators and the main inter-dependencies of temperature and q profile, detailed features are of lesser importance for the present work.
Rapid inter-discharge simulations and optimizations have guided the development of a reliable and reproducible early heating strategy.In a first step, post-shot RAPTOR simulations have allowed to understand the onset of 3/2 tearing modes, triggered during the high-β phase and having a significant detrimental effect on confinement.Our simulations indicated that the NBI heating was initiated too early during ramp-up, leading to the formation of an off-axis minimum of the q profile while q min > 1.5.The subsequent decrease of the q profile eventually leads to q min ∼ 1.5, with low magnetic shear at the rational surface, which was found to correlate with the onset of the corresponding tearing mode.
Since the q profile dynamics are highly sensitive to the timing of heating and fueling actuators, finding the onset timing of the various actuators gives rise to a delicate balancing act, benefiting from model-based pre-shot predictions.In a second step, the neutral beam heating onset time has been optimized, by a time interval found in simulations.Discharges with the optimized NBI trace were successful in avoiding these modes.Maintaining the relatively cold L-mode plasma for a longer time, q decreases below q = 1.5 earlier, at low β pol and before the q profile is strongly shaped.When the off-axis q min is formed, q min < 1.5, allowing to maintain non-zero shear through the rational 3/2 surface.
To reach the stationary state of the elevated q profile early during the pulse, a robust access strategy needs to be complemented with optimized flat-top heating and current drive deposition, such that the sum of ohmic and auxiliary current are consistent with a relaxed current diffusion solution.In the third step of our ramp-up studies, a non-linear optimization scheme has allowed to find the optimum ECCD deposition radii to maintain an elevated q min stationary, both for q min < 1.5 and q min > 1.5, whilst reaching the stationary high-performance phase early during flat-top.The proposed NBI trace and ECCD radii have been tested in experiment and resulted in a series of shots without confinementdegrading NTMs during ramp-up, reaching a stationary q min > 1.35 around 1.5 s, immediately after ramp-up, leaving a long flat-top time window for physics studies, e.g.regarding the improved ion heat confinement.
However, further experiments, aiming for q min > 1.5, have highlighted limitations of the present feedforward, modelbased control approach in the presence of shot-to-shot variations that are not included in the applied model.Interestingly, pedestal confinement was found to be an important physics driver for the central q profile.The key impact of T e | ρ=0.8 has been demonstrated: the redistribution of the plasma current due to an increase or reduction of the pedestal temperature significantly impacts the central q profile, as a transport effect leads to small changes of j par at mid radius and since the enclosed plasma current at small radii is relatively small for this scenario.Extending this work with a more advanced pedestal model, like proposed in [67,68], would allow extrapolation to different scenarios, as well as different machines, hence generalizing the applicability.
Finally, real-time model-based control was proposed as an interesting avenue to complement the feedforward optimization techniques applied here and compensate model-reality mismatches.If the pedestal is well-diagnosed, a linearized model around the nominal non-linear RAPTOR simulation could be leveraged to rapidly calculate the change in ECCD radii in real-time, in order to maintain q min > 1.5.

Appendix B. Feasibility study for replacement current-drive NBI source for ASDEX Upgrade advanced scenario
The beamline geometry of the neutral beam injectors on ASDEX Upgrade allows to distinguish two current drive sources that are best suited to drive off-axis current, namely sources 6 and 7.The strong off-axis current driven by source 6 and 7 is important for the advanced scenario described in this paper, which aims for an off-axis peak in the current density profile.
By the end of 2021, a leak in the ion dump impeded further use of NBI source 7 in the ongoing AUG campaign.To assess whether the scenario could still be run with the remaining sources, the modeling frameworks in ASTRA and RAPTOR were used to make a quantitative assessment [52].
To obtain source profiles for the individual NBI sources, interpretative ASTRA [50,51] runs have been performed for the various NBI sources, applying the RABBIT code [22] to calculate the deposition profiles.Figure B1 shows the current density contributions of the individual NBI sources at 3.5 s.
Summing different combinations of these source profiles, assuming the absence of non-linear effects, the impact of the loss of an individual source can rapidly be explored in RAPTOR.Firstly, it has been checked that the RAPTOR simulation using the NBI sources present in 39342 (NBI 3, 6 and 7) properly recovers the RAPTOR post-shot simulation making use of the total NBI deposition profiles from IDE (not shown here).
The next step is to assess the impact of the loss of NBI source 7 on the ability to maintain a stationary elevated q profile.A RAPTOR simulation has been performed with a newly proposed NBI heating mix, combining NBI sources 2, 6, 8, while maintaining the original mix NBI 3, 6 before 1.27 s.For this NBI heating mix, only minor differences in q profile evolution are observed with respect to the initial NBI set-up, as illustrated in figure B2.For comparison, we have also simulated the impact in case of additional loss of NBI source 6, the source with the largest off-axis current drive.Note that the original set-up relies on source 6 from 0.7 s onwards.Replacing source 6 by source 4, we can see that the impact of losing source 6 would have been much more significant compared to the loss of source 7. The reduction in the neutral beam driven current I nb leads to a significant drop in the q profile.It is however interesting to note that even in this case maintaining q min > 1.2 seems feasible.Table C1.Parameters of the AUG advanced scenario discharges discussed in this paper.MHD activity during t < 2.5 s is summarized, with boldface indicating a clear impact on core Te, and underline indicating a subsequent disruption.The timing is chosen to select the best performance phase (for 40188 at 1.6 s, the indicated performance is achieved transiently due to a β pol overshoot, and later recovered after a pre-programmed β pol ramp).

Figure 1 .
Figure 1.An illustration is given of the post-shot processing algorithm that is used to extract EC deposition radii and widths for individual EC sources from TORBEAM, applied for AUG discharge 39342.(a) Injected power per EC source; (b) deposition radius ρ dep for each EC source, based on a Gaussian fit of the TORBEAM heat deposition profile; (c) deposition width w dep for each EC source, based on a Gaussian fit of the TORBEAM heat deposition profile; (d) total driven current by all EC sources according to TORBEAM (blue solid) and RAPTOR (red dashed); (e)-(f ) at 2 s (e) and 4 s (f ), the TORBEAM electron heat deposition profile for each individual EC source is given in the color of the corresponding trace in (a)-(c); the respective Gaussian fits are shown as black dashed curves; the total EC electron heat source profile resulting in RAPTOR is shown as a shaded grey area.Note that a minimum ρ dep min = 0.08 is imposed to avoid the formation of a current hole for ρ ∼ 0.

Figure 2 .
Figure 2.For a range of different AUG advanced scenario discharges, the total EC driven current calculated by TORBEAM (blue solid line) is compared with the current evaluated with the current drive efficiency formula inside the RAPTOR code (red dashed line), with a fixed current drive efficiency parameter c cd = −11 (fixed in time and same for all discharges).

6 )Figures 5
Figures 5 and 6 illustrate the quality of the scaling law fit.The linear fit of the logarithmic quantities has a coefficient of determination R 2 = 1 − Σ(y − y scaling ) 2 /Σ(y − ȳ) 2 = 0.96 (with ȳ the mean of y) and a root mean square error of RMSE = 0.22 (normalized by ȳ).For n e T e | ρ=0.8 = exp(y), this metric yields R 2 exp = 1 − Σ(exp(y) − exp(y scaling )) 2 /Σ(exp(y) − exp(y)) 2 = 0.89 and RMSE exp = 0.16 (normalized by exp(y)).In figure5, the IDA data points used to derive the scaling law are shown: the x-axis indicates the IDA measurements of n e T e | ρ=0.8 for all IDA time points; the y-axis indicates the corresponding scaling law prediction.Clustering of the data points along the black solid line gives an indication of the quality of the fit.

Figure 4 .
Figure 4. RAPTOR post-shot simulation vs IDA/IDE reconstruction for advanced scenario on AUG (discharge 40029) with counter-ECCD.The Te evolution is accurately captured, until a 3/2 mode degrades confinement after 3 s.(a) Ip (blue solid), Ip of CHEASE equilibria used for RAPTOR metrics (black dots), NBI power P nb (green dotted), EC power Pec (red dashed); (b) Te at ρ = [0.20.4 0.6 0.8],for IDA (blue solid) and RAPTOR (red dashed); (c) q at ρ = [0.40.6 0.8] and q 95 , for IDE (blue solid) and RAPTOR (red dashed), q 95 of CHEASE equilibria (black dots); (d) Te profile at various times, for IDA (blue solid) and RAPTOR (red dashed), the IDA error bar is indicated with a blue shaded area; (d) q profile at various times, for IDE (blue solid) and RAPTOR (red dashed), the IDE error bar is indicated with a blue shaded area.

Figure 5 .
Figure 5.A scaling law for neTe| ρ=0.8 is presented, with input variables Ip, Paux + P oh and n el , based on IDA data available for 8 AUG advanced scenario discharges (considering the full discharge, sampling the time trace neTe| ρ=0.8 obtained from IDA, with a time step of 10 ms).The colored dots show the available data points based on which the scaling law presented in this section has been derived: for the different data points, the scaling law prediction is shown on the y-axis, versus the IDA measurement on the x-axis.Data points close to the black solid line indicate good agreement between the IDA reconstruction and the scaling law prediction.The black dash-dotted lines indicate an over-/under-prediction of 20%.

Figure 6 .
Figure 6.The scaling law presented in figure 5 is applied to a set of 6 AUG advanced scenario discharges: based on the input variables Ip, Paux + P oh and n el , a prediction of the neTe| ρ=0.8 trace is provided.To illustrate the performance we compare the scaling law prediction (black) to the IDA data (color).

Figure 7 .
Figure 7. Pre-and post-shot RAPTOR simulations of discharge 40398 are compared.The blue solid lines show the experimental reconstruction from IDA/IDE.In red dash-dotted lines, the pre-shot simulation is presented: the Te| ρ=0.8 boundary condition is set with a scaling law derived based on earlier discharges; the density is adapted from the IDA reconstruction of discharge 39342, taking into account a delay due to the later onset of NBI.The post-shot simulation, utilizing the IDA data for Te| ρ=0.8 and ne, is shown in green dashed lines.(a) Te traces at ρ = [0.20.4 0.6 0.8]; (b) ne traces at ρ = 0 and ρ = 0.8 (the black dotted lines show the IDA reconstruction for discharge 39342); (c) q traces at ρ = [0.40.6 0.8].The vertical dotted line indicates the onset time of NBI.

Figure 8 .
Figure 8.Comparison of RAPTOR post-shot simulation for late heating discharge 36087 (blue lines) and early heating discharge 39342 (red lines).(a) Ip (black dotted), neutral beam power P nb ; (b) EC power P ec ; (c) minimum of q for ρ > 0.25; (d)-(e) q profile at 1.2 s (solid line), 1.8 s (dashed lines) and 4.5 s (dash-dotted lines), for 36087 and 39342.Note that for the late heating discharge 36087 the off-axis minimum of q is formed while locally q < 1.5, while for the late heating discharge 39342 the off-axis minimum of q is formed while locally q > 1.5, leading to a time point after 1.8 s when q = 3/2 and the local magnetic shear is zero.

Figure 11 .
Figure 11.Three RAPTOR simulations: (1) a post-shot simulation of the early heating attempt 39342, with NBI onset time t nb = 0.5 s in blue solid lines; (2) a predictive simulation assessing the impact of an NBI onset delay of 200 ms in red dash-dotted lines; (3) a post-shot simulation of discharge 40192 implementing the NBI delay proposed by RAPTOR in green dotted lines.(a) NBI power P nb ; (b) EC power Pec and Ip (black dash-dotted line); (c) neutral beam driven current I nb ; (d) Te at ρ = 0.8; (e) n e0 ; (f ) minimum of q for ρ > 0.25; (g)-(i) q profile at 1.2 s, 1.6 s and 2.0 s.For the plots with time traces, the delay of the NBI has been indicated with a grey shaded area and the end of the pre-programmed β ramp has been indicated with a vertical dotted line.

Figure 12 .
Figure 12.Magnetic spectrograms for discharge 40030, quiescent prior to the premature termination of the discharge due to a false locked mode alarm.

Figure 13 .
Figure 13.Magnetic spectrograms for discharge 40192.A 3/2 mode is triggered at 1.1 s, disappearing around 1.8 s, without a clear impact on confinement.Later during the flat-top phase, sporadic occurrences of (short) time windows with 3/2 modes (likely triggered by a growing n = 3 mode) indicate that the scenario continues to be metastable regarding NTMs.

Figure 14 .
Figure 14.Results of forcing fast access to stationary q min > 1.35, by automatically optimizing EC deposition radii.For the initial (blue solid) and optimized (red dashed) state evolution in RAPTOR: (a) q profile at t = 3.5 s (the grey shaded area highlights q < 1.35); (b) time evolution of the minimum of q for ρ > 0.25 (the grey shaded area highlights t > 1.5 s, where the cost function is active); (c) EC power deposition profile; (d) EC power and Te| ρ=0.8 boundary condition (from scaling law).

Figure 15 .
Figure 15.The optimum state evolution forcing q min > 1.35, presented in figure 14, is shown in red dashed lines.Sensitivity of the modeled state evolution to a 20% increase or reduction of the Te| ρ=0.8 assumption is shown in respectively black dotted and black dash-dotted lines.The post-shot simulation of discharge 40398 that implements the actuator traces proposed by the optimizer is shown in green solid lines: (a) q profile at 3.5 s; (b) minimum of q for ρ > 0.25; (c) EC power deposition profile; (d) Te| ρ=0.8 .
(i) late heating discharge 36087; (ii) previous attempt for an early heating scenario 39342, triggering a confinement-degrading NTM around 1.7 s; (iii) optimized early heating discharge 40398.

Figure 18 .
Figure 18.The time traces of the safety factor q at various ρ = [0.350.45 . . .0.95], for the IDE reconstruction (blue solid lines) and the RAPTOR post-shot simulation (red dashed lines), for the late heating discharge 36087 on the left-hand side and for the optimized early heating discharge 40398 on the right-hand side.For discharge 40398, also an ASTRA post-shot simulation is included, in green dotted lines.

Figure 19 .
Figure 19.Results of forcing fast access to stationary q min > 1.55, by automatically optimizing EC deposition radii.For the initial (blue solid) and optimized (red dashed) state evolution in RAPTOR: (a) q profile at t = 3.5 s (the grey shaded area highlights q < 1.55); (b) time evolution of the minimum of q for ρ > 0.25 (the grey shaded area highlights t > 1.5 s, where the cost function is active); (c) EC power deposition profile; (d) EC power and Te| ρ=0.8 boundary condition (from scaling law).

Figure 20 .
Figure 20.The optimum state evolution forcing q min > 1.55, presented in figure 19, is shown in red dashed lines.Sensitivity of the modeled state evolution to a 20% increase or reduction of the Te| ρ=0.8 assumption is shown in respectively black dotted and black dash-dotted lines.The post-shot simulation of discharge 41102 that implements the actuator traces proposed by the optimizer is shown in green solid lines: (a) q profile at 3.5 s; (b) minimum of q for ρ > 0.25; (c) EC power deposition profile; (d) Te| ρ=0.8 .

Figure 21 .
Figure 21.Left-hand panel: (a) q profile at 3.5 s of the RAPTOR simulation of the optimum state evolution (red dashed) forcing q min > 1.55, presented in figure 19; post-shot simulation of discharge 41102 that implements the actuator traces proposed by the optimizer (blue); post-shot simulation of discharge 41102 with T e ped taken from the pre-shot simulation; post-shot simulation of discharge 41102 with T e ped , ne and NBI deposition taken from the pre-shot simulation.Right hand panels, for the optimum state evolution (red dashed lines) and the post-shot simulation of discharge 41102 (blue solid lines): (b) Te; (c) ne; (d) neutral beam current density j nb ; (e) bootstrap current density j bs , at 3.5 s.

Figure 23 .
Figure 23.For the RAPTOR optimum state evolution forcing q min > 1.55 (magenta dashed line) obtained in section 4.3, the sensitivity of the q profile with respect to variations in the EC deposition radii is studied.The profiles corresponding to 50 perturbed simulations (adding a random number to each of the optimization variables) are plotted in light grey, together with a mean, an upper and a lower standard deviation for the resulting q profile.(a) q profile; (b) EC driven current density; (c) integrated enclosed EC driven current.

Figure 24 .
Figure 24.The optimum state evolution forcing q min > 1.55, presented in figure 19, is shown in red dashed lines.Sensitivity of the modeled state evolution to a 20% increase or reduction of the Te| ρ=0.8 assumption is shown in respectively black dotted and black dash-dotted lines, illustrating the coupling between the current in the pedestal and the central q profile.(a) Te| ρ=0.8 ; (b) q profile at 3.5 s; (c) plasma current integrated over ρ ∈ [0.8 1]; (d) Te at 3.5 s; (e) plasma current integrated over ρ ∈ [0 0.4]; (f ) parallel current density jpar at 3.5 s.

Figure 25 .
Figure 25.Inductive and non-inductive contributions to the current density and the enclosed plasma current at t = 3.5 s for the optimum RAPTOR simulation forcing q min > 1.55, presented in figure 19.(a) contributions to jpar and total jpar; (b) contributions to enclosed plasma current; (c) q profile.

Figure A1 .
Figure A1.Data of four RAPTOR post-shot simulations are presented to illustrate the application of formula equation (A.1) to estimate the ohmic power P oh MW , with plasma current I p MA and total auxiliary heating P aux MW as input variables.Left hand side: ohmic power P oh MW versus I p MA for all time points of post-shot RAPTOR simulations for four AUG discharges, during the ohmic heating phase (crosses) and the auxiliary heating phase (dots).During the ohmic heating phase, a linear relation P oh MW = c Ip I p MA captures well the dependency of the ohmic power predicted in RAPTOR on the plasma current (indicated with black dash-dotted line).Right hand side: for these four shots, the ohmic power time trace found by evaluating formula equation (A.1) is compared to the data points in the RAPTOR simulation.

Figure B1 .
Figure B1.Individual current density contributions of the different NBI sources, calculated by RABBIT in an interpretative ASTRA run of discharge 39342 at 3.5 s.NBI sources 6 and 7 provide the most significant off-axis current drive contribution.

Figure B2 .
Figure B2.Using the deposition profiles of individual neutral beams, as obtained from interpretative ASTRA-RABBIT simulations, RAPTOR simulations for different combinations of NBI sources have been run, testing the impact on the q profile evolution.(a) neutral beam driven current I nb ; (b) neutral beam driven current density j nb at 3.5 s; (c) q at ρ = 0.4; (d) q profile at 3.5 s.

Figure C2 .
Figure C2.The time traces H 98y,2 (left-hand panels), β pol (middle panels) and β N (right-hand panels) for the late heating discharges 36087 and the early heating discharges 39342, 40188, 40192, 40398, 41102 and 41400.For H 98y,2 , the TTR post-processed diagnostic signal is used, applying a RABBIT simulation to subtract the fast ion energy.

Table 1 .
Setting of tuning parameters in the RAPTOR transport and Gaussian EC model.Note that these model parameters are kept constant for all simulations performed in this paper.