A hybrid physics/data-driven logic to detect, classify, and predict anomalies and disruptions in tokamak plasmas

The amplitude of the locked mode has traditionally been a workhorse for triggering mitigation actions in many devices and on JET in particular (see [16–18, 50] and references therein). In this work, the detection of anomalous locked mode signals is implemented with a new approach that combines physics understanding and experimental data. Indeed, at beginning of a new experimental campaign or of a new tokamak, there are not specific data and therefore the predictors must be trained with physical information provided by other tokamaks, other campaigns, theory and/or numerical information.

Large locked mode signals are correlated with the size of the magnetic islands due to small toroidal number (usually n = 1) macroscopic instabilities [1, 48]. It is now clear that when the velocity of these modes decreases in the laboratory frame of reference (finally locking to the wall), the amplitude of the radial perturbations increases exponentially and usually leads to a disruption. From a theoretical point of view, the threshold for the locked mode amplitude triggering a disruption is interpreted in terms of the Chirikov parameter [51]. However, in real time the correct computation of this quantity is quite challenging and in the last decades empirical models have been developed to understand what is the best mode locking threshold to be used in multimachine applications. A multimachine scaling in power law form for the mode locking, interpreting the various quantities in terms of the magnetic island maximum width, was formulated in [48]. Unfortunately, such a law does not perform very well for JET with the ITER like wall [47] and therefore alternatives have been proposed. In [47] a simple scaling law, derived with a machine learning based methodology, was proposed. It links the locked mode amplitude with the internal inductance. In this work, we have implemented this approach to find an empirical equation that correlates the normalised locked mode signal (the locked mode amplitude divided by the magnetic poloidal field at the geometrical minor radius, calculated through the plasma current, LM_N ) with the internal inductance (l_i ):

$$\begin{equation}{\text{L}}{{\text{M}}_{N,{\text{threshold}}}} = al_i^b.\end{equation} \tag{ 2.1 }$$

The equation has been particularised using the same campaigns of the cited article [47] with exactly the same setup. To this end, first the disruption probability is modelled with the sigmoid function reported in equation (2.2):

$$\begin{equation}{p_{{\text{LM}}}} = {e^{{d_1}\left( {{\text{L}}{{\text{M}}_N} - {\text{L}}{{\text{M}}_{N,{\text{ threshold}}}}} \right) + {d_2}}}/\left( {1 + {e^{{d_1}\left( {{\text{L}}{{\text{M}}_N} - {\text{L}}{{\text{M}}_{N,{\text{ threshold}}}}} \right) + {d_2}}}} \right).\end{equation} \tag{ 2.2 }$$

At the beginning of a new experimental campaign, no new data is available and therefore the locked mode predictor is trained to reproduce equation (2.2). The predictor consists of a simple feed forward neural network (3 hidden layers of 20 neurons each with tanh transfer functions plus one neuron output layer), which returns the probability (q_LM) of an anomalous locked mode amplitude. The cross-entropy loss function has been used to train the neural network:

$$\begin{equation}{\text{los}}{{\text{s}}_{{\text{physics}}}} = - {\mathop \sum \nolimits }{p_{{\text{LM}},1}}\text{ln} \left( {{q_{{\text{LM}},1}}} \right) - {\mathop \sum \nolimits }{p_{{\text{LM}},0}}\text{ln} \left( {{q_{{\text{LM}},0}}} \right)\end{equation} \tag{ 2.3 }$$

where the subscript 0 indicates a safe discharge, the subscript 1 a disruptive shot. Moreover, p_LM,0 is the expected probability that the mode locking is not anomalous (calculated as 1—p_LM and derived from equation (2.2); q_LM and q_LM,0 are the predicted probabilities (the neural network predicts q_LM and q_LM0 is 1—q_LM).

Once the new device or campaign starts, pulse by pulse, the training set is meant to be updated (see section 3) and the neural network retrained adding a new data-driven term to the loss function:

$$\begin{equation}{\text{loss}} = \frac{{\left( {{\text{los}}{{\text{s}}_{{\text{physics}}}} + \alpha {\text{ los}}{{\text{s}}_{{\text{data}}}}} \right)}}{{1 + \alpha }}\end{equation} \tag{ 2.4 }$$

where:

$$\begin{equation}{\text{los}}{{\text{s}}_{{\text{data}}}} = - \sum\limits_{{\text{class}} = 1}^{} \text{ln} \left( {{q_{{\text{LM}},1}}} \right) - \sum\limits_{{\text{class}} = 0}^{} \text{ln} \left( {{q_{{\text{LM}},0}}} \right).\end{equation} \tag{ 2.5 }$$

And α is a weight parameter. In this work, α has been set to a constant value equal to 1. However, it is possible to adjust its value as a function of the amount of the training set data and their reliability.

In terms of interpretation, this predictor returns the probability that the time slice has an anomalous disruptive mode locked amplitude. Since the locking to the wall of the magnetic perturbation is the last event before a disruption, p_LM can also be conceived as the probability of the plasma disrupting. This is due to the fact that in general, once a macroscopic mode has locked to the wall, there is not enough time to undertake any other action except mitigation. As discussed in the next subsections, this correspondence, between probability of an anomaly and probability of disruption, does not apply to the outputs of the other predictors for avoidance and prevention.

The two electron temperature anomalies considered in this work are EC and core hollowness. Both anomalies leave a clear signature in the electron temperature profiles. In the case of hollowness, the core electron temperature is not peaked in the centre but off-axis, becoming bimodal. In the case of EC, a strong decrease of the peripheral temperature is observed [12, 41, 42, 52]. The two indicators proposed in [52] to quantify these anomalies are the ones used in this work. They are briefly described in the following and they are calculated using the data of JET high resolution Thomson scattering.

Regarding the electron temperature hollowness (ETH), the 'Gaussian fitted hollowness' indicator has been implemented. The indicator is derived by fitting a bimodal symmetric Gaussian function to the electron temperature profile. If the profile is unimodal (no hollowness), the separation of the two Gaussians is in principle zero. On the contrary, for hollow profiles, the separation of the two Gaussians increases with and it is directly proportional to the degree of hollowness. The separation is quantified by the Bhattacharya indicator particularised for symmetric Gaussians. Then, by using a sigmoid function, the probability to have an ETH is computed (the coefficients d_1ETH and d_2ETH have been calculated using numerical cases reported in the reference paper [52]):

$$\begin{align}{p_{{\text{ETH}}}} & = {e^{{d_{1{\text{ETH}}}}\left( {{\text{ETH}} - {\text{ET}}{{\text{H}}_{{\text{threshold}}}}} \right) + {d_{2{\text{ETH}}}}}}\nonumber\\ & \quad / \left( {1 + {e^{{d_{1{\text{ETH}}}}\left( {{\text{ETH}} - {\text{ET}}{{\text{H}}_{{\text{threshold}}}}} \right) + {d_{2{\text{ETH}}}}}}} \right).\end{align} \tag{ 2.6 }$$

EC is detected using the cumulative distribution function (CDF) of the normalised electron temperature profile and by calculating an 'edge radius' ρ₉₈ as the radius containing the 98% of the CDF. Then, the lowest is the radius, the more severe is the EC. Again, by using a sigmoid function, the probability of the plasma being in an EC dangerous situation can be computed (the coefficients d_1EC and d_2EC have been calculated using numerical cases shown in [52]):

$$\begin{equation}{p_{{\text{EC}}}} = {e^{{d_{1{\text{EC}}}}\left( {{\text{EC}} - {\text{E}}{{\text{C}}_{{\text{threshold}}}}} \right) + {d_{2{\text{EC}}}}}}/\left( {1 + {e^{{d_{1{\text{EC}}}}\left( {{\text{EC}} - {\text{E}}{{\text{C}}_{{\text{threshold}}}}} \right) + {d_{2{\text{EC}}}}}}} \right).\end{equation} \tag{ 2.7 }$$

To interpret equations (2.6) and (2.7) two comments are in place. First, it should be noticed that both anomalies deform the temperature profile independently from the absolute value of this quantity. The relevant information resides in the shape of the profile. Indeed, for example, what really matters for the control system is how hollow a temperature profile is and this can be well quantified by a simple form factor. Consequently, the detection of electron temperature anomalies does not require any machine learning method (no training set, etc), because the critical values can be determined on the basis of previous experiments and do not need to be adjusted during the campaigns. Secondly, the probabilities expressed by equations (2.6) and (2.7) have to be interpreted differently from the one of the mode locking, described in the previous subsection. Indeed, p_ETH and p_EC quantify how anomalous the temperature profile is, not directly the plasma tendency to disrupt. Indeed, typically the mode locked also must become anomalous before the plasma configuration collapses. As will be discussed in section 2.4, therefore p_ETH and p_EC can be utilised to provide information about the propensity of the plasma to developing macroscopic magnetic instabilities, not directly about the plasma propensity to disrupt.

Contrary to mode locking and electron temperature anomalies, radiation anomalies do not lend themselves easily to a straightforward interpretation for various reasons. First, the phenomena influencing the radiation emission in a tokamak are extremely complex. Consequently neither physical nor empirical equations of radiative disruptive thresholds have been found, which can be applied in real time and easily transferred from one device to another. Moreover, contrary to temperature hollowness or EC, radiation cannot simply be eliminated. Therefore, an alternative physics-guided semi-supervised methodology has been developed to detect anomalous radiation patterns.

Another difficulty resides in the fact that the plasma radiation emission in tokamaks is usually estimated with the inversion of bolometric signals, utilising quite sophisticated tomographic algorithms. This approach allows obtaining high spatial resolution information and the advanced method based on the maximum likelihood can also associate confidence intervals to the results [53–56]. Unfortunately, the available tomographic algorithms are too time consuming and their deployment in real-time is, at the moment, unrealistic [57]. The workaround solution of calculating radiation peaking factors, as developed in [42], can be useful but, their interpretation can be delicate in some situations and they lack sufficient spatial resolution to identify all the potentially disruptive emission patterns detected on JET. A systematic analysis of the peaking factor limitations can be found in [58]. A possible alternative to estimating the local emissivity is the use of deep learning [59, 60]. However, being these methodologies based on supervised techniques, they may not return reliable results in anomalous situations. Moreover, transfer learning to new tokamaks should be developed to render this technology applicable from day 0 of a new device.

The alternative adopted in the present work consists of trading space resolution for time resolution. The emission is calculated for macro-pixels derived by macro-views, obtained by summing the corresponding lines of sight of the bolometers [58, 61]. Given the topology of JET bolometric cameras, it is possible to divide the projections into three vertical and three horizontal macro views, whose intersections identify the eight regions reported in figure 1. Obtaining the emission from these macro/pixels requires solving a system of six equations with eight degrees of freedom. By imposing that the emissivity must be positive, this system of equations can be solved by making recourse to non-negative least square minimisation methods [62]. Given the types of emissivities encountered on JET, the iterative routine implemented always converges in less than 50 μs, more than adequate since the cycle time of JET real time network is 2 ms. The absolute values of the radiation emission from the macro-pixels have been checked with the most accurate tomographic method available, the one based on the maximum likelihood [55]. Since it has been found that the average error is around 15%–20% in most cases, the proposed simplified inversion method has the spatial resolution, the time resolution and the accuracy to detect all the most important radiation anomalies potentially leading to disruptions on JET.

Zoom In Zoom Out Reset image size

Then, the eight macro-pixels are used to calculate the radiation in four macro-regions of interest:

• Core: this region is exactly the one shown in figure 1 left. This region allows detecting anomalous core radiation due to heavy impurity accumulation and is usually correlated with the onset of ETH.
• Divertor (Div): also this region is exactly the one shown in figure 1 left. This region allows evaluating excessive radiation in the region of the divertor.
• High-field: calculated as the sum of high-field left, high-field top and top. This region is used to detect MARFE and MARFE-like radiation patterns, phenomena that sometimes precede electron temperature edge reductions at the edge.
• Low-field: calculated as the sum of low-field bottom, low-field right and low-field top, this region is useful to monitor the radiation emitted by the frequently formed radiative blobs on the low-field side, which sometimes precede core radiation and hollowness or EC.

To predict the onset of radiation anomalies, the simple emission of each macro-pixel is not enough though. The really relevant quantity is the ratio of the radiated power divided by the local internal energy:

$$\begin{equation}{\Lambda _i} = \frac{{{P_{{\text{rad}}}}\left( {\rho ,\theta ,t} \right)}}{{{E_p}\left( {\rho ,\theta ,t} \right)}}\end{equation} \tag{ 2.8 }$$

where ρ and θ are the coordinates of the corresponding macro-pixel barycentre. The Λ_i quantities have the units of s⁻¹ and can therefore be interpreted as the inverse of cooling factors. How to calculate the denominator of equation (2.8), the average plasma energy of a macro-pixel, is explained in appendix A.

On a shot-to-shot basis, the training sets of 'safe' and 'anomalous' plasmas are updated and the cooling factors recalculated (see training set update logic section). However, a typical classification algorithm of classification cannot be used because:

• 1.  
  The anomalous training set would be very unbalanced, containing a lot (most) of time slices when the plasma is not anomalous.
• 2.  
  The safe training set may contain some time slices where the plasma is anomalous but the plasma does not disrupt for various reasons (for example because the ramp down has already started and the plasma dynamics is slow enough to allow the safe termination of the discharge).

This implies that the classification training set could be inaccurate and classification algorithms may fail. Therefore, the following simple but more robust approach has been developed. The probability density function of the various cooling factors is calculated in the two cases of safe and disruptive discharges. So, for each region, the pdfs of 'safe' and '(potentially) anomalous' Λ_i are derived from the available examples. What it is expected is that the pdfs would diverge at higher Λ_i values as shown pictorially in figure 2(a) (the pdf of anomalous data should contain more cases at larger Λ_i ). An example for the core radiation anomaly is shown in figure (b, left). From the pdfs the corresponding CDFs are calculated and the potential true positive and false alarm (in percentage) trends are obtained (figures (a) and (b), right). The true positive percentage (or sensitivity) trend is calculated as one minus the anomalous CDF, while false alarm percentage (or 1—specificity) is one minus the safe CDF.

Zoom In Zoom Out Reset image size

Once the false alarms and sensitivity trends are derived (figure 2(b) right), the Λ_threshold,i is evaluated by using the following rules:

• 1.  
  Search for the threshold that ensures a false alarm rate lower than 0.05% and a true positive larger than the false alarm rate.
• 2.  
  If no threshold satisfies the previous rule, the false alarm limit is increased to 0.1%.
• 3.  
  If no threshold satisfies the previous rules, the false alarm limit is increased to 1%.

Once the threshold has been determined on the basis of the CDFs, the probability to have a local radiative anomaly is calculated with a sigmoid function:

$$\begin{equation}{p_\Lambda } = {e^{{d_{1\Lambda }}\left( {{\Lambda _i} - {\Lambda _{{\text{threshold}},i}}} \right) + {d_{2\Lambda }}}}/\left( {1 + {e^{{d_{1\Lambda }}\left( {{\Lambda _i} - {\Lambda _{{\text{threshold}},i}}} \right) + {d_{2\Lambda }}}}} \right).\end{equation} \tag{ 2.9 }$$

Note that a Λ_threshold,i is found for each macro region (Core, HF, LF and Divertor) and an alarm is raised if one of the p_Λ exceeds 0.5.

As in the case of p_ETH and p_EC, also equation (2.8) is to be interpreted as the probability of the radiation being anomalous in a certain region, compared to the value assumed in safe discharges. Therefore, as discussed in the next subsection, p_Λ provides information about the probability of anomalies developing in the temperature profile not directly about the propensity of the plasma to disrupt.

By using the anomaly indicators described in the previous subsections, the algorithm is able to detect an anomaly and to classify it in one of the classes: mode locking, EC, hollowness, core radiation, HF radiation, LF radiation, divertor radiation. As mentioned, in the case of mode locking, the associated probability can be interpreted as the propensity of the plasma to disrupt. For the other classes, the probabilities quantify the how severe the corresponding anomaly is with respect to safe discharges. This is the basic information that the control system requires to take the proper countermeasures (mitigation, prevention, and avoidance). However, in order to accurately control the plasma and optimise feedback schemes, it would be very useful to have an idea of how much time is available to undertake a certain action. For this reason, three specific predictors have been developed that predict the time to hollowness, time to EC and time to mode locking (t_{r ,LM}, t_{r ,ETH} and t_{r ,EC}).

The task is a typical supervised regression problem and for each predictor a feed-forward neural network is used. The architecture used is a three hidden layers neural network with 20, 10 and 5 neurons respectively. The output, i.e. the time to the anomaly, are predicted in logarithm scale in the range from 1 ms to 2 s (2 s is considered the time horizon). The loss function is a weighted MSE aiming at weighting more shorter times, to compensate for the fact that in shorter intervals there are less data (for example, from 1 ms to 10 ms there are 10 points per pulse, while from 100 ms to 1 s there are 900 points):

$$\begin{equation}{\text{los}}{{\text{s}}_{W,{\text{MSE}},j}} = \frac{{\sum\limits{W_i}{{\left( {{{\text{log} }_{10}}{t_{j,{ }i}} - {{\text{log} }_{10}}{y_{j,i}}} \right)}^2}}}{{\sum\limits{W_i}}}\end{equation} \tag{ 2.10 }$$

where:

$$\begin{equation}{W_i} = \frac{1}{{{y_{j,i}}}}.\end{equation} \tag{ 2.11 }$$

While ${t_{j,\,i}}$ and ${y_{j,i}}$ are the predicted and the effective time to anomaly j of the ith observation. These times are calculated continuously for the entire discharge providing each network with the following signals as inputs: the inputs are the plasma current (I_p ), the magnetic toroidal field (B_t ), the dimensionless locked mode (LM_N ), the internal inductance (l_i ), the electron temperature and density in the core, middle and edge regions, the electron temperature profile anomaly indicators (ETH and EC), the Λ_i in the four regions (Core, HF, LF and Div) and the Λ_heating.

The electron temperature and density in the core, middle and edge regions are the average values calculated considering the radii proposed in [12] (core for R₀ < R < 3.4, middle for 3.4 < R < 3.6, and edge for 3.6 < R < 3.8), while Λ_heating is the total input power normalised to the plasma energy (note that analogously to radiative Λ, this value is the inverse of a time that represents, in a certain sense, a characteristic heating time).

The mode locking predictor is based on the locked mode signal and it typically needs to intervene at the end of the pre-disruptive phase, to trigger the measures of last resort, the mitigation actions. Consequently, the training set must be update only in clear cases:

• 1.  
  If the previous pulse is disruptive and no remedial actions have been undertaken (no mitigation, prevention or avoidance), the last 100 ms are inserted in the training set and classified as anomalous (class = 1).
• 2.  
  If the previous pulse is safe and no alarms have been raised, the pulse is stable and 100 random points are sampled sequentially, added to the training set sequentially (respecting the time ordering) and classified as safe (class = 0).

In the case of radiation anomaly, the training set update is less restrictive, and it follows the rules:

• 1.  
  If a mode locking or an electron temperature anomaly is detected in the last pulse, the previous 2000 points are added to the training set and labelled as anomalous (class = 1). It has indeed proved more effective to train the predictor with the information about the anomalies caused by abnormal radiation than with the disruptive time slices.
• 2.  
  If the last pulse triggered no alarms and it did not disrupt (safe pulse), all the points of the pulse are added in the training set and labelled as not anomalous (class = 0).

The training set, to estimate the time remaining before an anomaly, is updated following the criteria:

For time to mode locking:

• 1.  
  If the mode locking detection algorithm predicts a locked mode anomaly, the 2000 points before the alarm are taken and the targets are given as a function of the actual difference between the alarm time and the actual beginning of the current quench.
• 2.  
  If no mode locking alarms are triggered and the pulse is not disruptive, 200 points are added to the training set with the target equal to the time horizon (2 s).

For Time to EC

• 1.  
  If the EC detection algorithm predicts an EC anomaly, the 2000 points before the alarm are taken and the targets are given as a function of the actual difference between the alarm time and the actual beginning of the current quench.
• 2.  
  If no EC alarms are triggered and the pulse is not disruptive, 200 points are added to the training set with the target equal to the time horizon (2 s).

For time to hollowness:

• 1.  
  If the hollowness detection algorithm predicts a hollowness anomaly, the 2000 points before the alarm are taken and the targets are given as a function of the actual difference between the alarm time and the actual beginning of the current quench.
• 2.  
  If no hollowness alarms are triggered and the pulse is not disruptive, 200 points are added to the training set with the target equal to the time horizon (2 s).

The dataset is randomly divided into training set (70%) and test set (30%).

As mentioned all these updates are implemented after each discharge. The computational times, being of the order of minutes, are fully compatible with the time between discharges.

Note that disruptive pulses are fewer than safe pulses and this may lead to unbalanced dataset. To alleviate this problem, different methodologies have been used. For mode locking classification, the cross-entropy loss function has been used. For time to anomaly prediction, on a safe pulse we sample fewer points than disruptive points and we also weighted the one closer to the disruption (equation (2.10)).

The database consists of 1683 JET pulses, spanning from the high-power DD campaign to the full tritium and DT campaign (C38, C39, C40, C41). The DB includes 1141 not disruptive pulses and 542 disruptive. The flat-top plasma current ranges from 1 MA to 3.5 MA, while the toroidal magnetic field goes from 1.7 T to 3.9 T. The maximum input power is 37.8 MW, where 31.6 MW were provided by the NBI and 5.1 MW by the ICRH. See table 1 for more details.

All disruptive shots have been included, except intentional disruptions and a few cases, in which essential data are missing. Note that, even if some of these pulses are in an operational regime where disruption would not involve danger for the tokamak, they have been included to test the capability of the predictor to perform well in a large range of plasma parameters. A graphical overview of the entire database is provided in figure 3. The bins cover 0.5 MA centred on the values indicated on the x axis (for example the bin centred at 2.5 MA includes discharges in the interval 2.25–2.75 MA). The absolute frequencies (figure 3(a)) are calculated dividing the number of safe (disruptive) cases in the bin e by the total number of safe (disruptive) pulses in the entire database. The relative frequencies (figure 3(b)) are calculated as the number of safe (disruptive) discharges in each bin divided by the total number of discharges in the same bin.

Zoom In Zoom Out Reset image size

In order to evaluate the performances of the proposed predictors, they have been run in conditions exactly the same as the ones they would have faced in real life. To simulate an online situation, it is important to define the performance metrics and the alarm priorities.

The alarm priorities are in order of decreasing priority: mitigation, followed by prevention and then avoidance. It means that, for example, if both mitigation and prevention alarms are raised at the same time, it is assumed that the pulse must be mitigated because there is no time left for anything else. Therefore, the performance of the predictors will be evaluated according to the following classification, in which mitigation time threshold indicates the minimum time needed for implementing mitigation actions, prevention time threshold indicates the minimum time needed for implementing prevention actions and avoidance time threshold the minimum time needed for undertaking avoidance strategies.

In the case of not disruptive pulses:

• 1.  
  If no alarms are triggered the pulse is classified as safe no false alarms.
• 2.  
  If the predictor launches only avoidance alarms, the pulse is classified as safe avoidance.
• 3.  
  If the predictor triggers only a prevention alarm, the pulse is classified as safe false prevention.
• 4.  
  If the predictor triggers only a mitigation alarm, the pulse is classified as safe false mitigation.
• 5.  
  If the predictor triggers both prevention and mitigation alarms, there are two possibilities. If the prevention alarm time is larger than the prevention time threshold, the pulse is classified as safe false prevention, otherwise it is classified as safe false mitigation.
• 6.  
  If the predictor triggers avoidance alarms and prevention and/or mitigation alarms, there are three possibilities. If the avoidance alarm time is larger than avoidance time threshold, the pulse is classified as safe avoidance, otherwise it is classified as safe false mitigation or safe false prevention according to the previous rules.

In the case of disruptive pulses:

• 1.  
  If no alarms are triggered the pulse is classified as missed.
• 2.  
  If only a mitigation alarm has been triggered and the mitigation warning time is smaller than the minimum time for mitigation, the pulse is classified as tardy mitigation.
• 3.  
  If only a mitigation alarm is triggered and the mitigation warning time is larger than the mitigation time threshold, the pulse is classified as good mitigation.
• 4.  
  If only a prevention alarm is triggered and the prevention warning time is smaller than the prevention time threshold, the pulse is classified as Tardy prevention.
• 5.  
  If only a prevention alarm is triggered and the prevention warning time is larger than the prevention time threshold, the pulse is classified as good prevention.
• 6.  
  If both mitigation and prevention alarms are triggered, if the prevention warning time is larger than the prevention time threshold, the pulse is classified as good prevention. Otherwise, if the mitigation warning time is larger than the mitigation time threshold, the pulse is classified as good mitigation, else it is a tardy mitigation.
• 7.  
  If an avoidance alarm is raised and the avoidance warning time is larger than the avoidance time threshold, the pulse is classified as good avoidance, otherwise the previous rules apply.

Figure 4 summarises the metrics and the main parameters implemented to define the metrics. Note that two avoidance levels have been considered in the scheme. They are just triggered as a function of the alarm. If a radiation anomaly is detected without an ETH alarm, the plasma should be in an operational region more stable than the operational region where hollowness is detected. Therefore, different avoidance schemes may be used. Consequently, with avoidance lv2 it is intended a radiation anomaly with no electron temperature anomalies, while anomaly lv1 is triggered when hollowness is detected (irrespective of the presence of radiation anomalies). All the details and a table about the evaluation criteria are provided in appendix B.

Zoom In Zoom Out Reset image size

In this section, an example of how the predictor evolves during the campaign, for a fixed value of the hyper-parameters, is provided. In section 5.4, how performances change as a function of the required warning times for prevention and avoidance is analysed. The results in this section have been obtained assuming that the minimum warning time to undertake avoidance and prevention actions are 100 ms. The minimum time required to undertake successful mitigation actions is assumed to be 10 ms, to reflect the capabilities of JET tools.

At the beginning of the campaign, it is assumed that no experimental data is available. Thus, the mode locking predictor is trained with the physics loss function (the operational space after the first training is shown in figure 5—top left). The hollowness and EC indicators are determined by the thresholds that are self-determined (no training is needed). The radiation anomalies cannot be detected because it is considered too difficult to extrapolate the impurity and radiation behaviour to new, possibly much larger devices (however, if future studies will find some reliable scaling laws for radiation anomalies, they may be implemented following a logic similar to the one used in this work for mode locking). It is assumed that also for the time-to-anomaly predictors there is not enough prior knowledge to train them without specific data of the new device.

Zoom In Zoom Out Reset image size

In the analysed database, the evolution of the experiments is such that only after 30 discharges it is meaningful to perform the first training with the new data. Up to that point (before the first training), 14 disruptive pulses have been observed and 6 were mitigated, 3 prevented, 5 avoided (ETH) without neither tardy alarms nor missed disruptions. Regarding the 16 safe pulses, 3 were mitigated and 3 were avoided. However, by analysing the 6 pulses where 'false alarms' have been triggered, it can be observed that the anomalies are actually present, even if they do not lead to a major disruption. A typical example is pulse 94 162, the 6th pulse in the database. In this case, a 'false avoidance' due to ETH and a following false mitigation alarms are raised. However, both anomalies actually occur during the discharge. Figure 6 shows the main signals and anomaly indicators for this pulse. Both core radiation anomaly and hollowness indicators increase from t = 13 s, with a quite deep hollow profile from 13.22 s to 14.042 s, when the n = 1 mode locking amplitude increases (triggering a mitigation alarm). However, in this case unusually the mode locking does not lead to a disruption. Moreover in this shot mitigation actions would not have had any major detrimental effect because the plasma current was already very low at the moment of the alarm.

Zoom In Zoom Out Reset image size

After 30 pulses, all the predictors (mode locking, temperature profile anomalies, radiation anomaly and the time-to-anomaly) have been retrained by combining the physics loss function and the data driven loss function, derived from the training set automatically extracted from the first 30 pulses of the campaign. From this discharge and on the shot to shot basis, the entire algorithm detects and classifies the various anomalies and is automatically updated, by taking decisions according to the rules described in the previous section. The evolution of the operational space of the LM predictor and the radiation thresholds for the core region are shown in figures 5 and 7 respectively. The trends of this macropixel are representative of those of the other 3 regions. Figure 7 must be interpreted as described in figure 2 and the relative text.

Zoom In Zoom Out Reset image size

The overall performances as a function of the progressive pulse number are shown in figure 8. At first, it can be observed that the sensitivity of the predictor is very high. At the end of the campaign, only three disruptive pulses were not correctly avoided, prevented, or mitigated (they were all tardy mitigations). Of the disruptive pulses, 190 have been mitigated, 81 prevented and 268 avoided. Of the three pulses for which the mitigation alarm was launched less than 10 ms before the beginning of the current quench, in two cases the warning time was 9 ms and in one case 2 ms. Given the length of the current quench (at JET, most of current quenches last for at least 30 ms, and usually they are between 50 and 150 ms), mitigation would have been potentially effective also in these three discharges [63, 64]. Moreover, current quenches should be longer for large size future tokamaks (in ITER it is expected that they will be at least one order of magnitude longer) [63, 64].

Zoom In Zoom Out Reset image size

For what concerns the specificity, it is important to discern between 'false anomaly detection' and 'wrong action'. 'False anomaly detection' occurs when an indicator detects an anomaly that is not actually present in the discharge. On the contrary, 'wrong action' is when a control/mitigation action is triggered but the plasma does not actually need it. For example, if the LM indicator detects an alarm and a locked mode is actually present, the outcome is a 'true anomaly detection'. However, if a mitigation action is triggered but the plasma would have survived without this action, this is considered a 'wrong action'.

With regard to the anomaly detection specificity, this value is extremely high. For example, every time the hollowness indicator detects a hollow profile the hollowness is actually present (this can be easily confirmed by visual inspection, and a detailed analysis is performed in [52]). Therefore, the specificity of the indicators is of the order of 100% for mode locking, EC and hollowness. Note that these statistics cannot be provided for radiation anomalies, since today there is not a physics-based definition of radiation anomaly. However, it is possible to estimate the reliability of radiation indicators through a comparison with the other anomaly indicators. For example, in 86% of safe pulses where core radiation anomaly has been detected also electron temperature anomaly and/or mode locking is observed. In the case of high-field and low-field radiation anomaly, this percentage drop to 62% and 71% respectively. However, this estimation is conservative since it is expected that radiation is only one of the drivers (other are input power and transport) of electron temperature anomalies.

Different is the case of 'wrong actions'. In this work, it is assumed that all actions by the control system would have been undertaken at a probability threshold of 50% for each anomaly. For example, if the mode locking amplitude probability exceeds 0.5, mitigation actions would have been triggered. However, it is clear that in a real-case, when different control strategies can be taken, performances can be improved. For example, there are several studies regarding the possibility of de-locking magnetic islands, recovering from this instability, without the need to mitigate the plasma [20, 65, 66]. In such a case, the pulse could have been recovered by the appropriate control action, without the need of mitigating. This would reduce significantly the amount of false mitigation and prevention actions. Moreover, it has to be considered that triggers of avoidance do not lead to 'hard' actions in the plasma but only 'soft' ones, i.e. control actions that should allow to run the pulse in a safer situation. Most of them would improve the plasma performances even in discharges not disrupting. In other words, a 'false avoidance' trigger is much less harmful than a 'false mitigation' trigger and must be interpreted differently.

Another important consideration is about the feedback actions actually implemented by JET real time control system. After years of work, on JET several controllers are operational to mitigate or early terminate the discharge [67–70]. This implies that the predictors developed in the present work may detect anomalies that did not lead to a disruption because of JET control system intervening with the appropriate remedial actions. In this part of the work, the statistics have been calculated assuming that safe pulses are stable pulses and that disruptive ones are not. Therefore, the statistical results previously reported are to be considered conservative, since many safe pulses, in which anomalies are detected, are classified as 'wrong actions' even if the detected anomaly is absolutely real (but the plasma does not disrupt for various reasons including JET real time system interventions). To investigate this aspect in more details, section 5.3 reports the statistics for a subset of pulses, which have been analysed manually, to assess in detail the performances of the developed predictors and how they fare compared with JET control/prevention/mitigation actions.

Figure 9 (left) shows the cumulative fraction of warning times for each indicator, while figure 9 (right) show the cumulative fraction of three different predictors: the first one uses only the mode locking indicators, the second implements also the electron temperature anomalies, while the third combines all the indicators (magnetic, kinetic profiles, radiation patterns). The mode locking-based predictor LM returns the typical warning time curve observed in other works, with a median warning time around 200 ms. The use of hollowness and EC indicators allow to increase the median warning time from ∼200 ms to ∼500 ms, where the strongest contribution in the large warning time is given by the hollowness (EC are quite fast phenomena, as also discussed by [12], with a warning time between EC and mode locking of the order of 50–100 ms). The implementation of radiative indicators allows to increase the median warning time to 900 ms, with at least 90% of pulse predicted at least 200 ms before the beginning of the current quench.

Zoom In Zoom Out Reset image size

Concerning the time to anomaly predictions, figure 10 (left) shows the median error committed by the predictors as a function to the actual time to anomaly. The errors decrease linearly in logarithmic scale (as expected, being the predictor trained in log-scale), showing the capability of the neural networks to estimate correctly the order of magnitude of the time remaining before any of the three anomalies is observed. It is worth noting that the time-to-mode-locking is the most accurate prediction since mode locking is normally preceded by the other anomalies (radiation and temperature) with a quite large margin, allowing a reliable evaluation of the plasma stability. On the contrary, prediction of time-to-hollowness is the most difficult because this modification of the temperature profile can appear very close to core radiation anomaly and it is not preceded by other (real-time) clear features on the signals. For this reason, the 'time horizon' of the time-to-hollowness prediction is very small and therefore the error of the prediction quite high. The EC error is in between, since it is usually preceded by a clear anomaly radiation on the high-field side (MARFE-like radiation) and/or low field side radiation.

Zoom In Zoom Out Reset image size

Figure 10 (right) shows the coefficient of determination R² calculated for each predicted pulse as a function of the pulse numbers [71]. Here, it can be seen that, on average, the prediction performances increase with the pulse number, and this is to be expected since more pulses are used to train the predictors and therefore they gain in accuracy. However, it can be observed that performances are quite good even just after the few pulses used for the training. It has to be highlighted that the time to anomaly prediction is based on a fully data-driven methodology, since theoretical/empirical equations deployable in real-time have not been developed yet. However, future works on this subject may achieve a better prediction of the time to anomaly from day-0, by exploiting physical or empirical models of these quantities.

In this section, some examples are reported to illustrate the functioning of the predictors, the various cases leading to recovering actions, and the amount of information that they are able to provide for the control of the plasma.

The first example is the pulse 96 486, a baseline scenario with a 3.5 MA flat top plasma current and ∼35 MW input power (30 MW from NBI, 4.2 MW from ICRH). The main signals are shown in figure 11. From t = 13 s, a second before the ramp-down phase, core radiation starts to slightly increase. At the beginning of the ramp down, due to the input power decrease, the ELMy mode terminates (figure 11, top-left) with an increase of impurity accumulation in the core. At t = 14.179 s, anomalous core radiation is detected, and at t = 14.246 s also the ETH alarm is triggered. The core radiation anomaly grows exponentially and at t = 15.47 the mode locking alarm is triggered. In this pulse, the only disruption prevention/mitigation action undertaken by JET control system was the activation of the disruption mitigation valve at t = 15.492. In figure 11 (bottom-left), the various alarms raised (following the logic described in the previous sections) are shown. In brief, the predictor raised an avoidance alarm 1.292 s in advance with respect to the mode locking. Then, the mitigation alarm is triggered 31 ms before the beginning of the current quench. It has to be mentioned that this is a pulse mitigated by JET control system; consequently, the actual warning time provide by the developed predictor may be larger than the one just quoted.

Zoom In Zoom Out Reset image size

The core-electron temperature anomaly behaviour can be observed also in figure 12 (right). It can be observed that at first the core radiation anomaly slightly increases. This leads to an increase of the ETH that amplifies the core radiation for two reasons. First the impurities have a radiation function which increases when the temperature goes down. Secondly to maintain pressure equilibrium, during the onset of the ETH the core density peaking is observed to increase, favouring the accumulation of impurities in the core. Both effects contribute to establishing a clear positive feedback loop that makes the thermal-radiative instability grow exponentially.

Zoom In Zoom Out Reset image size

Figure 13 shows the time to anomaly predictions. When the plasma in a stable operational regime, far away from the anomaly onset, the time to anomaly prediction return a very high value, comparable with the time horizon (between 1 s and 2 s). When core radiation starts increasing (t > 13 s), a drop of the time to hollowness ETH signals is observed and when an anomalous core radiation is detected (t = 14.179), the time to ETH is around 35 ms. After 67 ms, the ETH anomaly is detected. The good agreement between predicted and actual time to mode locking is shown in figure 13 (right). As anticipated in the previous section, the prediction of the time to mode locking anomaly is much more reliable than the time to electron temperature anomaly since the mode locking is preceded by clear and measurable effects (radiation, temperature).

Zoom In Zoom Out Reset image size

The second example is pulse 94 650. This is another baseline scenario at 3 MA plasma current and 25 MW input power (22.3 MW from NBI and 2.4 MW from ICRH).

Figure 14 shows the most relevant signals for pulse 94 650. At first, two spikes are observed in the edge regions (HF and LF) at t = 6.8 s and t = 13.45 s, probably caused by some impurity influxes. From t = 15.5, an increase of edge radiation is observed and at t = 15.8 s an intense radiation activity is observed on the high field side (MARFE-like radiation, also confirmed by a visible camera). After 175 ms, the edge radiation leads to an EC of the electron temperature profile (detected a t = 16.142 s) and a consequent mode locking at t = 16.515 s. The pulse disrupts at t = 16.611 s, with a plasma current equal to 1.73 MA. Therefore, in this case, the avoidance warning time (with respect to the prevention alarm) lasts for 175 ms, followed by a prevention warning time (with respect to the mitigation alarm) of 373 ms, and a mitigation warning time of 96 ms.

Zoom In Zoom Out Reset image size

Figure 15 reports the signal in the space HF radiation indicator—EC indicator. Also in this case, it is possible to observe that the radiation and temperature anomalies increase together, with an exponential growth of the EC signals, since the local decrease of the electron temperature causes a local increase of the electron density (pressure balance), with a consequent amplification of the local radiation (confirming the positive feedback loop).

Zoom In Zoom Out Reset image size

In figure 16, the time to anomaly predictions for pulse 94 650 are shown. In this case, both the time to mode locking and time to EC predictions follow quite well the actual evolution of the respective anomalies, when the pulse is in the unstable phase.

Zoom In Zoom Out Reset image size

In this subsection, a subset of the entire database is considered, to compare the results of the predictors with the already implemented real-time JET controllers and prevention/mitigation systems. The pulses included in this analysis belong to the DD campaign for DT preparation (the so-called C38) and they were performed in the framework of to the two main JET tasks of developing the 'ITER baseline scenario' and the 'hybrid scenario'. In total, this data set comprises 75 disruptive and 211 safe discharges. These are discharges in which the interventions of JET control system was very clear and properly documented.

On JET, different indicators are already implemented to detect plasma anomalies during a discharge and trigger some actions to mitigate or prevent disruptions.

Regarding the mitigation actions, on JET there are disruption mitigation valves (DMV) that inject gas inside the plasma. This gas causes the plasma to irradiate, loosing energy, and so reducing the loads on the plasma facing components, since most of the plasma energy is lost isotropically by the radiation channel [72, 73].

A second relevant action, that can be considered a prevention action, is the so called 'jump to termination' or JTT; it is used to terminate the plasma discharge early in a controlled way when some anomalous radiation events occur (e.g. tungsten impurity accumulation).

Detailed descriptions of the JET real-time implemented algorithms and control actions can be found in the literature [67–70].

Table 2 shows the comparison between the present predictor and the JET termination actions in the dataset considered and the main findings are reported below.

All the pulses where the JET control fired the disruption mitigation valve are detected significantly earlier by the predictors. In 28 pulses, the avoidance alarms are triggered with sufficient warning time to make avoidance strategies feasible, while 5 of them would require prevention and 3 mitigation. The median anticipation between the avoidance alarm and the DMV is 689 ms, the average is 1.115 s, and the minimum value is 135 ms. For prevention, the median is 118 ms, the average 135 ms, and the minimum warning time is 80 ms. Finally, for mitigation there is a median value of 165 ms, an average of 481 ms and a minimum value of 25 ms.

In 19 disruptive pulses a JTT procedure was undertaken and in 18 of them successful avoidance alarms would have been triggered by our predictors, whereas in one only the mitigation alarm would have been triggered. The median avoidance warning time (with respect to JTT) is 120 ms, with an average of 350 ms. The minimum value is negative (−0.45 ms), since two pulses triggered an avoidance alarm slightly delayed with respect to JTT.

In 20 disruptive discharges, the JET system did not trigger neither a JTT nor DMV action. However, all of these disruptions occur below the JET plasma current threshold (2 MA) and therefore actions might not have been triggered for this reason. In any case, of these 5 pulses would have satisfied the criterion for mitigated, 4 for prevention, and 11 for avoidance.

Out of all 33 safe pulses where JET triggered a JTT action, the mitigation predictor triggers seven times: it has been checked that mode locking was actually present in the ramp down, below the plasma current threshold for DMV triggering. In ten discharges the avoidance predictor properly launches and alarm. In all safe pulses where JET did not take any actions, also the predictors trigger no alarms. This confirms that the algorithm is very specific, and that the false alarms are mostly triggered in situations where the plasma is actually unstable even if it does not disrupt.

The results reported in the previous subsection have been obtained assuming that the warning times required for successful avoidance and prevention are both 100 ms (and 10 ms for mitigation). However, since there is still some research on control schemes for avoidance and prevention, it is clear that the minimum warning times may be different in future experiments. It has also to be considered that different warnings time may be required as a function of the anomaly typology and intensity. For these reasons, a parametric study of how the performances vary as a function of the required warning times has been performed.

Figure 17 shows how the performances in disruptive pulses change as a function of the required avoidance warning times (left) and prevention warning times (right). As expected, avoided disruptive pulse percentages decrease with increasing warning times. However, this decrease is limited and does not affect the missed/tardy pulses, but simply implies that the numbers of mitigated/prevented pulses increase. Similar results are observed for the increase of required prevention warning times. If the required warning increases, less pulses will be prevented, but mitigation actions will suffice to avoid unmitigated disruptions. It has to be note that in this case, since prevention priority is higher than avoidance (section 4.2—figure 4), the pulses, in which avoidance alarms are triggered, remain the same.

Zoom In Zoom Out Reset image size

The results shown in figure 17 clearly demonstrate that the adaptive training logic has a good reliability and that all the predictors work in such a way as to ensure that disruptive pulses are at least mitigated, favouring actions that aim at triggering prevention/mitigation actions.