New Identiﬁcation Approach and Methods for Plasma Equilibrium Reconstruction in D-Shaped Tokamaks

: The paper deals with the identiﬁcation of plasma equilibrium reconstruction in D-shaped tokamaks on the base of plasma external magnetic measurements. The methods of such identiﬁcation are directed to increase their speed of response when plasma discharges are relatively short, like in the spherical Globus-M2 tokamak (Ioffe Inst., St. Petersburg, Russia). The new approach is ﬁrst to apply to the plasma discharges data the off-line equilibrium reconstruction algorithm based on the Picard iterations, and obtain the gaps between the plasma boundary and the ﬁrst wall, and the second is to apply new identiﬁcation methods to the gap values, producing plasma shape models operating in real time. The inputs for on-line robust identiﬁcation algorithms are the measurements of magnetic ﬂuxes on magnetic loops, plasma current, and currents in the poloidal ﬁeld coils measured by the Rogowski loops. The novel on-line high-performance identiﬁcation algorithms are designed on the base of (i) full-order observer synthesized by linear matrix inequality (LMI) methodology, (ii) static matrix obtained by the least square technique, and (iii) deep neural network. The robust observer is constructed on the base of the LPV plant models which have the novelty that the state vector contains the gaps which are estimated by the observer, using input and output signals. The results of the simulation of the identiﬁcation systems on the base of experimental data of the Globus-M2 tokamak are presented.


Introduction
Tokamaks [1], toroidal vessels with magnetic coils (Figure 1), originated at the I.V. Kurchatov Institute of Atomic Energy in the USSR and spread around the world to solve the problem of controlled thermonuclear fusion: obtaining energy from the fusion of the light elements nuclei.The most promising devices for solving this problem are vertically elongated tokamaks with increased gas-kinetic pressure (D-shaped tokamaks) (Figure 1).Plasma (the fourth state of matter) vertically elongated by an external magnetic field is unstable in the vertical direction, and it is necessary to use automatic feedback control systems to keep it near the first tokamak wall.
In this work, the new plasma equilibrium reconstruction algorithms are to be inserted into the plasma position, current, and shape feedback control system of the Globus-M2 tokamak.In Figure 2, one can see the digital model of that system [16] where the plasma equilibrium reconstructed algorithm is to be identified by the new methods proposed in the paper.
< l a t e x i t s h a 1 _ b a s e 6 4 = " f 0 3 z 2 t C 5 i 2 M M u e B 6 K C 4

s C P B A E = " >
R k l 2 k j t w o n j 2 H 4 I P S a y M C J a s v u R f e v 7 + N r 2 Q X R + N M d p F 4 p u y Z 0 A L x J v R o p o h m q 7 8 O l 3 Y p p G I D T l R K m W 5 y Y 6 y I j U j H I Y 2 X 6 q I C F 0 Q H r Q M l S Q C F S Q T V 4 a 4 S O j d H A 3 l q a E x h P 1 9 0 R G I q W G U W g 6 z Y l 9 N e + N x f + 8 V q q 7 5 0 H G R J J q E H S 6 q J t y r G M 8 z g d 3 m A S q + d A Q Q i U z t 2 L a J 5 J Q b V K 0 T Q j e / M u L p F 4 u e a c l 7 6 p c r O B Z H H l 0 g A 7 R M f L Q G a q g S 1 R F N U T R A 3 p C L + j V e r S e r T f r f d q a s 2 Y z + + g P r I 9 v 0 m y c N Q = = < / l a t e x i t >  Z R < l a t e x i t s h a 1 _ b a s e 6

I c
< l a t e x i t s h a 1 _ b a s e 6 4 = " e H L U a N V a w 7 R j N m D z h q O B 7 3 c B o S g = " > A A A C V n i c b V B N T 9 w w E P W m p d D 0 K 4 U j F 6 u 7 l a o e V g m H w h G p l 3 K j i F 1 Q N 1 H k O J P F w r E j e 4 J Y W f l n / S H t t T d E f w P C u 0 Q q H 3 2 S p e f 3 Z j z j V z R S W I z j 3 4 P g 2 f O 1 F + s b L 8 N X r 9 + 8 f R e 9 3 5 x a 3 R o O E 6 6 l N q c F s y C F < l a t e x i t s h a 1 _ b a s e 6 4 = " 8 g e 4 7 J I x q y d M 6 5 u s

Reconstructing Plasma Equilibria from External Magnetic Measurements
A tokamak is an axially symmetrical device, so the tokamak plasma equilibrium is described in the poloidal plane (r, z), typically in terms of the poloidal magnetic flux distribution Ψ(r, z), which is defined as the flux of the magnetic field vector B through a surface S bounded by the line (r = const, z = const): The magnetic field lines, along which plasma particles move, lie on the flux surfaces Ψ(r, z) = const; therefore, the boundary of the magnetically confined plasma can be found as the largest closed flux surface.
The toroidal current density J ϕ in the tokamak is connected with the poloidal flux through the linear second-order partial differential equation [1]: The boundary conditions for the equation are obtained from the definition and the physical meaning of the poloidal flux: When the right-hand side of the Equation ( 2) is known, it can be solved with the standard numerical methods, for example, using the corresponding Green's function G [17]: where K and E are the elliptic integrals of the first and the second kind respectively, and In practice, the plasma current distribution and the induced currents in the conductive Vacuum Vessel (VV) of the tokamak are often not available for real-time reconstruction and must be identified together with the poloidal flux distribution from the external magnetic measurements, which include coil currents I 1 , . . ., I N c , total plasma current I p measured by Rogowski coils and poloidal flux values Ψ at finite number of points (r 1 , z 1 ), . . ., (r N l , z N l ) by magnetic loops outside the plasma.Hence, the plasma equilibrium reconstruction problem is to find plasma area S p , plasma current distribution J p and induced current density J v such that: where σ p and σ j are uncertainties of the plasma current and poloidal flux at j th magnetic loop, Ψ(r, z) is the solution of the Equation ( 2) with boundary conditions (3) and the right-hand side: S k and S v are the area occupied by the k th coil and the VV, respectively, N k is the number of turns of the k th coil.Optionally, coil current measurements may also be considered uncertain and accounted in the functional (4) by terms The functional may also include other measurements that can be expressed in terms of the currents and the magnetic flux.Finally, as the plasma equilibrium reconstruction problem is ill-posed in the sense of Hadamard, the functional may include a regularization term.
To find the plasma shape in the Globus-M2 tokamak (Figure 3), the flux-current distribution identification (FCDI) code was used [14].The FCDI code applies the following expression for the plasma toroidal current density, obtained from the plasma force balance equations [1,14]: where p is plasma pressure and F is poloidal current defined analogous to poloidal flux (1): The Picard iteration method is used to find the poloidal flux distribution.Since F and p depend only on poloidal flux, on each iteration, the FCDI code approximates the plasma current density by polynomials of the poloidal flux from the previous iteration: Similarly, the VV currents are approximated as a linear combination of some basis functions, for example, orthogonal VV current modes [18]: The coefficients of the J p polynomials and the J v basis function regression are found then by minimizing the error functional (4) which can be written in the matrix form: Here, c is the N × 1 column-vector of the coefficients A is the M × N matrix, where M is the number of magnetic measurements used, and b is the N × 1 column-vector.To regularize the problem, the SVD truncation method is used to minimize the quadratic functional [19].After the coefficients are determined, the corresponding poloidal flux distribution is calculated, which is used for the polynomials construction on the next iteration.The iterations are continued until the error χ 2 is sufficiently small or the maximal number of iterations is reached.

Experimental Data
The FCDI code was applied to 50 discharges of the Globus-M2 tokamak.For each discharge, there are magnetic measurements available y (1) , y (2) , . . ., y (50) , which include currents in the 8 control coils (Horizontal Field Coil, Vertical Field Coil, Central Solenoid, Poloidal Field Coil 1, upper and lower sections of the Poloidal Field Coil 2, Poloidal Field Coil 3 and Correcting Coil) (Figure 4), poloidal magnetic flux from 21 loops, vertical dipole magnetic flux (difference between magnetic flux above and below plasma), horizontal dipole magnetic flux (difference between magnetic flux on the left and on the right of the plasma), and quadrupole magnetic flux (expressed as ψ(L 1 ) − ψ(L 2 ) + ψ(L 3 ) − ψ(L 4 ), with location of loops L 1 -L 4 shown in Figure 4) so that y (i) ∈ R 33×s i , i ∈ [1; 50], where each s i = T i /τ, T i is the duration of the discharge, τ is the discretization step.Here, the discretization step is the time step between the reconstructed off-line equilibria.It is constrained only by the discretization time of the experimental measurements.From these data, the FCDI code obtains plasma current distribution and plasma boundary coordinates for the divertor phases of the discharges.The calculated plasma shapes are represented by the positions of 2 strike points (g 1 , g 2 ) on the VV and values of 4 gaps (g 3 -g 6 ) between plasma and VV (Figure 4) g (1) , g (2) , . . ., g (50) ; g (i) ∈ R 6×s i .The strike points are points of intersection of the poloidal flux isoline, which bounds the plasma and the VV.Their coordinates g 1 and g 2 are calculated as the distance from point P 6 in Figure 4 along the VV.The gap g 3 is calculated as the distance between point P 3 on the VV outer wall and the plasma boundary on the horizontal line, g 4 is the distance between P 4 and the plasma boundary on the 45 • line, g 5 is the distance between P 5 and the plasma boundary along the vertical line, and g 6 is the distance between point P 6 on the VV inner wall and the plasma along the horizontal line.The g 1 -g 6 values describe plasma shape in the LSND (lower single null divertor) configuration, typical for the Globus-M2 tokamak.Other configurations may require different sets of descriptors, but the identification methods described below are applicable all the same.

Plasma Model
The plasma dynamics is described by Faraday's law equations: and force balance equation The measured fluxes and the plasma shape are determined by currents in the tokamak: Here, I = [I T c , I T v , I p ] T , Φ, R, and U are respectively the column-vector of currents, column-vector of magnetic flux, diagonal matrix of electrical resistance and column-vector of the voltage applied to the control coils, VV, and plasma F is the force acting on the plasma, g is the column vector of strike points positions on the VV and the gaps between the plasma and VV, Ψ is the column vector of the fluxes measured by the tokamak diagnostics.The plasma mass is neglected.
The magnetic flux vector can be expressed as Φ(J p , I) = M(J p )I, where M is the inductance matrix.The dependence of the inductance matrix M, force F and plasma shape g on plasma current distribution J p is nonlinear but for the small deviations from the reconstructed equilibrium, the linearized model is sufficient.Assuming that plasma can rigidly move in vertical and radial directions, the linearized Equations ( 5)-( 7) take form: where r p is the radius-vector r p = [r p , z p ] T of plasma center of mass, δ denotes deviation from the scenario value.
Introducing state vector x = δI = [δI T c , δI T v , δI p ] T , input vector u = δU and output vector of plasma and coil currents, gaps, and fluxes deviations y = [δI p , δI T c , δΨ T , δg T ] T , the LPV (linear parameter varying) model takes the standard state-space form: The reconstructed plasma current distributions J p are used to calculate series of linear models {A, B, C} nm describing plasma dynamics in each considered discharge.Here, index n denotes time moment t n for which the model is obtained where t 1 , . . ., t N m correspond to the time points of the divertor phase of the m th tokamak discharge with the time step of 1 ms and index m denotes the serial number of discharge.This represents the LPV model ( 8) as an array of LTI (linear time invariant) models.During modeling each discharge, a linear interpolation is performed between time points from t 1 to t N m .
The models have 24 states, 8 inputs, and 39 outputs.Each obtained model has a single real positive pole.
Although the models include expressions for the gaps as the outputs, the gaps are not directly measured on the tokamak, so it may be convenient to apply state-space coordinate transformation, replacing any 6 currents with gaps in the state vector and removing gaps from the outputs.Furthermore, use the ZOH (zero-order hold) for discretization with sample time T s = 0.1 ms such that The final array of discrete-time models in the state-space form is obtained The models have 8 inputs u = δU, 24 states x = [δg T , δ ÎT ] T consisting of 6 gaps and truncated to 18 elements current vector Î, 33 outputs y = [δI p , δI T c , δΨ T ] T directly corresponding to the values measured by the diagnostics at Globus-M2 tokamak: plasma current, 8 currents in control coils, poloidal magnetic flux from 21 loops, quadrupole magnetic flux, and vertical and horizontal dipole magnetic flux.The inclusion of the gaps in the state vector is convenient for some applications, one of which is described in the next section.

Plasma Shape Identification by Robust Observer Synthesized by LMI
The idea of gap estimation with a robust discrete state observer is as follows.Using the FCDI code, a series of LPV models for a series of plasma discharges is computed.The gaps are included in the state vector of all linear models, and the output vector includes the signals measured by the magnetic diagnostics system of the tokamak.Then, using the LMI method, a unified state observer is synthesized, which provides minimal error between states and state estimates for a series of LPV models.
The synthesized observer can be used in a real experiment, with experimental signals connected to its input as shown in Figure 5.
The unified observer for an array of linear models of the plant ensures the minimum error between the state vectors and state estimation and consequently between the gaps values and gaps estimation over the entire discharge duration.This further guarantees the robust behavior of the synthesized plasma shape control system.

Observer
< l a t e x i t s h a 1 _ b a s e 6 4 = " n X x < l a t e x i t s h a 1 _ b a s e 6 4 = " s t 6 M Y P 9 B B a J t J X a q n I c p 7 X q 2 J F 9 g 1 R F 3 Z n Z e P w A C x t i 5 T f 4 B P 4 C p + 1 A W 4 5 k 6 e i c + / I J E s E N u O 6 3 s 7 S 8 s r q 2 X t g o b m 5 t 7 + y W 9 v b r R q W a M p 8 q o X Q z The state equation of the full-order discrete-time observer [20] is given as follows where x is the state estimation vector of discrete-time state-space plant model {A d , B d , C d }, T s is the sample time and L is the matrix of the observer.
Then it is necessary to perform the transition to the error equation of the observer where e = x − x is the error between the states and state estimations.The matrix inequalities systems for the observer synthesis are obtained using the generalized Lyapunov theorem [21]   where the symbol "⊗" denotes the Kronecker product.The poles of the observer are placed in the D-region formed by the disk with the characteristic function The choice of this D-region is due to the need, on the one hand, to provide shorter transition times in the observer compared to the plant model, and on the other hand, the D-region should not be too small; otherwise, it would be impossible to find a solution of the LMI system for the array of plant models.
The synthesizable observer should qualitatively estimate the states for each LTI model from (9), which is obtained from the LPV model (8) for the m th plasma discharge.In addition, the same observer should qualitatively estimate the states for several LPV models corresponding to several discharges.In this approach, the robust performance of the synthesized observer is achieved.
The LMI system for obtaining the observer matrix for the array of models in state-space (9) with the replacement of V = XL is as follows where n = 1, . . ., N m , m = 1, . . ., M and r = 1, . . ., N m M.
The LMI system (11) includes N m M + 1 LMIs, and it must be solved with respect to the two unknown matrices, X and V.The matrix of the observer is defined as Finally, the gap estimation vector δ g is obtained as follows where S g (Figure 5) is the gaps estimation selection matrix from the state vector estimation S g = I 6 0 6,18 , where I 6 is the identity matrix and 0 6,18 is the zeros matrix of the appropriate size.
The comparison of the gaps variations δg derived by LPV model obtained from the FCDI code and the gaps variations estimation δ g obtained by the robust observer synthesized via LMIs for plasma discharge #37263 is shown in Figure 6.Reconstructed plasma shape Estimated plasma shape This is done to obtain matrix K and the base gap values using the data of 50 discharges of the Globus-M2 tokamak.The calculated base gap values are g1 = 0.5235 m, g2 = 0.6264 m, g3 = 0.0217 m, g4 = 0.1058 m, g5 = 0.1590 m, g6 = 0.0278 m. ( The obtained matrix K and base gap values are tested on the discharge #37712 that was not used for identification (Figure 7).The mean squad error (MSE) of all gaps estimation at all moments of time is 1.5 × 10 −5 m 2 .(cm) Reconstructed plasma shape Estimated plasma shape

Neural Network for Plasma Shape Identification
In this section, the identification system based on an artificial neural network is proposed.It is assumed that the input and the output data can be linked using some unknown function f .Neural networks are well known for their ability to approximate unknown functions [22].Attempts to apply them to plasma research in tokamaks began as early as the 1990s.Several major results have been achieved, including the tasks of plasma equilibrium reconstruction [23][24][25][26].However, the vast majority of studies use feed-forward neural networks with multiple hidden layers to approximate the unknown mapping function, which have not shown good results in this problem in the area of generalization to various unknown discharges.To improve this ability, this paper proposes an approach using an encoder-decoder network structure [27].
Neural networks are based on the concept of artificial neurons.The first concept was proposed by Rosenblat [28], called perceptron.It receives inputs (X 1 , X 2 , .., X n ) and sums it with weights (W 1 , W 2 , .., W n ).Then the special function, named the transfer function, is applied to this sum product.The result of the transfer function is the output of the neuron.The most simple neural network, called multilayer perceptron, consists of three layers of neurons: the first one gets the input data, the second one is hidden and processes this data, and the third one is an output layer (Figure 8).
To approximate an unknown function f , a neural network needs to be trained on some given data.The better approximation is achieved by adjusting weights of network's neurons to minimize the value of the loss function, which is computed between the network output and groundtruth values.
In this work, the input and output data are represented as time sequences.Each point in time corresponds to a vector of features, so the data can be described by the matrix with dimensions of time intervals and parameters.The dimensionality in time is equal to 4110, i.e., there are 4110 data vectors at each time point of discharge.The time step is 6.38 × 10 −5 s, so the entire signal is 0.262154 s.The variation in the absolute values of the various parameters, particularly the coil currents, is quite large, so the normalized values are used.They are obtained by subtracting the average value over the entire time sequence from each value of a particular parameter, and then dividing the result by the standard deviation.
However, the time dimension length is not fixed, as the plasma shape parameters are only determined during the divertor phase, when the strike points and corresponding values g 1 and g 2 exist.The start time and duration of this phase are not known from the available for the real-time reconstruction diagnostic (currents and magnetic fluxes) and require further determination.In general, the plasma shape identification problem is dynamic, i.e., the gap values at some point in time during the divertor phase depend not only on the values of magnetic fluxes and coil currents at the same point in time, but also on the values at previous time steps.Therefore, to determine the required parameters, it is advisable to use recurrentbased neural networks.However, it is not practical to use the entire input signal in such a network for several reasons.First, the longer the sequence of the data fed to the recurrent network input, the longer the training and prediction processes, which are important factors in real-time identification.Second, the coordinate values are only significantly affected by data over a relatively small time range.Based on this, the task can be divided into two subtasks.The first one is to determine whether a given moment in time is a divertor phase.The second one is calculation of the required parameters during the already known divertor phase.The first subtask can be solved using a simple feed forward network without the use of recurrence blocks because it is a classification problem, not a regression problem, unlike the second one.In addition, the first subtask is only necessary to limit the length of the input signal to the recurrent network and achieve a simultaneous increase in system speed and improved positioning accuracy.
The values of magnetic fluxes through the loops and coil currents are fed into the network separately.Each input is processed by a densely connected layer, whose outputs are then concatenated.The merged result is fed into two densely connected layers with a dropout between them.This solution is designed to combat overtraining, which has a significant impact in this task because the signals provided for training have a similar structure.The output of the network is the probability that the current time moment belongs to the divertor phase (Figure 9).where Θ is the neural network parameters, x is the network's output value, and y is a label.
The sigmoid function is taken as a transfer function of the output neuron and RelU for the neural network hidden layer ones The Adam [29] optimization algorithm with learning rate α = 0.0001 is used to minimize the loss function.Learning takes place on 50 discharges and the remaining one is left for tests.To measure how often output values match with groundtruth values, the binary accuracy function is utilized.The obtained accuracy of the identification of the divertor phase for all time points equals 0.986 (Figure 10).The second subtask is to determine the required gap values during the divertor phase.As mentioned above, the recurrent neural network based on an encoder-decoder architecture [27] can combat this problem.This type of networks allows to capture temporal dependencies both in input and output data and build mapping between them.The first major block-encoder-created state describing input signal and the second block decoder is responsible for mapping the data into an output sequence.Both the encoder and decoder consist of GRU cells [30].
Figure 11 shows the encoder-decoder network schematic diagram.The network input is divided into two parts: encoder input and decoder input.An input signal is a sequence of vectors with the values of magnetic fluxes and currents, and it is applied to the encoder input.The decoder input can vary.Therefore, it is best to set the decoder input to 0, which will make it work with the dependencies passed to it by the encoder.
where gi are the network's estimation of gaps and coordinates, g i are the groundtruth values, and N is the number of values.This network is also trained on 50 signals and tested on the remaining one.Figure 12 shows the results for the required plasma parameters during the divertor phase of the discharge #37270.
The deviation is calculated for all values of each gap using the MSE.The results obtained have the order of 10 −5 .Reconstructed plasma shape Estimated plasma shape

Real Time Simulation of Identification Systems
To develop and realize plasma control systems for tokamaks, it is effective to apply so-called digital twins with real and digital control systems (Figure 13).This idea is used intensively in the industry because it gives a lot of advantages for the design, modeling, and application of control systems in real time.The digital twin is the interface between the digital and real world because it can have the ability to link physical and virtual worlds in real time, which provides more a realistic and holistic measurement of unforeseeable and unpredictable scenarios [31].
All signals from the magnetic diagnostics system of the tokamak are analog, which are then digitized by passing through an analog-digital converter (ADC).We can simulate the signal digitization process on our real-time test bed (Figures 14 and 15).In Figure 13, one can see the digital twin containing the real dynamical plant with the real feedback controller in the real space and the virtual dynamical plant with the real feedback controller in the virtual space.Between these spaces, there is a feedback of information flows that offers the opportunity to use the results obtained on the digital feedback system for the real control system, and vice versa.The data flows between an existing physical object and a digital object are fully integrated in both directions, which one might refer to as a digital twin [32].The digital twin in the paper consists of the spherical Globus-M2 tokamak with a plasma feedback control system and a test bed with a digital controlled plant model and a feedback controller.The test bed was created by Lomonosov Moscow State University, Trapeznikov Institute of Control Sciences (Moscow), and Ioffe Institute (St.Petersburg).A photo of the test bed for Globus-M2 that is operating in real time is given in Figure 14.In Figure 15, one can see the test bed scheme in detail consisting of the digital plant model and the digital controller with an internal plant model for controller tuning.The digital plant model contains the plasma model in the tokamak and a set of feedback loops for plasma horizontal and vertical position control, for currents control in the poloidal field coils.The digital controller contains plasma equilibrium reconstruction algorithm as well as plasma current and shape controller.
This connection of the two real-time target machines is real and reliable.The realtime test bed is away from sources of powerful electromagnetic radiation, and all of its components have high-quality protection by means of shielding and grounding.
The two identification algorithms described in this paper are applied on the real-time testbed.Figure 16 shows the real-time running of a robust observer synthesized via LMIs.Figure 17 shows the real-time running of the identification algorithm with the static matrix.Real-time simulations are performed with a sample time of 0.1 ms.
These signals demonstrate the workability of two new approaches to reconstruct plasma equilibrium in real time on the test bed.That means important value of these signals in Figures 16 and 17

Comparison of Identification Algorithms
Table 1 shows the comparison results of the different gap identification algorithms.For each gap g and an estimate of this gap g, the value of the MSE (20) is calculated.For each algorithm, the value of the TET (task execution time) in the real-time simulation is given.The FCDI code has an execution time of approximately 25 ms, which is too slow for real-time applications at Globus-M2 tokamak.To apply a plasma shape identification algorithm in real time, the algorithm must have an execution time of less than 1 ms, preferably less than 0.1 ms.All algorithms in Table 1 satisfy that criteria.
The MSE of the robust observer is 100 times smaller than other algorithms.This advantage is due to the fact that the observer is the dynamic model in the state-space form with 24 states.It contains 24 integrators with the help of which the error between the states of the plant model and the estimates of the states at the observer's output is fast minimized.The time it takes to minimize the error is determined by the location of the observer's poles.The observer's poles are defined by the D-region (10).
The disadvantages of using the observer include the fact that it requires the use of scenario values of currents and fluxes, i.e., values relative to which the deviations from gaps are calculated.Other algorithms use the full values of experimental signals as inputs.The synthesis of the observer is possible only in the presence of linear models of the plasma in the tokamak as (9).Linear plant models can be derived only for deviations of currents and fluxes from the scenario values.
The fastest of these estimation algorithms of the gaps between the plasma boundary and the first wall is the static matrix algorithm with a TET of 6.3 µs because it is the simplest and requires only matrix-vector multiplication.The neural network algorithm is attractive because it can be adapted to a large number of discharges during experiments.

Discussion
In this work, the authors developed the original direction of plasma equilibrium reconstruction in D-shaped tokamaks using the magnetic measurements outside the hot plasma [33].The basic criteria of this development are speed of response and accuracy.In practice, there is a set of such approaches, mentioned in the Introduction, which are used on working D-shaped tokamaks all over the world.Some of them use Picard iterations or current filaments methods.However, they rely only on the measurements outside plasma and most do not use the information from the database of the previous plasma discharges.If one uses this information, it may be possible to increase the speed of plasma equilibrium reconstruction.Moreover, when the history information of the plasma discharges is used, one can apply various reconstruction approaches from very simple ones, such as approximation with static matrices, to complex ones, such as artificial neural networks, which can be adjusted by and learn from dynamic processes.It gives a chance to improve not only the process of plasma identification on-line, but to understand the patterns of plasma processes from the experiment.These patterns cannot be deduced from the theory of high-temperature plasma physics because the plasma in a magnetic field is an extremely complicated object.These patterns represent the relationships between the gaps, which are the outputs of the plasma equilibrium reconstruction algorithm applied off-line to a set of plasma discharges and the inputs of this algorithm.The input signals are the measured fluxes on the magnetic loops, the currents in the CS/PF coils, and the plasma current.Then, one can use these patterns to apply any plasma reconstruction algorithms with the highest speed of response, e.g., state observers, static matrices, neural networks, and others.This activity is similar to machine learning techniques, where the search process is automated on big data [34].In future, we can use these patterns for effective plasma control systems design with the fast plasma equilibrium reconstruction algorithms in the feedback with the usage of our new testbed for the installation of plasma control systems in real time on operating tokamaks [16].

Conclusions
The development of the fusion problem is moving forward but not very quickly because the plasma in tokamaks is an extremely complicated plant.In spite of that, the technologies in this field have had great progress and new technologies are appearing.One of these directions is plasma diagnostics to which our research belongs.The algorithms of plasma equilibrium reconstruction, such as ones using static matrix, state observer, and artificial neural network, can be included into the feedback of plasma shape control.The first real-time test of these algorithms is done on the digital model of the plasma shape control system (Figures 2 and 15).After that, the control system can be used in a real experiment on the Globus-M2 tokamak by means of a controller based on the third machine of the digital complex shown in Figures 14 and 15.The third machine will be connected to the tokamak as the control unit of the real control system, like in Figure 13.The real control system will interact with the virtual control system shown in Figures 14 and 15, realizing the concept of the digital twin shown in Figure 13.This approach is in line with the digital twins which are applied in Industry 4.0 [35].
In any case, there is a critical point in this new identification approach.The point is that this approach greatly increases the response rate of plasma equilibrium reconstruction, but the estimation accuracy may not be as high as, for example, in the filaments (current rings) approach.So, the designer of the magnetic plasma control system should choose what is more adequate for the specific control problem since the plasma equilibrium reconstruction algorithm is included in the feedback (Figure 15).
4 = " O j j d 8 s z n P / D p o f v / R T V c Y H o a T 3 4 = " > A A A B 7 H i c b Z C 7 S g N B F I b P x l u M t 1 U L C 5 v B K F i F X Q u 1 M 5 B G u w h u E k i W M D u Z T Y b M z i 4 z s 0 J Y 8 g w 2 F o r Y + g q + h J X Y + A i + g a W T S 6 G J P w x 8 / P 8 5 z D k n S D h T 2 n E + r d z C 4 t L y S n 6 1 s L a + s b l l b + / U V J x K Q j 0 S 8 1 g 2 A q w o Z 4 J 6 m m l O G 4 m k O A o 4 r Q f 9 y i i v 3 1 K p W C x u 9 C C h f o S 7 g o W M Y G 0 s 7 6 q d k W H b L j o l Z y w 0 D + 4 U i h f f b 6 8 f e 5 W v a t t + b 3 V i k k Z U a M K x U k 3 X S b S f Y a k Z 4 X R Y a K W K J p j 0 c Z c 2 D Q o c U e V n 4 2 G H 6 M g 4 H R T G 0 j y h 0 d j 9 3 Z H h S K l B F J j K C O u e m s 1 G 5 n 9 Z M 9 X h u Z 8 x k a S a C j L 5 K E w 5 0 j E a b Y 4 6 T F K i + c A A J p K Z W R H p Y Y m J N v c p m C O 4 s y v P Q + 2 k 5 J 6 W 3 G u n W D 6 E i f K w D w d w D C 6 c Q R k u o Q o e E G Bw B w / w a A n r 3 n q y n i e l O W v a s w t / Z L 3 8 A E g W k 2 Q = < / l a t e x i t > t e x i t s h a 1 _ b a s e 6 4 = " O j j d 8 s z n P / D p o r b d 5 a 8 F a z B z C H 1 j v P 7 b 5 l P 4 = < / l a t e x i t > U ref < l a t e x i t s h a 1 _ b a s e 6 4 = " O j j d 8 s z n P / D p o

Figure 2 .
Figure 2. Structure scheme of the plasma position, current and shape control system of Globus-M2 tokamak.
r 9 a b 9 T 4 r X b P m P U e w A O v j F / h 9 m o 4 = < / l a t e x i t > S g < l a t e x i t s h a 1 _ b a s e 6 4 = " v S + I x v a Y / U E 3 m e w A y a a G S W j Z z q s s 0 4 2 j D m r g y 3 I 1 K S T M o 2 G z 8 6 S R m y f V B 2 T 8 s + 5 d e 6 W y X j J E n 2 2 S H 7 B O f H J E z c k E q p E o 4 u S M P 5 J E 8 O f f O s / P i v I 5 b c 8 7 v z B a Z g P P 2 A / 3 8 o e o = < / l a t e x i t > g < l a t e x i t s h a 1 _ b a s e 6 4 = " 1 F

Figure 5 .
Figure 5. Robust observer synthesized via LMIs for use in a real experiment.The red vector signal includes experimental signals obtained by the magnetic diagnostic system.The blue vector signal includes voltages on the poloidal coils.The yellow signal contains states estimation, which includes gaps estimation.

Figure 6 .
Figure 6.Comparison of gaps variations δg derived by LPV model obtained from the FCDI code (blue line) and estimation of gap variations δ g obtained from robust observer synthesized by LMIs (red line).Globus-M2 discharge #37263.

Figure 7 .
Figure 7. Estimation of the gap values in the discharge #37712 of the Globus-M2 spherical tokamak.

Figure 9 .
Figure 9. FNN model.After this, the binary crossentropy is used as the loss function to measure the difference between the network output and training data BCS(Θ, y) = −y log(x) + (1 − y) log(1 − x),

Figure 11 .
Figure 11.Encoder-decoder model.The MSE as loss function and linear function as transfer function for the output neuron have the best performance for this regression problem.

Figure 12 .
Figure 12.Neural network estimation of the gap values.Discharge #37270 of the Globus-M2 spherical

FeedbackFigure 13 .
Figure 13.Digital twin containing a real dynamical plant with real feedback controller and virtual dynamical plant with a real feedback controller closed by the feedback of information flows: real-time data and algorithms, commands, adaptations, and recommendations.

Figure 14 .
Figure 14.Real-time test bed for plasma control in tokamaks.The test bed consists of two Speedgoat performance real-time target machines that are connected in feedback: one computer plays the role of the controlled plant model and the other one is the MIMO controller (https://www.ipu.ru/presscenter/62866,accessed on 21 December 2021).

Figure 15 .
Figure 15.Scheme of digital twin of Globus-M2.The block with plasma equilibrium reconstruction algorithm is marked by red color. .

Figure 16 .
Figure 16.Real-time simulation of gaps variations δg derived by LPV model obtained from the FCDI code (blue line) and estimation of gap variations δ g obtained from robust observer synthesized by LMIs (red line).Discharge #37263.

Figure 17 .
Figure 17.Real-time simulation of gaps g derived by LTI model obtained from the FCDI code (blue line) and estimation of gaps g obtained from static matrix (red line).Discharge #37263.

Table 1 .
Comparison of identification results.