Detection of changes in the dynamics of thermonuclear plasmas to improve the prediction of disruptions

In particular circumstances, nonlinear systems can collapse suddenly and abruptly. Anomalous detection is therefore an important task. Unfortunately, many phenomena occurring in complex systems out of equilibrium, such as disruptions in tokamak thermonuclear plasmas, cannot be modelled from first principles in real-time compatible form and therefore data-driven, machine learning techniques are often deployed. A typical issue, for training these tools, is the choice of the most adequate examples. Determining the intervals, in which the anomalous behaviours manifest themselves, is consequently a challenging but essential objective. In this paper, a series of methods are deployed to determine when the plasma dynamics of the tokamak configuration varies, indicating the onset of drifts towards a form of collapse called disruption. The techniques rely on changes in various quantities derived from the time series of the main signals: from the embedding dimensions to the properties of recurrence plots and to indicators of transition to chaotic dynamics. The methods, being mathematically completely independent, provide quite robust indications about the intervals, in which the various signals manifest a pre-disruptive behaviour. Consequently, the signal samples, to be used for supervised machine learning predictors, can be defined precisely, on the basis of the plasma dynamics. This information can improve significantly not only the performance of machine learning classifiers but also the physical understanding of the phenomenon.


Introduction
Many complex systems, either human-made or natural, sometimes can collapse abruptly and in quite dramatic fashion. The spectrum of these failures includes the most diverse examples such as extreme weather events, crashing markets, bankruptcy of companies, technological disasters, epidemics, revolutionary political changes and many others. Such failures constitute a challenge for the relatively new field of complexity theory, whose objective is precisely studying the behaviour of these systems. The post-mortem study of collapses is very important. Some of them have been understood quite well and can be reliably attributed to some specific properties such as group behaviour, evolutionary processes, nonlinear dynamics and network effects, just to name a few [1]. The appropriate analysis of these aspects is crucial since ill-founded statistical assumptions about the occurrence of extreme events have led to catastrophic results or jeopardized their prevention. Consequently, the goal of real-time prediction remains probably the most important and the most challenging.
Nuclear fusion is the process, by which light nuclei coalesce into heavier ones. The fusion of the hydrogen isotopes and He is the most exo-energetic reaction in the known universe. Nuclear fusion is indeed the process of powering the stars and the sun. The potential advantages of fusion as a source of energy are very appealing [2]. The fuel is almost inexhaustible, and the reactors would not produce radioactive waste and could not explode. Unfortunately, to come close enough for the strong forces to take over and fusion to occur with sufficiently high probability, the nuclei have to overcome the Coulomb barriers. This requires the fuel to reach very high temperatures of hundreds of millions of degrees (higher than in the core of the sun). The most credible solution for achieving controlled thermonuclear fusion in the laboratory relies on confining high-temperature plasmas with the help of strong magnetic fields. In this area of magnetic confinement fusion (MCNF), the most advanced configuration is the tokamak [3]. In tokamak devices, the plasma is confined in toroidal vacuum chambers. High currents, of the order of several millions of Amperes, are generated inductively; the experiments therefore consist of discharges with a repetition rate typically of about one hour in the present largest devices.
Tokamak plasmas are very complex systems, from both a technological and physical point of view. They are kept well out of equilibrium by continuous injection of material and megawatts of power. They are very delicate and poorly accessible for measurement. They are also strongly nonlinear and plagued by a series of instabilities of various consequences.
One of the major issues on the route of a commercial tokamak reactor is the occurrence of macroscopic instabilities called disruptions, because they cause a complete loss of confinement and the abrupt extinction of the discharge. The consequences of disruptions are high thermal loads on the plasma-facing components, strong forces on the electromagnetic structures and the generation of beams of runaway electrons in the MeV range. Unfortunately, all these effects scale badly with the dimensions of the devices, and therefore, in the demonstrative reactor, even a single disruption could force the shutting down of the machine. Indeed, already on ITER (International Thermonuclear Experimental Reactor, Cadarache, France, a 500 MW tokamak fusion device designed to prove the feasibility of fusion as a large-scale source of energy), full current disruptions could cause forces on the structure comparable to the weight of an Airbus. The thermal loads could even exceed by an order of magnitude the most resistant materials available nowadays.
Disruptions are very complex and nonlinear phenomena. Therefore, to this day, no first principle model has been devised to reliably predict their occurrence in real-time compatible form. An example of the signal evolution leading to a disruption is reported in Fig. 1. From a simple inspection of the Fig. 1 Time evolution of most relevant global signals in the interval preceding the beginning of a disruption. The time instance when the disruption occurs is marked by the vertical red bar experimental signals, it is clear that significant dynamical changes occur in the plasma dynamics preceding a disruption. One strategy to prepare for ITER relies on the development of data-driven predictors based on machine learning. Some of them have been very performing, reaching success rates of the order of 98% and false alarm rates of a few per cent [4][5][6]. Three were indeed installed in the JET tokamak real-time network, and one, called APODIS, was also used in the feedback loop [7]. For a review on machine learning methods for disruption prediction, the reader is referred to [8]. The main drawback of this first generation of predictors resided in their requirements of examples for the training. Hundreds of disruptive discharges were necessary, and this is a luxury that the new net generation of devices will not be able to afford. Therefore, optimizing the training process, by improving the choice of the examples, is of paramount importance.

The database and the main signals for disruption prediction
The examples analysed in the rest of the paper are discharges of the Joint European Torus (JET), the largest tokamak ever operated in the world. Located near Culham in Oxfordshire, JET was the device of the European Community until the end of 2021. For about 10 years, it has been operated with the new ITER-like wall (ILW), made of the same plasma-facing components as those of the next-generation device ITER: W in the divertor and Be in the main chamber. JET is also the only operational device in the world, which can use the fuel of the reactor, a 50/50 mixture of deuterium and tritium. As can be deduced by a simple inspection of the plots of Fig. 1, before the actual onset of the disruption, many important quantities exhibit signs of anomalies, called precursors. The disruption proper starts with the thermal quench, during which the plasma loses almost all its kinetic energy on time scales of the order of 1 ms. The thermal quench is then followed by the current quench, when the plasma is extinguished on time scales of hundreds of milliseconds on JET with the ILW.
The time series most important for mitigation are known to be: The first two time series, as the present study shows, present morphological properties, which can provide rich information about the dynamical changes and they are available routinely. The calculations have been performed for pulses from various experiments, as shown in Appendix 1, in an attempt to remove bias. It should be noted that the disruption times listed in Appendix 1 are provided by the JET operator. The disruption time is conventionally determined as the beginning of the current quench. This is an inherent property of tokamak plasmas accepted by all the specialists.
3 Overview of the tools for the detection of changes in the plasma dynamics Critical transitions in complex systems often occur unexpectedly. Consequently, they are usually difficult to predict and to manage. Therefore, the identification of early warning signals to predict significant changes in dynamical systems, bifurcations and tipping points has been a growing research subject in the last years, covering a wide range of topics such as climate, paleoclimate and ecology studies [9][10][11], psychotherapeutic phase transitions [12], financial markets [11,13,14] and behavioural science [15]. However, the applicability of the methods for predicting critical transitions is frequently restricted to specific applications. Consequently, the limitations of the different methods to identify approaching critical transitions should not be neglected [16]. The choice of the early warning signals and the associated detection techniques should be carefully adjusted to specific applications. In the case of fusion plasma disruption predictions, the data analysis tools, adapted and deployed to investigate the evolution of the disruption precursors, are overviewed in this section. Section 3.1 describes an established method of calculating the embedding dimension. The indicators derived from the recurrence plots are covered in Sect. 3.2. The technique aimed at detecting the onset of chaos is the subject of Sect. 3.3.

Embedding dimension
The extraction of information about the features of the phase space from time series is allowed by the Takens' theorem [17]. The unfolding of a system attractor from observations of a single dynamical variable is retrieved by the time delay method: where m is the embedding dimension, s is the time delay and e j ! are the unit vectors spanning an orthogonal coordinate system. Takens' theorem ensures that the original and the reconstructed attractor have the same properties, provided that m [ D 2 þ 1, where D 2 is the correlation dimension of the attractor. The embedding delay can be computed by mutual information or auto-correlation [18], while the embedding dimension can be efficiently computed by using the method of the ''False Neighbours'' [19]. This method searches for the optimal embedding dimension by iteratively increasing the value of m. If, at a certain step of this iterative process, m is too small to unfold the attractor, then not all the points that lie close to each other are real neighbours. Part of them look like neighbours only because the structure of the attractor has been projected down to an excessively small dimensional space. When the percentage of false nearest neighbours is negligible, the correct embedding dimension is reached. The use of the embedding dimension changes, as indicators of approaching transitions, has been proposed in [15], and it was interpreted as a quantifier of the change in the complexity of the context influencing the behaviour of the dynamical system.

Recurrence plots
In the last decade, the understanding of the underlying dynamics and the detection of regime changes in the evolution of dynamical systems have increasingly relied on recurrence-based approaches, as recurrences encapsulate the properties of dynamical systems. The concept of recurrences has its origin in the famous Poincaré recurrence theorem [20], which states that, under certain conditions, any dynamical system returns arbitrarily close to a neighbourhood of each of its former points at least once. A significant advance in the recurrence-based analysis came with the proposal of the recurrence plots (RP), in the famous paper of Eckman et al. [21]. RP encodes the pairwise recurrences of time series values and creates a visual representation of the recurrence dynamics. The reader is referred to [22] and [23] for reviews of RPs.
Considering two state vectors on a trajectory of a system in phase space, they are considered to be recurrent if the second vector falls into the neighbourhood (defined as a -radius sphere) of the first vector. RPs depict the collection of pairs of times i; j ð Þ at which recurrences occur, which means that the trajectory returns sufficiently close to the same place. RPs are based on the following matrix representation: where H Á ð Þ is the Heaviside function, Á k k is the Euclidian norm, d ij is the Kronecker delta symbol, N is the trajectory length and is usually an arbitrarily selected small number. The choice of the parameter may have a strong influence on the structure of the RP. Various methods have been proposed for the selection of this parameter [24][25][26][27]. In the present approach, has been set at 5% of the phase space dimension.
The graphical representation of an RP is an N Â N grid of points, which are encoded as black for 1 and white for 0. A black point in the RP means that the system returns to an e-neighbourhood of the corresponding point in phase space.
Besides the visualization of recurrences, RPs allow also deriving very important properties of a system phase space. Several measures of complexity, for quantifying the small-scale structures in RPs, have been introduced by Zbilut and Webber Jr. [28]. This ensemble of measures, known as recurrence quantification analysis (RQA), are based on the recurrence point density and the diagonal, vertical and horizontal line structures. In the last decades, RQA has been extended with new measures of complexity [29,30].
In the last years, the transformation of time series in complex networks has become a frequent approach to the analysis of complex dynamical systems. Significant examples are visibility graphs [31,32] and recurrence networks [33][34][35][36]. From this point of view, RPs can be interpreted as the adjacency matrix of an underlying complex network C. This allows the study of time series dynamics by means of the organization of networks, using various specific measures.
In the present approach, we have used two measures, the average edge overlap and the interlayer mutual information, for assessing the similarity of the evolution of different time series.
The average edge overlap, usually noted as x, has been introduced in [37], and it quantifies the simultaneous presence of edges in complex networks fC k g, which are described by the adjacency matrices given by the recurrence plots fRP k g: where d ij is the Kronecker delta symbol and k is indexing the recurrence plots. The interlayer mutual information (IMI) has been introduced also in [37], and it quantifies the presence of correlations between the degree distributions in two different complex networks C a and C b : where Pðk a Þ and Pðk b Þ are the degree distributions of the networks whose adjacent matrices are the recurrence plots RP a andRP b . P k a ; k b À Á is the joint probability to find a node having degree k a in the complex network C a and degree k b in the network C b . The degree distribution PðkÞ of a network is defined as the fraction of nodes in the network with degree k. For a network with N nodes in total, n k of them having the degree k, P k ð Þ ¼ n k N . IMI is indeed a similarity measure which quantifies the information flow between networks. Important insights into the dynamics of a system, when only one realization is available, may be achieved also by using recurrence time statistics [38]. RPs allow the evaluation of the recurrence time (RT). As shown in [39], the white vertical lines in an RP, which corresponds to the vertical non-recurrence points, are an estimator of RT. Various measures of complexity based on recurrence time statistics proved to be powerful tools to detect non-stationary time series and dynamical transitions from regular to chaotic behaviour (see, e.g. [38,40,41]). In the present approach, the changes in the structure of individual RPs have been monitored by two recurrence time statistics measures, the recurrence time spectrum (RTS) and the recurrence time entropy (RTE). RTE is the Shannon entropy of the recurrence times contained in the recurrence matrix R, ruling out the lines shorter than a certain value (here lmin¼ 2) and all the points inside the Theiler window [42]. For perfectly periodic signals RTE = 0, while the maximum value is reached for a pure stochastic signal.

Detection of chaos onset
Transition to chaotic dynamics can be a route to collapse. It is therefore interesting to investigate also this aspect of the dynamics to understand disruptions. The possible onset of the chaotic regime has been evaluated using the 0-1 test for chaos proposed by Gottwald and Melbourne [43][44][45]. A slightly modified version, able to cope with moderate amounts of measurement noise, has been proposed in [46], while a detailed methodology and a MATLAB computer code have been provided in [47]. The 0-1 chaos test has the advantage to work directly on the time series, avoiding the space reconstruction and the calculation of Lyapunov exponents.
The studied time series / is used to construct the following system of equations: where c 2 0; 2p ½ is a random number and is indexing the time series values. For a certain value of c, it is possible to derive the following relations: Gottwald and Melbourne found that for a regular time series /, the motion of p and q is bounded, while if / is chaotic then p and q exhibit an asymptotic Brownian motion. Also, they showed that the timeaveraged mean square displacement of p and q is given by the following relation: where g n 2 À 1 2 ; 1 2 Â Ã is a uniformly distributed random variable and r is the noise level.
The growth rate of M c can be evaluated by means of a correlation coefficient: Usually, K c is calculated for 100 different values of c 2 ½0; 2p and the median value K represents the final output. K will approach 1 for chaotic dynamical systems, while K approaching 0 is specific of periodic systems.

Exploring time series
All the indicators described above have been calculated using a sliding window. As disruptions are an irreversible process in the plasma evolution, the size of the sliding window has been defined based on the concept of time series irreversibility (see, e.g. [48]).
A simple method to characterize (ir)reversibility is based on the third-order statistical moment, suitably standardized in the form of a skewness coefficient [49]. Besides detecting irreversibility, it is important also to determine at which time scale it occurs. Therefore, a temporal multi-scale representation is required. Based on the process U i at scale 1, it is possible to define a process U ðkÞ i , at scale k [ 1 by the relation: where k is a positive integer. Then, the marginal third moment and the skewness coefficient are given by the following equations: where k denotes the scale. The calculations have been performed on the differenced version of the time series in order to ensure stationarity. For differenced time series, the first-order moment is always zero, while the secondorder moment has only positive values, being unable to provide any indication about the time asymmetry. On the contrary, a high magnitude skewness value indicates a strong irreversibility, while values close to zero are related to reversible processes. Figure 2 presents the evolution of the skewness coefficient at different time scales for the LM signal, for all the analysed pulses. In most cases, the irreversible processes take place at scales below 100 ms (which corresponds to approximately 500 time steps in the time series). Therefore, these results suggest that the use of a sliding window of 100 ms, to find distinctive features related to disruptions, is a good choice. Moreover, the use of a window size of 100 ms is also justified by physical reasons. Indeed, this choice provides a sufficient resolution, compared with the characteristic time scales of JET evolution. As some of the indicators could be used in online disruption prediction, a size as small as possible is needed.
4 Summary of the results for the investigated database As described above, all the indicators described in the previous section have been calculated using a sliding window of 500 points, equivalent to 100 ms, circulated with shifts of 1 point over the time series. It should also be mentioned that the time series are composed of several tens of thousands of samples, which ensures the reliability of the third-order moment estimation [50]. Significant changes in the embedding dimension occur before the disruptions for all the analysed pulses, in the case of the plasma current and the locked mode amplitude time series. Figure 3 shows the results obtained for a representative JET discharge (#92045, hybrid scenario). In this figure, the position of the sliding window is indexed by its right margin (so a point in a plot corresponds to calculations performed for the previous 500 points). A very sharp increase in the embedding dimension appears in both I and ML signals, for all pulses, more than 300 ms before the beginning of the current quench. The increased dimensionality persists until the moment when the disruption begins. An overview of the changes in the embedding dimension for the analysed database is provided in Appendix 2.
A significant dynamical change has been recorded also for the chaotic regime indicator, in the case of the ML signal. It shows the onset of the chaotic regime also with more than 300 ms before the disruption. A typical result is presented in Fig. 4 for JET discharge #81159. It is interesting to notice that the chaotic regime ends just before the disruption. The magnitude of the time interval between the end of the chaotic regime and the disruption varies among discharges in the range of 21-354 ms. The duration of the chaotic regime varies between 38 and 395 ms. All the details are presented in Table 3 in Appendix 3. For the I signal, the chaos test does not provide positive results for all pulses. When it occurs, the duration of the chaotic regime is significantly shorter, in the range of  Table 4 of Appendix 3. The existence of a time interval, before the disruption, characterized by significant dynamical variations in the signals, is supported also by the results obtained for the indicators, which characterize the changes in the structure of the recurrence plots. The evolution of the recurrence time spectra shows the existence of three time intervals with distinct behaviour. During a time interval roughly synchronous with the time interval of the chaotic regime onset, a sharp shrinking of the spectra is recorded. This behaviour is illustrated in Fig. 6 for the JET discharge #91570. The shrinking tendency is increasing even more when approaching the moment of the disruption, after the end of the chaotic regime. A similar behaviour is revealed by the evolution of the recurrence time entropy, as illustrated in Fig. 7 for the same JET discharge. Again, three distinct time intervals can be identified. The time series used for calculating RTS and RTE in Figs. 6, 7 is presented in Fig. 8.  The indicators evaluating the similarity between recurrence plots corresponding to ML and I signals show a clear depression in the same time interval when the other indicators have a significant variation. An illustration is shown in Fig. 9, which shows the evolution of the average edge overlap and of the interlayer mutual information for JET discharge #92226. The complete results obtained for these two indicators are given in Tables 5 and 6 in Appendix 4.
As already mentioned, the indicators provide coherent results in identifying the time intervals when significant changes in the signal morphology occur. This is illustrated in Fig. 10, which shows graphically how each indicator identifies the interesting time intervals in the case of the JET discharge #81985, for the ML signal. A complete comparison, between the time when dynamical changes occur, is shown in Fig. 11 for the entire database analysed.

Discussion and conclusion
In this work, various techniques to detect regime variables have been applied to the signals of tokamak diagnostics. The goal of the analysis consists of determining in an objective and statistically sound way when disruptions start leaving a signature in the acquired measurements. The deployed methods are mathematically completely independent and therefore the fact that they provide very coherent results is quite reassuring and gives confidence in the conclusions. The analysis of several discharges in different experiments reveals that the first signs of an approaching disruption are present in the plasma current and locked mode amplitude signals about 300 ms before the beginning of the collapse. Earlier warning or alarms, launched by predictors using the signals considered in the present work, should therefore be considered false alarms. However, 300 ms is a comfortable time interval to undertake mitigation actions. The techniques presented should therefore allow building significantly more reliable training sets for the predictors already providing quite satisfactory performance.
With regard to the plasma dynamics, very interesting new information has emerged. The deployment of these advanced signal processing techniques has revealed that the plasma evolution leading to disruptions typically consists of two phases. First, the dynamics becomes more chaotic or at least less coherent, for a couple of hundreds of ms. Then, in the interval closer to the beginning of the collapse, the signal becomes again more coherent and less chaotic. This behaviour should be further investigated also because it shows analogies with other anomalous events, namely the onset of epilepsy [51].
Unfortunately, most of the analysis techniques deployed are quite demanding in terms of computational resources. The only fast one is the 0-1 chaos test. On usual computers, the routine, provided by the DynamicalSystems.jl software package written in Julia language [52], runs in about 1.2 ms, already compatible with real-time applications (on JET the cycle time of the real-time network is 2 ms). All the other methods require from about 500 ms to 1 s to run. The potential of acceleration measures and implementation using GPUs will be a matter of future investigations.  Data availability Data will be available on reasonable request.

Declarations
Conflict of interest The authors have no conflicts of interest to declare that are relevant to the content of this article.

Appendix 1: Database
This appendix provides some details (see Table 1) about the characteristics of the discharges analysed in  the paper. The choice of the individual discharges has been motivated by the intention to cover a wider range of experimental conditions.