Introduction

In the context of the design of nuclear fusion technologies and devices, Tokamak Energy Ltd. (TE) is conceiving compact fusion power plants with the aim of building a first demonstrative plant by the 2030’s1,2, based on two promising technologies: the Spherical Tokamak (ST) and High Temperature Superconductors. TE is presently operating ST401,2,3, a high field ST, with the mission to advance the spherical tokamak physics basis towards fusion pilot plants/power plants. Recently, ST40 has achieved plasmas with: major radius R0 = 0.4–0.6 m, aspect ratio A = 1.7–2.0, plasma current Ipla = 0.8 MA, toroidal field BT = 2.6 T and elongation k = 2.3. ST40 has the highest current density among the spherical tokamaks, which makes it ideal for studying plasma stability and the possible reasons which might lead to off-normal termination of the scenario (i.e., disruptions). The abrupt termination of the plasma scenario leads to several critical issues such as thermal and electromechanical loads on the structures both whenever it is triggered by intended or unintended and/or unpredictable plasma perturbations. A plasma disruption4,5,6, often occurring as a consequence of transient plasma perturbations7, is a complex phenomenon involving plasma instabilities. A complete plasma thermal and magnetic energy transfer to the surrounding structures occurs on very short timescales, potentially inflicting severe damages to components. The study of the machine disruptive operational space for the evaluation of the disruptions-related effects are among the main drivers of the design of tokamaks components, to ensure that the machines reach the projected lifetime and fulfil their function.

Plasma disruptions typically proceed in two steps. First, a sudden loss of the plasma confinement occurs determining the generation of high thermal loads on the plasma facing components. This phase, known as Thermal Quench (TQ), is characterised by a rapid cooling of the plasma to sub-keV temperatures within very short time scale, and by a global contamination of the plasma after the plasma-wall interactions. The combination of these factors results in a significant increase of the plasma resistivity, leading to a rapid decay of the plasma current, known as Current Quench (CQ) phase. During CQ, the abrupt drop of the plasma current generates high currents in the nearby conducting structures, producing substantial electromagnetic forces on the components. In JET, a plasma current of several MA typically decades on a time scale of ~ 5–50 ms, resulting in forces of several MN on the vacuum vessel6. Disruptions can be clustered in several classes depending on the chain of physical events driving the plasma towards the thermal and the current collapse8. Several approaches have been used to automatically classify disruptions9,10,11 with the aim of developing avoidance strategies by controlling or supressing the disruption precursor for each class. Through this analysis, two general classes of disruptions have been identified4,5: Major Disruptions (MDs) and Vertical Displacement Events (VDEs), depending on the role of the plasma vertical instability within the chain of events. MDs occur with the plasma close to its neutral position (central disruptions) and a loss of vertical position control is experienced after the TQ. On the other hand, VDEs occur when the loss of vertical position control destabilizes the plasma before any appreciable cooling of the plasma core, typical of elongated plasma scenarios. During a VDE, the plasma moves vertically; it goes away from its equilibrium position towards the wall, and reduces progressively its area with a little change in the total plasma current. Then, when the plasma starts to significatively interact with the structures or when plasma reaches a critical value of the safety factor (q95)8, conditions for a rapid growth of Magneto-HydroDynamics activities inevitably arise and the TQ occurs.

As an example of a VDE, Fig. 1 shows some snapshots taken from the ST40 internal visible light camera, where the plasma separatrix (blue solid line) inside the inner vacuum chamber highlights upwards movement of the plasma column during the pulse #7153.

Fig. 1
figure 1

Image of a VDE taken from an internal visible light camera of ST40 with the plasma separatrix snapshot (blue solid line) inside the inner vacuum chamber.

This paper adopts the disruptions classification reported in5 and, following the classification approach proposed in7, proposes a data driven algorithm for automatic plasma pulses clustering, based on the deviation of the vertical position with respect to the reference one. Thus, a DataBase (DB) of disruptions and regularly terminated pulses selected from ST40 2021–2022 experimental campaign was built to map the ST40 operational space using a Machine Learning (ML) algorithm. The mapping of the operational space is then used to extrapolate unexplored plasma configurations, whose reliability is investigated by a suitably calibrated ST40 numerical model. Finally, the new pulses of interest are reproduced with a time evolutionary equilibrium code. In this way, it is possible to perform preliminary predictive simulations aimed at investigating the plasma vertical displacement, to provide lessons to be learnt for the next ST40 experimental campaign and for the design of future ST devices.

This paper is structured as follows: first, the criteria and the main features used to build and to populate the ST40 experimental DB, alongside the algorithms for the automatic classification and the subsequent database population are reported. Then, the results in terms of ST40 operational space mapping with the thresholds to achieve a regularly terminated pulse or disrupted pulses are shown. The ST40 numerical model is described as well as the numerical validation of the classification in MAXFEA environment12,13. Finally, the projection of the multi-dimensional experimental database into a 2D map by using ML algorithms is presented. The 2D map allows to visualize a clustering of the experimental dataset and to perform preliminary predictive simulations on the plasma pulse termination trajectory.

Preliminary analysis for populating the ST40 experimental database

The 2021–2022 ST40 experimental campaign under analysis was based on Low-Confinement Mode pulses with limiter-type plasmas14. Within the campaign, a DB of disruptions and regularly terminated pulses was built with the aim of selecting MDs and VDEs for studying the dynamic of the plasma vertical displacement and for predicting its consequences in the most dangerous cases. To this purpose, a data driven algorithm to automatically classify disruptions in MDs and VDEs is proposed in this paper to provide a well-defined identification criterion and to avoid manual classification for the disruption studies on future campaigns. The first issue to be faced in developing a data driven algorithm to automatically classify disruptions is the creation of a reliable target dataset needed to optimize the decision-making metrics. To begin this study, we created a “target database” where we manually classified pulses into two classes depending on whether the plasma position control is kept controlled or not before the end of the flat-top phase.

First, a preliminary list of relevant pulses was prepared manually. The execution of a manual classification is not a trivial task when dealing with raw experimental data and when outliers and incorrect or missing measurements may affect the discrimination values. The pulses selection was performed according to the plasma current evolution and plasma column movement by checking the agreement between their experimental evolutions and the preprogrammed signals, taken as reference. To this end, it is first necessary to identify the flat-top phase according to the pulse termination behaviour (disrupted or non-disrupted). For a regularly terminated pulse, the plasma current flat-top duration is defined by two times: the Start of Flat-top (SOF) time (tsof), corresponding to the end of the ramp up phase (in ST40, it is the Merging Compression start-up15), and the End of Flat-top (EOF) time (teof), corresponding to the beginning of the Ramp Down phase (RD). In a disrupted pulse, i.e., when a disruption occurs during the flat-top, its ending time is the disruption time (td), which is assumed within the TQ phase. In Fig. 2, two examples of the characteristic times detected for a disrupted and a non-disrupted pulse are reported. A semi-automatic identification of the characteristic times was performed. The following relevant times were collected for each disrupted pulse (see Fig. 2a), besides the flat-top start and end times:

  • tEQ: an equilibrium time during the flat-top phase when the plasma is still in its central position and the effects of the plasma start-up are negligible.

  • tCQstart: current quench starting time.

  • tbtd: time taken shortly after the disruption time, useful to evaluate the plasma parameters variations experienced at the beginning of the CQ phase.

As regards the non-disrupted (or regularly terminated) pulses, the semi-automatic analysis collects the following relevant time (see Fig. 2b), besides the flat-top start and end times:

  • tEQ: an equilibrium time during the flat-top phase approximatively near the half of its duration when the effects of the plasma start-up are negligible.

Fig. 2
figure 2

Examples of characteristic times detected respectively for a disrupted (a) and a non-disrupted plasma pulse (b). Red solid lines report the plasma current (Ipla). Magenta and green solid lines represent the actual plasma vertical (Zc) and radial (Rc) positions, respectively. Purple dashed lines indicate the SOF time (tsof). Green dashed line indicates the considered equilibrium time during the flat-top (tEQ). Cyan dashed line indicates the EOF time (teof) and the disruption time  (td) for non-disrupted and disrupted pulses, respectively. Dark red dashed line reports the Current Quench start time (tCQstart) and blue dashed line reports a time during the CQ taken shortly after its starting time (both for disrupted pulses only).

Following the definition of the characteristic times, the preliminary analysis allowed us to identify, with high confidence, the following two classes:

  • VDE: it is the so-called hot VDE, where a large plasma vertical displacement (several centimetres) induces a disruption (Fig. 2a from shot #9054).

  • Non-Disrupted or regularly terminated pulses (ND): plasma pulses where there is no disruptive event or there is a possible minor disruption that does not lead to an anomalous termination (Fig. 2b from shot #9501).

A preliminary list of 432 relevant shots was generatedc starting from a set of over 2000 experimental pulses. Among them, 95 pulses show a large deviation of the vertical position with respect to the reference and were classified as VDEs; 73 non-disrupted pulses were selected to represent the class of pulses where the position control is maintained until the end of the plasma current flat-top. For the remaining 264 pulses, the preliminary analysis did not result into an unambiguous classification because of the lack of evident features typical of a disruption or of a regularly terminated pulse. The 95 VDE pulses and 73 non-disrupted pulses were then used to optimize an automatic classifier, used to distinguish among the remaining ambiguous 264 pulses.

ST40 experimental database automatic analysis

The automatic classifier consists of a cascade of two automatic binary classifiers. The first one is a 2D binary classifier able to distinguish between disrupted and non-disrupted pulses by thresholding the values of the plasma current time derivative and the plasma energy variation. After selecting disrupted pulses, a 1D binary classifier automatically distinguishes between MDs and VDEs by thresholding a feature based on the deviation of the plasma vertical position. The first classification is performed by considering the pulse termination behaviour, depicted by the two clustering features, in a reference time window. The reference time window starts at tCQstart and teof for disrupted and non-disrupted pulses, respectively, and ends when the plasma current reaches the value of 30% of its flat-top value, in both cases. In this paper, the automatic clustering of pulses into non-disrupted and disrupted categories is accomplished by thresholding the plasma current time derivative and the plasma energy variation (2D binary classifier). To this purpose, two threshold values for each feature were optimized by minimizing the classification error of VDEs and non-disrupted pulses previously classified manually. When developing an automatic clustering algorithm, the manual classification reflects on the classification criterion and its performance. Therefore, only the unambiguous pulses classified from the manual analysis can be used to define the automatic classification criterion. VDEs and non-disrupted pulses allow to define two threshold values for each clustering features, which separate the disrupted and non-disrupted pulses with 100% of confidence. The lowest values of the plasma current time derivatives and plasma energy variations of the VDE class define the thresholds that the ambiguous pulses must exceed to be classified as disrupted. Since the manual analysis did not classify these pulses as VDEs, they are labelled as disruption with fast CQ. The highest values of plasma current time derivatives and plasma energy variations within the non-disrupted class give the thresholds not to be exceeded for classifying a pulse as non-disrupted or with a late disruption (defined as a disruption occurring after the plasma current drops below 30% of its flat-top value). Additionally, pulses characterised by the plasma current time derivative and the plasma energy variation falling between the two thresholds for each feature are classified as disruptions, but specifically characterised by a slow CQ. Therefore, the developed algorithm allows to classify pulses as: non-disrupted, disrupted with fast CQ, and disrupted with slow CQ. Note that, in the scientific literature, the classification of disruption events as fast or slow is distinctly related to the duration of the plasma current drop phase, which is influenced by the size of the machine, as shown in4,8. However, the label “fast” and “slow” are used here to differentiate between the two types of disruptions identified by the automatic algorithm described above.

Using the automatic classification, the 432 pulses were classified as: 101 regularly terminated pulses and 300 disruptions, while 31 pulses were discarded because of kinetic energy data lack during the CQ. Considering the 300 disrupted pulses in more detail: 204 were classified as fast CQ and 96 as slow CQ. To distinguish between MDs and VDEs, a new variable was defined, DevZc,max, which is the maximum deviation of the plasma vertical position, Zc, in the flat-top time window with respect to Zc at SOF. Pulses where DevZc,max exceeds some threshold value, DevZc,thr, are classified as VDEs by the 1D binary classifier, otherwise, they are classified as MD (i.e., VDEs have DevZc,maxDevZc,thr and MDs have DevZc,max < DevZc,thr). The optimum threshold was found to be 2.874 cm. Using this classification we found: 95 VDEs from 73 non-disrupted pulses in the target dataset. Only three non-disrupted pulses were misclassified, leading to a confidence level of 98.21%. The MDs vs. VDEs classifier was applied to the 300 disrupted pulses populating the ST40 DB. Figure 3 shows the deviation of the actual plasma vertical position Zc for the 300 disruptions; the horizontal black line indicates the optimized threshold value, the pulses below the threshold (blue circles) are classified as MDs whilst the pulses above the threshold (red circles) are classified as VDEs. Finally, the ST40 DB collects: 118 VDEs (39.3% of the disruptions) and 182 MDs (60.7% of the disruptions), and 101 non-disrupted pulses.

Fig. 3
figure 3

Disruptions classification into MD and VDE for the ST40 DB based on DevZc,max: VDEs—red dots (118 pulses), MDs—blue dots (182 pulses). The optimized threshold (2.874 cm) is indicated by the horizontal black line.

Vertical displacement characterisation of ST40 pulses

After populating the ST40 DB, a vertical displacement analysis of ST40 pulses was implemented considering the plasma internal parameters variations. This is aimed at defining a preliminary safe operating space and the correlation between the vertical displacement \(\:\left(\varDelta\:{Z}_{\text{c}}\right)\) and the poloidal beta variation \(\:(\varDelta\:{\beta\:}_{\text{p}})\)6 during a prefixed time window, when it is possible to highlight differences between disrupted and regularly terminated pulses. This study is the starting point for the comparative simulations reported in the next section. Figure 4 reports a scatter plot between \(\:\varDelta\:{Z}_{\text{c}}\:\)and \(\:\varDelta\:{\beta\:}_{\text{p}}\), together with a qualitative operational space clusterization for each considered pulse type.

Fig. 4
figure 4

Scatter plot with correlation between \(\:\varDelta\:{Z}_{\text{c}}\:\)and \(\:\varDelta\:{\beta\:}_{\text{p}}\) in the time windows of interest, and qualitative clusterization of regions of the operational space: ND region—green, MD region—blue, VDE region—red. Cross-shaped markers report the position of three numerically reconstructed plasma pulses.

During a disruption event, the magnetic diagnostics of interest turn out to be incomplete and characterised by a reconstruction that stops before the disruption time in 59% of the pulses. This failure is due to the intensity of the event itself which saturates the magnetic measurements. Indeed, these missing equilibrium times were mostly found for VDEs (more than 90% of the pulses were discarded), resulting in the unavailability of the \(\:\varDelta\:{Z}_{\text{c}}\) and \(\:\varDelta\:{\beta\:}_{\text{p}}\) signals for 90% of the samples. In addition to a preliminary cleaning focusing on the goodness of the signal reconstruction, a statistical analysis on \(\:\varDelta\:{Z}_{\text{c}}\) and \(\:\varDelta\:{\beta\:}_{\text{p}}\) in the mentioned time windows was performed to identify and exclude outliers. Specifically, the mean value and the lower and upper bounds of the distribution indicating the 25th and 75th percentiles, respectively, were calculated. Samples exceeding the percentiles were considered outliers and were excluded from the analysis (a total of 25 pulses were excluded). As expected, a high concentration of ND pulses around a \(\:\varDelta\:{\beta\:}_{\text{p}}\) close to zero is observed in Fig. 4. In this way, it is possible to establish a qualitative range of variation for the poloidal beta and the plasma centroid vertical displacement to separate the NDs region from MDs and VDEs regions.

ST40 numerical model calibration via experimental pulses from ST40 database

To perform preliminary predictive analyses of the plasma displacement, a ST40 numerical model was set up in MAXFEA environment and then validated by comparing the numerically reconstructed equilibrium and the subsequent plasma time evolution with experimental data provided by the diagnostics. MAXFEA is a 2D electromagnetic equilibrium code that solves the Grad-Shafranov non-linear problem with a Finite Element Method (FEM) approach12. Starting from an initial reference equilibrium, MAXFEA calculates the plasma evolution considering electromagnetic coupling with both passive structures and active circuits. The data used for the equilibrium reconstruction and dynamic simulations come from signals of plasma internal parameters and the evolving currents within active circuits, as provided by the ST40 magnetic diagnostics. The ST40 main components considered in the model are: the Inner and Outer Vacuum Chambers (IVC and OVC) and BVL, BVM, BVT and DIV poloidal field copper coils3,13, as shown in Fig. 5.

Fig. 5
figure 5

Overview of ST40 main components and PF coils (a, b) and detailed view of divertor system components (c)4,13.

More specifically, OVC is a Stainless Steel (SS) 316 structure housing the PF and TF coils, and the IVC, which is made in SS314 with the central tube in Inconel 718. E-rings made of SS316 are placed on the OVC3. Gas Baffles (GASBFL) made of Inconel 718 are also located inside the IVC3. Recently, two up-down symmetric divertors were introduced in ST40 (Fig. 6a) to enable Double Null Diverted (DND) operations14. They are complex structures consisting of: (i) passive stabilization copper plates, namely the Divertor Stabilizing Plates (DIVPSR red areas in Fig. 6b) and the High Field Side Stabilizing Plates (HFSPSR blue areas in Fig. 6b); (ii) castellated graphite plasma facing components on the divertor plates connected with the IVC through SS supports (Fig. 6b); (iii) poloidal limiters on the stabilizing plates and (iv) steel support structures. Since MAXFEA is a 2D axisymmetric numerical procedure, the divertor 3D structure was simplified and suitably calibrated by modelling the DIVPSR and the HFSPSR as full toroidal continuous structures. Moreover, their equivalent conductivity was calculated in such a way that the numerically computed eddy currents are the closest possible to the experimental signal given by GASBFL Rogowski coils.

Fig. 6
figure 6

Overview of ST40 geometrical model: (a) ST40 top divertor system. (b) ST40 poloidal section model with indication of the main passive structures and layouts of active circuits (PF coils).

The model so calibrated was used to first reproduce equilibrium configurations at the equilibrium time reported in the ST40 DB. As an example, Fig. 7 shows a comparison between the experimental and numerically reconstructed plasma Last Closed Flux Surfaces (LCFS) for three plasma pulses. A good agreement can be found not only in terms of plasma shape and position but also in terms of plasma internal parameters (Ipla, \(\:{\beta\:}_{\text{p}}\), li), reproduced within a tolerance lower than 1%.

Fig. 7
figure 7

Reference starting equilibria reproduced in MAXFEA—blue LCFS; and experimental plasma boundary—red boundary: (a) VDE, (b) MD, (c) ND.

The plasma electromagnetic time evolutions were modelled referring to the corresponding experimental events depicted in Table 1.

Table 1 Electromagnetic modelling strategies of disrupted and regularly terminated pulses.

Both disrupted and regularly terminated pulses were simulated according to the modelling strategy in Table 1, as shown in Figs. 8a, 9a and 10a for each pulse considered. Comparisons between experimental and numerical data are reported in terms of plasma current and \(\:{\beta\:}_{\text{p}}\) variations used to drive the numerical analysis. Numerical simulation results are shown in Figs. 8b, 9b and 10b in terms of plasma trajectory (Rc and Zc) during the scenario, alongside a comparison with experimental reconstructed signals for the disrupted and the regularly terminated pulses. In addition, the poloidal projection of the plasma centroid position time evolution and plasma separatrices at relevant times are reported in Figs. 8c, 9c and 10c for VDE, MD and ND, respectively.

Fig. 8
figure 8

Comparison between numerical and experimental data for VDE #9054: (a) numerical and experimental plasma current– continuous blue and dashed magenta lines respectively, and \(\:{\beta\:}_{\text{p}}\)—continuous cyan and dashed green lines respectively; (b) experimental (dashed lines) and numerical (continuous lines) plasma current centroid position time evolution—blue and magenta lines for radial displacement, cyan and green lines for vertical displacement; (c) plasma current evolution and trajectory of the centroid on the poloidal plane—magenta line, together with the LCFS evolution at different times—in black at tEQ, in blue at the First-wall Touch (FT) time, in red at TQ time, in green and in cyan at 70% and 30% of Ipla during the CQ respectively.

Fig. 9
figure 9

Comparison between numerical and experimental data for MD #9820: (a) numerical and experimental plasma current—continuous blue and dashed magenta lines respectively, and \(\:{\beta\:}_{\text{p}}\)—continuous cyan and dashed green lines respectively; (b) experimental (dashed lines) and numerical (continuous lines) plasma current centroid position time evolution—blue and magenta lines for radial displacement, cyan and green lines for vertical displacement; (c) plasma current evolution and trajectory of the centroid on the poloidal plane—magenta line, together with the LCFS evolution at different times—in black at tEQ, in blue at the First-wall Touch (FT) time, in red at TQ time, in green and in cyan at 70% and 30% of Ipla during the CQ respectively.

Fig. 10
figure 10

Comparison between numerical and experimental data for ND #9501: (a) numerical and experimental plasma current—continuous blue and dashed magenta lines respectively, and \(\:{\beta\:}_{\text{p}}\)—continuous cyan and dashed green lines respectively; (b) experimental (dashed lines) and numerical (continuous lines) plasma current centroid position time evolution—blue and magenta lines for radial displacement, cyan and green lines for vertical displacement; (c) plasma current evolution and trajectory of the centroid on the poloidal plane—magenta line, together with the LCFS evolution at different times—in black at tEQ, in blue at the First-wall Touch (FT) time, in red at EOF time.

Despite all the assumptions for each pulse, Figs. 8, 9 and 10 show a good agreement in terms of plasma current centroid trajectory (especially for the vertical displacement). The electromagnetic model proves to be reliable and well calibrated from the comparison of the numerical results with experimental data. In fact, the poloidal beta drop perturbation produces a vertical displacement in a good agreement with the experimental observation shown in Fig. 4. Each plasma pulse falls into the corresponding region of the operational space and is very close to the experimental data point (the non-perfect matching is due to limitations in MAXFEA numerical procedure). In particular:

  • In pulse #9054 (VDE), the imposed \(\:\varDelta\:{\beta\:}_{\text{p}}\) is \(\:\sim\)0.328 and \(\:\varDelta\:{Z}_{\text{c}}\sim-17.3\:\text{c}\text{m}\) was obtained in the time window from tEQ to tbtd;

  • In pulse #9820 (MD), the imposed \(\:\varDelta\:{\beta\:}_{\text{p}}\) is \(\:\sim\)-0.370 and \(\:\varDelta\:{Z}_{\text{c}}\sim-0.54\:\text{c}\text{m}\) was obtained in the time window from tEQ to tbtd;

  • In pulse #9501 (ND), the imposed \(\:\varDelta\:{\beta\:}_{\text{p}}\) is \(\:\sim\)-0.034 and \(\:{\Delta\:}{Z}_{\text{c}}\:\sim0.61\:\text{c}\text{m}\) was obtained in the time window from tEQ to teof.

Machine learning clustering applications to ST40 DB for preliminary predictive simulations

Machine Learning is a branch of Artificial Intelligence encompassing classes of mathematical models that can be classified into three main learning paradigms: supervised learning, unsupervised learning, and reinforcement learning16. These paradigms differ both in the tasks they can accomplish and in how the data is presented. Reinforcement learning is used to teach agents to learn a given task by interacting with an external environment; no fixed set of training examples is required16. On the other hand, algorithms performing pattern recognition, classification, prediction and dimensionality reduction tasks learn from a training dataset by using two different approaches: supervised and unsupervised17.

Supervised learning implies a priori knowledge of the input labels; it can be used to classify new and unseen datasets and predict outcomes. Unsupervised learning uses unlabelled input data; it is usually adopted to better understand the datasets by identifying unknown patterns and trends, or to cluster similar data into a specific number of groups. Among ML algorithms, mapping methods are employed to learn or uncover underlying patterns embedded in the data manifold of multidimensional input spaces. In addition, 2D or 3D visualizations of the multidimensional input space can be achieved through dimensionality reductions. Unsupervised mappings algorithms were already used in fusion research to map tokamaks operational spaces in the framework of disruption prediction and classification18,19,20,21. Self-Organizing Map (SOM)22 is an unsupervised ML algorithm which performs both clustering and mapping of a multidimensional input space into another one with reduced dimensions (dimensionality reduction), according to the similarity between the data. The input space dimension is reduced while preserving the mutual topological position of the clusters, which are commonly displayed in a 2D or 3D grid, for visualization purposes. A 2D SOM consist of a single-layer neural network of \(\:M\) neurons arranged in a \(\:\left(n\times\:l\right)\) grid. Let \(\:X\) be a \(\:D-dimensional\:\) input space of data, \(\:X=\left\{{x}_{1},\dots\:,\:{x}_{N}\right\}\), where \(\:N\) is the number of the data sample; every neuron in the grid is defined by a \(\:D\)-\(\:dimensional\) weight vector in the space \(\:W=\left\{{w}_{1},\:.,\:{w}_{M}\right\}\). The weight vector space \(\:W\) is randomly initialized and trained by a competitive unsupervised learning. In each training step, the \(\:{i}^{th}\) sample vector \(\:{x}_{i}\) from the input space is chosen randomly and a similarity measure is calculated between it and all the \(\:M\) weight vectors of the grid. Then, the closest neuron is assigned to the input vector as the winning neuron. The weight of the winning neuron and the ones contained in its neighbourhoods are updated to minimize the distance from \(\:{x}_{i}\), according to Eq. 1:

$$\:{\varvec{w}}_{j}(k+1)\:=\:{\varvec{w}}_{j}\left(k\right)\:+\:\alpha\:{h}_{jc}\:\left(r\right)\left[d\left(x,\:{w}_{j}\left(k\right)\right)\right]$$
(1)

where \(\:h\) is a kernel function defining the neighbourhood (it is required to be symmetrical and to decrease with the distance to the winning neuron), \(\:\alpha\:\) is the learning rate, \(\:d\) is the similarity metric, and \(\:r\) is neighbourhood radius22. The similarity metric and the kernel neighbourhood function usually adopted are the Euclidean distance and the gaussian function23,24, respectively, as in this paper.

To reveal underlying similarities among the pulses into the DB and to perform preliminary predictive analyses of plasma termination and vertical displacement, a SOM was trained to map the ST40 operational space defined by 24 plasma parameters for 401 pulses. Each pulse is represented by the sample at the equilibrium time tEQ and is treated as a vector of the 24D input space, defined by the following components:

  • Plasma current and current in the 6 active coil circuits.

  • Eddy current in the main in-vessel components (4 divertor plates).

  • Poloidal beta, internal inductance, toroidal field at the machine radius, and radial and vertical coordinates of the magnetic axis on the poloidal plane.

  • Radial and vertical position of 4 plasma-wall gaps.

Despite the ST40 experimental campaign under analysis was based on limiter-type plasmas, there is to specify that Single Null Diverted (SND) configurations were observed in most of the VDEs, in some MDs and in a very limited number of ND pulses. These configurations resulted from high eddy currents on the divertor HFS plates during some pulses due to the low resistivity and their mutual inductance with the plasma. These undesired effects seem to lead to these unintended SND configurations that are poorly controlled by the vertical stabilization system, leading to VDE-type disruptions in most cases. Therefore, a preliminary filtering of the database was performed, keeping the SND equilibria only for those pulses evolving in VDEs, and the limiter equilibria only for those pulses regularly terminating. As for MDs, both types of plasma configuration were considered. The regular lattice in Fig. 11 gives a 2D representation of the 24D clustering performed during the training phase, where the clusters colour (namely, clusters labelling) is given a posteriori, depending on their composition in terms of pulse classes. Red, blue and green clusters contain only one class of pulses: VDE, MD, and ND respectively, grey clusters collect samples of at least two classes, whilst white clusters are empty. The clusterization procedure allows to extrapolate a new plasma configuration for each cluster, not contained in the experimental dataset, and defined by the cluster centroid. Since each plasma centroid is defined by the 24 features defining the SOM input space, these features were assumed as new plasma configurations (in what follows, centroid equilibria) from which to predict the plasma evolution of a VDE, MD, and a non-disrupted pulse. The choice of the clusters is general and any of them may be selected for the purpose. However, those belonging to a macro-aggregation of clusters of the same type and highly density populated are to be preferred given the higher statistical representativeness of their class. Following this criterion, the clusters with yellow diamonds shown in Fig. 11 were selected for extrapolating the new ND, MD, and VDE equilibria. Fig. 12 shows the corresponding plasma separatrices.

Fig. 11
figure 11

SOM trained with equilibria taken at tEQ of MD, VDE and ND pulses collected into the ST40 DB. The clusters colour depends on its composition in terms of sample classes: only VDEs are red, only MDs are blue and only NDs are green, while empty clusters are white, and clusters containing more than one class are grey.

Fig. 12
figure 12

LCFS (in blue) reconstructed for the centroids equilibrium of a VDE cluster (a), MD cluster (b) and ND cluster (c). Respective clusters are indicated with yellow diamond in Fig. 11. Main plasma-wall gaps are marked with red crosses.

ST40 preliminary predictive simulations

Preliminary predictive transient simulations can be performed once the three starting equilibria were obtained from the centroid features, and the following weighted average parameters were calculated:

  • for regularly terminated pulses: duration of the time window Δt = tEQ÷teof, βp variation in the time window Δt, and RD duration;

  • for disrupted pulses: duration of the time window Δt = tEQ÷td, βp variation in the time window Δt, q95 at td, and TQ and CQ duration.

Note that the calculation of the weighted average value considers that the centroid position is influenced not only by the pulses within that cluster but also by all pulses belonging to the same class. Therefore, the weights are chosen as the reciprocal of the Euclidean distance of the centroid from all the pulses belonging to the same class. Finally, predictive analyses were carried out by modelling the plasma evolution, as reported in Table 2:

Table 2 Electromagnetic modelling strategies of predictive simulations for disrupted and regularly terminated pulses.

In the disrupted clusters examples, the \(\:{\beta\:}_{\text{p}}\) variation is set to \(\:\sim\)0.199 in about 70 ms and \(\:\sim\)0.061 in about 21 ms for the VDE and MD, respectively. Afterwards, TQ is imposed when the trigger condition is met (q95 limit is equal to \(\:\sim\)4.7 and \(\:\sim\)5.6 for the VDE and the MD, respectively) and a current spike is found due to the flux conservation and a subsequent redistribution of the plasma current, as observed experimentally. In the regularly terminated cluster example, it is possible to observe a quasi-flat plasma current evolution with no disruptions after the \(\:{\beta\:}_{\text{p}}\) perturbation of \(\:\sim\)0.074 in about 52 ms, as expected.

Figures 13, 14 and 15 report the \(\:{\beta\:}_{\text{p}}\) time evolution and the transient analyses results in terms of plasma current evolution for the clusters selected in the SOM shown in Fig. 11. Specifically, plasma current, \(\:{\beta\:}_{\text{p}}\) and radial and vertical plasma centroid position time evolution are reported in Figs. 13a, 14a and 15a for VDE, MD and ND, respectively. Poloidal projection of the plasma centroid position time evolution and plasma separatrices at relevant times are reported Figs. 13b, 14b and 15b for VDE, MD and ND, respectively.

Fig. 13
figure 13

Centroid VDE: (a) reconstructed plasma current—blue, radial and vertical plasma centroid position—red and green lines respectively, and \(\:{\beta\:}_{\text{p}}\) imposed—magenta; (b) plasma current evolution and trajectory of the centroid on the poloidal plane—magenta line, together with snapshots of LCFSs and at different times—in black at tEQ, in blue at the First-wall Touch (FT) time, in red at TQ time, in green and in cyan at 70% and 30% of Ipla during the CQ respectively.

Fig. 14
figure 14

Centroid MD: (a) reconstructed plasma current—blue, radial and vertical plasma centroid position—red and green lines respectively, and \(\:{\beta\:}_{\text{p}}\) imposed—magenta; (b) plasma current evolution and trajectory of the centroid on the poloidal plane—magenta line, together with snapshots of LCFSs and at different times—in black at tEQ, in blue at the First-wall Touch (FT) time, in red at TQ time, in green and in cyan at 70% and 30% of Ipla during the CQ respectively.

Fig. 15
figure 15

Centroid ND: (a) reconstructed plasma current—blue, radial and vertical plasma centroid position —red and green lines respectively, and \(\:{\beta\:}_{\text{p}}\) imposed—magenta; (b) plasma current evolution and trajectory of the centroid on the poloidal plane—magenta line, together with snapshots of LCFSs and at different times—in black at tEQ, in blue at the First-wall Touch (FT) time, in red at EOF time.

The results from the preliminary predictive analysis perfectly fall within the respective operational region shown in Fig. 4. In fact, a large downward vertical displacement (~ 5.5 cm) is expected by the plasma for the VDE after the \(\:{\beta\:}_{\text{p}}\) perturbation of about 0.199, and a much smaller vertical displacement (~ 0.03 cm) is expected for the MD after the \(\:{\beta\:}_{\text{p}}\) perturbation of about 0.061, before the start of the TQ. On the other hand, a small vertical displacement of the plasma centroid (~ 0.73 cm) is expected before the ramp down phase for the ND pulse because of the \(\:{\beta\:}_{\text{p}}\) perturbation of about 0.074.

This analysis can provide preliminary predictive information on the evolution of plasma pulses not yet executed. In fact, the exploitation of the “similarity” among the pulses belonging to the same cluster allows to classify new equilibria by observing where they fall within the SOM. Then, the performance of numerical analyses starting from new equilibria and the input \(\:{\beta\:}_{\text{p}}\) variation calculated for that cluster’s centroid allows estimating the subsequent \(\:\varDelta\:{Z}_{\text{c}}\), providing a new numerical data point to be inserted into the operational space shown in Fig. 4. Therefore, the study can provide first indications about how the unexecuted plasmas are expected to evolve from that equilibrium configuration. It provides preliminary predictions on the plasma behaviour and allows the scientific coordinator to redesign it before its execution whenever the pulse does not meet the design specifications or terminates with a disruption.

Conclusions

In order to support the ST40 Tokamak operations, preliminary predictive analyses aimed at understanding the evolution of plasma displacement and the possibility of ending in disruptions were performed in this paper. The analyses were supported by the definition, construction, and populating of an experimental database based on the plasma pulses performed during the 2021–2022 experimental campaign. The DB collects 101 non-disrupted and 300 disrupted experimental pulses deriving from more than 2000 experiments, showing a plasma disruptivity of about 12.3%. Their classification into MDs and VDEs was performed by using data-driven algorithms based on the deviation of the plasma centroid vertical position with respect to the reference one. The automatic procedure resulted in 118 VDEs (39.3% of the disruptions) and 182 MDs (60.7% of the disruptions). Based on this classification, a qualitative operational space was built in terms of \(\:\:{\beta\:}_{\text{p}}\) range of variation and plasma centroid vertical position to separate the NDs region from MDs and VDEs regions. The aim is to provide first indications of the maximum allowable variations to avoid the plasma pulse ending in a disruption. The application of the Self-Organized Map Machine Learning algorithm allowed the identification of hidden similarities among the pulses contained in the database, through a clusterization based on common features and the calculation of a representer (namely, the centroid) for each cluster. The centroid represents the cluster composition and provides a new configuration extrapolated from the experimental ones.

In order to support this research, an electromagnetic FEM model of ST40 was developed in MAXFEA environment and was validated against experimental data measured by the magnetic diagnostics. The numerical model is able to provide preliminary predictions on the evolution of unexecuted plasma pulses starting from the new equilibrium condition and the experimental perturbation of \(\:{\beta\:}_{\text{p}}\), extrapolated from experimental observations in a prefixed time window specific for the cluster to which the equilibrium belongs. In this way it is possible to observe whether the numerical transient analysis results in a disruption (central MD or VDE) or in a non-disruptive pulse. Without lacking generality, the study of an equilibrium was performed, starting from the centroid of a cluster selected from the SOM for each pulse class, based on its numerosity. The numerical transient analysis confirmed the class of belonging and provided preliminary information on the evolution of that pulse before its eventual execution. The analysis presented in this work is preliminary since a number of other variables characterising the pulse behaviour were not accounted due to the lack of reliable signal diagnostics in most cases.

It should be stressed that this methodology is general and can be applied to any other existing machine. Indeed, the exploitation of the common characteristic among the pulses belonging to the same cluster allows the classification of new proposed equilibria from future experimental campaigns according to the SOM. In this way, the performance of new numerical analyses starting from future experimental equilibria and the input \(\:{\beta\:}_{\text{p}}\) variation calculated for the centroid of these clusters allows the estimation of \(\:\varDelta\:{Z}_{\text{c}}\), providing new numerical data points to be inserted into the operational space. In conclusion, the study can provide a predictive tool that can estimate the plasma evolution in terms of the displacement that is triggered in the presence of specific parameters variations, giving indications on the plasma behaviour and the presence or absence of possible disruptive occurrences. In the future, it is planned to extend the population of the DB and to validate this methodology with data coming from plasma pulses executed under new ST40 experimental campaigns, after a proper recalibration of the electromagnetic model, necessary due to the modification of the divertor plates that took place in 2023.