Concepts of dynamical systems theory had been frequently used to improve neural network performance [1,2,3] but Neural Ordinary Differential Equations (NODE) appear to have been formally introduced by Chen et al. [4]. NODE provide an explicit connection between deep feed-forward neural networks and dynamical systems and are considered to provide a bridge between modern deep learning and classical mathematical/numerical modelling. NODE provide a flexible trade-off between efficiency, memory costs and accuracy. The approximation capabilities [5,6] of NODE are particularly useful for time-series modelling [4,7,8], generative models for continuous normalizing flows [4,9] and modeling/controlling physical environments see, e.g., [10].

Neural ODEs are trained by minimizing a least-squares quadratic scalar-valued “loss function” by computing its gradients with respect to the weights to be optimized using a first-order optimizer such as “stochastic gradient descent” [11,12]). Since ODE solvers (e.g., Runge–Kutta solvers) perform differentiable algebraic operations, the gradients of the loss function can be calculated by the so-called “direct method” which directly backpropagates through the operations performed by the ODE solver. However, when the dynamics are complex, the “direct method” can lead to an arbitrarily large number of function evaluations for adaptive solvers while storing all of the intermediate activations during the “solving” process, so the “direct method” becomes prohibitively memory intensive. A NODE-training method which is less memory intensive is the so-called “adjoint method” [13,14,15], which solves an ODE (related to the original NODE) backwards in time. The direct method is faster but is more memory intensive than the adjoint method. The one-dimensional definite integrals, which appear when computing gradients via the “adjoint method” are traditionally evaluated by solving them as differential equations, which considerably slows down the training process. Evaluating these one-dimensional definite integrals by using Gauss–Legendre quadrature (rather than solving them as ODEs) has been shown [16] to be faster than ODE-based methods while retaining memory efficiency, thus speeding up the training of NODE.

The gradients of the loss function are often called “sensitivities” in the literature on neural nets and aspects of the optimization/training procedure are occasionally called “sensitivity analysis.” But the “loss function” is of interest only for the “training” phase of the NODE and the “sensitivities of the loss function” are driven towards the ideal zero-values by the minimization process while optimizing the NODE weights/parameters. Furthermore, after the NODE is optimized to reproduce the underlying physical system as closely as possible, the responses of interest for the NODE-modeled system is no longer a “loss function” but are various functions of the NODE’s “decoder”-output. Since the physical system being modeled by the NODE comprises itself parameters that stem from measurements or computations, they are not perfectly well-known, but are afflicted by uncertainties that stem from the respective experiments and/or computations. Hence, it is important to quantify the uncertainties induced in the NODE decoder-output by the uncertainties that afflict the parameters/weights underlying the physical system modeled by the NODE. The quantification of the uncertainties in the NODE-decoder and derived results (i.e., “NODE-responses”) of interest require the computation of the sensitivities of the NODE-decoder with respect to the optimized NODE-weights/parameters. However, a “NODE sensitivity analysis” methodology for computing efficiently exact expressions of decoder-sensitivities with respect to the optimized parameters/weights, including with respect to the initial conditions/encoder, does not seem to be available in the literature.

The scope of this work is to present a novel methodology for computing all of the first-order sensitivities, exactly and exhaustively, of the responses of the post-training optimized NODE-decoder with respect to the optimized/trained weights involved in the NODE’s decoder, hidden layers, and encoder. The general mathematical representation of the NODE-network considered in this work is presented in Section 2. As a specific illustrative paradigm application, Section 3 presents the NODE conceptual representation of the Nordheim-Fuchs phenomenological reactor dynamics/safety model [17,18]. This paradigm illustrative model has been chosen because it is representative of typical NODE-applications while admitting closed-form analytical solutions for the quantities of interest, including the functions describing the hidden layers, encoder, decoder, and sensitivities of decoder responses. Section 4 presents the Mathematical Framework of the novel “First-Order Comprehensive Adjoint Sensitivity Analysis Methodology for Neural Ordinary Differential Equations (1^st-CASAM-NODE).” Section 5 illustrates the application of the 1^st-CASAM-NODE methodology to compute all of the first-order sensitivities of Nordheim-Fuchs model responses with respect to the underlying parameters. Specifically, Subsections 5.1 through 5.4, respectively, illustrate the application of the 1^st-CASAM-NODE methodology for computing the first-order sensitivities with respect to the underlying model parameters and initial conditions of the following responses: (i) the reactor’s flux; (ii) the reactor’s energy release; (iii) the reactor’s temperature; and (iv) the reactor’s thermal conductivity. Using the “energy-released” response as a paradigm, Subsection 5.5 illustrates an alternative path for computing first-order sensitivities by applying the “First-Order Feature Adjoint Sensitivity Analysis Methodology for Nonlinear Systems (1^st-FASAM-N)” [19], which is the most efficient procedure for computing first-order sensitivities, but which may require the construction of a dedicated neural net for this purpose.

A general mathematical representation of a NODE-network is provided by the following system of so-called “augmented” equations:

where:

The weights of the NODE are adjusted/calibrated by “training” the NODE, using gradients of a scalar loss functional, denoted as $L\left[h\left(t\right);\theta ;t\right]$ which is designed to represent the deviations/discrepancies between the responses/outputs of the NODE and the “true” values obtained from measurements (or by other means, independently of the NODE). There are several methods for accomplishing this “training,” all of which require that the functions underlying the NODE, i.e., $h\left(t\right)$, $f\left[h\left(t\right);\theta ;t\right]$, $h_e\left(x,w\right)$ and $h_d\left[h\left(t\right);\phi \right]$ be differentiable with respect to their arguments. For complex systems, involving many parameters, the so-called “adjoint method” [13,14,15] offers an optimal compromise between memory requirements and computational intensity. This method computes the required gradients of the loss function by evaluating the following integral: where the so-called “adjoint function” $a\left(t\right)\triangleq {\left[a_1\left(t\right),\mathrm{...},a_{TH}\left(t\right)\right]}^{†}$ satisfies the following “adjoint equation” computed backwards in time:

After the “training” of the NODE has been accomplished, the various “weights” will have been assigned “optimal” values which will have minimized the chosen loss functional $L\left[h\left(t\right);\theta ;t\right]$. These “optimal” values will be denoted using a superscript “zero” as follows: ${\theta }^0\triangleq {\left[{\theta }_1^0,\mathrm{...},{\theta }_{TW}^0\right]}^{†}$ and $w^0\triangleq {\left[w_1^0,\mathrm{...},w_{TEW}^0\right]}^{†}$. These optimal values are used to compute the optimal values for the system responses, which will be denoted as $r_n^0\left[h\left(t_f\right);\phi \right]$. However, since the physical parameters and the initial conditions underlying the actual physical system (which is represented by the optimized NODE) are not known exactly (because they are actually subject to uncertainties), it follows that the optimal values obtained for the weights are actually just nominal values that are used to compute the nominal/optimal response-values $r_n^0\left[h\left(t_f\right);\phi \right]$. The uncertainties in the various weights and initial conditions will induce uncertainties in the system responses, which can be computed deterministically by using the well-known “propagation of errors” methodology, originally proposed by Tuckey [20] and subsequently extended to sixth-order by Cacuci [21].

The Nordheim-Fuchs phenomenological model [17,18] describes a short-time self-limiting power transient in a nuclear reactor system having a negative temperature coefficient in which a large amount of reactivity is suddenly inserted, either intentionally or by accident. The response of such a reactor system can be estimated by considering that the reactivity insertion is sufficiently large and the time-span of the transient phenomena under consideration is of the order of the life-time of prompt-neutrons. For such short times, the effects of delayed neutrons and the local spatial variations of the neutron distribution in the reactor can be neglected, and the heat generated during the transient remains within the reactor. Using the notation of Lamarsh [17], the Nordheim-Fuchs paradigm model describing such a self-limiting power transient comprises the following balance equations:

where $l_p$ denotes the prompt-neutron lifetime, $k\left(t\right)$ denotes the reactor’s multiplication factor, and ${\phi }_0$ denotes the initial (i.e., extant flux) prior to initiating the transient at time $t=0$.

where $\gamma$ denotes the recoverable energy per fission; ${\Sigma }_f\triangleq {\sigma }_f{N}_f$ denotes the reactor’s effective macroscopic fission cross section, where ${\sigma }_f$ denotes the reactor’s equivalent microscopic fission cross section while $N_f$ denotes the reactor’s equivalent atomic number density.

where $E\left(t\right)$ denotes the total energy released (per cm³) at time $t$ in the reactor since the onset of reactivity change; $c_p$ denotes the specific heat (per cm³) of the reactor.

Equations (8)‒(12) can be transformed into the following system of nonlinear differential equations:

The Nordheim-Fuchs model described by Eqs. (13)‒(15) can be solved analytically to obtain exact closed-form expression for the state functions $\phi \left(t\right)$, $E\left(t\right)$, and $T\left(t\right)$, as follows:

The closed-form solution of Eq. (17) has the following form: where:

The typical results of interest (called “model response”) for the Nordheim-Fuchs model are as follows:

where $\delta \left(t-\tau \right)$ denotes the Dirac-delta functional.

Comparing the structure of the Nordheim-Fuchs model, cf. Eqs. (13)‒(15), to the generic structure of a NODE, cf. Eqs. (1) and (2), indicates the following correspondences:

The actual values of the components of the vectors $\theta$ and $x$ are unknown even after having trained the NODE, since the actual values of the parameters underlying the Nordheim-Fuchs model are experimentally-measured and are thus subject to uncertainties. However, the nominal values of these parameters are considered to be known, and are considered to be exactly reproducible by the “trained” NODE; these nominal values will be denoted using a superscript “zero,” as follows:

Consequently, the exact values of the functions $h\left(t\right)\triangleq {\left[h_1\left(t\right),h_2\left(t\right),h_3\left(t\right)\right]}^{†}\triangleq {\left[\phi \left(t\right),E\left(t\right),T\left(t\right)\right]}^{†}$ are unknown but their nominal values $h^0\left(t\right)\triangleq {\left[h_1^0\left(t\right),h_2^0\left(t\right),h_3^0\left(t\right)\right]}^{†}\triangleq {\left[{\phi }^0\left(t\right),E^0\left(t\right),T^0\left(t\right)\right]}^{†}$ are known after having solved Eqs. (13)‒(15) at the nominal values $\left({\theta }^0,\hspace{0.17em}\hspace{0.17em}{x}^0\right)$.

The NODE-representations, cf. Eq. (4), of the responses considered in Eqs. (23)‒(24)have the following expressions, respectively:

To illustrate the efficient computation of responses involving decoders having their own parameters/weights, the thermal conductivity of the conceptual material of the Nordheim-Fuchs reactor model will be considered to be a “decoder” response having the following expression:

At the optimal/nominal parameter values, the optimal/nominal solution $h^0\left(t\right)$ will satisfy the following forms of Eqs. (1) and (2):

, Furthermore, the vector of optimal/nominal response will have components that are obtained by using the nominal values for the respective functions and parameters, i.e.:

The known nominal values $x^0$ of the initial conditions will differ from the true but unknown values $x$ of the initial conditions by variations denoted as $\delta x\triangleq x-x^0$. Furthermore, the known nominal values $w^0$ of the weights characterizing the encoder will differ from the true but unknown values $w$ of the respective weights by variations denoted as $\delta w\triangleq w-w^0$. Similarly, the nominal values ${\theta }^0$ and ${\phi }^0$, respectively, will differ by variations $\delta \theta \triangleq \theta -{\theta }^0$ and $\delta \phi \triangleq \phi -{\phi }^0$, respectively, from the corresponding true but unknown values $\theta$ and $\phi$. Since the forward state functions $h\left(t\right)$ are related to the weights and initial conditions through Eqs. (1) and (2), it follows that the variations in these weights and initial conditions will induce corresponding variations $v^{\left(1\right)}\left(t\right)\triangleq {\left[\delta h_1\left(t\right),\dots ,\delta h_{TH}\left(t\right)\right]}^{†}$ around the nominal solution $h^0\left(t\right)$. In turn, the variations $\delta \phi$ and $v^{\left(1\right)}\left(t\right)$ will induce variations $\delta r_n\left(h^0;{\phi }^0;v^{\left(1\right)};\delta \phi \right)$ in the system’s response.

The 1^st-CASAM-NODE methodology for computing the first-order sensitivities of the response with respect to the model’s weights and initial conditions will be established by following the same principles as those underlying the 1^st-CASAM-N methodology [22], which commence by noting that Cacuci [23] has shown that the most general definition of the sensitivity of an operator-valued model response $R\left(e\right)$ with respect to variations $\delta e$ in the model parameters and state functions in a neighborhood around the nominal functions and parameter values $e^0$, is given by the 1^st-order Gateaux- (G-) variation, which will be denoted as $\delta R\left(e^0;\delta e\right)$ and is defined as follows: for a scalar $\epsilon$ and for all (i.e., arbitrary) vectors $\delta e$ in a neighborhood $\left(e^0+\epsilon \delta e\right)$ around $e^0$. The G-variation $\delta R\left(e^0;\delta e\right)$ is an operator defined on the same domain as $R\left(e\right)$ and has the same range as $R\left(e\right)$. The G-variation $\delta R\left(e^0;\delta e\right)$ satisfies the relation: $R\left(e^0+\epsilon \delta e\right)-R\left(e^0\right)=\delta R\left(e^0;\delta e\right)+\Delta \left(\delta e\right)\hspace{0.33em},$ with $\underset{\epsilon \to 0}{\lim }\left[\Delta \left(\epsilon \delta e\right)\right]/\epsilon =0$. When the G-variation $\delta R\left(e^0;\delta e\right)$ is linear in the variation $\delta e$, it can be written in the form $\delta R\left(e^0;\delta e\right)={\left\{\partial R/\partial e\right\}}_{e^0}\delta e$, where ${\left\{\partial R/\partial e\right\}}_{e^0}$ denotes the first-order G-derivative of $R\left(e\right)$ with respect to $e$ evaluated at $e^0$.

Applying the definition provided in Eq. (38) to Eq. (4) yields the following expression for the first-order G-variation $\delta r_n\left(h^0;{\phi }^0;v^{\left(1\right)};\delta \phi \right)$ of the response $r_n\left(h;\phi \right)$: where $v^{\left(1\right)}\triangleq {\left[v_1^{\left(1\right)}\left(t\right),\mathrm{...},v_{TH}^{\left(1\right)}\left(t\right)\right]}^{†}$ and:

This, the quantity ${\left\{\partial h_n^d\left[h\left(t\right);\phi \right]/\partial \phi \right\}}_{\left(h^0;{\phi }^0\right)}$ in Eq. (40) denotes the partial G-derivatives of the response $h_n^d\left[h\left(t\right);\phi \right]$ with respect to the decoder weights $\phi \triangleq {\left[{\phi }_1,\mathrm{...},{\phi }_{TD}\right]}^{†}$, evaluated at the nominal values $\left(h^0;{\phi }^0\right)$. The quantity ${\left\{\delta r_n\left(h^0;{\phi }^0;\delta \phi \right)\right\}}_{dir}$ is called the “direct-effect term” because it arises directly from parameter variations $\delta \phi$ and can be computed directly using the nominal values $\left(h^0;{\phi }^0\right)$. The quantity ${\left\{\delta r_n\left(h^0;{\phi }^0;\delta h;\delta \phi \right)\right\}}_{ind}$ is called the “indirect-effect term” because it arises indirectly, through the variations $v^{\left(1\right)}\left(t\right)$ in the hidden state functions $h\left(t\right)$. The indirect-effect term can be quantified only after having determined the variations $v^{\left(1\right)}\left(t\right)$, which are caused by the variations $\delta x$, $\delta w$ and $\delta \theta$.

The first-order relationships between the variations $v^{\left(1\right)}\left(t\right)$, $\delta x$, $\delta w$ and $\delta \theta$ are obtained by computing the first-order G-variation of Eqs. (1) and (2), which are obtained, by definition, as follows:

Carrying out the operations indicated in Eqs. (42) and (43) yields the following system of equations: where:

The system comprising Eqs. (44) and (45) is called the “1^st-Level Variational Sensitivity System” (1^st-LVSS), and its solution, $v^{\left(1\right)}\left(t\right)$, is called the “1^st-level variational sensitivity function.” Note that the 1^st-LVSS would need to be solved anew for each component of the variations $\delta x$, $\delta w$ and $\delta \theta$, which would be prohibitively expensive computationally.

The need for solving the 1^st-LVSS can be avoided if the indirect-effect term defined in Eq. (41) could be expressed in terms of a “right-hand side” that does not involve the function $v^{\left(1\right)}\left(t\right)$. This goal can be achieved by expressing the right-side of Eq. (41) in terms of the solutions of the “1^st-Level Adjoint Sensitivity System (1^st-LASS),” the construction of which requires the introduction of adjoint operators. Adjoint operators can be defined in Banach spaces but are most useful in Hilbert spaces. For the NODE considered in this work, the appropriate Hilbert space is defined on the domain ${\Omega }_t\triangleq \left[t_0,t_f\right]$ and will be denoted as $H_1\left({\Omega }_t\right)$, so that $v^{\left(1\right)}\left(t\right)\in H_1\left({\Omega }_t\right)$. In $H_1\left({\Omega }_t\right)$, the inner product of two vectors in $u^{(a)}\left(t\right)\in H_1\left({\Omega }_t\right)$ and $u^{(b)}\left(t\right)\in H_1\left({\Omega }_t\right)$ will be denoted as ${〈u^{(a)},u^{(b)}〉}_1$, and is defined as follows: where the “dot” indicates the “scalar product of two vectors” defined as follows: $u^{(a)}\left(t\right)·u^{(b)}\left(t\right)\hspace{0.17em}\hspace{0.17em}\triangleq {\left[u^{(a)}\left(t\right)\right]}^{†}{u}^{(b)}\left(t\right)\hspace{0.17em}\hspace{0.17em}\triangleq {\displaystyle \sum _{i=1}^{TH}{u}_i^{(a)}\left(t\right)u_i^{(b)}}\left(t\right)={\left[u^{(b)}\left(t\right)\right]}^{†}{u}^{(a)}\left(t\right)$.

The next step is to form the inner product of Eq. (44) with a vector $a^{(1)}\left(t\right)\triangleq \left[a_1^{\left(1\right)}\left(t\right),\dots ,a_{TD}^{\left(1\right)}\left(t\right)\right]\in H_1\left({\Omega }_t\right)$, where the superscript “(1)” indicates “1^st-Level”, to obtain the following relationship:

Using the definition of the adjoint operator in $H_1\left({\Omega }_t\right)$, the left-side of Eq. (51) is transformed as follows, after integrating by parts over the independent variable $t$:

The last term on the right-side of Eq. (52) is now required to represent the “indirect-effect” term defined in Eq. (41), which is achieved by requiring that the 1^st-level adjoint function $a^{(1)}\left(t\right)$ satisfy the following relation:

The definition of the 1^st-level adjoint sensitivity function $a^{(1)}\left(t\right)$ is now completed by requiring it to satisfy (adjoint) “boundary conditions at the final time” $t=t_f$ so as to eliminate the term containing the unknown values $v^{\left(1\right)}\left(t_f\right)$ in Eq. (52). This aim is achieved by requiring that

The system of equations comprising Eqs. (53) and (54) constitute the “1^st-Level Adjoint Sensitivity System (1^st-LASS)” for the 1^st-level adjoint function $a^{(1)}\left(t\right)$. Evidently, the 1^st-LASS is independent of parameter variations and needs to be solved just once to obtain the 1^st-level adjoint function $a^{(1)}\left(t\right)$. Notably, the 1^st-LASS has the same form as the “adjoint equations” used for training the NODE, cf. Eqs. (6) and (7), but with the “response” $\partial h_n^d\left[h\left(t\right);\phi \right]/\partial h\left(t\right)\delta \left(t-t_f\right)$ being the “source” for the 1^st-LASS, whereas the “source” in the “training” of the NODE was the “loss functional” $L\left[h\left(t\right);\theta ;t\right]/\partial h\delta \left(t-t_f\right)$. Evidently, the 1^st-level adjoint sensitivity function $a^{(1)}\left(t\right)$ is the counterpart of the “adjoint function” $a\left(t\right)$ in the “training” of the NODE.

Using the results represented by Eqs. (53), (54), (51), and (41) in Eq. (52) yields the following alternative expression for the “indirect-effect” term, which does not involve the 1^st-level variational sensitivity function $v^{\left(1\right)}\left(t\right)$ but involves the 1^st-level adjoint function $a^{(1)}\left(t\right)$:

Using in Eq. (55) the expression provided for $v^{\left(1\right)}\left(t_0\right)$ in Eq. (45) yields the following expression for the “indirect-effect” term:

Replacing the expression obtained in Eq. (55) for the “indirect-effect term” together with the expression of the direct-effect term provided by Eq. (40) into Eq. (39) yields the following expression for the first-order G-variation $\delta r_n\left(h^0;{\phi }^0;v^{\left(1\right)};\delta \phi \right)$ of the response $r_n\left(h;\phi \right)$:

As indicated by the right-side of Eq. (57), the (partial) sensitivities of the response $r_n\left(h;\phi \right)$ are provided by the following expressions, all of which are to be evaluated at the nominal values of all functions and parameters/weights:

The application of the 1^st-CASAM-NODE methodology to compute the first-order sensitivities of the responses $r_1\left(h\right)$, $r_2\left(h\right)$, $r_3\left(h\right)$ and $r_4\left(h\right)$ with respect to the Nordheim-Fuchs model’s parameters and initial conditions will be presented below in Subsections 5.1 through 5.4, respectively. Using the “energy-released” response $r_2\left(h\right)=E\left(t_f\right)$ as a paradigm, Subsection 5.5 will illustrate an alternative path for computing first-order sensitivities by applying the “First-Order Feature Adjoint Sensitivity Analysis Methodology for Nonlinear Systems (1^st-FASAM-N)” [24], which is the most efficient procedure for computing sensitivities, but which may require the construction of a dedicated neural net for this purpose.