Development of Data Assimilation System for Temperature and Density Control in Helical Plasmas

Abstract

We develop a real-time adaptive predictive control system based on data assimilation (DA) for the temperature and density of helical fusion plasmas. The DA-based control approach enables the harmonious integration of measurement, heating, fueling, and simulation and can provide a flexible platform for adaptive model predictive control. The core part of the control system, ASTI, is built upon the integrated simulation code TASK3D and a data assimilation framework DACS. DACS integrates adaptation of the predictive model (digital twin) to the actual system using real-time measurements and control estimation that is robust against model and observation uncertainties. We perform numerical experiments using ASTI to control the electron temperature profile and density of a virtual plasma generated by TASK3D. The results demonstrate that ASTI can effectively drive the virtual plasma state toward the target state while bridging the gap between the digital twin and the virtual plasma. Furthermore, the numerical experiments clarify the effects of hyperparameters in the DA-based control approach on control performance.

1.  Introduction

Future fusion reactors will require state estimation of the fusion plasma based on limited observations and harmonious control of various actuators. One powerful approach to address this challenge is digital twin control, which captures the internal state of the system and estimates the optimal control input for the target state by leveraging a virtual model of a physical object (digital twin) [1–3]. However, there are several problems in achieving a digital twin control for fusion plasmas. First, it is inherently difficult to model all the components of fusion plasma and their interactions with sufficient accuracy. Second, it is difficult to predict the behavior of the entire plasma in real time with high accuracy. Furthermore, nonlinear model predictive control generally requires iterations of predictive calculations. Therefore, the actual prediction must be sufficiently faster than real time. To overcome these challenges, we need to develop a system in which measurements, actuators (or technologies for individual control problems), and the digital twin, based on physics knowledge and accumulated experimental data, are closely coupled and complement each other.

We have introduced a control approach based on data assimilation (Data Assimilation and Control System, DACS) [4] to achieve a digital twin control system for fusion plasmas. Data assimilation (DA) is a statistical method for estimating the state vector, which consists of the variables in a numerical model, based on observation data [5–7]. It can make the behavior of the model similar to that of the actual system. DACS is a DA framework that integrates digital twin adaptation using real-time measurements and control estimation that is robust to model and observation uncertainties. We are developing a data assimilation system, ASTI [4, 8, 9], based on the DACS framework and an integrated simulation code, TASK3D [10, 11]. ASTI approximates the probability distribution of the state vector (state distribution) with many ensemble members (simulations with slightly different conditions) and performs time evolution and assimilation calculations. The effectiveness of the DA-based control was demonstrated through numerical experiments [4] and a simple control experiment in Large Helical Device (LHD) [8].

In recent years, research on data-driven approaches for plasma control (e.g., control of instability and avoidance of sudden phenomena) has been active [12–15]. Our approach has both physics-based and data-driven aspects to build a comprehensive control system in which many observations and actuators are harmoniously combined. Observations compensate for the digital twin’s imperfections to enhance the prediction and control performance, and the digital twin interconnects various observations to help the state estimation. In ASTI, physical knowledge and control constraints can be easily taken into account by incorporating them into the state vector and the digital twin. ASTI can be applied to complex control with multiple variables, considering both observable and unobservable variables.

In this study, we extend ASTI for simultaneous predictive control of temperature and density, considering electron cyclotron heating (ECH), neutral beam injection (NBI) heating, and gas-puff as actuators. Real-time prediction of heat transport in ECH plasmas has already been achieved in the previous study [8]. In this study, we accelerate the computation of particle transport by developing a surrogate model to evaluate the particle source due to neutral particles. We perform numerical experiments using ASTI to simultaneously control the electron temperature profile and density of a virtual LHD plasma. In addition, we investigate the control performance by varying the hyperparameters and varying the gaps between the digital twin and the actual system.

This paper is organized as follows. Section 2 describes the DACS framework and the real-time prediction model of helical plasma (TASK3D). Section 3 explains the numerical experiments to demonstrate and investigate the control performance. The results of the numerical experiments are shown in Sec. 4. The conclusions are given in Sec. 5.

2.  DA-Based Control for Helical Plasmas

2.1  Data assimilation framework for adaptive model predictive control

Consider a control problem where control estimation and measurement of the system state are performed at time interval $\mathrm{\Delta }t$. We use the discrete time $t_i=t_0+i\mathrm{\Delta }t$, where $t_0$ denotes the initial time. The DACS framework is based on the following state space model:

  

$${\mathbf{x}}_{i+1}=f_{i+1}({\mathbf{x}}_i,{\mathbf{v}}_{i+1}),$$(1)

  

$${\mathbf{z}}_i=h_i^{\mathbf{z}}({\mathbf{x}}_i,{\mathbf{w}}_i^{\mathbf{z}}),$$(2)

  

$${\mathbf{u}}_i^{*}=H^{\mathbf{u}}{\mathbf{x}}_i+{\mathbf{w}}_i^{\mathbf{u}},$$(3)

  

$${\mathbf{y}}_i=h_i^{\mathbf{y}}({\mathbf{x}}_i,{\mathbf{w}}_i^{\mathbf{y}}).$$(4)

The state vector at time $t_i$, ${\mathbf{x}}_i$, is defined by

  

$${\mathbf{x}}_i=\left(\begin{array}{c}\hfill {\tilde{\mathbf{x}}}_i\hfill \\ \hfill {\mathbf{u}}_i\hfill \end{array}\right),$$(5)

where ${\tilde{\mathbf{x}}}_i$ is part of the state vector containing the system state and model parameters, which is used as the state vector in typical data assimilation. Vector ${\mathbf{u}}_i$ is the control input that determines the time evolution of the system from time $t_{i-1}$ to $t_i$.

Equation (1) is the system model that describes the time evolution of the system, ${\mathbf{x}}_i\to {\mathbf{x}}_{i+1}$. Here, $f_{i+1}$ is a time evolution operator considering the system noise ${\mathbf{v}}_{i+1}$ and corresponds to a digital twin. In this study, we assume that ${\mathbf{u}}_i$ is constant in the prediction interval, and Eq. (1) can be broken down as

  

$${\mathbf{u}}_{i+1}={\mathbf{u}}_i+{\mathbf{v}}_{i+1}^{\mathbf{u}},$$(6)

  

$${\tilde{\mathbf{x}}}_{i+1}={\tilde{f}}_{i+1}({\tilde{\mathbf{x}}}_i,{\mathbf{u}}_{i+1},{\tilde{\mathbf{v}}}_{i+1}).$$(7)

We add the system noise for control input, ${\mathbf{v}}_{i+1}^{\mathbf{u}}$, to ${\mathbf{u}}_i$ by Eq. (6) before the time evolution calculation by Eq. (7). In Eq. (7), ${\tilde{\mathbf{v}}}_{i+1}$ is the system noise for $\tilde{\mathbf{x}}$ and is added to $\tilde{\mathbf{x}}$ before or after the time evolution. Equation (2) represents the relationship between the state vector ${\mathbf{x}}_i$ and target-state vector ${\mathbf{z}}_i$ considering the target-state noise ${\mathbf{w}}_i^{\mathbf{z}}$. In the same manner, Eqs. (3) and (4) represent the relationships between ${\mathbf{x}}_i$ and optimal control input ${\mathbf{u}}_i^{*}$ with control-input noise ${\mathbf{w}}_i^{\mathbf{u}}$, and that between ${\mathbf{x}}_i$ and observation vector ${\mathbf{y}}_i$ with observation noise ${\mathbf{w}}_i^{\mathbf{y}}$, respectively. The relation between ${\mathbf{x}}_i$ and ${\mathbf{u}}_i$ is written by a matrix $H^{\mathbf{u}}$ from Eq. (5). Note that Eqs. (2), (3) and (4) are referred to as “target-state model”, “control-input model”, and “observation model”, respectively. In this study, we assume that the system noise ${\mathbf{v}}_i$ follows a Gaussian distribution with zero mean and covariance matrix $Q_i$, i.e., ${\mathbf{v}}_i\sim N(\mathbf{0},Q_i)$. Similarly, ${\mathbf{w}}_i^{\mathbf{z}}$, ${\mathbf{w}}_i^{\mathbf{u}}$, and ${\mathbf{w}}_i^{\mathbf{y}}$ are assumed to follow the probability distributions $N(\mathbf{0},R_i^{\mathbf{z}})$, $N(\mathbf{0},R_i^{\mathbf{u}})$, and $N(\mathbf{0},R_i^{\mathbf{y}})$, respectively. These covariance matrices are important hyper-parameters of the control system that determine the control capability, as discussed in Sec. 3.2.

The state distribution $p(\mathbf{x})$ is approximated by an ensemble consisting of ${\{{\mathbf{x}}^{(k)}\}}_{k=1}^{k=n}$, where $k$ and $n$ are the index and the number of ensemble members. The computations for the ensemble in the DACS framework are divided into two main processes: prediction and filtering (assimilation). The former performs a time-evolution calculation for each ensemble member based on the system model (1), and the ensemble approximating the state distribution $\mathrm{\Delta }t$ ahead is obtained. In the latter process, information imposition on a state distribution is computed for system state estimation (adaptation) and control estimation. The ensemble approximating the state distribution reflecting the information $\mathbf{d}$, $p(\mathbf{x}|\mathbf{d})$, is obtained, where $\mathbf{d}=\mathbf{z}\mathbf{,}{\mathbf{u}}^{\mathbf{*}}$, or $\mathbf{y}$. This process is performed based on Bayesian filters (e.g., the ensemble Kalman filter [16] and the particle filter [17]) for the information ($\mathbf{z}$, ${\mathbf{u}}^{\mathbf{*}}$, $\mathbf{y}$) to be imposed (assimilated) using the models (2)–(4).

The control estimation is executed by assimilating the target state ${\mathbf{z}}_i$ into the state distribution $p(\mathbf{x})$ to obtain $p({\mathbf{x}}_i|{\mathbf{z}}_i)$ (z-filter). We obtain the probability distribution of the control input ${\mathbf{u}}_i$ required to produce the target state ${\mathbf{z}}_i$, $p({\mathbf{u}}_i|{\mathbf{z}}_i)$, by marginalizing the z-filtered distribution $p({\mathbf{x}}_i|{\mathbf{z}}_i)$ for ${\tilde{\mathbf{x}}}_i$. In this study, the optimal control input ${\mathbf{u}}_i^{*}$ is defined as the expected value of ${\mathbf{u}}_i$ for the distribution $p({\mathbf{u}}_i|{\mathbf{z}}_i)$. The estimated control input ${\mathbf{u}}_i^{*}$ is reflected to the state distribution by assimilating ${\mathbf{u}}_i^{*}$ into the predicted distribution (u-filter). The u-filtered distribution $p({\mathbf{x}}_i|{\mathbf{u}}_i^{*})$ represents the predicted distribution when ${\mathbf{u}}_i^{*}$ is input. In addition to control estimation, the system state, including the model parameters, can be estimated by assimilating the observation ${\mathbf{y}}_i$ into the latest distribution (y-filter). The DACS framework realizes adaptive predictive control by combining these operations on the state distribution: prediction, z-filter, u-filter, and y-filter. See Ref. [4] for the detail of the DACS framework and Sec. 3.2 for the specific control algorithms employed in this study.

2.2  Real-time prediction of helical plasmas

We employ the integrated transport simulation code TASK3D as a digital twin of helical plasma (the operator $f$ in Eq. (1)). This section gives an overview of TASK3D and the extensions made to achieve real-time prediction of temperature and density. TASK3D treats the radial heat and particle transport issues in a torus-shaped plasma as one-dimensional (1D) problems for the normalized minor radius $\rho$. In this study, we assume that the electron density and ion density profiles are the same, i.e., $n=n_{\mathrm{e}}=n_{\mathrm{i}}$. The particle and heat transport equations for each electron and ion species are solved:

  

$$\frac{\partial }{\partial t}\left(n{𝒱}^{\prime }\right)=-\frac{\partial }{\partial \rho }{𝒱}^{\prime }(⟨|\mathrm{\nabla }\rho |⟩nV-⟨|\mathrm{\nabla }\rho {|}^2⟩D\frac{\partial n}{\partial \rho })+S{𝒱}^{\prime },$$(8)

  

$$\begin{array}{ccc}\frac{\partial }{\partial t}\left(\frac{3}{2}n{T}_s{{𝒱}^{\prime }}^{5/3}\right)& =& -{{𝒱}^{\prime }}^{2/3}\frac{\partial }{\partial \rho }{𝒱}^{\prime }[⟨|\mathrm{\nabla }\rho |⟩n{T}_s(V_{K_s}+\frac{3}{2}V)\hfill \\ & & -⟨|\mathrm{\nabla }\rho {|}^2⟩\frac{3}{2}{\mathit{DT}}_s\frac{\partial n}{\partial \rho }-⟨|\mathrm{\nabla }\rho {|}^2⟩n{\chi }_s\frac{{\mathit{\partial T}}_s}{\partial \rho }]+P_s{{𝒱}^{\prime }}^{5/3},\hfill \end{array}$$(9)

where $n$ and $T_s$ are the density and temperature of the s-species, and $⟨⟩$ represents the magnetic flux surface average. In addition, $𝒱$ is the plasma volume, and ${𝒱}^{\prime }=\mathrm{d}𝒱/\mathrm{d}\rho$. The equilibrium magnetic field is computed by the VMEC code [18], which calculates the three-dimensional magnetohydrodynamics equilibrium. Here, we use a typical LHD magnetic configuration: the major radius of the magnetic axis at the vacuum is 3.6 m, and the magnetic field strength at the plasma center is 2.85 T.

Coefficients $D$ and ${\chi }_s$ are the particle and thermal diffusivities, and the convection velocities $V$ and $V_{K_s}$ are the particle and heat pinch velocities, respectively. The $S$ and $P_s$ terms are the particle and heat source, respectively. The heating power source $P_s$ comprises external heating, power exchange between species, and the loss term due to interaction with neutrals. For real-time computation of ECH, we employ the following simple model:

  

$$P_{\mathrm{e}}^{\mathrm{ECH}}(\rho )=\sum _{\xi }^{}{p}_{\xi }(\rho ),$$(10)

  

$$p_{\xi }(\rho )=A_{\xi }\exp \left(-\frac{1}{2}\frac{({\mu }_{\xi }-\rho {)}^2}{{\sigma }_{\xi }^2}\right),$$(11)

where $\xi$ is the index of gyrotron, and ${\mu }_{\xi }$ and ${\sigma }_{\xi }$ are the corresponding parameters determining the radial profile of the heat deposition. The coefficient $A_{\xi }$ is determined from the input power ${𝒫}_{\xi }$ by

  

$${𝒫}_{\xi }={\int }_0^1{p}_{\xi }(\rho ){𝒱}^{\prime }(\rho )\mathrm{d}\rho .$$(12)

For real-time computation of NBI heating, we employ the FIT3D-RC module [19]. This module evaluates the beam ion birth profile using the Gaussian process regression model applied to precomputation results by the Monte Carlo simulation. It calculates the heat deposition and particle source profiles based on the simple analytical solution of the Fokker-Planck equation. The particle source $S$ is primarily determined by the ionization of neutral particles evaluated by the AURORA module [11, 20] of TASK3D. The AURORA module calculates the component of $S$ using the plasma profiles and the density and temperature of neutral particles at the plasma edge.

Real-time calculations are already possible for heat transport considering ECH and NBI, but they have not yet been achieved for particle transport. In TASK3D, the AURORA module accounts for most of the computational time for the particle transport. Although the AURORA module is essential for computing the time evolution of the density profile, it cannot be used for real-time prediction because it is based on the Monte Carlo method, which is computationally expensive. To achieve real-time particle transport computation, we have developed a neural network-based surrogate model of AURORA. The input variables are the plasma profiles ($T_{\mathrm{e}}$, $T_{\mathrm{i}}$, $n$) and the density and temperature of neutral particles at the plasma edge ($n_{\mathrm{n}}$, $T_{\mathrm{n}}$). Each profile has 60 radial grid points; thus, the dimension of the input vector is 182. The output is the particle source due to the ionization of neutral particles, $S^{\mathrm{neu}}$, (60 radial points). The surrogate model is built by using a multilayer perceptron with four hidden layers (182, 1,000, 1,000, 180 units). We employ ReLU as the activation function and Adam as the optimizer. This structure was selected from a range of configurations with up to five hidden layers and up to 1,000 units per layer. It was confirmed that increasing the number of parameters beyond this structure does not significantly improve prediction accuracy.

To train the model, 500,000 input-output data sets were used (400,000 as training data and 100,000 as test data). These data sets were generated by the AURORA computation from randomly generated input vectors. The input radial profiles of $T_{\mathrm{e}}$, $T_{\mathrm{i}}$, and $n$ were generated using the following function [21]:

  

$$A(\rho )=c_1\tanh [-c_2(\rho -c_3)]+c_4.$$(13)

The parameters $c_2$ and $c_3$, which determine the shape of the radial profile, are randomly generated from uniform distributions ($2<c_2<15$ and $0.3<c_4<1$). The parameters $c_1$ and $c_4$ are determined by the values $A(0)$ and $A(1)$, which are also randomly generated from uniform distributions. The range of each variable was set as follows:

  

$$\begin{array}{ccc}\hfill 0.1<& \hfill T_{\mathrm{e}}(0),T_{\mathrm{i}}(0)\hfill & <10(\mathrm{keV}),\hfill \\ \hfill 0.1<& \hfill T_{\mathrm{e}}(1),T_{\mathrm{i}}(1)\hfill & <1(\mathrm{keV}),\hfill \\ \hfill 0.01<& \hfill n(0)\hfill & <1({10}^{20}{\mathrm{m}}^{-3}),\hfill \\ \hfill 0.01<& n(1)& <0.1({10}^{20}{\mathrm{m}}^{-3}),\\ \hfill 0<& \hfill n_{\mathrm{n}}\hfill & <2({10}^{17}{\mathrm{m}}^{-3}),\hfill \\ \hfill 1<& \hfill T_{\mathrm{n}}\hfill & <50(\mathrm{eV}).\hfill \end{array}$$

Here, profiles where $A(0)<A(1)$ were not used for the training. Furthermore, as a quarter of the input density profiles, hollow profiles generated by the following function were used:

  

$$n_{\mathrm{h}}(\rho )=w_1\tanh [-6(\rho -w_2)]+w_3\tanh [-w_4(\rho -1)]+w_5.$$(14)

The parameters $w_1$, $w_3$, and $w_5$ are determined by the parameter $H=({\max }_{\rho }n-n(0))/n(0)$ and randomly generated $n(0)$ and $n(1)$. The parameters $w_2$, $w_4$, and $H$, which determine the shape of the hollow profile, are randomly selected from a set of reasonable values, $(w_2,w_4,H)$ ∈ {(0.3, 2, 0.2), (0.3, 2, 0.4), (0.5, 4, 0.2), (0.5, 4, 0.4), (0.7, 6, 0.2), (0.7, 6, 0.4)}.

Figure 1 shows a test result of the surrogate model for $S^{\mathrm{neu}}$ ($\rho$ = 0.9). The coefficient of determination is 0.995, which shows that the surrogate model can reproduce the AURORA simulation results with high accuracy. Figure 2 shows the TASK3D simulation results using the AURORA module and the surrogate model. In the simulations, we assumed $D(\rho )$ = 0.8 m²s⁻¹ and the typical LHD magnetic configuration, with a major radius of 3.6 m and a magnetic field strength at the plasma center of 2.85 T. It can be seen that the results of the integrated simulation using AURORA are reproduced with high accuracy, both for the density profile as shown in Fig. 2(a) and for the particle source as shown in Fig. 2(b).

Fig. 1.  Comparisons between the AURORA simulation results (“Simulation”) and the surrogate model prediction (“NN prediction”) of $S^{\mathrm{neu}}$ (10²⁰ m⁻³s⁻¹) at $\rho$ = 0.9.

Fig. 2.  Integrated simulation results of (a) the density and (b) the particle source calculated by using the AURORA module (“Simulation”) and the surrogate model (“NN”) in three cases: (i) $n_{\mathrm{n}}$ = 1 × 10¹⁶ m⁻³, (ii) $n_{\mathrm{n}}$ = 5 × 10¹⁶ m⁻³, and (iii) $n_{\mathrm{n}}$= 1 × 10¹⁷ m⁻³.

The TASK3D simulation using the surrogate model takes 0.04 seconds to perform a particle transport calculation for one second, while the simulation using the AURORA module takes 9.5 seconds. Here, this calculation was performed in a vector computer (NEC, SX-Aurora TSUBASA), and the calculation of the $S^{\mathrm{nue}}$ term in TASK3D was performed every 10 ms. The time required for particle transport calculations has been significantly reduced. Even when coupling heat transport considering ECH and NBI, the transport simulation could be performed in 130 ms. This surrogate model has sufficient computational speed and accuracy for our DA-based control system. It enables real-time predictive control including the plasma density.

3.  Numerical Experiments to Control Virtual LHD Plasma

We perform numerical experiments to control a virtual plasma (hydrogen plasma) in LHD to demonstrate the effectiveness of the DA-based control system and investigate the dependence of the control performance on the hyperparameters.

3.1  Control problem

Consider simultaneous control of electron temperature profile and density. As the actuators, we assume gas-puff to control the density and ECH with two separate heating positions (around $\rho$ = 0.1 and $\rho$ = 0.4) to control the radial profile of the electron temperature. The ECH parameters in Eq. (11) are set as shown in Table 1 based on the five gyrotrons implemented in LHD. Here, we introduce the ECH power to around $\rho$ = 0.1, ${𝒫}_{0.1}={𝒫}_1+{𝒫}_2$, and the power to around $\rho$ = 0.4, ${𝒫}_{0.4}={𝒫}_3+{𝒫}_4+{𝒫}_5$. In addition, we introduce the parameter $𝒢$ to simulate the gas-puff. We assume that the neutral density at the plasma edge, $n_{\mathrm{n}}$, is determined by the linear model, $n_{\mathrm{n}}=0.6\times 𝒢$ (10¹⁶ m⁻³).

Table 1. ECH parameters employed in the numerical experiments.

$\xi$	${\mu }_{\xi }$	${\sigma }_{\xi }$	available ${𝒫}_{\xi }$ (MW)
1	0.1	0.015	0, 0.15, 0.3
2	0.1	0.02	0, 0.2, 0.4
3	0.4	0.025	0, 0.2, 0.4
4	0.4	0.03	0, 0.3, 0.6
5	0.4	0.03	0, 0.3, 0.6

Suppose that the electron temperature and density radial profiles are available using the real-time Thomson scattering measurement system in LHD [8, 22]. Table 2 shows the state, target, and observation variables used in the numerical experiments. The radial profile in the state vector is defined on 11 grid points ($\rho$ = 0, 0.1, 0.2, …, 1), while the profile in TASK3D is defined on 60 grid points. We put the parameters $c_{\mathrm{e}}$, $c_{\mathrm{i}}$, $d$, and $v$ into the state vector to consider the uncertainties in the thermal and particle diffusivities and particle convection velocity [9]. Denoting the coefficients used in ordinary simulations by variables with ′, the transport parameters used in ASTI are defined by ${\chi }_{\mathrm{e}}=c_{\mathrm{e}}{\chi }_{\mathrm{e}}^{\prime }$, ${\chi }_{\mathrm{i}}=c_{\mathrm{i}}{\chi }_{\mathrm{i}}^{\prime }$, and $D=d{D}^{\prime }$, $V=V^{\prime }+v$. The base models for transport parameters in ASTI are assumed as follows:

Table 2. State, target, and observation variables with their dimensions in the vectors ($M_j$).

Variable			$M_j$
$\tilde{\mathbf{x}}$	$T_{\mathrm{e}}$	Electron temperature	11
	$T_{\mathrm{i}}$	Ion temperature	11
	$n$	Density	11
	$c_{\mathrm{e}}$	Factor for electron thermal diffusivity	11
	$c_{\mathrm{i}}$	Factor for ion thermal diffusivity	11
	$d$	Factor for particle diffusivity	11
	$v$	Additional particle convection velocity	10
$\mathbf{u}$	${𝒫}_{0.1}$	ECH input power for $\rho$ = 0.1	1
	${𝒫}_{0.4}$	ECH input power for $\rho$ = 0.4	1
	$𝒢$	Gas-puff input	1
$\mathbf{z}$	$T_{\mathrm{e},\rho =0}$	Electron temperature at $\rho$ = 0	1
	$T_{\mathrm{e},\rho =0.25}$	Electron temperature at $\rho$ = 0.25	1
	$n_{\rho =0.25}$	Electron density at $\rho$ = 0.25	1
$\mathbf{y}$	$n$	Density	11
	$T_{\mathrm{e}}$	Electron temperature	11

  

$${\chi }_{\mathrm{e}}^{\prime }(\rho )={\chi }_{\mathrm{e}}^{\mathrm{const}},$$(15)

  

$${\chi }_{\mathrm{i}}^{\prime }(\rho )=C_{\mathrm{i}}\left(\frac{T_{\mathrm{i}}}{eB}\right)\left(\frac{{\rho }_{\mathrm{i}}}{a}\right)\left(\frac{\mathrm{\nabla }{T}_{\mathrm{i}}}{T_{\mathrm{i}}}a\right),$$(16)

  

$$D^{\prime }(\rho )=D^{\mathrm{const}},$$(17)

  

$$V^{\prime }(\rho )=0,$$(18)

where $B$, ${\rho }_{\mathrm{i}}$, and $a$ are the magnetic field strength, the ion Larmor radius, and the plasma minor radius, respectively. We set the constant parameters as ${\chi }_{\mathrm{e}}^{\mathrm{const}}$ = 2 (m²/s), $C_{\mathrm{i}}$ = 0.57, and $D^{\mathrm{const}}$ = 1 (m²/s).

3.2  Control algorithm

We have constructed an algorithm of adaptive model predictive control using the DACS framework for the numerical experiments. The control and observation interval $\mathrm{\Delta }t$ is set to 0.3 seconds, considering the computation time of the prediction and the filtering, the time delay of the Thomson scattering measurement, and the communication time. The algorithm has been constructed so that the measured information is most quickly reflected in the control estimation. It comprises the following steps on the state distributions.

• Prediction

  

$$p\left({\mathbf{x}}_i\right|{\mathbf{y}}_{0:i-1},{\mathbf{u}}_{0:i}^{*})\to p\left({\mathbf{x}}_{i+1}\right|{\mathbf{y}}_{0:i-1},{\mathbf{u}}_{0:i}^{*}),$$(19)

• y-filter

  

$$p\left({\mathbf{x}}_{i+1}\right|{\mathbf{y}}_{0:i-1},{\mathbf{u}}_{0:i}^{*})\to p\left({\mathbf{x}}_{i+1}\right|{\mathbf{y}}_{0:i},{\mathbf{u}}_{0:i}^{*}),$$(20)

• z-filter

  

$$p\left({\mathbf{x}}_{i+1}\right|{\mathbf{y}}_{0:i},{\mathbf{u}}_{0:i}^{*})\to p\left({\mathbf{u}}_{i+1}\right|{\mathbf{y}}_{0:i},{\mathbf{u}}_{0:i}^{*},{\mathbf{z}}_{i+1}),$$(21)

• u-filter

  

$$p\left({\mathbf{x}}_{i+1}\right|{\mathbf{y}}_{0:i},{\mathbf{u}}_{0:i}^{*})\to p\left({\mathbf{x}}_{i+1}\right|{\mathbf{y}}_{0:i},{\mathbf{u}}_{0:i+1}^{*}).$$(22)

Here the subscript $t_1:t_2$ denotes all the timings in $[t_1,t_2]$. Given the distribution $p({\mathbf{x}}_i|{\mathbf{y}}_{0:i-1},{\mathbf{u}}_{0:i}^{*})$, the prediction step computes the time evolution of the state distribution to $p({\mathbf{x}}_{i+1}|{\mathbf{y}}_{0:i-1},{\mathbf{u}}_{0:i}^{*})$ based on the system model. This step can be performed by computing the time evolution of each ensemble member approximating the distribution $p({\mathbf{x}}_i|{\mathbf{y}}_{0:i-1},{\mathbf{u}}_{0:i}^{*})$. After the next observation ${\mathbf{y}}_i$ is obtained, ASTI assimilates the observation to the predicted distribution using the y-filter and optimizes the state vector (adaptation). After the y-filtering, ASTI assimilates the target ${\mathbf{z}}_{i+1}$ to the y-filtered distribution and determines the next control input by

  

$${\mathbf{u}}_{i+1}^{*}=E\left[{\mathbf{u}}_{i+1}\right|{\mathbf{y}}_{0:i},{\mathbf{u}}_{0:i}^{*},{\mathbf{z}}_{i+1}].$$(23)

This control input is sent to the actuators immediately after the z-filtering. This determined control input information is reflected in the y-filtered distribution using the u-filter. This u-filtered distribution $p({\mathbf{x}}_{i+1}|{\mathbf{y}}_{0:i},{\mathbf{u}}_{0:i+1}^{*})$ is the distribution $\mathrm{\Delta }t$ ahead from the given distribution $p({\mathbf{x}}_i|{\mathbf{y}}_{0:i-1},{\mathbf{u}}_{0:i}^{*})$; thus this control algorithm can proceed to the next prediction step. In this study, we implement all filters using the EnKF and take 390 ensemble members. This is the maximum number of ensemble members that can be computed by the actual DA-based control system in LHD [8].

We perform only observation assimilation with a fixed control input ${\mathbf{u}}_{\mathrm{init}}=({𝒫}_1,{𝒫}_2\mathcal{,}𝒢)=(0.3,0.4,1)$ during the first phase $t$ < 2.1 seconds to get the digital twin close to the actual system (virtual LHD) at the beginning of control (2.1 seconds). During this initial adaptation phase, the system noise is not applied to the control input. ASTI assimilates the observations up to 1.8 seconds, after which the system noise is added to the control input, and the control estimation begins. Subsequently, the control estimation starts at 2.1 seconds to produce the target state.

The covariance matrices $Q_i$, $R_i^{\mathbf{z}}$, $R_i^{\mathbf{u}}$, and $R_i^{\mathbf{y}}$ are the key parameters determining the overall control performance. The covariance matrix to generate the initial ensemble members, $V_{\mathrm{init}}$, is also required. In the numerical experiments, diagonal matrices are used for these covariance matrices. The standard deviations of $V_{\mathrm{init}}$, $Q_i$, and $R_i^{\mathbf{u}}$ are set as shown in Table 3. We introduce the parameter ${\sigma }_Q$ to vary the system noise level for the model parameters. System noise for $\mathbf{u}$ is an important parameter determining the range of control inputs considered in a single control estimation and the change rate of ${\mathbf{u}}^{*}$. The standard deviations of the noise for $\mathbf{u}$ are set to sufficiently large values to track the target state. The control input noise ${\mathbf{w}}^{\mathbf{u}}$, which represents the uncertainty in the actually applied control input, affects the assimilation of ${\mathbf{u}}^{*}$ into the state distribution (u-filter). In many cases, the variance can be set sufficiently small within the range that ensures stable data assimilation, as shown in Table 3. The covariance matrix $R^{\mathbf{y}}$ determines the effect of the observed information on the state distribution and affects the adaptation performance. The standard deviation of the observation noise is assumed to be proportional to the difference between the observation and mean of the state distribution [4, 9],

Table 3. Standard deviations of the initial distribution (${\sigma }_{\mathrm{init}}$), system noise (${\sigma }_{\mathrm{sys}}$), and control-input noise (${\sigma }_{\mathit{Ru}}$). The values with % as the unit represent the rate for determining the standard deviation in proportion to the state distribution mean.

Variable		${\sigma }_{\mathrm{init}}$	${\sigma }_{\mathrm{sys}}$	${\sigma }_{\mathit{Ru}}$
$\tilde{\mathbf{x}}$	$T_{\mathrm{e}}$	15%	5%
	$T_{\mathrm{i}}$	15%	5%
	$n$	15%	5%
	$c_{\mathrm{e}}$	0.1	${\sigma }_Q$
	$c_{\mathrm{i}}$	0.1	${\sigma }_Q$
	$d$	0.1	${\sigma }_Q$
	$v$	0.1 m/s	0.1 m/s
$\mathbf{u}$	${𝒫}_{0.1}$	0	0.3 MW	0.03 MW
	${𝒫}_{0.4}$	0	0.3 MW	0.03 MW
	$𝒢$	0	0.3	0.05

  

$${\left(R_i^{\mathbf{y}}\right)}_{ll}=r_y^2{({\mathbf{y}}_i-H^{\mathbf{y}}{\hat{\mathbf{x}}}_i)}_l^2,$$(24)

where ${\hat{\mathbf{x}}}_i$ is the mean of the ensemble approximating $p({\mathbf{x}}_i|{\mathbf{y}}_{0:i-1},{\mathbf{u}}_{0:i}^{*})$. The subscripts $({)}_{ll}$ and $({)}_l$ denote the l-th diagonal component and the l-th element of the vector.

The parameters $r_y$ in Eq. (24), the standard deviation of the system noise for the model parameters (${\sigma }_Q$), and the standard deviation of the target noise (${\sigma }_{Rz}$) are set or varied according to numerical experiments. Note that the same value of ${\sigma }_{Rz}$ is used for the target variables $T_{\mathrm{e},\rho =0}$, $T_{\mathrm{e},\rho =0.25}$, and $n_{\rho =0.25}$ (unit is keV or 10¹⁹ m⁻³). The initial ensemble means of $T_{\mathrm{e}}$, $T_{\mathrm{i}}$, and $n$ are set to the steady-state radial profiles calculated by the TASK3D simulation for the initial input ${\mathbf{u}}_{\mathrm{init}}$. The initial mean profiles of $c_{\mathrm{e}}$, $c_{\mathrm{i}}$, and $d$ are set to 1, and that of $v$ is set to 0.

3.3  Virtual LHD plasma

The virtual LHD plasma is generated in numerical space using TASK3D and controlled by ASTI, which implements the control algorithm described in the previous section. To create gaps between the digital twin and the actual system (virtual LHD), transport models different from those employed in ASTI are used in the virtual LHD plasma. In this paper, we consider two models of thermal diffusivity depending on $T_{\mathrm{e}}$, ${\chi }_{\mathrm{e}}(T_{\mathrm{e}})$, as shown in Fig. 3. Model 1 is a model in which the thermal diffusion increases linearly with the electron temperature, and model 2 is the opposite. Both models make the diffusion term nonlinear. Furthermore, the particle diffusivity in the virtual plasma is defined by $D=(2/3){\chi }_{\mathrm{e}}$. For simplicity, other conditions are aligned with those used in ASTI. When observing the virtual plasma, ASTI obtains an observation vector $\mathbf{y}$, which is created by adding white noise as the measurement error to the radial profiles of the virtual plasma. The white noise is sampled from a Gaussian distribution with a standard deviation of 3% of the simulation value.

Fig. 3.  Thermal diffusivity models employed in the virtual LHD plasma.

4.  Results

4.1  Control performance to control the virtual plasma

First of all, we show the results of the numerical experiments with ${\sigma }_Q$ = 0.15, ${\sigma }_{Rz}$ = 0.1, and $r_y$ = 0.6. We set the target state as ${\mathbf{z}}_i=(T_{\mathrm{e},\rho =0},T_{\mathrm{e},\rho =0.25},n_{\rho =0.25})$ = (2.5, 2.0, 1.5). Figure 4 shows the control results of the virtual plasma with model 1. ASTI successfully adapts to the virtual plasma by the initial adaptation phase up to 2.1 seconds. After the start of the control, oscillations in $T_{\mathrm{e}}$ and $n$ are observed up to 4 seconds in Figs. 4(a)–(c). This is caused by a change in the $T_{\mathrm{e}}$ radial profile, which changes the diffusivities and widens the gap between the digital twin and the virtual plasma. This gap reduces the control estimation performance and creates oscillations in the control inputs, as shown in Figs. 4(d) and (e). This gap is gradually narrowed by the assimilation of observations, and the virtual plasma state converges to the target state with good accuracy from 4 seconds.

Fig. 4.  Results of a numerical experiment for the virtual plasma (model 1). (a–c) Control results of $T_{\mathrm{e}}(\rho =0)$, $T_{\mathrm{e}}(\rho =0.25)$, and $n(\rho =0.25)$, respectively. The plotted values labeled “Prediction” correspond to the expected values of the predicted distribution for t ≤ 2.1 seconds and those of the u-filtered distributions for t > 2.1 seconds. The shaded areas represent a single standard deviation of the distributions. The plotted values labeled “Observation” are those obtained from the virtual plasma. (d) Estimated control input of the ECH powers ${𝒫}_1$ (labeled “axis”) and ${𝒫}_2$ (labeled “$\rho$ = 0.4”). (e) Estimated control input of the gas-puff parameter $𝒢$.

Figure 5 shows the radial profiles of $T_{\mathrm{e}}$, ${\chi }_{\mathrm{e}}$, $n$, and $D$ at 6 seconds. ASTI predicts the $T_{\mathrm{e}}$ and $n$ profiles with high accuracy through optimization of ${\chi }_{\mathrm{e}}$ and $D$. While ASTI prediction of the thermal diffusivity is in good agreement with that of the virtual plasmas as shown in Fig. 5(b), the reproducibility of the particle diffusivity at $\rho$ < 0.5 is low in Fig. 5(d). This is because the density gradient near the center is close to zero, which makes the estimation of the diffusivity difficult. This may need to be addressed for accurate state estimation but is not a major problem in terms of prediction.

Fig. 5.  Radial profiles of (a) electron temperature, (b) electron thermal diffusivity, (c) density, and (d) particle diffusivity at 6 seconds in a numerical experiment for the virtual plasma (model 1). The plotted values labeled “Prediction” represent the expected values of the u-filtered distribution. The shaded areas around the radial profiles represent a single standard deviation. The profiles labeled “Virtual LHD” represent true profiles in the virtual LHD plasma without measurement error. The lines labeled “w/o DA” in (b) and (d) are the diffusivities in ASTI without the DA optimization.

Table 4 shows the computation time required for $\mathrm{\Delta }t$ = 0.3 seconds on the vector computer in the ASTI-centered control system built on LHD (NEC SX-Aurora TSUBASA, 16VE: 128 parallel processes). This computation was performed for 384 (∼390) ensemble members; thus, one process computed 3 ensemble members (TASK3D simulations). All filter calculations (EnKF) were performed in a single process. The total computation time is 134 milliseconds, less than half the real time (300 milliseconds). Using a larger computer that can assign one ensemble member to one process can reduce the computation time required for the prediction step to approximately 40 ms. ASTI can perform real-time predictive control for control problems in the transport time scale, even with the current computer.

Table 4. Computation time required for one time step ($\mathrm{\Delta }t$ = 300 ms).

Step	computation time (ms)
Prediction	120
y-filter	6
z-filter	4
u-filter	4

Next, we show the control results for the virtual plasma with model 2. In this model, a higher electron temperature results in a lower diffusivity. The control results are shown in Fig. 6 in the same manner as Fig. 4. As observed in the case of model 1, oscillations emerge at the start of the control process. However, the virtual plasma also converges to the target state after 4 seconds by compensating the gap through the observation assimilation. Figure 7 shows the radial profiles of $T_{\mathrm{e}}$, ${\chi }_{\mathrm{e}}$, $n$, and $D$ at 6 seconds. It can be seen that the radial profiles are also predicted (reproduced) with high accuracy in model 2.

Fig. 6.  As Fig. 4, but for the virtual plasma with model 2. (a and b) Control results of $T_{\mathrm{e}}(\rho =0)$, and $n(\rho =0.25)$, respectively. (c and d) Estimated control inputs of the ECH powers and the gas-puff parameter.

Fig. 7.  As Fig. 5, but for the virtual plasma with model 2.

The control results for the target ($T_{\mathrm{e},\rho =0}$, $T_{\mathrm{e},\rho =0.25}$, $n_{\rho =0.25}$) = (2.5, 2.0, 1.0) are shown in Fig. 8, and those for the target (2.5, 2.0, 2.0) are shown in Fig. 9. Model 1 is employed for the diffusivity model of the virtual plasma in these numerical experiments. In both cases, the virtual plasma state successfully approaches the target state. The control seems to be more stable when controlling toward $n_{\rho =0.25}$ = 2.0 than toward $n_{\rho =0.25}$ = 1.0. This is because the electron temperature becomes more sensitive to changes in ECH power as the density decreases. A larger electron temperature response increases the gap between the diffusivity in ASTI and that in the virtual plasma, reducing the accuracy of the control estimate. Thus, the control of one variable can be destabilized by the control of other variables. This destabilization can be mitigated by setting either target-state noise or observation noise large, as discussed in Sec. 4.2.

Fig. 8.  Control results of (a) $T_{\mathrm{e}}(\rho =0)$ and (b) $n(\rho =0.25)$ for the target $(T_{\mathrm{e},\rho =0},T_{\mathrm{e},\rho =0.25},n_{\rho =0.25})$ = (2.5, 2.0, 1.0).

Fig. 9.  Control results of (a) $T_{\mathrm{e}}(\rho =0)$ and (b) $n(\rho =0.25)$ for the target $(T_{\mathrm{e},\rho =0},T_{\mathrm{e},\rho =0.25},n_{\rho =0.25})$ = (2.5, 2.0, 2.0).

Figure 10 is the scatter plots of the ensemble members at a z-filter step. The ensemble of the y-filtered distribution (labeled “Prediction”) is colored by the log-likelihood $p({\mathbf{z}}_i|{\mathbf{x}}_i)$. We can see that ASTI can capture regions of the control input with high likelihood and find a proper control input for each target state. The numerical experiments using standard parameter values demonstrate the validity of the DA-based control approach for the adaptive predictive control of electron temperature profile and density.

Fig. 10.  Scatter plots of the ensemble members in the z-filtering step at 5.1 s in numerical experiments for the target density of (a) 1.0 × 10¹⁹ m⁻³, (b) 1.5 × 10¹⁹ m⁻³, and (c) 2.0×10¹⁹ m⁻³. The ensembles labeled “Prediction” represent the y-filtered ensembles. Each control input to which ensemble members are assigned is colored by the median of the log-likelihoods for the target state. The ensembles labeled “Filtered” represent the z-filtered ensembles, and the x points represent the mean of the z-filtered ensemble.

4.2  Dependence of control performance on hyperparameters

In this subsection, we discuss the hyperparameter dependence of the control accuracy of ASTI using numerical experiments. The settings of the experiments are the same as in Sec. 4.1, except for ${\sigma }_{Rz}$, $r_y$, and ${\sigma }_Q$. We set the target state as ${\mathbf{z}}_i=(T_{\mathrm{e},\rho =0},T_{\mathrm{e},\rho =0.25},n_{\rho =0.25})$ = (2.5,2.0,1.5) and use model 1 as the diffusivity of the virtual plasma.

 4.2.1 Target-state noise and observation noise

The target-state noise affects the performance of the z-filter and determines the accuracy and speed with which the actual system state approaches the target state. A large value of ${\sigma }_{Rz}$ weakens the z-filter’s force to attract the prediction distribution (y-filtered distribution in this control algorithm) to the target state, which results in a gradual change in the actual system state. The observation noise determines the performance of the y-filter (adaptation). A large value of $r_y$ slows the update of model parameters and makes the adaptation robust to outlier observations. To investigate the dependence of the control performance on the target-state noise and the observation noise, we have performed numerical experiments varying ${\sigma }_{Rz}$ and $r_y$. The system noise parameter ${\sigma }_Q$ was fixed at 0.15.

Figure 11 shows the ${\sigma }_{Rz}$- and $r_y$-dependence of the control performance as root mean square percentage error (RMSPE) between the observation and target state of $T_{\mathrm{e}}(\rho =0)$. Each line is smoothed by a moving average of window size 0.2. The RMSPEs tend to increase for large ${\sigma }_{Rz}$. This trend is expected, as ${\sigma }_{Rz}$ determines the accuracy in achieving the target state. It should be noted that the trend weakens when the parameter $r_y$ is small ($r_y$ = 0.2 and 0.4). This may be because a strong y-filter reduces the uncertainty in the digital twin, making the control estimation by the z-filter more accurate. Moreover, large values of ${\sigma }_{Rz}$ slow down the rate of change of the plasma state ($T_{\mathrm{e}}$) and help keep the gap between the digital twin and the virtual plasma small. Thus, while large ${\sigma }_{Rz}$ helps the digital twin adapt to observations, it temporarily degrades the control performance.

Fig. 11.  RMSPEs between the observation and the target state of $T_{\mathrm{e}}(\rho =0)$ in numerical experiments with various ${\sigma }_{Rz}$ and $r_y$ values.

While the small ${\sigma }_{Rz}$ and small $r_y$ impose higher accuracy on the control estimation and the adaptation, they have the risk of destabilizing the control. Figure 12 shows the $r_y$-dependence of the control performance for small target-state noise (${\sigma }_{Rz}$ = 0.02, 0.06, 0.1). The RMSPEs increase significantly when both ${\sigma }_{Rz}$ and $r_y$ are small. The control results for the case of ${\sigma }_{Rz}$ = 0.06 and $r_y$ = 0.1 are shown in Fig. 13. The oscillation that occurs at the start of control is not damped until the end. The z-filter works to expand the gap between the digital twin and the actual plasma when the model parameters (e.g., $c_{\mathrm{e}}$) depend on the controlled variables ($T_{\mathrm{e}}$ and $n$ in this case). This instability arises because the z-filter, which tries to quickly get the plasma state close to the target state, and the y-filter, which tries to quickly adapt the digital twin to the actual plasma, clash and work to maintain (expand) the gap. We can prevent this control instability by making the target-state noise or the observation noise larger. In other words, we should weaken either the control estimation or the adaptation.

Fig. 12.  RMSPEs between the observation and the target state for small ${\sigma }_{Rz}$.

Fig. 13.  Results of a numerical experiment for the virtual plasma (model 1) using ${\sigma }_{Rz}$ = 0.06 and $r_y$ = 0.1. (a) Control results of $T_{\mathrm{e}}(\rho =0)$. (b) Estimated control input of the ECH powers ${𝒫}_1$ (labeled “axis”) and ${𝒫}_2$ (labeled “$\rho$ = 0.4”). (c) Estimated control input of the gas-puff parameter $𝒢$.

 4.2.2 System noise for model parameters

The system noise for the model parameters controls the uncertainty of the digital twin in ASTI and affects the control estimation and the adaptation performance. We have performed numerical experiments varying ${\sigma }_Q$ to investigate the dependence of the control performance on the system noise. The other noise parameters ${\sigma }_{Rz}$ and $r_y$ are fixed at 0.1 and 0.6, respectively.

Figure 14 shows the ${\sigma }_Q$-dependence of the control performance as the RMSPEs of $T_{\mathrm{e}}(\rho =0)$, and Fig. 15 shows a comparison of the radial profiles of the electron thermal diffusivity in ASTI and those in the virtual plasma. The RMSPE is small for 0.1 ≤ ${\sigma }_Q$ ≤ 0.4, as shown in Fig. 14. For a small ${\sigma }_Q$ case, the uncertainties considered in ASTI become small, and ASTI can not adequately track changes in the diffusivity of the virtual plasma, as shown in Fig. 15(a). The loss of control accuracy is related not only to the system noise added to the model parameters but also to the system noise added to the control inputs. When the system noise to the control inputs is significantly larger in comparison to that to the model parameters, the relationship between the model parameters (e.g., $c_{\mathrm{e}}$) and the observation variables (e.g., $T_{\mathrm{e}}$) in the ensemble becomes difficult to identify at the y-filter.

Fig. 14.  RMSPEs between the observation and the target state of $T_{\mathrm{e}}(\rho =0)$ in the numerical experiments with various ${\sigma }_Q$ values.

Fig. 15.  Radial profiles of electron thermal diffusivity at 5.1 s in the numerical experiments with (a) ${\sigma }_Q$ = 0.05, (b) ${\sigma }_Q$ = 0.15, and (c) ${\sigma }_Q$ = 0.6.

On the other hand, the control accuracy also decreases when ${\sigma }_Q$ is large, as shown in Fig. 14. It can be seen that the adaptive capacity is also reduced when the system noise is large from Fig. 15(c). The y-filter updates the latest predicted distribution by assimilating the observation to the joint distribution of the state variables at the latest time and the observation time. Therefore, large system noise added between the two timings can break the relationship of the state variables between the two timings and reduce the adaptation capacity. Furthermore, a large system noise to the model parameters can reduce the control estimation (z-filter) performance, just as a large system noise to the control input can reduce the adaptation (y-filter) performance. Figure 16 shows the scatter plots of the ensemble members at a z-filter step. In Figs. 16(a) and (b), ASTI can capture regions of the control input with high likelihood and find a proper control input. However, in Fig. 16(c), ASTI is missing the proper control input because the digital twin has too large uncertainties.

Fig. 16.  As Fig. 10, but for the numerical experiments with (a) ${\sigma }_Q$ = 0.05, (b) ${\sigma }_Q$ = 0.15, and (c) ${\sigma }_Q$ = 0.6 at 5.1 s.

These numerical experiments indicate that the system noise must be set considering the balance between the noise magnitude to the model parameters and that to the control input and the gap between the digital twin and the actual system. However, there is no need to be nervous about setting the noise parameters. The DA-based control is generally stable over a wide range of noise parameters (e.g., Fig. 14). Unless excessively large or small noise is required, the adaptation aids the control estimation, and the actual system state is expected to converge to the target state. Methods can be developed to adjust the noise intensities according to the situation and the nature of the target system. This is an issue to be addressed in the future.

5.  Conclusion

In this study, we developed a real-time predictive control system based on DA for heat and particle transport in helical fusion plasmas and clarified the fundamental properties of the DA-based control approach through numerical experiments. First, we extended the integrated simulation code TASK3D to implement a digital twin that can compute the thermal and particle transport of helical plasmas in real time. To enable real-time calculation of the particle source from neutral particles, we built a surrogate model for the AURORA code using a neural network. This surrogate model accurately reproduces the AURORA calculations and operates sufficiently faster than real time. TASK3D can now simulate the time evolution of the helical plasma, including ECH, NBI heating, and gas-puff fueling, faster than real time. A real-time computational model for pellet injection is currently under development.

To demonstrate and investigate the control performance of the DA-based control system, we conducted numerical experiments to control a virtual LHD plasma using the real-time TASK3D as a digital twin of ASTI. In these numerical experiments, the electron temperature radial profile and the density were controlled in the presence of electron temperature and electron density measurements, assuming an actual situation of the LHD plasma control. We demonstrated, for several target states and virtual plasma models, that ASTI can bring the virtual plasma state close to the target state while bridging the gap between the digital twin and the virtual plasma. It was shown that the time required for ASTI to perform real-time predictive control was less than half of the real time. This indicates that simultaneous control of temperature profile and density in LHD is already feasible with the current ASTI. Furthermore, the numerical experiments clarified the effects of target-state noise, observation noise, and system noise on control performance. The relationship between the noise parameters and the control accuracy obtained in this study is an important finding for enhancing the stability and robustness of the DA-based control.

Throughout the numerical experiments, we implemented all filters using the EnKF. When the state variables are strongly coupled through nonlinear interactions and the state distribution deviates significantly from a Gaussian distribution, other filters (e.g., the PF) should be considered for implementation. The PF can address control problems involving strong nonlinearities, such as phase transitions. Moreover, the PF is expected to suppress the oscillations observed in the numerical experiments and facilitate faster convergence of the plasma state to the target state. The DA-based control using the PF will be discussed in a future paper.

The DA-based control approach enables the harmonious integration of measurement, heating, fueling, and simulation, providing a flexible platform for digital twin control of future fusion reactors. The ASTI extended in this study has already been incorporated into the LHD control system, and experimental plans for advanced plasma control in LHD using various actuators are underway. Although ASTI has been developed for helical plasmas, we are also working to extend it for tokamak plasma control.

 Acknowledgments

This work was supported by the NIFS Collaborative Research Program (NIFS22KAPT008 and NIFS23KIPT010), the ISM Cooperative Research Program (2022-ISMCRP-2026), JSPS KAKENHI (Grant Numbers JP23K19033, JP21K13901, and JP24K00609), the ENEOS TonenGeneral Research Encouragement & Scholarship Foundation, and QST Research Collaboration for Fusion DEMO.

References

• [1]   M.  Liu et al., J. Manuf. Syst. 58, 346 (2021).
• [2]   F.  Emmert-Streib et al., PNAS nexus 3, pgae456 (2024).
• [3]   Z.  Lei et al., Nat. Commun. 14, 5604 (2023).
• [4]   Y.  Morishita et al., J. Comput. Sci. 72, 102079 (2023).
• [5]   A.  Gettelman et al., Sci. Adv. 8, eabn3488 (2022).
• [6]   M.  Li et al., Nat. Geosci. 1, 1299 (2024).
• [7]   G.  Ueno et al., Sola 3, 5 (2007).
• [8]   Y.  Morishita et al., Sci. Rep. 14, 137 (2024).
• [9]   Y.  Morishita et al., Comput. Phys. Commun. 274, 108287 (2022).
• [10]   S.  Murakami et al., Plasma Phys. Control. Fusion 57, 054009 (2015).
• [11]   A.  Sakai et al., Plasma Fusion Res. 10, 3403048 (2015).
• [12]   J.  Seo et al., Nature 626, 746 (2024).
• [13]   J.  Degrave et al., Nature 602, 414 (2022).
• [14]   J.  Kates-Harbeck et al., Nature 568, 526 (2019).
• [15]   A.  Murari et al., Nucl. Fusion 61, 036027 (2021).
• [16]   G.  Evensen, Ocean Dyn. 53, 343 (2003).
• [17]   G.  Kitagawa, J. Comput. Graph. Stat. 5, 1 (1996).
• [18]   S.  Hirshman et al., Comput. Phys. Commun. 43, 143 (1986).
• [19]   Y.  Morishita et al., J. Fusion Energy 41, 1 (2022).
• [20]   M.  Hughes and  D.  Post, J. Comput. Phys. 28, 43 (1978).
• [21]   M.  Honda, Comput. Phys. Commun. 231, 94 (2018).
• [22]   H.  Funaba et al., Sci. Rep. 12, 15112 (2022).