Illustrative Application of the 2^nd-Order Adjoint Sensitivity Analysis Methodology to a Paradigm Linear Evolution/Transmission Model: Reaction-Rate Detector Response

Dan Gabriel Cacuci

Center for Nuclear Science and Energy, Department of Mechanical Engineering, University of South Carolina, Columbia SC, USA.

DOI: 10.4236/ajcm.2020.103020                

Abstract

This work continues the illustrative application of the “Second Order Comprehensive Adjoint Sensitivity Analysis Methodology” (2^nd-CASAM) to a benchmark mathematical model that can simulate the evolution and/or transmission of particles in a heterogeneous medium. The model response considered in this work is a reaction-rate detector response, which provides the average interactions of particles with the respective detector or, alternatively, the time-average of the concentration of a mixture of substances in a medium. The definition of this model response includes both uncertain boundary points of the benchmark, thereby providing both direct and indirect contributions to the response sensitivities stemming from the boundaries. The exact expressions for the 1^st- and 2^nd-order response sensitivities to the boundary and model parameters obtained in this work can serve as stringent benchmarks for inter-comparing the performances of all (deterministic and statistical) sensitivity analysis methods.

Keywords

Second-Order Adjoint Comprehensive Sensitivity Analysis Methodology (2^nd-CASAM), Evolution Benchmark Model, Exact and Efficient Computation of First- and Second-Order Response Sensitivities

1. Introduction

This work continues to illustrate the application of the general second-order adjoint sensitivity analysis methodology (2^nd-CASAM) presented in [1] by using the evolution/transmission mathematical benchmark model introduced in [2], but considering a “reaction-rate” detector response, as opposed to the “point-detector” response considered in [2]. As in [2], the mathematical model considered in this work could represent [3] [4] the time-evolution of the concentration of a substance in a homogeneous mixture of materials or, alternatively, it could represent [4] [5] [6] the transmission/attenuation of the flux of uncollided particles (e.g., photons) travelling through a one-dimensional homogenized multi-material slab of imprecisely known thickness.

Although simple, the model comprises a large number of model parameters, thereby involving a correspondingly large number of sensitivities (i.e., functional derivatives) of the model’s responses to the model parameters. Furthermore, the model has been deliberately designed so that a large number of relative response sensitivities display identical values. The application of the 2^nd-CASAM [1] yields the exact expressions of the 1^st- and 2^nd-order sensitivities of the reaction-rate response to the uncertain model and boundary parameters. These exact expressions can be used to benchmark any other statistical or deterministic software used for computing sensitivities.

This work is organized as follows: Section 2 presents the mathematical formulation and expression of a reaction-rate detector for the paradigm evolution/transmission model. Section 3 illustrates the application of the 2^nd-CASAM [1] for obtaining efficiently the exact closed-form expressions of the first- and second-order sensitivities of the reaction-rate detector response to both model and boundary parameters. Section 4 offers concluding remarks.

2. Mathematical Modeling of a Paradigm Evolution/Transmission Benchmark Problem

The general methodology presented in Part I (Cacuci, 2020) is applied in this work to a simple paradigm model, admitting a closed-form analytic solution for convenient verification of all results to be obtained, which simulates a typical evolution or attenuation of a quantity that will be denoted as $\rho \left(t\right)$, satisfying the following linear conservation equation:

$\frac{\text{d}\rho \left(t\right)}{\text{d}t}+\rho \left(t\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}=0,0\le {\beta }_{\mathcal{l}}\le t\le {\beta }_u<\infty ,$ (1)

$\rho \left({\beta }_{\mathcal{l}}\right)={\rho }_{in},\text{at}t={\beta }_{\mathcal{l}}.$ (2)

As has been discussed in [2], Equations (1) and (2) occur in the mathematical modeling of many physical systems, including the evolution of the concentration of a substance in a homogeneous mixture of N materials (from an imprecisely known initial quantity, denoted as ${\rho }_{in}$, measured at an initial-time value $t={\beta }_{\mathcal{l}}$ towards an imprecisely known final-time value $t={\beta }_u$ ) or the mono-directional propagation of the flux of uncollided particles travelling through a one-dimensional homogenized multi-material slab of imprecisely known thickness $\left({\beta }_u-{\beta }_{\mathcal{l}}\right)$ in a direction parallel to the t-coordinate.

An important typical response of interest for the physical problem modeled by Equations (1) and (2) is a “reaction rate” detector response, which will be denoted as $R_2\left(\rho ;\alpha ,\beta \right)$, and which is represented mathematically by the following functional:

$R_2\left(\rho ;\alpha ,\beta \right)={\Sigma }_d{\displaystyle \underset{{\beta }_{\mathcal{l}}}{\overset{{\beta }_u}{\int }}\rho \left(t\right)}\text{d}t,$ (3)

where ${\Sigma }_d$ represents the detector’s reaction cross section. In Equation (3), the vector $\alpha$ denotes the “vector of model parameters” and defined as follows:

$\alpha \triangleq {\left({\alpha }_1,\cdots ,{\alpha }_{N_{\alpha }}\right)}^{†}\triangleq {\left(n_1,\cdots ,n_N,{\sigma }_1,\cdots ,{\sigma }_N,{\rho }_{in},{\Sigma }_d\right)}^{†}.$ (4)

Similarly, the vector $\beta$ denotes the “vector of boundary parameters” and is defined as follows:

$\beta \triangleq {\left({\beta }_{\mathcal{l}},{\beta }_u\right)}^{†}.$ (5)

Throughout this work, the symbol “ $\triangleq$ ” is used to denote “is defined as” or “is by definition,” while the “dagger” $(†)$ superscript is used to denote “transposition.”

For subsequent verification of the expressions that will be obtained for various response sensitivities, the closed-form solution of Equations (1) and (2) is provided below:

$\rho \left(t\right)={\rho }_{in}\exp \left[\left({\beta }_{\mathcal{l}}-t\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}\right].$ (6)

Although the model parameters ${\rho }_{in}$, $n_i$, ${\sigma }_i$, ${\Sigma }_d$, together with the boundary parameters ${\beta }_{\mathcal{l}}$ and ${\beta }_u$ are considered to be imperfectly known and subject to uncertainties, the actual probability distributions of these parameters are not known in practice. Usually, only the “nominal” (or “mean”) values and the respective variations from the nominal values (e.g., standard deviations) of the respective components are known. The nominal values will be denoted using the superscript “zero” so that the vector comprising the nominal values of the model parameters, denoted as ${\alpha }^0$, will be defined for the system under consideration as follows:

${\alpha }^0\triangleq {\left({\alpha }_1^0,\cdots ,{\alpha }_{N_{\alpha }}^0\right)}^{†}\triangleq {\left(n_1^0,\cdots ,n_N^0,{\sigma }_1^0,\cdots ,{\sigma }_N^0,{\rho }_{in}^0,{\Sigma }_d^0\right)}^{†}.$ (7)

Similarly, the vector comprising the nominal values of the boundary parameters is denoted as ${\beta }^0$ and is defined for the system under consideration as follows:

${\beta }^0\triangleq {\left({\beta }_{\mathcal{l}}^0,{\beta }_u^0\right)}^{†}$. (8)

In practice, the nominal solution, denoted as ${\rho }^0\left(t\right)$, is computed by solving numerically Equations (1) and (2) using the nominal values for the model and boundary parameters. For this illustrative example, the nominal solution of Equations (1) and (2) has the following expression:

${\rho }^0\left(t\right)={\rho }_{in}^0\exp \left[\left({\beta }_{\mathcal{l}}^0-t\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i^0{\sigma }_i^0}\right].$ (9)

The closed-form expression of $R_2\left(\rho ;\alpha ,\beta \right)$ is readily obtained by inserting into Equation (3) the expression of $\rho \left(t\right)$ provided in Equation (6) and performing the respective integration, which yields the following expression:

$R_2\left(\rho ;\alpha ,\beta \right)=\frac{{\rho }_{in}{\Sigma }_d}{{\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}}\left\{1-\exp \left[\left({\beta }_{\mathcal{l}}-{\beta }_u\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}\right]\right\}.$ (10)

As indicated be the expression in Equation (10), even though the forward function $\rho \left(t\right)$ is independent of the boundary point ${\beta }_u$, the response $R_2\left(\rho ;\alpha ,\beta \right)$ evidently depends on ${\beta }_u$, so its sensitivities with respect to ${\beta }_u$ will not vanish.

3. Application of the 2^nd-CASAM for Computing Exactly and Efficiently the 1^st- and 2^nd-Order Response Sensitivities of $R_2\left(\rho ;\alpha ,\beta \right)$ with Respect to the Uncertain Model and Boundary Parameters

The variations between the true and the nominal values of the model and boundary parameters will be considered to constitute the components of the vectors $\delta \alpha$ and $\delta \beta$, respectively, defined as follows:

$\delta \alpha \triangleq \left(\delta {\alpha }_1,\cdots ,\delta {\alpha }_{N_{\alpha }}\right),\delta {\alpha }_i\triangleq {\alpha }_i-{\alpha }_i^0,$ (11)

$\delta \beta \triangleq {\left(\delta {\beta }_{\mathcal{l}},\delta {\beta }_u\right)}^{†},\delta {\beta }_{\mathcal{l}}\triangleq {\beta }_{\mathcal{l}}-{\beta }_{\mathcal{l}}^0,\delta {\beta }_u\triangleq {\beta }_u-{\beta }_u^0.$ (12)

Since the state function is related to the model and boundary parameters $\alpha$ and $\beta$ through Equations (1) and (2), it follows that the variations and $\delta \beta$ in the model and boundary parameters will cause a corresponding variation in the state function $\rho \left(t\right)$ around the nominal solution ${\rho }^0\left(t\right)$. In turn, these variations will cause variations in the responses $R_1\left(\rho ;\alpha ,\beta \right)$ around the respective nominal response values. For subsequent derivations, it is convenient to use the compact notation $e\triangleq \left(\rho ;\alpha ,\beta \right)$, with the corresponding nominal values denoted as $e^0\triangleq \left({\rho }^0;{\alpha }^0,{\beta }^0\right)$.

3.1. Computing the 1^st-Order Sensitivities $R_2\left(\rho ;\alpha ,\beta \right)$ Using the 1^st-LASS

The first-order total sensitivity of the response $R_2\left(\rho ;\alpha ,\beta \right)$ defined in Equation (3) is provided by the first-order G-differential of this response, which is obtained by G-differentiating Equation (3) as shown in Part I (Cacuci, 2020), to obtain the following expression:

$\delta R_2\left(e^0;\delta \rho ;\delta \alpha ,\delta \beta \right)={\left(\delta R_2\right)}^{ind}+{\left(\delta R_2\right)}^{dir}$ (13)

where the indirect-effect term ${\left(\delta R_2\right)}^{ind}$ and the direct-effect term ${\left(\delta R_2\right)}^{dir}$, respectively, are defined as follows:

${\left(\delta R_2\right)}^{ind}\triangleq {\Sigma }_d^0{\displaystyle \underset{{\beta }_{\mathcal{l}}^0}{\overset{{\beta }_u^0}{\int }}\delta \rho \left(t\right)}\text{d}t,$ (14)

${\left(\delta R_2\right)}^{dir}\triangleq \left(\delta {\Sigma }_d\right){\displaystyle \underset{{\beta }_{\mathcal{l}}^0}{\overset{{\beta }_u^0}{\int }}\rho \left(t\right)}\text{d}t+{\Sigma }_d^0{\left\{\rho \left(t\right)\right\}}_{t={\beta }_u^0}\delta {\beta }_u-{\Sigma }_d^0{\left\{\rho \left(t\right)\right\}}_{t={\beta }_{\mathcal{l}}^0}\delta {\beta }_{\mathcal{l}}.$ (15)

The direct effect term, ${\left(\delta R_2\right)}^{dir}$, can be computed directly by using in Equation (15) the nominal solution provided in Equation (9) to obtain:

$\begin{array}{c}{\left(\delta R_2\right)}^{dir}\triangleq \left(\delta {\Sigma }_d\right)\frac{{\rho }_{in}^0}{{\displaystyle \underset{i=1}{\overset{N}{\sum }}{w}_i^0{\sigma }_i^0}}\left\{1-\exp \left[\left({\beta }_{\mathcal{l}}^0-{\beta }_u^0\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{w}_i^0{\sigma }_i^0}\right]\right\}\\ +\left(\delta {\beta }_u\right){\Sigma }_d^0{\rho }_{in}^0\exp \left[\left({\beta }_{\mathcal{l}}^0-{\beta }_u^0\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{w}_i^0{\sigma }_i^0}\right]-\left(\delta {\beta }_{\mathcal{l}}\right){\Sigma }_d^0{\rho }_{in}^0.\end{array}$ (16)

The indirect effect term, ${\left(\delta R_2\right)}^{ind}$, depends on the unknown variation $\delta \rho \left(t\right)$, of the state function $\rho \left(t\right)$, which is the solution of the following First-Level Forward Sensitivity System (1^st-LFSS) obtained by G-differentiating Equations (1) and (2) around the nominal parameter values:

$\frac{\text{d}\left[\delta \rho \left(t\right)\right]}{\text{d}t}+\delta \rho \left(t\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i^0{\sigma }_i^0}=-{\rho }^0\left(t\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}\left(n_i^0\delta {\sigma }_i+{\sigma }_i^0\delta n_i\right)}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{},0\le {\beta }_{\mathcal{l}}^0\le t\le {\beta }_u^0<\infty ,$ (17)

$\delta \rho \left({\beta }_{\mathcal{l}}^0\right)+\delta {\beta }_{\mathcal{l}}{\left\{\frac{\text{d}\rho \left(t\right)}{\text{d}t}\right\}}_{t={\beta }_{\mathcal{l}}^0}=\delta {\rho }_{in},\text{at}t={\beta }_{\mathcal{l}}^0.$ (18)

The need for performing the many large-scale computations for obtaining for all possible variations in the model and boundary parameters can be avoided by applying the 2^nd-CASAM presented in [1]. In order to apply the 2^nd-CASAM, the function $\delta \rho \left(t\right)$ is considered to be an element of a Hilbert space $H^{\left(1\right)}\left({\Omega }_t\right)$, ${\Omega }_t\triangleq \left({\beta }_{\mathcal{l}}^0,{\beta }_u^0\right)$, endowed with the following inner product, denoted as $\langle {\rho }_1\left(t\right),{\rho }_2\left(t\right)\rangle$, between two (square-integrable) functions ${\rho }_1\left(t\right)\in H^{\left(1\right)}\left({\Omega }_t\right)$ and ${\rho }_2\left(t\right)\in H^{\left(1\right)}\left({\Omega }_t\right)$ :

$\langle {\rho }_1\left(t\right),{\rho }_2\left(t\right)\rangle \triangleq {\displaystyle \underset{{\beta }_{\mathcal{l}}^0}{\overset{{\beta }_u^0}{\int }}{\rho }_1\left(t\right){\rho }_2\left(t\right)}\text{d}t.$ (19)

The indirect effect term, ${\left(\delta R_2\right)}^{ind}$, is computed by applying the general 2^nd-CASAM presented in [1] to Equation (17), which commences by using the definition of the inner product provided in Equation (19) to form the inner product of Equation (17) with a square-integrable function ${\theta }^{\left(1\right)}\left(t\right)\in H^{\left(1\right)}\left({\Omega }_t\right)$ and integrating the left-side of the resulting equation by parts once, so as to transfer the differential operation from $\delta \rho \left(t\right)$ onto ${\theta }^{\left(1\right)}\left(t\right)$, to obtain:

$\begin{array}{l}{\displaystyle \underset{{\beta }_{\mathcal{l}}^0}{\overset{{\beta }_u^0}{\int }}{\theta }^{\left(1\right)}\left(t\right)\left[\frac{\text{d}\left[\delta \rho \left(t\right)\right]}{\text{d}t}+\delta \rho \left(t\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{w}_i^0{\sigma }_i^0}\right]}\text{d}t\\ ={\displaystyle \underset{{\beta }_{\mathcal{l}}^0}{\overset{{\beta }_u^0}{\int }}\delta \rho \left(t\right)\left[-\frac{\text{d}{\theta }^{\left(1\right)}\left(t\right)}{\text{d}t}+{\theta }^{\left(1\right)}\left(t\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{w}_i^0{\sigma }_i^0}\right]}\text{d}t+{\theta }^{\left(1\right)}\left({\beta }_u^0\right)\delta \rho \left({\beta }_u^0\right)-{\theta }^{\left(1\right)}\left({\beta }_{\mathcal{l}}^0\right)\delta \rho \left({\beta }_{\mathcal{l}}^0\right).\end{array}$ (20)

The following sequence of operations is performed next using Equation (20):

1) Require that the first term on the right-side of Equation (20) be identical with the indirect effect ${\left(\delta R_2\right)}^{ind}$ defined in Equation (14).

2) Use the right-side of Equation (17) to replace the term multiplying ${\theta }^{\left(1\right)}\left(t\right)$ on the left-side of Equation (20).

3) Eliminate the unknown quantity $\delta \rho \left({\beta }_u^0\right)$ on the right-side of Equation (20) by imposing the condition ${\theta }^{\left(1\right)}\left({\beta }_u^0\right)=0$.

4) Insert the boundary condition provided in Equation (18) into the resulting expression.

The result of the above sequence of operations is the following expression for the indirect-effect term ${\left(\delta R_2\right)}^{ind}$ :

$\begin{array}{l}{\left(\delta R_2\right)}^{ind}=-{\displaystyle \underset{i=1}{\overset{N}{\sum }}\left(n_i^0\delta {\sigma }_i+{\sigma }_i^0\delta n_i\right)}{\displaystyle \underset{{\beta }_{\mathcal{l}}^0}{\overset{{\beta }_u^0}{\int }}{\theta }^{\left(1\right)}\left(t\right){\rho }^0\left(t\right)}\text{d}t\\ +{\theta }^{\left(1\right)}\left({\beta }_{\mathcal{l}}^0\right)\left[\delta {\rho }_{in}-\delta \beta {\left\{\frac{\partial \rho \left(t\right)}{\partial t}\right\}}_{t={\beta }_{\mathcal{l}}^0}\right],\end{array}$ (21)

where the first-level adjoint function ${\theta }^{\left(1\right)}\left(t\right)$ appearing in Equation (21) is the solution of the following First-Level Adjoint Sensitivity System (1^st-LASS):

$-\frac{\text{d}{\theta }^{\left(1\right)}\left(t\right)}{\text{d}t}+{\theta }^{\left(1\right)}\left(t\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i^0{\sigma }_i^0={\Sigma }_d^0},0\le {\beta }_{\mathcal{l}}^0\le t\le {\beta }_u^0<\infty ,$ (22)

${\theta }^{\left(1\right)}\left({\beta }_u^0\right)=0$. (23)

Since the 1^st-LASS does not depend on any parameter or boundary variation, a single large-scale computation for obtaining the adjoint function ${\theta }^{\left(1\right)}\left(t\right)$ suffices for computing exactly and efficiently, using quadrature formulas, all of the $2N+4$ sensitivities of the response $R_2\left(\rho ;\alpha ,\beta \right)$ with respect to all model and boundary parameters. Solving the 1^st-LASS defined by Equations (22) and (23) yields the following expression for the first-level adjoint function ${\theta }^{\left(1\right)}\left(t\right)$ :

${\theta }^{\left(1\right)}\left(t\right)=\frac{{\Sigma }_d^0}{{\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i^0{\sigma }_i^0}}\left\{1-\exp \left[\left(t-{\beta }_u^0\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i^0{\beta }_i^0}\right]\right\}.$ (24)

Adding Equations (21) and (15), and identifying the quantities that multiply the respective parameter variations yields the following expressions for the first-order sensitivities of $R_2\left(\rho ;\alpha ,\beta \right)$ in terms of the first-level adjoint function ${\theta }^{\left(1\right)}\left(t\right)$ :

${\left(\frac{\partial R_2}{\partial {\rho }_{in}}\right)}_{e^0}={\theta }^{\left(1\right)}\left({\beta }_{\mathcal{l}}^0\right),$ (25)

${\left(\frac{\partial R_2}{\partial {\sigma }_i}\right)}_{e^0}=-n_i^0{\displaystyle \underset{{\beta }_{\mathcal{l}}^0}{\overset{{\beta }_u^0}{\int }}{\theta }^{\left(1\right)}\left(t\right){\rho }^0\left(t\right)\text{d}t},i=1,\cdots ,N,$ (26)

${\left(\frac{\partial R_2}{\partial n_i}\right)}_{e^0}=-{\sigma }_i^0{\displaystyle \underset{{\beta }_{\mathcal{l}}^0}{\overset{{\beta }_u^0}{\int }}{\theta }^{\left(1\right)}\left(t\right){\rho }^0\left(t\right)}\text{d}t,i=1,\cdots ,N,$ (27)

${\left(\frac{\partial R_2}{\partial {\Sigma }_d}\right)}_{e^0}={\displaystyle \underset{{\beta }_{\mathcal{l}}^0}{\overset{{\beta }_u^0}{\int }}\rho \left(t\right)}\text{d}t,$ (28)

${\left(\frac{\partial R_2}{\partial {\beta }_u}\right)}_{e^0}={\Sigma }_d^0{\left\{\rho \left(t\right)\right\}}_{t={\beta }_u^0},$ (29)

${\left(\frac{\partial R_2}{\partial {\beta }_{\mathcal{l}}}\right)}_{e^0}=-{\Sigma }_d^0{\left\{\rho \left(t\right)\right\}}_{t={\beta }_{\mathcal{l}}^0}-{\theta }^{\left(1\right)}\left({\beta }_{\mathcal{l}}^0\right){\left\{\frac{\partial \rho \left(t\right)}{\partial t}\right\}}_{t={\beta }_{\mathcal{l}}^0}.$ (30)

The results in Equations (26) and (27) indicate that the relative sensitivities of the response to the parameters ${\sigma }_i$, $\left[\partial R_2/\partial {\sigma }_i\right]\left[{\sigma }_i/R_2\right]$, are identical to the relative sensitivities $\left[\partial R_2/\partial n_i\right]\left[n_i/R_2\right]$ of the response to the parameters $n_i$, for each $i=1,\cdots ,N$, for all of these 2N model parameters, since

$\frac{\partial R_2}{\partial {\sigma }_i}\frac{{\sigma }_i}{R_2}=\frac{{\sigma }_i{n}_i}{R_2}{\displaystyle \underset{0}{\overset{t_f}{\int }}\left[-{\theta }^{\left(1\right)}\left(t\right)\rho \left(t\right)\right]}\text{d}t\equiv \frac{\partial R_2}{\partial n_i}\frac{n_i}{R_2},i=1,\cdots ,N.$ (31)

As indicated in Equations (28) and (29), the first-order sensitivities of $R_2\left(\rho ;\alpha ,\beta \right)$ with respect to the model/detector parameter ${\Sigma }_d$ and the boundary parameter ${\beta }_u$ stem exclusively from the direct effect term, ${\left(\delta R_2\right)}^{dir}$, and can therefore be computed directly using the solution $\rho \left(t\right)$, provided in Equation (6), of the forward model [cf., Equations (1) and (2)]. Thus, inserting the expression of $\rho \left(t\right)$ and performing the operations indicated in Equations (28) and (29), respectively, yields the following expressions:

${\left(\frac{\partial R_2}{\partial {\Sigma }_d}\right)}_{e^0}=\frac{{\rho }_{in}^0}{{\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i^0{\sigma }_i^0}}\left\{1-\exp \left[\left({\beta }_{\mathcal{l}}^0-{\beta }_u^0\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i^0{\sigma }_i^0}\right]\right\},$ (32)

${\left(\frac{\partial R_2}{\partial {\beta }_u}\right)}_{e^0}={\Sigma }_d^0{\rho }_{in}^0\exp \left[\left({\beta }_{\mathcal{l}}^0-{\beta }_u^0\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i^0{\sigma }_i^0}\right].$ (33)

On the other hand, as indicated in Equations (25)-(27), the first-order sensitivities of $R_2\left(\rho ;\alpha ,\beta \right)$ with respect to the model parameters ${\rho }_{in}$, ${\sigma }_i$, and $n_i$ stem exclusively from the indirect effect term ${\left(\delta R_2\right)}^{ind}$ and can therefore be computed after having obtained the first-level adjoint function ${\theta }^{\left(1\right)}\left(t\right)$ by solving the 1^st-LASS, namely Equations (22) and (23). Finally, the sensitivity of $R_2\left(\rho ;\alpha ,\beta \right)$ with respect to the boundary parameter ${\beta }_{\mathcal{l}}$ stems from contributions arising from both the direct and the indirect effect terms, as indicated by Equation (30).

Using the result for the first-level adjoint function ${\theta }^{\left(1\right)}\left(t\right)$ obtained in Equation (24) in Equations (25)-(27) and performing the respective operations yields the following results:

${\left(\frac{\partial R_2}{\partial {\rho }_{in}}\right)}_{e^0}=\frac{{\Sigma }_d^0}{{\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i^0{\sigma }_i^0}}\left\{1-\exp \left[\left({\beta }_{\mathcal{l}}^0-{\beta }_u^0\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i^0{\sigma }_i^0}\right]\right\},$ (34)

${\left(\frac{\partial R_2}{\partial {\sigma }_i}\right)}_{e^0}=-\frac{{\rho }_{in}^0{\Sigma }_d^0{n}_i^0}{{\left({\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i^0{\sigma }_i^0}\right)}^2}\left\{1+\left[\left({\beta }_{\mathcal{l}}^0-{\beta }_u^0\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i^0{\sigma }_i^0}-1\right]\exp \left[\left({\beta }_{\mathcal{l}}^0-{\beta }_u^0\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i^0{\sigma }_i^0}\right]\right\},$ (35)

${\left(\frac{\partial R_2}{\partial n_i}\right)}_{e^0}=-\frac{{\rho }_{in}^0{\Sigma }_d^0{\sigma }_i^0}{{\left({\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i^0{\sigma }_i^0}\right)}^2}\left\{1+\left[\left({\beta }_{\mathcal{l}}^0-{\beta }_u^0\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i^0{\sigma }_i^0}-1\right]\exp \left[\left({\beta }_{\mathcal{l}}^0-{\beta }_u^0\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i^0{\sigma }_i^0}\right]\right\}.$ (36)

Finally, the sensitivity of $R_2\left(\rho ;\alpha ,\beta \right)$ with respect to the boundary parameter ${\beta }_{\mathcal{l}}$ is computed by using in Equation (30) the expression of first-level adjoint function ${\theta }^{\left(1\right)}\left(t\right)$ obtained in Equation (24) together with the expression of ${\rho }^0\left(t\right)$ provided in Equation (9), to obtain the following result:

${\left(\frac{\partial R_2}{\partial {\beta }_{\mathcal{l}}}\right)}_{e^0}=-{\Sigma }_d^0{\rho }_{in}^0\exp \left[\left({\beta }_{\mathcal{l}}^0-{\beta }_u^0\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i^0{\sigma }_i^0}\right].$ (37)

Notably for the result obtained in Equation (37), the contribution from the indirect effect term involving ${\theta }^{\left(1\right)}\left(t\right)$ partially offsets the contribution from the direct effect term involving the expression of $\rho \left(t={\beta }_{\mathcal{l}}^0\right)$.

3.2. Computing the 2^nd-Order Sensitivities of the Response $R_2\left(\rho ;\alpha ,\beta \right)$ Using Second-Level Adjoint Sensitivity Systems (2nd-LASS)

As indicated by the general methodology presented in [1], the 2^nd-order sensitivities of $R_2\left(\rho ;\alpha ,\beta \right)$ are obtained by computing successively the 1^st-order G-differentials of the 1^st-order sensitivities obtained in Equations (25)-(30). The detailed steps will be illustrated by computing the 2^nd-order sensitivities corresponding to defined in Equation (26), which is representative of the steps that would be repeated for the computations of the other 2^nd-order sensitivities of $R_2\left(\rho ;\alpha ,\beta \right)$, which would be derived from the 1^st-order G-differentials of the expressions presented in Equations (25) and (27)-(30). Omitting, for notational simplicity, the superscript “zero” (which denotes “nominal values” in this work), the 1^st-order G-differential of the expression provided in Equation (26) is obtained as follows:

$\begin{array}{l}{\left\{\delta \left(\frac{\partial R_2}{\partial {\sigma }_i}\right)\right\}}_{e^0}\\ \triangleq -\frac{\text{d}}{\text{d}\epsilon }{\left\{\left(n_i+\epsilon \delta n_i\right){\displaystyle \underset{{\beta }_{\mathcal{l}}^0+\epsilon \delta {\beta }_{\mathcal{l}}}{\overset{{\beta }_u^0+\epsilon \delta {\beta }_u}{\int }}\left[{\theta }^{\left(1\right)}\left(t\right)+\epsilon \delta {\theta }^{\left(1\right)}\left(t\right)\right]\left[\rho \left(t\right)+\epsilon \delta \rho \left(t\right)\right]\text{d}t}\right\}}_{\left(e=e^0,\epsilon =0\right)}\\ ={\left\{{\left[\delta \left(\frac{\partial R_2}{\partial {\sigma }_i}\right)\right]}_{dir}\right\}}_{e^0}+{\left\{{\left[\delta \left(\frac{\partial R_2}{\partial {\sigma }_i}\right)\right]}_{indir}\right\}}_{e^0},\end{array}$ (38)

where

$\begin{array}{l}{\left\{{\left[\delta \left(\frac{\partial R_2}{\partial {\sigma }_i}\right)\right]}_{dir}\right\}}_{e^0}\\ \triangleq -\left(\delta n_i\right){\displaystyle \underset{{\beta }_{\mathcal{l}}^0}{\overset{{\beta }_u^0}{\int }}{\theta }^{\left(1\right)}\left(t\right)\rho \left(t\right)\text{d}t}+\left(\delta {\beta }_{\mathcal{l}}\right)n_i{\theta }^{\left(1\right)}\left({\beta }_{\mathcal{l}}^0\right)\rho \left({\beta }_{\mathcal{l}}^0\right),\text{for}i=1,\cdots ,N;\end{array}$ (39)

$\begin{array}{l}{\left\{{\left[\delta \left(\frac{\partial R_2}{\partial {\sigma }_i}\right)\right]}_{indir}\right\}}_{e^0}\\ \triangleq {\left\{-n_i{\displaystyle \underset{{\beta }_{\mathcal{l}}^0}{\overset{{\beta }_u^0}{\int }}\delta {\theta }^{\left(1\right)}\left(t\right)\rho \left(t\right)\text{d}t}-n_i{\displaystyle \underset{{\beta }_{\mathcal{l}}^0}{\overset{{\beta }_u^0}{\int }}{\theta }^{\left(1\right)}\left(t\right)\delta \rho \left(t\right)\text{d}t}\right\}}_{e^0},\text{}\text{}\text{}\text{}\text{for}i=1,\cdots ,N.\end{array}$ (40)

The direct-effect term defined by Equation (39) can be computed immediately, since the 1^st-level adjoint function ${\theta }^{\left(1\right)}\left(t\right)$ and the forward function $\rho \left(t\right)$ are known. However, the indirect-effect term defined by Equation (40) contains the variation $\delta {\theta }^{\left(1\right)}\left(t\right)$ in the 1^st-level adjoint function and, respectively, the variation $\delta \rho \left(t\right)$ in the forward function, both of which depend on parameter variations and neither of which is immediately available. The variation $\delta {\theta }^{\left(1\right)}\left(t\right)$ of the 1^st-level adjoint function ${\theta }^{\left(1\right)}\left(t\right)$ is related to the parameter variations through the G-differential of the 1^st-LASS, which is derived by applying the definition of the G-differential to Equations (22) and (23) to obtain the following equations, evaluated at the nominal parameter values:

$\begin{array}{l}-\frac{\text{d}\left[\delta {\theta }^{\left(1\right)}\left(t\right)\right]}{\text{d}t}+\delta {\theta }^{\left(1\right)}\left(t\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}\\ =\left(\delta {\Sigma }_d\right)-{\theta }^{\left(1\right)}\left(t\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}\left(n_i\delta {\sigma }_i+{\sigma }_i\delta n_i\right)},0\le {\beta }_{\mathcal{l}}^0\le t\le {\beta }_u^0<\infty ,\end{array}$ (41)

$\begin{array}{l}{\left\{\delta {\theta }^{\left(1\right)}\left(t\right)+\left(\delta \beta {}_u\right)\frac{\text{d}{\theta }^{\left(1\right)}\left(t\right)}{\text{d}t}\right\}}_{t={\beta }_u^0}\\ =\delta {\theta }^{\left(1\right)}\left({\beta }_u^0\right)=0.\end{array}$ (42)

Since Equations (41) and (42) depend on the parameter variations, solving them would be just as impractical computationally as solving the 1^st-LFSS. Therefore, the need for solving these equations will be circumvented by expressing the indirect-effect term defined in Equation (40) in an alternative way so as to eliminate the appearance of $\delta {\theta }^{\left(1\right)}\left(t\right)$ and $\delta \rho \left(t\right)$. For this purpose, we introduce another Hilbert space, denoted as $H^{\left(2\right)}\left({\Omega }_t\right),{\Omega }_t\triangleq \left({\beta }_{\mathcal{l}}^0,{\beta }_u^0\right)$, which comprises, as elements, two-component vectors of the form

${\psi }_i^{\left(2\right)}\left(t\right)\triangleq \left[{\psi }_{i1}^{\left(2\right)}\left(t\right),{\psi }_{i2}^{\left(2\right)}\left(t\right)\right]$, with square-integrable functions ${\psi }_{ij}^{\left(2\right)}\left(t\right),j=1,2$. The inner product between two elements ${\psi }_i^{\left(2\right)}\left(t\right)\triangleq \left[{\psi }_{i1}^{\left(2\right)}\left(t\right),{\psi }_{i2}^{\left(2\right)}\left(t\right)\right]\in H^{\left(2\right)}(\; \Omega \; t\; )$

and ${\phi }_i^{\left(2\right)}\left(t\right)\triangleq \left[{\phi }_{i1}^{\left(2\right)}\left(t\right),{\phi }_{i2}^{\left(2\right)}\left(t\right)\right]\in H^{\left(2\right)}\left({\Omega }_t\right)$ in the Hilbert space $H^{\left(2\right)}\left({\Omega }_t\right)$ will be denoted as ${\langle {\psi }_i^{\left(2\right)}\left(t\right),{\phi }_i^{\left(2\right)}\left(t\right)\rangle }_2$ and is defined as follows:

${\langle {\psi }_i^{\left(2\right)}\left(t\right),{\phi }_i^{\left(2\right)}\left(t\right)\rangle }_2\triangleq {\displaystyle \underset{j=1}{\overset{2}{\sum }}{\displaystyle \underset{{\beta }_{\mathcal{l}}^0}{\overset{{\beta }_u^0}{\int }}{\psi }_{ij}^{\left(2\right)}\left(t\right){\phi }_{ij}^{\left(2\right)}\left(t\right)\text{d}t}}.$ (43)

Using the definition given in Equation (43), construct the inner product of Equations (41) and (17) with a square integrable two-component function ${\theta }_1^{\left(2\right)}\left(t\right)\triangleq \left[{\theta }_{11}^{\left(2\right)}\left(t\right),{\theta }_{12}^{\left(2\right)}\left(t\right)\right]\in H^{\left(2\right)}\left({\Omega }_t\right)$ to obtain the following relation:

$\begin{array}{l}{\displaystyle \underset{{\beta }_{\mathcal{l}}^0}{\overset{{\beta }_u^0}{\int }}\left[{\theta }_{11}^{\left(2\right)}\left(t\right),{\theta }_{12}^{\left(2\right)}\left(t\right)\right]}\left(\begin{array}{cc}-\frac{\text{d}}{\text{d}t}+{\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}& 0\\ 0& \frac{\text{d}}{\text{d}t}+{\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}\end{array}\right)\left(\begin{array}{c}\delta {\theta }^{\left(1\right)}\left(t\right)\\ \delta \rho \left(t\right)\end{array}\right)\text{d}t\\ ={\displaystyle \underset{{\beta }_{\mathcal{l}}^0}{\overset{{\beta }_u^0}{\int }}\left[{\theta }_{11}^{\left(2\right)}\left(t\right),{\theta }_{12}^{\left(2\right)}\left(t\right)\right]}\left(\begin{array}{c}\left(\delta {\Sigma }_d\right)-{\theta }^{\left(1\right)}\left(t\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}\left(n_i\delta {\sigma }_i+{\sigma }_i\delta n_i\right)}\\ -\rho \left(t\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}\left(n_i\delta {\sigma }_i+{\sigma }_i\delta n_i\right)}\end{array}\right)\text{d}t.\end{array}$ (44)

Integrating by parts the left-side of Equation (44) so as to transfer the differential operations on $\delta {\theta }^{\left(1\right)}\left(t\right)$ and $\delta \rho \left(t\right)$ to differential operations on ${\theta }_{11}^{\left(2\right)}\left(t\right)$ and ${\theta }_{12}^{\left(2\right)}\left(t\right)$ yields the following result:

$\begin{array}{l}{\displaystyle \underset{{\beta }_{\mathcal{l}}^0}{\overset{{\beta }_u^0}{\int }}{\theta }_{11}^{\left(2\right)}\left(t\right)\left[-\frac{\text{d}}{\text{d}t}\delta {\theta }^{\left(1\right)}\left(t\right)+\delta {\theta }^{\left(1\right)}\left(t\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}\right]\text{d}t}\\ +{\displaystyle \underset{{\beta }_{\mathcal{l}}^0}{\overset{{\beta }_u^0}{\int }}{\theta }_{12}^{\left(2\right)}\left(t\right)\left[\frac{\text{d}}{\text{d}t}\delta \rho \left(t\right)+\delta \rho \left(t\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}\right]\text{d}t}\\ =-{\theta }_{11}^{\left(2\right)}\left({\beta }_u^0\right)\delta {\theta }^{\left(1\right)}\left({\beta }_u^0\right)+{\theta }_{11}^{\left(2\right)}\left({\beta }_{\mathcal{l}}^0\right)\delta {\theta }^{\left(1\right)}\left({\beta }_{\mathcal{l}}^0\right)+{\theta }_{12}^{\left(2\right)}\left({\beta }_u^0\right)\delta \rho \left({\beta }_u^0\right)\\ -{\theta }_{12}^{\left(2\right)}\left({\beta }_{\mathcal{l}}^0\right)\delta \rho \left({\beta }_{\mathcal{l}}^0\right)+{\displaystyle \underset{{\beta }_{\mathcal{l}}^0}{\overset{{\beta }_u^0}{\int }}\delta {\theta }^{\left(1\right)}\left(t\right)\left[\frac{\text{d}{\theta }_{11}^{\left(2\right)}\left(t\right)}{\text{d}t}+{\theta }_{11}^{\left(2\right)}\left(t\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}\right]\text{d}t}\\ +{\displaystyle \underset{{\beta }_{\mathcal{l}}^0}{\overset{{\beta }_u^0}{\int }}\delta \rho \left(t\right)\left[-\frac{\text{d}{\theta }_{12}^{\left(2\right)}\left(t\right)}{\text{d}t}+{\theta }_{12}^{\left(2\right)}\left(t\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}\right]\text{d}t}.\end{array}$ (45)

The last two terms on the right-side of Equation (45) are now required to represent the indirect-effect term defined in Equation (40) by imposing the following relations:

$\frac{\text{d}{\theta }_{11}^{\left(2\right)}\left(t\right)}{\text{d}t}+{\theta }_{11}^{\left(2\right)}\left(t\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}=-n_i\rho \left(t\right)=-n_i{\rho }_{in}\exp \left[\left({\beta }_{\mathcal{l}}-t\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}\right],$ (46)

$\begin{array}{l}-\frac{\text{d}{\theta }_{12}^{\left(2\right)}\left(t\right)}{\text{d}t}+{\theta }_{12}^{\left(2\right)}\left(t\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}\\ =-n_i{\theta }^{\left(1\right)}\left(t\right)=-\frac{n_i{\Sigma }_d}{{\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}}\left\{1-\exp \left[\left(t-{\beta }_u\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}\right]\right\}.\end{array}$ (47)

The boundary conditions for Equations (46) and (47) are established by requiring the contributions involving the unknown quantities $\delta {\theta }^{\left(1\right)}\left({\beta }_{\mathcal{l}}\right)$ and $\delta \rho \left({\beta }_u\right)$ in Equation (45) to vanish, which can be accomplished by imposing the following conditions:

${\theta }_{11}^{\left(2\right)}\left({\beta }_{\mathcal{l}}\right)=0,{\theta }_{12}^{\left(2\right)}\left({\beta }_u\right)=0.$ (48)

The system of equations comprising Equations (46)-(48), which is independent of parameter variations, constitutes the 2^nd-Level Adjoint Sensitivity System (2^nd-LASS) for the two-component 2^nd-level adjoint function ${\theta }_1^{\left(2\right)}\left(t\right)\triangleq \left[{\theta }_{11}^{\left(2\right)}\left(t\right),{\theta }_{12}^{\left(2\right)}\left(t\right)\right]\in H^{\left(2\right)}\left({\Omega }_t\right)$.

Replacing the left-side of Equation (45) by the right-side of Equation (44) and taking into account Equations (46)-(48) yields the following expression for the indirect effect term defined in Equation (40):

$\begin{array}{l}{\left[\delta \left(\frac{\partial R_2}{\partial {\sigma }_i}\right)\right]}_{indir}\\ ={\displaystyle \underset{{\beta }_{\mathcal{l}}^0}{\overset{{\beta }_u^0}{\int }}\text{d}t\left[{\theta }_{11}^{\left(2\right)}\left(t\right),{\theta }_{12}^{\left(2\right)}\left(t\right)\right]}\left(\begin{array}{c}\left(\delta {\Sigma }_d\right)-{\theta }^{\left(1\right)}\left(t\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}\left(n_i\delta {\sigma }_i+{\sigma }_i\delta n_i\right)}\\ -\rho \left(t\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}\left(n_i\delta {\sigma }_i+{\sigma }_i\delta n_i\right)}\end{array}\right)\\ +{\theta }_{11}^{\left(2\right)}\left({\beta }_u^0\right)\delta {\theta }^{\left(1\right)}\left({\beta }_u^0\right)+{\theta }_{12}^{\left(2\right)}\left({\beta }_{\mathcal{l}}^0\right)\delta \rho \left({\beta }_{\mathcal{l}}^0\right),\text{}\text{}\text{}\text{}\text{for}i=1,\cdots ,N.\end{array}$ (49)

Using the conditions given in Equations (18) and (42) in the last terms on the right side of Equation (49) yields the following expression for the indirect-effect term:

$\begin{array}{l}{\left[\delta \left(\frac{\partial R_2}{\partial {\sigma }_i}\right)\right]}_{indir}\\ ={\displaystyle \underset{{\beta }_{\mathcal{l}}^0}{\overset{{\beta }_u^0}{\int }}\left[{\theta }_{11}^{\left(2\right)}\left(t\right),{\theta }_{12}^{\left(2\right)}\left(t\right)\right]}\left(\begin{array}{c}\left(\delta {\Sigma }_d\right)-{\theta }^{\left(1\right)}\left(t\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}\left(n_i\delta {\sigma }_i+{\sigma }_i\delta n_i\right)}\\ -\rho \left(t\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}\left(n_i\delta {\sigma }_i+{\sigma }_i\delta n_i\right)}\end{array}\right)\text{d}t\\ +{\theta }_{12}^{\left(2\right)}\left({\beta }_{\mathcal{l}}\right)\left\{\delta {\rho }_{in}-\delta {\beta }_{\mathcal{l}}{\left[\frac{\text{d}\rho \left(t\right)}{\text{d}t}\right]}_{t={\beta }_{\mathcal{l}}^0}\right\},\text{}\text{}\text{}\text{}\text{for}i=1,\cdots ,N.\end{array}$ (50)

Adding the direct-effect term defined in Equation (39) to Equation (50) and identifying in the resulting expression the coefficients multiplying the variations $\delta {\sigma }_i$, $\delta n_i$, $\delta {\rho }_{in}$, $\delta t_d$, $\delta {\beta }_{\mathcal{l}}$ and $\delta {\beta }_u$ yields the following expression for the respective 2^nd-order sensitivities of the response $R_2\left(\rho ;\alpha ,\beta \right)$ :

$\frac{{\partial }^2{R}_2}{\partial {\sigma }_j\partial {\sigma }_i}=-n_j{\displaystyle \underset{{\beta }_{\mathcal{l}}^0}{\overset{{\beta }_u^0}{\int }}{\theta }_{11}^{\left(2\right)}\left(t\right){\theta }^{\left(1\right)}\left(t\right)\text{d}t}-n_j{\displaystyle \underset{{\beta }_{\mathcal{l}}^0}{\overset{{\beta }_u^0}{\int }}{\theta }_{12}^{\left(2\right)}\left(t\right)\rho \left(t\right)}\text{d}t;i,j=1,\cdots ,N;$ (51)

$\begin{array}{l}\frac{{\partial }^2{R}_2}{\partial n_j\partial {\sigma }_i}=-{\sigma }_j{\displaystyle \underset{{\beta }_{\mathcal{l}}^0}{\overset{{\beta }_u^0}{\int }}{\theta }_{11}^{\left(2\right)}\left(t\right){\theta }^{\left(1\right)}\left(t\right)\text{d}t}-{\sigma }_j{\displaystyle \underset{{\beta }_{\mathcal{l}}^0}{\overset{{\beta }_u^0}{\int }}{\theta }_{12}^{\left(2\right)}\left(t\right)\rho \left(t\right)}\text{d}t\\ -{\delta }_{ij}{\displaystyle \underset{{\beta }_{\mathcal{l}}^0}{\overset{{\beta }_u^0}{\int }}{\theta }^{\left(1\right)}\left(t\right)\rho \left(t\right)}\text{d}t;i,j=1,\cdots ,N;{\delta }_{ij}=\{\begin{array}{l}1,i=j;\\ 0,i\ne j;\end{array}\end{array}$ (52)

$\frac{{\partial }^2{R}_2}{\partial {\rho }_{in}\partial {\sigma }_i}={\theta }_{12}^{\left(2\right)}\left({\beta }_{\mathcal{l}}\right);i=1,\cdots ,N;$ (53)

$\frac{{\partial }^2{R}_2}{\partial {\Sigma }_d\partial {\sigma }_i}={\displaystyle \underset{{\beta }_{\mathcal{l}}^0}{\overset{{\beta }_u^0}{\int }}{\theta }_{11}^{\left(2\right)}\left(t\right)}\text{d}t;i=1,\cdots ,N;$ (54)

$\frac{{\partial }^2{R}_2}{\partial {\beta }_u\partial {\sigma }_i}=-{\theta }_{11}^{\left(2\right)}\left({\beta }_u\right){\left[\frac{\text{d}{\theta }^{\left(1\right)}\left(t\right)}{\text{d}t}\right]}_{t={\beta }_u^0};i=1,\cdots ,N;$ (55)

$\frac{{\partial }^2{R}_2}{\partial {\beta }_{\mathcal{l}}\partial {\sigma }_i}=n_i{\theta }^{\left(1\right)}\left({\beta }_{\mathcal{l}}\right)\rho \left({\beta }_{\mathcal{l}}\right)-{\theta }_{12}^{\left(2\right)}\left({\beta }_{\mathcal{l}}\right){\left[\frac{\text{d}\rho \left(t\right)}{\text{d}t}\right]}_{t={\beta }_{\mathcal{l}}^0};i=1,\cdots ,N.$ (56)

The 2^nd-order sensitivities shown in Equations (51)-(56) can be computed after having determined the 2^nd-level adjoint function ${\theta }_1^{\left(2\right)}\left(t\right)\triangleq \left[{\theta }_{11}^{\left(2\right)}\left(t\right),{\theta }_{12}^{\left(2\right)}\left(t\right)\right]$ by solving the 2^nd-LASS comprising Equations (46)-(48) using the nominal parameter values (the superscript “zero,” which indicates “nominal values,” has been omitted, for simplicity). Since the model parameters $n_i$ depend on the index $i=1,\cdots ,N$, it follows that the right-sides of Equations (46)-(48) also depend on this index. Strictly speaking, therefore, the 2^nd-level adjoint sensitivity function ${\theta }_1^{\left(2\right)}\left(t\right)\triangleq \left[{\theta }_{11}^{\left(2\right)}\left(t\right),{\theta }_{12}^{\left(2\right)}\left(t\right)\right]$ is a function of the index $i=1,\cdots ,N$. Hence, in the most unfavorable situation, the 2^nd-LASS, comprising Equations (46)-(48) would need to be solved numerically for each distinct value $n_i$, for a total of N-times. Even in such a “worse-case scenario,” however, only the right sides (i.e., “sources”) of Equations (46) and (47) would need to be modified, which is relatively easy to implement computationally. The left-sides of these equations remain unchanged since they are independent of the index $i=1,\cdots ,N$. Notably, the change of the dependent variables

${\theta }_{11}^{\left(2\right)}\left(t\right)=n_i{u}_{11}^{\left(2\right)}\left(t\right)=n_i{\phi }_{11}^{\left(2\right)}\left(t\right),{\theta }_{12}^{\left(2\right)}\left(t\right)=n_i{u}_{12}^{\left(2\right)}\left(t\right),$ (57)

would transform (46)-(48) into the following form:

$\frac{\text{d}{u}_{11}^{\left(2\right)}\left(t\right)}{\text{d}t}+u_{11}^{\left(2\right)}\left(t\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\beta }_i}=-\rho \left(t\right),$ (58)

$-\frac{\text{d}{u}_{12}^{\left(2\right)}\left(t\right)}{\text{d}t}+u_{12}^{\left(2\right)}\left(t\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\beta }_i}=-{\theta }^{\left(1\right)}\left(t\right),$ (59)

$u_{11}^{\left(2\right)}\left({\beta }_{\mathcal{l}}\right)=0,u_{12}^{\left(2\right)}\left({\beta }_u\right)=0.$ (60)

The above (alternative) 2^nd-LASS, comprising Equations (58)-(60) is independent of the index $i=1,\cdots ,N$, and would need to be solved (numerically or analytically) only once, to obtain the following expressions for the functions $u_{11}^{\left(2\right)}\left(t\right)$ and $u_{12}^{\left(2\right)}\left(t\right)$ :

$u_{11}^{\left(2\right)}\left(t\right)={\rho }_{in}\left({\beta }_{\mathcal{l}}-t\right)\exp \left[\left({\beta }_{\mathcal{l}}-t\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}\right]={\phi }_{11}^{\left(2\right)}\left(t\right),$ (61)

$u_{12}^{\left(2\right)}\left(t\right)=\frac{{\Sigma }_d}{{\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}}\left\{\left[{\beta }_u-t+{\left({\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}\right)}^{-1}\right]\exp \left[\left(t-{\beta }_u\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}\right]-{\left({\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}\right)}^{-1}\right\}.$ (62)

The components of the 2^nd-level adjoint function ${\theta }_1^{\left(2\right)}\left(t\right)$ are obtained by multiplying the functions $u_{11}^{\left(2\right)}\left(t\right)$ and $u_{12}^{\left(2\right)}\left(t\right)$ by the respective model parameters $n_i$, to obtain the following expressions for the components of the 2^nd-level adjoint function ${\theta }_1^{\left(2\right)}\left(t\right)\triangleq \left[{\theta }_{11}^{\left(2\right)}\left(t\right),{\theta }_{12}^{\left(2\right)}\left(t\right)\right]$ :

${\theta }_{11}^{\left(2\right)}\left(t\right)=n_i{\rho }_{in}\left({\beta }_{\mathcal{l}}-t\right)\exp \left[\left({\beta }_{\mathcal{l}}-t\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}\right],$ (63)

${\theta }_{12}^{\left(2\right)}\left(t\right)=-\frac{n_i{\Sigma }_d}{{\left({\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}\right)}^2}\left\{1+\left[\left(t-{\beta }_u\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}-1\right]\exp \left[\left(t-{\beta }_u\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}\right]\right\}.$ (64)

Using Equations (63) and (64) in Equations (51)-(56) and performing the respective operations yields the following results for the partial 2^nd-order sensitivities:

$\begin{array}{l}\frac{{\partial }^2{R}_2}{\partial {\sigma }_j\partial {\sigma }_i}\\ =\frac{2n_i{n}_j{\rho }_{in}{\Sigma }_d}{{\left({\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}\right)}^3}\left\{1+\left[\left({\beta }_{\mathcal{l}}-{\beta }_u\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}-1\right]\exp \left[\left({\beta }_{\mathcal{l}}-{\beta }_u\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}\right]\right\}\\ -\frac{n_i{n}_j\left({\beta }_{\mathcal{l}}-{\beta }_u\right){\rho }_{in}{\Sigma }_d}{{\left({\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}\right)}^2}\left[\left({\beta }_{\mathcal{l}}-{\beta }_u\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}-1\right]\exp \left[\left({\beta }_{\mathcal{l}}-{\beta }_u\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}\right]\\ -\frac{n_i{n}_j\left({\beta }_{\mathcal{l}}-{\beta }_u\right){\rho }_{in}{\Sigma }_d}{{\left({\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}\right)}^2}\exp \left[\left({\beta }_{\mathcal{l}}-{\beta }_u\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}\right],i,j=1,\cdots ,N;\end{array}$ (65)

$\begin{array}{l}\frac{{\partial }^2{R}_2}{\partial n_j\partial {\sigma }_i}\\ ={\rho }_{in}{\Sigma }_d{\left({\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}\right)}^{-3}\left[2n_i{\sigma }_j-{\delta }_{ij}{\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}\right]\\ \times \left\{1+\left[\left({\beta }_{\mathcal{l}}-{\beta }_u\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}-1\right]\exp \left[\left({\beta }_{\mathcal{l}}-{\beta }_u\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}\right]\right\}\\ -\frac{n_i{\sigma }_j\left({\beta }_{\mathcal{l}}-{\beta }_u\right){\rho }_{in}{\Sigma }_d}{{\left({\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}\right)}^2}\left[\left({\beta }_{\mathcal{l}}-{\beta }_u\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}-1\right]\exp \left[\left({\beta }_{\mathcal{l}}-{\beta }_u\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}\right]\\ -\frac{n_i{\sigma }_j\left({\beta }_{\mathcal{l}}-{\beta }_u\right){\rho }_{in}{\Sigma }_d}{{\left({\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}\right)}^2}\exp \left[\left({\beta }_{\mathcal{l}}-{\beta }_u\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}\right],i,j=1,\cdots ,N;\end{array}$ (66)

$\begin{array}{l}\frac{{\partial }^2{R}_2}{\partial {\rho }_{in}\partial {\sigma }_i}\\ =-\frac{n_i{\Sigma }_d\left\{1+\left[\left({\beta }_{\mathcal{l}}-{\beta }_u\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}-1\right]\exp \left[\left({\beta }_{\mathcal{l}}-{\beta }_u\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}\right]\right\}}{{\left({\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}\right)}^2},i=1,\cdots ,N;\end{array}$ (67)

$\begin{array}{l}\frac{{\partial }^2{R}_2}{\partial {\Sigma }_d\partial {\sigma }_i}\\ =-\frac{n_i{\rho }_{in}\left\{1+\left[\left({\beta }_{\mathcal{l}}-{\beta }_u\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}-1\right]\exp \left[\left({\beta }_{\mathcal{l}}-{\beta }_u\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}\right]\right\}}{{\left({\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}\right)}^2},i=1,\cdots ,N;\end{array}$ (68)

$\frac{{\partial }^2{R}_2}{\partial {\beta }_u\partial {\sigma }_i}=n_i{\rho }_{in}{\Sigma }_d\left({\beta }_{\mathcal{l}}-{\beta }_u\right)\exp \left[\left({\beta }_{\mathcal{l}}-{\beta }_u\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}\right];i=1,\cdots ,N;$ (69)

$\frac{{\partial }^2{R}_2}{\partial {\beta }_{\mathcal{l}}\partial {\sigma }_i}=-n_i{\rho }_{in}{\Sigma }_d\left({\beta }_{\mathcal{l}}-{\beta }_u\right)\exp \left[\left({\beta }_{\mathcal{l}}-{\beta }_u\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{n}_i{\sigma }_i}\right],i=1,\cdots ,N.$ (70)

As before, the right-sides of expressions shown in Equations (65)-(70) are to be evaluated at the nominal values for the parameters, but the superscript “zero,” which indicates “nominal values,” has been omitted, for notational simplicity.

The expressions of the remaining 2^nd-order sensitivities of $R_2\left(\rho ;\alpha ,\beta \right)$ can be derived following the same procedure as outlined above, and as also outlined in [2] since these derivations would not illustrate any new concepts, they will not be presented in this work.

4. Conclusions

The results obtained by applying the general 2^nd-CASAM presented in [1] to the paradigm evolution/transmission benchmark analyzed in this work indicate the following major characteristics of this powerful methodology for computing exactly and efficiently the 1^st- and 2^nd-order sensitivities of model responses with respect to model and boundary parameters:

1) For a model comprising $N_{\alpha }$ distinct but uncertain model parameters and $N_{\beta }$ distinct but uncertain distinct boundary parameters, a single adjoint computation, to solve the 1^st-LASS, is necessary for computing exactly all of the $N_{\alpha }+N_{\beta }$ 1^st-order response sensitivities to the model’s uncertain parameters and boundaries. In contradistinction, $N_{\alpha }+N_{\beta }$ computations would be needed to obtain approximate values for the 1^st-order sensitivities if forward and/or statistical methods were used.

2) By considering each 1^st-order sensitivity as a response and developing a corresponding 2^nd-level adjoint sensitivity system (2^nd-LASS) for computing the respective 2^nd-order sensitivities, the application of the 2^nd-CASAM yields exact expressions/values for all (mixed and unmixed) 2^nd-order sensitivities. For each 1^st-order sensitivity, the solution of each of the 2^nd-LASS is at most a two-component (vector) 2^nd-level adjoint sensitivity function of the form ${\psi }_i^{\left(2\right)}\left(t\right)\triangleq \left[{\psi }_{i1}^{\left(2\right)}\left(t\right),{\psi }_{i2}^{\left(2\right)}\left(t\right)\right],i=1,\cdots ,N_{\alpha }+N_{\beta }$. Hence, although these 2^nd-level adjoint sensitivity functions are independent of any parameter variations, the 2^nd-LASS would need to be solved for $N_{\alpha }+N_{\beta }$ distinct right-sides (i.e., “source terms”), in the most unfavorable situation. Even in this “worse-case scenario,” only the right sides of the 2^nd-LASS would need to be modified in computational computer codes, which is a relatively easy programming task. The left-sides of the 2^nd-LASS equations (which contain differential operators, and which would therefore involve “solvers” that would be much more difficult to modify computationally) remain unchanged.

3) In many practical situations, it is possible to reduce drastically the number of computations involving the 2^nd-LASS. Occasionally, the solutions of some of the 2^nd-LASS can be written down by inspection, without actually solving the corresponding 2^nd-LASS. For the paradigm evolution problem analyzed in this chapter, for example, the 2^nd-LASS would need to be solved only four times, to compute the 2^nd-level adjoint functions ${\psi }_j^{\left(2\right)}\left(t\right)\triangleq \left[{\psi }_{j1}^{\left(2\right)}\left(t\right),{\psi }_{j2}^{\left(2\right)}\left(t\right)\right],j=1,2,3,4$, which would suffice for obtaining all of the $\left(N_{\alpha }+N_{\beta }\right)\left(N_{\alpha }+N_{\beta }+1\right)/2$ distinct 2^nd-order response sensitivities. Of course, such a very large reduction in the number of large-scale computation cannot be expected in every practical problem, but in most cases, the number of computations required for computing the complete set of 2^nd-order response sensitivities if far less than the number, $\left(N_{\alpha }+N_{\beta }\right)$, of model parameters.

4) The specific 2^nd-order sensitivities of interest can be selected “a priori,” based on the magnitude/importance of the 1^st-order sensitivities.

5) As has been generally shown by Cacuci [1] [2] [3], the mixed 2^nd-order sensitivities are obtained twice, stemming from distinct 2^nd-LASS. This fact enables the 2^nd-CASAM to provide an inherent independent verification of the correctness and accuracy of the 2^nd-level adjoint sensitivity functions that are used to compute the respective mixed 2^nd-order sensitivities.

6) The un-mixed 2^nd-order sensitivities of the form ${\partial }^2\rho \left(t_1\right)/\partial {\alpha }_i^2$ are obtained only once. Therefore, they can be independently verified either by solving the 2^nd-LFSS, which would yield their exact values, or they can be computed approximately (rather than exactly) by using finite difference formulas in conjunction with re-computations, e.g., ${\partial }^{2}R/\partial {\alpha }_i^2\cong \left[R\left({\alpha }_i^0+\delta {\alpha }_i\right)-2R\left({\alpha }_i^0\right)+R\left({\alpha }_i^0-\delta {\alpha }_i\right)\right]{\left|\delta {\alpha }_i\right|}^{-2}$.

7) Notably, contributions to the response sensitivities with respect to the uncertain boundary parameters can arise from either the model’s boundary conditions, from the definition of the model’s response or from both. It has been shown that in some cases, the contributions from the model’s boundary conditions may cancel partially the “direct-effect” contributions stemming from the response’s definition.

Ongoing work aims at extending the 2^nd-CASAM to include the consideration of additional responses (e.g., ratios of functionals), as well as consideration of coupled physical systems having common imprecisely known internal interfaces in phase-space.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Cacuci, D.G. (2020) Second-Order Adjoint Sensitivity Analysis Methodology for Computing Exactly Response Sensitivities to Uncertain Parameters and Boundaries of Linear Systems: Mathematical Framework. AJCM, 10, 329-354. https://doi.org/10.4236/ajcm.2020.103018
[2]	Cacuci, D.G. (2020) Illustrative Application of the 2nd-Order Adjoint Sensitivity Analysis Methodology to a Paradigm Linear Evolution/Transmission Model: Point-Detector Response, AJCM, 10, 355-381. https://doi.org/10.4236/ajcm.2020.103019
[3]	Cacuci, D.G., Badea, A.F., Badea, M.C. and Peltz, J.J. (2016) Efficient Computation of Operator-Type Response Sensitivities for Uncertainty Quantification and Predictive Modeling: Illustrative Application to a Spent Nuclear Fuel Dissolver Model. International Journal for Numerical Methods in Fluids, 83, 149-174. https://doi.org/10.1002/fld.4258
[4]	Cacuci, D.G. and Badea, M.C. (2018) Second-Order Adjoint Sensitivity Analysis of a Generic Evolution System with Applications to Radiation Transport and Fuel Reprocessing. PHYSOR-2018, Cancun, Mexico.
[5]	Cacuci, D.G. (2017) Inverse Predictive Modeling of Radiation Transport through Optically Thick Media in the Presence of Counting Uncertainties. Nuclear Science and Engineering, 186, 199-223. https://doi.org/10.1080/00295639.2017.1305244
[6]	Cacuci, D.G. and Favorite, J.A. (2018) Second-Order Sensitivity Analysis of Uncollided Particle Contributions to Radiation Detector Responses. Nuclear Science and Engineering, 190, 105-133. https://doi.org/10.1080/00295639.2018.1426899