
\section*{Acknowledgements}
This material is based upon work supported by the National Science Foundation Graduate Research Fellowship Program under Grant No. DGE1745016. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation.
\section{Appendix}
\label{appendix01}
\subsection{Additional Dynamical Systems Used}
\subsubsection{Ballistic System}
\label{ss:ballistic}
We also predict trajectories for ballistics: an object is shot out of a cannon in the presence of air resistance. It has a mass and drag coefficient and follows a nearly parabolic trajectory. This system has a two-dimensional state space (altitude $y$ and horizontal range $x$) and 2 parameters (mass and drag coefficient), which we reduce down to one: terminal velocity $v_t$. It is a second order system of differential equations which we reduce to first order using the standard refactoring.

The system is given by
\begin{equation}
\label{eq:ballistics}
\begin{aligned}
\dot{x} &= \dot{x} & 
\ddot{x} &= -\frac{g \dot{x}}{v_t}\\
\dot{y} &= \dot{y} & 
\ddot{y} &= -g \left(1 + \frac{\dot{y}}{v_t}\right)\\
\end{aligned}
\end{equation}
(here, $g$ is the constant of gravitational acceleration).

\subsubsection{Cartpole System}
\label{ss:cartpole}
In section \ref{sec:control} we discuss experiments on a modified Cartpole system with randomly sampled parameters. Here we give a full delineation of that system as defined in \cite{brockman2016openai} with our modifications.

This system has three parameters, one control, and four state variables. The parameters are $l$, the length of the pole, $m_c$, the mass of the cart, and $m_p$, the mass of the pole. The control is $F$, the horizontal force applied to the cart. In the current setup, $u$ can be one of $\pm 10$. The state variables are the lateral position $x$, the angle of the pole from vertical $\theta$, and their first derivatives, $\dot{x}$ and $\dot{\theta}$.

The system is given by the following equations
\begin{equation}
    \label{eq:cartpole}
    \begin{aligned}
    \dot{\theta} &= \dot{\theta} &
    \ddot{\theta} &= \frac{g\sin\theta + \cos\theta\left(\frac{-F - m_p l\dot{\theta}^2 \sin\theta}{m_c + m_p}\right)}{l\left(\frac{4}{3} - \frac{m_p\cos^2 \theta}{m_c + m_p}\right)}\\
    \dot{x} &=\dot{x} &
    \ddot{x} &= \frac{F + m_p l (\dot{\theta}^2 \sin\theta - \ddot{\theta} \cos\theta)}{m_c + m_p}
    \end{aligned}
\end{equation}

We give a reward of $+1$ for every timestep that $|\theta| \leq \frac{\pi}{15}$ and $|x| \leq 2.4$. Initial conditions are uniformly sampled from $[-0.05, 0.05]^4$ at test time and at training $x\sim U([-2, 2])$, $\dot{x}\sim U([-1, 1])$, $\theta \sim U([-0.2, 0.2])$ (which is slightly wider than the $\frac{\pi}{15}$ threshold used at test), and $\dot{\theta}\sim U([-0.1, 0.1])$. The parameters are also uniformly sampled at train and test with $l\sim U([0.6, 1.2])$, $m_c\sim U([0.5, 2])$ and $m_p \sim U([0.03, 0.2])$.

For the partial NDS for Cartpole, we remove the equations corresponding to $\dot{x}$ and $\ddot{x}$.

\subsection{Hyperparameters}

We trained using Adam with a learning rate of $3\times 10^{-3}$ and an ODE relative and absolute tolerance when applicable of $10^{-3}$. This wide tolerance was basically a function of training time as with tighter tolerances experiments took a long time to run and we were more concerned with sample complexity than the tightness of the integration. Hyperparameters were predominantly chosen by trial and error.

Over the course of figuring out how this works and then evaluating the models we were evaluating, we ran some form of this code 347 times. The typical setup was either a 1080Ti GPU and 6 CPU cores or 7 CPU cores for roughly a day per experiment. We found that most of these trainings were marginally faster on the GPU ~1.5x speedup, but weren't religious about the GPU as we had many CPU cores and many experiments to run in parallel.

\paragraph{Lorenz Initial Conditions and Parameters}
For our Lorenz system, we sampled $\rho \sim U([15, 35])$, $\sigma\sim U([9, 12])$, $\beta \sim U([1, 3])$, $x_0\sim U([0, 5])$, $y_0\sim U([0, 5])$, $z_0\sim U([0, 5])$.

\paragraph{Ballistic Initial Conditions and Parameters}
For our Ballistic System, we sampled masses $m \sim U([1, 100])$, drag coefficients $c_d\sim U([0.4, 3])$, $x_0 \sim U([-100, 100])$, $y_0\sim U([0, 200])$. We then use $v_t = \frac{mg}{c_d}$ to recover the terminal velocity used in our model.
\subsection{Fusion Model Parameters}
\label{fusion_model_parameters}
$n_i$ is ion density (which we approximate as a constant value of $5\times 10^{19}$ deuterium ions per $m^3$), $m_i$ is ion mass (which we know since our dataset contains deuterium shots and the mass of a deuterium ion is $3.3436\times 10^{-27}$ kg), and $R_0$ is the tokamak major radius of $1.67$ m. 
\subsection{Fusion Data Preprocessing}
Data is loaded and partially processed within the OMFIT framework \cite{meneghini_integrated_2015}. A new module, ``SIGNAL\_PROCESSING", has been developed for this task and is publicly available on the ``profile\_prediction\_data\_processing" branch of the OMFIT source code. The rest of the processing is done on Princeton's Traverse computing cluster, and is available in the GitHub sourcecode for this project (https://github.com/jabbate7/plasma-profile-predictor). 

DIII-D shots from 2010 through the 2019 campaign are collected from the MDS+ database. Shots with a pulse length less than 2s, a normalized beta less than 1, or a non-standard topology are excluded from the start. A variety of non-standard data is also excluded, including the following situations:
\begin{enumerate}
    \item during and after a dudtrip trigger (as determined by the pointname "dustripped") 
    \item during and after ECH (pointname "pech" greater than .01) activation, since ECH is not currently included as an actuator
    \item whenever density feedback is off (as determined by the pointname "dsifbonoff")
    \item during and after non-normal operation of internal coil, with an algorithm described by Carlos Paz-Soldan
\end{enumerate}


All signals are then put on the same 50ms time base by averaging all signal values available between the current time step and 50ms prior. If no data is available in the window, the most recent value is floated. Only time steps during the "flattop" of the plasma current are included. 
The flattop period is determined by DIII-D's ``t\_ip\_flat" and ``ip\_flat\_duration" PTdata pointnames. 

All profile data is from Zipfit \cite{meneghini_integrated_2015}. The profiles are linearly interpolated onto 33 equally spaced points in normalized toroidal flux coordinates, denoted by $\rho$, where $\rho=0$ is the magnetic axis and $\rho=1$ is the last closed flux surface.

\subsection{Computation of Noise}
\label{computation_of_noise}
In our previous experiments, we did not add noise to the trajectories generated by the synthetic systems. We generate noise added to the trajectories by first sampling a set of 100 trajectories and computing for the $i$th component of the state space the RMS value $c_i$ across our 100 trajectory samples. Then we sample noise $\mathcal{N}(0, rc_i)$ from a normal distribution where the variance is $c_i$ scaled by our `relative noise' term $r$. We vary $r$ to control the amount of noise added to a system in a way that generalizes across problems and across components of a problem's state space.
\subsection{Robustness to irregular timesteps}
\label{a:jitter}
We also explored how these models respond to data sampled in intervals that are not regularly spaced in a fashion discussed in \ref{computation_of_jitter}. 
In Figure \ref{fig:synthetic_jitter}, we see that the ODE-based models are able to handle the irregularly spaced data much better than the others. Like with noise, we are focusing on relative performance but even so, the full NDS even does substantially better under large jitter settings. As we have discussed, this makes sense because these models natively handle time. We conjecture that the Full NDS neural networks may learn something slightly more general by this `domain randomization' given that they are correctly specified models that receive all the information of the differing timesteps. The symbolic regression method fails under jitter, presumably because it relies heavily on finite differencing. 
\begin{figure}[t]
    \centering
    \includegraphics[width=0.48\columnwidth]{figures/ballistic_jitter.pdf}
    \hfill
    \includegraphics[width=0.48\columnwidth]{figures/lorenz_jitter.pdf}
    \caption{\textbf{Change in performance when data are irregularly sampled in time.} The Full NDS performs better than any comparison model under added jitter.}
    \label{fig:synthetic_jitter}
\end{figure}
\subsection{Robustness to added noise}
\label{noise}
We evaluate our experiments on additive noise scaled relative to the system as discussed in \ref{computation_of_noise}. As the NDS does substantially better than the other models on the synthetic data, we look at the effect of noise relative to the original performance to focus on the properties of the new model. When a large amount of noise is added to the system, NDS's performance degrades faster than that of the other models, though it still outperforms them on an absolute basis. We think full and partial NDS might be unstable here due to errors propagating through the prior knowledge model. The other models all otherwise perform fairly similarly under added noise.
\begin{figure}[t]
    \centering
    \includegraphics[width=.48\columnwidth]{figures/ballistic_noise.pdf}
    \hfill
    \includegraphics[width=.48\columnwidth]{figures/lorenz_noise.pdf}
    \caption{\textbf{Change in performance under added relative noise $r$.} The NDS seem to have more variance in performance under high amounts of noise as the comparison models but it still greatly outperforms them on an absolute basis.}
   
    \label{fig:synthetic_noise}
\end{figure}
\subsection{Computation of Jitter}
\label{computation_of_jitter}
Though, as we will discuss, our fusion data in section \ref{fusion_experiments} has been postprocessed to be regularly spaced, in practice, data on tokamaks and in many other systems come in irregularly and are sometimes missing. As both the fully connected neural ODE and NDS models are integrated in continuous time, they can handle arbitrarily spaced data. For the LSTM and Fully Connected models, we concatenated the times of the datapoints to the associated data and fed this into the model. For each batch of training and test data and some value of `jitter', $j$, we create a new time series $\{t_i + j_i\}_{i=1}^T$, where $j_i \sim U([-j, j])$. Since our timestep for synthetic experiments is $0.5$s we try values of $j$ of $0.1$s and $0.25$s. We then generate the batch of trajectories by integrating our systems to the new timesteps.

\begin{table}[]
\centering
\resizebox{\textwidth}{!}{
\begin{tabular}{|c|cccc|cccc|}
\hline
System & \multicolumn{4}{c}{Lorenz} & \multicolumn{4}{c}{Ballistic} \\\hline
Samples &   100,000  &  25,000   &   5,000  &  1,000  &  100,000   &   25,000  & 5,000    & 1,000  \\\hline
Full NDS &   $36.4\pm 2.8$  & $38.3\pm 2.8$& $\mathbf{40.3 \pm 4.9}$ & $\mathbf{57.8\pm 5.5}$ & $\mathbf{1110\pm 497}$  & $1267.81 \pm 531$    & $\mathbf{1308.9\pm 394}$ & $\mathbf{2349\pm 278}$   \\
Partial NDS & n/a &  n/a  & n/a   & n/a &    $1249\pm 324$  &    $\mathbf{1194\pm293}$  & $1434.8\pm426$  & $2378\pm 593$   \\
NDS0 & $\mathbf{32.5 \pm3.8}$ & $\mathbf{35.4\pm 4.3}$ & $\mathbf{41.3 \pm 4.3}$ & $\mathbf{55.4\pm 6.4}$ & $\mathbf{1034 \pm434}$ & $\mathbf{1128 \pm 453}$ & $1624\pm 484$ & $2833 \pm 233$\\
GBO &  $233.8\pm 15$   & n/a &  n/a  &  n/a &  $160990 \pm 1261$   &   n/a  & n/a    &  n/a \\\hline
\end{tabular}
}
\caption{Parameter Error Results as discussed in Section \ref{sec:synthetic_experiments} and Figure \ref{fig:parameters}.}
\end{table}
\begin{table}[]
\centering
\begin{tabular}{|c|cc|cc|}
\hline
System & \multicolumn{2}{c}{Lorenz} & \multicolumn{2}{c}{Ballistic} \\
\hline
Jitter Level (s) &    0.1      &      0.25    &   0.1    &    0.25  \\\hline
FC &     $0.104$      &    $0.109$   & $0.106$          &  $0.125$        \\
FC NODE & $0.008$     & $0.012$    &  $0.075$    &  $0.077$   \\
Full NDS &   $-0.010$   &  $\mathbf{-0.066}$    &    $\mathbf{-0.012}$ & $\mathbf{-0.096}$ \\
Partial NDS& $\mathbf{-0.011}$  &  $-0.012$   & $0.024$  &    $0.016$\\
LSTM  & $0.032$  & $0.031$  &  $0.130$  &   $0.132$  \\
SR &   $2.69$   & $2.17$  &   $0.239$   &  $0.342$   \\
GBO &  $-0.002$ &  $0.058$ & $0.002$  &  $0.013$ \\\hline    
\end{tabular}
\caption{The effect of jitter on a given algorithm for a given level of jitter as computed in \ref{computation_of_jitter}. The effect is calculated by dividing the performance with jitter by the mean performance without and subtracting $1$.}
\end{table}

\begin{table}[]
\centering
\begin{tabular}{|c|c|}\hline
Model & $L2$ Error on the Fusion Test Set \\\hline
FC &  $1.942\times 10^{11}$\\
FC NODE & $8.288 \times 10^{10}$ \\
NDS with Approximate Dynamics &  $\mathbf{4.837 \times 10^{10}}$\\
NDS0 & $8.944 \times 10^{10}$\\
LSTM &  $2.527\times 10^{11}$\\
SR & $2.546 \times 10^{12}$ \\
GBO & $6.452\times 10^{11}$ \\\hline
\end{tabular}
\caption{The performance of our comparison models on the nuclear fusion problem, as discussed in Section \ref{fusion_experiments}.}
\end{table}
\begin{table}[]
\centering
\begin{tabular}{|c|cc|cc|}
\hline
 & \multicolumn{2}{c|}{MSE of Model} & \multicolumn{2}{c|}{Performance of MPC} \\
 \hline
\# of Training Examples & 100K & 25K &100K & 25K \\\hline
FC & $\mathbf{0.016\pm 0.004}$  &$0.022\pm 0.006$ &$\mathbf{86\pm 4}$ & $76\pm 4$ \\
FC NODE & $\mathbf{0.015\pm 0.008}$ &$\mathbf{0.019 \pm 0.007}$ &  $\mathbf{87\pm 4}$  & $\mathbf{78\pm 4}$  \\
LSTM &   $0.081\pm 0.023$ & $0.092 \pm 0.025$   & $33\pm 6$ & $28\pm 5$        \\
Full NDS &  $\mathbf{0.015 \pm 0.007}$ & $\mathbf{0.017 \pm 0.002}$  & $\mathbf{86\pm 3}$ &$\mathbf{82\pm 5}$  \\
Partial NDS & $\mathbf{0.016 \pm 0.011}$  & $\mathbf{0.017 \pm 0.005}$     & $79\pm 6$ & $78 \pm 5$ \\
NDS0 & $\mathbf{0.013\pm 0.009}$ & $\mathbf{0.017\pm0.008}$ & $76 \pm 5$ & $74 \pm 6$\\
SR &   $0.032\pm 0.023$  &     $0.034 \pm 0.017$& $68\pm 4$ & $66\pm 3$ \\
\hline
\end{tabular}
\caption{Modeling and Control on the Cartpole system with large amounts of data.}
\label{t:cartpole_big_data}
\end{table}

\section{Conclusion}
\label{conclusion}
\vspace{-1mm}
In conclusion, we give a framework that merges theoretical dynamical system models with deep learning by backpropagating through a numerical ODE solver. This framework succeeds even when there is a partial or approximate model of the system. We show there is an empirical reduction in sample complexity and increase in accuracy on two synthetic systems and on a real 
nuclear fusion dataset.
In the future, we wish to expand upon our work to make more sophisticated models in the nuclear fusion setting as we move toward practical use. We also hope to explore applications of this framework in other area which have ODE-based models of systems.



\section{Experiments}
\label{experiments}







In the following experiments, we aim to show that our methods improve predictions of physical systems by including prior dynamical knowledge. These improvements hold even as we vary the configurations between structured and fluid settings. We show that our models learn from less data and are more accurate, that they handle irregularly spaced data well, and that they learn the appropriate parameters of the prior knowledge systems even when they only ever see trajectories. 




We use L2 error as our evaluation measure for predictive accuracy as given by Equation \ref{eq:loss}, though in cases where the absolute errors aren't very interpretable, we normalize for ease of comparison.
We also evaluate our model's ability to predict the system parameters by computing the L2 error, i.e. $\sum_{i=1}^n ||\hat{\phi}_i - \phi_i||_2^2$.
For synthetic examples, we consider the Lorenz system  in \eqref{eq:lorenz} and a simple Ballistic system modeling projectile motion under a variety of drag conditions. 
We learn over trajectories $\{(x_{t_i}, u_{t_i}, t_i)\}_{i=1}^T$ where the $x_{t_i}$ are generated by numerically integrating $\dot{x}_\phi(x, u, t)$ using scipy's odeint function \cite{virtanen_scipy_2019}, with $x_0$ and $\phi$ uniformly sampled from $\mathcal{X}$ and $\Phi$, and $u_{t_i}$ given. Note that $\phi$ remains fixed throughout a single trajectory. Details on the ranges of initial conditions and parameters sampled are in the appendix. We evaluate the synthetic experiments on a test set of 20,000 trajectories that is fixed for a particular random seed generated in the same way as the training data.
We use a timestep of $0.5$ seconds for the synthetic trajectories.
On the Ballistic system this allows us to see trajectories that do not reach their peak and those that start to fall.
Since the Lyapunov exponent of the Lorenz system is less than $3$, in 16 predicted timesteps we get both predictable and unpredictable data \cite{froyland_lyapunov-exponent_1984}).
We believe it is important to look at the progress of the system across this threshold to understand whether the NDS model is robust to chaotic dynamics --- since the Lorenz system used for structure is itself chaotic, we want to make sure that the system does not blow up over longer timescales. 

We note that ReLU activations were chosen for all feedforward and recurrent architectures, while in the Neural-ODE-based architectures, we follow the recommendations of \cite{chen_neural_2018} and use the Softplus. The sizes and depths of the baselines were chosen after moderate hyperparameter search.


We can view the Partial NDS and NODE as ablations of the Full NDS model which remove some and all of the prior knowledge, respectively.
Each model takes 32 timesteps of state 
and control information as input and are trained on predictions for the following 16 timesteps. The ODE-based models are integrated from the initial conditions of the last given state.
All neural networks are all trained with a learning rate of $3\times 10^{-3}$, which was seen to work well across models. We generated a training set of 100,000 trajectories, test set of 20,000 trajectories, and validation set of 10,000 trajectories. Training was halted if validation error did not improve for 3 consecutive epochs.
\vspace{-2mm}
\subsection{Comparison Methods}
\label{a:comparisons}
We compare our models with other choices along the spectrum of structured to flexible models from both machine learning and system identification.
In our paper, we compared the following methods in our experiments:
\begin{itemize}
    \item \textbf{Full NDS}: A Neural Dynamical System with the full system dynamics for the problem being analyzed. The full construction of this model is given by Equation \ref{eq:full_dynamics}. For the functions $h_{\theta}, c_{\vartheta}, d_{\tau}$, we use fully connected networks with 2 layers, Softplus activations, 64 hidden nodes in each layer, and batch normalization.
    \item \textbf{Partial NDS}: A Neural Dynamical System with partial system dynamics for the problem being analyzed. These follow Equation \ref{eq:partial_dynamics} as applied to Equation \ref{eq:full_dynamics}.
    For the Ballistic system, we only provide equations for $\dot{x}$ and $\ddot{x}$, excluding the information about vertical motion from our network. For the Lorenz system, we only provide equations for $\dot{x}$ and $\dot{y}$, excluding information about motion in the $z$ direction. For the Cartpole system, we only provide information about $\dot{\theta}$ and $\ddot{\theta}$. These equations were chosen somewhat arbitrarily to illustrate the partial NDS effectiveness. We use similar neural networks here as for the Full NDS.
    \item \textbf{NDS0}: A Full NDS with residual terms removed. This serves as an ablation which shows the use of the residual terms.
    \item \textbf{Fully Connected} (FC): A Fully-Connected Neural Network with $4$ hidden layers containing $128$ nodes with ReLU activations and batch normalization.
    \item \textbf{Fully Connected Neural ODE} (FC NODE): A larger version of the Neural ODE as given in \cite{chen_neural_2018}, we use $3$ hidden layers with $128$ nodes, batch norm, and Softplus activations for $\dot{x}$. This can be interpreted as a version of our NDS with no prior knowledge, i.e. $g(x) = 0$.
    \item \textbf{LSTM}: A stacked LSTM with 8 layers as in \cite{graves_generating_2013}. The data is fed in sequentially and we regress the outputs of the LSTM against the true values of the trajectory.
    \item \textbf{Gray Box Optimization} (GBO): We use MATLAB's gray-box system identification toolbox \cite{ljung_developments_2009} along with the prior knowledge ODEs to fit the parameters $\hat{\phi}$ as an alternative to using neural networks. This algorithm uses trust-region reflective nonlinear least squares with finite differencing \cite{coleman_interior_1996} to find the parameter values which minimize the error of the model rollouts over the observed data.
   
    \item \textbf{Sparse Identification of Nonlinear Systems} (SR): We use the method from \cite{brunton_discovering_2015} to identify the dynamical systems of interest. This method uses sparse symbolic regression to learn a linear mapping from basis functions of the state $x_t$ and control $u_t$ to the derivatives $\dot{x}_t$ computed by finite differences. Our synthetic systems are in the span of the polynomial basis that we used.
    \item \textbf{APHYNITY}: We use the method from \cite{guen2021augmenting}, which fits a min-error parameter but then has an additional neural network component to model unknown dynamics.
\end{itemize}
\subsection{Synthetic Experiments}
\label{sec:synthetic_experiments}
\vspace{-1mm}
We first present results on a pair of synthetic physical systems where the data is generated in a noiseless and regularly spaced setting. 


\paragraph{Sample Complexity and Overall Accuracy}
In order to test sample complexity in learning or fitting, we generated data for a full training dataset of 100,000 trajectories. We then fit our models on different fractions of the training dataset: $1, 0.25, 0.05, 0.01$. We repeated this process with $5$ different random seeds and computed the $L2$ error of the model over a the various dataset splits seen by the model in Table \ref{table_of_results}. The error regions are the standard error of the errors over the various seeds.

We also see that with small amounts of data, the NDS models greatly outperform the Neural ODE, but with the full dataset, their performances get closer. This makes sense as the Neural ODE is likely able to infer the structure of the system with large amounts of data. Also, the Fully Connected Neural ODE outperforms the other baselines, which we posit is due to the fact that it implicitly represents that this system as a continuous time dynamical process and should change in a continuous fashion. From a sample-complexity perspective it makese sense that the better initialization of NDS should matter most when data is limited. A table of the full results of all experiments can be seen in Table \ref{table_of_results}.

We notice that the NDS0 slightly outperforms the NDS with higher variance on these systems. Since it has a perfect model of the system, the residual components aren't necessary for the model to perform well, however, there is no way the network can `correct' for a bad estimate. 

Curiously, we see on the ballistic system that the partial NDS slightly outperforms the full NDS in the small data setting, but they converge to similar performance with slightly more data.  
A potential explanation for this is that errors propagate through the dynamical model when the parameters are wrong, while the partial systems naturally dampen errors since, for example, $\dot{z}$ only depends on the other components through a neural network. Concretely this might look like a full NDS predicting the wrong Rayleigh number $\sigma$ which might give errors to $y$ which would then propagate to $x$ and $y$. Conversely, this wouldn't happen as easily in a partial NDS because there are neural networks intermediating the components of the system. We also see APHYNITY (which fits a min-error parameter and is otherwise very similar to NDS) performs worse than NDS in this setting.


\paragraph{Parameter Learning without Explicit Supervision}
\label{a:parameters}
For experiments in Figure 2, we stored the parameter estimates $\hat{\phi}$ for the NDS and gray box models and compared them to the true values to see how they perform in identification rather than prediction. None of these models were ever supervised with the true parameters. We see in Figure \ref{fig:parameters} that the NDS is better able to estimate the parameter values than the gray-box method for both systems tested.  We believe this is  because  our  method  is  able  to  leverage  many trajectories to infer  the  parameters whereas the gray-box method only uses the single trajectory.
\begin{figure}
    \includegraphics[width=0.48\columnwidth]{figures/ballistic_params.png}
    \hfill
    \includegraphics[width=0.48\columnwidth]{figures/lorenz_params.png}
    
    \caption{\textbf{$L2$ distance between $\phi$ and $\hat{\phi}$.} As the NDS are trained under the usual $L2$ supervision, the parameters $\hat{\phi}$ of the system approach the correct values.} 
   
\label{fig:parameters}
\end{figure}


\begin{table*}[t]
\centering
\resizebox{\textwidth}{!}{
\begin{tabular}{|c|cccc|cccc|}
\hline
System & \multicolumn{4}{c}{Lorenz} & \multicolumn{4}{c}{Ballistic} \\
\hline
Samples & 100,000 & 25,000 & 5,000 & 1,000  & 100,000 & 25,000 & 5,000 & 1,000 \\\hline
FC &  $1.06 \pm 0.002$ & $1.39 \pm 0.003$ & $1.694 \pm 0.002$ &  $4.692 \pm 0.6$ & $7.8 \pm 1.4$   & $8.3 \pm 1.6$  & $9.2\pm 2.5$ &  $13.4\pm 3.4$ \\
FC NODE & $1.205 \pm 0.01$ & $1.233 \pm 0.01$ & $1.27 \pm 0.01$ & $1.917 \pm 0.12$ & $2.53 \pm 0.2$ & $2.6 \pm 0.4$ & $4.8 \pm 0.7$ & $9.7 \pm 1.3$\\
Full NDS & $\mathbf{1.004\pm 0.06}$ & $\mathbf{1.087\pm 0.06}$ & $1.14 \pm 0.05$ & $1.42\pm .05$ & $1.2 \pm 0.1$ & $1.4\pm 0.2$ & $1.5\pm0.2$  & $4.23\pm 0.3$  \\
Partial NDS & $1.036\pm 0.03$  &$\mathbf{1.064 \pm 0.1}$ &  $\mathbf{1.12 \pm 0.06}$   & $\mathbf{1.39\pm 0.04}$  & $\mathbf{1.05 \pm 0.1}$ & $\mathbf{1 \pm 0.03}$ & $\mathbf{1.48\pm 0.02}$ & $\mathbf{1.9\pm 0.15}$ \\
NDS0 & $\mathbf{1 \pm 0.03}$ & $\mathbf{1.075\pm0.1}$ & $\mathbf{1.13\pm 0.11}$ & $\mathbf{1.36 \pm 0.17}$ & $\mathbf{1.2\pm 0.06}$ & $\mathbf{1.35 \pm 0.14}$ & $\mathbf{1.45\pm 0.18}$ & $4.3\pm 0.6$\\
APHYNITY & $1.08 \pm 0.03$ & $1.17 \pm 0.05$ & $1.23 \pm 0.04$ & $1.61 \pm 0.07$ &$1.7\pm 0.1$ & $1.75 \pm 0.13$ & $1.94\pm 0.18$ & $6.2\pm 0.4$\\
LSTM &  $4.98\pm 0.01$   & $5.98 \pm 0.3$ & $5.99 \pm 0.3$ & $6.13 \pm 0.4$  &   $8.5\pm 1.6$  &  $9.2 \pm 1.5$  & $10.1 \pm 2.1$ & $14.9\pm 1.9$ \\
SR & $2.3 \pm 0.6$ & n/a & n/a & n/a & $3.5\pm 0.3$  & n/a & n/a & n/a \\
GBO & $2.8\pm0.4$ & n/a & n/a & n/a & $2.94 \pm 0.3$ & n/a & n/a & n/a  \\
\hline
\end{tabular}
}
\caption{Sample Complexity Results as discussed in Section \ref{sec:synthetic_experiments}. Here, the values are normalized by the smallest reported value for comparison purposes.}
\label{table_of_results}
\vspace{-9mm}
\end{table*}

\vspace{-2mm}

   
   
   
   
   
   
   
   
   






\subsection{Fusion Experiments}
\label{fusion_experiments}
We explored the concept of approximate system dynamics in a simplified fusion system. 
We predict the state of the tokamak as summarized by its stored energy and rotational frequency given the time series of control input in the form of injected power and torque. As mentioned in Section \ref{methods}, we have a simplified physical model given by Equation \ref{eq:tokamak_model} that approximately gives the dynamics of these quantities and how they relate to one another through time. 
Though there is a lot of remaining work to apply this model in a real experiment, approaches merging theoretical models with data to make useful predictions can be embedded into useful controller designs and improve the state of fusion.








Our full dataset consisted of 17,686 trajectories, which we randomly partitioned into 1000 as a test set and 16,686 as a training set.\footnote{Data is loaded and partially processed within the OMFIT framework \cite{meneghini_integrated_2015}. We used the ``SIGNAL\_PROCESSING" module which has recently been developed for this task and is publicly available on the ``profile\_prediction\_data\_processing" branch of the OMFIT source code. Details of the preprocessing are in the Appendix.}
The data are measured from the D-IIID tokamak via magnetic reconstruction \cite{ferron_real_1998} and charge-exchange recombination spectroscopy \cite{haskey_active_2018}. Similar to our synthetic experiments, we cut each trajectory into overlapping 48 timestep sections and train on 32 timesteps to predict 16 timesteps.  We compare with the same models as in the previous section, but using our Fusion Neural Dynamical System as described in Equation \ref{eq:full_dynamics} with $g$ given by Equation \ref{eq:tokamak_model}. As we discussed above, the dynamics in this equation are approximate.
To illustrate this, we have included the accuracy of the naive dynamics with no learning on our data with fixed confinement times $\tau_e = \tau_m = 0.1$s as the Nominal Fusion Model in Table \ref{fusion_table}. We use a larger fully connected network with 6 layers with 512 hidden nodes to attempt to capture the added complexity of the problem.


\begin{table}[]
\begin{tabular}{|c|c|}\hline
Model & $L2$ Error on the Fusion Test Set \\\hline
FC &  $4.02 \pm 0.27$\\
FC NODE & $1.71 \pm 0.11$ \\
Nominal Fusion Model & $2.89$\\
NDS with Approximate Dynamics &  $\mathbf{1 \pm 0.06}$\\
NDS0 & $1.85 \pm 0.09$\\
APHYNITY & $1.21 \pm 0.15$\\
LSTM &  $5.23 \pm 0.43$\\
SR & $5.26 \pm 0.35$ \\
GBO & $2.98$ \\\hline
\end{tabular}
\caption{The performance of our comparison models on the nuclear fusion problem, as discussed in Section \ref{fusion_experiments}. We again normalize by the smallest value for ease of comparison. We use the same data but different random seeds for the standard errors.}
\label{fusion_table}
\end{table}



\paragraph{Sample Complexity and Overall Accuracy}



When comparing our NDS models, the machine learning baselines, the system ID baselines, and a nominal model from \cite{boyer_real-time_2019}, we see that the Fusion NDS model performs best by a large margin. Although the fully connected neural ODE performs competitively, it fails to reach the same performance. We speculate that the dynamical model helps with generalization whereas the fully connected network may overfit the training data and fail to reach good performance on the test set. Here the NDS0 is unable to perform well compared to the NDS, as the approximate dynamics mean that the model error is somewhat catastrophic for predictions. We see however that the NDS0 outperforms the Nominal Physics Model as it is able to estimate the parameters for each example rather than fixing values of the parameters for the whole dataset. Similarly, APHYNITY outperforms the nominal model but underperforms the NDS0, suggesting that the online parameter estimation component is necessary for good performance. 

We see these results as highly encouraging and will continue exploring uses of NDS in fusion applications.
\subsection{Control Experiment}
\label{sec:control}
We also explored the use of these models for control purposes using model-predictive control \cite{camacho2013model}. For this purpose, we modified the Cartpole problem from \cite{brockman2016openai} so that there are a variety of parameter values for the weight of the cart and pole as well as pole length. Typically, a `solved' cartpole environment would imply a consistent performance of 200 from a control algorithm. However, there are three factors that make this problem more difficult. First, in order to allow each algorithm to identify the system in order to make appropriate control decisions, we begin each rollout with 8 random actions. The control never fails at this point but would certainly fail soon after if continued. Second, the randomly sampled parameters per rollout make the actual control problem more difficult as the environment responds less consistently to control. For example, MPC using the typical Cartpole environment as a model results in rewards of approximately 37. Third, all training data for these algorithms uses random actions with no exploration, which has been seen to degrade the performance of most model-based RL or control algorithms \cite{random_exploration}.

We then trained dynamics models on this `EvilCartpole' enviroment for each of our comparison algorithms on datasets of trajectories on the environment with random actions. At that point, we rolled out trajectories on our EvilCartpole environment using MPC with control sequences and random shooting with 1,000 samples and a horizon of 10 timesteps. The uncertainties are standard errors over 5 separately trained models.

As shown in Table \ref{t:cartpole}, the NDS algorithms outperform all baselines on the cartpole task for both the modeling and control objectives. We see that all algorithms degrade in performance as the amount of data is limited. We notice however that with larger amounts of data (we performed other experiments with 25,000 and 100,000 samples) the Fully Connected and Neural ODE models perform as well as the NDS models. We hypothesize that this is due to the fact that the cartpole dynamics are ultimately not that complicated and with sufficient data unstructured machine learning algorithms can learn the appropriate dynamics to reach a modestly performing controller as well as NDS. 



\begin{table}[]
\centering
\begin{tabular}{|c|cc|cc|}
\hline
 & \multicolumn{2}{c|}{MSE of Model} & \multicolumn{2}{c|}{MPC Returns} \\
 \hline
\# Train &  5K & 1K & 5K & 1K\\\hline
FC &  $0.031 \pm 0.009$    & $0.058\pm 0.018$ & $52\pm 3$ & $41\pm 4$\\
FC NODE &  $0.028\pm 0.011$ & $0.049 \pm 0.013$ &   $55\pm 4$&  $46\pm 3$  \\
LSTM &   $0.081\pm 0.023$ & $0.092 \pm 0.025$ & $23\pm 6$      & $25 \pm 8$       \\
Full NDS & $\mathbf{0.020 \pm 0.006}$ & $\mathbf{0.029 \pm0.007}$   &  $\mathbf{72\pm 4}$& $\mathbf{60\pm 3}$  \\
Partial NDS & $\mathbf{0.022 \pm 0.009}$   &$\mathbf{0.033\pm 0.011}$     & $\mathbf{69\pm 8}$     & $55\pm 6$     \\
NDS0 & $\mathbf{0.023 \pm 0.013}$ & $\mathbf{0.028 \pm 0.014}$ & $\mathbf{71 \pm 11}$ & $\mathbf{57 \pm 8}$\\
SR &  $0.037 \pm 0.023$ & $0.041\pm 0.015$  &  $65 \pm 4$ & $56 \pm 4$ \\
GBO & $0.046 \pm 0.019$    &    n/a  &  $49\pm 5$ & n/a   \\
\hline
\end{tabular}
\caption{Modeling and Control on the EvilCartpole system.}
\label{t:cartpole}
\end{table}




\section{Introduction}
\label{introduction}
\vspace{-2mm}


The use of function approximators for dynamical system modeling has become increasingly common. This has proven quite effective when a substantial amount of real data is available relative to the complexity of the model being learned \cite{chua_deep_2018, janner_when_2019, chen_non-linear_1990}. These learned models are used for downstream applications such as model-based reinforcement learning \cite{nagabandi_neural_2017, ross_agnostic_2012} or model-predictive control (MPC) \cite{wang_exploring_2019}.

Model-based control techniques are exciting as we may be able to solve new classes of problems with improved controllers. Problems like dextrous robotic manipulation \cite{nagabandi_deep_2019}, game-playing \cite{schrittwieser_mastering_2019}, and nuclear fusion are increasingly being approached using model-based reinforcement learning techniques. However, learning a dynamics model using, for example, a deep neural network can require large amounts of data. 
This is especially problematic when trying to optimize real physical systems, where data collection can be expensive.
As an alternative to data-hungry machine learning methods, there is also a long history of fitting models to a system using techniques from system identification, some of which include prior knowledge about the system drawn from human understanding \cite{nelles_nonlinear_2013, ljung_developments_2009, sohlberg_grey_2008}. 
These models, especially in the gray-box setting, are typically data-efficient and often contain interpretable model parameters. However, they are not well suited for the situation where the given prior knowledge is approximate or incomplete in nature. They also do not generally adapt to the situation where trajectories are drawn from a variety of parameter settings at test time. This is an especially crucial point as many systems of interest exhibit path-dependent dynamics, 
which we aim to recover on the fly.

In total, 
system identification methods are sample efficient but inflexible given changing parameter settings and incomplete or approximate knowledge. Conversely, deep learning methods are more flexible at the cost of many more samples. 
In this paper, we aim to solve both of these problems by biasing the model class towards our physical model of dynamics. Physical models of dynamics are often given in the form of systems of ordinary differential equations (ODEs), which are ubiquitious and may have 
free parameters that specialize them to a given physical system.
We develop a model that uses neural networks to predict the free parameters
of an ODE system from the previous timesteps as well as residual terms added to each component of the
system. To train this model, we integrate
over the ODE and backpropagate gradients from
the prediction error. This particular combination of prior knowledge
and deep learning components is effective in quickly learning the dynamics and allows us to adjust system behavior in response to a wide variety of dynamic parameter settings. Even when the dynamical system is partially understood and only a subset of the ODEs are known, we find that our method still enjoys these benefits.
We apply our algorithm to learning models in three synthetic settings: a generic model of ballistics, the Lorenz system \cite{lorenz_deterministic_1963}, and a generalized cartpole problem, which we use for control as well. We also learn a high-level model of plasma 
dynamics for a fusion tokamak from real data.


The contributions of this paper are
\begin{itemize}
    \item We introduce Neural Dynamical Systems (NDS), a new class of model for learning dynamics that can incorporate prior knowledge about the system.
    \item We show that these models naturally handle the issue of partial or approximate prior knowledge, irregularly spaced data, and system dynamics that change across instantiations, which generalizes the typical system identification setting. We also show that these advantages extend to control settings. 
   
   
    \item We demonstrate this model's effectiveness on a real dynamics problem relevant to nuclear fusion and on synthetic problems where we can compare against a ground truth model.
   
\end{itemize}
\section{Methods}
\label{methods}
\vspace{-1.5mm}

We build up the description of our proposed method by first describing the two methods that inspire it: gray box system identification through optimization \cite{ljung_developments_2009}, and using a Neural ODE \cite{chen_neural_2018} to predict future states in a dynamical system.

To apply grey box optimization \cite{ljung_developments_2009} to a dynamical system $\dot{x} = f_\phi(x, u, t)$ for problem data $\{x_{t_i}\}_{t=0}^{T'}$, we would use nonlinear least squares \cite{coleman_interior_1996} to find 
\begin{equation}\label{eq:gray_box}
\begin{aligned}
\hat{\phi} = \argmin_{\hat{\phi}}\sum_i \left|\left| \int_0^{t_i}f_{\phi}(x, u(t), t)\,dt -\right.\right.\\\left.\left. \int_0^{t_i}f_{\hat{\phi}}(x, u(t), t)dt\right|\right|.
\end{aligned}\end{equation}
This makes good use of the prior knowledge component of our system but is prone to compounding errors through integration and does not leverage data that may have come from alternate system parameters.

A data driven approach would be to  minimize the same objective with a fully connected neural ODE \cite{chen_neural_2018} $h_\theta$ in place of $f$.
However, we find that this procedure requires large amounts of training data and doesn't leverage any prior knowledge we might have, though it is flexible to classes of dynamical systems.

We define a Neural Dynamical System by taking the advantages of both these methods in the setting where we know the full and correct ODEs and then show how to generalize it to situations where only some ODEs are known or they are approximate. Specifically, a \emph{Neural Dynamical System} (NDS) is a class of dynamical systems where a neural network predicts some part of $f_\phi(x, u, t)$, usually parameters $\phi$ or a term which is added to $f$.
























\paragraph{NDS with Full System Dynamics}

Consider a class of dynamical systems as defined in Section \ref{problem_setting} where $x\in \mathbb{R}^n$, $u\in \mathbb{R}^m$, $\phi\in \mathbb{R}^{d_p}$, $d_h, d_c \in \mathbb{N}$ and let $\theta$, $\vartheta$, $\tau$ be trainable neural network weights. Let $h_\theta(x_{t_{1:T'}}, u_{t_{1:T}})$ be a neural net mapping state history and control sequence to the $d_p$ parameters of the system $\hat{\phi}$ and an embedding $b_h \in \mathbb{R}^{d_h}$. Also let $c_\vartheta(x_t, u_t)$ be a similar network taking a single state and control that outputs an embedding $b_c \in \mathbb{R}^{d_c}$. Finally, let $d_\tau(b_h, b_c)$ be a network which takes the two output embeddings from the previous network and outputs residual terms $\hat{r}$.
Intuitively, we would like to use the observed history to estimate our system parameters, and some combination of the observed history and current observation to estimate residuals, which influences the design of our model, the neural dynamical system 
(a visualization of which is shown in Figure~\ref{fig:arch}), written
\begin{equation}
\begin{aligned}
\dot{x} & = \underbrace{g_{\hat{\phi}}(x_t, u_t, t)}_{\text{Prior knowledge}} + \hat{r} &
\hat{\phi}, b_h &= \underbrace{h_\theta(x_{t_{1:T'}}, u_{t_{1:T}})}_{\text{History encoder}} & \\
b_c &= 
\underbrace{c_\vartheta(x_t, u_t)}_{\text{Context encoder}} & 
\hat{r} &= 
\underbrace{d_\tau(b_h, b_c)}_{\text{Residual prediction}} & 
\\
\end{aligned}
\label{eq:full_dynamics}
\end{equation}


where $g$ are domain-specific ODEs which are the input `domain knowledge' about the system being modeled. 
Note that if the prior knowledge $g$ is identically zero, this method reduces to the Neural ODE predictions we discussed at the beginning of this section. We also study an ablated model, NDS0, which lacks the residual component $\hat{r}$ and context encoder network $d_\tau$. We note here that the context encoder is intended to potentially correct for model misspecification and noise but in the noiseless case with a model which is perfect, it may not be necessary. We explore this throughout Section \ref{experiments}.



\emph{Example 1: Lorenz system.}
To illustrate the full construction, we operate on the example of the the Lorenz system: a chaotic dynamical system originally defined to model atmospheric processes \cite{lorenz_deterministic_1963}. The system has 3-dimensional state (which we'll denote by $x, y, z$), 3 parameters, $\rho$, $\sigma$, and $\beta$, and no control input. The system is given by
\begin{equation}
\label{eq:lorenz}
\begin{aligned}
\dot{x} &= \sigma (y - x) & 
\dot{y} &= x (\rho - z) - y & 
\dot{z} &= xy - \beta z.
\end{aligned}
\end{equation}
For a given instantiation of the Lorenz system, we have values of $\phi = [\beta, \sigma, \rho]$ that are constant across the trajectory. So, we can instantiate  $h_\theta$ which outputs $\hat{\phi} = [\hat{\beta}, \hat{\sigma}, \hat{\rho}]$. We use the DOPRI5 method \cite{dormand_family_1980} to integrate the full neural dynamical system in Equation \ref{eq:full_dynamics}, with $g$ given by the system in Equation \ref{eq:lorenz} using the adjoint method of \cite{chen_neural_2018}. We use the state $x_{T'}$ as the initial condition for this integration. This gives a sequence $\{\hat{x}_t\}_{t=T'}^T$, which we evaluate and supervise with a loss of the form
\begin{equation}
\label{eq:loss}
\mathcal{L}_{\theta, \vartheta, \tau}(\{\hat{x}_{t_i}\}_{i=T' + 1}^T, \{x_{t_i}\}_{t = T' + 1}^T) = \sum_{t=T' + 1}^T||x_{t_i} - \hat{x}_{t_i}||_2^2.
\end{equation}
Because of the way we generate our input data, this is equivalent to Equation \ref{eq:prediction_problem}. We assume in our setting with full dynamics that the true dynamics lie in the function class established in Equation \ref{eq:full_dynamics}. By the method in \cite{chen_neural_2018} we can backpropagate gradients through this loss into the parameters of our NDS. Then algorithms in the SGD family will converge to a local minimum of our loss function. 

\paragraph{NDS with Partial System Dynamics}
Suppose we only had prior knowledge about some of the components of our system and none about others.
We can easily accomodate this incomplete information by simply `zeroing out' the function.
This looks like
\begin{equation}
\label{eq:partial_dynamics}
g_{\phi}(x, u, t) =\begin{bmatrix} \begin{cases}
g_{\phi}^{(i)}(x, u, t) & \text{if }g_{\phi}^{(i)}\text{ is known,}\\
0&\text{else.}
\end{cases}
\end{bmatrix}
\end{equation}
substituted into equation \ref{eq:full_dynamics}. In this setup, the residual term essentially makes the unknown dimensions unstructured Neural ODEs, which still can model time series well \cite{portwood_turbulence_2019}.




\paragraph{NDS with Approximate System Dynamics}
For Neural Dynamical Systems to be useful, they must handle situations where the known model is approximate. This is transparently handled by our formulation of Neural Dynamical Systems: the parameters of the approximate model $\hat{\phi}$ are predicted by $h_\theta(x_{1:T'}, u_{1:T'})$ and the residuals $\hat{r}$ are predicted by $d_\tau (b_h, b_c)$. This is the same as in the case where we have the correct dynamics, but we remove the assumption of a perfect model.

\emph{Example 2: Nuclear Fusion System.}
In this paper, we apply this technique to plasma dynamics in a tokamak. In a tokamak, two quantities of interest are the stored energy of the plasma, which we denote $E$ and its rotational frequency, $\omega$. The neutral beams and microwave heating allow us to add power ($P$) and torque ($T$) to the plasma. The plasma also dissipates energy and rotational momentum via transport across the boundary of the plasma, radiative cooling, and other mechanisms. While the detailed evolution of these quantities is described by turbulent transport equations, for the purposes of control and design studies, physicists often use reduced, volume-averaged models. The simple linear model (up to variable parameters) used for control development in \cite{boyer_feedback_2019} is used in this work.
\begin{equation}
\label{eq:tokamak_model}
\begin{aligned}
\dot{E} &= P - \frac{E}{\tau_e} & 
\dot{\omega} &= \frac{T}{n_im_iR_0} - \frac{\omega}{\tau_m}
\end{aligned}
\end{equation}

Here, $n_i$ is ion density, $m_i$ is ion mass, and $R_0$ is the tokamak major radius. We use the constant known values for these. $\tau_e$ and $\tau_m$ are the confinement times of the plasma energy and angular momentum, which we treat as variable parameters (because they are!). These are predicted by the neural network in our model.
We again use the model in Equation \ref{eq:full_dynamics} to give us a neural dynamical system which can learn the real dynamics starting from this approximation in Section \ref{fusion_experiments}.










\section{Problem Setting}
\vspace{-2mm}
\label{problem_setting}


Typically, a dynamical system $\dot{x} = f_\phi(x, u, t)$ with some parameters $\phi$ is the conventional model for system identification problems.
Here, state is $x\in \mathcal{X}$, control is $u\in \mathcal{U}$, and time is $t \in \mathbb{R}$. The objective is to predict future states given past states, past and future controls, and prior knowledge of the form of $f$.
We denote $x(\phi, t,{\bf u}, x_0 ) = x_0 + \int_0^t f_\phi(x, u, t)dt$ as the state obtained by integrating our dynamical system around $f$ to time $t$.




We consider in this work a more general setting and address the problem of prediction and control over a \textit{class of dynamical systems}, which we define as the set $\{\dot{x} = f_\phi(x, u, t)\mid \phi \in \Phi\}$
, where $\Phi$ is the space of parameters for the dynamical system (e.g. spring constant or terminal velocity). We can generate a trajectory from a class by sampling a $\phi\sim P(\Phi)$ for some distribution $P$ and choosing initial conditions and controls. In real data, we can view nature as choosing, but not disclosing, $\phi$. For a particular example $j$, we sample $\phi \sim P(\Phi)$ and $x_0 \sim P(X_0)$ and are given controls $\bf u$ indexed as $u(t)$ and input data $\{x(\phi, t_i, {\bf u}, x_0)\}_{i=0}^{T}$ during training. At test time, we give a shorter, prefix time series $\{x(\phi, t_i, {\bf u}, x_0)\}_{i=0}^{T'}$ but assume access to future controls.
Then the prediction objective for a class of systems for $N$ examples for timesteps $\{t_i\}_{T' + 1}^T$ is
   

   
\begin{equation*}
    \label{eq:prediction_problem}
    \hat{x} = \argmin_{\hat{x}}\underset{\underset{\phi \sim P(\Phi)} {x_0 \sim P(X_0)}}{\mathbb{E}}\left[\sum_{i=T' + 1}^T\left|\left|x(\phi, t_i, {\bf u}, x_0) - \hat{x}_{t_i}\right|\right|_2^2\right].
\end{equation*}


This objective differs from the traditional one in that implicitly, identifying $\phi$ for each trajectory needs to be done from the problem data in order to be able to predict the data generated by $f_\phi$.




Similarly, the control problem i
\begin{equation*}
    \label{eq:control_problem}
    \begin{aligned}
      \mathbf{u} = \min_{{\mathbf u}} \mathbb{E}_{\phi\sim P(\Phi), x_0 \sim P(X_0)}
    \left[ \int_0^t c(u(t), x(t))dt\right],\\
    \text{ s.t. }
    x(t) = x_0 + \int_0^{t}f_\phi(x, u, t)\,dt
    \end{aligned}
\end{equation*}
for some cost functional $c$.
We will primarily explore the prediction problem in this setting, but as secondary considerations, we explore robustness to noise, the ability to handle irregularly spaced input data, and the ability to recover the parameters $\phi$ which generated the original trajectories. We will also consider the control problem in a simple setting.







\section{Related Work}
\vspace{-1.5mm}
\label{relatedwork}

There is a long tradition of forecasting physical dynamics with machine learning. \cite{cressie_statistics_2015} lays out ideas from classical statistics for predicting spatiotemporal data. In the deep learning world, recurrent neural networks and long short-term memory networks have been used for a long time to address sequential problems \cite{hochreiter_long_1997}. However, these models struggle with continuous-time data and do not allow for the introduction of physical priors. Recently, there has been work in learning provably stable dynamics models by constraining the neural network to have a stable and jointly learned Lyapunov function \cite{manek_learning_2020}. We see opportunities to use these techniques in follow up work.

\paragraph{Neural Ordinary Differential Equations}

As most numerical ODE solvers are algorithms involving differentiable operations, it has always been possible in principle to backpropagate through the steps of these solvers dating back to at least \cite{runge_ueber_1895}. 
However, since each step of the solver involves calling the derivative function, na\"ive backpropagation incurs an $O(n)$ memory cost in the number of ODE solution steps, where $n$ is the number of derivative calls made by the solver. 
\cite{chen_neural_2018} demonstrated that by computing gradients through the adjoint sensitivity method, the memory complexity of backpropagating through a family of ODE solvers can be reduced to $O(1)$ for a fixed network, as opposed to the naive $O(n)$.
However, this work only used generic neural networks as the derivative function and did not consider dynamics. They also provide a PyTorch package which we have built off of in our work.


There has been some work using neural ordinary differential equations to solve physical problems. \cite{portwood_turbulence_2019} used a fully-connected neural ODE with an RNN encoder and decoder to model Navier-Stokes problems. \cite{rudy_deep_2019} used a neural network integrated with a Runge-Kutta method for noise reduction and irregularly sampled data. There has also been work learning the structure of dynamical systems, first with a convolutional-deconvolutional warping scheme inspired by the solutions to advection-diffusion PDEs \cite{de_bezenac_deep_2018}, then with a Neural ODE which was forced to respect boundary conditions and a partial observation mechanism \cite{ayed_learning_2019}. There was also work on constraining a neural network model to respect separable conservative Hamiltonian dynamics, 
but like all of the aforementioned methods, this does not incorporate information about the dynamics of a particular system \cite{chen_symplectic_2019}.
In general, none of these methods incorporate prior knoweledge as explicitly as including appropriate equations in the statistical model. Furthermore, many of these methods focus on a specific problem, whereas we give a way to apply specific knowledge about a variety of problems. 

An interesting approach that builds on a previous version of this paper \cite{mehta_neural_2020} is developed in \cite{guen2021augmenting}, where it is shown that iteratively solving a Lagrangian formulation of the problem for a fixed set of parameters and varying Lagrange multiplier results in predictions which maximally use the physical model. Our problem setting is a generalization of theirs in that we allow the parameters to vary among rollouts and predict them on the fly.



\paragraph{Machine Learning for Nuclear Fusion:}
As far back as 1995, \cite{van_milligen_neural_1995} showed that by approximating the differential operator with a (single-layer, in their case) neural network, one could fit simple cases of the Grad-Shafranov equation for magnetohydrodynamic equilibria. Recently, work has shown that plasma dynamics are amenable to neural network prediction.
In particular, \cite{kates-harbeck_predicting_2019} used a convolutional and LSTM-based architecture to predict possible plasma disruptions (when a plasma instability grows large and causes a loss of plasma containment and pressure).

There has also been work in the field of plasma control: a neural network model of the neutral beam injection for the DIII-D tokamak has been deployed in order to diagnose the effect of controls on shots conducted at the reactor \cite{boyer_real-time_2019}. 
Additionally, \cite{boyer_feedback_2019} used classic control techniques and a simpler model of the dynamics to develop a controller that allows characteristics of the tokamak plasma to be held at desired levels.
