Plasma current tomography for HL-2A based on Bayesian inference

1.   Introduction

In a magnetic confinement fusion device, the temperature of the plasma can exceed 100 million degrees Celsius. Under such a high temperature, it is difficult to directly obtain control parameters of the device, such as the plasma shape, magnetic axis, X-point, etc. In practice, plasma parameters are indirectly reconstructed through diagnostics measurements to provide control information for plasma. On HL-2A, the plasma parameters are provided from Equilibrium FITting (EFIT) [1]. EFIT assumes that the plasma current satisfies the magnetohydrodynamic equilibrium, electromagnetic measurement is used as the boundary condition, and the optimal solution of the plasma current distribution is obtained based on the Picard linearization scheme. EFIT is based on specific physical models, and it is highly robust and has been widely used in various devices, such as JET [2], DIII-D [3], NSTX [4], EAST [5], KSTAR [6], etc.

With the popularity of machine learning, there are more and more successful cases based on machine learning in the field of fusion, which provides new ideas for the study of fusion. The most famous cases include plasma disruption prediction based on deep learning [7] and plasma profile control based on reinforcement learning [8]. Machine learning can mainly be divided into two categories: traditional machine learning and deep learning. Deep learning is generally based on big data, and it uses neural network algorithms for training and learning. Data-driven deep learning algorithms are very suitable for handling complex non-linear problems in fusion. Migration between different devices is difficult because data from the new device are still required for retraining [9]. Bayesian inference is a classic algorithm in machine learning, which is mainly used to solve the inverse problem. Unlike deep learning, it does not require large amounts of data for training. In the Bayesian framework, each piece of information will be given a possible probability density function. When a new piece of information is added, a new joint probability density function will be obtained through Bayes’ theorem. Since the algorithm includes physical processes, such as the relationship between electromagnetic measurements and plasma current in plasma current inversion, it is easy to migrate to different devices. In existing studies, plasma profile inversion for multiple diagnostics using Bayesian inference has been investigated, such as the radiation emissivity profile [10, 11], carbon impurity radiation emissivity profile [12], electron temperature profile [13], electron density profile [14], plasma current profile [15–18], etc.

As shown in equation (1), the Bayesian inference algorithm is based on the Bayesian formula. The formula contains four items: prior probability $P(f)$, likelihood probability $P(d|f)$, posterior probability $P(f|d)$, and the evidence item $P(d)$. In essence, the prior probability $P(f)$ is the probability assigned to the unknown parameters $f$ before considering any observed data. It represents our belief in the unknown parameters $f$ in the absence of new information. The likelihood probability $P(d|f)$ is the probability of observing data $d$ given that $f$ is known, which describes how well $f$ explains or predicts the observed data $d$. The evidence term $P(d)=\int P(d|f) P(f) {\rm{d}} f$ represents the normalization constant, which does not impact the peak value of the posterior probability and is often disregarded. The posterior probability $P(f|d)$ is the probability of $f$ after taking into account the observed data $d$. It is the main goal of the Bayesian inference and represents our belief in $f$ after considering new data $d$. When it reaches its maximum, it corresponds to the optimal value of the unknown parameter $f$.

In previous work, the plasma current density profile was reconstructed for EAST based on Bayesian inference [17, 18]. In this work, the algorithm will be applied to the HL-2A device to invert the plasma current density. In addition, two different prior probabilities, the Conditional AutoRegressive (CAR) prior and the Advanced Squared Exponential (ASE) kernel prior will be used in this work. The rest of the paper is organized as follows: in section 2, the electromagnetic measurement diagnostics on HL-2A are introduced; in section 3, the tomographic method is developed and verified; in section 4, the tomography method is applied to the analysis of experimental data; in section 5, this paper is concluded.

2.   Electromagnetic measurements on HL-2A

HL-2A is a tokamak nuclear fusion experimental device with an advanced divertor configuration and non-circular cross-section. The major radius is 1.65 m, and the minor radius is 0.4 m. Its total plasma current can reach 480 kA and the discharge time is about 5 s. Its goal is to conduct improved confinement experiments under high-parameter plasma conditions and utilize its unique large-volume closed divertor to investigate various cutting-edge physics and engineering topics in the field of nuclear fusion. This research aims to lay a technological foundation for the study of next-generation fusion reactors in China. The HL-2A is equipped with six sets of coils, which are used for plasma generation, heating, and confinement. As shown in formulas (2), (3), and (4), the model incorporates diagnostics related to plasma current, including pickup coils, magnetic flux loops, and Rogowski coils [19].

where $B_{R}$ and $B_{Z}$ are the $R$ and $Z$ directions of the magnetic field at the position with angle $\theta$. $\psi$ is the magnetic flux at $(R,\ Z)$. $I_{R, Z}$ represents the plasma current at $(R,\ Z)$. Pickup coils are employed to measure local magnetic fields, while magnetic flux loops are used to measure the magnetic flux. As we know, the magnetic field and flux are both functions of the plasma current. Additionally, the plasma total current is measured by Rogowski coils. In the HL-2A device, there are 18 pickup coils, four magnetic flux loops, and one Rogowski coil, and they are distributed around the vacuum chamber. Based on these 23 diagnostic signals, the internal plasma current distribution is reconstructed. Previous studies have demonstrated that using only external magnetic diagnostics cannot accurately reconstruct the internal plasma current distribution [20], and this is why this study focuses mainly on the reconstruction of the plasma boundary.

The high-temperature plasma is confined within the first wall, also called the vacuum chamber, so the region for model inversion is established within the first wall. It is assumed that plasma is symmetric in the toroidal direction, and the three-dimensional (3D) problem can be transformed into the two-dimensional (2D) problem. Then, plasma region is divided by rectangular grid beams, and the current in any grid is assumed to be uniform, as illustrated in figure 1.

Figure  1.  The plasma area within the first wall is marked with green pixels. Each pixel is scaled by 0.93 cm × 1.49 cm and carries a uniform current. There are 4205 pixels in total.

DownLoad: Full-Size Img PowerPoint

3.   Bayesian model

3.1   Likelihood probability

As mentioned earlier, likelihood probability links unknown parameters and the observations. In this study, the observations are electromagnetic measurements, and the unknown parameters are the plasma currents at grid points. The electromagnetic measurements include pickup coils, flux loops, and Rogowski loops, which can be represented by the vector $\bar{D}^{\rm{Mag}}$. The plasma currents at grid points can be represented by the vector $\bar{I}$. According to the Biot–Savart law, there is a linear relationship between the diagnostics and plasma currents, as shown in equation (5):

$\bar{\bar{R}}$ is the response matrix. $\bar{C}$ is an integral part of the electromagnetic measurement signal, which is determined by the six sets of coil current. Since the coil current can be directly measured, $\bar{C}$ can be directly calculated as a constant. $\bar{\varepsilon}^{\rm{Mag}}$ is the error of electromagnetic measurements, and it satisfies the Gaussian distribution. Assigning probabilities to the distribution of errors, the likelihood probability based on electromagnetic measurements can be expressed in equation (6):

3.2   Gaussian prior probability

3.2.1   CAR prior.

The CAR prior associates the current of any grid point with the current of its four adjacent grid points [17]. It is assumed that the current at the center position is equal to the average of the currents at the four adjacent grid points. To capture the relationship of currents on grid points, the covariance matrix $\bar{\bar{Q}}$ of the Gaussian distribution is constructed:

where $\bar{\bar{1}}$ is the identity matrix and $\bar{\bar{W}}$ is the adjacency matrix. If $n$ and $m$ are adjacent, then $\bar{\bar{W}}^{n, m} = 1$; otherwise, $\bar{\bar{W}}^{n, m} = 0$. It is noteworthy that for different grid points, the values of $\tau$ are the same. Therefore, $\tau$ is a stationary hyperparameter, and it can be selected through Bayesian Occam’s razor [17]. To maintain the unbiasedness and translation invariance of the Gaussian prior probability, the mean of the Gaussian prior is set to be $\bar{0}$. Thus, the prior probability density can be represented as follows:

where $N_I$ denotes the number of grid points.

3.2.2   ASE kernel prior.

The ASE kernel prior is developed based on the squared exponential (SE) kernel prior [18]. The difference is that the reference discharge information is added. Compared with the CAR prior, the SE kernel function is more flexible, and does not need to satisfy the consistency condition [17]. The SE kernel prior can integrate more information by adjusting the hyperparameters. The introduction of reference discharge can ensure that the distribution form of plasma current is similar, thus eliminating more impossible values. The SE kernel function determines the covariance matrix of the prior probability, as shown in equation (9):

where ${\bar{x}^i}$ and ${\bar{x}^j}$ represent the different locations of the grid points. $\sigma$ is called the standard deviation, which determines the size of the average distance of the function away from its mean. $\ell$ is the characteristic length. $\sigma$ and $\ell$ are also hyperparameters. For the ASE prior, it is assumed that $\sigma$ is proportional to the plasma current of the reference discharge, as shown in equation (10). For the HL-2A device, once the operation scenarios are determined, the characteristic performance becomes quite predictable, and the difference between similar scenarios should be relatively small.

where $\bar{I}'_i$ and $\bar{I}'_j$ denote the currents on the grids at different locations. $k$ is a constant to be optimized by Bayesian Occam’s razor. Since the currents on different grid points are generally different, the current at the core position is large, and the current at the boundary is relatively small. Based on this, the spatial inhomogeneity can be achieved by only adjusting $k$. Therefore, the hyperparameters of the ASE prior are non-stationary. With this kernel function, the covariance matrix $\bar{\bar{\Sigma}}_{I}$ is constructed as follows:

For the mean of the Gaussian prior, it is also set to $\bar{0}$ like the CAR prior. In this way, the prior probability of the plasma current density has been constructed:

3.3   Posterior probability

In the previous two sections, the likelihood probability and the prior probability have been constructed. Then, according to the Bayesian formula, the posterior probability can be obtained, as shown in equation (13):

Bringing the Gaussian prior probability and likelihood probability into equation (13), it can be found that the posterior probability is still subject to Gaussian distribution.

For the CAR prior, where

For the ASE prior, where

$\bar{m}$ denotes the mean of the new Gaussian distribution after coupling the prior probability and the likelihood probability. It is well known that when the function value is equal to its mean value, the Gaussian distribution achieves the maximum probability. Therefore, the optimal current profile $\bar{I}$ is equal to $\bar{m}$. $\bar{\bar{\Sigma}}$ denotes the covariance of the posterior probability, and it can quantify the uncertainty of the entire inversion process.

3.4   Synthetic test

The known discharge data are used to evaluate the algorithm, and the process is called the synthetic test. The known discharge data are designed by the EFIT code in a semi-fixed boundary calculation model [21]. EFIT has been used to reconstruct plasma equilibrium and for plasma discharge control for many years, so the designed discharge is representative. In this mode, EFIT is used to calculate a desired plasma equilibrium with given plasma boundary. At the same time, the required diagnostic signals are also obtained. In addition, the data include six sets of coil currents, plasma current distribution, and magnetic flux distribution. Among them, the plasma current distribution, magnetic flux distribution, and plasma boundary are shown in figure 2.

Figure  2.  (a) The plasma current contour map of the design discharge. (b) The magnetic flux contour map of the design discharge, where the red solid line represents the plasma boundary.

DownLoad: Full-Size Img PowerPoint

3.4.1   CAR model.

After inputting the electromagnetic measurement signals into the current tomograph based on the CAR prior, the reconstructed results are shown in figure 3. In figures 3(b) and (e), $\xi^I$ and $\xi^{\psi}$ represent the relative errors of the plasma current and magnetic flux, respectively. Their defining formulas are expressed as $\xi_{i}^{I}=\dfrac{\left|\bar{I}_{i}^{\rm{rec}}-\bar{I}_{i}\right|}{\max \{\bar{I}\}}$ and $\xi_{i}^{\psi}=\dfrac{\left|\bar{\psi}_{i}^{\rm{rec}}-\bar{\psi}_{i}\right|}{\max \{\bar{\psi}\}}$, where $\bar{I}^{\rm{rec}}$ and $\bar{\psi}^{\rm{rec}}$ represent the recon-structed plasma current and magnetic flux, respectively. From figures 3(b) and (e), it can be found that the maximum relative errors of the plasma current and magnetic flux are about 12% and 1.2%, respectively. From figure 3(a), it can be found that the plasma current profile agrees with our experience that the core is high and the edge is low. The uncertainty in figure 3(c) shows that there is the least information at the core because the diagnostics are far away from the core but close to the plasma boundary. Figure 3(f) shows that the reconstructed boundary is consistent with the designed boundary. The error of the main control point (C1–C6 and X) is kept within 1 cm. The numerical values of the errors are listed in table 1.

Figure  3.  (a) The contour map of the reconstructed plasma current based on Bayesian inference. (b) The relative error distribution map between the reconstructed plasma current based on Bayesian inference and the designed discharge, where the black point represents the position of the maximum relative error. (c) The plasma current along the radius, where the blue solid line shows the current obtained from Bayesian reconstruction, and the red dotted line shows the designed current. The gradually lighter gray level shows the confidence interval (also called uncertainty) of the Gaussian distribution $[-\sigma,\ \sigma]$, $[-2\sigma,\ 2\sigma]$, and $[-3\sigma,\ 3\sigma]$, and their probabilities are 68.26%, 95.44%, and 99.74%, respectively. (d) The magnetic flux contour map based on the reconstructed current, where the red solid line shows the plasma boundary. (e) The relative error of the magnetic flux, where the black dot shows the position of the maximum error. (f) The reconstructed plasma boundary (blue solid line) and the design discharge boundary (red dotted line). C1–C6 represent the boundary errors from directions of 0, $0.25 {\text{π}}$, $0.75\mathrm{ {\text{π}} }$, ${\text{π}}$, $1.25 {\text{π}}$, and $1.75 {\text{π}}$, respectively. X represents the point X of the divertor configuration, and the black point represents the position of the maximum error (Maximum error). The error values of the boundary are shown in table 1.

DownLoad: Full-Size Img PowerPoint

Table  1.  The plasma boundary error from the CAR prior (cm).

${\rm{C1}}$	${\rm{C2}}$	${\rm{C3}}$	${\rm{C4}}$	${\rm{C5}}$	${\rm{C6}}$	${\rm{X}}$	Maximum error
0.0118	0.276	0.292	0.148	0.0596	0.0943	0.0209	0.521

 | Show Table

DownLoad: CSV

Table 1 shows that the maximum error of the boundary is less than 1 cm. As mentioned in section 3.1, the size of the grid is cm $0.93 \times 1.49$ cm, and the maximum length in the grid is 1.76 cm (diagonal length $\sqrt{0.93^{2}+1.49^{2}}=1.76$). The boundary error caused by interpolation should be dominant, so the error is below it, which is acceptable.

To simulate the real discharge environment, noises are added to the input electromagnetic measurement signals to test the model’s robustness. The electromagnetic measurements are accurate diagnostics, and the maximum error is no more than 3%. According to the characteristic of Gaussian distribution, when the random variable falls within $[-3\sigma,\ 3\sigma]$, the probability is 99.74%. Here, $\sigma$ is set to 1%, and 100 sets of diagnostic data with 3% random noise are obtained. The errors of boundary reconstruction based on these diagnostics are presented in figure 4. The boundary errors of the seven positions considered in this study are all controlled at about 1 cm and below the threshold of 1.76 cm, which meets the requirements of control.

Figure  4.  The plasma boundary errors based on the CAR model.

DownLoad: Full-Size Img PowerPoint

3.4.2   ASE model.

Figure 5 shows the data of reference discharge, and the ASE kernel prior described in section 3.2.2 is utilized to reconstruct the plasma current distribution. It has been demonstrated in a previous study that the ASE kernel prior is tolerant to the reference discharge, even if there are certain differences between the reference discharge and true discharge [18]. So, the reference discharge is randomly selected here, but it is necessary to ensure that plasma current distribution of the discharge is similar. It can be seen from the figure that there is a difference between the reference discharge current and the discharge current that needs to be reconstructed, and the largest difference is at the position of the core. Meanwhile, their boundaries are slightly different. Based on this reference discharge, the reconstruction results are shown in figure 6. As shown in figure 6(e), the largest current error appears at the core, and figure 6(c) shows that the uncertainty is the largest at this position. This is consistent with the location of the largest difference between the reference discharge and the design discharge. Figures 6(b) and (e) show that the maximum relative errors of plasma current and magnetic flux are only 6% and 0.35%, respectively. The maximum error of plasma current is near the boundary. The main reason is that the difference between the reference discharge and the real discharge is the largest at the boundary, as shown in figure 7. Compared with the results of CAR reconstruction, the reconstruction accuracy under the ASE model is higher. As mentioned in section 2, using only external magnetic diagnostics cannot accurately reconstruct the internal plasma current distribution. The CAR prior only smoothes the distribution of the plasma current, while the ASE kernel prior couples more plasma current information, extracts information about current distribution, and excludes more impossible results. Thus, the reconstruction results from the ASE kernel prior are better. From figure 6(f), the reconstructed plasma boundary agrees well with that of the designed discharge, and the boundary errors are shown in table 2. The errors of the main control points (C1–C6, X) and the global maximum boundary error are all less than 1 cm.

Figure  5.  The left panel is the radial distribution of the plasma current, and the right panel is the plasma boundary of the reference shot and the test synthetic. The red dashed line represents the synthetic test, and the blue line represents the reference discharge.

DownLoad: Full-Size Img PowerPoint

Figure  6.  (a) The plasma current contour map, (b) the relative error distribution map, (c) the plasma current along the radius, (d) the magnetic flux contour map, (e) the relative error of the magnetic flux, (f) the plasma boundary.

DownLoad: Full-Size Img PowerPoint

Figure  7.  Relative error between reference discharge and true discharge; the value of the white part exceeds 100%.

DownLoad: Full-Size Img PowerPoint

Table  2.  The plasma boundary error from ASE (cm).

${\rm{C1}}$	${\rm{C2}}$	${\rm{C3}}$	${\rm{C4}}$	${\rm{C5}}$	${\rm{C6}}$	${\rm{X}}$	Maximum error
0.0497	0.0782	0.120	0.0996	0.0734	0.0495	0.0270	0.450

 | Show Table

DownLoad: CSV

To test the robustness of the ASE prior, the diagnostics with 3% random noise are input to the model for profile inversion. The reconstruction results from 100 sets of diagnostics are counted and presented in figure 8. The boundary errors of the seven positions are all below 1 cm.

Figure  8.  The plasma boundary errors based on the ASE model.

DownLoad: Full-Size Img PowerPoint

3.4.3   Comparative analysis of two models.

As introduced in the previous two sections, the plasma current tomography of the two prior models can reconstruct the plasma current, and the calculation results under the two models are listed in table 3. The Max $\xi^I$ and Max $\xi^{\psi}$ in table 3 represent the maximum relative errors of the plasma current and flux, respectively. The ${\rm{RMSD}}^{I}$ and ${\rm{RMSD}}^{\psi}$ are respectively the root mean square deviation (RMSD) of the plasma current and flux, and they are expressed as follows: ${\rm{RMSD}}^{I}= \sqrt{\dfrac{\sum_{i=1}^{n}(\bar{I}_{i}^{\rm{rec}}-\bar{I}_{i})^{2}}{n}}$ and ${\rm{RMSD}}^{\psi}=\sqrt{\dfrac{\sum_{i=1}^{n}\left(\bar{\psi}_{i}^{\rm{rec}}- \bar{\psi}_{i}\right)^{2}}{n}}$. It can be seen from table 3 that the ASE prior coupling reference discharge contains more information on the plasma current, so both the relative error and the RMSD are lower than those of the CAR prior. Generally, under the conditions of electromagnetic measurement, the Bayesian model based on the ASE prior is more advantageous.

Table  3.  Results from different models.

${\rm{Model}}$	${\rm{Max}}$ ${\xi^{I}}$ ${(\%)}$	${{\rm{RMSD}}^{I}}$ ${({\rm{A}})}$	${\rm{Max}}$ ${{\xi}^{\psi}}$ ${(\%)}$	${{\rm{RMSD}}^{\psi}}$ ${({\rm{Wb}})}$
${\rm{CAR}}$	12.76%	6.0305	1.21%	0.0008
${\rm{ASE}}$	5.77%	2.0699	0.36%	0.0004

 | Show Table

DownLoad: CSV

4.   Experimental validation

The CAR and ASE prior plasma current tomography algorithms are applied to the experimental data. The shot number is 35000, and the time is 500 ms. The reconstruction results are illustrated in figure 9. The boundaries (red lines) reconstructed by the plasma tomography algorithms based on the two priors are almost consistent with that (the white line) reconstructed by EFIT, and the boundary differences are shown in table 4. It can be seen from the table that the maximum error between the boundary of the CAR model and the boundary of EFIT reconstruction is 1.6 cm, while the maximum error between the boundary of the CAR model and the boundary of EFIT reconstruction is only 0.8 cm. From figure 9(c), it can be found that the uncertainty under the CAR model is higher and remains at the same level. This is because the CAR model uses stationary parameters, and $\tau$ adopts the same value. In the inversion process, to match all the grids, the hyperparameters are set larger. But the ASE model couples the reference discharge to achieve non-stationary hyperparameters, so the uncertainty is lower.

Figure  9.  (a)–(c) Results from CAR model, (d)–(f) results from ASE model. The white line in each figure represents the boundary reconstructed by EFIT, and the red line represents the boundary reconstructed by plasma current tomography. (a) and (d) The magnetic flux distribution, where the black solid line represents the contour line. (b) and (e) The distribution of the plasma current, where the black solid line represents the contour line. (c) and (f) The standard deviation $\sigma$ distribution of the plasma current.

DownLoad: Full-Size Img PowerPoint

Table  4.  The plasma boundary difference compared to EFIT (cm).

${\rm{Model}}$	${\rm{C1}}$	${\rm{C2}}$	${\rm{C3}}$	${\rm{C4}}$	${\rm{C5}}$	${\rm{C6}}$	${\rm{X}}$	Maximum error
CAR	0.930	1.37	0.138	0.920	0.460	0.590	0.0146	1.60
ASE	0.280	0.450	0.700	0.700	0.420	0.050	0.0894	0.800

 | Show Table

DownLoad: CSV

5.   Conclusion

In this study, the plasma current tomography methods based on two priors (the CAR prior and the ASE kernel prior) are applied to the HL-2A device. The 2D plasma current profile is successfully reconstructed, and then the magnetic flux distribution, plasma boundary, etc. are identified. Different from the CAR prior, the ASE prior couples more information, adopts non-stationary hyperparameters, and achieves better reconstruction accuracy. Compared with the existing profile algorithm EFIT on the HL-2A device, the feasibility and accuracy of the proposed algorithm are demonstrated. This algorithm can not only invert the plasma parameters but also give the uncertainty in the inversion process. This study demonstrates that current tomography based on Bayesian inference can be applied to different devices and provide technical support for the equilibrium reconstruction of future fusion devices.

Acknowledgment

The authors are grateful to the HL-2A team for providing the experimental data. This work is supported by the National MCF Energy R&D Program of China (Nos. 2018YFE0301105, 2022YFE03010002 and 2018YFE0302100), the National Key R&D Program of China (Nos. 2022YFE03070004 and 2022YFE03070000) and National Natural Science Foundation of China (Nos. 12205195, 12075155 and 11975277)