Real-time control of NBI fast ions, current-drive and heating properties

Figure 1 shows time-traces of the first experiment (#38308) to control the NBI current-drive (NBCD). The scenario is based on a low-density, highly non-inductive scenario. We have programmed an NBCD trajectory which features ramps and a sharp step at the end of the plasma discharge. In between the ramps, there is a longer constant phase, where we want to demonstrate that the controller is able to keep the NBCD constant even under changing background plasma conditions. The latter are emulated by a pre-programmed drop of ECRH power at 4.2 s. Our first attempt, #38308, features some oscillations of the actually achieved NBCD trajectory. This is due to non-optimal controller gains (see table 1), because the beams drove more current than expected. A second attempt with optimized gains (#38310, figure 2) shows then a very good behavior and the target trajectory is met well. This could be automated in the future by calculating the current-drive efficiency with RABBIT, and scaling the controller gains accordingly.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To cross-check the real-time results we have carried out offline simulations both with RABBIT (in time-dependent mode) and TRANSP-NUBEAM. The offline inputs are more precise as there are much more diagnostics available than in real-time, and even the NBI power signals are less accurate in real-time (see section 3.2).

It should also be noted that the current real-time implementation of RABBIT in the DCS makes the assumption $\left \langle {1}/{R} \right\rangle \approx \sqrt{\left \langle {1}/{R^2} \right\rangle}$, where $\left\langle \ldots \right\rangle$ means flux-surface average, because the current version of the real-time equilibrium code JANET provides only $\left \langle {1}/{R^2} \right\rangle$ and not $\left \langle {1}/{R} \right\rangle$. $\left \langle {1}/{R} \right\rangle$ is needed to calculate the area elements [13] which are needed to integrate the driven current density into an absolute current. Above approximation hence only affects the current drive output of RABBIT and it should be a quite good approximation. It will be removed with the next upgrade of JANET.

The agreement between offline-RABBIT and TRANSP-NUBEAM is excellent, underlining the high fidelity of RABBIT. The offline calculations are done in time-dependent and the real-time calculation is done in steady-state mode. One can see that the steady-state assumption does not cause strong deviations: A delayed response to the large jump at 7 s is visible, but otherwise, the qualitative temporal behavior is very similar. Overall, both offline calculations are a bit lower than the real-time result (here mostly due to inaccuracies in the real-time T_e profile). It is planned to revisit the real-time T_e data evaluation and calibration. Since the raw data is coming from the same diagnostic (ECE) it should be possible to further improve the real-time data quality and hence the agreement. But also in the current state, the agreement is still reasonable and hence we can conclude that control of NBCD has been demonstrated successfully.

The fast-ion content has a strong influence on the plasma behavior. Too many fast-ions can trigger MHD instabilities such as Alfvén eigenmodes [4, 5]. Control of these Alfvén eigenmodes, e.g. by reducing their drive by lowered NBI power, is a research topic of its own [21], and the neutron emission calculation by RABBIT could help to determine how harmful the Alfvén eigenmodes are, e.g. by estimating how strong the measured neutrons fall short with respect to the expectation based on a neo-classical fast-ion distribution. On the other side, fast ions can stabilize turbulence [6] if their density and temperature gradient is appropriate. Advanced control schemes might help to tailor the fast-ion profiles to achieve and sustain these high confinement regimes more readily (at least in present-day machines, where most fast ions stem from external heating sources).

For a first controller test, several scalar quantities would be possible, e.g. the volume integrated fast-ion density (potentially also inside a certain radius). Here we chose the fast-ion stored energy, or more precisely its ratio with respect to the total plasma energy, i.e. $W_\text{fi} / W_\text{mhd}$. This has the advantage of being a dimensionless quantification of the fast-ion content, while perhaps being more descriptive than e.g. the ratio of particles. W_mhd is estimated by the real-time equilibrium reconstruction based on function parametrization (FP) [22, 23], since JANET does currently not output W_mhd.

Figure 3 shows ASDEX Upgrade discharge #39654. It features a similar requested trajectory with ramps, steps and a constant phase which is well achieved by the controller. However, the response of the controller seems to be a bit slower than in the second NBCD case (figure 2). This is most visible during the large step at 7 s, but also e.g. during the downward ramp (e.g. around 6 s). This suggests to use slightly more aggressive (i.e. larger) controller gains in the future (see table 1 for the used values).

Zoom In Zoom Out Reset image size

An offline RABBIT calculation with more precise inputs (i.e. W_mhd and equilibrium from CLISTE, T_e and n_e from IDA integrated data analysis [24]) shows a small deviation but overall good agreement. In this case, the deviation seems to originate from missing CXRS real-time measurements: An offline simulation without CXRS data, making the same assumptions of $T_\text{e} = T_\text{i}$ and zero rotation as in real-time, shows excellent agreement with the real-time simulation. During most of this discharge, core T_i is less than half of T_e, such that our real-time assumption of $T_\text{e} = T_\text{i}$ is not very good. Nevertheless, since T_i has only a very weak effect on the fast-ion distribution function, the deviations are still quite small and the control experiment can be considered successful.

The NBI-induced fast-ions slow down in collisions with the background plasma and heat it accordingly. Whether these collisions are happening with electrons or ions determines how strong the electron and ion heating of NBI is. RABBIT calculates radial profiles of both heating densities. For the control experiments, we focus on the volume-integrated heating of ions.

The ion heating controller has been tested for the first time piggy-back at the end of discharges #38134 and #38137. In the first, two power levels of ion heating (with a step) have been programmed into the discharge program and in the second it is a ion heating power ramp. As can be seen in figure 4, both for the ion heating power steps and ramp the target values are matched very well.

Zoom In Zoom Out Reset image size

In discharge #38164, we have also tested steeper ramps and a phase where the controller should keep the ion heating constant, while the background plasma changes. As in previously described experiments, this is done by lowering the pre-programmed ECRH power (here at 6.7 s). In figure 5(b) one can see that the real-time measured electron temperature T_e clearly drops during that time window. T_e has a very direct impact on the distribution of ion and electron heating, because the critical energy E_c, above which fast ions collide dominantly with electrons instead of thermal ions, is given by $E_\text{c} \approx 18.6 T_\text{e}$ in a D plasma with D beams [25–27]. In figure 5(a) we can clearly see that the controller reacts to the dropped T_e by increasing P_NBI, thus keeping the ion heating constant as requested.

Zoom In Zoom Out Reset image size

The controller output u_k (in the kth time-step) will then be a vector of both heating powers. To this end we use the general form of a MIMO (Multi-Input-Multi-Output) controller, which can be written as:

$$\begin{align} u_k & = C x_k + D e_k + u_{\text{ff}, k} \end{align} \tag{ 1 }$$

$$\begin{align} x_{k+1} & = A x_k + B e_k + K \left(u_{k,\text{sat}} - u_k\right).\end{align} \tag{ 2 }$$

Here, e_k is the error vector in the kth time-step, $u_{\text{ff}, k}$ is the feed-forward input vector. x_k is the controller state vector, which has no direct physical meaning but can be defined such to match the control goals. A, B, C, D and K are matrices with constants.

In this general form, the controller had already been implemented in the DCS. The values of the vectors and matrices can be configured flexibly by only modifying an xml configuration file, and new controllers of the same type with other parameters can be added flexibly as well.

In our case, we choose the vectors and matrices to be:

$$\begin{align} \begin{pmatrix} P_\text{NBI}^{k} \\[2pt] P_\text{ECRH}^{k} \end{pmatrix} = & \begin{pmatrix} 1 & 0 \\[2pt] 0 & 0 \end{pmatrix} \begin{pmatrix} x_\text{NBI}^{k} \\[2pt] x_\text{zero}^{k} \end{pmatrix} + \begin{pmatrix} K_\text{p} & 0 \\[2pt] 0 & 1 \end{pmatrix} e_k + u_{\text{ff}, k} \end{align} \tag{ 3 }$$

$$\begin{align} \begin{pmatrix} x_\text{NBI}^{k+1} \\[2pt] x_\text{zero}^{k+1} \end{pmatrix} & = \begin{pmatrix} 1 & 0 \\[2pt] 0 & 0 \end{pmatrix} \begin{pmatrix} x_\text{NBI}^{k} \\[2pt] x_\text{zero}^{k} \end{pmatrix} + \begin{pmatrix} K_\text{i} \cdot \Delta t & 0 \\[2pt] 0 & 0 \end{pmatrix} e_k \nonumber\\ & \quad + \begin{pmatrix} 0.99 & 0 \\[2pt] 0 & 0 \end{pmatrix} \left(u_{k\text{,sat}} - u_k\right) \end{align} \tag{ 4 }$$

$$\begin{align} \text{with: } e_k = \begin{pmatrix} P_\text{ion, reference}^k - P_\text{ion, RABBIT}^k \\[4pt] P_\text{total, reference}^k - P_\text{NBI, actual}^k \end{pmatrix} . \end{align} \tag{ 5 }$$

This means that effectively we setup a conventional PI controller for P_NBI in the first row, while the second row simply evaluates the boundary condition $P_\text{ECRH} = P_\text{total, reference} - P_\text{NBI, actual}$. $P_\text{NBI, actual}$ can be the raw injected NBI power signal $P_\text{NBI, inj}$, or the NBI heating power as calculated by RABBIT $P_\text{NBI, heat}$. The latter accounts for NBI power losses (e.g. due to shine-through and promptly unconfined orbits) to ensure that the total effective heating power (i.e. after subtraction of losses) is kept constant. We have tested both settings in the following. Δt is equal to the time step, i.e. the cycle time of the DCS which is 1.5 ms in all discharges of this publication. The last term in equation (4) is there to counter-act integral windup in the case of saturating actuators (only then $u_{k,\text{sat}}-u_k\neq0$, since $u_{k,\text{sat}} = \text{min}(\text{max}(u_k, 0), \text{max. Power})$), but it is irrelevant for our purposes because it is always zero during the experiments in this publication.

To test this new control scheme, we chose a very low density scenario (core densities around $2 \cdot 10^{19} \,\text{m}^{-3}$) which was designed to study LH transitions. This scenario is interesting, because we can connect our control experiments with physics of the LH transition, namely the established theory that the LH transition is triggered by a threshold in (surface integrated) ion heat flux Q_i, rather than total heating power [7–9].

The (volumetric) ion heat-flux density q_i is given by the ion heating density plus the equipartition power density (i.e. collisional heat transfer between electrons and ions):

$$\begin{equation} q_\text{i} = p_\text{ion} + 0.00246 \cdot \sum_j n_j \frac{Z_j^2}{A_j} \cdot n_\text{e} \text{ln} \Lambda \cdot \frac{T_\text{e} - T_j}{T_\text{e}^{3/2}} .\end{equation} \tag{ 6 }$$

Hereby, temperatures are in keV, densities in $10^{19}\,\text{m}^{-3}$, $\text{ln} \Lambda$ is the Coulomb logarithm, Z_j and A_j are charge and mass numbers of the ion species j, respectively, and the result is in MW m⁻³ [28–30]. Q_i and ion heating are then volume integrals of q_i and p_ion, respectively, and are both measured in Watts.

In low density plasmas, the coupling between electrons and ions (i.e. the equipartition term) is weak and hence one can modify the ion heat flux by changing the ion heating. In a high density plasma, this would not be possible because even for pure electron heating there would be a considerable ion heat-flux due to collisional heat transfer from the electrons to the ions.

Hence, it should be possible to design discharges that start with 100% electron heating still in L-mode and will transit into H-mode while ramping up the ion heating ratio at constant total power. A few setup discharges have been carried out to identify the parameter window where this is possible.

Figure 6 shows the first test of such a discharge. The plasma is initially heated with two gyrotrons, staying in L-mode. At 3.0 s, an ion heating ramp is programmed while keeping the total heating power constant and the LH transition is triggered at 4.3 s. The controller achieves the ion heating ramp almost perfectly. The total power is, however not fully constant.

Zoom In Zoom Out Reset image size

In phase (b), this is simply due to saturation. We have demanded too much ion heating, and the ECRH power cannot become negative. The controller here decides to match the ion heating time-trace, sacrificing to keep the total power constant. The opposite behavior would be better and more stable, i.e. to just keep the total power constant with 100% NBI. This can be programmed easily in future experiments by clipping the NBI power to the maximum total power. Such a clipping feature was not available in the DCS at the time of the experiment, but has been included in the mean time. Even more elegantly, the same could also be achieved in the future by setting matrix K in equation (4) to

$$\begin{align} K = \begin{pmatrix} 0.99 & -1 \\ -1 & 0 \end{pmatrix}, \end{align} \tag{ 7 }$$

ensuring that saturation in ECRH power will be compensated with opposite sign in the NBI power request and vice versa.

In phase a) the total power is not constant due to uncertainties in the real-time NBI power signals during pulse-width modulation. As one can see in figure 7, the NBI system reacts with a delay to the modulation requests sent by DCS. The delay is longer for switch-on than for switch-off, resulting in a different duty cycle percentage. To solve this, new signals for all eight NBI beams with the actual NBI power have been routed and made available in the DCS.

Zoom In Zoom Out Reset image size

Figure 8 shows a subsequent repeat of the experiment. $P_\text{NBI, actual}$ is now set to $P_\text{NBI, RABBIT}$, where as in the previous attempt it was set to $P_\text{NBI, raw}$. The ion heating target is again matched very well, and the constancy of the sum of ECRH power and effective NBI heating power (i.e. after subtraction of RABBIT-calculated losses) is better than the constancy of ECRH + NBI power in the previous attempt. However, it is still not perfect and as one can see in figure 9. This time, it seems to be mainly caused by uncertainties in the real-time ECRH power. For example, the actuator manager of the DCS assumes that each gyrotron has the same power of 700 kW whereas in this discharge the real power of the gyrotrons was in the range of 500–800 kW. For example, gyrotron 2 and 3, which were active from 3–5.8 s delivered both 800 kW while the later replacements 4 and 5 delivered 800 kW and 500 kW, respectively (comp. figure 9(b)). These deviations from the constancy of power could be improved further in future experiments, by configuring the actuator manager with the expected power from each gyrotron.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To validate the real-time calculations of RABBIT we have again carried out offline calculations with the TRANSP-NUBEAM code. They are drawn with orange curves in figures 6(b) and 8(b) and the agreement with real-time RABBIT is good in #38307 and excellent in #39653 (which had the improved NBI power signals). This indicates that we controlled the ion heating correctly and that the improved NBI power signals have also improved the controller accuracy. In summary, control of the ion heating ratio at constant total heating power was demonstrated successfully.

Now we want to move one step further and try to control the ion heat-flux Q_i. As mentioned in equation (6), it is given by the sum of the ion heating and the equipartition term. This control scheme is more challenging, because the latter term is rather self-organized by the plasma and cannot be controlled from outside.

For the previous discharges, we have calculated Q_i already and plotted it with yellow curves in figures 6 and 8. One can see that while controlling the ion heating, we have also modified the ion heat flux in a similar fashion: Surely, Q_i cannot go to zero, because even at $P_\text{ion} = 0$, there is a residual energy transfer from the electrons to the ions. But this Q_i baseline is low enough and eventually, the ion heating ramp coincides with a Q_i ramp. This means that the equipartition term is weak enough (due to the low density plasma scenario) to allow an attempt to control Q_i.

To evaluate the equipartition term, we need profiles of n_e, T_e and T_i. For the first two, measurements are readily available in real-time at ASDEX Upgrade with interferometry and electron cyclotron emission (respectively). The real-time profiles are then estimated by RAPTOR (taking into account the measurements and a transport model) [31]. For T_i, unfortunately no real-time diagnostic (e.g. CXRS) is available, and the RAPTOR assumption of $T_\text{i} = T_\text{e}$ is too simplified for our application. To overcome this, we estimate the energy of the thermal ions by taking the plasma energy W_mhd as determined by the real-time equilibrium function parametrization (FP) [22, 23] and subtracting all other plasma species:

$$\begin{align}\kern-7pt W_\text{ion} = W_\text{mhd} - \int \left(\tfrac 3 2 n_\text{e} T_\text{e} + \tfrac 3 2 w_{\text{fi}, \parallel} + \tfrac 3 4 w_{\text{fi}, \perp} + \,2 w_\text{rot}\right) \text{d} V .\end{align} \tag{ 8 }$$

Here, w stands for energy densities and the energy W is obtained by volume integration of the w-profiles. The parallel and perpendicular fast-ion energy densities $w_{\text{fi}, \parallel}$ and $w_{\text{fi}, \perp}$ are obtained by the RABBIT code and the rotational energy density w_rot is assumed simply zero (due to lack of real-time rotation measurements). This is a good approximation especially in cases of low NBI power (and hence torque input), e.g. in #39653 the value of $W_\text{rot} / W_\text{mhd}$ does not exceed 0.7%.

Then, a mock-up T_i profile is calculated to match the ion energy $W_\text{ion} = \int \tfrac 3 2 (n_\text{i,th} + n_\text{imp}) T_\text{i} \text{d} V$. Here, we assume one main ion species (i.e. deuterium) with thermal ion density $n_\text{i,th}$ and one impurity species (i.e. carbon) with density n_imp. The latter approach is justified because high-Z impurities (such as tungsten) do not contribute significantly to the plasma dilution (due to their very low concentrations), and carbon is a suitable representative for the various low-Z impurities. The mock-up T_i is created as follows:

$$\begin{equation} T_\text{i} = \begin{cases} T_\text{e} + \textrm{const} \cdot \left(\rho_\text{tor} - 0.8\right) / \left(-0.8\right) & \text{for}\,\,\, \rho_\text{tor} \lt 0.8 \\ T_\text{e} & \text{for} \,\,\,\rho_\text{tor} \unicode{x2A7E} 0.8 \end{cases} \end{equation} \tag{ 9 }$$

ρ_tor refers to the square-root of normalized toroidal flux. Outside of $\rho_\textrm{tor} = 0.8$, we assume $T_\text{i} = T_\text{e}$ and consequently zero equipartition, because otherwise, e.g. when using a simpler approach such as $T_\text{i} = \textrm{const} \cdot T_\text{e}$, the equipartition will get very high in this region and the Q_i value can become arbitrary as a consequence. Inside of $\rho_\textrm{tor} = 0.8$ we use T_e as a basis and subtract a linear function where the constant prefactor is calculated such to match W_ion.

The results of this approach for discharge #39656 are shown in figures 10 and 11. It can be seen that in this discharge, it works surprisingly well and values of core T_i are in very good agreement with offline CXRS data, capturing correctly even the large increase in T_i at 5.5 s. Also the profile shape estimates are reasonable, although in situations where $T_\text{i} \ll T_\text{e}$ a nonphysical reversal of the T_i-gradient appears due to our approach. Even though this has no strong implications for our application, this can be improved in future experiments.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Q_i is then estimated based on equation (6) and a subsequent volume integration (up to $\rho_\textrm{tor} = 1.0$). Hereby, we assume the Coulomb Logarithm, $\text{ln}\Lambda$, to be 16. In the sum over ion species j we consider one thermal ion species (deuterium) and one impurity species (carbon). The impurity density is given by the RAPTOR Z_eff estimate (which is hardwired to 1.5 due to lack of diagnostics—this is a reasonable assumption since AUG is a tungsten device). The thermal ion density is given by subtracting impurities and fast-ions (as calculated by RABBIT) from n_e. Figure 12 shows a comparison of the real-time Q_i estimate with an offline TRANSP calculation (i.e. including CXRS measurements for T_i) for discharge #39656 and the agreement is very good.

Zoom In Zoom Out Reset image size

Figure 12 shows results of the first experiments to control the ion heat flux Q_i. The programmed target trajectories and plasma scenario are similar as for the ion heating case. The controller worked quite well in the initial phase of the ramp. At the beginning (3 s), the controller just stays at 100% ECRH heating, because it cannot lower Q_i anyhow. After ≈3.5 s the controller reacts to the increasing target-Q_i by exchanging ECRH power with NBI power. This works quite well for about a second, when oscillations start to appear that get worse over time, and also lead once again to an over-shoot of NBI power, hence also violating the constancy of total power. This leads to a strong density increase which increases the equipartition term, and hence makes further Q_i control impossible.

We make the following suggestions for improvement in subsequent experiments: First, one should request less Q_i. Second, the recently added feature of clipping the NBI power should be activated here too, to ensure constancy of total power also in situations where the controller is at 100% NBI heating. This would prohibit also the density increase which leads to a plasma state where Q_i control is no longer possible. To mitigate the oscillations, one could thirdly implement low-pass filters on the quantities entering into the equipartition term (i.e. T_e, T_i and n_e.) In our analysis we could not identify a single source for the oscillations. T_i does have the strongest oscillations in the phase around 4.6 s, but e.g. also T_e data shows oscillations, which could partly be induced by the controller oscillating between high ECRH and high NBI heating. Lastly, improvements in the T_i estimate would be very beneficial. In attempts to implement the first two of the above suggestions and repeat the experiment, the T_i estimate was less accurate than in the example above (without changing anything about the T_i estimate). Hence it is currently not in a very reliable state, and it is unclear if the current method can be robustified, because W_mhd comes with uncertainties, which get worse as we subtract all other plasma species from it: We essentially subtract two large numbers from each other, and the remainder determines then T_i, consequently with a high uncertainty. The ideal solution to this issue would be a real measurement of T_i in real time, e.g. through real-time CXRS.

In conclusion we can however state that the proof of principle of a Q_i controller was already successful. With additional work, it could be robustified and made ready for experimental applications.