Investigation of the J-TEXT plasma events by k-means clustering algorithm

1.   Introduction

Supervised learning and semi-supervised learning algorithms have been widely applied to prediction and fitting in fusion plasma research [1–8], like disruption prediction [1–3], inversion reconstruction [4], simulation [5, 6], and event recognition [7, 8]. These algorithms learn from abundant diagnostic data marked with plasma events, like plasma disruption or L-H transition, and then repeat similar behaviors by given data, avoiding complicated or incomprehensible dynamical methods. The hierarchical clustering and k-mean method was implemented in the National Spherical Torus Experiment to identify characteristic evolution patterns of an edge localized mode event [8], indicating that clustering analysis is a powerful tool for classification and identification of events.

The k-means clustering algorithm [9], as a prototype-based cluster analysis method, can classify samples into diverse clusters by the prototype points (like weighted center points). Samples of the same type, gathering around center points, are divided into isolated clusters, which is similar to classification and identification of plasma events which occur in diverse parameter regions. Hence, this algorithm is expected to identify and label various plasma activities to avoid using the scientist's time to label abundant and various plasma events. Based on the above thinking, a semi-supervised learning algorithm, the k-means cluster analysis method, is utilized to survey the diagnostic data of the J-TEXT plasma events. The k-means clustering algorithm can cluster plasma events having similar features, which can be labelled by manually marking a few data. Applying the results of the clustering analysis, massive diagnostic data can be labelled with plasma events, which is beneficial to physics research on macroscopic magnetohydrodynamics (MHD) activities.

In this work, the k-means algorithm is utilized to conduct diagnosis of signals on the J-TEXT device, and some interesting results have been obtained that several types of plasma events can be classified by this algorithm, and then a simple analytical model is proposed to describe these behaviors. This article is organized as follows: section 2 induces the k-means algorithm, section 3 depicts the classification and a simple analytical model for identification of plasma events, and the last section is a summary of this work.

2.   Methodology

Assume a dataset X = {x₁, x₂, …, x_i, …, x_n}, where x_i stands for a multi-dimensional sample, is a collection R = {R₁, R₂, …, R_i, …, R_k}. R_i is a cluster made up of certain data in X, and it can represent certain plasma events, like the rotating tearing mode (RTM), sawtooth oscillations (ST) or locked mode (LM). The k-means algorithm can classify X into k clusters by comparing the distances between x_i and c_j in the corresponding center collection C = {c₁, c₂, …, c_j, …, c_k} of R, and the center c_j is determined by

Here, we employed the Euclidean distance to estimate the distance between x_i and c_j,

where x_{i, l} and c_{j, l} are the sub-parameters of L-dimensional x_i and c_j, respectively. The criterion is that x_i belongs to the cluster whose center is the nearest to x_i.

The center collection C can be determined by the iterative method: (a) randomly or manually select k elements in X as cluster centers in C; (b) determine distances between each point and k cluster centers by equation (2), divide all points into their nearest cluster and then update the centers by equation (1); (c) repeat step (b) until the centers have no change or the sum of squared error (SSE) reaches its minimum. Here, SSE is written as,

When the iteration terminates, we can finally derive the center collection C and event collection R. Given a few marked samples, the clusters in R can be labeled with diverse plasma events on the J-TEXT tokamak. J-TEXT is a medium size tokamak with a major radius R = 1.05 m and minor radius a = 0.255 m [10], and the J-TEXT plasma, with a circular cross section in typical discharges, is performed in the L-mode regime with ohmic heating and 500 kW ECRH.

The database employs 236 available discharges during #1061280 to #1061750 and those discharges with no necessary diagnostic signals and failed plasma currents are eliminated. Three types of signals are selected as sample data, i.e. core soft x-ray emission intensity (I_sxr) from the soft x-ray diagnostic (SXR) system [11], edge toroidal rotation velocity (V_φ) from the edge rotation measurement (ERD) [12], and the magnetic signal amplitude (b_θ) from the Mirnov probe (MP) [13], remarkably reflecting the features of diverse plasma events. The raw signals of SXR and MP have the sampling rate of 500 kHz while the signal of ERD has that of 10−50 Hz. The magnetic signal b is calculated by the amplitude of raw magnetic signal x in time interval of about 1 ms,

where l equals 250. l should be greater than the ratio of RTM time scale (usually greater than 1.5 kHz in J-TEXT) to the sampling time interval while overlarge l can reduce the temporal resolution of classification. Therefore, l should be within the region of 175 to 250 in J-TEXT. Soft x-ray emission intensity and the magnetic signal are down sampled with 1 kHz while edge toroidal rotation velocity of 1 kHz is derived by the linear interpolation.

In comparison with applying Z-score and Sigmoid-function normalization methods, with employing min-max normalization the results of cluster analysis have the best accuracy and then this method is employed to uniform the dataset,

where $X_{i}^{*}$ is the diagnostic signal of one kHz and X_i denotes normalized signals I_sxr, V_φ or b_θ, constituting the sample collection X. The database is randomly divided into two parts: 80% (189 discharges) as the train set and 20% (47 discharges) as the test set. Table 1 lists variate descriptions and their parameter ranges in the database. The line-averaged electron density n_e is measured by far-infrared HCN interferometer [14]. The database contains very broad discharge parameters on J-TEXT, covering experimental research on plasma disruption, modulation of electron density, MHD research on applying resonant magnetic perturbations (RMPs) and electron cyclotron resonance heating (ECRH), and so on. In addition, to assess the results of cluster analysis, the total data is marked manually: identifying ST (RTM) by sawtooth (quasi-sinusoidal) oscillations in diagnostic signals, and LM by diverse effects such as increasing V_φ or radial magnetic field b_r, decreasing I_sxr and even LM unlocking phenomena.

Table  1.  Plasma parameter range of sample data.

Signal description	Minimum	Maximum	Variable name
Plasma current	115	200	Ip (kA)
Edge safety factor	2.6	5.4	qa (a.u.)
Soft x-ray emission intensity	0.06	1.95	Isxr (a.u.)
Edge toroidal rotation velocity	−25	35	Vφ (km s−1)
Mirnov signal amplitude	0.04	2.40	bθ (a.u.)
Line-averaged electron density	0.6	4.84	ne (1019 m−3)

 | Show Table

DownLoad: CSV

Figure 1 shows a typical discharge applying the RMP system [15] in J-TEXT. In this discharge with plasma current I_p of 175 kA and toroidal magnetic field B_t of 1.74 T, RTM emerged spontaneously with b_θ of obvious oscillations in figure 1(d). With applying the RMP system, RTM was locked at 0.3195 s, accompanied with increasing V_φ in figure 1(e). Afterwards, LM was maintained until 0.3824 s and then started to rotate. The re-rotating mode slowly reduced and vanished at 0.402 s. Subsequently, ST appeared in figure 1(c) until the discharge end.

Figure  1.  A typical discharge applying RMP: (a) plasma current, (b) RMP current, (c) soft x-ray emission intensity, (d) the Mirnov probe signal and (e) edge toroidal rotation velocity. Rotating tearing mode (RTM), sawtooth oscillations (ST) and locked mode (LM) emerged in this discharge.

DownLoad: Full-Size Img PowerPoint

3.   Survey of plasma events by k-means algorithm

In the k-means algorithm, the cluster number k needs to be set initially. There are several ways to take an appropriate k value. In this work, the 'elbow' method of the k-SSE curve is employed. Figure 2 shows k-SSE curve when scanning k from 2 to 9. As it is known, larger cluster numbers can reduce SSE and hence SSE is decreasing with increasing k. However, too large k may not conform to the actual situation. The appropriate k would be at the 'elbow' point. In figure 2, SSE has sharp decrease when k = 3 while those of k > 3 descend slowly. Therefore, k is set to 3 in this work.

Figure  2.  SSE of cluster analysis with k = 2−9.

DownLoad: Full-Size Img PowerPoint

The k-means clustering analysis output the event classification and cluster centers on table 2 with an accuracy rate of 98.09% in comparison with manually labelled categories. The confusion matrix between the manually labeled categories and train set categories is listed on table 3. 427 samples of LM are labelled as ST while 242 and 148 samples of ST labelled as RTM and LM, respectively. The failed marked samples occur mainly during an event transition process, like excitation of locked mode by RMP. The error is likely due to temporal resolution (the time-integrated signal, with time sampling interval of 10−100 ms) of ERD and the judgment method of the event boundaries between cluster algorithm and physical understanding: the cluster analysis prefers to divide plasma events by their middle distance while we can define RTM (LM) strictly by its amplitude (rotation velocity). When two clusters have different sizes, their boundary has unequal distances away from cluster centers, leading to discrepancy between the realistic and middle-line decision boundaries. In the k-means algorithm, the boundary problem is inevitable due to the decision method (the nearest distance) although the fuzzy k-means algorithm [16] maybe deal with the sample data near the boundary by distribution probability.

Table  2.  The centers list of normalized sample data.

Cluster	Isxr	Vφ	bθ	Label
1	0.1427	0.3419	0.5824	RTM
2	0.2583	0.1905	0.0671	ST
3	0.0760	0.6649	0.0555	LM

 | Show Table

DownLoad: CSV

Table  3.  Confusion matrix between the manually labeled and train set categories.

Cluster manual	RTM	ST	LM
RTM	12973	0	50
ST	242	27849	148
LM	54	427	6450

 | Show Table

DownLoad: CSV

The clustering result is verified in test set. By comparing distances between sample x_i and three centers on table 2 based on equation (3), data in test set is identified and labelled as RTM, ST and LM with the accuracy rate of 98.05%. The results of identification are shown in figure 3. The confusion matrix between the manually labeled categories and predicted categories is listed on table 4. 172 samples of LM are falsely labelled as ST although RTM is identified accurately. The error is attributed to the above-mentioned boundary definition and temporal resolution of diagnosis.

Figure  3.  (a) Classification of samples in test set, and (b) the responding I_sxr, V_φ and b_θ in test set. The magenta '+' sign denotes test error.

DownLoad: Full-Size Img PowerPoint

Table  4.  Confusion matrix between the manually labeled and predicted categories.

Test manual	RTM	ST	LM
RTM	2491	0	6
ST	0	7553	41
LM	0	172	1396

 | Show Table

DownLoad: CSV

The boundary of two clusters can be determined by solving the plane where each point has the same distance from their centers, and hence the boundary is actually on the middle plane D_{i, j} between their centers c_i and c_j. The boundary of the ith cluster is made of the middle planes D_{i, j} and D_{i, k}. D_{i, j} can be derived by solving the equation,

where d denotes any point $\left\{d_{1}, d_{2}, d_{3}\right\}$ in plane D_{i, j}, and c_i (c_j) denotes event center $\left\{c_{i}^{1}, c_{i}^{2}, c_{i}^{3}\right\}\left(\left\{c_{i}^{1}, c_{j}^{2}, c_{j}^{3}\right\}\right)$. We substitute d, c_i and c_j back into equation (7), and then obtain the plane D_{i, j},

According to the centers on table 2, the functions of middle planes are

It should be noted that a portion of D_{1, 2} extends into the third cluster, and hence this part should be excluded by computing their distances from the centers. By the above method, the realistic boundaries (the green faces) are calculated and shown in figure 4 depicting plasma events distributed in train set spaces vividly. It is found in figure 4 that RTM occupies the space with larger b_θ while ST and LM have a smaller b_θ. Intuitively, there should be one simple method to identify RTM events: comparing critical value, instead of computing the distances from three centers. Here, the critical value c₁ can be estimated by the b_θ value $\left(C_{i}^{3}\right)$ of three centers.

Figure  4.  The respective spatial regimes of diverse plasma events. Magenta dots locate at the centers. Green faces denote their boundaries.

DownLoad: Full-Size Img PowerPoint

Clearly, one simple formula can be utilized to recognize RTM,

ST and LM are separated by plane D_{2, 3}, and hence equation (10) can be available to distinguish them. In equation (10), the third term on the left side can be ignored due to its coefficient being much less than the others. It is clear that equation (10) can be simplified as,

By equation (14), we can approximately distinguish ST and RTM, and hence we can obtain,

where c₂ is $0.3842 I_{\mathrm{sxr}}+0.362$ Formulas (13), (15) and (16) assemble an analytical model of identifying plasma events rapidly, as shown in figures 5 and 6. The criterion is consistent with the physical understanding: b_θ will have obvious oscillations of a large amplitude when TM exists, while edge V_φ would change in the direction of plasma current when LM was triggered. Besides an intuitive explanation of the analytical model for its classification, the analytical model can be still utilized to identify the plasma event in specific scenarios of missing certain signals: RTM only by equation (13) and b_θ; LM only by V_φ > c₃, and ST only by I_sxr > c₄. c₃ = 0.5138 and c₄ = 0.4014 can be derived by substituting the maximum I_sxr (0.3951) and maximum V_φ (0.5162) in clusters of ST and LM into equations (15) and (16), respectively. Hence, the case of V_φ > 0.5138 or I_sxr > 0.4014 can be still accurately identified as LM or ST, though missing I_sxr or V_φ, which can be also roughly verified in figure 5.

Figure  5.  Plasma events distributed in (b_θ, V_φ) and (b_θ, I_sxr) spaces.

DownLoad: Full-Size Img PowerPoint

Figure  6.  Plasma events distributed in (V_φ, I_sxr) space. LM and ST events were divided by the dashed magenta line.

DownLoad: Full-Size Img PowerPoint

The analytical model is extremely close to the clustering method with an accuracy rate of 99.78%. Therefore, the model can be used to label plasma event in the actual discharge #1062247. In this discharge with ramping I_p from 160 to 200 kA and q_a from 3.0 to 2.6, RMP is applied during the period of 0.16–0.426 s, and three events exist successively. It should be noted that V_φ calculated during plasma current ramp would be in great errors and hence its acquisition time is from 0.1 s to ~ 0.55 s in general, corresponding to the plasma flattop phase.

According to the analytical model (c₁ and c₂) and normalization coefficients in table 1, we can compute the reference values for the discharge parameter, C_RTM, and C_LM, to identify plasma events,

Then, we can obtain C_RTM = 0.799 and C_LM = 12.23I_sxr − 5.239. The rotating tearing mode has a larger amplitude than C_RTM and then is identified as RTM (the dark line in figure 7(d)). Afterward, the tearing mode is locked at 0.29 s by RMP, accompanied by increasing V_φ larger than C_LM (the blue line in figure 7(e)), leading to the plasma event labelled as LM. After RMP is off, LM reduces its amplitude and then disappears. When V_φ decreases, sawtooth oscillations emerge (the red line in figure 7(e)) and are labelled as ST. In this discharge, plasma events are identified with the accuracy rate of 95%, and the error comes from the event transition process, probably attributed to boundary problems and sampling time intervals (20 ms) of ERD, as previously mentioned.

Figure  7.  Recognition of plasma events in the fresh discharge #1062247: (a) plasma current, (b) RMP current, (c) soft x-ray emission intensity, (d) the Mirnov signal amplitude and (e) edge toroidal rotation velocity (dark squares). The cyan curves represent C_RTM = 0.799 and C_LM = 12.23I_sxr − 5.239.

DownLoad: Full-Size Img PowerPoint

It should be noted that attempts of adding more variables had been performed, and the results of clustering analysis revealed that V_φ, b_θ and I_sxr (n_e) played a dominant impact, probably owing to them including event characteristics. When n_e was absorbed into sample data, the results of clustering analysis were similar to the case without it, possibly due to n_e having a similar evolution to I_sxr [17]. Besides, too superfluous variables would not increase the accuracy, probably due to superfluous variates dimming the features of meaningful ones. Therefore, it is better to choose variables which are capable of significantly representing event characteristics.

4.   Summary

In this work, the k-means clustering algorithm is utilized to classify plasma events in the J-TEXT plasma, like the rotating tearing mode, sawtooth oscillations, and locked mode, and then one simple analytical model is proposed to identify these plasma events with a high accuracy rate by basic diagnosis signals having obvious features of plasma events. The analytical model can be applied to data markers of massive diagnostic data with plasma events, especially for research on controlling the tearing mode by RMP and RMP penetration in J-TEXT.

In this algorithm, Euclidean distance is employed to measure the difference between plasma events, and the distance to centers of plasma events decides the type of events. Therefore, it is difficult for this method to exactly calculate the event borders or event transition processing. At present, this model can be appropriate only for the recognition of independent plasma events but struggles to distinguish the coexisting MHD activities, like sawtooth and edge located mode in H-mode plasma. In J-TEXT, large rotating and locked mode coexists rarely with sawtooth oscillations [18], except for precursor oscillations of sawteeth. In the future, much more research on k-means and even other cluster analysis algorithms may solve these two problems.