First-Order Comprehensive Adjoint Sensitivity Analysis Methodology for Neural Ordinary Differential Equations: Mathematical Framework and Illustrative Application to the Nordheim-Fuchs Reactor Safety Model

1. Introduction

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.

2. Neural Ordinary Differential Equations (NODE): Basic Properties and Uses

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

$${\displaystyle \frac{dh\left(t\right)}{dt}}=f\left[h\left(t\right);\theta ;t\right],\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}t>0$$

(1)

$$h\left(t_0\right)=h_e\left(x,w\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}at\hspace{0.17em}\hspace{0.17em}t=t_0$$

(2)

$$r\left(t_f\right)=h_d\left[h\left(t_f\right);\phi \right],\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}at\hspace{0.17em}\hspace{0.17em}t=t_f$$

(3)

where:

(i)

The quantity $t$ is a time-like independent variable which parameterizes the dynamics of the hidden/latent neuron units; the initial value is denoted as $t_0$ (which can be considered to be an initial measurement time) while the stopping value is denoted as $t_f$ (which can be considered to be the next measurement time).

(ii)

The $TH$-dimensional vector-valued function $h\left(t\right)\triangleq {\left[h_1\left(t\right),\mathrm{...},h_{TH}\left(t\right)\right]}^{†}$ represents the hidden/latent neural networks. In this work, all vectors are considered to be column vectors and the dagger “$†$” symbol will be used to denote “transposition.” The symbol “$\triangleq$” signifies “is defined as” or, equivalently, “is by definition equal to.”

(iii)

The$TH$-dimensional vector-valued nonlinear function $f\left[h\left(t\right);\theta ;t\right]\triangleq {\left[f_1\left(h;\theta ;t\right),\mathrm{...},f_{TH}\left(h;\theta ;t\right)\right]}^{†}$ models the dynamics of the latent neurons with learnable scalar adjustable weights represented by the components of the vector $\theta \triangleq {\left[{\theta }_1,\mathrm{...},{\theta }_{TW}\right]}^{†}$, where $TW$ denotes the total number of adjustable weights in all of the latent neural nets.

(iv)

The $TH$-dimensional vector-valued function $h_e\left(x,w\right)\triangleq {\left\{h_1^e\left(x,w\right),\mathrm{...},h_{TH}^e\left(x,w\right)\right\}}^{†}$ represents the “encoder” which is characterized by “inputs” $x\triangleq {\left[x_1,\mathrm{...},x_{TI}\right]}^{†}$and “learnable” scalar adjustable weights $w\triangleq {\left[w_1,\mathrm{...},w_{TEW}\right]}^{†}$, where $TI$ denotes the total number of “inputs” and $TEW$ denotes the total number of “learnable encoder weights” that define the “encoder.”

(v)

The $TR$-dimensional vector-valued function $r\left(t_f\right)\triangleq {\left\{r_1\left[h\left(t_f\right);\phi \right],\mathrm{...},r_{TR}\left[h\left(t_f\right);\phi \right]\right\}}^{†}=h_d\left[h\left(t_f\right);\phi \right]$ represents the vector of “system responses.” The vector-valued function $h_d\left[h\left(t_f\right);\phi \right]\triangleq {\left\{h_1^d\left[h\left(t_f\right);\phi \right],\mathrm{...},h_{TR}^d\left[h\left(t_f\right);\phi \right]\right\}}^{†}$ represents the “decoder” with learnable scalar adjustable weights, which are represented by the components of the vector $\phi \triangleq {\left[{\phi }_1,\mathrm{...},{\phi }_{TD}\right]}^{†}$, where $TD$ denotes the total number of adjustable weights that characterize the “decoder.” Each component $r_n\left[h\left(t_f\right);\phi \right]$ can be epresented in integral form as follows:

$$r_n\left(h;\phi \right)={\displaystyle {\int }_{t_0}^{t_f}{h}_n^d\left[h\left(t\right);\phi \right]}\hspace{0.17em}\delta \left(t-t_f\right)dt;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}n=1,\mathrm{...},TR.$$

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:

$${\displaystyle \frac{\partial L}{\partial \theta }}={\displaystyle {\int }_{t_0}^{t_f}{\left[a\left(t\right)\right]}^{†}\frac{\partial f\left[h\left(t\right);\theta ;t\right]}{\partial \theta }}\hspace{0.17em}dt$$

(5)

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:

$${\displaystyle \frac{da\left(t\right)}{dt}}=-{\left[a\left(t\right)\right]}^{†}{\displaystyle \frac{\partial f\left[h\left(t\right);\theta ;t\right]}{\partial h}},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}t>0$$

(6)

$$a\left(t\right)={\displaystyle \frac{\partial L\left[h\left(t\right);\theta ;t\right]}{\partial h}},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{at}\hspace{0.17em}\hspace{0.17em}t=t_f$$

(7)

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].

3. Illustrative Paradigm Application: NODE Conceptual Modeling of the Nordheim-Fuchs Phenomenological Reactor Dynamics/Safety Model

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:

• The time-dependent neutron balance (point kinetics) equation for the neutron flux $\phi \left(t\right)$:

$${\displaystyle \frac{d\phi \left(t\right)}{dt}}={\displaystyle \frac{k\left(t\right)-1}{l_p}}\phi \left(t\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}t>0$$

(8)

$$\phi \left(0\right)={\phi }_0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}t=0$$

(9)

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$.

2.

The energy production equation:

$$\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}E\left(t\right)=\gamma {\Sigma }_f{\displaystyle \underset{0}{\overset{t}{\int }}\phi \left(x\right)\hspace{0.17em}}dx$$

(10)

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.

3.

The energy conservation equation:

$$c_p\left[T\left(t\right)-T_0\right]=E\left(t\right)$$

(11)

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.

4.

The reactivity-temperature feedback equation: $k\left(t\right)=k_0-{\alpha }_T{k}_0\left[T\left(t\right)-T_0\right]$, where $k_0\triangleq k\left(0\right)\ge 1$ denotes the changed multiplication factor following the reactivity insertion at $t=0$, ${\alpha }_T$ denotes the magnitude of the negative temperature coefficient, $T\left(t\right)$ denotes the reactor’s temperature, and $T_0$ denotes the reactor’s initial temperature at time $t=0$. For illustrating the application of the 1st-FASAM methodology, it suffices to consider the special case of a “prompt critical transient”, when the reactor becomes prompt critical after the reactivity insertion, i.e., when $k_0=1$, so that the reactivity-temperature feedback equation takes on the following particular form:

$$k\left(t\right)=1-{\alpha }_T\left[T\left(t\right)-T_0\right]$$

(12)

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

$${\displaystyle \frac{d\phi \left(t\right)}{dt}}=-{\displaystyle \frac{{\alpha }_T}{l_p{c}_p}}E\left(t\right)\phi \left(t\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}t>0. \phi \left(0\right)={\phi }_0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}t=0$$

(13)

$${\displaystyle \frac{dE\left(t\right)}{dt}}=\gamma {\sigma }_f{N}_f\phi \left(t\right),\hspace{0.17em}\hspace{0.17em}E\left(0\right)=0$$

(14)

$${\displaystyle \frac{dT\left(t\right)}{dt}}={\displaystyle \frac{\gamma {\sigma }_f{N}_f}{c_p}}\phi \left(t\right);\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}T\left(0\right)=T_0$$

(15)

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:

Eliminating the function $\phi \left(t\right)$ from Eqs. (13) and (14) yields a nonlinear differential equation which can be integrated directly to obtain the following relation:

$$\phi \left(t\right)=-{\displaystyle \frac{{\alpha }_T}{2l_p{c}_p\gamma {\sigma }_f{N}_f}}{E}^2\left(t\right)+{\phi }_0$$

(16)

(ii)

Using Eq. (16) in Eq. (14) yields the following nonlinear equation for the released energy $E\left(t\right)$:

$${\displaystyle \frac{dE\left(t\right)}{dt}}=-{\displaystyle \frac{{\alpha }_T}{2l_p{c}_p}}{E}^2\left(t\right)+{\phi }_0\gamma {\sigma }_f{N}_f,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}E\left(0\right)=0$$

(17)

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

$$E\left(t\right)=K_1\left(\alpha \right)\tanh \left[t{K}_2\left(\alpha \right)\right]$$

(18)

where:

$$K_1\left(\alpha \right)\triangleq {\left[{\displaystyle \frac{2{\phi }_0\gamma {\sigma }_f{N}_f{l}_p{c}_p}{{\alpha }_T}}\right]}^{1/2};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{K}_2\left(\alpha \right)\triangleq {\left[{\displaystyle \frac{{\alpha }_T{\phi }_0\gamma {\sigma }_f{N}_f}{2l_p{c}_p}}\right]}^{1/2}.$$

(19)

(iii)

Replacing Eq. (18) into Eq. (16) yields the following closed-form expression for $\phi \left(t\right)$:

$$\phi \left(t\right)={\phi }_0\left\{1-{\tanh }^2\left[t{K}_2\left(\alpha \right)\right]\right\}={\displaystyle \frac{{\phi }_0}{{\cosh }^2\left[t{K}_2\left(\alpha \right)\right]}}$$

(20)

(iv)

Replacing Eq. (18) into Eq. (11) yields the following closed-form expression for $T\left(t\right)$:

$$T\left(t\right)=T_0+{\displaystyle \frac{K_1\left(\alpha \right)}{c_p}}\tanh \left[t{K}_2\left(\alpha \right)\right]$$

(21)

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

(i)

The neutron flux $\phi \left(\tau \right)$ in the reactor at a “final time” instance denoted as $t=\tau$, after the initiation at $t=0$ of the prompt-critical power transient, which can be defined mathematically as follows:

$$\phi \left(\tau \right)={\displaystyle {\int }_0^{\tau }\phi \left(t\right)\delta \left(t-\tau \right)dt}$$

(22)

(ii)

The total energy per cm³, $E\left(\tau \right)$, released at a user-chosen “final time” instance denoted as $t=\tau$, after the initiation at $t=0$ of the prompt-critical power transient, which can be defined mathematically as follows:

$$E\left(\tau \right)={\displaystyle {\int }_0^{\tau }E\left(t\right)\delta \left(t-\tau \right)dt}$$

(23)

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

(iii)

The reactor’s temperature $T\left(\tau \right)$ at a “final time” instance denoted as $t=\tau$ after the initiation at $t=0$ of the prompt-critical power transient, which can be defined mathematically as follows:

$$T\left(\tau \right)={\displaystyle {\int }_0^{\tau }T\left(t\right)\delta \left(t-\tau \right)dt}$$

(24)

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:

$$h\left(t\right)\triangleq {\left[h_1\left(t\right),\mathrm{...},h_{TH}\left(t\right)\right]}^{†}\triangleq {\left[\phi \left(t\right),E\left(t\right),T\left(t\right)\right]}^{†};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}TH=3;$$

(25)

$$\theta \triangleq {\left[{\theta }_1,\mathrm{...},{\theta }_{TW}\right]}^{†}\triangleq {\left({\alpha }_T,l_p,c_p,\gamma ,{\sigma }_f,N_f\right)}^{†};\hspace{0.17em}\hspace{0.17em}x\triangleq {\left[x_1,x_2\right]}^{†}\triangleq {\left({\phi }_0,T_0\right)}^{†};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}TW=6,\hspace{0.17em}\hspace{0.17em}TI=2$$

(26)

$$f_1\left(h;\theta ;t\right)\triangleq -{\displaystyle \frac{{\alpha }_T}{l_p{c}_p}}E\left(t\right)\phi \left(t\right)\hspace{0.17em}\triangleq -{\displaystyle \frac{{\theta }_1}{{\theta }_2{\theta }_3}}{h}_1\left(t\right)h_2\left(t\right)$$

(27)

$$f_2\left(h;\theta ;t\right)\triangleq \gamma {\sigma }_f{N}_f\phi \left(t\right)\hspace{0.17em}\triangleq {\theta }_4{\theta }_5{\theta }_5{h}_1\left(t\right)$$

(28)

$$f_3\left(h;\theta ;t\right)\triangleq {\displaystyle \frac{\gamma {\sigma }_f{N}_f}{c_p}}\phi \left(t\right)\triangleq {\displaystyle \frac{{\theta }_4{\theta }_5{\theta }_6}{{\theta }_3}}{h}_1\left(t\right)$$

(29)

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:

$${\theta }^0\triangleq {\left[{\theta }_1^0,\mathrm{...},{\theta }_6^0\right]}^{†}\triangleq {\left({\alpha }_T^0,l_p^0,c_p^0,{\gamma }^0,{\sigma }_f^0,N_f^0\right)}^{†}\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{x}^0\triangleq {\left[x_1^0,x_2^0,x_3^0\right]}^{†}\triangleq {\left[{\phi }_0^0,0;T_0^0\right]}^{†}$$

(30)

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:

$$r_1\left(h\right)={\displaystyle {\int }_{t_0=0}^{t_f}{h}_1}\left(t\right)\hspace{0.17em}\delta \left(t-t_f\right)dt=\phi \left(t_f\right)$$

(31)

$$r_2\left(h\right)={\displaystyle {\int }_{t_0=0}^{t_f}{h}_2}\left(t\right)\hspace{0.17em}\delta \left(t-t_f\right)dt=E\left(t_f\right)$$

(32)

$$r_3\left(h\right)={\displaystyle {\int }_{t_0=0}^{t_f}{h}_3\left(t\right)\hspace{0.17em}\delta \left(t-t_f\right)dt}=T\left(t_f\right)$$

(33)

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:

$$\begin{array}{l}{r}_4\left(h;\phi \right)={\displaystyle {\int }_{t_0}^{t_f}{h}_4^d\left[h\left(t\right);\phi \right]}\hspace{0.17em}\delta \left(t-t_f\right)dt;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\\ h_4^d\left[h\left(t\right);\phi \right]\triangleq k(T)={\phi }_1+{\phi }_2{h}_3\left(t\right)\hspace{0.17em}+{\phi }_3{h}_3^2\left(t\right)\hspace{0.17em}={\phi }_1+{\phi }_{2}T\left(t\right)\hspace{0.17em}+{\phi }_3{T}^2\left(t\right)\hspace{0.17em}\hspace{0.17em}.\end{array}$$

(34)

4. First-Order Comprehensive Adjoint Sensitivity Analysis Methodology for Neural Ordinary Differential Equations (1^st-CASAM-NODE): Mathematical Framework

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):

$${\displaystyle \frac{d{h}^0\left(t\right)}{dt}}=f\left[h^0\left(t\right);{\theta }^0;t\right],\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}t>0$$

(35)

$$h^0\left(t_0\right)=h_e\left(x^0,w^0\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}at\hspace{0.17em}\hspace{0.17em}t=t_0$$

(36)

, 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.:

$$r_n^0\left(h^0;{\phi }^0\right)={\displaystyle {\int }_{t_0}^{t_f}{h}_n^d\left[h^0\left(t\right);{\phi }^0\right]}\hspace{0.17em}\delta \left(t-t_f\right)dt;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}n=1,\mathrm{...},TR.$$

(37)

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:

$$\delta R\left(e^0;\delta e\right)\triangleq {\left\{{\displaystyle \frac{d}{d\epsilon }}\left[R\left(e^0+\epsilon \delta e\right)\right]\right\}}_{\epsilon =0}\triangleq \underset{\epsilon \to 0}{\lim }{\displaystyle \frac{R\left(e^0+\epsilon \delta e\right)-R\left(e^0\right)}{\epsilon }}$$

(38)

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)$:

$$\begin{array}{l}\delta r_n\left(h^0;{\phi }^0;v^{\left(1\right)};\delta \phi \right)={\left\{{\displaystyle \frac{d}{d\epsilon }}{\displaystyle {\int }_{t_0}^{t_f}{h}_n^d\left[h^0\left(t\right)+\epsilon v^{\left(1\right)}\left(t\right);{\phi }^0+\epsilon \delta \phi \right]}\hspace{0.17em}\delta \left(t-t_f\right)dt;\hspace{0.17em}\right\}}_{\epsilon =0}\hspace{0.17em}\hspace{0.17em}\\ ={\left\{\delta r_n\left(h^0;{\phi }^0;\delta \phi \right)\right\}}_{dir}+{\left\{\delta r_n\left(h^0;{\phi }^0;v^{\left(1\right)}\right)\right\}}_{ind}\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}n=1,\mathrm{...},TR.\end{array}$$

(39)

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:

$${\left\{\delta r_n\left(h^0;{\phi }^0;\delta \phi \right)\right\}}_{dir}\triangleq {{\displaystyle {\int }_{t_0}^{t_f}\delta \left(t-t_f\right)\left\{\frac{\partial h_n^d\left[h\left(t\right);\phi \right]}{\partial \phi }\right\}}}_{\left(h^0;{\phi }^0\right)}\delta \phi \hspace{0.17em}\hspace{0.17em}dt\hspace{0.17em}$$

(40)

$${\left\{\delta r_n\left(h^0;{\phi }^0;v^{\left(1\right)}\right)\right\}}_{ind}\triangleq {{\displaystyle {\int }_{t_0}^{t_f}\delta \left(t-t_f\right)\left\{\frac{\partial h_n^d\left[h\left(t\right);\phi \right]}{\partial h\left(t\right)}\right\}}}_{\left(h^0;{\phi }^0\right)}{v}^{\left(1\right)}\left(t\right)\hspace{0.17em}dt\hspace{0.17em}$$

(41)

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:

$${\left\{{\displaystyle \frac{d}{d\epsilon }}\left[{\displaystyle \frac{d}{dt}}\left(h^0+\epsilon v^{\left(1\right)}\right)\right]\right\}}_{\epsilon =0}={\left\{{\displaystyle \frac{d}{d\epsilon }}f\left[h^0+\epsilon v^{\left(1\right)};{\theta }^0+\epsilon \delta \theta ;t\right]\right\}}_{\epsilon =0}$$

(42)

$${\left\{{\displaystyle \frac{d}{d\epsilon }}\left[h^0\left(t_0\right)+\epsilon v^{\left(1\right)}\left(t_0\right)\right]\right\}}_{\epsilon =0}={\left\{{\displaystyle \frac{d}{d\epsilon }}\left[h_e\left(x^0+\epsilon \delta x,w^0+\epsilon \delta w\right)\right]\right\}}_{\epsilon =0}$$

(43)

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

$${\displaystyle \frac{d{v}^{\left(1\right)}\left(t\right)}{dt}}-{\left\{{\displaystyle \frac{\partial f\left(h;\theta \right)}{\partial h}}\right\}}_{\left(h^0,{\theta }^0\right)}{v}^{\left(1\right)}\left(t\right)={\left\{{\displaystyle \frac{\partial f\left(h;\theta \right)}{\partial \theta }}\right\}}_{\left(h^0,{\theta }^0\right)}\delta \theta$$

(44)

$$v^{\left(1\right)}\left(t_0\right)={\left\{{\displaystyle \frac{\partial h_e\left(x,\mathbf{w}\right)}{\partial \mathbf{x}}}\right\}}_{\left(x^0,w^0\right)}\delta x+{\left\{{\displaystyle \frac{\partial h_e\left(x,w\right)}{\partial w}}\right\}}_{\left(x^0,w^0\right)}\delta w$$

(45)

where:

$${\displaystyle \frac{\partial f\left(h;\theta \right)}{\partial h}}\triangleq {\left(\begin{array}{ccc}\partial f_1/\partial h_1& ·& \partial f_1/\partial h_{TH}\\ ·& ·& ·\\ \partial f_{TH}/\partial h_1& ·& \partial f_{TH}/\partial h_{TH}\end{array}\right)}_{TH\times TH}$$

(46)

$${\displaystyle \frac{\partial f\left(h;\theta \right)}{\partial \theta }}\triangleq {\left(\begin{array}{ccc}\partial f_1/\partial {\theta }_1& ·& \partial f_1/\partial {\theta }_{TW}\\ ·& ·& ·\\ \partial f_{TH}/\partial {\theta }_1& ·& \partial f_{TH}/\partial {\theta }_{TW}\end{array}\right)}_{TH\times TW}$$

(47)

$${\displaystyle \frac{\partial h_e\left(x,w\right)}{\partial x}}\triangleq {\left(\begin{array}{ccc}\partial h_1^e/\partial x_1& ·& \partial h_1^e/\partial x_{TI}\\ ·& ·& ·\\ \partial h_{TH}^e/\partial x_1& ·& \partial h_{TH}^e/\partial x_{TI}\end{array}\right)}_{TH\times TI}$$

(48)

$${\displaystyle \frac{\partial h_e\left(x,w\right)}{\partial w}}\triangleq {\left(\begin{array}{ccc}\partial h_1^e/\partial w_1& ·& \partial h_1^e/\partial w_{TEW}\\ ·& ·& ·\\ \partial h_{TH}^e/\partial w_1& ·& \partial h_{TH}^e/\partial w_{TEW}\end{array}\right)}_{TH\times TEW}$$

(49)

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:

$${〈u^{(a)},u^{(b)}〉}_1\triangleq \hspace{0.17em}{\left\{{\displaystyle {\int }_{t_0}^{t_f}{u}^{(a)}\left(t\right)·u^{(b)}\left(t\right)dt}\right\}}_{\left(x^0;{\theta }^0;w^0;{\phi }^0\right)}$$

(50)

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:

$${〈a^{(1)}\left(x\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em} {\displaystyle \frac{d{v}^{\left(1\right)}\left(t\right)}{dt}}-{\left[{\displaystyle \frac{\partial f\left(h;\theta \right)}{\partial h}}\right]}_{\left(h^0,{\theta }^0\right)}{v}^{\left(1\right)}\left(t\right)〉}_1={〈a^{(1)}\left(x\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\left[{\displaystyle \frac{\partial f\left(h;\theta \right)}{\partial \theta }}\right]}_{\left(h^0,{\theta }^0\right)} \delta \theta 〉}_1.$$

(51)

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$:

$$\begin{array}{l}{\displaystyle {\int }_{t_0}^{t_f}{a}^{(1)}\left(t\right)· \frac{d{v}^{\left(1\right)}\left(t\right)}{dt}dt}-{\displaystyle {\int }_{t_0}^{t_f}{a}^{(1)}\left(t\right)·\hspace{0.17em}{\left[\frac{\partial f\left(h;\theta \right)}{\partial h}\right]}_{\left(h^0,{\theta }^0\right)}{v}^{\left(1\right)}\left(t\right)dt}=a^{(1)}\left(t_f\right)· v^{\left(1\right)}\left(t_f\right)\\ -a^{(1)}\left(t_0\right)· v^{\left(1\right)}\left(t_0\right)+{\displaystyle {\int }_{t_0}^{t_f}{v}^{\left(1\right)}\left(t\right)·\left\{-\frac{d{a}^{(1)}\left(t\right)}{dt}-{\left[\frac{\partial f\left(h;\theta \right)}{\partial h}\right]}_{\left(h^0,{\theta }^0\right)}^{†}{a}^{(1)}\left(t\right)\right\} dt}\hspace{0.17em}.\end{array}$$

(52)

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:

$$-{\displaystyle \frac{d{a}^{(1)}\left(t\right)}{dt}}-{\left[{\displaystyle \frac{\partial f\left(h;\theta \right)}{\partial h}}\right]}_{\left(h^0,{\theta }^0\right)}^{†}{a}^{(1)}\left(t\right)={\left\{{\displaystyle \frac{\partial h_n^d\left[h\left(t\right);\phi \right]}{\partial h\left(t\right)}}\right\}}_{\left(h^0;{\phi }^0\right)}\delta \left(t-t_f\right)$$

(53)

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

$$a^{(1)}\left(t_f\right)=0\hspace{0.17em}.$$

(54)

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)$:

$${\left\{\delta r_n\left(h^0;{\phi }^0;v^{\left(1\right)}\right)\right\}}_{ind}={〈a^{(1)}\left(x\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\left[{\displaystyle \frac{\partial f\left(h;\theta \right)}{\partial \theta }}\right]}_{\left(h^0,{\theta }^0\right)} \delta \theta 〉}_1+a^{(1)}\left(t_0\right)· v^{\left(1\right)}\left(t_0\right)$$

(55)

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:

$$\begin{array}{l}{\left\{\delta r_n\left(h^0;{\phi }^0;v^{\left(1\right)}\right)\right\}}_{ind}={〈a^{(1)}\left(x\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\left[{\displaystyle \frac{\partial f\left(h;\theta \right)}{\partial \theta }}\right]}_{\left(h^0,{\theta }^0\right)} \delta \theta 〉}_1\\ +a^{(1)}\left(t_0\right)·{\left\{{\displaystyle \frac{\partial h_e\left(x,w\right)}{\partial x}}\right\}}_{\left(x^0,w^0\right)}\delta x+a^{(1)}\left(t_0\right)·{\left\{{\displaystyle \frac{\partial h_e\left(x,w\right)}{\partial w}}\right\}}_{\left(x^0,w^0\right)}\delta w \end{array}$$

(56)

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)$:

$$\begin{array}{l}\delta r_n\left(h^0;{\phi }^0;v^{\left(1\right)};\delta \phi \right)={\left\{{\displaystyle \frac{\partial h_n^d\left[h\left(t_f\right);\phi \right]}{\partial \phi }}\right\}}_{\left(h^0;{\phi }^0\right)}\delta \phi +a^{(1)}\left(t_0\right)·{\left\{{\displaystyle \frac{\partial h_e\left(x,w\right)}{\partial w}}\right\}}_{\left(x^0,w^0\right)}\delta w\\ +\hspace{0.17em}{a}^{(1)}\left(t_0\right)·{\left\{{\displaystyle \frac{\partial h_e\left(x,w\right)}{\partial x}}\right\}}_{\left(x^0,w^0\right)}\delta x\hspace{0.17em}+{\displaystyle {\int }_{t_0}^{t_f}{a}^{(1)}\left(t\right)·\hspace{0.17em}\left\{{\left[\frac{\partial f\left(h;\theta \right)}{\partial \theta }\right]}_{\left(h^0,{\theta }^0\right)} \delta \theta \right\}}\hspace{0.17em}\hspace{0.17em}dt\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}n=1,\mathrm{...},TR.\end{array}$$

(57)

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:

$${\displaystyle \frac{\partial r_n}{\partial {\phi }_i}}={\displaystyle \frac{\partial h_n^d\left[h\left(t_f\right);\phi \right]}{\partial {\phi }_i}};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},TD;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}n=1,\mathrm{...},TR;$$

(58)

$${\displaystyle \frac{\partial r_n}{\partial w_i}}=a^{(1)}\left(t_0\right)·{\displaystyle \frac{\partial h_e\left(x,w\right)}{\partial w_i}};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},TEW;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}n=1,\mathrm{...},TR;$$

(59)

$${\displaystyle \frac{\partial r_n}{\partial x_i}}=a^{(1)}\left(t_0\right)·{\displaystyle \frac{\partial h_e\left(x,w\right)}{\partial x_i}}\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},TI;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}n=1,\mathrm{...},TR;$$

(60)

$${\displaystyle \frac{\partial r_n}{\partial {\theta }_i}}={\displaystyle {\int }_{t_0}^{t_f}{a}^{(1)}\left(t\right)·\hspace{0.17em}\frac{\partial f\left(h;\theta \right)}{\partial {\theta }_i}}\hspace{0.17em}\hspace{0.17em}dt\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},TW;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}n=1,\mathrm{...},TR.$$

(61)

5. Illustrative Application of the 1^st-CASAM-NODE Methodology to Compute First-Order Sensitivities of Nordheim-Fuchs Model Responses with respect to the Underlying Parameters

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.

5.1. First-Order Sensitivities of the Flux Response $r_1\left(h\right)=\phi \left(t_f\right)$

The first-order sensitivity of the response $r_1\left(h\right)=\phi \left(t_f\right)$ is provided by the first-order G-differential of the expression in Eq. (31), which is, by definition, obtained as follows:

$$\delta r_1\left(h;\delta h\right)={\left\{{\displaystyle \frac{d}{d\epsilon }}\left[{\displaystyle {\int }_0^{t_f}\left[{\phi }^0\left(t\right)+\epsilon \delta \phi \left(t\right)\right]}\hspace{0.17em}\delta \left(t-t_f\right)dt\right]\right\}}_{\epsilon =0}={\displaystyle {\int }_0^{t_f}\delta \phi \left(t\right)\delta \left(t-t_f\right)dt}$$

(62)

The variation $\delta \phi \left(t\right)$ is the solution of the “First-Level Variational Sensitivity System (1^st-LVSS) which is obtained by G-differentiating Eqs. (13)‒(15), which yields the following expressions:

$${\left\{{\displaystyle \frac{d}{d\epsilon }}\left[{\displaystyle \frac{d}{dt}}\left({\phi }^0+\epsilon \delta \phi \right)\right]\right\}}_{\epsilon =0}=-{\left\{{\displaystyle \frac{d}{d\epsilon }}\left[{\displaystyle \frac{{\alpha }_T^0+\epsilon \delta {\alpha }_T}{\left(l_p^0+\epsilon \delta l_p\right)\left(c_p^0+\epsilon \delta c_p\right)}}\left(E^0+\epsilon \delta E\right)\left({\phi }^0+\epsilon \delta \phi \right)\right]\right\}}_{\epsilon =0}\hspace{0.17em},$$

(63)

$${\left\{{\displaystyle \frac{d}{d\epsilon }}\left[{\displaystyle \frac{d}{dt}}\left(E^0+\epsilon \delta E\right)\right]\right\}}_{\epsilon =0}={\left\{{\displaystyle \frac{d}{d\epsilon }}\left[\left({\gamma }^0+\epsilon \delta \gamma \right)\left({\sigma }_f^0+\epsilon \delta {\sigma }_f\right)\left(N_f^0+\epsilon \delta N_f\right)\left({\phi }^0+\epsilon \delta \phi \right)\right]\right\}}_{\epsilon =0},$$

(64)

$${\left\{{\displaystyle \frac{d}{d\epsilon }}\left[{\displaystyle \frac{d}{dt}}\left(T^0+\epsilon \delta T\right)\right]\right\}}_{\epsilon =0}={\left\{{\displaystyle \frac{d}{d\epsilon }}\left[{\displaystyle \frac{\left({\gamma }^0+\epsilon \delta \gamma \right)\left({\sigma }_f^0+\epsilon \delta {\sigma }_f\right)\left(N_f^0+\epsilon \delta N_f\right)}{\left(c_p^0+\epsilon \delta c_p\right)}}\left({\phi }^0+\epsilon \delta \phi \right)\right]\right\}}_{\epsilon =0},$$

(65)

$${\left\{{\displaystyle \frac{d}{d\epsilon }}{\left[{\phi }^0\left(t\right)+\epsilon \delta \phi \left(t\right)\right]}_{t=0}\right\}}_{\epsilon =0}={\left\{{\displaystyle \frac{d}{d\epsilon }}\left({\phi }_0^0+\epsilon \delta {\phi }_0\right)\right\}}_{\epsilon =0}$$

(66)

$${\left\{{\displaystyle \frac{d}{d\epsilon }}{\left[E^0\left(t\right)+\epsilon \delta E\left(t\right)\right]}_{t=0}\right\}}_{\epsilon =0}=0$$

(67)

$${\left\{{\displaystyle \frac{d}{d\epsilon }}{\left[T^0\left(t\right)+\epsilon \delta T\left(t\right)\right]}_{t=0}\right\}}_{\epsilon =0}={\left\{{\displaystyle \frac{d}{d\epsilon }}\left(T_0^0+\epsilon \delta T_0\right)\right\}}_{\epsilon =0}$$

(68)

Performing the operations involving the scalar $\epsilon$ in Eqs. (63)‒(68) yields the following expression for the 1^st-LVSS:

$$\begin{array}{l}{\displaystyle \frac{d}{dt}}\delta \phi \left(t\right)+{\displaystyle \frac{{\alpha }_T^0}{l_p^0{c}_p^0}}{E}^0\left(t\right)\delta \phi \left(t\right)+{\displaystyle \frac{{\alpha }_T^0}{l_p^0{c}_p^0}}{\phi }^0\left(t\right)\delta E\left(t\right)\\ =\left[-{\displaystyle \frac{\delta {\alpha }_T}{l_p^0{c}_p^0}}+{\displaystyle \frac{{\alpha }_T^0}{{\left(l_p^0\right)}^2{c}_p^0}}\delta l_p+{\displaystyle \frac{{\alpha }_T^0}{l_p^0{\left(c_p^0\right)}^2}}\delta c_p\right]E^0\left(t\right){\phi }^0\left(t\right)\hspace{0.17em}\hspace{0.17em},\end{array}$$

(69)

$${\displaystyle \frac{d}{dt}}\delta E\left(t\right)-{\gamma }^0{\sigma }_f^0{N}_f^0\delta \phi \left(t\right)=\left[\left({\sigma }_f^0{N}_f^0\right)\delta \gamma +\left({\gamma }^0{N}_f^0\right)\delta {\sigma }_f+\left({\gamma }^0{\sigma }_f^0\right)\delta N_f\right]{\phi }^0\left(t\right),$$

(70)

$$\begin{array}{l}{\displaystyle \frac{d}{dt}}\delta T\left(t\right)-{\displaystyle \frac{{\gamma }^0{\sigma }_f^0{N}_f^0}{c_p^0}}\delta \phi \left(t\right)\\ =\left[\left({\sigma }_f^0{N}_f^0\right)\delta \gamma +\left({\gamma }^0{N}_f^0\right)\delta {\sigma }_f+\left({\gamma }^0{\sigma }_f^0\right)\delta N_f-{\displaystyle \frac{{\gamma }^0{\sigma }_f^0{N}_f^0}{c_p^0}}\delta c_p\right]{\displaystyle \frac{{\phi }^0\left(t\right)}{c_p^0}},\end{array}$$

(71)

$${\left[\delta \phi \left(t\right)\right]}_{t=0}=\delta {\phi }_0$$

(72)

$${\left[\delta E\left(t\right)\right]}_{t=0}=0$$

(73)

$${\left[\delta T\left(t\right)\right]}_{t=0}=\delta T_0$$

(74)

The 1^st-LVSS comprising Eqs. (69)‒(74) represents the specific form taken on by the general NODE-representation of the 1^st-LVSS provided by Eqs. (44) and (45) for the Nordheim-Fuchs model. Comparing Eqs. (69)‒(74) to Eqs. (44) and (45) indicates the following correspondences:

$${\displaystyle \frac{\partial f\left(h;\theta \right)}{\partial h}}\triangleq \left(\begin{array}{ccc}-{\displaystyle \frac{{\alpha }_T^0{E}^0\left(t\right)}{l_p^0{c}_p^0}}& -{\displaystyle \frac{{\alpha }_T^0{\phi }^0\left(t\right)}{l_p^0{c}_p^0}}& 0\\ {\gamma }^0{\sigma }_f^0{N}_f^0& 0& 0\\ {\displaystyle \frac{{\gamma }^0{\sigma }_f^0{N}_f^0}{c_p^0}}& 0& 0\end{array}\right);\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{v}^{\left(1\right)}\left(t\right)\triangleq \left(\begin{array}{c}\delta \phi \left(t\right)\\ \delta E\left(t\right)\\ \delta T\left(t\right)\end{array}\right);\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}$$

(75)

$${\displaystyle \frac{\partial f\left(h;\theta \right)}{\partial \theta }}\triangleq \left(\begin{array}{cccccc}{\displaystyle \frac{\partial f_1}{\partial {\theta }_1}}& {\displaystyle \frac{\partial f_1}{\partial {\theta }_2}}& {\displaystyle \frac{\partial f_1}{\partial {\theta }_3}}& 0& 0& 0\\ 0& 0& 0& {\displaystyle \frac{\partial f_2}{\partial {\theta }_4}}& {\displaystyle \frac{\partial f_2}{\partial {\theta }_5}}& {\displaystyle \frac{\partial f_2}{\partial {\theta }_6}}\\ 0& 0& {\displaystyle \frac{\partial f_3}{\partial {\theta }_3}}& {\displaystyle \frac{\partial f_3}{\partial {\theta }_4}}& {\displaystyle \frac{\partial f_3}{\partial {\theta }_5}}& {\displaystyle \frac{\partial f_3}{\partial {\theta }_6}}\end{array}\right);\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\delta \theta \triangleq \left(\begin{array}{c}\delta \alpha \\ \delta l_p\\ \delta c_p\\ \delta \gamma \\ \delta {\sigma }_f\\ \delta N_f\end{array}\right);$$

(76)

$$\begin{array}{l}{\displaystyle \frac{\partial f_1}{\partial {\theta }_1}}\triangleq -{\displaystyle \frac{E^0\left(t\right){\phi }^0\left(t\right)}{l_p^0{c}_p^0}};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\displaystyle \frac{\partial f_1}{\partial {\theta }_2}}\triangleq {\displaystyle \frac{{\alpha }_T^0{E}^0\left(t\right){\phi }^0\left(t\right)}{{\left(l_p^0\right)}^2{c}_p^0}};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\displaystyle \frac{\partial f_1}{\partial {\theta }_3}}\triangleq {\displaystyle \frac{{\alpha }_T^0{E}^0\left(t\right){\phi }^0\left(t\right)}{l_p^0{\left(c_p^0\right)}^2}};\\ {\displaystyle \frac{\partial f_2}{\partial {\theta }_4}}\triangleq {\sigma }_f^0{N}_f^0{\phi }^0\left(t\right);\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\displaystyle \frac{\partial f_2}{\partial {\theta }_5}}\triangleq {\gamma }^0{N}_f^0{\phi }^0\left(t\right);\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\displaystyle \frac{\partial f_2}{\partial {\theta }_6}}\triangleq {\gamma }^0{\sigma }_f^0{\phi }^0\left(t\right);\\ \hspace{0.17em}{\displaystyle \frac{\partial f_3}{\partial {\theta }_3}}\triangleq -{\displaystyle \frac{{\gamma }^0{\sigma }_f^0{N}_f^0{\phi }^0\left(t\right)}{{\left(c_p^0\right)}^2}};\hspace{0.17em}\hspace{0.17em}{\displaystyle \frac{\partial f_3}{\partial {\theta }_4}}\triangleq {\displaystyle \frac{{\sigma }_f^0{N}_f^0{\phi }^0\left(t\right)}{c_p^0}};\hspace{0.17em}\hspace{0.17em}{\displaystyle \frac{\partial f_3}{\partial {\theta }_5}}\triangleq {\displaystyle \frac{{\gamma }^0{N}_f^0{\phi }^0\left(t\right)}{c_p^0}};\hspace{0.17em}\hspace{0.17em}{\displaystyle \frac{\partial f_3}{\partial {\theta }_6}}\triangleq {\displaystyle \frac{{\gamma }^0{\sigma }_f^0{\phi }^0\left(t\right)}{c_p^0}}.\end{array}$$

(77)

$${\displaystyle \frac{\partial h_e\left(x,w\right)}{\partial x}}\triangleq \left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right)\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\delta x\triangleq \left(\begin{array}{c}\delta {\phi }_0\\ 0\\ \delta T_0\end{array}\right);\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\displaystyle \frac{\partial h_e\left(x,w\right)}{\partial w}}\triangleq \left(\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right)\hspace{0.17em}\hspace{0.17em}.$$

(78)

It is evident that the 1^st-LVSS would need to be solved repeatedly in order to compute the 1^st-level variational function $v^{\left(1\right)}\left(t\right)\triangleq {\left[\delta \phi \left(t\right),\delta E\left(t\right),\delta T\left(t\right)\right]}^{†}$ for every possible variations $\delta \theta$ in the model parameters and variations $\delta x$ in the initial conditions (‘encoder”). This computationally expensive path can be avoided by applying the concepts of the 1^st- CASAM-NODE previously outlined in Subsection 4.1, as follows:

• Consider that the 1^st-level variational function $v^{\left(1\right)}\triangleq {\left[\delta \phi \left(t\right),\delta E\left(t\right),\delta T\left(t\right)\right]}^{†}\in H_1\left({\Omega }_t\right)$, is an element in a Hilbert space denoted as $H_1\left({\Omega }_t\right)$, $\hspace{0.17em}{\Omega }_t\triangleq \left(0,t_f\right)$, comprising elements of the form $u^{(a)}\left(t\right)\triangleq {\left[u_1^{(a)}\left(t\right),u_2^{(a)}\left(t\right),u_3^{(a)}\left(t\right)\right]}^{†}$, $u^{(b)}\left(t\right)\triangleq {\left[u_1^{(b)}\left(t\right),u_2^{(b)}\left(t\right),u_3^{(b)}\left(t\right)\right]}^{†}$, and being endowed with the inner product ${〈u^{(a)},u^{(b)}〉}_1$ introduced in Eq. (50), which takes on the following particular form for the Nordheim-Fuchs model:

$${〈u^{(a)},u^{(b)}〉}_1\triangleq \hspace{0.17em}{\displaystyle {\int }_0^{t_f}{u}^{(a)}\left(t\right)·u^{(b)}\left(t\right)dt}={\displaystyle \sum _{i=1}^3{\displaystyle {\int }_0^{t_f}{u}_i^{(a)}\left(t\right)u_i^{(b)}\left(t\right)dt}}$$

(79)

2.

Use Eq. (79) to form the inner product of Eqs. (69)‒(71) with a yet undefined function $a^{(1)}\left(t\right)\triangleq {\left[a_1^{(1)}\left(t\right),a_2^{(1)}\left(t\right),a_3^{(1)}\left(t\right)\right]}^{†}\in H_1\left({\Omega }_t\right)$, to obtain the following relation, which is the particular form taken on by Eq. (51) for the Nordheim-Fuchs model:

$$\begin{array}{l}{\displaystyle {\int }_0^{t_f}{a}_1^{(1)}\left(t\right)\left[\frac{d}{dt}\delta \phi \left(t\right)+\frac{{\alpha }_T^0}{l_p^0{c}_p^0}{E}^0\left(t\right)\delta \phi \left(t\right)+\frac{{\alpha }_T^0}{l_p^0{c}_p^0}{\phi }^0\left(t\right)\delta E\left(t\right)\right]dt}\\ +{\displaystyle {\int }_0^{t_f}{a}_2^{(1)}\left(t\right)}\left[{\displaystyle \frac{d}{dt}}\delta E\left(t\right)-{\gamma }^0{\sigma }_f^0{N}_f^0\delta \phi \left(t\right)\right]dt\\ +{\displaystyle {\int }_0^{t_f}{a}_3^{(1)}\left(t\right)}\left[{\displaystyle \frac{d}{dt}}\delta T\left(t\right)-{\displaystyle \frac{{\gamma }^0{\sigma }_f^0{N}_f^0}{c_p^0}}\delta \phi \left(t\right)\right]dt\\ ={\displaystyle {\int }_0^{t_f}{a}_1^{(1)}\left(t\right)\left[-\frac{\delta {\alpha }_T}{l_p^0{c}_p^0}+\frac{{\alpha }_T^0}{{\left(l_p^0\right)}^2{c}_p^0}\delta l_p+\frac{{\alpha }_T^0}{l_p^0{\left(c_p^0\right)}^2}\delta c_p\right]E^0\left(t\right){\phi }^0\left(t\right)\hspace{0.17em}}\hspace{0.17em}dt\\ +{\displaystyle {\int }_0^{t_f}{a}_2^{(1)}\left(t\right)}\left[\left({\sigma }_f^0{N}_f^0\right)\delta \gamma +\left({\gamma }^0{N}_f^0\right)\delta {\sigma }_f+\left({\gamma }^0{\sigma }_f^0\right)\delta N_f\right]{\phi }^0\left(t\right)dt\\ +{\displaystyle {\int }_0^{t_f}{a}_3^{(1)}\left(t\right)}\left[\left({\sigma }_f^0{N}_f^0\right)\delta \gamma +\left({\gamma }^0{N}_f^0\right)\delta {\sigma }_f+\left({\gamma }^0{\sigma }_f^0\right)\delta N_f-{\displaystyle \frac{{\gamma }^0{\sigma }_f^0{N}_f^0}{c_p^0}}\delta c_p\right]{\displaystyle \frac{{\phi }^0\left(t\right)}{c_p^0}}dt.\end{array}$$

(80)

3.

Integrating by parts the terms on the left-side of Eq. (80) yields the following relation

$$\begin{array}{l}{\displaystyle {\int }_0^{t_f}{a}_1^{(1)}\left(t\right)\left[\frac{d}{dt}\delta \phi \left(t\right)+\frac{{\alpha }_T^0}{l_p^0{c}_p^0}{E}^0\left(t\right)\delta \phi \left(t\right)+\frac{{\alpha }_T^0}{l_p^0{c}_p^0}{\phi }^0\left(t\right)\delta E\left(t\right)\right]dt}\\ +{\displaystyle {\int }_0^{t_f}{a}_2^{(1)}\left(t\right)}\left[{\displaystyle \frac{d}{dt}}\delta E\left(t\right)-{\gamma }^0{\sigma }_f^0{N}_f^0\delta \phi \left(t\right)\right]dt\\ +{\displaystyle {\int }_0^{t_f}{a}_3^{(1)}\left(t\right)}\left[{\displaystyle \frac{d}{dt}}\delta T\left(t\right)-{\displaystyle \frac{{\gamma }^0{\sigma }_f^0{N}_f^0}{c_p^0}}\delta \phi \left(t\right)\right]dt=a_1^{(1)}\left(t_f\right)\delta \phi \left(t_f\right)-a_1^{(1)}\left(0\right)\delta \phi \left(0\right)\\ +a_2^{(1)}\left(t_f\right)\delta E\left(t_f\right)-a_2^{(1)}\left(0\right)\delta E\left(0\right)+a_3^{(1)}\left(t_f\right)\delta T\left(t_f\right)-a_3^{(1)}\left(0\right)\delta T\left(0\right)\\ +{\displaystyle {\int }_0^{t_f}{v}^{\left(1\right)}\left(t\right)·{\left\{A^{\left(1\right)}\left(h;\theta \right)a^{(1)}\left(t\right)\right\}}_{\left(h^0;{\theta }^0\right)}dt}\hspace{0.17em},\end{array}$$

(81)

where:

$$A^{\left(1\right)}\left(h;\theta \right)a^{(1)}\left(t\right)\triangleq -{\displaystyle \frac{d{a}^{(1)}\left(t\right)}{dt}}-{\left[{\displaystyle \frac{\partial f\left(h;\theta \right)}{\partial h}}\right]}_{\left(h^0,{\theta }^0\right)}^{†}{a}^{(1)}\left(t\right)$$

(82)

With

$${\left[{\displaystyle \frac{\partial f\left(h;\theta \right)}{\partial h}}\right]}_{\left(h^0,{\theta }^0\right)}^{†}\triangleq \left(\begin{array}{ccc}-{\alpha }_T^0{E}^0\left(t\right)/\left(l_p^0{c}_p^0\right)& {\gamma }^0{\sigma }_f^0{N}_f^0& {\gamma }^0{\sigma }_f^0{N}_f^0/c_p^0\\ -{\alpha }_T^0{\phi }^0\left(t\right)/\left(l_p^0{c}_p^0\right)& 0& 0\\ 0& 0& 0\end{array}\right)\hspace{0.17em}.$$

(83)

The relation obtained in Eq. (81) is the particular form taken on by Eq. (52) for the Nordheim-Fuchs model.

4.

The definition of the function $a^{(1)}\left(t\right)$ is now completed by requiring that: (i) the integral term on the right-side of Eq.(81) represent the G-differential $\delta r_1\left(h;\delta h\right)$ defined in Eq. (62), and (ii) the appearance of the unknown values of the components of $v^{\left(1\right)}\left(t_f\right)$ be eliminated from appearing in Eq. (81). These requirements will be satisfied if the function $a^{(1)}\left(t\right)\triangleq {\left[a_1^{(1)}\left(t\right),a_2^{(1)}\left(t\right),a_3^{(1)}\left(t\right)\right]}^{†}\in H_1\left({\Omega }_t\right)$ is the solution of the following “1^st-Level Adjoint Sensitivity System (1^st-LASS)”:

$$A^{\left(1\right)}\left(h;\theta \right)a^{(1)}\left(t\right)\triangleq -{\displaystyle \frac{d{a}^{(1)}\left(t\right)}{dt}}-{\left[{\displaystyle \frac{\partial f\left(h;\theta \right)}{\partial h}}\right]}_{\left(h^0,{\theta }^0\right)}^{†}{a}^{(1)}\left(t\right)={\left[\delta \left(t-t_f\right),0,0\right]}^{†}$$

(84)

$$a^{(1)}\left(t_f\right)\triangleq {\left[a_1^{(1)}\left(t_f\right),a_2^{(1)}\left(t_f\right),a_3^{(1)}\left(t_f\right)\right]}^{†}={\left[0,0,0\right]}^{†}$$

(85)

It is important to note that if the vector-valued function $f\left(h;\theta \right)$ is linear in $h\left(t\right)$ (in which case the NODE would be linear), then the 1^st-level adjoint sensitivity function $a^{(1)}\left(t\right)$ would not depend on $h\left(t\right)$, so the “forward solution path” would not need to be stored in order to compute $a^{(1)}\left(t\right)$. Otherwise, however, the “forward solution path” $h\left(t\right)$ would need to be stored in order to compute $a^{(1)}\left(t\right)$.

5.

Using Eqs. (84), (85), (80), (62), (72), (73) and (74) in Eq. (81) yields the following expression for the first G-differential $\delta r_1\left(h;\delta h\right)$ of the response under consideration:

$$\begin{array}{l}\delta r_1\left(h;\delta h\right)=\delta \phi \left(t_f\right)={\displaystyle {\int }_0^{t_f}{a}_1^{(1)}\left(t\right)\left[-\frac{\delta {\alpha }_T}{l_p^0{c}_p^0}+\frac{{\alpha }_T^0}{{\left(l_p^0\right)}^2{c}_p^0}\delta l_p+\frac{{\alpha }_T^0}{l_p^0{\left(c_p^0\right)}^2}\delta c_p\right]E^0\left(t\right){\phi }^0\left(t\right)\hspace{0.17em}}\hspace{0.17em}dt\\ +{\displaystyle {\int }_0^{t_f}{a}_2^{(1)}\left(t\right)}\left[\left({\sigma }_f^0{N}_f^0\right)\delta \gamma +\left({\gamma }^0{N}_f^0\right)\delta {\sigma }_f+\left({\gamma }^0{\sigma }_f^0\right)\delta N_f\right]{\phi }^0\left(t\right)dt\\ +{\displaystyle {\int }_0^{t_f}{a}_3^{(1)}\left(t\right)}\left[\left({\sigma }_f^0{N}_f^0\right)\delta \gamma +\left({\gamma }^0{N}_f^0\right)\delta {\sigma }_f+\left({\gamma }^0{\sigma }_f^0\right)\delta N_f-{\displaystyle \frac{{\gamma }^0{\sigma }_f^0{N}_f^0}{c_p^0}}\delta c_p\right]{\displaystyle \frac{{\phi }^0\left(t\right)}{c_p^0}}dt\\ +a_1^{(1)}\left(0\right)\delta {\phi }_0+a_3^{(1)}\left(0\right)\delta T_0\hspace{0.17em}.\end{array}$$

(86)

It follows from Eq. (86) that the first-order sensitivities of the response $\phi \left(t_f\right)$ with respect to the parameters and initial conditions underlying the Nordheim-Fuchs model have the following expressions, all of which are to be evaluated at the nominal values of the respective parameters and functions (but the superscript “zero” is omitted to simplify the notation):

$${\displaystyle \frac{\partial \phi \left(t_f\right)}{\partial {\alpha }_T}}=-{\displaystyle \frac{1}{l_p{c}_p}}{\displaystyle {\int }_0^{t_f}{a}_1^{(1)}\left(t\right)E\left(t\right)\phi \left(t\right)\hspace{0.17em}}\hspace{0.17em}dt$$

(87)

$${\displaystyle \frac{\partial \phi \left(t_f\right)}{\partial l_p}}={\displaystyle \frac{{\alpha }_T}{{\left(l_p\right)}^2{c}_p}}{\displaystyle {\int }_0^{t_f}{a}_1^{(1)}\left(t\right)E\left(t\right)\phi \left(t\right)\hspace{0.17em}}\hspace{0.17em}dt$$

(88)

$${\displaystyle \frac{\partial \phi \left(t_f\right)}{\partial c_p}}={\displaystyle \frac{{\alpha }_T}{l_p{\left(c_p\right)}^2}}{\displaystyle {\int }_0^{t_f}{a}_1^{(1)}\left(t\right)E\left(t\right)\phi \left(t\right)\hspace{0.17em}}\hspace{0.17em}dt-{\displaystyle \frac{\gamma {\sigma }_f{N}_f}{{\left(c_p\right)}^2}}{\displaystyle {\int }_0^{t_f}{a}_3^{(1)}\left(t\right)}\hspace{0.17em}\phi \left(t\right)dt$$

(89)

$${\displaystyle \frac{\partial \phi \left(t_f\right)}{\partial \gamma }}={\sigma }_f{N}_f{\displaystyle {\int }_0^{t_f}\left[a_2^{(1)}\left(t\right)+\frac{1}{c_p}{a}_3^{(1)}\left(t\right)\right]}\hspace{0.17em}\phi \left(t\right)dt$$

(90)

$${\displaystyle \frac{\partial \phi \left(t_f\right)}{\partial {\sigma }_f}}=\gamma N_f{\displaystyle {\int }_0^{t_f}\left[a_2^{(1)}\left(t\right)+\frac{1}{c_p}{a}_3^{(1)}\left(t\right)\right]}\hspace{0.17em}\phi \left(t\right)dt$$

(91)

$${\displaystyle \frac{\partial \phi \left(t_f\right)}{\partial N_f}}=\gamma {\sigma }_f{\displaystyle {\int }_0^{t_f}\left[a_2^{(1)}\left(t\right)+\frac{1}{c_p}{a}_3^{(1)}\left(t\right)\right]}\hspace{0.17em}\phi \left(t\right)dt$$

(92)

$${\displaystyle \frac{\partial \phi \left(t_f\right)}{\partial {\phi }_0}}=a_1^{(1)}\left(0\right);\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\displaystyle \frac{\partial \phi \left(t_f\right)}{\partial E\left(0\right)}}=0;\hspace{0.17em}\hspace{0.17em}{\displaystyle \frac{\partial \phi \left(t_f\right)}{\partial T_0}}=a_3^{(1)}\left(0\right)$$

(93)

5.2. First-Order Sensitivities of the Energy Released Response $r_2\left(h\right)=E\left(t_f\right)$

The first-order G-differential of the response $r_2\left(h\right)=E\left(t_f\right)$ defined in Eq. (32) is obtained as follows:

$$\delta r_2\left(h;\delta h\right)={\left\{{\displaystyle \frac{d}{d\epsilon }}\left[{\displaystyle {\int }_0^{t_f}\left[E^0\left(t\right)+\epsilon \delta E\left(t\right)\right]}\hspace{0.17em}\delta \left(t-t_f\right)dt\right]\right\}}_{\epsilon =0}={\displaystyle {\int }_0^{t_f}\delta E\left(t\right)\delta \left(t-t_f\right)dt}$$

(94)

where the variation $\delta E\left(t\right)$ is the solution of the “First-Level Variational Sensitivity System (1^st-LVSS) defined by Eqs. (69)‒(74).

The sensitivities of the response $r_2\left(h\right)=E\left(t_f\right)$ are determined by following the same procedure as has been outlined in Section 5.2, using an adjoint function denoted as ${\chi }^{(1)}\left(t\right)\triangleq {\left[{\chi }_1^{(1)}\left(t\right),{\chi }_2^{(1)}\left(t\right),{\chi }_3^{(1)}\left(t\right)\right]}^{†}\in H_1\left({\Omega }_t\right)$. Following the same steps as in Section 5.2 (which are omitted here to avoid undue repetition) leads to the following 1^st-LASS for the 1^st-level adjoint sensitivity function ${\chi }^{(1)}\left(t\right)$:

$$A^{\left(1\right)}\left(h;\theta \right){\chi }^{(1)}\left(t\right)\triangleq -{\displaystyle \frac{d{\chi }^{(1)}\left(t\right)}{dt}}-{\left[{\displaystyle \frac{\partial f\left(h;\theta \right)}{\partial h}}\right]}_{\left(h^0,{\theta }^0\right)}^{†}{\chi }^{(1)}\left(t\right)={\left[0,\delta \left(t-t_f\right),0\right]}^{†}$$

(95)

$${\chi }^{(1)}\left(t_f\right)\triangleq {\left[{\chi }_1^{(1)}\left(t_f\right),{\chi }_2^{(1)}\left(t_f\right),{\chi }_3^{(1)}\left(t_f\right)\right]}^{†}={\left[0,0,0\right]}^{†}$$

(96)

The sensitivities of $E\left(t_f\right)$ with respect to the model parameters and initial conditions have the same formal expressions as shown in Eqs. (87)‒(93), but with the components of the 1^st-level adjoint sensitivity function ${\chi }^{(1)}\left(t\right)$ replacing the components of $a^{(1)}\left(t\right)\triangleq {\left[a_1^{(1)}\left(t\right),a_2^{(1)}\left(t\right),a_3^{(1)}\left(t\right)\right]}^{†}.$

5.3. First-Order Sensitivities of the Temperature Response $r_3\left(h\right)=T\left(t_f\right)$

The first-order G-differential of the response $r_3\left(h\right)=T\left(t_f\right)$ defined in Eq. (33) is obtained as follows:

$$\delta r_3\left(h;\delta h\right)={\left\{{\displaystyle \frac{d}{d\epsilon }}\left[{\displaystyle {\int }_0^{t_f}\left[T^0\left(t\right)+\epsilon \delta T\left(t\right)\right]}\hspace{0.17em}\delta \left(t-t_f\right)dt\right]\right\}}_{\epsilon =0}={\displaystyle {\int }_0^{t_f}\delta T\left(t\right)\delta \left(t-t_f\right)dt}$$

(97)

where the variation $\delta T\left(t\right)$ is the solution of the “First-Level Variational Sensitivity System (1^st-LVSS) defined by Eqs. (69)‒(74).

The sensitivities of the response $r_3\left(h\right)=T\left(t_f\right)$ are determined by following the same procedure as has been outlined in Section 5.1, using an adjoint function denoted as${\xi }^{(1)}\left(t\right)\triangleq {\left[{\xi }_1^{(1)}\left(t\right),{\xi }_2^{(1)}\left(t\right),{\xi }_3^{(1)}\left(t\right)\right]}^{†}\in H_1\left({\Omega }_t\right)$. Following the same steps as in Section 5.1 (which are omitted here to avoid undue repetition) leads to the following 1^st-LASS for the 1^st-level adjoint sensitivity function ${\xi }^{(1)}\left(t\right)$:

$$A^{\left(1\right)}\left(h;\theta \right){\xi }^{(1)}\left(t\right)\triangleq -{\displaystyle \frac{d{\xi }^{(1)}\left(t\right)}{dt}}-{\left[{\displaystyle \frac{\partial f\left(h;\theta \right)}{\partial h}}\right]}_{\left(h^0,{\theta }^0\right)}^{†}{\xi }^{(1)}\left(t\right)={\left[0,0,\delta \left(t-t_f\right)\right]}^{†}$$

(98)

$${\xi }^{(1)}\left(t_f\right)\triangleq {\left[{\xi }_1^{(1)}\left(t_f\right),{\xi }_2^{(1)}\left(t_f\right),{\xi }_3^{(1)}\left(t_f\right)\right]}^{†}={\left[0,0,0\right]}^{†}$$

(99)

The sensitivities of $T\left(t_f\right)$ with respect to the model parameters and initial conditions have the same formal expressions as shown in Eqs. (87)‒(93), but with the components of the 1^st-level adjoint sensitivity function ${\xi }^{(1)}\left(t\right)\triangleq {\left[{\xi }_1^{(1)}\left(t\right),{\xi }_2^{(1)}\left(t\right),{\xi }_3^{(1)}\left(t\right)\right]}^{†}$ replacing the components of $a^{(1)}\left(t\right)\triangleq {\left[a_1^{(1)}\left(t\right),a_2^{(1)}\left(t\right),a_3^{(1)}\left(t\right)\right]}^{†}.$

5.4. First-Order Sensitivities of the Thermal Conductivity Response $r_4\left(h;\phi \right)=k(T_f;\phi )$

The first-order G-differential of the response $r_4\left(h;\phi \right)=k(T_f;\phi )$ defined in Eq. (34) is obtained as follows:

$$\begin{array}{l}\delta r_4\left(h;\phi ;\delta h;\delta \phi \right)=\delta k(T;\phi ;\delta T;\delta \phi )={\left\{{\displaystyle \frac{d}{d\epsilon }}{\displaystyle {\int }_0^{t_f}\left[{\phi }_1^0+\epsilon \delta {\phi }_1\right]}\hspace{0.17em}\delta \left(t-t_f\right)dt\right\}}_{\epsilon =0}\\ +{\left\{{\displaystyle \frac{d}{d\epsilon }}{\displaystyle {\int }_0^{t_f}\left[\left({\phi }_2^0+\epsilon \delta {\phi }_2\right)\left(T^0+\epsilon \delta T\right)\hspace{0.17em}+\left({\phi }_3^0+\epsilon \delta {\phi }_3\right){\left(T^0+\epsilon \delta T\right)}^2\right]}\hspace{0.17em}\delta \left(t-t_f\right)dt\right\}}_{\epsilon =0}\\ ={\left\{\delta k(T;\phi ;\delta \phi )\right\}}_{dir}+{\left\{\delta k(T;\phi ;\delta T)\right\}}_{ind},\end{array}$$

(100)

where the direct-effect and the indirect-effect terms, respectively, are defined as follows:

$${\left\{\delta k(T;\phi ;\delta \phi )\right\}}_{dir}\triangleq \delta {\phi }_1+\delta {\phi }_2{\displaystyle {\int }_0^{t_f}{T}^0\left(t\right)}\hspace{0.17em}\delta \left(t-t_f\right)+\hspace{0.17em}\hspace{0.17em}\delta {\phi }_3{\displaystyle {\int }_0^{t_f}\left[\hspace{0.17em}{T}^0\left(t\right)\right]}{\hspace{0.17em}}^2\delta \left(t-t_f\right)$$

(101)

$${\left\{\delta k(T;\phi ;\delta T)\right\}}_{ind}\triangleq {\displaystyle {\int }_0^{t_f}\left[{\phi }_2^0\hspace{0.17em}+2{\phi }_3^0{T}^0\left(t\right)\right]\delta T\left(t\right)}\hspace{0.17em}\delta \left(t-t_f\right)dt$$

(102)

The direct-effect term yields the following sensitivities which can be evaluated immediately:

$${\displaystyle \frac{\partial k(T_f)}{\partial {\phi }_1}}=1;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\displaystyle \frac{\partial k(T_f)}{\partial {\phi }_2}}=T^0\left(t_f\right);\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\displaystyle \frac{\partial k(T_f)}{\partial {\phi }_3}}={\left[\hspace{0.17em}{T}^0\left(t_f\right)\right]}^2$$

(103)

The indirect-effect term can be evaluated only after determining the variational function $\delta T\left(t\right)$, which is the solution of the 1^st-LVSS defined by Eqs. (69)‒(74). The need for solving (repeatedly) the 1^st-LVSS can be circumvented by applying the principles of the 1^st-CASAM-NODE, as previously outlined. Thus, following the same procedure as detailed in Section 5.1 leads to the following 1^st-LASS for the 1^st-level adjoint sensitivity function, denoted as ${\psi }^{(1)}\left(t\right)\triangleq {\left[{\psi }_1^{(1)}\left(t\right),{\psi }_2^{(1)}\left(t\right),{\psi }_3^{(1)}\left(t\right)\right]}^{†}\in H_1\left({\Omega }_t\right)$, for computing the sensitivities stemming from the indirect-effect term ${\left\{\delta k(T;\phi ;\delta T)\right\}}_{ind}$:

$$A^{\left(1\right)}\left(h;\theta \right){\psi }^{(1)}\left(t\right)={\left[\left[{\phi }_2^0\hspace{0.17em}+2{\phi }_3^0{T}^0\left(t\right)\right]\delta \left(t-t_f\right),0,0\right]}^{†}$$

(104)

$${\psi }^{(1)}\left(t_f\right)\triangleq {\left[{\psi }_1^{(1)}\left(t_f\right),{\psi }_2^{(1)}\left(t_f\right),{\psi }_3^{(1)}\left(t_f\right)\right]}^{†}={\left[0,0,0\right]}^{†}$$

(105)

It is important to note that all of the following 1^st-Level Adjoint Sensitivity Systems, enumerated in items (i) through (iv), below:

the 1^st-LASS defined by Eqs. (84) and (85), which are solved for obtaining the corresponding 1^st-level adjoint sensitivity function needed for computing the sensitivities of the component $h_1\left(t\right)\triangleq \phi \left(t\right)$ of the state function $h\left(t\right)$;

the 1^st-LASS defined by Eqs. (95) and (96), which are solved for obtaining the corresponding 1^st-level adjoint sensitivity function needed for computing the sensitivities of the component $h_2\left(t\right)\triangleq E\left(t\right)$ of the state function $h\left(t\right)$];

the 1^st-LASS defined by Eqs. (98) and (99), which are solved for obtaining the corresponding 1^st-level adjoint sensitivity function needed for computing the sensitivities of the component $h_3\left(t\right)\triangleq T\left(t\right)$ of the state function $h\left(t\right)$, and

the 1^st-LASS defined by Eqs. (104) and (105), which are solved for obtaining the corresponding 1^st-level adjoint sensitivity function needed for computing the sensitivities stemming from the indirect-effect term ${\left\{\delta k(T;\phi ;\delta T)\right\}}_{ind}$;

…have the same structures/operators on their left sides, and the respective adjoint sensitivity function all satisfy the same final-time conditions; only the source terms on the right-sides of the respective 1^st-LASS differ from each other. Consequently, the same numerical procedures and/or neural nets can be used for computing the respective 1^st-level adjoint sensitivity functions.

Since the NODE is a first-order ODE, the corresponding 1^st-LASS is solved “backwards” in time, starting at the final time-step $t=t_f$, as indicated by the general 1^st-CASAM-NODE methodology presented in Section 4. If the NODE is linear in the state function (dependent variable) $h\left(t\right)$, then the 1^st-LASS will be independent of $h\left(t\right)$, so the “forward solution path” would not need to be stored in order to compute the 1^st-level adjoint sensitivity functions. In contradistinction, if the NODE is nonlinear in the state function (dependent variable) $h\left(t\right)$, then the 1^st-LASS will depend on $h\left(t\right)$, so the “forward solution path” would need to be stored in order to compute the respective 1^st-level adjoint sensitivity functions.

Furthermore, the same formal expressions are obtained for the sensitivities of the responses considered. Thus, the respective 1^st-level adjoint sensitivity functions differ from each other according to the response considered, but the quadrature-schemes needed to evaluate the integrals defining the respective sensitivities are the same. Therefore, the same numerical procedures and/or neural nets can be used for computing the respective integrals that define the 1^st-order sensitivities, while using the appropriate/corresponding 1^st-level adjoint sensitivity functions. If the decoder-response depends on parameters/weights, additional sensitivities arise from the respective non-vanishing “direct-effect term.”

If simple relations can be obtained among the responses of interest, such as Eqs. (11) and (16) for the illustrative paradigm example, then the sensitivities of the various responses can be obtained by using these relationships, but this is seldom the case in practice.

5.5. Most Efficient Computation of First-Order Sensitivities: Application of the 1^st-FASAM-N

In most, if not all, practical situations, the equations modeling the physical system under consideration can be recast to suit the computation of the response under consideration and, consequently, the computation of the response sensitivities with respect to the underlying model parameters. For example, the response $r_2\left(h\right)=E\left(t_f\right)$ involves just the function $E\left(t\right)$; hence, this response would be ideally computed, together with its sensitivities to parameters, by using an equation containing as few as possible dependent variables other than the ones [e.g., $E\left(t\right)$] needed for computing the response. Such an equation was obtained in Eq. (17), which contains just the dependent variable $E\left(t\right)$, so it would be more advantageous to us it for the sensitivity analysis of $r_2\left(h\right)=E\left(t_f\right)$ rather than use the entire system of equations underlying the Nordheim-Fuchs model, as was done, for illustrative purposes, in Section 5.1. Furthermore, the form of Eq. (17) indicates that the “features” (i.e., functions) of model parameters characterizing this balance equation can be chosen as follows:

$$F_1\left(p\right)\triangleq {\displaystyle \frac{{\alpha }_T}{2l_p{c}_p}};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{F}_2\left(p\right)\triangleq {\phi }_0\gamma {\sigma }_f{N}_f;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}F\left(p\right)\triangleq {\left[F_1\left(p\right),F_2\left(p\right)\right]}^{†}$$

(106)

where the vector of primary model parameters is defined as follows:

$$p\triangleq {\left[p_1,\mathrm{...},p_7\right]}^{†}\triangleq {\left[{\alpha }_T,\hspace{0.17em}\hspace{0.17em}{l}_p,\hspace{0.17em}\hspace{0.17em}{c}_p,\gamma ,\hspace{0.17em}\hspace{0.17em}{\sigma }_f,\hspace{0.17em}\hspace{0.17em}{N}_f,\hspace{0.17em}\hspace{0.17em}{\phi }_0\hspace{0.17em}\right]}^{†}$$

(107)

Note that the vector $p$ includes the initial condition ${\phi }_0$.

In terms of the “feature function” $F\left(p\right)\triangleq {\left[F_1\left(p\right),F_2\left(p\right)\right]}^{†}$, Eq. (17) can alternatively be written as follows:

$${\displaystyle \frac{dE\left(t\right)}{dt}}=-F_1\left(p\right)E^2\left(t\right)+F_2\left(p\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}E\left(0\right)=0$$

(108)

In terms of the feature function $F\left(p\right)\triangleq {\left[F_1\left(p\right),F_2\left(p\right)\right]}^{†}$, the solution of Eq. (108) has the following form:

$$E\left(t\right)={\left[{\displaystyle \frac{F_2\left(p\right)}{F_1\left(p\right)}}\right]}^{1/2}\tanh \left[tG\left(p\right)\right];\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}G\left(p\right)\triangleq \sqrt{F_1\left(p\right)F_2\left(p\right)}$$

(109)

Of course, a specific NODE would need to be constructed to model Eq. (108).

The form of Eq. (108) is suitable for applying the “n^th-Order Features Adjoint Sensitivity Analysis Methodology for Nonlinear Systems (n^th-FASAM-N)” [24], which is the most efficient methodology for computing sensitivities, particularly for sensitivities of second- and higher-order. This methodology considers the specific “features” of model parameters, such as the function $F\left(p\right)\triangleq {\left[F_1\left(p\right),F_2\left(p\right)\right]}^{†}$, to compute sensitivities with respect to model parameters more efficiently than by considering directly the respective primary parameters.

For the computation of 1^st-order sensitivities, the 1^st-FASAM-N commences by constructing the 1^st-Level Variational Sensitivity System (1^st-LVSS) for the variational function $\delta E\left(t\right)$ by applying the definition of the first-order G-differential to Eq. (108), which yields:

$${\displaystyle \frac{d}{d\epsilon }}\left\{{\displaystyle \frac{d\left[E^0\left(t\right)+\epsilon \delta E\left(t\right)\right]}{dt}}+\left[F_1^0+\epsilon \delta F_1\right]{\left[E^0+\epsilon \delta E\right]}^2-\left[F_2^0+\epsilon \left(\delta F_2\right)\right]\right\}\epsilon =0=0$$

(110)

$${\displaystyle \frac{d}{d\epsilon }}{\left\{\hspace{0.17em}{\left[E^0\left(t\right)+\epsilon \delta E\left(t\right)\right]}_{t=0}\right\}}_{\epsilon =0}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}=0$$

(111)

Performing the operations indicated in Eqs. (110) and (111) yields the following expression for the 1^st-LVSS satisfied by the variational function $\delta E\left(t\right)$:

$$\left[{\displaystyle \frac{d}{dt}}+2F_{1}E\left(t\right)\right]\delta E\left(t\right)=-\delta F_1{E}^2\left(t\right)+\delta F_2,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}t>0$$

(112)

$$\delta E\left(0\right)=0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}t=0.$$

(113)

The 1^st-LVSS represented by Eq. (112) is to be solved at the nominal values for the parameters and the state function $E\left(t\right)$ but the superscript “0” (which indicates “nominal values”) has been omitted to simplify the notation.

Numerically, the 1^st-LVSS would need to be solved anew for the various variations $\delta F_1$, $\delta F_2$, in the components of the feature function $F\left(p\right)$. This need for repeatedly solving the 1^st-LVSS can be avoided by constructing the corresponding 1^st-Level Adjoint Sensitivity System (1^st-LASS). The Hilbert space appropriate for the construction of the 1^st-LASS corresponding to Eq. (112) is endowed with the following particular form of Eq. (79):

$${〈u^{(a)}\left(t\right),u^{(b)}\left(t\right)〉}_1\triangleq \hspace{0.17em}{\displaystyle {\int }_0^{t_f}{u}^{(a)}\left(t\right)u^{(b)}\left(t\right)dt}$$

(114)

Using Eq. (114) to form the inner product of Eq. (112) with a yet undefined function ${\omega }^{(1)}\left(t\right)$ yields the following relation:

$${\displaystyle {\int }_0^{t_f}{\omega }^{(1)}\left(t\right)\left[\frac{d}{dt}+2F_{1}E\left(t\right)\right]\delta E\left(t\right)\hspace{0.17em}dt}=-\left(\delta F_1\right){\displaystyle {\int }_0^{t_f}{\omega }^{(1)}\left(t\right)E^2\left(t\right)\hspace{0.17em}dt}+\left(\delta F_2\right){\displaystyle {\int }_0^{t_f}{\omega }^{(1)}\left(t\right)\hspace{0.17em}dt}$$

(115)

Integrating by parts the left side of Eq. (115) yields the following relation:

$$\begin{array}{l}{\displaystyle {\int }_0^{t_f}{\omega }^{(1)}\left(t\right)\left[\frac{d}{dt}+2F_{1}E\left(t\right)\right]\delta E\left(t\right)\hspace{0.17em}dt}={\omega }^{(1)}\left(\tau \right)\delta E\left(\tau \right)\hspace{0.17em}-{\omega }^{(1)}\left(0\right)\delta E\left(0\right)\\ +{\displaystyle {\int }_0^{t_f}\delta E\left(t\right)\left[-\frac{d{\omega }^{(1)}\left(t\right)}{dt}+2F_{1}E\left(t\right){\omega }^{(1)}\left(t\right)\right]\hspace{0.17em}dt}.\end{array}$$

(116)

Identifying the integral on the right-side of Eq. (116) with the G-differential $\delta E\left(\tau \right)$ of the response $E\left(\tau \right)$ obtained in Eq. (32) and eliminating the unknown value $\delta E\left(\tau \right)\hspace{0.17em}$ from the right-side of Eq. (116) by setting ${\omega }^{(1)}\left(\tau \right)=0$ yields the following 1^st-Level Adjoint Sensitivity System (1^st-LASS) for the 1^st-level adjoint sensitivity function ${\omega }^{(1)}\left(t\right)$:

$$\left[-{\displaystyle \frac{d}{dt}}+2F_{1}E\left(t\right)\right]{\omega }^{(1)}\left(t\right)=\delta \left(t-t_f\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}t>0$$

(117)

$${\omega }^{(1)}\left(t_f\right)=0,\hspace{1em}t=t_f$$

(118)

The 1^st-LASS represented by Eqs. (117) and (118) is independent of variations in the feature functions (and/or parameters) so it would need to be solved only once, numerically. In the present case, the 1^st-LASS can be solved analytically to obtain the following closed-form expression for the 1^st-level adjoint sensitivity function ${\omega }^{(1)}\left(t\right)$:

$${\omega }^{(1)}\left(t\right)=H\left(t_f-t\right){\left\{{\displaystyle \frac{\cosh \left[tG\left(p\right)\right]}{\cosh \left[t_{f}G\left(p\right)\right]}}\right\}}^2$$

(119)

where $H\left(t-t_f\right)$ denotes the Heaviside functional.

Using Eqs. (116)‒(118) in Eq. (115) yields the following expression for the first-order total G-differential $\delta E\left(t_f\right)$ of the response $E\left(t_f\right)$ in terms of the 1^st-level adjoint function ${\omega }^{(1)}\left(t\right)$:

$$\delta E\left(t_f\right)=-\left(\delta F_1\right){\displaystyle {\int }_0^{t_f}{\omega }^{\left(1\right)}\left(t\right)}{E}^2\left(t\right)dt+\left(\delta F_2\right){\displaystyle {\int }_0^{t_f}{\omega }^{\left(1\right)}\left(t\right)dt}$$

(120)

It follows from Eqs. (120), (119) and (109) that the two sensitivities of the response $E\left(t_f\right)$ with respect to the two components of the feature function $F\triangleq {\left(F_1,F_2\right)}^{†}$ have the following expressions:

$${\displaystyle \frac{\partial E\left(t_f\right)}{\partial F_1}}=-{\displaystyle {\int }_0^{t_f}{\omega }^{\left(1\right)}\left(t\right)E^2\left(t\right)\hspace{0.17em}}dt={\displaystyle \frac{1}{2}}{\left[{\displaystyle \frac{F_2\left(p\right)}{F_1\left(p\right)}}\right]}^{1/2}\left\{{\displaystyle \frac{t_f}{{\cosh }^2\left[t_{f}G\left(p\right)\right]}}-{\displaystyle \frac{\tanh \left[t_{f}G\left(p\right)\right]}{G\left(p\right)}}\right\};$$

(121)

$${\displaystyle \frac{\partial E\left(t_f\right)}{\partial F_2}}={\displaystyle {\int }_0^{t_f}{\omega }^{\left(1\right)}\left(t\right)dt}={\displaystyle \frac{1}{2G\left(p\right)}}\tanh \left[t_{f}G\left(p\right)\right]+{\displaystyle \frac{t_f}{2{\cosh }^2\left[t_{f}G\left(p\right)\right]}}$$

(122)

The above expressions are to be evaluated at the nominal parameter values but the superscript “zero” has been omitted, for simplicity. The expressions obtained in Eqs. (121) and (122) can be verified by differentiating the expression provided in Eq. (109), evaluated at a user-chosen time $t=t_f$ within the interval $0<t_f<\infty$.

The sensitivities of the response $E\left(t_f\right)$ with respect to the model parameters and initial condition are obtained by using the following “chain-rule” relationship:

$${\displaystyle \frac{\partial E\left(t_f;F_1;F_2\right)}{\partial p_i}}={\displaystyle \frac{\partial E\left(t_f\right)}{\partial F_1}}{\displaystyle \frac{\partial F_1\left(p\right)}{\partial p_i}}+{\displaystyle \frac{\partial E\left(t_f\right)}{\partial F_2}}{\displaystyle \frac{\partial F_2\left(p\right)}{\partial p_i}};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},7.$$

(123)

The explicit expressions for the specific sensitivities of the response $E\left(t_f\right)$ with respect to the parameters underlying the feature functions are obtained using Eq. (123) in conjunction with Eqs. (121) and (122) while recalling the definitions of the feature functions $F_1\left(p\right)$ and $F_2\left(p\right)$ defined in Eq. (106). The detailed expressions of these sensitivities are as follows:

$${\displaystyle \frac{\partial E\left(t_f\right)}{\partial {\alpha }_T}}={\displaystyle \frac{\partial E\left(t_f\right)}{\partial F_1}}{\displaystyle \frac{\partial F_1}{\partial {\alpha }_T}}+{\displaystyle \frac{\partial E\left(t_f\right)}{\partial F_2}}{\displaystyle \frac{\partial F_2}{\partial {\alpha }_T}}={\displaystyle \frac{1}{2l_p{c}_p}}{\displaystyle \frac{\partial E\left(t_f\right)}{\partial F_1}}$$

(124)

$${\displaystyle \frac{\partial E\left(t_f\right)}{\partial l_p}}={\displaystyle \frac{\partial E\left(t_f\right)}{\partial F_1}}{\displaystyle \frac{\partial F_1}{\partial l_p}}+{\displaystyle \frac{\partial E\left(t_f\right)}{\partial F_2}}{\displaystyle \frac{\partial F_2}{\partial l_p}}=-{\displaystyle \frac{{\alpha }_T}{2{\left(l_p\right)}^2{c}_p}}{\displaystyle \frac{\partial E\left(t_f\right)}{\partial F_1}}$$

(125)

$${\displaystyle \frac{\partial E\left(t_f\right)}{\partial c_p}}={\displaystyle \frac{\partial E\left(t_f\right)}{\partial F_1}}{\displaystyle \frac{\partial F_1}{\partial c_p}}+{\displaystyle \frac{\partial E\left(t_f\right)}{\partial F_2}}{\displaystyle \frac{\partial F_2}{\partial c_p}}=-{\displaystyle \frac{{\alpha }_T}{2{\left(c_p\right)}^2{l}_p}}{\displaystyle \frac{\partial E\left(t_f\right)}{\partial F_1}}$$

(126)

$${\displaystyle \frac{\partial E\left(t_f\right)}{\partial \gamma }}={\displaystyle \frac{\partial E\left(t_f\right)}{\partial F_1}}{\displaystyle \frac{\partial F_1}{\partial \gamma }}+{\displaystyle \frac{\partial E\left(t_f\right)}{\partial F_2}}{\displaystyle \frac{\partial F_2}{\partial \gamma }}={\phi }_0{\sigma }_f{N}_f{\displaystyle \frac{\partial E\left(t_f\right)}{\partial F_2}}$$

(127)

$${\displaystyle \frac{\partial E\left(t_f\right)}{\partial {\sigma }_f}}={\displaystyle \frac{\partial E\left(t_f\right)}{\partial F_1}}{\displaystyle \frac{\partial F_1}{\partial {\sigma }_f}}+{\displaystyle \frac{\partial E\left(t_f\right)}{\partial F_2}}{\displaystyle \frac{\partial F_2}{\partial {\sigma }_f}}={\phi }_0\gamma N_f{\displaystyle \frac{\partial E\left(t_f\right)}{\partial F_2}}$$

(128)

$${\displaystyle \frac{\partial E\left(t_f\right)}{\partial N_f}}={\displaystyle \frac{\partial E\left(t_f\right)}{\partial F_1}}{\displaystyle \frac{\partial F_1}{\partial N_f}}+{\displaystyle \frac{\partial E\left(t_f\right)}{\partial F_2}}{\displaystyle \frac{\partial F_2}{\partial N_f}}={\phi }_0\gamma {\sigma }_f{\displaystyle \frac{\partial E\left(t_f\right)}{\partial F_2}}$$

(129)

$${\displaystyle \frac{\partial E\left(t_f\right)}{\partial {\phi }_0}}={\displaystyle \frac{\partial E\left(t_f\right)}{\partial F_1}}{\displaystyle \frac{\partial F_2}{\partial {\phi }_0}}+{\displaystyle \frac{\partial E\left(t_f\right)}{\partial F_2}}{\displaystyle \frac{\partial F_2}{\partial {\phi }_0}}=\gamma {\sigma }_f{N}_f{\displaystyle \frac{\partial E\left(t_f\right)}{\partial F_2}}$$

(130)

Notably, the application of the 1^st-FASAM-N requires one “large-scale” computation to solve the 1^st-LASS, cf. Eq. (117) and (118), which is a single ODE, to obtain the 1^st-level adjoint function ${\omega }^{(1)}\left(t\right)$, which is a scalar-valued function. However, solving the forward model, cf. Eq. (17), and the corresponding 1^st-LASS, comprising Eq. (117) and (118), would require the construction of a separate (albeit simpler) NODE. The 1^st-level adjoint function ${\omega }^{(1)}\left(t\right)$is subsequently used in performing two integrals (quadrature) for obtaining the two sensitivities of the response $E\left(t_f\right)$ with respect to the two components $F_1\left(p\right)$ and $F_2\left(p\right)$ of the feature function $F\left(p\right)\triangleq {\left(F_1,F_2\right)}^{†}$. Subsequently, all of the response sensitivities with respect to the model’s primary parameters are obtained analytically by using the chain-rule to differentiate the components of the feature function with respect to the underlying model parameters and initial conditions.

In contradistinction, if one wishes to compute directly the sensitivities of the response with respect to the model parameters and initial conditions, it has been shown in Subsections 5.1‒5.4 that the original NODE can be used to solve (backward in time) the 1^st-LASS, which comprises a system of three coupled ODEs (rather than a single ODE if the 1^st-FASAM is used) for obtaining the 1^st-level adjoint function, which is a vector-valued function comprising three components, cf. ${\chi }^{(1)}\left(t\right)\triangleq {\left[{\chi }_1^{(1)}\left(t\right),{\chi }_2^{(1)}\left(t\right),{\chi }_3^{(1)}\left(t\right)\right]}^{†}$ for the response $E\left(t_f\right)$. The respective vector-valued 1^st-level adjoint function is subsequently used in computing six (rather than two, if the 1^st-FASAM is used) integrals (quadrature) for obtaining the six sensitivities of the respective response with respect to the six model parameters.

Equations similar to Eq. (17) can be derived for the reactor-flux and reactor temperature responses, so the 1^st-FASAM can be applied in a similar fashion to compute the first-order sensitivities of these responses. Using the sensitivities with respect to the reactor temperature response would readily provide the first-order sensitivities of the reactor thermal conductivity response. However, corresponding to each of these responses, a specific NODE would need to be constructed. Of course, any of these specific NODE would have much simpler structures than the NODE for solving simultaneously the system of coupled ODEs presented in Subsections 5.1 through 5.4.

6. Discussion and Conclusions

This work has introduced the mathematical framework of the novel “First-Order Comprehensive Adjoint Sensitivity Analysis Methodology for Neural Ordinary Differential Equations (1^st-CASAM-NODE)” which yields exact expressions for the first-order sensitivities of NODE decoder-responses to the NODE parameters, including encoder initial conditions, while enabling the most efficient computation of these sensitivities. The application of the 1^st-CASAM-NODE has been illustrated by using the Nordheim-Fuchs reactor dynamics/safety phenomenological model, which is representative of physical systems that would be modeled by NODE while admitting exact analytical solutions for all quantities of interest (hidden states, decoder outputs, sensitivities with respect to all parameters and initial conditions, etc.). It has also been shown that if the equations underlying the physical model can be re-arranged so as to group the parameters/weights into functional “features” of several parameters, then the “First-Order Feature Adjoint Sensitivity Analysis Methodology for Nonlinear Systems (1^st-FASAM-N)” can be advantageously applied to compute the response sensitivities with respect to the feature functions (which are by definition fewer than the number of parameters). The response sensitivities with respect to the primary parameters are subsequently obtained analytically by using the chain-rule to differentiate the components of the feature function with respect to the underlying model parameters and initial conditions. Applying the 1^st-FASAM-N, however, would require the construction of a specific NODE for this purpose.

This work has also laid the foundation for the ongoing work on conceiving the “Second-Order Comprehensive Adjoint Sensitivity Analysis Methodology for Neural Ordinary Differential Equations (2^nd-CASAM-NODE)” which aims at yielding exact expressions for the second-order sensitivities of NODE decoder-responses to the NODE parameters and initial conditions while enabling the most efficient computation of these sensitivities.