Representation learning algorithms for inferring machine independent latent features in pedestals in JET and AUG

The standard high-performance tokamak scenarios are based on the high-confinement mode (H-mode) operation with a self-organized transport barrier and plasma pedestal at the edge.^1,2 As the heat transport in tokamak core plasmas is generally observed to be stiff, showing strong increase with temperature gradients exceeding a critical value, the overall achievable plasma performance becomes directly dependent on the boundary condition provided by the pedestal, such as the pedestal top pressure, $p ped$⁠. For example, simulations for ITER plasmas indicate fusion power to scale as $p ped 2$⁠.³ Furthermore, in future fusion reactors, these high-performance pedestals must be compatible with the necessary power exhaust measures to guarantee sufficiently long component duty cycles.⁴ However, due to the different scalings of the pedestal and power exhaust physics from present-day fusion devices to reactor-scale facilities, a present-day fusion device cannot conclusively demonstrate reactor-relevant integration of high-performance pedestal with a power exhaust solution. Therefore, accurate and fast predictions for pedestal plasmas are needed to design and operate the future tokamak fusion reactors.

Due to the multiple relevant physical processes, and spatial and temporal scales, predicting pedestal plasma states is extremely challenging, and pedestal transport and stability are very active topics of research.^5–9 Hence, reduced models assuming MHD-constrained pedestal pressure profiles provide presently the standard approach for projecting pedestal performance between scenarios and devices.^10–13 Even though these models have been very successful in predicting pedestal pressures in many presently operating tokamaks, such a reduced model approach does increase the prediction risk when extrapolating to a new domain, such as was observed when the carbon wall at JET was changed to the ITER-like wall,² or when scaling the device size toward the reactor-scale devices and integrating radiative power exhaust to the pedestal solution. On the other hand, a more complete, high physics fidelity predictive model for pedestal and core-edge integration would be expected to be computationally prohibitively demanding for agile design optimization tasks or real-time and control applications. Advancement of deep learning algorithms for high-performance computers has opened a path to a data-driven approach to overcome these gaps in predictive capability and model throughput.^14–16 

In this work, representation learning algorithms are explored for pedestal profile measurements from the JET and ASDEX Upgrade (AUG) tokamaks.¹⁷ Representation learning algorithms aim to learn a useful representation that characterizes the training data. In the context of magnetic confinement fusion devices, a useful representation learning algorithm could be considered to compress the high-dimensional experimental observations to a lower-dimensional space that represents the actual degrees of freedom of the plasma scenarios, such that the high-dimensional observations can be reconstructed with little loss from the low-dimensional representation. More generally, state representation learning (SRL) algorithms aim to learn representations that are in low dimension, evolve in time, and are influenced by actions of an agent.¹⁷ There could be several attractive applications for such representations, including reactor and scenario design by supplementing reduced model assumptions in systems codes with fast, high-fidelity integrated scenario predictions and allowing computationally efficient optimization of feed forward scenario trajectories. Furthermore, such approaches could flexibly leverage previous and current information from the system to essentially open a pathway toward plasma scenario digital twins with plasma state awareness and enable agility for timely and informed decision, similar to Kalman filtering in control applications.¹⁸ In this work, the representation learning algorithms are based on variational autoencoders (VAE).¹⁹ Similar approaches have been previously discussed, for example, in the context of divertor detachment models,²⁰ pedestals,²¹ and disruption predictions.²² 

Building on top of the previous study by Kit et al.,²¹ the capabilities of these models to infer machine size scaling in large databases of experimentally observed pedestal plasmas in JET and AUG are investigated. One of the primary challenges of data-driven algorithms is that generally these are not expected to extrapolate beyond the data that were used for training the model. This work aims to address this challenge by investigating methods to disentangle the learned representation of the latent variable models into machine independent and dependent features. Theoretically, the machine independent features can extract the information that extrapolates between fusion devices of various scales. A related study by Kit et al.²³ investigates state representation learning algorithms that aim to capture the dynamical evolution of the plasmas by training a forward model together with the representation.

The data used to explore the deep learning models in this work originate from experimental measurements of outer midplane (OMP) electron density, n_e, and temperature, T_e, profiles. For JET, these are obtained from the high-resolution Thomson scattering (HRTS) diagnostic,²⁴ and for AUG, the profiles produced by the integrated data analysis (IDA) were used.²⁵ The HRTS measurements are used for JET without including profile fits and dedicated instrument function deconvolution procedures, such that a certain level of profile smearing is present in the data.²⁶ Only ITER-like wall (ILW) JET pulse numbers (JPN) are included in this work, which are all obtained after the improvement of the optical design of the laser input system. Following the improvement of the optical design, the profile smearing caused by the instrument function was reduced to a level of a relatively small correction rather than substantially changing the profile gradient and width.²⁶ The individual profile measurements are assumed to have negligibly small measurement error relative to the plasma fluctuation between individual time slices. This is expected to be a well-justified assumption for the plasmas in this study.

It is also quite usual in the HRTS measurements to see relatively high scatter in T_e values in the scrape-off layer (SOL) as measured n_e values are low. To avoid feeding these large scatter T_e values to the deep learning algorithms, a heuristic approach was taken such that for any HRTS measurements, where n_e is lower than $5 × 10 18$ m⁻³ and T_e is higher than 50 eV, the T_e value is replaced by $T e , corr = ( n e × 10 − 18 m − 3 ) × 10$ eV. Such a procedure biases the data at low T_e values in the SOL. However, for the time being this was considered a better approach than confusing the model with a large scatter ranging between a few eV and values beyond keV, when it is known from SOL physics that the expected values are significantly below keV. In future, the plan is to improve these SOL corrections with information from SOL diagnostics and boundary models.

It is also well known that the magnetic equilibrium reconstruction is often not sufficiently accurate to locate the separatrix within the profile sufficiently well considering the profile gradient scale lengths. To obtain a common separatrix alignment procedure for the investigated plasmas, the separatrix is shifted to the location where T_e equals 100 eV for all plasmas for both JET and AUG, which is of the order of the expected values for conduction-limited SOL transport in H-mode conditions. However, this is an oversimplification and will be reevaluated in future studies, as all dependencies on power densities, SOL conditions, and device sizes are neglected by this approach. Nevertheless, such a procedure does align the separatrix approximately to the location where the radially steep gradients of the pedestal T_e begin and provides a pragmatic approach to align 10⁴ profiles with a batch algorithm. Reaching such level of robustness with a more sophisticated separatrix alignment algorithms is expected to require substantially more work than was allocated here. Furthermore, varying $T e , SEP$ within 20%–30%, as observed in SOL simulations of deuterium and nitrogen injection scans at JET,²⁷ does not affect the determination of the key pedestal parameters significantly.²⁸ It is acknowledged that by doing this, the original uncertainty in the equilibrium reconstruction is replaced with the physics uncertainty of the separatrix realignment procedure. However, the latter is considered more attractive here as it provides a controlled bias to the data that is consistent between the data points, whereas the equilibrium reconstruction uncertainty magnitude and direction is expected to be more uncertain between the data points when combining thousands of measured profiles from more than one fusion devices. All profiles were mapped to a common normalized poloidal flux, $Ψ N$⁠, grid of 30 points uniformly spaced between 0.85 and 1.05.

Only deuterium plasmas without significant impurity seeding, resonant magnetic perturbations, kicks, or pellets were considered in this study. For JET, the plasma discharge numbers are selected from the JET pedestal database,² whereas for AUG, the algorithm selects plasma discharges that have been flagged as useful in the operating journal, which in practice only filters out plasma discharges that failed. Furthermore, only plasma time slices with core T_e exceeding 2 keV are selected. The total resulting number of plasma discharges in the set of data used for training the deep learning algorithms is 1280 for JET and 3724 for AUG. From these plasmas, time periods where the total heating power exceeds 3 MW for JET and 2 MW for AUG were selected for analysis. To reduce the size of the dataset, the total database of measured profiles was downsampled by randomly selecting ten time slices from each investigated time window as well as choosing the three time slices with the highest pedestal pressure within each time window. The intention of the former approach is to approximately retain the statistical distribution of the full data while downsampling, whereas the intention of the latter is to make sure that the highest-pressure profiles right before edge localized modes (ELMs) are retained in the training database. The downsampling procedure also helps to balance the dataset as there are significantly more profiles in a given AUG time window relative to a given JET time window. Since there is no effort to filter ELMs or identify ELM cycles in this study, the resulting scatter in the profile data is relatively large and encoded as stochastic variation by the model in this study. In future studies, methods to identify fractions of ELM cycles automatically from the data will be investigated. The dataset reduction is expected to help with agile algorithm testing. The total resulting number of plasma time slices is 55 087, as the average number of time slices retained for a plasma discharge is a bit more than 11.

The following input machine control parameters were selected for the model in this work (Table I): (1) toroidal magnetic field, $B T$ (2) safety factor at 95% flux surface, q₉₅, (3) total heating power normalized by the Martin LH-threshold scaling,^29, $P TOT / P LH$⁠, (4) total deuterium injection rate, $D TOT$⁠, (5) major radius, R, (6) aspect ratio, $A = R / a$⁠, where a is the minor radius, (7) elongation, κ, (8) upper triangularity, δ^u, and (9) lower triangularity, δ^l, of the plasma. Even though the Martin LH-threshold scaling might not provide a fully accurate LH-threshold value for all plasmas considered in this study, the scaling provides a pragmatic approach to normalize the total heating power as two devices of different scales are considered. As can be seen in Table I, the average plasma in the investigated database is an H-mode at about 2.4 T with a moderately low q₉₅ of 3.5–5.0 and moderately low upper triangularity of 0.1–0.2. Figure 1 illustrates the distributions of the n_e and T_e values near the pedestal top in this dataset. The dataset used for training the model is split into training, validation, and test sets in proportions of 70%, 20%, and 10%. The data are further normalized by subtracting the mean and dividing by the standard deviation of the training set. In addition to this, two representative plasmas that have been documented in previous publications^13,30 are specifically selected to be retained outside the dataset used for training the model and applied for prediction testing. In the remainder of this manuscript, these will be called the holdout plasmas. For AUG, the plasma discharge number 33 173 was selected, which is further discussed in Ref. 13. For JET, the plasma discharge number 96 202 was selected, which is further discussed in Ref. 30. The representative input machine control parameters for the two holdout plasma discharges are shown in Table II. These plasmas represent relatively standard operational points and, therefore, test the model capability in interpolating between previously observed points within the data distribution.

Algorithms to learn a representation with tokamak size regression are explored. All approaches are based on the variational autoencoder (VAE),¹⁹ which encodes the representation for the observed plasma profile information in this work (Fig. 2). In Sec. III A, a standard VAE is first trained with the dataset to demonstrate that the representation learning model has sufficient capacity to reconstruct the profiles well. Once this is demonstrated, Sec. III B explores an approach similar to domain invariant variational autoencoder (DIVA)³¹ to connect the machine control configurations with the learned representation, building on top of the previous work by Kit et al.²¹ However, as will be discussed in Sec. III B, when connecting the device size together with the machine control parameters to the same auxiliary learning objective, the resulting device size regression is mixed with changes of the other parameters. This leads to large reconstruction error for the device size, indicating that the device size-dependent features are not well disentangled from the size independent features. Therefore, in Sec. III C, a separate latent space is dedicated to the device size regression, which is observed to facilitate disentangled encoding of the size-dependent features into this latent variable.

A VAE is a generative model, defining a joint distribution over the observed data, $x$⁠, and latent variables, $z$⁠, as $p θ ( x , z ) = p ( z ) p θ ( x | z )$⁠. Learning representation of the data means to learn the posterior distribution $p θ 1 ( z | x )$ that maps the observations to distributions of latent variables. The generative side of the model aims to model the dataset, $x$⁠, as a conditional distribution given the latent representation, $z , p θ 2 ( x | z )$⁠. In this work, this is implemented using encoder, $q Φ ( z | x )$⁠, and decoder, $p θ 2 ( x | z )$⁠, distribution, parameterized by convolutional neural networks (Fig. 2). The $q Φ ( z | x )$ represents an amortized variational approximation of the intractable true posterior, $p θ 1 ( z | x )$⁠. The VAE model architecture in this work is shown in Table III. The encoder consists of convolutional layers that take the radial n_e and T_e profiles as input. The output from the second convolutional layer is transformed into mean and standard deviation of 12 latent variables through two dense layers. The decoder consists of transposed convolutional layers. The last transposed convolution layer provides two output channels, which are transformed into two radial profiles of 30 points each through a dense layer.

The second term constraints the model to learn a latent representation that resembles a Gaussian hypersphere. This constraint aims to regularize the learned representation toward continuity. The $E q Φ ( z | x )$ stands for expectation of the loss integrated over the amortized variational approximations of the posterior. In practice, the reconstruction error is obtained by sampling a reconstruction from the VAE latent posterior distribution and computing the mean squared error. $β prior$ multiplier is used to control the strength of the KL constraint, as motivated by the β-VAE by Higgins et al.³² Similarly, $α obs .$ multiplier is added in front of the reconstruction error to also adjust the strength of the loss term. In addition to the standard objective in Eq. (1), a physics regularization term is added to the objective function $∝ β press | | n e T e − n e ̂ T e ̂ | | 2$ to add a penalty if the reconstructed static pressure, $n e ̂ T e ̂$⁠, deviates from the measured static pressure, $n e T e$⁠. The intention is that this term encourages reconstruction of physically mutually consistent n_e and T_e profiles. Such a model was trained for 100 epochs with the adam optimizer³³ using a learning rate of 0.01 and batch size of 1024. The loss multipliers are shown in Table IV. These hyperparameters were manually selected based on manual testing of various hyperparameter values. In future, a more careful hyperparameter optimization studies will be performed. Figure 3 illustrates the reconstruction quality for the holdout plasmas and their representative locations within a 2D cut of the learned 12-dimensional latent representation. The overall reconstruction errors are very small throughout the profiles, which indicates that the present model setup has the necessary capacity to represent the profiles in the investigated dataset (Table V).

Since there are no other learning objectives than the reconstruction error and the overall shape of the latent representation through the KL constraint, the information may be entangled between the latent variables in this learned representation. To generate a useful representation, a desired property is organization of the information into disentangled latent variables, enabling interpretation of the latent variables as representations of certain, semantically meaningful characteristics of the dataset.^31,34–37 To achieve disentangled latent variables, semi-supervision of the learning task is explored.^31,35,37 By introducing auxiliary learning objectives to a selection of the latent variables, the learning algorithm is given an incentive to encode the semantic information related to the learning objective to those latent variables. Attaching the auxiliary learning objective to the machine parameter configuration of the plasma discharge, a connection is also established between the machine configuration and the state of the plasma. The first approach, discussed in Sec. III B, follows broadly the structure of the domain invariant variational autoencoder (DIVA).³¹ In Sec. III C, this approach is extended by separating a dedicated latent space to represent the scaling with the machine size. The latter approach is actually closer to the original DIVA publication, where the machine size can be considered to represent a class label and the machine parameter configuration a domain label. However, the machine size in this work is treated as a continuous rather than discrete variable.

The original DIVA publication by Ilse et al. discusses splitting the latent space of a VAE to three sub-spaces: (1) domain, (2) class, and (3) residual variations.³¹ The motivation is to disentangle the latent representation to domain variant and invariant features, such that the domain invariant features could potentially be applicable to new, previously unseen domains. In this work, this approach is applied for JET and AUG experimental data with the aim to learn a representation for which certain latent features are independent of the device size. Fundamentally, the idea is that if it is possible to disentangle the features that are dependent and independent of the device size, these tools could be used to infer fundamental device size scalings in the experimental databases and potentially also to build machine learning models that are transferable to new device sizes with relatively small amounts of training data at scale.

This learning objective aims to encode information from the machine parameter configuration into the latent variable, $z y$⁠, together with the observed n_e and T_e profile information. As the neural network generated prior distribution depends conditionally on the machine parameter configuration, the algorithm is given incentive to encode such a latent representation that also conditionally depend on the machine parameter configuration. Furthermore, the machine parameter reconstruction objective also incentives the algorithm to encode such a latent representation that the information content is sufficient to reconstruct the machine parameter configuration. Therefore, correlations with plasma state and machine parameters would be expected to be preferentially encoded to $z mp$⁠, while any residual variations would be expected to be encoded to the residual stochastic part of the latent space, $z s$⁠. To retain the same latent space dimensionality as in the initial VAE example, $z mp = 9$ and $z s = 3$ were selected in this study. For the original VAE model, the only change is that for those nine latent dimensions, the prior distribution is given primarily by the neural network model rather than the static $N ( 0 , 1 )$ prior. For the remaining three latent dimension, the prior is simply the static $N ( 0 , 1 )$⁠. The chosen latent space dimensionality is small relative to the dimensionality of the profiles, but not relative to the dimensionality of the machine parameters. This is intended, as the machine parameters are supposed to provide the independent control parameters of the system. If this would be exactly the case, then one would not be able to compress the information any more than is the dimensionality of the control parameter vector. In practice, this is not the necessarily the case as the experimental databases have cross-correlations within them.

The loss multipliers are shown in Table VII. The $β prior mp$ stands for a small static prior loss on $N ( 0 , 1 )$ applied for the machine parameter latent space. This is a regularizing term in this model setup. The multiplier value is set low as this loss should not dominate over the prior regression loss for the model to work as intended. Overall, these hyperparameters resulted from simple manual scan and represent values that lead to reasonable model performance. A more systematic hyperparameter optimization would probably improve the model performance further but was not considered in this study.

Figure 4 illustrates prediction performance of the trained model for the holdout plasmas. By providing the machine control parameters of the holdout plasmas through the prior regression module, a conditional distribution is obtained in the latent representation. This conditional distribution is showed with stars and error bars in Fig. 4. From this conditional distribution, full profile samples can be generated with the decoder of the VAE. As can be seen, the generated profiles are broadly in line with the observations, although the detailed profile shapes are not necessarily fully reproduced and the prediction uncertainties are quite large. One of the leading reasons for both of these shortcomings of the model is presumably the fact that ELM cycles are averaged by the model in this work. Therefore, stochastic variation below the ELM fluctuations would not be expected through the generative model with this model setup. In future studies, focus will be given to more directly encode ELM information into the latent representation, which is expected to improve the prediction quality.

For a more quantitative performance test, mean absolute errors for reconstructions and conditional generation or prediction for the test set were calculated (Tables VIII and IX). There is no noticeable loss of reconstruction accuracy in the VAE model even though some of the VAE capacity is allocated for the additional training objectives (Tables V and VIII). The conditional generation or prediction error is about a factor of 5–10 larger than the reconstruction error (Table VIII), which is expected as solely based on the machine parameters, there is no information about the fraction of the ELM cycle at which the target profile is measured. Considering this challenge, the prediction error is also relatively small. Actually, compared to Fig. 4, it can be seen that the prediction error is of the order of the standard deviation of the predicted distribution, which is expected. Therefore, when predicting an individual time slice with this model, one would not expect the model to be more accurate than this. Same arguments apply for the reconstructions of the machine parameters (Table IX). In future, further development of the model to encode information about the ELM cycle is expected to reduce both the prediction error and the machine parameter reconstruction error.

As a test of the capability of the generative model to predict the impact of R on the pedestal n_e and T_e, the following experiment was conducted. The machine parameter inputs were set to representative values near the mean of the training set: $| B T | = 2.25$ T, $q 95 = 4.7$⁠, $P TOT / P LH = 2.4 , D TOT = 10 22$ e/s, A = 3.2, $κ = 1.67$⁠, $δ u = 0.14$⁠, $δ l = 0.4$⁠. Then, R was scanned from 1.6 to 2.9 m (Fig. 5). It is observed that the model does encode a scaling of reduced pedestal density as R is increased, which is consistent with the scaling of the empirical Greenwald density limit in tokamaks,³⁸  $n GW = I p π a 2 ∝ B T R$⁠, where the latter relation assumes constant plasma shape and q₉₅ when scaling the device size. However, as can be observed from the reconstructed R, the internal representation of the model is not fully consistent with the scan of the major radius (Fig. 5). As there are more AUG time slices than JET time slices in the training database, it seems that the average representation is weighted toward smaller device sizes. As in this model construction, all the machine parameters are encoded and reconstructed through a common fully connected neural network, the machine parameters end up being entangled with each other. Therefore, even though the model does learn the correlation of reduction of average pedestal top n_e with device size, this correlation is mixed with changes of other parameters, such that the internal state of the model is entangled. As a result, when the machine parameters are reconstructed from latent distribution, the machine parameter reconstruction model mixes the original drivers of the posterior latent distribution, leading to large reconstruction error on the device size. Therefore, although the results are already quite promising in terms of connecting machine parameters to the learned representation, the information organization is not yet quite as disentangled as would be needed to fully disentangle the machine size dependence from the other dependencies in the dataset. Therefore, in Sec. III C, a separate latent space is dedicated for the machine size.

A dedicated latent space of dimensionality of one was allocated for the machine size regression model. To retain the same latent space dimensionality as previously, the latent space dimensions were set as $z mp = 8 , z size = 1$⁠, and $z s = 3$⁠. The machine parameter and size regression architectures are similar to that presented in Table VI. The primary difference is that the input, output, and latent space dimensionality of the machine parameter regression are reduced to eight, and the representative dimensionality for the size regression is one. The loss is similar to Eq. (2) with the machine size reconstruction loss and KL constraint separated from the other machine parameters. The loss multipliers are shown in Table X. As before, these hyperparameters are obtained through a manual search and a more systematic hyperparameter optimization could potentially improve the model performance further.

A clear qualitative feature that can be seen with the model setup is that the R scaling is now well organized into latent dimension eight (Fig. 6). As the latent dimension eight is scanned from 4 to −4, the inferred tokamak major radius increases from mid-sized 1.5 m to large 3.0 m. Also, the Greenwald like density scaling can now clearly be seen when comparing the upper and lower latent space contour plots in Fig. 6, where the smaller device size is associated with higher pedestal n_e values.

As the model becomes more regularized through splitting of the latent space, there is a small, less than factor of 2, increase in the VAE reconstruction error (Tables VIII and XI). The prediction errors remain almost the same between the models with and without the separate latent space for device size (Tables VIII and XI). For most of the machine parameters, reconstruction errors remain the same with the exception that the reconstruction error on R is reduced by a factor of 2, as might be expected (Tables IX and XII).

Finally, the test on the capability of the generative model to predict dependence of pedestal n_e and T_e on machine size shows that the internal representation now captures the size dependence significantly better (Figs. 5 and 6). As the input R approaches JET sized values, the internal state also approaches these numbers, which was not the case when R was mixed in the same latent space with the other machine parameters (Figs. 5 and 6). Furthermore, similar tests were done by scanning the size from the latent representation of each of the holdout plasma discharges. When scanning from the AUG case, the dependence was very similar to that observed in Fig. 7. From the JET case, the dependence was significantly weaker, which indicates that further research is needed to fully understand how the size dependence emerges globally in the latent space.

Variational autoencoder-based representation learning algorithms are explored for datasets of thousands of experimentally observed pedestal n_e and T_e profiles at JET and AUG tokamaks. Representation learning aims to establish a useful representation that characterizes the dataset. In the context of magnetic confinement fusion devices, a useful representation could be considered to map the high-dimensional observations to a lower-dimensional manifold that represents the actual degrees of freedom of the plasma scenarios. In this work, the capability of these algorithms to infer tokamak size dependence in multimachine databases is investigated.

Using VAE-based representation learning models and following broadly the architecture of DIVA to connect the machine parameter configuration with the learned representation, the model learns joint probability distributions that associate a certain device size and machine control configuration with a certain pedestal n_e and T_e. Since ELM fluctuations are treated as stochastic variations of the profiles and not taken into account in the model design, the prediction variances for profiles and machine parameters remain quite large. Given this shortcoming, the performances for predicting pedestal n_e and T_e profiles based on the machine control configuration and device size within the test sets are quite good.

It was observed that when connecting the device size together with the machine control parameters to a common regression objective and latent space, the model ended up entangling the machine parameters within each other. As a result, conditional latent distributions with a certain device size were mixed with other device size and machine parameter configurations. By separating the device size to a dedicated regression objective and latent space, similar to separation of domain and class labels in the original DIVA publication,³¹ the device size information becomes well organized and disentangled from the other machine parameters. This highlights the importance of properly designing the learning algorithm structure appropriately to appropriately encode the information into semantically meaningful latent variables.

For these multimachine datasets, the representation does encode density scaling with a device size that is qualitatively consistent with Greenwald density limit scaling.³⁸ Since the training dataset consists of data from AUG and JET, it should be noted that there are no training observations at values in between the two extremes. Therefore, for a good fraction of the scan in Fig. 7, the model is interpolating in a domain where there are no training data, and this interpolation seems quite stable, which is not necessarily given when operating with overparameterized models capable for overfitting. The features that the model learns to be independent of the device size could in principle extrapolate beyond the training domain. In future studies, such interpolation and extrapolation performance could be explored by training the model on JET and Tokamak à Configuration Variable (TCV), for example, and testing the interpolation performance for AUG, or extrapolation performance by training on TCV and AUG and extrapolating to JET. How much real data are needed at scale and whether simulations can be used to supplement the real data at scale are interesting research questions to be addressed in future.

This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No. 101052200—EUROfusion). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them.

The work of Aaro Järvinen and Adam Kit was partially supported by the Research Council of Finland Grant No. 355460.

The authors wish to acknowledge CSC-IT Center for Science, Finland, for computational resources. The relevant CSC Project No. 2005083.

The authors have no conflicts to disclose.

Aaro E. Järvinen: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Funding acquisition (equal); Investigation (lead); Methodology (lead); Project administration (equal); Resources (equal); Software (lead); Supervision (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (lead). Adam Kit: Conceptualization (equal); Data curation (equal); Formal analysis (supporting); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (supporting); Writing – original draft (supporting); Writing – review & editing (supporting). Yoeri Poels: Formal analysis (supporting); Methodology (supporting); Software (supporting); Writing – original draft (supporting); Writing – review & editing (supporting). Sven Wiesen: Funding acquisition (lead); Project administration (lead); Resources (lead); Writing – original draft (supporting); Writing – review & editing (supporting). Vlado Menkovski: Conceptualization (supporting); Formal analysis (supporting); Methodology (supporting); Writing – original draft (supporting); Writing – review & editing (supporting). Lorenzo Frassinetti: Data curation (supporting); Writing – original draft (supporting); Writing – review & editing (supporting). Mike Dunne: Data curation (supporting); Writing – original draft (supporting); Writing – review & editing (supporting).

The experimental data used for training the deep learning models in this work are stored at the data storage facilities of AUG and JET, and the authors do not have permission to make these data publicly available. However, the training data preparation routines will be provided on request such that anyone with access to the data can regenerate the training dataset. The Python code encompassing the deep learning algorithms is available in GitHub at https://github.com/DIGIfusion/latent-state-modeling (Ref. 39).