Introducing the Second-Order Features Adjoint Sensitivity Analysis Methodology for Neural Ordinary Differential Equations. I: Mathematical Framework

1. Introduction

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

Neural ODEs are trained by minimizing a least-squares quadratic scalar-valued “loss function” which is usually meant to represent the discrepancy between a “reference solution” and the output produced by the NODE-decoder. The minimization procedure requires the computation of the gradients of the loss function with respect to the weights to be optimized. These gradients can be computed by using a first-order optimizer such as “stochastic gradient descent” [10,11], employing either the so-called “direct method” or the “adjoint method” [12,13,14]. The one-dimensional definite integrals, which appear when computing gradients via the “adjoint method” can be efficiently evaluated by using Gauss–Legendre quadrature, which has been shown [15] to be faster than ODE-based methods while retaining memory efficiency.

In the literature on neural nets, the “gradients of the loss function” are often called “sensitivities” and various aspects of the optimization/training procedure are occasionally called “sensitivity analysis.” However, this optimization procedure is not a bona-fide “sensitivity analysis,” since the “loss function” being minimized is of interest only for the “training” phase of the NODE-net, and the “sensitivities of the loss function” are driven towards the ideal zero-values by the minimization process while optimizing the NODE weights/parameters.

After the NODE-net is optimized to reproduce the underlying physical system as closely as possible, the subsequent responses of interest become various functionals of the NODE’s “decoder” output rather than some “loss function.” The physical system modeled by the NODE-net comprises parameters that stem from measurements and/or computations. Such parameters are not perfectly well known but are subject to 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-net. 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. Cacuci [16] has recently presented the “First-Order Comprehensive Adjoint Sensitivity Analysis Methodology for Neural Ordinary Differential Equations (1^st-CASAM-NODE)”, which is a pioneering sensitivity analysis 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.

It is well known that equations modeling real-life phenomena comprise not only scalar-valued model parameters but also functions of such scalar model parameters. For example, conservation equations modeling thermal-hydraulics phenomena depend on correlations, Reynolds, Nusselt numbers, etc., which are functions of various primary scalar parameters. Similarly, the neutron transport and/or diffusion equations in nuclear reactor engineering (e.g., reactor physics and shielding) involve group-averaged macroscopic cross section, which are scalar-valued functions of various sums of scalar-valued isotopic number densities and microscopic cross sections. It is convenient to refer to such scalar-valued functions as “features of primary model parameters.” Cacuci [17] has recently introduced the “n^th-Order Features Adjoint Sensitivity Analysis Methodology for Nonlinear Systems (n^th-FASAM-N),” which enables the most efficient computation of the exact expressions of arbitrarily high-order sensitivities of model responses with respect to the model’s “features.” Subsequently, the sensitivities of the responses with respect to the primary model parameters are determined, analytically and trivially, by applying the “chain-rule” to the expressions of the sensitivities with respect to the model’s “features/functions of parameters.”

The material presented in this work extends the 1^st-CASAM-NODE [16] methodology by using the concepts underlying the n^th-FASAM-N methodology [17] to introduce the newly developed “First-Order Features Adjoint Sensitivity Analysis Methodology for Neural Ordinary Differential Equations (1^st-FASAM-NODE)” and “Second-Order Features Adjoint Sensitivity Analysis Methodology for Neural Ordinary Differential Equations (2^nd-FASAM-NODE)”. The 1^st-FASAM-NODE methodology is introduced in Section 2 while the 2^nd-FASAM-NODE methodology is introduced in Section 3. The discussion presented in Section 4 concludes this work, noting that the 1^st-FASAM-NODE methodology enables the computation of exactly-determined first-order sensitivities of decoder response with respect to the NODE-parameters with unparalleled efficiency, while the 2^nd-FASAM-NODE methodology enables the computation of exactly-determined second-order sensitivities of decoder response with respect to the NODE-parameters with unparalleled efficiency. The efficiency of the 1^st-FASAM-NODE and the 2^nd-FASAM-NODE methodologies will be illustrated in the accompanying “Part II” [18] by means of heat and energy transfer responses in the Nordheim-Fuchs phenomenological model for reactor safety [19,20].

2. First-Order Features Adjoint Sensitivity Analysis Methodology for Neural Ordinary Differential Equations (1^st-FASAM-NODE)

The general representation of a NODE-network that comprises “features of primary model parameters/weights” is provided by the following system of “augmented” equations:

$${\displaystyle \frac{dh\left(t\right)}{dt}}=f\left[h\left(t\right);F\left(\theta \right);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$” will be used to denote “is defined as” or, equivalently, “is by definition equal to.”

(iii)

The$TH$-dimensional vector-valued nonlinear function $f\left[h\left(t\right);F\left(\theta \right);t\right]\triangleq {\left[f_1\left(h;F\left(\theta \right);t\right),\mathrm{...},f_{TH}\left(h;F\left(\theta \right);t\right)\right]}^{†}$ models the dynamics of the latent neurons. The components of the vector $\theta \triangleq {\left[{\theta }_1,\mathrm{...},{\theta }_{TW}\right]}^{†}$ represents learnable scalar adjustable weights, where $TW$ denotes the total number of adjustable weights in all of the latent neural nets. The components of the vector-valued function $F\left(\theta \right)\triangleq {\left[F_1\left(\theta \right),\mathrm{...},F_{TF}\left(\theta \right)\right]}^{†}$ represent the ”feature” functions of the respective weights, which are considered to be “primary parameters;” the quantity $TF$ denotes the “total number of feature/functions of the primary model parameters” comprised in the NODE. Evidently, the total number of feature functions must necessarily be smaller than the total number of primary parameters (weights), i.e., $TF<TW$. In the extreme case when there are no feature functions, it follows that $F_i\left(\theta \right)\equiv {\theta }_i$, for all $i=1,\mathrm{...},TW\equiv TF$.

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

(4)

After the “training” of the NODE has been accomplished (using the “adjoint” or other methods), the various “weights” will have been assigned “optimal” values which will have minimized the user-chosen loss functional, which is usually meant to represent the discrepancy between a “reference solution” and the output produced by the NODE-decoder. 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]}^{†}$. Using these optimal/nominal parameter values to solve the NODE-system will yield the optimal/nominal solution $h^0\left(t\right)$, which 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);F\left({\theta }^0\right);t\right],\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}t>0,$$

(5)

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

(6)

Furthermore, the optimal/nominal solution $h^0\left(t\right)$ will be used to obtain the components $r_n^0\left(h^0;{\phi }^0\right)$ of the vector of optimal/nominal decoder-responses; in view of Eq. (3), these components can be written in the following integral form:

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

(7)

However, the 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. 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$, while 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$. 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$. The forward state functions $h\left(t\right)$ are related to the weights and initial conditions through Eqs. (1)‒(3); consequently, 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-FASAM-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-FASAM-N [17] methodology. The fundamental concept for defining the sensitivity of an operator-valued quantity $R\left(e\right)$ with respect to variations $\delta e$ in a neighborhood around the nominal values $e^0$, has been shown by Cacuci [21] to be provided by the 1^st-order Gateaux- (G-) variation $\delta R\left(e^0;\delta e\right)$ of $R\left(e\right)$, which 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 }},$$

(8)

for a scalar $\epsilon \in F$ 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. (8) 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}$$

(9)

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]}^{†}\triangleq \delta h\left(t\right)\triangleq {\left[\delta h_1\left(t\right),\mathrm{...},\delta h_{TH}\left(t\right)\right]}^{†}$ and where the following definitions were used:

$${\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},$$

(10)

$${\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}.$$

(11)

The quantity ${\left\{\delta r_n\left(h^0;{\phi }^0;\delta \phi \right)\right\}}_{dir}$ in Eq. (10) 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 among the variations $v^{\left(1\right)}\left(t\right)$, $\delta x$, $\delta w$ and $\delta \theta$ are obtained from the first-order G-variations 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)};F\left({\theta }^0\right)+\epsilon \delta F;t\right]\right\}}_{\epsilon =0},$$

(12)

$${\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}.$$

(13)

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

$${\displaystyle \frac{d{v}^{\left(1\right)}\left(t\right)}{dt}}=J\left(h;F\right)v^{\left(1\right)}\left(t\right)+Q^{\left(1\right)}\left(h;F\right)\delta F,$$

(14)

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

(15)

where the various matrices and vectors are defined as follows:

$$\begin{array}{l}J\left(h;F\right)\triangleq {\left(\begin{array}{ccc}{J}_{11}\left(h;F\right)& ·& J_{1,TH}\left(h;F\right)\\ ·& ·& ·\\ J_{TH,1}\left(h;F\right)& ·& J_{TH,TH}\left(h;F\right)\end{array}\right)}_{TH\times TH}\triangleq {\displaystyle \frac{\partial f\left(h;F\right)}{\partial h}};\\ \hspace{0.17em}{J}_{ij}\left(h;F\right)\triangleq {\displaystyle \frac{\partial f_i}{\partial h_j}};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i,j=1,\mathrm{...},TH;\end{array}$$

(16)

$$\begin{array}{l}{Q}^{\left(1\right)}\left(h;F\right)\triangleq {\left(\begin{array}{ccc}{q}_{11}^{\left(1\right)}\left(h;F\right)& ·& q_{1,TF}^{\left(1\right)}\left(h;F\right)\\ ·& ·& ·\\ q_{TH,1}^{\left(1\right)}\left(h;F\right)& ·& q_{TF,TF}^{\left(1\right)}\left(h;F\right)\end{array}\right)}_{TH\times TF}\triangleq {\displaystyle \frac{\partial f\left(h;F\right)}{\partial F}};\\ q_{ij}^{\left(1\right)}\left(h;F\right)\triangleq {\displaystyle \frac{\partial f_i\left(h;F\right)}{\partial F_j}};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},TH;j=1,\mathrm{...},TF;\end{array}$$

(17)

$${\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};$$

(18)

$${\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}.$$

(19)

The system comprising Eqs. (14) and (15) 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. The 1^st-LVSS is to be satisfied at the nominal values $\left(h^0,F^0,x^0,w^0\right)$ for the respective functions and parameters, but this fact had not been indicated explicitly in order to simplify the notation. Note that the 1^st-LVSS would need to be solved anew for each component of the variations $\delta x$, $\delta w$ and $\delta F$, 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. (11) 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. (11) in terms of the solutions of the “1^st-Level Adjoint Sensitivity System (1^st-LASS),” to be constructed next, 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 will be defined on the domain ${\mathrm{\Omega }}_t\triangleq \left[t_0,t_f\right]$ and will be denoted as $H_1\left({\mathrm{\Omega }}_t\right)$, so that $v^{\left(1\right)}\left(t\right)\in H_1\left({\mathrm{\Omega }}_t\right)$. The inner product of two vectors ${\chi }^{(1)}\left(t\right)\triangleq {\left[{\chi }_1^{\left(1\right)}\left(t\right),\dots ,{\chi }_{TH}^{\left(1\right)}\left(t\right)\right]}^{†}\in H_1\left({\mathrm{\Omega }}_t\right)$ and ${\eta }^{(1)}\left(t\right)\triangleq {\left[{\eta }_1^{\left(1\right)}\left(t\right),\dots ,{\eta }_{TH}^{\left(1\right)}\left(t\right)\right]}^{†}\in H_1\left({\mathrm{\Omega }}_t\right)$ will be denoted as ${〈{\chi }^{(1)}\left(t\right),{\eta }^{(1)}\left(t\right)〉}_1$, and is defined as follows:

$${〈{\chi }^{(1)}\left(t\right),{\eta }^{(1)}\left(t\right)〉}_1\triangleq {\displaystyle {\int }_{t_0}^{t_f}{\left[{\chi }^{(1)}\left(t\right)\right]}^{†}{\eta }^{(1)}\left(t\right)dt}=\hspace{0.17em}{\displaystyle \sum _{i=1}^{TH}{\displaystyle {\int }_{t_0}^{t_f}{\chi }_i^{\left(1\right)}\left(t\right){\eta }_i^{\left(1\right)}\left(t\right)dt}}.$$

(20)

The scalar product ${〈{\chi }^{(1)}\left(t\right),{\eta }^{(1)}\left(t\right)〉}_1$ is required to hold in a neighborhood of the nominal values $\left(x^0;F^0;w^0;{\phi }^0\right)$.

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

$${〈a^{(1)}\left(t\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\displaystyle \frac{d{v}^{\left(1\right)}\left(t\right)}{dt}}-J\left(h;F\right)v^{\left(1\right)}\left(t\right)〉}_1={〈a^{(1)}\left(t\right),\hspace{0.17em}\hspace{0.17em}{Q}^{\left(1\right)}\left(h;F\right)\delta F〉}_1.$$

(21)

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

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

(22)

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

$${\displaystyle \frac{d{a}^{(1)}\left(t\right)}{dt}}=-{\left[J\left(h;F\right)\right]}^{†}{a}^{(1)}\left(t\right)-{\displaystyle \frac{\partial h_n^d\left[h\left(t\right);\phi \right]}{\partial h\left(t\right)}}\delta \left(t-t_f\right).$$

(23)

The definition of the 1^st-level adjoint sensitivity function $a^{(1)}\left(t\right)$ is now completed by requiring that it 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. (22), which is achieved by requiring that

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

(24)

It is evident from Eqs. (23) and (24) that to each response $h_n^d\left[h\left(t\right);\phi \right]$, there will correspond a distinct 1^st-level adjoint sensitivity function $a^{(1)}\left(n;t\right)$, $n=1,\mathrm{...},TR$. Since it is important to highlight this characteristic, the correspondence/dependence of the 1^st-level adjoint sensitivity function to/on the individual response under consideration will be made explicit by re-writing Eqs. (23) and (24) as follows:

$${\displaystyle \frac{d{a}^{(1)}\left(n;t\right)}{dt}}=-{\left[J\left(h;F\right)\right]}^{†}{a}^{(1)}\left(n;t\right)-{\displaystyle \frac{\partial h_n^d\left[h\left(t\right);\phi \right]}{\partial h\left(t\right)}}\delta \left(t-t_f\right);$$

(25)

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

(26)

The system of equations comprising Eqs. (25) and (26) constitute the “1^st-Level Adjoint Sensitivity System (1^st-LASS)” for the 1^st-level adjoint function $a^{(1)}\left(n;t\right)$, for the response $r_n\left(h;\phi \right)$. Evidently, the 1^st-LASS is linear in $a^{(1)}\left(n;t\right)$ and is independent of parameter variations, which implies that it needs to be solved just once to obtain the 1^st-level adjoint function $a^{(1)}\left(n;t\right)$. Notably, the left-side of the 1^st-LASS has the same as form as left-side of the “adjoint equations” used for training the NODE, but the “source” $\partial h_n^d\left[h\left(t\right);\phi \right]/\partial h\left(t\right)\delta \left(t-t_f\right)$ on the right-side of the 1^st-LASS stems from the “response” under consideration, while the “source” on the right-side of the “adjoint equations for training the NODE” stems from “loss functional” used in training the NODE. In component form, the 1^st-LASS has the following structure for $i=1,\mathrm{...},TH$:

$${\displaystyle \frac{d{a}_i^{\left(1\right)}\left(n;t\right)}{dt}}=-{\displaystyle \sum _{j=1}^{TH}\frac{\partial f_j}{\partial h_i}{a}_j^{\left(1\right)}\left(n;t\right)}-{\displaystyle \frac{\partial h_n^d\left[h\left(t\right);\phi \right]}{\partial h_i\left(t\right)}}\delta \left(t-t_f\right);$$

(27)

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

(28)

Using the results represented by Eqs. (23), (24), (21), and (11) in Eq. (22) 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(n;t\right)$:

$${\left\{\delta r_n\left(h^0;{\phi }^0;v^{\left(1\right)}\right)\right\}}_{ind}={〈a^{(1)}\left(n;t\right),\hspace{0.17em} Q^{\left(1\right)}\left(h;F\right)\delta F〉}_1+{\left[a^{(1)}\left(n;t_0\right)\right]}^{†} v^{\left(1\right)}\left(t_0\right).$$

(29)

Using in Eq. (29) the expression provided for $v^{\left(1\right)}\left(t_0\right)$ in Eq. (15) 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(n;t\right),\hspace{0.17em}\hspace{0.17em}{Q}^{\left(1\right)}\left(h;F\right)\delta F〉}_1\\ +{\left[a^{(1)}\left(n;t_0\right)\right]}^{†}{\displaystyle \frac{\partial h_e\left(x,w\right)}{\partial x}}\delta x+{\left[a^{(1)}\left(n;t_0\right)\right]}^{†}{\displaystyle \frac{\partial h_e\left(x,w\right)}{\partial w}}\delta w .\end{array}$$

(30)

Replacing the expression obtained in Eq. (29) for the “indirect-effect term” together with the expression of the direct-effect term provided by Eq. (10) into Eq. (9) yields the following expression (which is to be evaluated at the nominal values of all functions and parameters/weights) 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)={\displaystyle \frac{\partial h_n^d\left[h\left(t_f\right);\phi \right]}{\partial \phi }}\delta \phi +{\left[a^{(1)}\left(n;t_0\right)\right]}^{†}{\displaystyle \frac{\partial h_e\left(x,w\right)}{\partial w}}\delta w\\ +\hspace{0.17em}{\left[a^{(1)}\left(n;t_0\right)\right]}^{†}{\displaystyle \frac{\partial h_e\left(x,w\right)}{\partial x}}\delta x\hspace{0.17em}+{\displaystyle {\int }_{t_0}^{t_f}{\left[a^{(1)}\left(n;t\right)\right]}^{†}{Q}^{\left(1\right)}\left(h;F\right)\delta F}\hspace{0.17em}dt\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}n=1,\mathrm{...},TR.\end{array}$$

(31)

As indicated by the right-side of Eq. (31), the 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 F_i}}={\displaystyle {\int }_{t_0}^{t_f}{\left[a^{(1)}\left(n;t\right)\right]}^{†}\frac{\partial f\left(h;F\right)}{\partial F_i}}\hspace{0.17em}dt\hspace{0.17em}={\displaystyle \sum _{j=1}^{TH}{\displaystyle {\int }_{t_0}^{t_f}{a}_j^{\left(1\right)}\left(n;t\right)\frac{\partial f_j\left(h;F\right)}{\partial F_i}}\hspace{0.17em}dt};\hspace{0.17em}i=1,\mathrm{...},TF;\hspace{0.17em}n=1,\mathrm{...},TR;$$

(32)

$${\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}\hspace{0.17em}n=1,\mathrm{...},TR;$$

(33)

$${\displaystyle \frac{\partial r_n}{\partial w_i}}={\left[a^{(1)}\left(n;t_0\right)\right]}^{†}{\displaystyle \frac{\partial h_e\left(x,w\right)}{\partial w_i}}={\displaystyle \sum _{k=1}^{TH}{a}_k^{\left(1\right)}\left(n;t_0\right)}{\displaystyle \frac{\partial h_k^e\left(x,w\right)}{\partial w_i}};\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},TEW;\hspace{0.17em}\hspace{0.17em}n=1,\mathrm{...},TR;$$

(34)

$${\displaystyle \frac{\partial r_n}{\partial x_i}}={\left[a^{(1)}\left(n;t_0\right)\right]}^{†}{\displaystyle \frac{\partial h_e\left(x,w\right)}{\partial x_i}}={\displaystyle \sum _{k=1}^{TH}{a}_k^{\left(1\right)}\left(n;t_0\right)}{\displaystyle \frac{\partial h_k^e\left(x,w\right)}{\partial x_i}}\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},TI;\hspace{0.17em}\hspace{0.17em}n=1,\mathrm{...},TR;$$

(35)

The expressions of the first-order sensitivities obtained in Eqs. (32)‒(35) correspond to the n^th-response, $r_n$, where $n=1,\mathrm{...},TR$. The first-order sensitivities of the response with respect to the primary parameters/weights $\theta \triangleq {\left[{\theta }_1,\mathrm{...},{\theta }_{TW}\right]}^{†}$ are obtained analytically by applying the “chain-rule” to the expression in Eq. (32), which yields:

$${\displaystyle \frac{\partial r_n}{\partial {\theta }_j}}={\displaystyle \sum _{i=1}^{TF}\frac{\partial r_n}{\partial F_i}\frac{\partial F_i}{\partial {\theta }_j}}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}j=1,\mathrm{...},TW;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}n=1,\mathrm{...},TR.$$

(36)

When there are no feature functions, it follows that $F_i\left(\theta \right)\equiv {\theta }_i$, for all $i=1,\mathrm{...},TW\equiv TF$, so the expression obtained in Eq. (32) yields directly the first-order sensitivities $\partial r_n/\partial {\theta }_j$, for all $j=1,\mathrm{...},TW;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}n=1,\mathrm{...},TR$.

3. Second-Order Features Adjoint Sensitivity Analysis Methodology for Neural Ordinary Differential Equations (2^nd-FASAM-NODE)

The second-order sensitivities of the response $r_n\left(h;\phi \right)$ defined in Eq. (4) are obtained by computing the first-order sensitivities of the expressions obtained in Eqs. (32)‒(35), which represent the first-order sensitivities of the response $r_n\left(h;\phi \right)$ with respect to the various parameters (weights) and initial conditions (input). In other words, the second-order sensitivities of $r_n\left(h;\phi \right)$ will be computed by conceptually using their basic definitions as being the “first-order sensitivities of the first-order sensitivities.”

3.1. Second-Order Sensitivities Stemming from $\partial r_n/\partial F_i$, $i=1,\mathrm{...},TF$;$n=1,\mathrm{...},TR$

The second-order sensitivities stemming from the first-order sensitivities $\partial r_n/\partial F_i$ are obtained from the first-order G-differential of Eq. (32), for $i=1,\mathrm{...},TF$;$n=1,\mathrm{...},TR$, as follows:

$$\begin{array}{l}\delta \left({\displaystyle \frac{\partial r_n}{\partial F_i}}\right)\triangleq {\left\{{\displaystyle \frac{d}{d\epsilon }}\left[{\displaystyle {\int }_{t_0}^{t_f}{\left[a^{(1,0)}\left(n;t\right)+\epsilon \delta a^{(1)}\left(n;t\right)\right]}^{†}\frac{\partial }{\partial F_i}f\left[h^0+\epsilon v^{\left(1\right)};F\left({\theta }^0\right)+\epsilon \delta F;t\right]}\hspace{0.17em}dt\right]\right\}}_{\epsilon =0}\\ \triangleq \delta {\left({\displaystyle \frac{\partial r_n}{\partial F_i}}\right)}_{dir}+\delta {\left({\displaystyle \frac{\partial r_n}{\partial F_i}}\right)}_{ind}.\end{array}$$

(37)

In Eq. (37), the direct-effect term $\delta {\left(\partial r_n/\partial F_i\right)}_{dir}$ comprises the variation $\delta F$ (stemming from the model parameters) and is defined as follows:

$$\delta {\left({\displaystyle \frac{\partial r_n}{\partial F_i}}\right)}_{dir}\triangleq {\displaystyle {\int }_{t_0}^{t_f}{\left[a^{(1)}\left(n;t\right)\right]}^{†}\frac{{\partial }^{2}f\left(h;F\right)}{\partial F\partial F_i}\delta F}\hspace{0.17em}dt={\displaystyle \sum _{j=1}^{TH}{\displaystyle {\int }_{t_0}^{t_f}{a}_j^{\left(1\right)}\left(n;t\right){\displaystyle \sum _{k=1}^{TF}\frac{{\partial }^2{f}_j\left(h;F\right)}{\partial F_k\partial F_i}}}\hspace{0.17em}\delta F_k\hspace{0.17em}dt}\hspace{0.17em};$$

(38)

while the indirect-effect term $\delta {\left(\partial r_n/\partial F_i\right)}_{ind}$ comprises the variations $v^{\left(1\right)}\left(t\right)$and $\delta a^{(1)}\left(t\right)$, and is defined as follows:

$$\begin{array}{l}\delta {\left({\displaystyle \frac{\partial r_n}{\partial F_i}}\right)}_{ind}\triangleq {\displaystyle {\int }_{t_0}^{t_f}{\left[\delta a^{(1)}\left(n;t\right)\right]}^{†}\frac{\partial f\left(h;F\right)}{\partial F_i}}\hspace{0.17em}dt+{\displaystyle {\int }_{t_0}^{t_f}{\left[a^{(1)}\left(n;t\right)\right]}^{†}\left[\frac{{\partial }^{2}f\left(h;F\right)}{\partial h\partial F_i}{v}^{\left(1\right)}\left(t\right)\right]}\hspace{0.17em}dt\\ ={\displaystyle \sum _{j=1}^{TH}{\displaystyle {\int }_{t_0}^{t_f}\delta a_j^{\left(1\right)}\left(n;t\right)\frac{\partial f_j\left(h;F\right)}{\partial F_i}dt}\hspace{0.17em}}+{\displaystyle \sum _{j=1}^{TH}{\displaystyle \sum _{k=1}^{TF}{\displaystyle {\int }_{t_0}^{t_f}{a}_j^{\left(1\right)}\left(n;t\right)\frac{{\partial }^2{f}_j\left(h;F\right)}{\partial h_k\partial F_i}{v}_k^{\left(1\right)}\left(t\right)}}}\hspace{0.17em}dt;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},TF.\end{array}$$

(39)

The expressions in Eqs. (38) and (39) are to be evaluated at the nominal values of the respective functions and parameters, but the respective indication (namely: the superscript “zero”) has been omitted in order to simplify the notation.

The direct-effect term can be evaluated at this time for all variations $\delta F$, but the indirect-effect term can be evaluated only after having determined the variations $v^{\left(1\right)}\left(t\right)$ and $\delta a^{(1)}\left(n;t\right)$. The variation $v^{\left(1\right)}\left(t\right)$ is the solution of the 1^st-LVSS defined by Eqs. (14) and (15). On the other hand, the variational function $\delta a^{(1)}\left(n;t\right)$ is the solution of the system of equations obtained by G-differentiating the 1^st-LASS. This differentiation is performed most transparently by applying the definition of the G-differential to the 1^st-LASS presented in component form in Eqs. (27) and (28), to obtain the following system of equations for $j=1,\mathrm{...},TH$:

$$\begin{array}{l}{\left\{{\displaystyle \frac{d}{dt}}{\displaystyle \frac{d}{d\epsilon }}\left[a_j^{\left(1,0\right)}\left(n;t\right)+\epsilon \delta a_j^{\left(1\right)}\left(n;t\right)\right]\right\}}_{\epsilon =0}=-{\displaystyle \sum _{k=1}^{TH}\left\{\frac{d}{d\epsilon }\frac{\partial f_k\left(h^0+\epsilon v^{\left(1\right)};F^0+\epsilon \delta F\right)}{\partial h_j}\right.}\\ {\left.\times \left[a_k^{\left(1,0\right)}\left(n;t\right)+\epsilon \delta a_k^{\left(1\right)}\left(n;t\right)\right]\right\}}_{\epsilon =0}-{\left\{{\displaystyle \frac{d}{d\epsilon }}{\displaystyle \frac{\partial h_n^d\left[h^0+\epsilon v^{\left(1\right)};{\phi }^0+\epsilon \delta \phi \right]}{\partial h_j\left(t\right)}}\right\}}_{\epsilon =0}\delta \left(t-t_f\right);\end{array}$$

(40)

$${\left\{{\displaystyle \frac{d}{d\epsilon }}\left[a_j^{\left(1,0\right)}\left(n;t_f\right)+\epsilon \delta a_j^{\left(1\right)}\left(n;t_f\right)\right]\right\}}_{\epsilon =0}=0;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}j=1,\mathrm{...},TH.$$

(41)

Performing the operations indicated in Eqs. (40) and (41) yields the following relations for $j=1,\mathrm{...},TH$:

$$\begin{array}{l}{\displaystyle \frac{d}{dt}}\delta a_j^{\left(1\right)}\left(n\hspace{0.17em};t\right)=-{\displaystyle \sum _{k=1}^{TH}\frac{\partial f_k\left(h;F\right)}{\partial h_j}}\delta a_k^{\left(1\right)}\left(n\hspace{0.17em};t\right)-{\displaystyle \sum _{k=1}^{TH}{a}_k^{\left(1\right)}\left(n;t\right){\displaystyle \sum _{\mathcal{l}=1}^{TH}\frac{{\partial }^2{f}_k\left(h;F\right)}{\partial h_{\mathcal{l}}\partial h_j}{v}_{\mathcal{l}}^{\left(1\right)}\left(t\right)}}\\ -{\displaystyle \sum _{k=1}^{TH}{a}_k^{\left(1\right)}\left(n\hspace{0.17em};t\right){\displaystyle \sum _{\mathcal{l}=1}^{TF}\frac{{\partial }^2{f}_k\left(h;F\right)}{\partial F_{\mathcal{l}}\partial h_k}}}\delta F_{\mathcal{l}}-{\displaystyle \sum _{\mathcal{l}=1}^{TH}\frac{{\partial }^2{h}_n^d\left(h;\phi \right)}{\partial h_{\mathcal{l}}\partial h_j\left(t\right)}{v}_{\mathcal{l}}^{\left(1\right)}\left(t\right)}\hspace{0.17em}\delta \left(t-t_f\right)\\ -{\displaystyle \sum _{\mathcal{l}=1}^{TD}\frac{{\partial }^2{h}_n^d\left(h;\phi \right)}{\partial {\phi }_{\mathcal{l}}\partial h_j\left(t\right)}\delta {\phi }_{\mathcal{l}}}\delta \left(t-t_f\right);\end{array}$$

(42)

$$\delta a_j^{\left(1\right)}\left(n;t_f\right)=0.$$

(43)

The summations in Eq. (42) are rearranged as follows, for each $j=1,\mathrm{...},TH$, and $n=1,\mathrm{...},TR$:

$${\displaystyle \frac{d}{dt}}\delta a_j^{\left(1\right)}\left(n;t\right)=-{\displaystyle \sum _{\mathcal{l}=1}^{TH}{g}_{\mathcal{l}j}^n}{v}_{\mathcal{l}}^{\left(1\right)}\left(t\right)-{\displaystyle \sum _{k=1}^{TH}\frac{\partial f_k\left(h;F\right)}{\partial h_j}}\delta a_k^{\left(1\right)}\left(n;t\right)-p_j^{\left(2,n\right)}\hspace{0.17em},$$

(44)

where the following definitions were used:

$$\begin{array}{l}{g}_{\mathcal{l}j}^n={\displaystyle \sum _{k=1}^{TH}{a}_k^{\left(1\right)}\left(n;t\right)\frac{{\partial }^2{f}_k\left(h;F\right)}{\partial h_{\mathcal{l}}\partial h_j}}+{\displaystyle \frac{{\partial }^2{h}_n^d\left(h;\phi \right)}{\partial h_{\mathcal{l}}\partial h_j\left(t\right)}}\delta \left(t-t_f\right)\\ p_j^{\left(2,n\right)}\triangleq {\displaystyle \sum _{k=1}^{TH}{a}_k^{\left(1\right)}\left(n\hspace{0.17em};t\right){\displaystyle \sum _{\mathcal{l}=1}^{TF}\frac{{\partial }^2{f}_k\left(h;F\right)}{\partial F_{\mathcal{l}}\partial h_j}}}\delta F_{\mathcal{l}}+{\displaystyle \sum _{\mathcal{l}=1}^{TD}\frac{{\partial }^2{h}_n^d\left(h;\phi \right)}{\partial {\phi }_{\mathcal{l}}\partial h_j\left(t\right)}\delta {\phi }_{\mathcal{l}}}\delta \left(t-t_f\right);\end{array}$$

(45)

In matrix-format, Eq. (44) can be written as follows:

$${\displaystyle \frac{d}{dt}}\left[\delta a^{(1)}\left(n;t\right)\right]=-G^{\left(n\right)}{v}^{\left(1\right)}\left(t\right)-J^{†}\delta a^{(1)}\left(n;t\right)-p^{\left(2,n\right)}\hspace{0.17em},$$

(46)

where the following definitions were used:

$$G^{\left(n\right)}\triangleq \left[\begin{array}{ccc}{g}_{11}^n& ·& g_{1,TH}^n\\ ·& ·& ·\\ g_{TH,1}^n& ·& g_{TH,TH}^n\end{array}\right]\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{p}^{\left(2,n\right)}\triangleq {\left[p_1^{\left(2,n\right)},\mathrm{...},p_{TH}^{\left(2,n\right)}\right]}^{†}.$$

(47)

Concatenating Eqs. (14) and (46) yields the following 2^nd-Level Variational Sensitivity System (2^nd-LVSS) for the 2^nd-level variational function $v^{\left(2\right)}\left(n;t\right)\triangleq {\left[v^{\left(1\right)}\left(t\right),\delta a^{(1)}\left(n;t\right)\right]}^{†}$:

$${\displaystyle \frac{d{v}^{\left(2\right)}\left(n;t\right)}{dt}}=V^{\left(2,n\right)}{v}^{\left(2\right)}\left(n;t\right)+q^{\left(2,n\right)}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}n=1,\mathrm{...},TR,$$

(48)

where the superscript “2” indicates “second-level” and where the following definitions were used:

$$V^{\left(2,n\right)}\triangleq \left(\begin{array}{cc}J& 0\\ -\hspace{0.17em}{G}^{\left(n\right)}& -J^{†}\end{array}\right);\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{q}^{\left(2,n\right)}\triangleq {\left[Q^{\left(1\right)}\delta F;\hspace{0.17em}\hspace{0.17em}-p^{\left(2,n\right)}\right]}^{†}.$$

(49)

The matrix “$0$” in Eq. (49) denotes a $TH\times TH$ matrix having zero-elements. The initial-time and final-time conditions satisfied by $v^{\left(2\right)}\left(n;t\right)\triangleq {\left[v^{\left(1\right)}\left(t\right),\delta a^{(1)}\left(n;t\right)\right]}^{†}$ are provided in Eqs. (15) and (43), which can be written in vector form as follows:

$$v^{\left(1\right)}\left(t_0\right)={\displaystyle \frac{\partial h_e\left(x,w\right)}{\partial x}}\delta x+{\displaystyle \frac{\partial h_e\left(x,w\right)}{\partial w}}\delta w\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\delta a^{(1)}\left(n\hspace{0.17em};t_f\right)=0.$$

(50)

It is impractical to solve the 2^nd-LVSS repeatedly in order to compute each 2^nd-level variational function $v^{\left(2\right)}\left(n;t\right)\triangleq {\left[v^{\left(1\right)}\left(t\right),\delta a^{(1)}\left(n;t\right)\right]}^{†}$ that would be needed for every component of the variations $\delta F$ and $\delta \phi$. The need for computing $v^{\left(2\right)}\left(n;t\right)$ repeatedly can be circumvented by applying the principles of the 2^nd-FASAM-N methodology [17] which comprises the following sequence of steps:

• Consider that $v^{\left(2\right)}\left(n;t\right)\triangleq {\left[v^{\left(1\right)}\left(t\right),\delta a^{(1)}\left(n;t\right)\right]}^{†}$ is an element in a Hilbert space denoted as $H_2\left({\mathrm{\Omega }}_t\right)$, ${\mathrm{\Omega }}_t\triangleq \left[t_0,t_f\right]$, comprising as elements 2-block vectors having the following structure: ${\chi }^{\left(2\right)}\triangleq {\left[{\chi }_1^{\left(2\right)}\left(t\right),{\chi }_2^{\left(2\right)}\left(t\right)\right]}^{†}$, with ${\chi }_1^{\left(2\right)}\left(t\right)\triangleq {\left[{\chi }_{1,1}^{\left(2\right)}\left(t\right),\dots ,{\chi }_{1,TH}^{\left(2\right)}\left(t\right)\right]}^{†}$ and ${\chi }_2^{\left(2\right)}\left(t\right)\triangleq {\left[{\chi }_{2,1}^{\left(2\right)}\left(t\right),\dots ,{\chi }_{2,TH}^{\left(2\right)}\left(t\right)\right]}^{†}$. The Hilbert space $H_2\left({\mathrm{\Omega }}_t\right)$ is considered to be endowed with an inner product denoted as ${〈{\chi }^{\left(2\right)},{\eta }^{\left(2\right)}〉}_2$, between two vectors ${\chi }^{\left(2\right)}\triangleq {\left[{\chi }_1^{\left(2\right)}\left(t\right),{\chi }_2^{\left(2\right)}\left(t\right)\right]}^{†}\in H_2\left({\mathrm{\Omega }}_t\right)$ and ${\eta }^{\left(2\right)}={\left[{\eta }_1^{\left(2\right)}\left(t\right),{\eta }_2^{\left(2\right)}\left(t\right)\right]}^{†}\in H_2\left({\mathrm{\Omega }}_t\right)$, with ${\eta }_1^{\left(2\right)}\left(t\right)\triangleq {\left[{\eta }_{1,1}^{\left(2\right)}\left(t\right),\dots ,{\eta }_{1,TH}^{\left(2\right)}\left(t\right)\right]}^{†}$, ${\eta }_2^{\left(2\right)}\left(t\right)\triangleq {\left[{\eta }_{2,1}^{\left(2\right)}\left(t\right),\dots ,{\eta }_{2,TH}^{\left(2\right)}\left(t\right)\right]}^{†}$, which is defined as follows:

  $${〈{\chi }^{\left(2\right)},{\eta }^{\left(2\right)}〉}_2\triangleq {\displaystyle {\int }_{t_0}^{t_f}{\left[{\chi }^{\left(2\right)}\left(t\right)\right]}^{†}{\eta }^{\left(2\right)}\left(t\right)dt}={\displaystyle \sum _{j=1}^{TH}{\displaystyle {\int }_{t_0}^{t_f}{\chi }_{1,j}^{\left(2\right)}\left(t\right){\eta }_{1,j}^{\left(2\right)}\left(t\right)dt}}+{\displaystyle \sum _{j=1}^{TH}{\displaystyle {\int }_{t_0}^{t_f}{\chi }_{2,j}^{\left(2\right)}\left(t\right){\eta }_{2,j}^{\left(2\right)}\left(t\right)dt}}.$$

  (51)

The scalar product ${〈{\chi }^{\left(2\right)},{\eta }^{\left(2\right)}〉}_2$ is required to hold in a neighborhood of the nominal values $\left(h^0;x^0;F^0;w^0;{\phi }^0\right)$.

2.

Use the definition of the inner product provided in Eq. (51) to form the inner product of Eq. (48) with a vector $a^{\left(2\right)}\left(n;t\right)\triangleq {\left[a_1^{\left(2\right)}\left(n;t\right),a_2^{\left(2\right)}\left(n;t\right)\right]}^{†}\in H_2\left({\mathrm{\Omega }}_t\right)$, with $a_1^{\left(2\right)}\left(n;t\right)\triangleq {\left[a_{1,1}^{\left(2\right)}\left(n;t\right),\dots ,a_{1,TH}^{\left(2\right)}\left(n;t\right)\right]}^{†}$ and $a_2^{\left(2\right)}\left(n;t\right)\triangleq {\left[a_{2,1}^{\left(2\right)}\left(n;t\right),\dots ,a_{2,TH}^{\left(2\right)}\left(n;t\right)\right]}^{†}$, where the superscript “(2)” indicates “2^nd-level”, to obtain the following relationship:

$${〈a^{(2)}\left(n;t\right),\hspace{0.17em}{\displaystyle \frac{d{v}^{\left(2,n\right)}\left(t\right)}{dt}}-V^{\left(2,n\right)}{v}^{\left(2,n\right)}\left(t\right)〉}_2={〈a^{(2)}\left(n;t\right),\hspace{0.17em}{q}^{\left(2,n\right)}〉}_2.$$

(52)

3.

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

$$\begin{array}{l}{\displaystyle {\int }_{t_0}^{t_f}{\left[a_1^{\left(2\right)}\left(n;t\right)\right]}^{†}\frac{d{v}^{\left(1\right)}\left(t\right)}{dt}dt}-{\displaystyle {\int }_{t_0}^{t_f}{\left[a_1^{\left(2\right)}\left(n;t\right)\right]}^{†}J{v}^{\left(1\right)}\left(t\right)dt}+{\displaystyle {\int }_{t_0}^{t_f}{\left[a_2^{\left(2\right)}\left(n;t\right)\right]}^{†}\frac{d}{dt}\left[\delta a^{(1)}\left(t\right)\right]dt}\\ +{\displaystyle {\int }_{t_0}^{t_f}{\left[a_2^{\left(2\right)}\left(n;t\right)\right]}^{†}{G}^{\left(n\right)}{v}^{\left(1\right)}\left(t\right)dt}+{\displaystyle {\int }_{t_0}^{t_f}{\left[a_2^{\left(2\right)}\left(n;t\right)\right]}^{†}{J}^{†}\delta a^{(1)}\left(t\right)dt}={\left[a_1^{\left(2\right)}\left(n;t_f\right)\right]}^{†} v^{\left(1\right)}\left(t_f\right)\\ -{\left[a_1^{\left(2\right)}\left(n;t_0\right)\right]}^{†} v^{\left(1\right)}\left(t_0\right)+{\displaystyle {\int }_{t_0}^{t_f}{\left[v^{\left(1\right)}\left(t\right)\right]}^{†}\left\{-\frac{d{a}_1^{\left(2\right)}\left(n;t\right)}{dt}-J^{†}{a}_1^{\left(2\right)}\left(n;t\right)+{\left[G^{\left(n\right)}\right]}^{†}{a}_2^{\left(2\right)}\left(n;t\right)\right\} dt}\\ +{\left[a_2^{\left(2\right)}\left(n;t_f\right)\right]}^{†} \delta a^{(1)}\left(n;t_f\right)-{\left[a_2^{\left(2\right)}\left(n;t_0\right)\right]}^{†} \delta a^{(1)}\left(n;t_0\right)\\ +{\displaystyle {\int }_{t_0}^{t_f}{\left[\delta a^{(1)}\left(n;t\right)\right]}^{†}\left\{-\frac{d{a}_2^{\left(2\right)}\left(n;t\right)}{dt}+J{a}_2^{\left(2\right)}\left(n;t\right)\right\} dt}.\end{array}$$

(53)

4.

The two integral terms on the right-side of Eq. (53) are now required to represent the “indirect-effect” term defined in Eq. (39), which is achieved by imposing the following requirements:

$$\begin{array}{l}-{\displaystyle \frac{d}{dt}}{a}_1^{\left(2\right)}\left(i;n;t\right)-J^{†}{a}_1^{\left(2\right)}\left(i;n;t\right)+{\left[G^{\left(n\right)}\right]}^{†}{a}_2^{\left(2\right)}\left(i;n;t\right)={\left[{\displaystyle \frac{{\partial }^{2}f\left(h;F\right)}{\partial h\partial F_i}}\right]}^{†}{a}^{(1)}\left(n;t\right);\\ \hspace{0.17em}i=1,\mathrm{...},TF;\end{array}$$

(54)

$$-{\displaystyle \frac{d}{dt}}{a}_2^{\left(2\right)}\left(i;n;t\right)+J{a}_2^{\left(2\right)}\left(i;n;t\right)={\displaystyle \frac{\partial f\left(h;F\right)}{\partial F_i}};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},TF.$$

(55)

Note that the right-sides of Eqs. (54) and (55) also depend on the index $i=1,\mathrm{...},TF$, which means that a distinct 2^nd-level adjoint sensitivity function will correspond to each value of the index $i=1,\mathrm{...},TF$. This fact has been explicitly indicated by including the index $i=1,\mathrm{...},TF$ in the list of arguments of the components of the 2^nd-level adjoint sensitivity function $a^{\left(2\right)}\left(i;n;t\right)\triangleq {\left[a_1^{\left(2\right)}\left(i;n;t\right),a_2^{\left(2\right)}\left(i;n;t\right)\right]}^{†}$.

5.

The definition of the 2^nd-level adjoint sensitivity function $a^{\left(2\right)}\left(i;n;t\right)\triangleq {\left[a_1^{\left(2\right)}\left(i;n;t\right),a_2^{\left(2\right)}\left(i;n;t\right)\right]}^{†}$ is now completed by requiring it to satisfy the following boundary conditions, which eliminates the respective unknown terms on the right-side of Eq. (53)

$$a_1^{\left(2\right)}\left(i;n;t_f\right)=0\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{a}_2^{\left(2\right)}\left(i;n;t_0\right)=0\hspace{0.17em}.$$

(56)

The system of equations comprising Eqs. (54)‒(56) constitutes the “2^nd-Level Adjoint Sensitivity System (2^nd-LASS)” for the 2^nd-level adjoint sensitivity function $a^{\left(2\right)}\left(i;n;t\right)\triangleq {\left[a_1^{\left(2\right)}\left(i;n;t\right),a_2^{\left(2\right)}\left(i;n;t\right)\right]}^{†}$. Evidently, the 2^nd-LASS is linear in $a^{\left(2\right)}\left(i;n;t\right)\triangleq {\left[a_1^{\left(2\right)}\left(i;n;t\right),a_2^{\left(2\right)}\left(i;n;t\right)\right]}^{†}$ and is independent of parameter variations. Notably, this system of equations does not need to be solved simultaneously, but can be solved sequentially, by first solving Eq. (55) subject to the initial condition given in Eq. (56) to determine the function $a_2^{\left(2\right)}\left(i;n;t\right)$, and subsequently using the function $a_2^{\left(2\right)}\left(i;n;t\right)$ in Eq. (54) to solve it subject to the “final-time” condition given in Eq. (56) to obtain the function $a_1^{\left(2\right)}\left(i;n;t\right)$. The 2^nd-LASS is to be solved using the nominal values of all functions and parameters/weights, but the superscript “zero” has been omitted for notational simplicity.

6.

Using the results obtained in Eqs. (52), (54)‒(56) in Eq. (53) yields the following alternative expression for the “indirect-effect” term $\delta {\left(\partial r_n/\partial F_i\right)}_{ind}$, which no longer involves the 2^nd-level variational sensitivity function $v^{\left(2,n\right)}\left(t\right)$ but involves the 2^nd-level adjoint sensitivity function $a^{\left(2\right)}\left(i;n;t\right)\triangleq {\left[a_1^{\left(2\right)}\left(i;n;t\right),a_2^{\left(2\right)}\left(i;n;t\right)\right]}^{†}$:

$$\begin{array}{l}\delta {\left({\displaystyle \frac{\partial r_n}{\partial F_i}}\right)}_{ind}={\left[a_1^{\left(2\right)}\left(i;n;t_0\right)\right]}^{†} v^{\left(1\right)}\left(t_0\right)+{〈a^{(2)}\left(i;n;t\right),\hspace{0.17em}{q}^{\left(2,n\right)}〉}_2\\ ={\left[a_1^{\left(2\right)}\left(i;n;t_0\right)\right]}^{†} v^{\left(1\right)}\left(t_0\right)+{\displaystyle {\int }_{t_0}^{t_f}{\left[a_1^{\left(2\right)}\left(i;n;t\right)\right]}^{†}{Q}^{\left(1\right)}\delta F}\hspace{0.17em}dt-{\displaystyle {\int }_{t_0}^{t_f}{\left[a_2^{\left(2\right)}\left(i;n;t\right)\right]}^{†}}\hspace{0.17em}{p}^{\left(2,n\right)}dt.\end{array}$$

(57)

7.

Using in Eq. (57) the expression provided for $v^{\left(1\right)}\left(t_0\right)$ in Eq. (15), and adding the resulting expression for the indirect-effect term $\delta {\left(\partial r_n/\partial F_i\right)}_{ind}$ to the expression for the direct-effect term $\delta {\left(\partial r_n/\partial F_i\right)}_{dir}$ provided in Eq. (38), yields the following expression for the total first-order G-differential $\delta \left(\partial r_n/\partial F_i\right)$, for each $n=1,\mathrm{...},TR$:

$$\begin{array}{l}\delta \left({\displaystyle \frac{\partial r_n}{\partial F_i}}\right)={\displaystyle {\int }_{t_0}^{t_f}{\left[a^{(1)}\left(i;n;t\right)\right]}^{†}\frac{{\partial }^{2}f\left(h;F\right)}{\partial F\partial F_i}\delta F\hspace{0.17em}}dt+{\left[a_1^{\left(2\right)}\left(i;n;t_0\right)\right]}^{†}{\displaystyle \frac{\partial h_e\left(x,w\right)}{\partial x}}\delta x\\ +{\left[a_1^{\left(2\right)}\left(i;n;t_0\right)\right]}^{†}{\displaystyle \frac{\partial h_e\left(x,w\right)}{\partial w}}\delta w +{\displaystyle {\int }_{t_0}^{t_f}{\left[a_1^{\left(2\right)}\left(i;n;t\right)\right]}^{†}{Q}^{\left(1\right)}\delta F}\hspace{0.17em}dt-{\displaystyle {\int }_{t_0}^{t_f}{\left[a_2^{\left(2\right)}\left(i;n;t\right)\right]}^{†}}\hspace{0.17em}{p}^{\left(2,n\right)}(t)dt\\ \triangleq {\displaystyle \frac{{\partial }^2{r}_n}{\partial F\partial F_i}}\delta F+{\displaystyle \frac{{\partial }^2{r}_n}{\partial x\partial F_i}}\delta x+{\displaystyle \frac{{\partial }^2{r}_n}{\partial w\partial F_i}}\delta w+{\displaystyle \frac{{\partial }^2{r}_n}{\partial \phi \partial F_i}}\delta \phi \\ ={\displaystyle \sum _{j=1}^{TF}\frac{{\partial }^2{r}_n}{\partial F_j\partial F_i}}\delta F_j+{\displaystyle \sum _{j=1}^{TI}\frac{{\partial }^2{r}_n}{\partial x_j\partial F_i}}\delta x_j+{\displaystyle \sum _{j=1}^{TEW}\frac{{\partial }^2{r}_n}{\partial w_j\partial F_i}}\delta w_j+{\displaystyle \sum _{j=1}^{TD}\frac{{\partial }^2{r}_n}{\partial {\phi }_j\partial F_i}}\delta {\phi }_j;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},TF.\end{array}$$

(58)

The expression shown in Eq. (58) is to be evaluated at the nominal values of all functions and parameters/weights but the superscript “zero” (which has been used to indicate this fact) has been omitted for notational simplicity.

Identifying in Eq. (58) the expressions that multiply the individual variations $\delta F$, $\delta x$, $\delta w$ and $\delta \phi$ yields the following expressions for the second-order sensitivities that stem from the first-order sensitivity $\partial r_n/\partial F_i$ for $n=1,\mathrm{...},TR$:

$$\begin{array}{l}{\displaystyle \frac{{\partial }^2{r}_n}{\partial F_j\partial F_i}}={\displaystyle \sum _{k=1}^{TH}{\displaystyle {\int }_{t_0}^{t_f}{a}_k^{\left(1\right)}\left(n;t\right)\frac{{\partial }^2{f}_k\left(h;F\right)}{\partial F_j\partial h_i}}dt}+{\displaystyle \sum _{k=1}^{TH}{\displaystyle {\int }_{t_0}^{t_f}{a}_{1,k}^{\left(2\right)}\left(i;n;t\right)\frac{\partial f_k\left(h;F\right)}{\partial F_j}dt}}\hspace{0.17em}\\ -{\displaystyle \sum _{m=1}^{TH}{\displaystyle {\int }_{t_0}^{t_f}{a}_{2,m}^{\left(2\right)}\left(i;n;t\right)}\hspace{0.17em}{\displaystyle \sum _{k=1}^{TH}{a}_k^{\left(1\right)}\left(n;t\right)\frac{{\partial }^2{f}_k\left(h;F\right)}{\partial F_j\partial h_m}}dt};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i,j=1,\mathrm{...},TF;\hspace{0.17em}\hspace{0.17em}n=1,\mathrm{...},TR;\end{array}$$

(59)

$$\begin{array}{l}{\displaystyle \frac{{\partial }^2{r}_n}{\partial x_j\partial F_i}}={\left[a_1^{\left(2\right)}\left(i;n;t_0\right)\right]}^{†}{\displaystyle \frac{\partial h_e\left(x,w\right)}{\partial x_j}}={\displaystyle \sum _{k=1}^{TH}{a}_{1,k}^{\left(2\right)}\left(i;n;t_0\right)}{\displaystyle \frac{\partial h_k^e\left(x,w\right)}{\partial x_j}}\hspace{0.17em};\\ \hspace{0.17em}\hspace{0.17em}i=1,\mathrm{....},TF;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}j=1,\mathrm{....},TI;\hspace{0.17em}\hspace{0.17em}n=1,\mathrm{...},TR;\end{array}$$

(60)

$$\begin{array}{l}{\displaystyle \frac{{\partial }^2{r}_n}{\partial w_j\partial F_i}}={\left[a_1^{\left(2\right)}\left(i;n;t_0\right)\right]}^{†}{\displaystyle \frac{\partial h_e\left(x,w\right)}{\partial w_j}} ={\displaystyle \sum _{k=1}^{TH}{a}_{1,k}^{\left(2\right)}\left(i;n;t_0\right)}{\displaystyle \frac{\partial h_k^e\left(x,w\right)}{\partial w_j}}\hspace{0.17em};\\ \hspace{0.17em}\hspace{0.17em}i=1,\mathrm{....},TF;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}j=1,\mathrm{....},TEW\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}n=1,\mathrm{...},TR;\end{array}$$

(61)

$$\begin{array}{l}{\displaystyle \frac{{\partial }^2{r}_n}{\partial {\phi }_j\partial F_i}}=-{\displaystyle {\int }_{t_0}^{t_f}{\left[a_2^{\left(2\right)}\left(i;n;t\right)\right]}^{†}}{\displaystyle \frac{{\partial }^2{h}_n^d\left(h;\phi \right)}{\partial {\phi }_j\partial h}}\delta \left(t-t_f\right)\hspace{0.17em}dt \\ =-{\displaystyle \sum _{k=1}^{TH}{a}_{1,k}^{\left(2\right)}\left(i;n;t_f\right)\frac{{\partial }^2{h}_n^d\left[h\left(t_f\right);\phi \right]}{\partial {\phi }_j\partial h_k\left(t_f\right)}};\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{....},TF;\hspace{0.17em}\hspace{0.17em}j=1,\mathrm{....},TD;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}n=1,\mathrm{...},TR.\end{array}$$

(62)

The second-order sensitivities of the responses with respect to the primary parameters/weights $\theta \triangleq {\left[{\theta }_1,\mathrm{...},{\theta }_{TW}\right]}^{†}$ are obtained analytically by using the results obtained in Eqs. (59)‒(62) in conjunctions with the “chain-rule”:

$${\displaystyle \frac{{\partial }^2{r}_n}{\partial {\theta }_k\partial {\theta }_j}}={\displaystyle \frac{\partial }{\partial {\theta }_k}}\left[{\displaystyle \sum _{i=1}^{TF}\frac{\partial r_n}{\partial F_i}\frac{\partial F_i}{\partial {\theta }_j}}\right]\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}j,k=1,\mathrm{...},TW;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}n=1,\mathrm{...},TR.$$

(63)

When there are no feature functions, the second-order sensitivities are obtained by setting $F_i\left(\theta \right)\equiv {\theta }_i$, for all $i=1,\mathrm{...},TW\equiv TF$ in Eqs. (59)‒(62).

3.2. Second-Order Sensitivities Stemming from $\partial r_n/\partial {\phi }_i$,$i=1,\mathrm{...},TD$;$n=1,\mathrm{...},TR$

The second-order sensitivities stemming from the first-order sensitivities $\partial r_n/\partial {\phi }_i$ are obtained from the first-order G-differential of Eq. (33), for $\hspace{0.17em}i=1,\mathrm{...},TD$;$n=1,\mathrm{...},TR$, as follows:

$$\begin{array}{l}\delta \left({\displaystyle \frac{\partial r_n}{\partial {\phi }_i}}\right)\triangleq {\left\{{\displaystyle \frac{\partial }{\partial {\phi }_i}}{\displaystyle \frac{d}{d\epsilon }}\left[{\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\right]\right\}}_{\epsilon =0}\\ \triangleq \delta {\left({\displaystyle \frac{\partial r_n}{\partial {\phi }_i}}\right)}_{dir}+\delta {\left({\displaystyle \frac{\partial r_n}{\partial {\phi }_i}}\right)}_{ind},\end{array}$$

(64)

where the direct-effect term $\delta {\left(\partial r_n/\partial {\phi }_i\right)}_{dir}$ comprises the variation $\delta \phi$ and is defined as follows:

$$\delta {\left({\displaystyle \frac{\partial r_n}{\partial {\phi }_i}}\right)}_{dir}\triangleq {\displaystyle {\int }_{t_0}^{t_f}\frac{{\partial }^2{h}_n^d\left(h;\phi \right)}{\partial \phi \partial {\phi }_i}\delta \phi }\hspace{0.17em}\delta \left(t-t_f\right)\hspace{0.17em}dt={\displaystyle \sum _{k=1}^{TH}{\displaystyle {\int }_{t_0}^{t_f}\frac{{\partial }^2{h}_n^d\left(h;\phi \right)}{\partial {\phi }_k\partial {\phi }_i}}\delta \left(t-t_f\right)\hspace{0.17em}\delta {\phi }_{k}dt}\hspace{0.17em},$$

(65)

and where the indirect-effect term $\delta {\left(\partial r_n/\partial {\phi }_i\right)}_{ind}$ comprises the variation $v^{\left(1\right)}\left(t\right)$and is defined as follows:

$$\delta {\left({\displaystyle \frac{\partial r_n}{\partial {\phi }_i}}\right)}_{ind}\triangleq {\displaystyle {\int }_{t_0}^{t_f}\frac{{\partial }^2{h}_n^d\left(h;\phi \right)}{\partial h\partial {\phi }_i}{v}^{\left(1\right)}\left(t\right)\delta \left(t-t_f\right)}\hspace{0.17em}dt={\displaystyle \sum _{k=1}^{TH}{\displaystyle {\int }_{t_0}^{t_f}\frac{{\partial }^2{h}_n^d\left(h;\phi \right)}{\partial h_k\partial {\phi }_i}{v}_k^{\left(1\right)}\left(t\right)\delta \left(t-t_f\right)dt}}.$$

(66)

The expressions shown in Eqs. (65) and (66) are to be evaluated at the nominal values of the respective functions and parameters, but the respective indication (namely: the superscript “zero”) has been omitted in order to simplify the notation.

The direct-effect term $\delta {\left(\partial r_n/\partial {\phi }_i\right)}_{dir}$ can be evaluated at this time for all variations $\delta \phi$, but the indirect-effect term $\delta {\left(\partial r_n/\partial {\phi }_i\right)}_{ind}$ can be evaluated only after having determined the variation $v^{\left(1\right)}\left(t\right)$, which is the solution of the 1^st-LVSS defined by Eqs. (14) and (15). The need for solving this 1^st-LVSS can be avoided if the indirect-effect term defined in Eq. (66) 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 applying the same concepts and steps as in Section 2, as follows:

• Using Eq. (20), form the inner product of Eq. (14) with a vector $b^{(2)}\left(t\right)\triangleq {\left[b_1^{\left(2\right)}\left(t\right),\dots ,b_{TH}^{\left(2\right)}\left(t\right)\right]}^{†}\in H_1\left({\mathrm{\Omega }}_t\right)$, where the superscript “(2)” indicates “2^nd-Level”, to obtain the following relationship:

  $${〈b^{(2)}\left(t\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\displaystyle \frac{d{v}^{\left(1\right)}\left(t\right)}{dt}}-J\left(h;F\right)v^{\left(1\right)}\left(t\right)〉}_1={〈b^{(2)}\left(t\right),\hspace{0.17em}\hspace{0.17em}{Q}^{\left(1\right)}\left(h;F\right)\delta F〉}_1.$$

  (67)
• Using the definition of the adjoint operator in $H_1\left({\mathrm{\Omega }}_t\right)$, the left-side of Eq. (21) is integrated by parts over the independent variable $t$ to obtain the following relation:

  $$\begin{array}{l}{\displaystyle {\int }_{t_0}^{t_f}{\left[b^{(2)}\left(t\right)\right]}^{†}\frac{d{v}^{\left(1\right)}\left(t\right)}{dt}dt}-{\displaystyle {\int }_{t_0}^{t_f}{\left[b^{(2)}\left(t\right)\right]}^{†}J\left(h;F\right)v^{\left(1\right)}\left(t\right)dt}={\left[b^{(2)}\left(t_f\right)\right]}^{†} v^{\left(1\right)}\left(t_f\right)\\ -{\left[b^{(2)}\left(t_0\right)\right]}^{†} v^{\left(1\right)}\left(t_0\right)+{\displaystyle {\int }_{t_0}^{t_f}{\left[v^{\left(1\right)}\left(t\right)\right]}^{†}\left\{-\frac{d{b}^{(2)}\left(t\right)}{dt}-{\left[J\left(h;F\right)\right]}^{†}{b}^{(2)}\left(t\right)\right\} dt}\hspace{0.17em}.\end{array}$$

  (68)
• The last term on the right-side of Eq. (68) is now required to represent the “indirect-effect” term defined in Eq. (66), which is achieved by requiring that the 2^nd-level adjoint sensitivity function $b^{(2)}\left(t\right)$ satisfy the following relation written in NODE-format:

  $${\displaystyle \frac{d}{dt}}{b}^{(2)}\left(i;n;t\right)=-{\left[J\left(h;F\right)\right]}^{†}{b}^{(2)}\left(i;n;t\right)-{\displaystyle \frac{{\partial }^2{h}_n^d\left(h;\phi \right)}{\partial h\partial {\phi }_i}}\delta \left(t-t_f\right).$$

  (69)

It is evident from Eq. (69) that to each response ${\partial }^2{h}_n^d\left(h;\phi \right)/\partial h\partial {\phi }_i$, there will correspond a distinct 2^nd-level adjoint sensitivity function $b^{(2)}\left(i;n;t\right)$, for each $i=1,\mathrm{...},TI$ and $n=1,\mathrm{...},TR$. This important fact has been highlighted by adding the indices $i=1,\mathrm{...},TI$ and $n=1,\mathrm{...},TR$ in the list of arguments of $b^{(2)}\left(i;n;t\right)$.

4.

The definition of the 2^nd-level adjoint sensitivity function $b^{(2)}\left(i;n;t\right)$ is now completed by requiring it to eliminate the term containing the unknown values $v^{\left(1\right)}\left(t_f\right)$ in Eq. (68), which is accomplished by requiring $b^{(2)}\left(i;n;t\right)$ to satisfy the following boundary condition at the final time $t=t_f$:

$$b^{(2)}\left(i;n;t_f\right)=0\hspace{0.17em}.$$

(70)

The system of equations comprising Eqs. (69) and (70) constitute the “2^nd-Level Adjoint Sensitivity System (2^nd-LASS)” for the function $b^{(2)}\left(i;n;t\right)$; this 2^nd-LASS is linear in $b^{(2)}\left(i;n;t\right)$ and is independent of parameter variations, which implies that it needs to be solved once to obtain the 2^nd-level adjoint sensitivity function $b^{(2)}\left(i;n;t\right)$, for each $i=1,\mathrm{...},TI$, $n=1,\mathrm{...},TR$. Notably, the left-side of this 2^nd-LASS has the same as form as left-side of the “adjoint equations” used for training the NODE, but the sources ${\partial }^2{h}_n^d\left(h;\phi \right)/\partial h\partial {\phi }_i$ on the right-side of this 2^nd-LASS stem from each of the “responses” under consideration, while the “source” on the right-side of the “adjoint equations for training the NODE” stems from the “loss functional” used in training the NODE.

5.

Using the results obtained in Eqs. (69), (70), and (67) in Eq. (68) 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 instead involves the 2^nd level adjoint function $b^{(2)}\left(i;n;t\right)$:

$$\begin{array}{l}\delta {\left(\partial r_n/\partial {\phi }_i\right)}_{ind}={〈b^{(2)}\left(i;n;t\right),\hspace{0.17em} {\displaystyle \frac{\partial f\left(h;F\right)}{\partial F}}\delta F〉}_1\\ +{\left[b^{(2)}\left(i;n;t_0\right)\right]}^{†}{\displaystyle \frac{\partial h_e\left(x,w\right)}{\partial x}}\delta x+{\left[b^{(2)}\left(i;n;t_0\right)\right]}^{†}{\displaystyle \frac{\partial h_e\left(x,w\right)}{\partial w}}\delta w .\end{array}$$

(71)

6.

Adding the expression obtained in Eq. (71) for the “indirect-effect term” together with the expression of the direct-effect term provided by Eq. (65) yields the following expression for the G-variation $\delta \left(\partial r_n/\partial {\phi }_i\right)$ defined in Eq. (64), which is to be evaluated at the nominal values of all functions and parameters/weights, for $i=1,\mathrm{...},TI;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}n=1,\mathrm{...},TR$:

$$\begin{array}{l}\delta \left({\displaystyle \frac{\partial r_n}{\partial {\phi }_i}}\right)={\displaystyle {\int }_{t_0}^{t_f}\frac{{\partial }^2{h}_n^d\left(h;\phi \right)}{\partial \phi \partial {\phi }_i}\delta \phi }\hspace{0.17em}\delta \left(t-t_f\right)\hspace{0.17em}dt+{\left[b^{(2)}\left(i;n;t_0\right)\right]}^{†}{\displaystyle \frac{\partial h_e\left(x,w\right)}{\partial x}}\delta x \\ +{\left[b^{(2)}\left(i;n;t_0\right)\right]}^{†}{\displaystyle \frac{\partial h_e\left(x,w\right)}{\partial w}}\delta w+{\displaystyle {\int }_{t_0}^{t_f}{b}^{(2)}\left(i;n;t\right)\frac{\partial f\left(h;F\right)}{\partial F}\delta F}\hspace{0.17em}dt\hspace{0.17em}\\ \triangleq {\displaystyle \frac{{\partial }^2{r}_n}{\partial \phi \partial {\phi }_i}}\delta \phi +{\displaystyle \frac{{\partial }^2{r}_n}{\partial x\partial {\phi }_i}}\delta x+{\displaystyle \frac{{\partial }^2{r}_n}{\partial w\partial {\phi }_i}}\delta w+{\displaystyle \frac{{\partial }^2{r}_n}{\partial F\partial {\phi }_i}}\delta F\\ ={\displaystyle \sum _{j=1}^{TF}\frac{{\partial }^2{r}_n}{\partial F_j\partial {\phi }_i}}\delta F_j+{\displaystyle \sum _{j=1}^{TI}\frac{{\partial }^2{r}_n}{\partial x_j\partial {\phi }_i}}\delta x_j+{\displaystyle \sum _{j=1}^{TEW}\frac{{\partial }^2{r}_n}{\partial w_j\partial {\phi }_i}}\delta w_j+{\displaystyle \sum _{j=1}^{TD}\frac{{\partial }^2{r}_n}{\partial {\phi }_j\partial {\phi }_i}}\delta {\phi }_j;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},TD.\end{array}$$

(72)

As indicated by the rightmost-side of Eq. (72), the 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:

$$\begin{array}{l}{\displaystyle \frac{{\partial }^2{r}_n}{\partial F_j\partial {\phi }_i}}={\displaystyle {\int }_{t_0}^{t_f}{\left[b^{(2)}\left(i;n;t\right)\right]}^{†}\frac{\partial f\left(h;F\right)}{\partial F_j}}\hspace{0.17em}dt={\displaystyle \sum _{k=1}^{TH}{\displaystyle {\int }_{t_0}^{t_f}{b}_k^{\left(2\right)}\left(i;n;t\right)\frac{\partial f_k\left(h;F\right)}{\partial F_j}}\hspace{0.17em}dt};\hspace{0.17em}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}j=1,\mathrm{...},TF;\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;\end{array}$$

(73)

$$\begin{array}{l}{\displaystyle \frac{{\partial }^2{r}_n}{\partial x_j\partial {\phi }_i}}={\left[b^{(2)}\left(i;n;t_0\right)\right]}^{†}{\displaystyle \frac{\partial h_e\left(x,w\right)}{\partial x_j}}={\displaystyle \sum _{k=1}^{TH}{b}_k^{\left(2\right)}\left(i;n;t_0\right)}{\displaystyle \frac{\partial h_k^e\left(x,w\right)}{\partial x_j}};\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}j=1,\mathrm{...},TI;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},TD;\hspace{0.17em}\hspace{0.17em}n=1,\mathrm{...},TR.\end{array}$$

(74)

$$\begin{array}{l}{\displaystyle \frac{{\partial }^2{r}_n}{\partial w_j\partial {\phi }_i}}={\left[b^{(2)}\left(i;n;t_0\right)\right]}^{†}{\displaystyle \frac{\partial h_e\left(x,w\right)}{\partial w_j}}={\displaystyle \sum _{k=1}^{TH}{b}_k^{\left(2\right)}\left(i;n;t_0\right)}{\displaystyle \frac{\partial h_k^e\left(x,w\right)}{\partial w_j}};\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}j=1,\mathrm{...},TEW;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},TD;\hspace{0.17em}\hspace{0.17em}n=1,\mathrm{...},TR;\end{array}$$

(75)

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

(76)

The first-order sensitivities of the response with respect to the primary parameters/weights $\theta \triangleq {\left[{\theta }_1,\mathrm{...},{\theta }_{TW}\right]}^{†}$ are obtained analytically by applying the “chain-rule” $\partial \left[\right]/\partial {\theta }_m={\displaystyle \sum _{\mu =1}^{TF}\left(\partial \left[\right]/\partial F_{\mu }\right)\left(\partial F_{\mu }/\partial {\theta }_m\right)}$ to the expression in Eq. (73). When there are no feature functions, then $F_i\left(\theta \right)\equiv {\theta }_i$, for all $i=1,\mathrm{...},TW\equiv TF$.

3.3. Second-Order Sensitivities Stemming from $\partial r_n/\partial w_i$, $i=1,\mathrm{...},TEW$;$n=1,\mathrm{...},TR$

The second-order sensitivities stemming from the first-order sensitivities $\partial r_n/\partial w_i$ are obtained from the first-order G-differential of Eq. (34), for $\hspace{0.17em}i=1,\mathrm{...},TEW$;$n=1,\mathrm{...},TR$, as follows:

$$\begin{array}{l}\delta \left({\displaystyle \frac{\partial r_n}{\partial w_i}}\right)\triangleq {\left\{{\displaystyle \frac{d}{d\epsilon }}\left[{\left[a^{(1,0)}\left(n;t_0\right)+\epsilon \delta a^{(1)}\left(n;t_0\right)\right]}^{†}{\displaystyle \frac{\partial h_e\left(x^0+\epsilon \delta x,w^0+\epsilon \delta w\right)}{\partial w_i}}\right]\right\}}_{\epsilon =0}\\ \triangleq \delta {\left({\displaystyle \frac{\partial r_n}{\partial w_i}}\right)}_{dir}+\delta {\left({\displaystyle \frac{\partial r_n}{\partial w_i}}\right)}_{ind},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},TEW;\hspace{0.17em}\hspace{0.17em}n=1,\mathrm{...},TR,\end{array}$$

(77)

where the direct-effect term $\delta {\left(\partial r_n/\partial w_i\right)}_{dir}$ comprises the variations $\delta x$ and $\delta w$, defined as follows:

$$\delta {\left({\displaystyle \frac{\partial r_n}{\partial w_i}}\right)}_{dir}\triangleq {\left[a^{(1)}\left(n;t_0\right)\right]}^{†}\left[{\displaystyle \frac{{\partial }^2{h}_e\left(x,w\right)}{\partial x\partial w_i}}\delta x+{\displaystyle \frac{{\partial }^2{h}_e\left(x,w\right)}{\partial w\partial w_i}}\delta w\right];$$

(78)

and where the indirect-effect term $\delta {\left(\partial r_n/\partial w_i\right)}_{ind}$ comprises the variation $\delta a^{(1)}\left(n;t\right)$ and is defined as follows:

$$\delta {\left({\displaystyle \frac{\partial r_n}{\partial w_i}}\right)}_{ind}\triangleq {\displaystyle {\int }_{t_0}^{t_f}{\left[\delta a^{(1)}\left(n;t\right)\right]}^{†}\frac{\partial h_e\left(x,w\right)}{\partial w_i}\delta \left(t-t_0\right)}\hspace{0.17em}dt.$$

(79)

The expressions shown in Eqs. (78) and (79) are to be evaluated at the nominal values of the respective functions and parameters, but the respective indication (namely: the superscript “zero”) has been omitted in order to simplify the notation.

The direct-effect term $\delta {\left(\partial r_n/\partial w_i\right)}_{dir}$ can be evaluated at this time for all variations $\delta x$ and $\delta w$, but the indirect-effect term $\delta {\left(\partial r_n/\partial w_i\right)}_{ind}$ can be evaluated only after having determined the variation $\delta a^{(1)}\left(n;t\right)$, which is the solution of the 2^nd-LVSS defined by Eqs. (48) and (50). The need for solving this 1^st-LVSS can be avoided if the indirect-effect term defined in Eq. (79) could be expressed in terms of a “right-hand side” that does not involve the function $\delta a^{(1)}\left(n;t\right)$. This goal can be achieved by applying the same concepts and steps as in Section 3.1, as follows:

• Use the definition of the inner product provided in Eq. (51) to form the inner product of Eq. (48) with a vector $c^{\left(2\right)}\left(n;t\right)\triangleq {\left[c_1^{\left(2\right)}\left(n;t\right),c_2^{\left(2\right)}\left(n;t\right)\right]}^{†}\in H_2\left({\mathrm{\Omega }}_t\right)$, with $c_1^{\left(2\right)}\left(n;t\right)\triangleq {\left[c_{1,1}^{\left(2\right)}\left(n;t\right),\dots ,c_{1,TH}^{\left(2\right)}\left(n;t\right)\right]}^{†}$ and $c_2^{\left(2\right)}\left(n;t\right)\triangleq {\left[c_{2,1}^{\left(2\right)}\left(n;t\right),\dots ,c_{2,TH}^{\left(2\right)}\left(n;t\right)\right]}^{†}$, where the superscript “(2)” indicates “2^nd-level”, to obtain the following relationship:

  $${〈c^{(2)}\left(n;t\right),\hspace{0.17em}{\displaystyle \frac{d{v}^{\left(2,n\right)}\left(t\right)}{dt}}-V^{\left(2,n\right)}{v}^{\left(2,n\right)}\left(t\right)〉}_2={〈c^{(2)}\left(n;t\right),\hspace{0.17em}{q}^{\left(2,n\right)}〉}_2.$$

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

  $$\begin{array}{l}{\displaystyle {\int }_{t_0}^{t_f}{\left[c_1^{\left(2\right)}\left(n;t\right)\right]}^{†}\frac{d{v}^{\left(1\right)}\left(t\right)}{dt}dt}-{\displaystyle {\int }_{t_0}^{t_f}{\left[c_1^{\left(2\right)}\left(n;t\right)\right]}^{†}J{v}^{\left(1\right)}\left(t\right)dt}\\ +{\displaystyle {\int }_{t_0}^{t_f}{\left[c_2^{\left(2\right)}\left(n;t\right)\right]}^{†}\frac{d}{dt}\left[\delta a^{(1)}\left(n;t\right)\right]dt}+{\displaystyle {\int }_{t_0}^{t_f}{\left[c_2^{\left(2\right)}\left(n;t\right)\right]}^{†}{G}^{\left(n\right)}{v}^{\left(1\right)}\left(t\right)dt}\\ +{\displaystyle {\int }_{t_0}^{t_f}{\left[c_2^{\left(2\right)}\left(n;t\right)\right]}^{†}{J}^{†}\delta a^{(1)}\left(n;t\right)dt}={\left[c_1^{\left(2\right)}\left(n;t_f\right)\right]}^{†} v^{\left(1\right)}\left(t_f\right)\\ -{\left[c_1^{\left(2\right)}\left(n;t_0\right)\right]}^{†} v^{\left(1\right)}\left(t_0\right)+{\left[c_2^{\left(2\right)}\left(n;t_f\right)\right]}^{†} \delta a^{(1)}\left(n;t_f\right)-{\left[c_2^{\left(2\right)}\left(n;t_0\right)\right]}^{†} \delta a^{(1)}\left(n;t_0\right)\\ +{\displaystyle {\int }_{t_0}^{t_f}{\left[v^{\left(1\right)}\left(t\right)\right]}^{†}\left\{-\frac{d{c}_1^{\left(2\right)}\left(n;t\right)}{dt}-J^{†}{c}_1^{\left(2\right)}\left(n;t\right)+{\left[G^{\left(n\right)}\right]}^{†}{c}_2^{\left(2\right)}\left(n;t\right)\right\} dt}\\ +{\displaystyle {\int }_{t_0}^{t_f}{\left[\delta a^{(1)}\left(n;t\right)\right]}^{†}\left\{-\frac{d{c}_2^{\left(2\right)}\left(n;t\right)}{dt}+J{c}_2^{\left(2\right)}\left(n;t\right)\right\} dt}.\end{array}$$

  (81)
• The two integral terms on the right-side of Eq. (81) are now required to represent the “indirect-effect” term defined in Eq. (79), which is achieved by imposing the following requirements on the components of the 2^nd-level adjoint sensitivity function $c^{\left(2\right)}\left(i;n;t\right)\triangleq {\left[c_1^{\left(2\right)}\left(i;n;t\right),c_2^{\left(2\right)}\left(i;n;t\right)\right]}^{†}$:

  $$-{\displaystyle \frac{d}{dt}}{c}_1^{\left(2\right)}\left(i;n;t\right)-J^{†}{c}_1^{\left(2\right)}\left(i;n;t\right)+{\left[G^{\left(n\right)}\right]}^{†}{c}_2^{\left(2\right)}\left(i;n;t\right)=0;$$

  (82)

  $$-{\displaystyle \frac{d}{dt}}{c}_2^{\left(2\right)}\left(i;n;t\right)+J{c}_2^{\left(2\right)}\left(i;n;t\right)={\displaystyle \frac{\partial h_e\left(x,w\right)}{\partial w_i}}\delta \left(t-t_0\right).$$

  (83)

Since the right-side of Eq. (83) also depends on the index $i=1,\mathrm{...},TF$, it follows that a distinct 2^nd-level adjoint sensitivity function will correspond to each distinct value of the index $i=1,\mathrm{...},TF$. This fact has been explicitly indicated by including the index $i=1,\mathrm{...},TF$ in the list of arguments of the components of the 2^nd-level adjoint sensitivity function $c^{\left(2\right)}\left(i;n;t\right)\triangleq {\left[c_1^{\left(2\right)}\left(i;n;t\right),c_2^{\left(2\right)}\left(i;n;t\right)\right]}^{†}$:

4.

The definition of the 2^nd-level adjoint sensitivity function $c^{\left(2\right)}\left(i;n;t\right)\triangleq {\left[c_1^{\left(2\right)}\left(i;n;t\right),c_2^{\left(2\right)}\left(i;n;t\right)\right]}^{†}$ is now completed by requiring that it satisfy the following boundary conditions, which eliminates the respective unknown terms on the right-side of Eq. (53):

$$c_1^{\left(2\right)}\left(i;n;t_f\right)=0\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{c}_2^{\left(2\right)}\left(i;n;t_0\right)=0\hspace{0.17em}.$$

(84)

The system of equations comprising Eqs. (82)‒(84) constitutes the “2^nd-Level Adjoint Sensitivity System (2^nd-LASS)” for the 2^nd-level adjoint sensitivity function $c^{\left(2\right)}\left(i;n;t\right)\triangleq {\left[c_1^{\left(2\right)}\left(i;n;t\right),c_2^{\left(2\right)}\left(i;n;t\right)\right]}^{†}$. This 2^nd-LASS is independent of parameter variations, is linear in $c^{\left(2\right)}\left(i;n;t\right)\triangleq {\left[c_1^{\left(2\right)}\left(i;n;t\right),c_2^{\left(2\right)}\left(i;n;t\right)\right]}^{†}$, and can be solved sequentially, by first solving Eq. (83) subject to the initial condition given in Eq. (84) to determine the function $c_2^{\left(2\right)}\left(i;n;t\right)$, and subsequently using the function $c_2^{\left(2\right)}\left(i;n;t\right)$ in Eq. (82) to solve it subject to the “final-time” condition given in Eq. (84) to obtain the function $c_1^{\left(2\right)}\left(i;n;t\right)$. The 2^nd-LASS is to be solved using the nominal values of all functions and parameters/weights, but the superscript “zero” has been omitted for notational simplicity.

5.

Using the results obtained in Eqs. (80), (82)‒(84) in Eq. (81) yields the following alternative expression for the “indirect-effect” term $\delta {\left(\partial r_n/\partial w_i\right)}_{ind}$, which no longer involves the 2^nd-level variational sensitivity function $v^{\left(2,n\right)}\left(t\right)$ but involves the 2^nd-level adjoint sensitivity function $c^{\left(2\right)}\left(i;n;t\right)\triangleq {\left[c_1^{\left(2\right)}\left(i;n;t\right),c_2^{\left(2\right)}\left(i;n;t\right)\right]}^{†}$:

$$\begin{array}{l}\delta {\left({\displaystyle \frac{\partial r_n}{\partial w_i}}\right)}_{ind}={\left[c_1^{\left(2\right)}\left(i;n;t_0\right)\right]}^{†}{\displaystyle \frac{\partial h_e\left(x,w\right)}{\partial x}}\delta x+{\left[c_1^{\left(2\right)}\left(i;n;t_0\right)\right]}^{†}{\displaystyle \frac{\partial h_e\left(x,w\right)}{\partial w}}\delta w \\ +{\displaystyle {\int }_{t_0}^{t_f}{\left[c_1^{\left(2\right)}\left(i;n;t\right)\right]}^{†}{Q}^{\left(1\right)}\delta F}\hspace{0.17em}dt-{\displaystyle {\int }_{t_0}^{t_f}{\left[c_2^{\left(2\right)}\left(i;n;t\right)\right]}^{†}}\hspace{0.17em}{p}^{\left(2,n\right)}dt.\end{array}$$

(85)

6.

Adding the expression obtained in Eq. (85) to the expression for the direct-effect term $\delta {\left(\partial r_n/\partial w_i\right)}_{dir}$ provided in Eq. (78) yields the following expression for the total first-order G-differential $\delta \left(\partial r_n/\partial w_i\right)$, for each $n=1,\mathrm{...},TR$:

$$\begin{array}{l}\delta \left({\displaystyle \frac{\partial r_n}{\partial w_i}}\right)={\left[a^{(1)}\left(n;t_0\right)\right]}^{†}\left[{\displaystyle \frac{{\partial }^2{h}_e\left(x,w\right)}{\partial x\partial w_i}}\delta x+{\displaystyle \frac{{\partial }^2{h}_e\left(x,w\right)}{\partial w\partial w_i}}\delta w\right]+{\left[c_1^{\left(2\right)}\left(i;n;t_0\right)\right]}^{†}{\displaystyle \frac{\partial h_e\left(x,w\right)}{\partial x}}\delta x\\ +{\left[c_1^{\left(2\right)}\left(i;n;t_0\right)\right]}^{†}{\displaystyle \frac{\partial h_e\left(x,w\right)}{\partial w}}\delta w +{\displaystyle {\int }_{t_0}^{t_f}{\left[c_1^{\left(2\right)}\left(i;n;t\right)\right]}^{†}{Q}^{\left(1\right)}\delta F}\hspace{0.17em}dt-{\displaystyle {\int }_{t_0}^{t_f}{\left[c_2^{\left(2\right)}\left(i;n;t\right)\right]}^{†}}\hspace{0.17em}{p}^{\left(2,n\right)}dt\\ \triangleq {\displaystyle \frac{{\partial }^2{r}_n}{\partial F\partial F_i}}\delta F+{\displaystyle \frac{{\partial }^2{r}_n}{\partial x\partial F_i}}\delta x+{\displaystyle \frac{{\partial }^2{r}_n}{\partial w\partial F_i}}\delta w+{\displaystyle \frac{{\partial }^2{r}_n}{\partial \phi \partial F_i}}\delta \phi \\ ={\displaystyle \sum _{j=1}^{TF}\frac{{\partial }^2{r}_n}{\partial F_j\partial F_i}}\delta F_j+{\displaystyle \sum _{j=1}^{TI}\frac{{\partial }^2{r}_n}{\partial x_j\partial F_i}}\delta x_j+{\displaystyle \sum _{j=1}^{TEW}\frac{{\partial }^2{r}_n}{\partial w_j\partial F_i}}\delta w_j+{\displaystyle \sum _{j=1}^{TD}\frac{{\partial }^2{r}_n}{\partial {\phi }_j\partial F_i}}\delta {\phi }_j;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},TF.\end{array}$$

(86)

The expression shown in Eq. (86) is to be evaluated at the nominal values of all functions and parameters/weights but the superscript “zero” (which has been used to indicate this fact) has been omitted for notational simplicity.

Identifying in Eq. (86) the expressions that multiply the individual variations $\delta F$, $\delta x$, $\delta w$ and $\delta \phi$ yields the following expressions for the second-order sensitivities that stem from the first-order sensitivity $\partial r_n/\partial F_i$ for $n=1,\mathrm{...},TR$:

$$\begin{array}{l}{\displaystyle \frac{{\partial }^2{r}_n}{\partial F_j\partial w_i}}={\displaystyle \sum _{k=1}^{TH}{\displaystyle {\int }_{t_0}^{t_f}{c}_{1,k}^{\left(2\right)}\left(i;n;t\right)\frac{\partial f_k\left(h;F\right)}{\partial F_j}dt-{\displaystyle \sum _{m=1}^{TH}{\displaystyle {\int }_{t_0}^{t_f}{c}_{2,m}^{\left(2\right)}\left(i;n;t\right)}\hspace{0.17em}{\displaystyle \sum _{k=1}^{TH}{a}_k^{\left(1\right)}\left(n;t\right)\frac{{\partial }^2{f}_k\left(h;F\right)}{\partial F_j\partial h_m}}dt}}}\\ \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}j=1,\mathrm{...},TF;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}n=1,\mathrm{...},TR;\end{array}$$

(87)

$$\begin{array}{l}{\displaystyle \frac{{\partial }^2{r}_n}{\partial x_j\partial w_i}}={\left[a^{(1)}\left(n;t_0\right)\right]}^{†}{\displaystyle \frac{{\partial }^2{h}_e\left(x,w\right)}{\partial x_j\partial w_i}}+{\left[c_1^{\left(2\right)}\left(i;n;t_0\right)\right]}^{†}{\displaystyle \frac{\partial h_e\left(x,w\right)}{\partial x_j}}\\ ={\displaystyle \sum _{k=1}^{TH}{a}_k^{\left(1\right)}\left(n;t_0\right)}{\displaystyle \frac{{\partial }^2{h}_k^e\left(x,w\right)}{\partial x_j\partial w_i}}+{\displaystyle \sum _{k=1}^{TH}{c}_{1,k}^{\left(2\right)}\left(i;n;t_0\right)}{\displaystyle \frac{\partial h_k^e\left(x,w\right)}{\partial x_j}}\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}j=1,\mathrm{....},TI;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}n=1,\mathrm{...},TR;\end{array}$$

(88)

$$\begin{array}{l}{\displaystyle \frac{{\partial }^2{r}_n}{\partial w_j\partial w_i}}={\left[a^{(1)}\left(n;t_0\right)\right]}^{†}{\displaystyle \frac{{\partial }^2{h}_e\left(x,w\right)}{\partial w_j\partial w_i}}+{\left[c_1^{\left(2\right)}\left(i;n;t_0\right)\right]}^{†}{\displaystyle \frac{\partial h_e\left(x,w\right)}{\partial w_j}} \\ ={\displaystyle \sum _{k=1}^{TH}{a}_k^{\left(1\right)}\left(n;t_0\right)}{\displaystyle \frac{{\partial }^2{h}_k^e\left(x,w\right)}{\partial w_j\partial w_i}}+{\displaystyle \sum _{k=1}^{TH}{c}_{1,k}^{\left(2\right)}\left(i;n;t_0\right)}{\displaystyle \frac{\partial h_k^e\left(x,w\right)}{\partial w_j}}\hspace{0.17em};\hspace{0.17em}i,j=1,\mathrm{....},TEW\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}n=1,\mathrm{...},TR;\end{array}$$

(89)

$$\begin{array}{l}{\displaystyle \frac{{\partial }^2{r}_n}{\partial {\phi }_j\partial w_i}}=-{\displaystyle {\int }_{t_0}^{t_f}{\left[c_2^{\left(2\right)}\left(i;n;t\right)\right]}^{†}}{\displaystyle \frac{{\partial }^2{h}_n^d\left(h;\phi \right)}{\partial {\phi }_j\partial h}}\delta \left(t-t_f\right)\hspace{0.17em}dt \\ =-{\displaystyle \sum _{k=1}^{TH}{c}_{1,k}^{\left(2\right)}\left(i;n;t_f\right)\frac{{\partial }^2{h}_n^d\left[h\left(t_f\right);\phi \right]}{\partial {\phi }_j\partial h_k\left(t_f\right)}};\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{....},TEW;\hspace{0.17em}\hspace{0.17em}j=1,\mathrm{....},TD;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}n=1,\mathrm{...},TR.\end{array}$$

(90)

The second-order sensitivities of the responses with respect to the primary parameters/weights $\theta \triangleq {\left[{\theta }_1,\mathrm{...},{\theta }_{TW}\right]}^{†}$ are obtained analytically by using the results obtained in Eqs. (87)‒(90) in conjunctions with the “chain-rule” $\partial \left[\right]/\partial {\theta }_m={\displaystyle \sum _{\mu =1}^{TF}\left(\partial \left[\right]/\partial F_{\mu }\right)\left(\partial F_{\mu }/\partial {\theta }_m\right)}$ to the expression in Eq. (87). In the absence of feature functions, the second-order sensitivities are obtained by setting $F_i\left(\theta \right)\equiv {\theta }_i$, for all $i=1,\mathrm{...},TW\equiv TF$.

3.4. Second-Order Sensitivities Stemming from $\partial r_n/\partial x_i$, $i=1,\mathrm{...},TI$;$n=1,\mathrm{...},TR$

The second-order sensitivities stemming from the first-order sensitivities $\partial r_n/\partial x_i$ are obtained from the first-order G-differential of Eq. (35), for $\hspace{0.17em}i=1,\mathrm{...},TI$;$n=1,\mathrm{...},TR$, as follows:

$$\begin{array}{l}\delta \left({\displaystyle \frac{\partial r_n}{\partial x_i}}\right)\triangleq {\left\{{\displaystyle \frac{d}{d\epsilon }}\left[{\displaystyle \frac{\partial h_e\left(x^0+\epsilon \delta x,w^0+\epsilon \delta w\right)}{\partial x_i}}{\displaystyle {\int }_{t_0}^{t_f}{\left[a^{(1,0)}\left(n;t\right)+\epsilon \delta a^{(1)}\left(n;t\right)\right]}^{†}\delta \left(t-t_0\right)}\hspace{0.17em}dt\right]\right\}}_{\epsilon =0}\\ \triangleq \delta {\left({\displaystyle \frac{\partial r_n}{\partial x_i}}\right)}_{dir}+\delta {\left({\displaystyle \frac{\partial r_n}{\partial x_i}}\right)}_{ind},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},TI;\hspace{0.17em}\hspace{0.17em}n=1,\mathrm{...},TR,\end{array}$$

(91)

where the direct-effect term $\delta {\left(\partial r_n/\partial x_i\right)}_{dir}$ comprises the variations $\delta x$ and $\delta w$, defined as follows:

$$\delta {\left({\displaystyle \frac{\partial r_n}{\partial x_i}}\right)}_{dir}\triangleq a^{(1)}\left(n;t_0\right)\left[{\displaystyle \frac{{\partial }^2{h}_e\left(x,w\right)}{\partial x\partial x_i}}\delta x+{\displaystyle \frac{{\partial }^2{h}_e\left(x,w\right)}{\partial w\partial x_i}}\delta w\right];$$

(92)

and where the indirect-effect term $\delta {\left(\partial r_n/\partial x_i\right)}_{ind}$ comprises the variation $\delta a^{(1)}\left(n;t\right)$ and is defined as follows:

$$\delta {\left({\displaystyle \frac{\partial r_n}{\partial x_i}}\right)}_{ind}\triangleq {\displaystyle \frac{\partial h_e\left(x,w\right)}{\partial x_i}}{\displaystyle {\int }_{t_0}^{t_f}{\left[\delta a^{(1)}\left(n;t\right)\right]}^{†}\delta \left(t-t_0\right)}\hspace{0.17em}dt.$$

(93)

The expressions shown in Eqs. (78) and (79) are to be evaluated at the nominal values of the respective functions and parameters, but the respective indication (namely: the superscript “zero”) has been omitted in order to simplify the notation.

The direct-effect term $\delta {\left(\partial r_n/\partial x_i\right)}_{dir}$ can be evaluated at this time for all variations $\delta x$ and $\delta w$, but the indirect-effect term $\delta {\left(\partial r_n/\partial x_i\right)}_{ind}$ can be evaluated only after having determined the variation $\delta a^{(1)}\left(n;t\right)$, which is the solution of the 2^nd-LVSS defined by Eqs. (48) and (50). The need for solving this 2^nd-LVSS can be avoided if the indirect-effect term defined in Eq. (93) could be expressed in terms of a “right-hand side” that does not involve the function $\delta a^{(1)}\left(n;t\right)$. This goal can be achieved by applying the same concepts and steps as in Section 3.3, which will not be shown in detail here in order to minimize repetitive derivations.

The final expressions of the second-order sensitivities stemming from $\partial r_n/\partial x_i$, $i=1,\mathrm{...},TI$;$n=1,\mathrm{...},TR$ are as follows:

$$\begin{array}{l}{\displaystyle \frac{{\partial }^2{r}_n}{\partial F_j\partial x_i}}={\displaystyle \sum _{k=1}^{TH}{\displaystyle {\int }_{t_0}^{t_f}{d}_{1,k}^{\left(2\right)}\left(i;n;t\right)\frac{\partial f_k\left(h;F\right)}{\partial F_j}dt-{\displaystyle \sum _{m=1}^{TH}{\displaystyle {\int }_{t_0}^{t_f}{d}_{2,m}^{\left(2\right)}\left(i;n;t\right)}\hspace{0.17em}{\displaystyle \sum _{k=1}^{TH}{a}_k^{\left(1\right)}\left(n;t\right)\frac{{\partial }^2{f}_k\left(h;F\right)}{\partial F_j\partial h_m}}dt}}}\\ \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}j=1,\mathrm{...},TF;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}n=1,\mathrm{...},TR;\end{array}$$

(94)

$$\begin{array}{l}{\displaystyle \frac{{\partial }^2{r}_n}{\partial x_j\partial x_i}}={\left[a^{(1)}\left(n;t_0\right)\right]}^{†}{\displaystyle \frac{{\partial }^2{h}_e\left(x,w\right)}{\partial x_j\partial x_i}}+{\left[d_1^{\left(2\right)}\left(i;n;t_0\right)\right]}^{†}{\displaystyle \frac{\partial h_e\left(x,w\right)}{\partial x_j}}\\ ={\displaystyle \sum _{k=1}^{TH}{a}_k^{\left(1\right)}\left(n;t_0\right)}{\displaystyle \frac{{\partial }^2{h}_k^e\left(x,w\right)}{\partial x_j\partial x_i}}+{\displaystyle \sum _{k=1}^{TH}{d}_{1,k}^{\left(2\right)}\left(i;n;t_0\right)}{\displaystyle \frac{\partial h_k^e\left(x,w\right)}{\partial x_j}}\hspace{0.17em};\\ \hspace{0.17em}\hspace{0.17em}i,j=1,\mathrm{....},TI;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}n=1,\mathrm{...},TR;\end{array}$$

(95)

$$\begin{array}{l}{\displaystyle \frac{{\partial }^2{r}_n}{\partial w_j\partial x_i}}={\left[a^{(1)}\left(n;t_0\right)\right]}^{†}{\displaystyle \frac{{\partial }^2{h}_e\left(x,w\right)}{\partial w_j\partial x_i}}+{\left[d_1^{\left(2\right)}\left(i;n;t_0\right)\right]}^{†}{\displaystyle \frac{\partial h_e\left(x,w\right)}{\partial w_j}} \\ ={\displaystyle \sum _{k=1}^{TH}{a}_k^{\left(1\right)}\left(n;t_0\right)}{\displaystyle \frac{{\partial }^2{h}_k^e\left(x,w\right)}{\partial w_j\partial x_i}}+{\displaystyle \sum _{k=1}^{TH}{d}_{1,k}^{\left(2\right)}\left(i;n;t_0\right)}{\displaystyle \frac{\partial h_k^e\left(x,w\right)}{\partial w_j}}\hspace{0.17em};\hspace{0.17em}\\ \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}j=1,\mathrm{....},TEW\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}n=1,\mathrm{...},TR;\end{array}$$

(96)

$$\begin{array}{l}{\displaystyle \frac{{\partial }^2{r}_n}{\partial {\phi }_j\partial x_i}}=-{\displaystyle {\int }_{t_0}^{t_f}{\left[d_2^{\left(2\right)}\left(i;n;t\right)\right]}^{†}}{\displaystyle \frac{{\partial }^2{h}_n^d\left(h;\phi \right)}{\partial {\phi }_j\partial h}}\delta \left(t-t_f\right)\hspace{0.17em}dt \\ =-{\displaystyle \sum _{k=1}^{TH}{d}_{1,k}^{\left(2\right)}\left(i;n;t_f\right)\frac{{\partial }^2{h}_n^d\left[h\left(t_f\right);\phi \right]}{\partial {\phi }_j\partial h_k\left(t_f\right)}};\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{....},TI;\hspace{0.17em}\hspace{0.17em}j=1,\mathrm{....},TD;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}n=1,\mathrm{...},TR.\end{array}$$

(97)

The 2^nd-level adjoint sensitivity function $d^{\left(2\right)}\left(i;n;t\right)\triangleq {\left[d_1^{\left(2\right)}\left(i;n;t\right),d_2^{\left(2\right)}\left(i;n;t\right)\right]}^{†}$, with $d_1^{\left(2\right)}\left(n;t\right)\triangleq {\left[d_{1,1}^{\left(2\right)}\left(n;t\right),\dots ,d_{1,TH}^{\left(2\right)}\left(n;t\right)\right]}^{†}$ and $d_2^{\left(2\right)}\left(n;t\right)\triangleq {\left[d_{2,1}^{\left(2\right)}\left(n;t\right),\dots ,d_{2,TH}^{\left(2\right)}\left(n;t\right)\right]}^{†}$, which appears in Eqs. (94)‒(97) is the solution of the following 2^nd-Level Adjoint Sensitivity System (2^nd-LASS):

$$-{\displaystyle \frac{d}{dt}}{d}_1^{\left(2\right)}\left(i;n;t\right)-J^{†}{d}_1^{\left(2\right)}\left(i;n;t\right)+{\left[G^{\left(n\right)}\right]}^{†}{d}_2^{\left(2\right)}\left(i;n;t\right)=0;$$

(98)

$$-{\displaystyle \frac{d}{dt}}{d}_2^{\left(2\right)}\left(i;n;t\right)+J{d}_2^{\left(2\right)}\left(i;n;t\right)={\displaystyle \frac{\partial h_e\left(x,w\right)}{\partial x_i}}\delta \left(t-t_0\right).$$

(99)

$$d_1^{\left(2\right)}\left(i;n;t_f\right)=0\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{d}_2^{\left(2\right)}\left(i;n;t_0\right)=0\hspace{0.17em}.$$

(100)

3.5. Discussion: Double-Computation of the Mixed Second-Order Sensitivities

It is evident from the derivations presented in Subsections 3.1 through 3.4 that the unmixed second-order sensitivities are computed just once, but the mixed second-order sensitivities are obtained/computed twice, having obtained two distinct expressions involving distinct 2^nd-level adjoint sensitivity functions for each mixed second-order sensitivity. The specific distinct expressions for the respective mixed second-order sensitivities are as follows:

• The mixed second-order sensitivities ${\partial }^2{r}_n/\partial x_j\partial F_i$ are obtained in Eq. (60) in terms of the 2^nd-level adjoint sensitivity functions $a_{1,k}^{\left(2\right)}\left(i;n;t\right)$,$k=1,\mathrm{....},TH;$$i=1,\mathrm{....},TF;$$j=1,\mathrm{....},TI;$$n=1,\mathrm{...},TR$. On the other hand, the mixed second-order sensitivities ${\partial }^2{r}_n/\partial F_j\partial x_i$ are obtained in Eq. (94) in terms of the 2^nd-level adjoint sensitivity functions $d_{1,k}^{\left(2\right)}\left(i;n;t\right)$ and $d_{2,k}^{\left(2\right)}\left(i;n;t\right)$, for $\hspace{0.17em}i=1,\mathrm{...},TI;$$j=1,\mathrm{...},TF;$$n=1,\mathrm{...},TR$. Due to the symmetry property of the mixed second-order sensitivities, the numerical results obtained for the corresponding mixed second-order sensitivities by computing them using Eq. (60), on the one hand, and using Eq. (94), on the other hand, provides a verification mechanism for assessing the accuracy of the computation of the 2^nd-level adjoint functions involved in the respective computations.
• The mixed second-order sensitivities ${\partial }^2{r}_n/\partial w_j\partial F_i$ are obtained in Eq. (61) in terms of the 2^nd-level adjoint sensitivity functions $a_{1,k}^{\left(2\right)}\left(i;n;t\right)$, for $k=1,\mathrm{....},TH;$$i=1,\mathrm{....},TF;$$j=1,\mathrm{....},TI;$$n=1,\mathrm{...},TR$. On the other hand, the mixed second-order sensitivities ${\partial }^2{r}_n/\partial F_j\partial w_i$ are obtained in Eq. (87) in terms of the 2^nd-level adjoint sensitivity functions $c_{1,k}^{\left(2\right)}\left(i;n;t\right)$ and $c_{2,k}^{\left(2\right)}\left(i;n;t\right)$, for $k=1,\mathrm{....},TH;$$\hspace{0.17em}i=1,\mathrm{...},TEW;$$j=1,\mathrm{...},TF;$$n=1,\mathrm{...},TR$. Due to the symmetry property of the mixed second-order sensitivities, the numerical results obtained for the corresponding mixed second-order sensitivities by computing them using Eq. (61), on the one hand, and using Eq. (87), on the other hand, provides a verification mechanism for assessing the accuracy of the computation of the 2^nd-level adjoint functions involved in the respective computations.
• The mixed second-order sensitivities ${\partial }^2{r}_n/\partial {\phi }_j\partial F_i$ are obtained in Eq. (62) in terms of the 2^nd-level adjoint sensitivity functions $a_{1,k}^{\left(2\right)}\left(i;n;t\right)$, for $k=1,\mathrm{....},TH;$$i=1,\mathrm{....},TF;$$j=1,\mathrm{....},TD;$$n=1,\mathrm{...},TR$. On the other hand, the mixed second-order sensitivities ${\partial }^2{r}_n/\partial F_j\partial {\phi }_i$ are obtained in Eq. (73) in terms of the 2^nd-level adjoint sensitivity functions $b_k^{\left(2\right)}\left(i;n;t\right)$, for $k=1,\mathrm{....},TH;$$\hspace{0.17em}i=1,\mathrm{...},TD;$$j=1,\mathrm{...},TF;$$n=1,\mathrm{...},TR$. Due to the symmetry property of the mixed second-order sensitivities, the numerical results obtained for the corresponding mixed second-order sensitivities by computing them using Eq. (62), on the one hand, and using Eq. (73), on the other hand, provides a verification mechanism for assessing the accuracy of the computation of the 2^nd-level adjoint functions involved in the respective computations.
• The mixed second-order sensitivities ${\partial }^2{r}_n/\partial x_j\partial w_i$ are obtained in Eq. (88) in terms of the adjoint sensitivity functions $a_k^{\left(1\right)}\left(n;t\right)$ and $c_{1,k}^{\left(2\right)}\left(i;n;t\right)$, for $k=1,\mathrm{....},TH;$$i=1,\mathrm{....},TEW;$$j=1,\mathrm{....},TI;$$n=1,\mathrm{...},TR$. On the other hand, the mixed second-order sensitivities ${\partial }^2{r}_n/\partial w_j\partial x_i$ are obtained in Eq. (96) in terms of the adjoint sensitivity functions $a_k^{\left(1\right)}\left(n;t\right)$ and $d_{1,k}^{\left(2\right)}\left(i;n;t\right)$, for $k=1,\mathrm{....},TH;$$\hspace{0.17em}i=1,\mathrm{...},TI;$$j=1,\mathrm{...},TEW;$$n=1,\mathrm{...},TR$. Due to the symmetry property of the mixed second-order sensitivities, the numerical results obtained for the corresponding mixed second-order sensitivities by computing them using Eq. (74), on the one hand, and using Eq. (97), on the other hand, provides a verification mechanism for assessing the accuracy of the computation of the adjoint functions involved in the respective computations.
• The mixed second-order sensitivities ${\partial }^2{r}_n/\partial x_j\partial {\phi }_i$ are obtained in Eq. (74) in terms of the 2^nd-level adjoint sensitivity functions $b_k^{\left(2\right)}\left(i;n;t\right)$, for $k=1,\mathrm{....},TH;$$i=1,\mathrm{....},TD;$$j=1,\mathrm{....},TI;$$n=1,\mathrm{...},TR$. On the other hand, the mixed second-order sensitivities ${\partial }^2{r}_n/\partial {\phi }_j\partial x_i$ are obtained in Eq. (97) in terms of the 2^nd-level adjoint sensitivity functions $d_{1,k}^{\left(2\right)}\left(i;n;t\right)$, for $k=1,\mathrm{....},TH;$$\hspace{0.17em}i=1,\mathrm{...},TI;$$j=1,\mathrm{...},TD;$$n=1,\mathrm{...},TR$. Due to the symmetry property of the mixed second-order sensitivities, the numerical results obtained for the corresponding mixed second-order sensitivities by computing them using Eq. (74), on the one hand, and using Eq. (97), on the other hand, provides a verification mechanism for assessing the accuracy of the computation of the 2^nd-level adjoint functions involved in the respective computations.
• The mixed second-order sensitivities ${\partial }^2{r}_n/\partial w_j\partial {\phi }_i$ are obtained in Eq. (75) in terms of the 2^nd-level adjoint sensitivity functions $b_k^{\left(2\right)}\left(i;n;t\right)$, for $k=1,\mathrm{....},TH;$$i=1,\mathrm{....},TD;$$j=1,\mathrm{....},TEW;$$n=1,\mathrm{...},TR$. On the other hand, the mixed second-order sensitivities ${\partial }^2{r}_n/\partial {\phi }_j\partial w_i$ are obtained in Eq. (90) in terms of the 2^nd-level adjoint sensitivity functions $c_{1,k}^{\left(2\right)}\left(i;n;t\right)$, for $k=1,\mathrm{....},TH;$$\hspace{0.17em}i=1,\mathrm{...},TEW;$$j=1,\mathrm{...},TD;$$n=1,\mathrm{...},TR$. Due to the symmetry property of the mixed second-order sensitivities, the numerical results obtained for the corresponding mixed second-order sensitivities by computing them using Eq. (75), on the one hand, and using Eq. (90), on the other hand, provides a verification mechanism for assessing the accuracy of the computation of the 2^nd-level adjoint functions involved in the respective computations.

4. Discussion and Conclusions

As has been discussed in this work, after the NODE-net is optimized to reproduce the underlying physical system as closely as possible, the subsequent NODE-responses of interest are various functionals of the NODE’s “decoder” output rather than some “loss function.” The physical system modeled by the NODE-net comprises parameters that stem from measurements and/or computations. Such parameters are not perfectly well known but are subject to uncertainties that stem from the respective experiments and/or computations. These uncertainties will induce uncertainties in the NODE decoder’s output/response, which can only be quantified if the sensitivities of the NODE decoder’s response with respect to the optimized NODE weights/parameters are available.

It is well known that equations modeling real-life phenomena comprise not only scalar-valued model parameters but also functions of such scalar model parameters. For example, conservation equations modeling thermal-hydraulics phenomena depend on correlations, Reynolds, Nusselt numbers, etc. Similarly, neutron transport and/or diffusion equations in nuclear reactor engineering (e.g., reactor physics and shielding) involve group-averaged macroscopic cross section, which are scalar-valued functions of various sums of scalar-valued isotopic number densities and microscopic cross sections. Such scalar-valued functions are conveniently called “features of primary model parameters.” This work has introduced the mathematical framework of the “First-Order Features Adjoint Sensitivity Analysis Methodology for Neural Ordinary Differential Equations (1^st-FASAM-NODE)”, which derives and computes most efficiently the exact expressions of all of the first-order sensitivities of NODE-decoder responses with respect to the features of parameters (weights, initial conditions, etc.) that characterize the NODE’s decoder, hidden layers, and encoder. It has been shown that a single large-scale computation to solve the 1^st-Order Adjoint Sensitivity System(1^st-LASS) suffices to obtain the first-order sensitivities of a decoder-response with respect to all of the feature functions of parameters. Subsequently, the sensitivities of the responses with respect to the primary model parameters are determined, analytically and trivially, by applying the “chain-rule” to the expressions of the sensitivities with respect to the model’s “features/functions of parameters.”

This work has also presented the “Second-Order Features Adjoint Sensitivity Analysis Methodology for Neural Ordinary Differential Equations (2^nd-FASAM-NODE),” which builds upon the 1^st-FASAM-NODE methodology to derive and compute most efficiently the exact expressions of all of the second-order sensitivities of NODE-decoder responses with respect to the features of parameters (weights, initial conditions, etc.) that characterize the NODE’s decoder, hidden layers, and encoder. It has been shown that, for each decoder response, the computation of all of the second-order sensitivities requires as many large-scale computations ( to solve the 2^nd-Level Adjoint Sensitivity System) as there are first-order sensitivities of the decoder-response with respect to the feature functions, initial conditions, and decoder-weights. The 2^nd-FASAM-NODE methodology computes the mixed second-order sensitivities twice, involving distinct 2^nd-level adjoint sensitivity functions in distinct expressions.

The concepts underlying the 1^st-FASAM-NODE, and the 2^nd-FASAM-NODE methodologies will be illustrated in the accompanying “Part II” [18] by considering the energy and heat transfer processes described by the well-known Nordheim-Fuchs phenomenological model [19,20] of reactor safety.