Multi-Timescale Dynamics Model Bayesian Optimization for Plasma Stabilization in Tokamaks
Authors: Rohit Sonker
Alexandre Capone
Andrew Rothstein
Hiro Josep Farre Kaga
Egemen Kolemen
Jeff Schneider
Venue: ICML 2025 poster (2025)
Keywords: Nuclear Fusion
Plasma Instabilities
Bayesian Optimization
Applied Machine Learning
Primary area: applications->chemistry_physics_and_earth_sciences
TL;DR: Multi Scale Bayesian Optimization which avoids instabilities in nuclear fusion with a success rate of 50% - marking a 117% improvement over past experiments

Abstract: Machine learning algorithms often struggle to control complex real-world systems. In the case of nuclear fusion, these challenges are exacerbated, as the dynamics are notoriously complex, data is poor, hardware is subject to failures, and experiments often affect dynamics beyond the experiment's duration. Existing tools like reinforcement learning, supervised learning, and Bayesian optimization address some of these challenges but fail to provide a comprehensive solution. To overcome these limitations, we present a multi-scale Bayesian optimization approach that integrates a high-frequency data-driven dynamics model with a low-frequency Gaussian process. By updating the Gaussian process between experiments, the method rapidly adapts to new data, refining the predictions of the less reliable dynamical model. We validate our approach by controlling tearing instabilities in the DIII-D nuclear fusion plant. Offline testing on historical data shows that our method significantly outperforms several baselines. Results on live experiments on the DIII-D tokamak, conducted under high-performance plasma scenarios prone to instabilities, shows a 50\% success rate — marking a 117\% improvement over historical outcomes.

Multi-Timeecaal Dynamics Model Bayesian Optimization for Plasma Stabilization in Tokamaks
Rohit Sonker
Alexandre Capone"
Andrew Rothstein Jeff Schneider
Hiro Josep Farre Kaga² Egemen Kolemen' 23
Abstract
are especially pronounced when controlling tokamaks, a form of nuclear fusion reactor. For tokamaks, good models are unavailable, safety paramount, and instabilities are notoriously hard control. These issues are further compli- cated by the fact that the dynamics may fluctuate strongly due to various reasons: tokamaks experience dozens of ex- periments any given day. some of which deliberately change its base configurations, e.g., by introducing impu- rities into the plasma, potentially affecting the dynamics of posterior experiments. Moreover, tokamaks frequently undergo hardware changes or outages, e.g., wall repairs or beam failures, making any single specific model unreliable. However, despite these challenges, designing good control policies for tokamaks is highly desirable due to their promise to generate abundant clean energy via nuclear fusion.
Machine learning algorithms often struggle to con- trol complex real-world systems. In the case of nuclear fusion, these challenges are exacerbated, as the dynamics are notoriously complex, data is poor, hardware is subject to failures, and exper- iments often affect dynamics beyond the experi- ment's duration. Existing tools like reinforcement learning, supervised learning, and Bayesian op- timization address some of these challenges but fail to provide comprehensive solution. To over- come these limitations, we present multi-scale Bayesian optimization approach that integrates a high-frequency data-driven dynamics model with a now-frequency Gaussian process. By updating the Gaussian process between experiments, the method rapidly adapts new data, refining the predictions of the less reliable dynamical model. We validate approach by controlling tearing instabilities in the DIII-D nuclear fusion plant. Offline testing on historical data shows that our method significantly outperforms several base- lines. Results on live experiments on the DIII- D tokamak, conducted under performefore plasma scenarios prone to instabilities, shows 50% success rate marking 117% improve- ment over historical outcomes.
In many real-world control applications, model-free rein- forcement learning promising solution and has seen successful applications (He al., 2024; Kumar et al., 2021; Lee et al., 2020). However, most of these methods rely on prohibitive amount of policy rollouts for training, which is typically only achievable with reliable simulation envi- ronments. In complex environments like tokamaks, this is particularly problematic, as operation costs typically only permit a handful of rollouts, and existing simulators do not reflect the true dynamics for many aspects of the plasma (Char al., 2023a). Offline RL seeks to overcome these issues by learning policy from offline data that conserva- tively stays within the bounds of the observed data (Levine et al., 2020). However, the performance of offline RL meth- ods depends crucially on high-quality expert data that con- tains advantageous states. these are not present, then offline RL can suffer from extrapolation errors (Fujimoto et al., 2019). This is major drawback for tokamak con- trol, where significant exploration and improvement are still required to achieve energy production.
1. Introduction
Controlling real-world systems is inherently challenging, even when powerful machine-learning tools are employed: nonlinearities are often pronounced, data is scarce, and safety issues impose severe limitations. These challenges
A Alternatively model-based reinforcement learning offers a solution where dynamics models are trained from historic data and rollouts from the model are then used for policy learning or plamning (Deisenroth & Rasmussen, 2011: Chua et al., 2018; Osinski et al., 2020). In the past, machine learn- ing algorithms have been used to directly model plasma dynamics (Char et al., 2023b: Abbate et al., 2021; Boyer
'Equal contribution 'Robotics Institute, Carnegie Mellon Uni- versity, Pittsburgh, PA, USA PPrnnceton University. Princeton, NJ, USA PPrinceton Plasma Physics Laboratory, Princeton, NJ, US. Correspondence to: Rohit Sonker <rsonkerr unndrw.ccm.edu..
Proceedings of the ndnd International Conference Machine Learning, Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

Multi-Timescale Dynamics Model Bayesian Optimization for Plasma Stabilization in Tokamaks
al., 2021). Reinforcement learning policies have also been trained in models trained solely on fusion data (Char et al., 2023a; Wakatsuki et al., 2023; Degrave et al., 2022). However, the performance of these approaches crucially hinges on the assumption that the data faithfully captures the model at test time. This is problematic in the case of tokamak dynamics, where time-dependent model changes cannot be neglected. Though this issue can be potentially ad- dressed by updating the model with new data, the scareity of experiments implies that too little data is typically produced to reliably update the model.
temperature and pressures for increased energy, hence, they are more unstable. Our method was able to find stabilizing ECH actuator values in four of eight experiments despite changes to other actuators, 1 117% improvement compared to historical experiments with the same configuration.
Our paper is structured as follows: first, we provide some necessary background to nuclear fusion, and define our prob- lem mathematically. Then we discuss our complete pipeline and methodology, followed by the results and analysis on offline historical data and live experiments on a Tokamak reactor. Finally, we provide conclusions and discuss oppor- tunities for future work. Additional details are provided in the Appendix.
In low-dimensional settings, the obstacles posed by conven- tional RL methods can potentially be addressed by Bayesian optimization (BO). BO is a data-efficient tool for optimizing black box functions (Garnett, 2023). By quantifying model uncertainty, BO achieves tradeoff between exploration and exploitation, leading to fast convergence in many practical settings (Shalloo et al., 2020; Shields et al., 2021). In the case of tokamak control, BO has been used, e.g., to control the rampdown of a real tokamak (Mehta et al., 2024), and to control neutral beams in a tokamak simulator (Char et al., 2019). However, the work of (Mehta et al., 2024) does not address critical plasma instabilities, whereas (Char et al., 2019) relies on simulator. Moreover, these methods use poorly specified prior and require an extensive amount of experiments to perform well.
2. Background and Problem Statement
In this section, we first provide some background on nuclear fusion and then present the formal problem statement. We also include more details in the Appendix section.
2.1. Nuclear Fusion
Nuclear fusion is seen as promising solution for clean, limitless energy, producing no high-level radioactive waste. Among the fusion technologies, tokamaks are the most ad- vanced, using magnetic fields to confine hot plasma to en- able fusion conditions. Many countries have invested in tokamak research facilities and currently more 35 nations are collaborating to build ITER, a global project aiming to demonstrate the viability of large-scale commercial fusion reactors (Mohamed et al., 2024; Shimada et al., 2007).
Motivated by the strengths and shortcomings of existing machine barning-based approaches for tokamak control, we design novel approach that combines a dynamic model predictor and Bayesian Optimization. Our approach em- ploys a multi-scale approach: a recurrent probabilistic neu- ral network models the high-frequency model dynamics, while a Gaussian process models the effect of low-frequency marginal statistics on the dynamics. After adequate pre- processing, we use historical data to train both models, where the dynamic model serves as a prior for the Gaus- sian process. Additionally, by leveraging physics-informet assumptions, we design a low-dimensional state space for the Gaussian process. This naturally leads to a contex- tual Bayesian optimization algorithm tailored to the task at hand, allowing it to find stabilizing actions in a highly data-efficient manner. Moreover, due to its ability to per- form fast updates, allows us to efficiently leverage small batches of data collected during experiments best inform new decisions on the fly.
One of the key challenges in tokamak development is plasma disruptions, which can cause severe damage to reactor walls and components, particularly in larger reactors like ITER (Schuller, 1995; Lehnen et al., 2015). These disruptions of- ten stem from tearing mode instabilities (or tearing modes), where magnetic islands form, leading to energy loss and instability. Prior work proposes avoiding tearing instability with predictive models using real-time control (Fu et al., 2020) and reinforcement learning (Seo et al., 2024). How- ever, these methods reduce neutral beam power and add torque to stabilize the plasma. This is undesirable, as reduc- ing beam power leads to lower confinement energy, decreas- ing the total energy output of the tokamak. On the other hand, adding torque to large tokamaks is itself a challenging issue.
We test our approach on large dataset from past tokamak experiments, where we can quickly identify stable configura- tions, outperforming a naive approach based exclusively on the recurrent neural network model. Furthermore, we apply our approach in a real-world on the DIII-D Tokamak to find stabilizing actions for high performing plasma scenario. High performing plasma scenarios need to maintain high
Electron Cyclotron Heating (ECH) has shown promise in counteracting tearing instabilities by driving localized cur- rents at the site of instability (Gantenbein et al., 2000; Kole- men et al., 2014). These and other findings have motivated the inclusion of gyrotrons capable of delivering ECH in fu- ture reactors to potentially control tearing instabilities, e.g., ITER will have over 40 gyrotrons. So far, the best results
2

---

Multi-Timescale Dynamics Model Bayesian Optimization for Plasma Stabilization in Tokamaks
for stabilizing instabilities with ECH have been achieved by keeping the ECH constant over time, as this minimizes the chance of plasma disruptions. However, how to best deploy ECH is still an open question, as researchers often struggle to find ECH profiles to suit their goals. An ECH profile represents the heating achieved by the gyrotrons across the cross section of the plasma. This can be seen in Fig 4.
og corresponds to gyrotron angles, operated at constant power, which we use to keep the tearing instability from occurring. The gyrotrons operate on the plasma by gener- ating an ECH profile af = o(a.). Unlike af and of, the number of gyrotrons, i.e., the dimension of af potentially changes between each individual experiment. This is due to various reasons, e.g., due to hardware issues or because some gyrotrons might be required for other tasks, such as elm suppression or density control (Hu et al., 2024; Ono et al., 2024).
In this work, we aim to find feedforward ECH profiles to avoid tearing instability (or modes) in high Jmin tokamak scenarios. EC waves from gyrotrons lead to EC heating (ECH) profiles and EC current drive (ECCD). The amount of current drive depends on angle of incidence on the plasma. In all experiments, we set angle to produce both ECH and ECCD. More details on effect of ECH on plasma are given in A. Pre-emptive suppression of tearing modes with ECCD has been previously shown in (Bardóczi et al., 2023), where ECCD was used at q=2 location. Our goal in this work is to optimize the size and locations of ECH profiles, conditioned on the state of the plasma. During our experiments we focus on the High Jmin plasma scenario. This is a scenario that supports long duration steady-state plasma operations, mak- ing crucial for future commercial fusion reactors. More details on High gm scenario are provided in the Appendix. We also focus our attention on 2-1 tearing instability, a type which is the most common and significantly disruptive.
This paper considers the case where the gyrotron angles of are kept fixed throughout each experiment rollout, i.e., of = a = = at =: ad, where T is the length of the rollout horizon. This a common operating mode and also a design choice, which we make because we need to search as efficiently as possible within the action space, impossible task if its dimension is too large. Thus, the feedforward actions of and the target normalized plasma pressure 3y. which defines the set-point for af, are specified beforehand and can change between rollouts. Our goal, is to then select a separately for each experiment such that the probability of encountering a tearing mode Té 1 is minimized over the full rollout horizon.
3. Methodology
2.2. Problem Statement
We now introduce our method, DynaBO, which aims to find stationary feedforward ECH profiles that mitigate tearing instabilities. Our complete pipeline is illustrated in Fig. 1. On high level, the process is as follows We model the system at two different time scales to inform the choice of actuator commands for each experiment. At smaller, more granular time scale. we use a recurrent probabilistic neu- ral network model (RPNN) to estimate the high-frequency behavior during each experiment. The coarser model corre- sponds to a Gaussian process model and is trained to predict the behavior of the system based on marginal statistics from experimental observations and RPNN predictions of the ob- jective function, which act as prior mean. In this case, the objective function is the time-to-tearing instability. Given the target normalized plasma pressure BN, we leverage the Gaussian process to select actions (ECH profiles) in low- dimensional space via Bayesian optimization. The desired profile is then converted gyrotron angles and applied to the tokamak. Finally, we update our model with the result- ing time-to-tearing instability and actual ECH profile and repeat the procedure. We update the model with the mea- sured ECH profile because it can diverge significantly from the desired one. In the following sections, we discuss the individual components of our method - the high- frequuryy RPNN and instability predictor, the Gaussian Process model, and action selection by Bayesian optimization.
We treat the tokamak dynamics an unknown discrete-timg stochastic system
81+1 ST
with states 8, E S and actions on E A. and the probability of tearing mode occurring follows a Bernoulli distribution, parameterized by the tokamak states and actions
T- ~ Bernoulli(p(s, a:)).
(2)
Of the state variables describing the plasma, the most im- portant for our approach is the normalized plasma pressure 3N E 8,. A full description of the state space is given in the appendix. The action vector can be decomposed into three different sub-vectors
at = [of,
[ af af
(3)
as follows. The actions af correspond to feedforward inputs specified before the experiment. These correspond, e.g., to gas flows, plasma density, and shape controls. They are typically picked manually based on the success of previous experiments. The actions of are part of a feedback control loop that aims to stabilize the normalized plasma pressure 3N. E S,, arguably one of the most important quantities since it measures the efficiency of plasma confinement rela- tive to the magnetic field strength. The third set of actions
3

---

Multi-Timescal Dynamics Model Bayesian Optimization for Plasma Stabilization in Tokamaks
Initial Conditions, Context (
Selected Action BB
Contextual Bayesian Optimization
Tokamak Reactor
ao
argmaxa 1 (a. (o ao (a.Pa 1)
Prior GP mean Eym
1. g(++
Observed Datapoint
(n+1) (R+11 CyM
Dynamics
Model
Figure 1. DynaBO Pipeline to generate feedforward trajectory actions. Initial conditions and feedforward actuators first used by the RPNN to generate rollouts through which we compute the prior mean of the objective function (time to tearing instability). Our Bayesian optimization algorithm uses this to optimize for actions (ECH). Noisy outputs from the Tokamak then used to update the Gaussian process model used for Bayesian optimization.
3.1. Recurrent Probabilistic Neural Networks and Instability Prediction
predictions often exhibit significant errors, largely due to the sim2real gap caused by lime-dependen fluctuations in the environment variables, e.g., due to maintenance or hardware changes provoked by previous experiments. Furthermore, retraining the RPNN between experiments and using it to select ECH profiles ag is virtually impossible because the newly collected data is too small and we only have few minutes between experiments.
We employ Recurrent Probabilistic Neural Network (RPNN) (Char et al., 2023b) to model the high-frequency behavior of the tokamak. An RPNN has a Gated Rectifier Unit (GRU) cell, which stores information about past states and actions. We use this network as state transition model which takes current state and action and outputs probabil- ity distribution the RPNN outputs a multivariate normal distribution with mean v) and variance D which we employ to approximate the system dynamics:
the next state. Given 8, and a, as inputs,
We address the above issues by employing Gaussian process (GP) model, a nonparametric model that is very data-efficient, especially low-dimensional spaces (Deisenroth & Rasmussen, 2011). A GP corresponds to an infinite collection of random variables, of which any finite number is jointly normally distributed. To fully leverage the strengths of GP models, we need to carefully summarize the information collected between experiments before training the GP. This is done as follows.
(.),
(4)
To bypass the issue that the number gyrotrons differs for each rollout the training dataset, we assume that the resulting heating profiles atah can be controlled directly, disregard af both in training and testing- approximate the
First, we assume the achieved normalized plasma pressure Oy is independent of the ECH profile og. This is a rea- sonable assumption because 3N largely determined by neutral beams, which are controlled through the feedback variables of. We then approximate the feedforward and feedback control actions of ....... of and af, ......... of by as- suming that they uniquely specified by the target normal- ized plasma pressure, denoted by 3r This choice is partly justified because the feedforward and feedback control ac- tions are often primarily informed by a target normalized plasma pressure. Finally, we employ the GP to predict the time-io-terring mode tTM. which we use as a proxy for the probability of tearing mode occurring. The rationale be- hind this choice is twofold. First, scenario where tearing instability occurs late implies higher degree of stability than a scenario where it occurs earlier. Moreover, this al- lows us to use the GP in a regression setting, where GPs are strongest and best understood. The GP inputs are thus B and ag. whereas the output is fTM-
allowing
Moreover, to facilitate optimization, ECH profile with stationary Gaussian curve Og. yielding the approximate ECH profile ag x af for all I. When carrying out experiments, we then project a, onto af. which can done for an arbitrary number of gyrotrons, i.e., for an arbitrary dimension of a?
In addition the RPNN, we train binary classifier, which call the tearing mode predictor h to predict the proba- bility of a tearing mode occurring based on input state and actuator.
11,8,, az) Bernoulli(y(se, at)).
3.2. Gaussian Process Model
Using RPNN rollouts exclusively for experimental design is challenging for various reasons. Although the RPNN accu- rately captures some of the tokamak behavior, the resulting

---

Multi-Timescale Dynamics Model Bayesian Optimization for Plasma Stabilization in Tokamaks
The GP is fully specified by a prior mean function i and kernel & that specifies the similarity between training inputs. In this work, we employ squared-exponential kernel k, which is appropriate for approximating most continuous functions. The prior mean m corresponds then to the aver- age Lime-to-tearing mode &TM predicted by autoregressive rollouts of the RPNN and tearing mode predictor:
where a balances exploration and exploitation. In con- ventional BO methods, the next step consists of setting dg (x+1)) = = measuring the time-to-tearing mode (TM and updating the GP accordingly. However, in our setting there is the added challenge that the target plasma 3N and the desired ECH profile corresponding to ag are not repro- duced exactly. This is due to various reasons, including a changing number of available gyrotrons, hardware failures, actuator noise, and unmodeled disturbances. To alleviate this issue, we measure the ECH profile and BN obtained dur- (n+1) ing the experiment and treat them as the true inputs dg and 3(x+1) used to update the GP model. Formally, this equivalent standard contextual BO where the GP inputs 0g and context By in equation 8 are perturbed by unknown noise.
(TM(3N, ag)= = E arg mint Te 2 0.5, Ti h(s),a))
81+1 ~ (((((es E)), (s), BN =
7 Σ
1
NN.
The dataset used to train the GP has the form
D, = ( 04 (MMi=,,......
(8) (0)
4. Results
We use it to compute the posterior distribution of (TM for arbitrary test inputs Bi: ag. which corresponds to normal distribution with mean and covariance
This section presents results from offline tests using his- torical data and results from experiments the General Atomics DIII-D Tokamak Fusion Facility. We use fixed RPNN in all experiments, trained using 15, 000 one-step state transition observations collected between 2010 and 2019 at the DIII-D tokamak. More details on the dataset can be found in A.2.
Pn (P (3) = (TM (5) + KT_ (K + -0 (6)
(BN, (g)) k.. - (K + ond ki + +o
KIT
respectively. Here on is the noise variance, [k =
k(3 ag. NN agg [K]
MG ai (i)
3W. aq (2)
Through our analysis of offline and online experiments, we aim to answer the following questions:
k.. = k( No Gg (..)). The vector [An]: TM
(i)
(TM( du contains the difference between the observed and the predicted Line-to-tearing mode. In practice, the posterior variance o² is typically small when evaluated in distribution and larger when out of distribution. Hence, in- tuitively, the posterior GP mean Hn, can be viewed as the predictive model, whereas of quantifies model uncertainty. This distinction is important for understanding Bayesian optimization, which we introduce in the next section.
1. How does DynaBO compare to other baselines? How robust is it in terms of the choice of kernel?
2. Can DynaBO find heating profiles that avoid tearing instabilities altogether using only handful of trials? If not, can it prolong the stable operation time of the plasma?
We address question 1 by conducting simulated experiments from offline data and comparing performance across all methods. We then address question 2 with results from live experiments on the DIII-D Tokamak. As we show in the following, both questions have an affirmative answer.
3.3. Contextual Bayesian Optimization with Noisy Inputs
Contextual Bayesian optimization is a data-efficient tool that leverages GPs to optimize black-box functions. Given a context that specifies the environment, it optimizes an acquisition function that carefully balances exploration ver- sus exploitation. By recursively updating the acquisition function after every observation, it gradually becomes more confident about its predictions, resulting in convergence. In every experiment, we treat the target normalized plasma pressure 3/1 g(n+1)), specified before the experiment, as the context and choose the BCH profile by optimizing the so- called upper confidence bound (UCB) acquisition function
4.1. Offline Data Analysis
This section employs historical data from the DIII-D toka- mak to compare DynaBO with several baselines. Specifi- cally, we employ data from 281 past experiments carried out at the DIII-D tokamak between 2012 and 2023. We selected the data points based their similarity with our live experiment, particularly the range of 3x and the high Jmin specification. Appendix A.2.2 provides a detailed description of the selection procedure.
BO- 42
arg max pm (ay, g(x++1)) + u0nnaa,,3 g(a+1)) 1), (8)
We employ the historical data to emulate our live experiment from Section 4.2. This is achieved as follows. At every
5

---

Multi-Timescale Dynamics Model Bayesian Optimization for Plasma Stabilization in Tokamaks
time step, we sample the target plasma pressure By from uniform distribution corresponding to the range of the historical data and condition DynaBO on 3y. We then select a subset from the historical data with plasma pressure values within the interval [3x c, By +J], where = 0.04, and use DynaBO to select ECH profile values corresponding to an element of that subset. After selecting an ECH profile, we treat the historical data point corresponding to that particular ECH profile new observation, which we use update our GP model.
kernel. DynaBO and DynaBO+time exhibit more variety in the queried ECH profiles than in the vanilla GP and the RPNN baseline. This corresponds to better exploration, resulting in lower regret for our approach. We observed similar trends using all other kernels except the linear one.
We additionally tested different acquisitions functions for our problem setup A.4. However, we finally chose Upper Confidence Bound (UCB) as our acquisition function be- cause it allows easy tuning of exploration exploitation during the actual experiment.
We compare our approach to four different baselines and ablations: the setting where we fully trust the RPNN to predict tearing modes without updating it, a vanilla GP with a zero-mean prior, a vanilla GP with a prior value which represents mean of the past data and our approach using deeeedeeendent kernel. The motivation for the latter ap- proach is that one could naively trust older data less than that closer to present-day experiments, attempting to model changes made due to deliberate hardware changes and re- pairs. In addition, we consider a linear kernel, Gaussian kernels, and Matérn kernels with different hyperparameter configurations to analyze our approach's robustness under varying model specifications.
4.2. DIIIID Tokamak Experiments
We tested algorithm at the DIII-D Tokamak hosted by General Atomics during two-hour time window allocated to us during the FY24 campaign. Each experiment run at DIII-D is known as a 'shot'. Each shot is then assigned a unique shot number.
To make the most of our time and make significant state- ments about results, we opted for pre-specified set of feedforward actuators oj that represents high performance experiments in high Jmin scenario.
In addition to our dynamics model which is trained on larger set of historical data, we further conditioned our GP on 125 historic high Qmin experiments. Appendix A.2.2 provides an overview of the data used train the GP. Our experiment consisted of 8 BO iterations with DynaBO.
In 2ig.2, we depict the cumulative regret
N
Cumulative Regret(NV)
Σ (7 max i=1
TM),
where mm = 10s is the maximal shot length. As can be seen, DynaBO and DynaBO with time dependency achieve the highest performance in all settings except the linear kernel setting. By contrast, the RPNN-based method and vanilla BO (with both zero and mean value priors) cannot consistently find good solutions despite performing more steps than the total number of data. This indicates that Dyn- aBO does not become overconfident and is robust to the choice of kernel and hyperparameters except when the ker- nel is clearly misspecified, c.g., when using a linear kernel. While including time as an input to the GP performs com- petitively, overall, the improvement seems only marginal. One possible explanation is that the reliability of the data depends on multiple factors, many of which cannot be ex- plained exclusively as a continuous function of time, e.g., sensor and actuator upgrades, particle absorption and re- lease by the tokamak wall, and the presence of impurities which may have been used in preceding experiments.
After each run of DynaBO, the selected heating profiles BO were converted to gyrotron angles and entered into the Plasma Control System, the interface that controls the tokamak. After a few seconds of maintaining the plasma, we ramped down the actuators and terminated the shot.
0g
We started our experiments by recreating a high-performing historical high gmin experiment with tearing instability, For this, we recreated the conditions in shot 180636, a plasma shot executed previously at DIII-D, Once this shot with positive tearing instability was recreated, we ran addi- tional shots where we varied the ECH using DynaBO while keeping the remaining settings identical. The details for each of the 8 shots carried out using DynaBO are shown in Table 1. As can be seen, our algorithm was able to avoid tear- ing instabilities in 4 out of 8 shots successfully. Moreover, we maintained a stable plasma for 3339 milliseconds on average. Although this number of shots is too low to be sta- tistically significant, we stress that the chosen configuration is very challenging. For reference, there were 61 historical experiments at DIII-D with very similar settings (high Jmin experiments that maintain high normalized plasma pressure NN > 3.0). Of those experiments, 47 reported tearing insta- bilities, corresponding to tearing instability rate of 77%. This set of experiments is selected as a baseline considering that these experiments also use ECH with current drive (sim-
We note that the vanilla GP does converge after more than 500 steps, i.e., after the plots in Fig. 2 end. However, such a long convergence time is unacceptable for our setting since fusion experiments are very costly, and we only get handful of experiments to explore.
To compare exploration. in Fig. 3, we display the ECH
profiles queried by the different baselines using a Gaussian
6

---

Multi-Timescale Dynamics Model Bayesian Optimization for Plasma Stabilization in Tokamaks
10 Cumulative Regret (Log Scale) Kernel RBF, length scale 0.01
Cumulative Regret (Log Scale) Kamel RBF, length scale 0.1
1d
Curulative Regret
Cumulative Regret
10
10
Daseline (RPNN only)
Baseline (RPNN only)
Zero Prior + GP
Zero Prior + GP
Mean value Prior GP
Mean value Prior + GP
RPNN Prior GP (Ours)
RPNN Prior + GP (Ours)
RPNN Prior (with Time Context]
RPNN Prior + GP (with Time Context)
10a
0
100
200
300
400
500
100
200
300
400
500
Iterations
Iterations
Cumulative Regret (Log Scale) Kernel RBF, length scale
Cumulative Regret (Log Scale) Kernel matern, nu 0.5
10
105
Ctvv Regret
Cumulative Regret
10a
10*
Baseline (RPNN only)
Baseline (RPNN only)
Prior
Prior
RPNN Prior (Ours)
RPNN Prior GP (Ours)
RPNN Prior GP (with Time Context)
RPNN Prior GP (with Time Context)
10a
10a
300
500
100
200
400
Iterations
Iterations
Cumulative Regret (Log Scale) Kernel matern, nu 2.5
Cumulative Regret (Log Scale)
Kemel
Linear
10a
Curmllatt laive Regret
Cumulative Regret
10
10+
Baseline (RPNN only)
Baseline (RPNN only)
Zero Prior GP
Zero Prior + GP
Mean value Prior + GP
Mean value Prior + GP
RPNN Prior (Ours)
RPNN Prior + GP (Ours)
RPNN Prior (with Time Context)
RPNN Prior + GP (with Time Context)
10a
0
100
200
300
400
500
100
200
300
400
500
Iterations
Iterations
Figure 2. Cumulative Regret (log scale) achieved by DynaBO (green), DynaBO with time-dependent kernel (gray), the RPNN only (bluc), vanilla GP with a zero-mean prior (red) and vanilla GP with mean value as prior (orange) using six different kernels.
ilar to our experiment), however, they use manual methods for selecting ECH profiles. The average time to tearing insta- bility in these experiments is 2424 milliseconds, well below our average of 3339 milliseconds. More details on identifi- cation of tearing mode instabilities in experiments from raw signals is shown in Section A.5.
shortcoming is that our current method is only applicable to feedforward control scenarios. This means the algorithm cannot adapt to unexpected real-time changes the plasma, such as MHD activity or impurity changes. In future work, aim extend our learning feedback control systems.
6. Conclusion
5. Limitations and Future Work
In this work, motivated by the challenges of tokamak control, we develop a multi-scale modeling approach for decision making on complex real world systems with limited avail- ability of data. Our pipeline leverages high highffrequency neural network model of the system dynamics and Gaus- sian process that makes predictions based on marginal statis- tics. Together, both models form Bayesian optimization
While our method is shown to preemptively suppress tearing instabilities, mainly data-driven, and potential improve- ments are possible by incorporating physics knowledge. One possible solution is to develop physics-infonen informed neu- ral network models, such incorporating elements of the Rutherford equation to improve interpretability. Another

---

Multi-Timescale Dynamics Model Bayesian Optimization for Plasma Stabilization in Tokamaks
ECH Profiles Queried by Different Methods
RPNN Prior+GP (Ours)
Mean Value Prior + GP
Prior+GP (Time Context)
Zero Prior GP
Baseline
3.0
3.0
3.0
3.0
10
2.5
2.5
2.5
2.5
2.5
ECHμ
ECHP
ECH
2.5
1.5
1.3
1.3
1.0
1.0
1.0
2.0
1.0
0.5
0.5
0
50
150
200
50
100
150
200
1.00
150
200
100
150
200
100
150
200
ECHu
ECHμ
ECHμ
ECHμ
ECMμ
Figure 3. ECH Profiles queried by different methods during simulated oftine using Gaussian kernel. We see that DynaBO and DynaBO with a time-dependent GP explore the most, highlighting the importance of our dynamic model prior mean.
Experiment ID
Target 3N
Tearing Instability
Stability Time
(Shotnumber)
Avoided
(ms)
199599
3.37
Yes
4566
199601
3.27
Yes
4632
199602
3.27
No
2107
199603
3.27
No
2149
199604
3.27
Yes
4592
199605
3.14
No
1512
199606
3.45
No
3512
199607
3.43
Yes
3654
Table 1. Results from two-hour experiments at DIII-D Tokamak DynaBO avoids tearing modes in 4/8 runs in a high-performance configuration with a historical rate of occurrence of the tearing instabilities of 77%. The time under stability with DynaBO is 3339 ms while the historical time under stability is 2424 ms corresponding to 914 ms improvement in stability duration. Generally, for stable experiments DIII-D, plasma stability is maintained for 4-5s.
algorithm tailored to the task at hand that can quickly iden- tify stabilizing control actions. This achieved by making decisions on the fly based on newly collected data. On historical data set, method outperforms vanilla BO and naive baseline that relies exclusively neural network predictions. This is mainly due to our approach having better exploration capability. Moreover, our method shows promise in live experiments on the DIII-D Fusion reactor. During the experiments, our approach successfully avoided tearing instability in 4/8 runs despite highly unstable condi- tions, representing an improvement of over 117% compared to past experiments.
and DE-SC0015480.
Disclaimer
This report prepared account of work sponsored by agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or im- plied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, represents that its would not infringe privately owned rights. Reference herein to any specific commercial product, process, or vice by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, rec- ommendation. or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.
Acknowledgments
We would like to thank the exceptional staff at the DIII-D National Fusion Facility that helped with the preparation and running of this experiment. This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Fusion Energy Sciences, using the DIIIID National Fusion Facility, a DOE Office of Science user facil- ity, under Awards DE-FC02-044R54698, DE-SC0024544
8

---

Multi-Timescale Dynamics Model Bayesian Optimization for Plasma Stabilization in Tokamaks
Impact Statement
Chua, K., Calandra, R., McAllister, R., and Levine, S. Deep reinforcement learning in handful of trials using proba- bilistic dynamics models. Advances in neural information processing systems, 31, 2018.
Our work illustrates the potential of combining complex high-frequency and low-frequency models address chal- lenges of complex real world systems like Tokamaks. Solv- ing nuclear fusion is one of the most important goals of our time and our work aims to utilize machine learning research methods to advance this field. The need for such methods in this area will increase in the future, as new and larger reactors such as ITER become operational, and a significant gap between existing and new models needs to be bridged with very little data. This is the case not only for the stabi- lization setting considered in this paper but also for settings such ramp-up design, where a different set of actuators is considered. Moreover, we believe this approach could be of interest to several other applications where the discrepancy between past and present data is considerable, e.g., complex physical systems with changing dynamics.
Degrave, J., Felici, F., Buchli, J., Neunert, M., Tracey, B., Carpanese, F., Ewalds, T., Hafner, R., Abdolmaleki, A., de Las Casas, D., et al. Magnetic control of tokamak plasmas through deep reinforcement learning. Nature, 602(7897):414 419, 2022.
Deisenroth, M. and Rasmussen, C. E. Pilco: A model-based and data-efficien approach to policy search. In Proceed- ings of the 28th International Conference on machine learning (ICML-11), pp. 465-472, 2011.
Fu, Y., Eldon, D., Erickson, K., Kleijwegt, K., Lupin- Jimenez, L., Boyer, M. D., Eidietis, N., Barbour, N., Izacard, O., and Kolemen, E. Machine learning control for disruption and tearing mode avoidance. Physics of Plasmas, 27(2), 2020.
References
Abbate, J., Conlin, R., and Kolemen, E. Data-driven pro- file prediction for diii-d. Nuclear Fusion, (44:0077,, 2021.
Fujimoto, S., Meger, D., and Precup, D. Off- Off-policy deep reinforcement learning without exploration. In Interna- tional conference on machine learning, pp. 2052-2062. PMLR, 2019.
Bardóczi, L., La Haye, R., Strait, E.,, Logan, N., Smith, S.,
Richner, N., and Callen, J. Direct preemptive stabilization of neoclassical tearing modes by electron cyclotron cur- rent drive in the diii-d low-torque iter baseline scenario. Nuclear Fusion, 63(9),096021, 2023.
Gantenbein, G., Zohm, H., Giruzzi, G., Günter, S., Leuterer, F., Maraschek, M., Meskat, J., Yu, Q., Team, A. U., al. Complete suppression of neoclassical tearing modes with current drive at the ceattron-cyclotron--ooeo frequency in asdex upgrade tokamak. Physical Review Letters, 8566)1244, 2000.
Boyer, M., Wai, J., Clement, M., Kolemen, E., Char, I., Chung, Y., Neiswanger, W., and Schneider, J. Machine learning for tokamak scenario optimization: combining accelerating physics models and empirical models. In APS Division of Plasma Physics Meeting Abstracts, vol- ume 2021, pp. PP11-164, 2021.
Garnett, R. Bayesian optimization. Cambridge University Press, 2023.
He, T., Zhang, C., Xiao, W., He, G., Liu, C., and Shi, G. Agile but safe: Learning collision-free high-speed legged locomotion. arXiv preprint ar arXiv:2401.1758, 2024.
Char, I., Chung, Y., Neiswanger, W., Kandasamy, K., Nel- son, A. O., Boyer, M., Kolemen, E., and Schneider, Offline contextual bayesian optimization. In Wallach, H., Larochelle, H., Beygelzimer, A., d'Alché-Buce Buc, F., Fox, E., and Garnett, R. (eds.), Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
Holcomb, C., Ferron, J., Luce, T., Petrie, T., Park, J., Turco, F., Zeeland, M. V., Okabayashi, M., Lasnier, C., Hanson, J., Politzter, P., In, Y., Hyatt, A., Haye, R. L.,, and Lanc- tot, M. Steady state scenario development with elevated minimum safety factor on dii-d. Nuclear Fusion, 54(9): 093009, aug 2014. doi: 10. 1080000009-511554/9909999999
Char, I., Abbate, J., Bardóezi, L., Boyer, M., Chung, Y., Con- lin, R., Erickson, K., Mehta, V., Richner, N., Kolemen, E., et al. Offline model-based reinforcement learning for tokamak control. In Learning for Dynamics and Control Conference, pp. 1357-1372. PMLR, 2023a.
Hu, Q. M., Logan, N. C., Yu, Q. and Bortolon, A. Ef- fects of edge-localized electron cyclotron current drive on edge-localized mode suppression by resonant magnetic perturbations in diii-d. Nuclear Fusion, 64(4):046027, 2024.
Char, I., Chung. Y., Abbate, J., Kolemen, E., and Schneider, J. Full shot predictions for the diii-d tokamak via deep recurrent networks. In APS Division of Plasma Physics Meeting Abstracts, volume 2023, pp. UP11 1P9169 2023b.
Kolemen, E., Welander, A., La Haye, R., Eidictis, N., Humphreys, D.., Lohr, J., Noraky, V., Penaflor, B., Prater, R., and Turco, F. State-of-ther the-art neoclassical tearing mode control in diii-d using real-time steerable electron

---

Multi-Timescale Dynamics Model Bayesian Optimization for Plasma Stabilization in Tokamaks
cyclotron current drive launchers. Nuclear Fusion, 54(7): 073020, 2014.
deep reinforcement learning. Nature, 626(8000):746-731 2024.
Kumar, A., Fu, Z., Pathak, D., and Malik, Rma: Rapid motor adaptation for legged robots. Robotics: Science and Systems XVII, 2021.
Shalloo, R., Dann, S., Gruse, J.-N., Underwood, C., An- toine, A., Arran, C., Backhouse, M., Baird, C., Balcazar, M., Bourgeois, N., et al. Automation and control of laser wakefield accelerators using bayesian optimization. Nature communications, 11(1):6355, 2020.
Lee, J., Hwangbo, J., Wellhausen, L., Koltun, V., and Hutter, M. Learning quadrupedal locomotion over challenging terrain. Science robotics, 5(47))eabc5986, 2020.
Shields, B. J., Stevens, J., Li, J., Parasram, M., Damani, F., Alvarado, M., Janey, M., Adams, R. P. and Doyle, A. G. Bayesian reaction optimization tool for chemical synthesis. Nature, 590(7844):89-96, 2021.
Lehnen, M., Aleynikova, K., Aleynikov, P., Camphell, D., Drewelow, P., Eidictis, N., Gasparyan, Y., Granetz, R., Gribov, Y., Hartmann, N., et al. Disruptions in iter and strategies for their control and mitigation. Journal of Nuclear materials, 463:39-44 2015.
Shimada, M., Campbell, D., Mukhovatov, V., Fujiwara, M.,
Kirneva, N., Lackner, K., Nagami, M., Pustovitov, V.,
Uckan, N., Wesley, J., al. Overview and summary. Nuclear Fusion, 47(6):S1, 2007.
Levine, S., Kumar, A., Tucker, G., and Fu, J. Offline rein- forcement learning: Tutorial, review. and Perspectives on Open Problems, 5, 2020.
Wakatsuki, T., Yoshida, M., Narita, E., Suzuki, T., and Hayashi, N. Simultaneous control of safety factor pro- file and normalized beta for jt-60sa using reinforcement learning. Nuclear Fusion, 63(7):076017, 2023.
Mehta, V., Barr, J., Abbate, J., Boyer, M. D., Char, I,, Neiswanger, W., Kolemen, E., and Schneider, J. Auto- mated experimental design of safe rampdowns via proba- bilistic machine learning. Nuclear Fusion, 64(4):046014 2024.
Meneghini, O., Smith, S., Lao, L., Izacard, O., Ren, Q., Park, J., Candy, J., Wang, Z., Luna, C., Izzo, V., Grierson, B., Snyder, P., Holland, C., Penna, J., Lu, G., Raum, P., McCubbin, A., Orlov, D., Belli, E., Ferraro, N., Prater, R., Osborne, T., Turnbull, A., and Staebler, Integrated modeling applications for tokamak experiments with om- fit. Nuclear Fusion, 5588):083008 2015.
Mohamed, M., Zakuan, D., Tengku Hassan, T. N. A., Lock, S. M., and Mohd Shariff, A. Global develop- ment and readiness of nuclear fusion technology the alternative source for clean energy supply. Sustainability, 16(10):4089, 2024.
Ono, M., Berkery, J. W., Bertelli, N., Shiraiwa, S., Delgado- Aparicio, L., Menard, J. E., Sánchez- Villar, Á., Shah, K,, Shevchenko, V. F., Idei, H., al. Efficient eced non-inductive plasma current start-up, ramp-up, and sus- tainment for an st fusion reactor. Nuclear Fusion, 2024.
Osinski, B., Finn, C.,, Erhan, D., Tucker, G., Michalewski, H., Czechowski, K., Kaiser, L. M., Babaeizadeh, M., Kozakowski, P., Milos, P., al. Model-based reinforce- ment learning for atari. ICLR, 1:2, 2020.
Schuller, F. Disruptions in tokamaks. Plasma Physics and Controlled Fusion, 37(1 A)::155, 1995.
Seo, J., Kim, S., Jalalvand, A., Conlin, R., Rothstein, A., Abbate, J., Erickson, K., Wai, J., Shousha, R., and Kole- men, E. Avoiding fusion plasma tearing instability with
10

---

Multi-Timescale Dynamics Model Bayesian Optimization for Plasma Stabilization in Tokamaks
A. Appendix
A.I. Additional Fusion Background
Effect of ECH on Plasma :
Electron cyclotron heating waves released by gyrotrons interact with the plasma by being absorbed by electrons whose gyrofrequency matches the frequency of the ECH wave. This absorption is highly localized, making electron eyclotron heating (ECH) a powerful tool for precise plasma control. ECH increases the electron temperature at the absorption point and often reduces electron density, phenomenon known density pump-out (Wang et al., 2017). As density decreases, plasma rotation tends to speed up due to reduced inertia. By adjusting the toroidal injection angle of the EC wave, different effects can be achieved. When injected perpendicular to the toroidal direction, EC waves primarily heat the plasma. When injected parallel, they drive electron acceleration in the direction of the plasma current, a process called electron cyclotron current drive (ECCD). While the injection angle can be switched between shots, it cannot be actively adjusted during a single shot. In this experiment, the toriodal angle is set, such that we get both ECH and ECCD profiles which are of similar shape.
Details on high qmin plasma scenario The high qrai scenario refers to group of scenarios related with elevated values of ymin the minimum value of the safety factor profile. Under this umbrella of scenarios, there are three main groups: Qmin = 1.4. 9min = 1.5 2, and 9min > 2. The lowest of the range with 9min = 1.4 has shown promise as being stable to TMs, but did not have the greatest confinement, while the highest of Jmin > 2 was stable to TMs but had lower energy confinement. The middle range of qman = 1.5 2 has very desirable confinement but is very susceptible to TMs (Holcomb et al., 2014). The purpose of this experiment to work in that middle range of (min; referred to as the elevated gmin scenario, to stabilize TMs and achieve higher confinement than either of the other similar scenario options. In the elevated Imin scenario, 2/1 TMs are the most prevalent mode because they require the least energy to perform magnetic reconnection and form magnetic island. Other lower order modes like 3/1 or 5/2 can sometimes occur but are significantly less frequent as they require more energy to form a magnetic island. Stabilizing TMSs in the elevated Jmin scenario would show path forward for this high-confinement scenario as a possible operating scenario for a fusion power plant.
A.2. Dataset
Plasma trajectories on Tokamak consists of three phases. The ramp-up phase, where the gases are heated and pressure is increased to generate the plasma state where fusion occurs. During this phase, the normalized plasma pressure A. rises. Then, we enter the flat-top phase, where the plasma pressure By is sustained, allowing fusion to occur. In this phase, BN is mostly constant and the aim to maintain this state without instabilities. Finally, the actuators are gradually ramped down and the plasma is safely terminated as the shot concludes. In this paper, we stay in the flat top phase and aim to stabilize it. To create our dataset, we hence use only flat top data and only control actuators during this phase of the experiment.
Our complete dataset consists of ~ 15000 plasma trajectories from historical experiments at DIII-D Tokamak. The data contains signals from different diagnostics have different dimensions and spatial resolutions, and the availability and target positions of each channel vary depending on the discharge condition. Therefore, the measured signals are preprocessed into structured data of the same dimension and spatial resolution using the profile reconstruction and equilibrium fitting (EFIT). These shots contain many different signals, some of which are described below. The dataset consists of scalar signals defined every timestep and profile signals which are defined along 33 or 200 points along the radius of the plasma cross-section. These consist of temperature, ion temperature, pressure, rotation, safety (Q) factor and density. For these signals we first convert them into PCA components and select the top components which are able to explain 99% of the variance in data. The Electron Cyclotron Heating (ECH) profile we choose to control, is also defined at 200 points along the plasma radius. PCA is unable to describe ECH profiles, however they can be described well by a Gaussian curve and are hence parameterized by the center, width and amplitude of the curve. These 3 parameters form our parameterization og of the ECH profiles. The model state space 8t is shown in table 2 while the actuator space o, is shown in table 3.
A.2.1. DATASET FOR DYNAMICS MODEL
For training the RPNN we utilize this data set with data points every 20 ms in time intervals with trajectories having an average length of 5 seconds. The RPNN is trained to predict 4 given (s). We add tearing mode labels to this dataset and train random forest classifier to predict the probability of tearing modes at every time step. We tried incorporating tearing mode predictions inside the RPNN; however, this did not yield good results.. This is likely due to the formation of spurious correlations and causality issues formed by introducing tearing modes into the dataset.

---

Multi-Timescale Dynamics Model Bayesian Optimization for Plasma Stabilization in Tokamaks
State Variables
Dimensions
Normalized Plasma Pressure 3N
Scalar
Line averaged density
Scalar
Loop voltage
Scalar
Confinement Energy
Scalar
Temperature Profile
Decomposed to 4 PCA components
Ion Temperature Profile
Decomposed to 4 PCA components
Density Profile
Decomposed to PCA components Decomposed PCA Components Decomposed PCA components
Rotation Profile
Pressure Profile
Profile (safety factor)
Decomposed to PCA components
Table 2. Plasma Features used as state space for RPNN model.
Actuator Variables
Dimensions
Power Injected
Scalar
Torque Injected
Scalar
Target Current
Scalar
Target Density
Scalar
Magnetic Field
Scalar
Gas Puffing
Scalar
Shape Controls
6 Scalars
ECH Profile
Decomposed to Gaussian curve with mean, stddev, amplitude (/4, a, w)
Table 3. Plasma Features as actuator space the RPNN model.
A.2.2. DATASET FOR GAUSSIAN PROCESS
To create the dataset for offline testing DH we first limit ourselves to High Jmin trajectories which achieve high normalized plasma pressure Ny > 3.0. This constraint follows our experiment conditions. This leaves us with 281 trajectories. We subsequently convert this data from a time-step scale trajectory level scale. We take average fyy of the flat-top phase of the trajectory. For ECH profile dy: take mean of all profiles in the flat-top phase of the experiment, This is the phase where the high-energy plasma state maintained. We thus get the dataset D where DMd consists of triplet (8x ay. fTM) ie. the normalized plasma pressure, parameterized ECH profile and the observed time-to-tearing mode. This dataset is used for offline testing.
For online testing, we subset this dataset further by only keeping whose ECH profiles are lie in the achievable parameter space as per experiment configuration, which leaves us with 125 training points. This used as a training set for the Gaussian Process before testing on the real tokamak.
A.3. Training Details for Recurrent Probabilistic Dynamics Model (RPNN)
Network Architecture:
Encoder:
Fully Connected (FC) layer: input_dim x 512
FC layer: 512 512
Memory Unit:
- Gated Recurrent Unit (GRU) block: 512 x 256
Decoder (with residual connections between FC layers):
FC layer: 256 x 512
FC layers: 512 512 (repeated times)
12

---

Miles-Timescale Dynamics Model Bayesian Optimization for Plasma Stabilization in Tokamaks
Tokamak Cioso section
Tokamak Top View
Gynotron
270"
230"
Power absorbed [Milim')
0.60
0.30
Driven current [MAim")
0.100
0.030
0.000
0.025
0.0
Figure 4. Gyrotron action on the Plasma inside the Tokamak. The bottom 2 curves indicate the power absorbed (heating profile) and current driven the plasma from the centre outer region of the plasma.
FC layer: 512 x 128
Output Heads:
Mean head: 128 x output dim
Log-varianee head: 128 x output.dim
The network predicts the parameters of a probability distribution, and is trained using a log-likelihood loss. We use the Adam optimizer with a learning rate x 10-4 and weight decay of 1 10-3 Early stopping applied with patience of 250 epochs based on performance validation set comprising 10% of the total data.
A.4. Comparing Different Acquisition Functions for Offline Experiments
We tested different acquisitions functions for this problem setup. We tested Expected Improvement, Thompson Sampling, and Upper Confidence Bound (UCB). There results of running the offline experiments shown below. All acquisitions functions gave similar order cumulative regret at the end of 500 epochs. We finally chose UCB our acquisition function because it allows easy tuning of exploration vs exploitation during the actual experiment.
Acq. Fun.
RBF (ls 0.1)
RBF (Is 1)
Matern (v 2.5)
Mean
STD
Mean
STD
Mean
STD
UCB
14418.2
8038.94
10726.4
3717.18
10223.0
3997.86
Thompson Sampling
11184.0
2835.95
15696.6
4091.07
18139.6
5818.07
EI
8201.0
10183.31
10342.6
9859.99
10231.8
9582.02
Table Cumulative regret values at end epochs under offline testing setup different acquisition functions across different kernels and length scales. The orders of cumulative regret are similar across all acquisitions functions.
A.5. Details about Online Experiment Results
In this section, we analyze the signals from our experiment runs at DIII-D. The results are shown Fig. and Fig. 6. We show the nlmms signal which measure the n=l magnetic perturbations. We also show the normalized plasma pressure 8x. a quantity which directly corresponds to energy confinement levels the plasma, For experiments 199606-199607, 199607, it tricky
13

---

Multi-Timescale Dynamics Model Bayesian Optimization for Plasma Stabilization in Tokamaks
detect a tearing instability, hence also add power injected. A sustained high value in rlrms along with drops in BN usually denotes tearing modes. However, (xx drops may vary depending severity.
n1rme 199599 (MHD)
n1mm 199601 (MHD)
beton
199599 (EFITO)
belan 1996C1 (CTIGC1)
690
2000
400
69C
(a)) Shots 199599 & 199601 No tearing modes were observed in these shots
ntrms (M (D)
nilrms 139603 (MHD)
Tanning Mode
Tearing
boton 180602 TO1)
boton 100603 (5PTOO)
Decrease plaama pressure
Decresse plasma
((b) Clear Tearing mode happens in 199602 which leads loss normalized plasma pressure Bx. In 199603 we see a tcaring mode form however its difficult to spot as the loss in By is minor.
n1rms 199604 (Mr)
n' 189605 (WHD)
Tearing Mode
HMA
belan 199304 (CF TO1)
beton 199605 (EFITO
Decrease in plasma pressune
(()) No tearing mode happens in 199604. A tearing mode happens in 199605
Figure 5. Identifying Tearing modes from Raw signal data. We use nlrms signal (denotes magnetic perturbations) and normalized plasma pr identify tearing modes. sustained high nlmms signal denotes learing modes. We label the drop in 8x due tearing mode formation. Note that in all experiments, By drops towards then end as power injected is dropped to safely end the experiment
14

---

Multi-Timescale Dynamics Model Bayesian Optimization for Plasma Stabilization in Tokamaks
n'rms 199808 (MHD)
Tearing Made
nimms 199607 (MHR)
batan 199606 (LH1001)
beton 199607 (ERTO))
pin.
199806 (NB)
pinj 199607 (NB)
100
1.00
50.1
MAR
2000
4000
635
Figure this figure we also include the power injected (pinj) along with nlrms and 8N. 199606, we see a very late tearing mode which occurs just before power injected is dropped. Very low loss in By is seen due to the tearing mode. Finally, in 199607 no tearing modes are seen.
A.6. Approximating the Prior
The historical data used to train the RPNN and the GP does contain the target normalized plasma pressure B Instead, it only contains the actions d, achieved during the shot. Similarly, the RPNN is trained exclusively on the actions, and not on Bx. hence direct mapping from fy does not take place in the RPNN. In the experiments, we address these issues as follows. In the historical data, we set By to be equal to the average normalized plasma pressure, i.e.,
B
$0 10
(10)
This is reasonable assumption since the target B is mostly achieved in practice. We then approximate the time-to-tcaring mode (BN, ag) predicted by the RPNN given By and ag as follows. We first use ag to compute the actions aeh We then compute the remaining actions of and of by sampling full rollouts from the historical data and setting of and a equal the corresponding actions. We then look the resulting average normalized plasma pressure and set it equal I(0) We do this for all ECH actions a, within 10 10 x 10 grid within the space of ECH parameters, specified by the historically largest and smallest parameter values in the historical data set. We then separate the results into bins that have the same value 3 up to margin of 0.04, and average over all tearing modes within that bin, yielding (ā a,). At test time, we project all points the closest point on the grid, both when performing queries and before updating the GP model.
A.7. Conversion of ECH Profile Gyrotron Angles
Even though we selected ECH profiles as our action space, the Plasma Control System (PCS) at DIID tokamak expects the output to be Gyrotron angles, which denote locations where they will be aimed. To make this conversion, used OMFIT software (Meneghini et al., 2015). We selected ECH profiles action space instead of gyrotron angles because experiment time one does not know how many gyrotrons available. With this choice of action space, we ensure our method is agnostic of number of gyrotrons.
15


Reviews:
[{"review_id": "oq6BMQIpUH", "reviewer": "Reviewer_zDdy", "summary": "- Machine learning algorithms face difficulties in controlling complex real-world systems, i.e., nuclear fusion.\n- Existing methods like reinforcement learning and Bayesian optimization fall short of fully addressing these challenges.\n- A new approach integrates a high-frequency data-driven model with a low-frequency Gaussian process to improve real-time adaptability.\n- This method, validated on the DIII-D nuclear fusion plant, shows significant improvement over baseline approaches.\n- Live experiments demonstrate a 50% success rate, marking a 117% improvement compared to historical results.", "questions": "", "limitations": "", "rating": 2, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you kindly for your review. We have addressed your questions and comments in detail below.\n\n**Definition of frequencies:** By frequencies, we mean the resolution of the input signal. The RPNN takes as input the full step-by-step actuator signals, which have a frequency of 50Hz, and makes predictions on a recursive, step-by-step basis, at the same frequency. By contrast, the Gaussian process model takes only the average ECH and betaN as input, which corresponds to a low-frequency smoothed version of the ground-truth actuator signals, and outputs the residual of time to tearing mode. (In other words, the RPNN is a high-resolution model, whereas the Gaussian process is a low-resolution model.) Based on this discussion, we have substituted the sentence you indicated with a clearer one.\n\n**Comparison with GP with marginal statistics:** As you suggested, we have added the relevant comparison with a GP with marginal statistics as a prior mean. The table below shows the cumulative regret value after 500 iterations under the setup mentioned in the offline experiments setup. We observed that a prior mean corresponding to the historical mean of the data yields similar performance to that of a zero mean prior, indicating that a constant prior is insufficient for the present setting.\n\n|| RBF ls.=0.01 | | RBF ls.=0.1 || RBF ls.=1 || Matern nu=0.5 | | Matern nu=2.5 |  | Linear ||\n| ----- | ----- | ----- | ---- | ---|----|----|-----|-----|-----| ------ | -------- | ------ | \n|  | Mean     | Std. | Mean     | Std.| Mean  | Std. | Mean  | Std.| Mean | Std. | Mean | Std.  | \n| Mean Prior  | 26899  | 30114.68 | 26971.6 | 30055.93 | 29228   | 28225.8  | 29478   | 28039.02 | 29417.8 | 28064    | 56309.4 | 6611.14  |\n| Zero Prior  | 27118 | 29939.2  | 27167.2 | 29895 | 30252 | 27432.5 | 29707.4 | 27845.35 | 29742.2 | 27798.6 | 51926.6 | 11134 |\n\n**Initial conditions and context:** The experiment configuration determines the initial conditions, which are agreed upon with DIII-D scientists before the experiment. The context for the GP is the target betaN, which is also decided before the experiment. We will make this clearer in the revised version. \n\n**Details of recurrent probabilistic neural network:** \nNetwork Architecture - \n\nEncoder : \n\n1. FC layer (input dim x512)\n2. FC layer (512x512)\n\nMemory Unit :\n1. GRU Block (512 x 256)\n\nDecoder : Residual connections used between FC layers\n1. FC Layer (256x512)\n2. FC layer (512x512) repeated 8 times\n3. FC layer (512x128)\n\nThis connects to two outputs\n1. Mean head (128 x output dim)\n2. Log Var head (128 x output dim)\n\nThe network predicts distribution parameters as output; hence, we train with a log likelihood loss. We use the Adam optimizer with a learning rate of 3e-4 and weight decay of 1e-3. We also use early stopping (patience = 250 epochs) on a validation set with 10% of the data points. We will include these details in the appendix of the revised manuscript.\n\n**Equation (5):** We train a binary classifier that outputs a probability of tearing mode at time t given the state s_t and action a_t. The output specifies a Bernoulli distribution (conditioned on s_t and a_t), which we sample to predict the probability of a tearing mode occurring. We will make this clearer in the revised manuscript.\n\n**Experiments with varying lengthscale:** The main goal of varying the lengthscale is to show that our method is robust to the choice kernel and hyperparameters, particularly compared to other approaches. We will state this more clearly in the revised manuscript.\n\n**Choice of acquisition function:** We ran more offline experiments with different acquisition functions:\n| Acq. fun. | RBF ls 0.1 || RBF ls 1 || Matern nu 2.5 ||\n| ----|---|-----|----|------|-----|-----| \n||Mean| STD | Mean | STD | Mean | STD |\n| UCB | 14418.2    | 8038.94  | 10726.4  | 3717.18 | 10223 | 3997.86 |\n| Thompson Sampling | 11184      | 2835.95  | 15696.6 | 4091.07 | 18139.6 | 5818.07 |\n| EI | 8201 | 10183.31 | 10342.6 | 9859.99 | 10231.8 | 9582.02 |\n\nThe numbers show the cumulative regret at the end of 500 iterations. We noticed similar performance across all acquisition functions. We chose UCB because it allows easy tuning of exploration vs exploitation during the actual experiment.\n\n**Lengthscale for Matern kernel:** As mentioned above, the offline experiments with varying lengthscales aim to analyze our approach's robustness under varying model specifications. With this in mind, we first observed the lengthscale obtained by log-likelihood maximization over the full data set, then chose lengthscales for the offline experiments that varied from the log-likelihood maximum by at most an order of magnitude. For the live experiments, we used the lengthscales that yielded the highest fit according to the log-likelihood loss.\n\n**GP data efficiency:** Thank you for pointing this out. We will revise the statement accordingly.", "claims_and_evidence": "- It solves an interesting application of Bayesian optimization to nuclear fusion.\n- Could you elaborate this sentence \"a recurrent probabilistic neural network models the high-frequency model dynamics, while a Gaussian process models the effect of low-frequency marginal statistics on the dynamics\"? In particular, what is the definition of frequencies here? What is the meaning of high and low frequencies?\n- For the setting of a GP prior mean, can you compare your method to a GP prior mean with simple basic statistics.  For example, you can set it as the arithmetic mean of historical data.", "methods_and_evaluation_criteria": "- How do the authors determine initials conditions and context?\n- Details of the recurrent probabilistic neural network are missing.  For example, neural network architecture, training loss, and training scheme are missing.\n- I don't understand Equation (5). How do you handle the Bernoulli distribution here?  Is it not a logistic function?\n- Why is GP-UCB selected for an acquisition function?", "theoretical_claims": "It is not applicable for this work.", "experimental_designs_or_analyses": "- In Figure 2, why did the authors vary a length scale?  Didn't the authors determine the length scale optimizing a GP model?\n- In Figure 2, how did the authors choose a length scale for the Matern kernel?", "relation_to_broader_scientific_literature": "It is related to broader scientific literature on Bayesian optimization.", "essential_references_not_discussed": "This paper seems to discuss essential references.", "other_strengths_and_weaknesses": "Please see above.", "other_comments_or_suggestions": "- In Page 4, the authors argue \"a Gaussian process (GP) model, a nonparametric model that is very data-efficient.\" I don't think GP is data-efficient.", "questions_for_authors": "Please see above.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "8ThI2LMgVR", "reviewer": "Reviewer_GHhS", "summary": "This paper introduces a multi-scale Bayesian optimization approach, termed DynaBO, specifically designed to control tearing instabilities in tokamaks. The approach integrates a high-frequency Recurrent Probabilistic Neural Network (RPNN) with a low-frequency Gaussian Process (GP), allowing rapid adaptation between experiments. Offline experiments on historical data from the DIII-D tokamak demonstrated significant improvement over baselines, while live experiments at DIII-D resulted in a 50% success rate in suppressing instabilities, showing a 117% improvement over historical outcomes.\n\nStrengths\n\nThe paper addresses a highly relevant, real-world application with significant practical implications for fusion energy.\nThe experiments conducted on the real DIII-D tokamak provide a compelling demonstration of the method's potential.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for reviewing our paper. If you feel that our responses to other reviewers improve the quality of the paper, we would be very thankful if you would consider changing your score accordingly.", "claims_and_evidence": "See above", "methods_and_evaluation_criteria": "See above", "theoretical_claims": "See above", "experimental_designs_or_analyses": "See above", "relation_to_broader_scientific_literature": "See above", "essential_references_not_discussed": "See above", "other_strengths_and_weaknesses": "see above", "other_comments_or_suggestions": "see above", "questions_for_authors": "see above", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "Ql9nVSga7h", "reviewer": "Reviewer_j9tX", "summary": "The authors propose a principled approach using Bayesian optimization to optimize stabilizing actions to efficiently stabilize a tokamak (a type of nuclear fusion reactor). The author’s proposed approach integrates both more-reliable observed data from the tokamak when specific actions are taken, with less-reliable information from a data-driven dynamics model, into the GP surrogate model, in order to improve optimization performance. Experimental results demonstrate that the proposed approach improves tokamak stability over other approaches. Notably, the authors demonstrate that their proposed approach achieves a 117% improvement over past methods in avoiding tearing instability.\n\n## update after rebuttal\nAs I stated in my comment below, the authors sufficiently answered all of my questions in their rebuttal, and I maintain that this work should be accepted for all of the reasons stated in my initial review.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you kindly for your detailed review. We have corrected all the typos you pointed out. Here are the answers to your questions:\n\n**Lessons:** During the experiments, we realized the importance of having a compressed representation of the actuator space. The high data efficiency of our approach would not have been possible if the search space had been much more granular. At he same time, having a high-fidelity model was particularly important, as we saw it clearly outperforms models with more naive and less sophisticated priors. During our experiment, we also observed that the suggestions for ECH provided by humans sometimes differed from our approach. We think investigating ways to incorporate these preferences on the fly would be a very exciting avenue of research.\n\n\n**Using self-attention:** We briefly experimented with an architecture that included self-attention while developing our model. Although we observed good predictive performance on the training data, the performance on the test data was poor, indicating overfitting. We speculate this is due to insufficient data and the overall tasks' difficulty (it has been shown that transformers sometimes only become good with a large amount of data (https://arxiv.org/abs/2106.04554)). For this reason, we stopped using transformers early during our research. Instead, we chose the RPNN architecture, which has been known to achieve state-of-the-art results on fusion data (https://arxiv.org/pdf/2404.12416). The paper cited looks into this question more deeply and analyzes GRU vs LSTM. Although we did not focus on this question in this work, we feel that designing appropriate attention-based architectures for nuclear fusion dynamics prediction is a very interesting research direction.", "claims_and_evidence": "Yes.", "methods_and_evaluation_criteria": "Yes.", "theoretical_claims": "There are no theoretical claims or proofs in this paper as far as I am aware.", "experimental_designs_or_analyses": "Yes, I checked all experimental design setups and experimental results provided in the paper. All are valid as far as I am aware.", "relation_to_broader_scientific_literature": "This paper focuses on one specific application: stabilization of Tokamaks. Scientists who work with Tokamaks may be able to apply the proposed method in practice to more effectively stabilize these reactors. This application is relevant to the broader scientific community because nuclear reactor stabilization is a very relevant and challenging real-world problem.", "essential_references_not_discussed": "There are no missing references as far as I am aware. However, as I noted below, I am not at all familiar with existing literature on stabilizing Tokamaks or nuclear reactors in general, so if there is missing relevant work in this area, I might not be aware of it.", "other_strengths_and_weaknesses": "Strengths: \n\n1: Application: This paper focuses on one specific application: stabilization of Tokamaks. This application is a strength of the paper because stabilizing these nuclear reactors is clearly an important and relevant real-world problem. \n \n2: Compelling experimental results: The experimental results demonstrate a compelling improvement in stabilizing Tokamaks compared to past approaches. In particular, the 117% improvement in successfully avoiding tearing instability provides very compelling evidence that the author’s proposed method is better than past approaches for stabilizing Tokamaks. \n\nWeaknesses: \n\nThe primary weakness of this paper is that there is a lack of novel machine learning methodology, as it combines ideas from well known methods in i.e. Bayesian optimization to create a principled approach for the specific problem of stabilizing Tokamaks. However, I think that the importance of the application and the compelling experimental results are strong enough that the novelty of the methods used are less important here. I therefore think that this is a very minor weakness and that the paper should likely be accepted.", "other_comments_or_suggestions": "Typos: \n\n1: In abstract: “Results on Live experiments…”. “Live” should not be capitalized here. \n\n2: Near line 119-120: “... are kept fixed throughout each experiment roullout…”. “Roullout” should be “rollout”.  \n\n3: Near lines 156-158: “In the following sections, we individual components our method -” \nShould instead be: “In the following sections, we discuss individual components of our method.” \n\n4: Line 150: “measure the n=1 magnetic pertubations”, “pertubations” should be “perturbations”\n\n5: Table 2 Row 7 “Decomponsed” should be “Decomposed”\n\n6: Line 568 “suspectible” should be “susceptible”", "questions_for_authors": "Question 1: \nIn your work on designing an approach specifically to improve tokamak optimization, what lessons have you learned that you think might be generally applicable to researchers trying to apply Bayesian optimization to other challenging problem domains? In particular, integrating high and low fidelity information in a single optimization run is a common problem in BO. In this case, the high-fidelity observations are the reliable ones from the Tokamak, and the low-fidelity/less-reliable information is that which we get from the dynamical model. What generally applicable lessons have you learned regarding this (if any)?\n\nQuestion 2: \nThe RPNN uses a GRU cell to store information about past states and actions. What motivated this particular choice? It’s possible that performance of this model could be improved by using a few attention layers instead of a GRU. Using attention typically outperforms other approaches including LSTM, GRU, CNN, etc. (https://arxiv.org/abs/1706.03762).", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: This paper uses Bayesian optimization applied to stabilization in Tokamak nuclear reactors. Given the new "Application-Driven Machine Learning" track at ICML--to which this paper was submitted--I think this paper is clearly above the bar for acceptance under that track. As far as I can tell, there aren't substantial (or any?) criticisms leveled at the paper by the one rejecting review, and the two other reviewers are in favor. The one rejecting review hasn't updated since the author feedback period. The authors seem to get some pretty cool results here on a very important (from a societal perspective) application domain, and I'm not seeing an obvious reason resubmission would be necessary on this track.
Decision: Accept (poster)