Particle accelerators are among the most complex engineering systems in the world. They are crucially important to the study of the elementary constituents of matter and the forces governing their interactions. Tuning and controlling accelerators is challenging and time consuming, but even marginal improvements can translate very efficiently into improved scientific yield for an experimental particle physics program, where integrated accelerator runtime imposes a substantial cost.
To date, the most common approach to accelerator systems control has largely consisted of hand-tuning by experts, and has been guided by physical models whenever possible. However, accelerator physics is complex and highly nonlinear. While we may model accelerator beams with impressive and improving precision [1] , it is often the case that real-world systems are too complex to admit practical analytical or Monte Carlo-based models, leaving optimization and control to the intuition of experienced experts. Even though this process has been successful to date, it is arduous and likely contains hidden inefficiencies. Recently, deep learning [2] , the training of neural networks consisting of many hidden layers, has proven itself useful for complex control problems [3] [4] [5] [6] [7] . Notably, processes within accelerator complexes occur * stjohn@fnal.gov between microsecond and millisecond timescales, much faster than human operators can react.
In this study, we present a real-time artificial intelligence (AI) control system for precisely regulating an important subsystem of the Fermilab Booster accelerator complex [8] . Our ultimate goal is to achieve a tenfold improvement in precision for the regulation system we describe in this paper. To realize this system, we use machine learning (ML) for two key features: to develop surrogate models [9] that reproduce the behaviors of the real-world system, and to train an online agent to take control actions in the system. The online agent is developed using reinforcement learning (RL) [10] , a framework in which an artificial agent learns by interacting with its environment. This online agent will ultimately control the actual accelerator system. The surrogate model provides an initial "safe" environment in which to train and evaluate control algorithms, and does so at higher rates and in a more controlled fashion than real-time training could allow. In this work, to demonstrate the complete methodology, we show how the control algorithms perform within the surrogate model.
Additionally, we propose a scheme for the first onboard ML control system at the Fermilab accelerator complex. Our approach utilizes field-programmable gate arrays (FPGAs) for real-time responses and includes full integration into the Booster controls system for highspeed data ingestion. On-board implementation is important because our latency requirements are more stringent than in most AI industrial controls applications 2 (e.g., Ref. [3] operated at a 5-minute update interval). In advance of commissioning the full system during online operations, we present important results from running a pre-trained, static RL model in FPGA test benches. Fully realizing this system will inevitably create an important and versatile set of tools that could find application in many areas in accelerator controls and monitoring, thereby paving the way for an ambitious control scheme that eventually encompasses the entire accelerator complex. While open questions remain as to whether data-intensive RL algorithms can be effectively used in this context without an independent, high-quality simulator, the results presented in this paper are encouraging.
In Section II we briefly highlight some previous publications at the intersection of machine learning and accelerator control. Section III describes the Fermilab Booster accelerator and surrounding complex in sufficient detail to understand the specific controls problem we address. In Section IV, we describe the accelerator data used to train a surrogate model, as well as the data processing used to support our ML workflow. Next, in Section V we describe the ML algorithms for the surrogate model and candidate agents. Then, in Section VI we briefly detail the process of deploying ML control algorithms to an FPGA. Section VII concludes with a discussion of our plans for full RL-based control.
The field of AI for accelerators is quite active in the development of new applications. AI algorithms are poised to play a salient role in accelerator control, tuning, diagnostics, and modeling. Here we provide a brief overview with a focus on recent work.
Predictive diagnostics are important precursors for control networks because they enable modeling of accelerator functions to a high degree of precision. The authors of Ref. [11] use ML-based diagnostics to predict the longitudinal phase space distribution in a particle accelerator. Ref. [12] offers an example of AI applications for fault classification in superconducting radio frequency cavities. Furthermore, in groundbreaking conceptual studies [13] [14] [15] [16] , the authors built predictive networks for accelerator modeling that could serve as predecessors to control networks, but did not implement these in hardware.
In Ref. [17] the authors use a neural network to demonstrate control of the longitudinal phase space of relativistic electron beams with very fine time resolution. However, model-based methods face serious challenges in regard to managing complexity. We seek to train, test, and deploy model-free RL control algorithms, but nevertheless must initiate that training using a surrogate model that provides a safe environment in which to evaluate algorithms before deployment. Surrogate models are an important ML-based tool for understanding accelerators and for training other ML models because it is impractical to train from scratch using real accelerator hardware due to the risk of downtime and the absence of precisely repeatable history. The authors of Ref. [18] use neural network (NN) based surrogate models for the Compact Linear Collider (CLIC) final-focus system. Additionally, ML algorithms for tuning have been studied at the Linac Coherent Light Source (LCLS) at SLAC National Accelerator Laboratory, and at the Swiss Free Electron Laser (FEL). Bayesian optimization of the Swiss FEL is explored in Ref. [19] . Likewise, in Ref. [20] , the authors use Bayesian optimization via Gaussian processes for fast tuning in response to the need to change beam configurations frequently. Although very large Gaussian process models tend to converge slowly on large datasets, this approach is notably fast and data efficient, which is extremely important at accelerator facilities where beam time for studies and tuning is in direct competition with beam for science. Furthermore, in Ref. [21] , the authors extend the Gaussian process-based methods to employ physics-based simulations in cases where archival data is not available. The challenge addressed in this paper does not involve an explicit physical model for use by such a simulation. Although Bayesian optimization is justifiably popular because of its effectiveness for a broad scope of problems, we instead pursue an RL approach as we believe it can better respond to rare events.
Ref. [22] summarizes a workshop that functions as a partial review of ML for accelerator controls and diagnostics and Ref. [23] provides an overview of beam dynamics studies using AI at CERN. In the most closely related previous work [24] , the authors approach the controls problem of maximization of the CERN Low Energy Ion Ring multi-turn injected intensity in a similar way to how we approach rapidly cycling magnet power supply stabilization here. Ref. [25] similarly studies RL to control the RF system of the KIT Karlsruhe Research Accelerator (KARA) storage ring and improve micro-bunching instability. As a proof of concept, the authors deploy a RL environment on an FPGA-ARM system for solving a classic cart-pole control problem [26] .
Relative to the body of prior work on AI for accelerators, we offer new developments in the integration of realistic RL algorithms into embedded systems (FPGAs). By embedding the algorithm in an FPGA, we may take advantage of the very stable, low-latency performance offered by that platform, which is critical for controls in a particle accelerator. Further, we indicate how this result may be extended through the use of a statistical ensembles of agent models, enabling increased learning rates, prediction stability, and robustness against potential mode collapse in individual agents. Given models of sufficiently small size, FPGA implementation of the models enables us to run an ensemble of models in parallel in order to optimize decision stability at no additional cost to overall latency. The Booster rapid-cycling synchrotron receives the 400 MeV (kinetic energy) beam from the Fermilab Linear Accelerator (Linac) via charge-exchange (H − to H + ). This beam is accelerated to 8 GeV by synchronously raising ("ramping") the Booster accelerator cavities' frequency and the magnetic field of the combined-function bending and focusing electromagnets known as gradient magnets, which are powered by the gradient magnet power supply (GMPS) [27, 28] . The beam is extracted at peak kinetic energy, after which the system is returned to the injection state. This complete cycle repeats, sinusoidally varying the GMPS magnet current between programmed current minimum and maximum at 15 Hz.
Meanwhile, other nearby high-current, high-power electrical loads are varying in time, causing unwanted fluctuations of the actual GMPS electrical current, and thus fluctuations of the magnetic field in the Booster gradient magnets. The role of the GMPS regulator is to calculate and apply small compensating offsets in the GMPS driving signal, improving the agreement of the resulting minimum and maximum currents with their set points. The present GMPS regulator system is a proportional-integral-derivative (PID) controller [29, 30] . The power supplies that provide the combined DC and AC components of the desired gradient magnet current consist of a pair of three-phase series-connected silicon controlled rectifier (SCR) [27] bridges fired at 720 Hz. An LC filter network at the output greatly reduces 720 Hz ripple, resulting in an output voltage that is proportional to the sinusoidal program provided by the GMPS regulator system. The series-connected electromagnet circuits and their cell capacitor banks are driven at resonance at 15 Hz and coupled via a distributed choke system for bypassing DC current. Figure 2 shows the power supply output voltage due to the cosine program.
For monitoring purposes, a special series-connected half-cell reference magnet located in the equipment gallery includes a pickup coil located between its poles to measure the time rate of change of the magnetic fielḋ B, which is an important input to the GMPS regulator. Powered with the other gradient magnets and housed in a low-radiation environment without charged particle beam passing through it, this reference magnet provides an accurate representation of the magnetic field under control throughout the accelerator.
Timing information derived from˙ B = 0 synchronizes the GMPS regulator system to the minimum and maximum values of the magnetic field and provides a transistor-transistor-logic (TTL) based 15 Hz master clock signal that drives the timing system for the GMPS regulator and indeed the rest of the accelerator complex. The high-frequency sampled measurements near the minimum and maximum values of the magnetic field are automatically fitted each cycle, and the finite impulse response (FIR) parameters are used as the primary feedback mechanism for the GMPS regulator system. Reducing the errors (the difference between the target and realized GMPS current especially at injection) of the GMPS system is of primary concern in the operational performance and efficiency of the Booster. The following Sections discuss the details of the present and proposed regulation systems.
The GMPS regulation system seeks to minimize the impact of disturbances due to environmental factors such as ambient temperature, nearby high-power pulsed RF systems, and ramping power supplies with inductive loads. Variations in the AC line frequency and amplitude are also significant sources of error, and are due in part to other particle accelerators in the complex changing currents in their own high-current electromagnets as part of their normal operations. Without regulation, the fitted minimum of the magnetic field may vary from the set point by as much as a few percent.
The existing regulator reduces GMPS regulation errors to roughly 0.1% of the set value by implementing a proportional-integral-derivative (PID) control scheme. Each cycle, the fitted minimum and maximum of the magnetic field reflect the combined influences of the set points, any compensation applied by the regulator for that cycle, and any new influence of other nearby electrical loads. Calculated estimates for the minimum and maximum values of the changing magnetic field of the previous 15 Hz cycle are used to adjust the power supply program to decrease the errors of the system. The calculation of the compensated minimum current B:VIMIN, used in the next cycle of the power supply program voltage, proceeds as follows at each time step t:
where B:IMINER(t) and B:VIMIN(t) are the corresponding control system readings at time t; I min fit is the measured (fitted) current minimum of the present cycle; I min set is the nominal set value of the current minimum; Γ is the integral gain; and α is the proportional gain. The running parameter β is effectively an integral of the error which rapidly downweights past measurements of the error. Typical values of the gain constants are Γ = 7.535 × 10 −5 and α = 8.5 × 10 −2 , adjustable by system experts. B:IMINER is proportional to the error in the system, and for our purposes they are synonymous. The environmental perturbations, discussed above, increase the distributed long-term steady-state errors of the GMPS system. A traditional PID regulation loop, given a sufficient amount of time and unchanging perturbations, will decrease the steady-state error of a system to zero. In reality, and within the timescale of the Booster beam cycle, the PID loop decreases the steady-state error to some nonzero error, typically of order 0.1% for the cycle minimum. See Figure 3 for a sample distribution of measured errors for the minimum value of the sinusoidally varying magnet current. Efforts to decrease the steady-state error further by adjusting the closedloop gains to be more responsive would come at the cost of reduced overall system stability. Therefore, a balance between the steady-state error and stability of the system has been struck.
We seek to improve upon this regulation performance in order to decrease injection losses and improve Booster efficiency. An RL model, whose inputs include all of the most important outside influences on the GMPS system, can infer more accurate cycle-by-cycle compensations on an FPGA, reliably producing better regulation against those influences. Creating such an algorithm begins with collecting time-series data from the accelerator complex in operation.
We collected a dataset for Booster GMPS regulation to provide cycle-by-cycle time series of readings and settings from the most relevant devices available in the Fermilab control system [31] . This data was drawn from a select subset of the roughly 200,000 entries that populate the device database of the accelerator control network [32] . Data was sampled at 15 Hz for 54 devices that pertain to the regulation of the GMPS during two separate periods of time: from June 3, 2019 to July 11, 2019 and from December 3, 2019 to April 13, 2020 when the accelerator was shut down in response to the COVID-19 pandemic. In this paper, we use a subset of this dataset's devices and a single day-March 10, 2020-for development and demonstration, as detailed in Section V B. For a more detailed overview, please see the corresponding datasheet available online [33] .
In this accelerator complex data-logging nomenclature, device parameters with the B: prefix are related to the Booster, whereas device parameters beginning with I: are related to the Main Injector. Additionally, "MDAT" denotes the accelerator (machine) data communication broadcast. This expert-chosen set of limited parameters was used to prototype simplified algorithms and train the Booster surrogate model in Sec. V as well as test lightweight RL agent models in hardware in Sec. VI. However, in order to more fully explore the utility of the full set of logged signals for future studies, a Granger causality study was performed using variables correlated or anticorrelated with B:IMINER, with absolute Pearson correlation coefficient |r| > 0.20. B:LINFRQ and I:IB were also considered on advice of system experts.
For a given pair of concurrent time series-one potentially causal, the other potentially responsive-Granger causality does not prove pure causation.
Instead, Granger causality suggests the response variable could be better predicted using both time series jointly rather than using the response variable's self-history alone. This test consists of creating and comparing two linear regression models: a self-model and a joint-model for both the response variable and the potentially causal variable, and calculating coefficients for each lag value (time difference) being tested, from one up to some predetermined maximum number of time steps [34] . If at least one coefficient is not zero in the joint model, then the other variable is said to be "Granger causal" with respect to the response variable, at the lag value being tested. We compared p-values to test statistical significance at each lag value up to 50 lags (approximately 3.33 s) as well as looked at the difference between the Bayesian information criterion (BIC) of the self and joint-models. As a result of these iterative calculations, we identified three additional variables-B:VIMAX, the compensated maximum GMPS current, B:VIPHAS, the GMPS ramp phase with respect to line voltage, and I:MXIB, the main injector dipole bend current-that will be considered in the next iteration of the surrogate model.
The initial selection of B:IMINER, B:LINFRQ, B:VIMIN, I:IB, and I:MDAT40 formed a subset of the top eight "causal" variables identified through the causality study. Table I briefly summarizes these parameters of interest, used in this iteration of the surrogate model. Although the devices were configured to write out reading and setting data at 15 Hz, actual timestamp intervals varied from this nominal frequency , and timestamps were not well synchronized across devices. Thus, for time alignment purposes, we made use of the recorded timestamps for Event0C, a broadcast accelerator control event which is synchronized to the fitted minimum of the periodically varying magnetic field in the gradient magnets described in Sec. III. This control event serves as a logical choice of reference time for GMPS-related parameters. Using Apache Spark-based [35] algorithms to process the data in parallel, we first calculate the maximum interval between successive timestamps for each device across all 176 days (necessarily excluding the five-month gap between our two data-taking periods). Then, we use the corresponding largest observed lag between recorded values within each device to place the upper limit on a look-back window from an Event0c timestamp in order to align the most recent device timestamp within that window to a corresponding event timestamp.
ML refers to the process by which we adjust the randomly initialized parameters of generic function approximators, termed "models," so as to minimize an appropriately chosen loss function, or conversely to maximize a reward function. ML model architectures specify arrangements of "nodes," usually co-evaluated in layers, where each node calculates the weighted sum of the inputs and a bias term, and outputs the value of a nonlinear "activation function." The outputs of one layer typically provide the inputs to the next layer; in the simple multilayer perceptron (MLP) architecture, all nodes in each layer send copies of their output values to all nodes in the next layer. Layers not at the input or output are termed "hidden" layers. There are useful variations on this simple layer stack architecture such as recurrent neural networks (RNN), wherein some outputs from a previous forward inference are taken as inputs. Our work 6 these actions and presenting new situations to the agent. 1 The environment also gives rise to rewards, special numerical values that the agent seeks to maximize over time through its choice of actions. More specifically, the agent and environment interact at each of a sequence of discrete time steps, t = 0, 1, 2, 3, . . .. 2 At each time step t, the agent receives some representation of the environment's state, S t 2 S, and on that basis selects an action, A t 2 A(s). 3 One time step later, in part as a consequence of its action, the agent receives a numerical reward , R t+1 2 R ⇢ R, and finds itself in a new state, S t+1 . 4 The MDP and agent together thereby give rise to a sequence or trajectory that begins like this:
In a finite MDP, the sets of states, actions, and rewards (S, A, and R) all have a finite number of elements. In this case, the random variables R t and S t have well defined discrete probability distributions dependent only on the preceding state and action. That is, for particular values of these random variables, s 0 2 S and r 2 R, there is a probability of those values occurring at time t, given particular values of the preceding state and action:
p(s 0 , r|s, a)
. For a given architecture with its activation functions, the weights and biases are the parameters being adjusted in the optimization or "training." The "learning rate" sets the proportionality of parameter adjustment to gradients of improved performance (based on lower loss or higher reward).
RL is the subfield of AI aimed at optimizing of control or planning of complex tasks based on feedback inputs from an environment, as explained below. The main components of RL are the environment and the agent, a ML model, as illustrated in Fig. 4 . RL trains a model over many time steps, and the resulting model may then be taken as fixed, to be deployed on new data; the model may be expected to perform similarly on new data as it did in training, if the dynamics of the new data were well represented in the training data. However RL extends naturally to continuous online learning, which would allow our regulator to adapt to changing environmental dynamics such as seasonality or new modes of accelerator complex operation, even though we would initially deploy a static or infrequently updated model out of prudence. We set out to use RL to train an optimal regulation policy dictating which action the GMPS regulator should take in any given state of the system.
The environment, usually formulated as a Markov decision process (MDP), is represented by a timeindependent, discrete system with which the RL agent interacts (e.g. the accelerator complex). For the regulation of the GMPS current minimum, the environment includes the time-varying outside influences for which the GMPS regulator makes compensating changes. At each time step t, the environment takes in the control action A t determined by the RL agent based on the current state S t , and provides the new system state S t+1 (e.g. settings and measured quantities) along with an associated reward R t+1 . Optimizing the agent's policy actions is defined to mean maximizing the long-term integrated reward. In this study, the reward is calculated from the error in the minimum value of the GMPS current, B:IMINER:
The larger the magnitude of B:IMINER, the lower the re-ward. The possible actions A t we consider correspond to adjusting the value of B:VIMIN, the lone control variable.
Recently, significant progress has been made in RL by combining it with advances in deep learning. Deep learning models are well suited to representing complex policies for high-dimensional problems such as regulation in a dynamically variable environment. The deep Q-network (DQN) [37, 38] approach, which we adopt for this study, involves using a deep neural network to learn the actionvalue function, or Q-value, and is usually deployed in environments that take discrete control actions. The optimal policy can then be derived by choosing the action that maximizes the expected Q-value.
More formally, a policy π is used by an agent to decide what actions A t = π(S t ) to take given a state S t at time t. An optimal policy π * maximizes the Q-value,
where E is the expectation value operator, R t ≡ R(S t , A t ) is the reward at time t, and γ is the discount factor that de-emphasizes future rewards relative to present ones. For this study, a value of γ = 0.85 was found to be performant. The Q-value is the sum of the expected discounted rewards from the current time t up to the horizon T . In practice, the optimal actionvalue function Q * is not known a priori, but it can be approximated iteratively because it satisfies the Bellman equation [39] ,
In a DQN, the Q-value is approximated using a deep neural network, or policy model, with parameters θ. In particular, the loss function at a time t is given by the mean squared error (MSE) in the Bellman equation,
where the (unknown) optimal target values are replaced by the approximate target values y t = R t + γ max At+1 Q(S t+1 , A t+1 ; θ − t ) using parameters θ − t derived from previous iterations.
Continuous action space environments, such as the compensating adjustments of our GMPS current regulator, can use the DQN algorithm by discretizing the action space directly. However, this approach scales poorly because the number of samples required to estimate the optimal policy (compensating setting of B:VIMIN) at some fixed accuracy grows exponentially with the number of input variables. To avoid this, we discretize instead the change of control signal B:VIMIN using steps of just a few different sizes, including the option for zero-size change. By doing so, we explicitly limit the control signal variation between time steps, while helping to minimize the tuning and resonance of the overall surrogate-RL training loop. Other RL algorithms, such as deep deterministic policy gradient (DDPG) [40] , proximal policy optimization (PPO) [41] , and twin delayed DDPG (TD3) [42] can provide continuous control suggestions which fit the GMPS system more naturally. However, the sample efficiency and tuning stability of these algorithms presents difficulties as discussed in various studies [43, 44] .
The iterative nature of RL training requires an environment that accommodates trial and error. As mentioned in Sec. II, most real-world, complex systems cannot afford such an approach due to the risk of system failure, therefore some level of offline training is essential. To facilitate offline training of a control agent, here we develop a surrogate model of the Booster in the accelerator complex, using the select collection of historical measurements described in Sec. IV.
In the following subsections, we discuss the workflow employed to develop and use an offline policy optimization, a prudent step prior to deploying the model in the real system. Sec. V B details the surrogate model that is developed to test the deep RL agent. Then, in Sec. V C, we describe how the agent itself is trained.
In order to train the RL policy model, we adopted a long short-term memory (LSTM) [45] architecture for a surrogate of Fermilab Booster accelerator in its operating environment using data from the Main Injector and Booster systems. Its purpose is to predict the impact of changing the control variable B:VIMIN on B:IMINER over the course of the RL episodes.
An LSTM model is a specific type of recurrent neural network (RNN), which is appropriate for modeling sequential data such as our accelerator data. The key feature of an RNN is the network recursion, which enables it to describe the dynamic performance of systems over sequential time steps. One difficulty of training RNNs is the vanishing gradient problem [46, 47] , which is a result of the gradient at a recursion step k being dominated by the product of partial derivatives of hidden-state vectors h,
where R is the number of recursions and w is a recurring weight. If w < 1, then this product tends to zero for large R. The design of the LSTM includes an internal memory state, which effectively mitigates the vanishing gradient problem, and allows the network to retain a longer memory of past inputs [45] . Given the complexity of the GMPS regulation environment, we selected the LSTM architecture in order to capture the multiple frequency modalities observed in the data.
At each time step, the input to the Booster ML model is the most recent 150 time steps from the five in- While the latter two values are not included in the environment state propagated to the RL agent, the additional output structure helps to regularize the training process and validate performance. The LSTM model uses this 150-step look-back window to "recall" the historical patterns of the time-series data, and thus achieve high accuracy in prediction. The LSTM surrogate model was developed and implemented using the Keras library [48] . We used the Adam optimizer [49] and a cost function of the mean squared error (MSE) of the predictions. The total number of data samples used for the analysis of the LSTM surrogate was 250,000 time steps from March 10, 2020, which we split into two non-overlapping data sets composed of 175,000 time steps for training and 75,000 time steps for testing. The training samples were further split using a K-fold cross-validation method: we defined five cross-validation folds that split the training and validation in a 80%/20% split. This technique was used in order to estimate how the model is expected to perform in general as well as to monitor over-fitting. The loss values from the validation sample were then used to determine if the learning rate should be reduced or if the model had stopped learning, as shown in Figure 5 . After more than 300 training epochs, the figure shows a bifurcation between the values of training loss and validation loss, suggesting some overfitting. Therefore, we used the model parameters prior to this bifurcation as the parameters of our surrogate model. On separate test data, the loss value for this surrogate model was determined to be 9 × 10 −4 which is consistent with the training data set prior to the bifurcation. Overlaid time series from the data and from LSTM predictions for a selected time window are shown in Fig. 6 . Based on the great similarity of these results, the surrogate model was deemed adequate to use for initial training of the RL policy model.
For this study, we formulated the problem as an episodic Markov decision process, where every episode contains 50 time steps. As in all Q-learning, the agent learns to maximize the reward within the time horizon of an episode. We developed our RL workflow based on a variant of DQN, the double DQN algorithm [37, 50, 51] , using Keras [48] to optimize B:VIMIN settings dynamically to minimize GMPS error B:IMINER. The double DQN explicitly decouples the target model, which is used to evaluate actions, from the policy model, which is used to select actions, although they take the same form.
We used the OpenAI gym package [52] to develop the environment that serves as a wrapper around the virtual accelerator complex model described above in Sec. V B to interact with the RL agent. The observation state space is defined by the aforementioned five variables in the virtual complex model section, shown to causally relate to B:IMINER. The action state space only contains one free parameter of control: adjustments to B:VIMIN. Th seven discrete control options relative to the previous B:VIMIN are 0 (no change), ±0.0001, ±0.005, and ±0.001. The choice of these values was based on the actual distribution of the changes in B:VIMIN observed in the data.
At the start of each episode, 150 time steps from the data are used to initialize the system state, as is required for the Booster surrogate model. The 150th step defines the observation state used by the agent. For each step thereafter, the agent provides a new action specified by the change in B:VIMIN, and the system state is updated.
The new system state is then used to predict the next B:IMINER. After the prediction, the system state is incremented to the next time step. The current state, reward, and status for each step is passed to the agent to be used for training the DQN policy model.
During training, event samples are placed into a buffer before calculating the loss. This memory buffer is sampled randomly in a process called experience replay [37] in order to remove instabilities found to arise from training on time-ordered samples. Once the memory buffer has sufficient experiences (32 experiences for this study) the active policy model begins training and continuously updating. We use the -greedy [53] method to control the agent's tradeoff between exploration (random choice of action) and exploitation (deterministic action dictated by the current policy), in which the optimal action according to the current policy is chosen with probability , while a random action is selected with probability 1 − . At the beginning of the training session we set = 1 with a decay factor of 0.9995, applied multiplicatively whenever an exploration action is selected, until a minimum value of = 0.0025 is reached. For this study, we use a multilayer perceptron (MLP) as the policy model architecture, and rectified linear unit (ReLU) activation functions [54] , as summarized in Table VI A.The active policy model is continuously updated during training by using randomly selected experiences from the memory buffer. At each training step the weights of the target model θ target are incrementally updated to reflect the weights of the active policy model θ policy ,
where we set τ = 0.5 [40] . RL algorithms learn from the reward provided by the environment, which in this study is given by Eqn. 4. Additionally, we impose a penalty term if B:VIMIN and B:IMINER go outside of the defined parameter space. The penalty term is proportional to the number of steps taken to ensure that the agent learns to perform the necessary steps in one episode. The result of the DQN MLP training, in terms of the total reward per episode versus the number of training episodes, is shown in Fig. 7 . In this proof-of-concept study, each episode is initialized with the same 150 time steps. The current controller reward, the black dashed line, was determined using the data by summing the values of R t = −|B:IMINER(t)| for the same 150 time steps used to initialize the RL loop, and thus given the fixed history does not vary as a function of the episode number. Additionally, the current controller reward value was compared to the reward associated with the actual controller action (B:VIMIN) within the RL loop, and was within 5% of this real value. The DQN MLP controller reward, the solid red line, shows steady improvements during training and converges to approximately a factor of 2 improvement over the current system.
Fast GMPS control electronics are required to collect information from the Booster environment, decide whether to apply a corrective action, and distribute the corresponding control signal, all with low latency. An FPGA is a natural choice to implement the corresponding circuit, accommodating latencies far below those achievable with a CPU or GPU while allowing reconfigurability impossible in a custom application-specific integrated circuit (ASIC) solution. The DQN MLP control model proposed in Section V C requires an efficient, but adjustable, implementation of NN algorithms, strongly suggesting an FPGA-based implementation. As a preliminary step, we take the offline-trained DQN MLP with weights fixed and deploy it in an FPGA.
The following subsections review the computational steps required for a single NN inference ( §VI A); describe the basic elements of an FPGA and how a deep NN calculation can be efficiently mapped to a corresponding circuit ( §VI B); present an implementation of the DQN MLP described above in Section V and the impact of various design choices ( §VI C); and lastly discuss possible extensions of the implementation to accommodate more complex algorithms which are of interest ( §VI D).
The structure of an MLP is a series of alternating linear and nonlinear transformations (layers), with the i th layer mapping a set of inputs x i (features) to a discrete list of outputs y i . In the present application, the features may include any measurements of the GMPS environment, such as digitized traces from the reference magnet system, line voltage frequency, and equipment gallery temperature. For the DQN MLP, the outputs y i are scores associated to a discrete set of possible actions, with the highest-scoring action being the one taken by the controller. An MLP layer f yielding m outputs may be written in terms of its action on a set of inputs
where w ij (the n × m weight matrix) and b j (the mdimensional bias vector) are configurable parameters of the linear translation and σ is an m-to-m nonlinear activation function. For each layer, the activation function is prescribed as a part of the model architecture while optimal values for the weights and biases are found through a training procedure. The DQN MLP utilizes the linear (identity) ReLU(x i ) = max(x i , 0) activation functions. The complete, k-layer NN is specified by an ordered composition of layers y = f (1) f (2) . . . f (k) (x). While the input and output dimensions are fixed by the set of features and actions, the dimensionality of intermediate layers is arbitrary. An FPGA consists of an array of logic gates that may be programmed to emulate any circuit (up to the physical resource constraints of the specific hardware device). This allows FPGA designs to profit from the many advantages of custom ASICs including massive parallelization and low power consumption while maintaining reconfigurability. However, a significant advantage of the FPGA architecture lies in the fact that it is not simply a homogeneous fabric of low-level gates (e.g. NAND gates). Rather, modern FPGAs are heterogeneous structures including more complex logical blocks, each specialized for a dedicated task, repeated many times. In this way, FPGA designs can simultaneously exploit both the flexibility of a programmable architecture and the performance of a dedicated printed circuit.
An efficient implementation of the NN model in firmware requires a design that exploits the FPGA's specialized computational units to perform each step of the NN calculation. Digital signal processor (DSP) slices are flexible circuits for addition, multiplication, wide-array bitwise logical operations, and more. DSPs may be further chained to accommodate more complex operations. In the Intel Arria 10 FPGA, to be deployed in the GMPS control system, DSP blocks may be configured to multiply and accumulate fixed-point numbers up to 27 bits, providing a solution for the linear component of Equation 9. The affine map from m to n dimensions requires mn scalar multiplications and sums that, in a fully parallelized design, may be accomplished with mn cascading DSPs. To evaluate an arbitrarily complex activation function in FPGAs, it is more efficient to store a precomputed table of values than to re-calculate the function many times per inference. This may be accomplished using block RAM (BRAM), embedded memory that is configurable for read/write access. BRAMs are available in segments of 20 kb in the Arria 10 to store, for example, a bank of 1024 function values at 20-bit precision. Registers are groups of flip flops used to record temporary numerical values or internal states, and to facilitate signal routing across the major computational blocks of the design. Finally, ALMs are lightweight, configurable modules of combinational logic elements, used throughout designs for basic operations such as simple arithmetic and logical operations.
The GMPS control model described in Section V must be converted to firmware in a manner that takes full advantage of the FPGA's architectural features described in Section VI B. This is accomplished through the translation of the Keras description of the NN function into high-level synthesis (HLS) code using the hls4ml [55] toolkit, whose functionality has been recently extended to Intel FPGAs [56] . The HLS design is converted to firmware using Intel Quartus [57] . The use of hls4ml brings the significant advantage of enabling a fast development cycle from model prototyping to implementation in firmware. Thus, the present work has focused not only on achieving an optimal design for the benchmark ML algorithm proposed in Section V C but also on more generic design-space exploration. Establishing scalable strategies such as FPGA implementation is critical for scaling up to more complex ML models that will inevitably become necessary as larger data sets allow for increasingly nuanced treatments of the control problem.
The conversion of the Keras model to firmware requires a number of design choices. Chief among these are the numerical precision to which the calculation is carried out and the degree of parallelism incorporated into the design. Fixed-point values with a specified number of total and integer bits are used to represent model inputs, weights, and all intermediate results of the calculation.
To determine the range of values to encode in the fixedpoint representation, the number of integer bits is set to be at least as large as the maximum weight value. The number of total bits, which sets the number of bits used to encode the fractional component of the weight (once the integer part is specified), is set to minimize the impact of quantization. Figure 8 displays a histogram of the weight values for the trained DQN MLP model. The total number of bits retained for each weight is selected by comparing the floating-point model inference with that of the fixedpoint model for a range of bit widths, scanning over a representative sample of input test data. Figure 9 shows that using 14 bits to encode the fractional component of all operands is sufficient to replicate the decision taken by the floating-point model for over 99.5% of the test data. This performance comparison of fixed and floating-point models shows that while nine fractional bits are sufficiently precise to represent the weights, additional precision in the representation of the intermediate sums in the NN calculation is necessary to achieve full performance. Following this, comparisons were performed using a 20bit representation for all internal fixed-point parameters (5 integer bits, 14 fractional bits, and a sign bit) where not explicitly varied.
The degree of design parallelization can be motivated by the number of operations required for each NN infer-2 10 address this is to assign a reuse factor to each NN layer, specifying the number of operations each physical MAC unit may contribute to the set of necessary computations required by that layer. Larger reuse factors will result in a design utilizing fewer FPGA resources at the expense of longer inference latency. For simplicity a single, modellevel reuse factor is considered in the following, where the per-layer reuse factor is given by the greatest common divisor of the reuse factor and the product of input and output multiplicity. Figure 10 demonstrates how the resources required to implement the NN algorithm are affected by the precision to which internal calculations are carried out, and the degree of parallelization specified by the selected reuse factor. In general, the required low-level resources such as ALMs and registers scale linearly with precision to accommodate the widths of increasing data paths. Conversely, a single DSP slice can accommodate a range of operand bit-widths, up to the limit of the design specification at which point a second DSP must be used per calculation. In the case of reuse factor of 128, the largest burden on FPGA resources comes from the required DSPs (either 4 or 8% of the Arria 10 total, depending on the necessary precision), while ALMs are the limiting factor (3-8%) in the case of reuse factor of 1568, where the design parallelization is at most a factor of four in each NN layer. While inference latency depends on the degree of parallelism directly through the reuse factor, it is essentially invariant under changes to the operand precision.
Table VI C compares several implementations for constant precision (20 total bits) and various reuse factors. In general the algorithm latency increases as a function of increasing reuse factor while the numbers of DSPs and BRAMs required are inversely proportional to the reuse factor. Variations in the required registers and ALMs are generally not significant by comparison. These results demonstrate a range of feasible firmware implementations of the algorithm that fit within the resource and latency budget of the GMPS control board. The ability to tune resource usage provides significant flexibility to accommodate potential future scenarios where the NN algorithm will inevitably coexist on a single FPGA with additional control logic that may present inflexible resource constraints of its own.
Up to this point, the discussion of the hardware implementation has centered around the three-hidden-layer MLP architecture found to be performant for the GMPS control problem in the context of RL studies described in Section V. However, the conclusions of the studies described above may be extended to more complex NN algorithms providing improved GMPS performance in tandem with the experience gained through future datataking campaigns.
The simplest extension to the single MLP solution, well-motivated in the context of RL studies, is to run inference with an ensemble of multiple copies of the network in parallel on the FPGA, to improve robustness of performance. Each NN may be programmed with a unique set of weights, allowing for disagreement among the models, where additional voter logic determines the final action to be taken by the control system. This is straightforward to achieve for models with similar complexity to the one studied in Section VI C. Achieving designs that consume ≤6% of all available resources suggests that an ensemble of O(10) models is feasible.
Alternatively, instead of an ensemble of relatively simple models, more complex networks can be pursued. The MLP architecture studied can be extended to additional layers and larger numbers of nodes per layer maintaining an acceptable footprint through corresponding ad-justment of the reuse factor. The theoretical scaling behavior was shown in the calculations of Section VI A and observed in the implementation using Quartus HLS. As an illustrative example, one could consider a refinement of the baseline architecture where the number of nodes per layer is uniformly increased by a scaling factor s. In this case, the number of required multipliers may be kept constant by simultaneously increasing the reuse factor by a factor of s 2 , at the expense of a corresponding increase in algorithm latency. More sophisticated architectures such as convolutional and recurrent NNs may also be considered, taking advantage of their representations as compositions of multiple dense sub-layers. A detailed study of such possibilities is left to future work.
In this paper, we have described a method for controlling the gradient magnet power supply (GMPS), an important subsystem of the Fermilab Booster accelerator, using machine learning models and demonstrated the feasibility of embedding such a model on a fieldprogrammable gate array (FPGA) for a high-uptime, low-latency implementation. We first developed a surrogate LSTM model, based on a recurrent neural network, to reproduce the behaviors of the real GMPS system in the context of the accelerator complex, establishing a safe environment for training reinforcement learning algorithms. Within this environment, we trained a deep Q-network, based on a multilayer perceptron, to choose an optimal action (adjustment of one control knob) to maximize the long-term reward, taken from the negative absolute value of the regulation error (difference between the set and observed values of the minimum GMPS current). We found this surrogate-trained network achieved a factor of 2 improvement over the existing controller in terms of the achieved rewards. Finally, we implemented this network on an Intel Arria 10 FPGA and found it reproduces the CPU-based model, consumes less than 6% of the total FPGA resources, and executes with a latency as low as 2.8 µs, which bodes well for future extensions.
Real-time and operations-hardened solutions will be critical for deploying this technology in an accelerator control context, but we believe a large number of other application spaces will be able to benefit from reinforcement learning on embedded systems. Surrogate models appear promising for supplying the large training data volumes required by reinforcement learning agents. This is particularly important for accelerator facilities where large-scale simulations of the entire complex are absent. Although many open questions remain, this proof-ofprinciple provides confidence to test our proposed concept on "live" hardware. The next steps of this work, including mechanisms for online training and model updates for systems operating with a running accelerator, will be the subject of a future report.
In general, the future for machine learning algorithms