First-Order Comprehensive Adjoint Sensitivity Analysis Methodology for Critical Points in Coupled Nonlinear Systems. I: Mathematical Framework

Abstract

This work presents the novel first-order comprehensive adjoint sensitivity analysis methodology for critical points (1st-CASAM-CP), which enables the exact and efficient computation of the first-order sensitivities of responses defined at critical points (maxima, minima, saddle points) of coupled nonlinear models of physical systems characterized by imprecisely known parameters underlying the models, boundaries, and interfaces between the coupled systems. Responses defined at critical points are important in many applications, including system optimization, safety analyses and licensing. For the design and licensing of nuclear reactors, such essentially important responses include the maximum temperatures of the fuel and cladding in hot channels. The 1st-CASAM-CP presented in this work makes it possible to determine, using a single large-scale “adjoint” computation, the first-order sensitivities of the magnitude of a response defined at a critical point of a function in the phase-space of the systems’ independent variables. In addition, the 1st-CASAM-CP enables the computation of the sensitivities of the location in phase-space of the critical point at which the respective response is located: one “adjoint” computation is required for each component of the respective critical point in the phase-space of independent variables. By enabling the exact and efficient computation of the sensitivities of responses and of their critical locations to imprecisely known model parameters, boundaries, and interfaces, the 1st-CASAM-CP significantly extends the practicality of analyzing crucially important responses for large-scale systems involving many uncertain parameters, interfaces, and boundaries.

1. Introduction

Computational models comprise independent and dependent variables, imprecisely known parameters, imprecisely known internal interfaces and external boundaries, and relations (equations) that relate all of these quantities to each other. The efficient computing of the exact functional derivatives (called “sensitivities”) of the result of interest (called “model response”) produced by the computational model has been of fundamental interest for predicting the changes in the response induced by changes in the parameters. The adjoint sensitivity analysis methodology for nonlinear systems introduced by Cacuci [1] has become the most efficient tool for sensitivity analysis [2,3,4,5] of large-scale systems involving a large number of parameters. Response sensitivities are also essential for computing the uncertainties in responses induced by uncertainties in the parameters, data assimilation [6,7,8], model calibration and predictive modeling [9].

Often, the result of interest produced by a computational model is located at a critical point (e.g., maximum, minimum, saddle point, bifurcation point, etc.) of the respective response in the phase-space of the model’s underlying independent variables. Such responses occur in optimization problems and also in safety-related considerations. For example, among the most important safety parameters for the licensing of nuclear reactors are the maximum fuel temperature, the maximum cladding temperature, the maximum admissible heat flux, etc. For such model responses, variations in the imprecisely known parameters and/or interfaces and/or boundaries will induce perturbations not only in the magnitude of the respective response but will also cause the critical point of the perturbed response to “move” in phase-space, i.e., to occur at another (perturbed) location in the phase-space of the model’s independent variables, as opposed to re-occurring at the original location of the critical point of the original, unperturbed model response. This fact has been recognized by Cacuci [10] in his pioneering work on developing an adjoint sensitivity analysis methodology for computing exactly and efficiently the sensitivities (functional derivatives) of the magnitude of such a response—as well as the sensitivities of the respective critical point(s)—with respect to the imprecisely known model parameters. The methodology conceived in [10] has been successfully applied [11] to a large-scale computational tool used for reactor physics and safety design, enabling the exact computation of the first-order sensitivities of the magnitude of the maximum clad temperature (and several other similar responses) to the many parameters underlying the computational model, as well as the sensitivities of the location of the respective maximum to the respective model parameters. Since the pioneering work described in [11], however, no other works have attempted to compute sensitivities of the locations of responses defined at critical points to uncertain model parameters. Furthermore, none of the works published thus far, including Cacuci’s original work [10], have the capabilities of computing sensitivities of critical points with respect to imprecisely known interfaces between coupled nonlinear subsystems and/or with respect to their imprecisely known external boundaries in phase space. This work alleviates this gap by extending the first-order comprehensive adjoint sensitivity analysis methodology (1st-CASAM) recently developed by Cacuci [12], to enable the efficient and exact computation of the first-order response sensitivities of model responses (of coupled nonlinear physical systems) defined at critical points with respect to the systems’ imprecisely known parameters, interfaces between systems and domain boundaries. Thus, the extended 1st-CASAM-CP presented in this work fundamentally generalizes all previously published works on this topic, and includes the following novel characteristics:

(i)

The system response considered in this work is a generic model response which is defined at a critical point in the phase-space of the coupled nonlinear models’ independent variables. The critical point is defined implicitly as the location where some (or all) of the first-order derivatives of the response with respect to some (or all) of the models’ independent variables vanish.

(ii)

The 1st-CASAM-CP enables the simultaneous, efficient, and exact computation of the 1st-order sensitivities of the magnitude of the model’s response to the imprecisely known parameters, interfaces, and domain boundaries of the coupled physical systems.

(iii)

The 1st-CASAM-CP also enables the simultaneous, efficient, and exact computation of the 1st-order sensitivities of the response’s critical point location (in the phase-space of independent variables) to the imprecisely known parameters, domain interfaces and boundaries of the coupled physical systems.

The application of the 1st-CASAM-CP developed in this work will be illustrated in an accompanying work [13] by considering a nuclear reactor thermal-hydraulics safety benchmark [14,15,16] which models heat conduction in a heated rod coupled to convection heat transport through a coolant that surrounded the rod. This benchmark [14,15,16] admits exact closed-form solutions for the sensitivities of the temperature distribution in the coupled rod/coolant system which can be used to benchmark thermal-hydraulics production codes. This benchmark is of fundamental importance both to the simulation of heat transport through reactor core channels, as well as to the design and operation of electrically heated rods in experimental thermal-hydraulics facilities. In particular, this benchmark was used to simulate the geometry of an advanced (“Generation-IV”) modular nuclear reactor [17] while verifying the numerical results produced by the commercially available FLUENT Adjoint Solver [18] and highlighting, in the process, various shortcomings of the current version of this code.

This work is structured as follows: Section 2 presents the mathematical definition of a response located at its critical point in the phase-space of the independent variables underlying two coupled generic nonlinear physical systems comprising imprecisely known parameters, interfaces, and boundaries. Section 3 presents the extended 1st-CASAM-CP for computing exactly and efficiently the first-order sensitivities of the magnitude of the response and of its phase-space location of its critical point to the physical systems’ imprecisely known parameters, interfaces, and boundaries. The discussion presented in Section 4 offers concluding remarks

2. Mathematical Modeling of a Response Defined at Critical Points of Two Generic Coupled Nonlinear Physical Systems

The mathematical model of the physical system considered in this work comprises two coupled nonlinear sub-systems which will be called “sub-system I” and, respectively, “sub-system II”. These sub-systems are considered to be coupled to each other across a common internal interface (boundary) in phase-space. The first subsystem is represented mathematically as follows:

$$N^{\left(I\right)}\left[u\left(x\right);\alpha \right]=Q^{\left(I\right)}\left(\alpha ;x\right)\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}x\in {\mathsf{\Omega }}_x\left(\alpha \right)$$

(1)

Bold letters will be used in this work to denote matrices and vectors. The superscript “I” will be used to denote quantities referring to “sub-system I.” Unless explicitly stated otherwise, the vectors in this work are considered to be column vectors. The second subsystem is represented mathematically as follows:

$$N^{\left(II\right)}\left[v\left(y\right);\alpha \right]=Q^{\left(II\right)}\left(\alpha ;y\right)\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}y\in {\mathsf{\Omega }}_y\left(\alpha \right)$$

(2)

The superscript “II” will be used to denote quantities referring to “sub-system II.” The sub-systems are coupled through interface conditions. Furthermore, when Equations (1) and (2) include differential operators, the interface conditions appear as internal boundary conditions, which may need to be supplemented by conditions on the sub-systems’ external boundaries, as well as initial/final conditions. Altogether, the interface, boundary and initial (or final time) conditions can be represented in operator form as follows:

$$B\left[u\left(x\right),v\left(y\right);\alpha ;x,y\right]=\hspace{0.17em}0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}x\in \partial {\mathsf{\Omega }}_x\left(\alpha \right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}y\in \partial {\mathsf{\Omega }}_y\left(\alpha \right)$$

(3)

The quantities (operators, vectors, scalars) appearing in Equations (1)–(3) are defined/described in Appendix A.

The system responses (i.e., quantities of interest) considered in previous works [12] were generic function-valued operators of the model’s dependent variables and parameters. In contradistinction to previous works [12], however, the specific response considered in this work is a scalar-valued functional defined at a critical point of the model’s dependent variable(s), including maxima, minima, saddle points or bifurcation points. A function defined on the phase-space domain of “Sub-system I” will be denoted as $F\left[u\left(x\right);\alpha ;x\right]$. The critical point of $F\left[u\left(x\right);\alpha ;x\right]$ will be denoted as $\xi \left(\alpha \right)\triangleq {\left[{\xi }_1\left(\alpha \right),\dots ,{\xi }_K\left(\alpha \right)\right]}^{†},\hspace{0.17em}\hspace{0.17em}K\le N_x$, and is defined implicitly as the point at which the first-order partial derivatives of $F$ with respect to some or all of the independent variables $x_k\hspace{0.17em},\hspace{0.17em}k=1,\dots ,K\le N_x$, vanish, i.e.,

$${\left\{\frac{\partial F\left[u\left(x\right);\alpha ;x\right]}{\partial x_k}\right\}}_{\zeta \left(\alpha \right)}=0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}k=1,\dots ,K\le N_x$$

(4)

The value (i.e., magnitude) of the response at the critical point $\xi \left(\alpha \right)\triangleq {\left[{\xi }_1\left(\alpha \right),\dots ,{\xi }_K\left(\alpha \right)\right]}^{†},\hspace{0.17em}\hspace{0.17em}K\le N_x$ of $F\left[u\left(x\right);\alpha ;x\right]$, which is defined implicitly by Equation (4), will be denoted as $R_x$ and is represented by the following functional:

$$R_x\triangleq {\displaystyle \prod _{i=1}^{N_x}{\displaystyle \underset{a_i\left(\alpha \right)}{\overset{b_i\left(\alpha \right)}{\int }}d{x}_i}}F\left[u\left(x\right);\alpha ;x,\right]{\displaystyle \prod _{k=1}^K\delta \left[x_k-{\xi }_k\left(\alpha \right)\right]}\hspace{0.17em}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}K\le N_x$$

(5)

where the following definition has been used:

$${\displaystyle \prod _{i=1}^{N_x}{\displaystyle \underset{a_i\left(\alpha \right)}{\overset{b_i\left(\alpha \right)}{\int }}d{x}_i}}\triangleq {\displaystyle \underset{a_1\left(\alpha \right)}{\overset{b_1\left(\alpha \right)}{\int }}d{x}_1\dots {\displaystyle \underset{a_{N_x}\left(\alpha \right)}{\overset{b_{N_x}\left(\alpha \right)}{\int }}d{x}_{N_x}}}\hspace{0.17em}$$

(6)

The phase-space coordinates ${\xi }_k\left(\alpha \right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}k=1,\dots ,K\le N_x,$ of the critical point $\xi \left(\alpha \right)\triangleq {\left[{\xi }_1\left(\alpha \right),\dots ,{\xi }_K\left(\alpha \right)\right]}^{†}$ of $F\left[u\left(x\right);\alpha ;x\right]$ are imprecisely known since they are functions of the uncertain values of the imprecisely known model parameters $\alpha \triangleq {\left({\alpha }_1,\dots ,{\alpha }_{Z_{\alpha }}\right)}^{†}\in {ℝ}^{Z_{\alpha }}$. Hence, the actual value/magnitude of the response $R_x$ is also unknown. The nominal value of the response at the critical point defined implicitly by Equation (4) will be denoted as $R_x^0$ and is determined by computing this response at the nominal values ${\alpha }^0$ and $u^0\left(x\right)$, i.e.,

$$R_x^0\triangleq {\displaystyle \prod _{i=1}^{N_x}{\displaystyle \underset{a_i\left({\alpha }^0\right)}{\overset{b_i\left({\alpha }^0\right)}{\int }}d{x}_i}}F\left[u^0\left(x\right);{\alpha }^0;x,\right]\hspace{0.17em}{\displaystyle \prod _{k=1}^K\delta \left[x_k-{\xi }_k\left({\alpha }^0\right)\right]}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}K\le N_x$$

(7)

Similarly, when the system response is a scalar-valued functional defined at a critical point of a function $G\left[v\left(y\right);\alpha ;y\right]$ pertaining to “Sub-system II,” the critical point of $G\left[v\left(y\right);\alpha ;y\right]$ will be denoted as $\eta \left(\alpha \right)\triangleq {\left[{\eta }_1\left(\alpha \right),\dots ,{\eta }_M\left(\alpha \right)\right]}^{†},\hspace{0.17em}\hspace{0.17em}M\le N_y$, and is defined implicitly as the point at which the first-order partial derivatives of $G$ with respect to some or all of the independent variables $y_j$ vanish, i.e.,:

$${\left\{\frac{\partial G\left[v\left(y\right);\alpha ;y\right]}{\partial y_m}\right\}}_{\eta \left(\alpha \right)}=0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}m=1,\dots ,M\le N_y$$

(8)

The imprecisely known value of the response at the critical point of $G\left[v\left(y\right);\alpha ;y\right]$, determined implicitly by Equation (8), will be denoted as $R_y$ and is defined as follows:

$$R_y\triangleq {\displaystyle \prod _{i=1}^{N_y}{\displaystyle \underset{c_i\left(\alpha \right)}{\overset{d_i\left(\alpha \right)}{\int }}d{y}_i}}\hspace{0.17em}\hspace{0.17em}G\left[v\left(y\right);\alpha ;y\right]{\displaystyle \prod _{m=1}^M\delta \left[y_m\hspace{0.17em}-{\eta }_m\left(\alpha \right)\right]}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}M\le N_y$$

(9)

The nominal value of the response at the critical point defined implicitly by Equation (9) will be denoted as $R_y^0$ and is determined by computing this response at the nominal values ${\alpha }^0$ and $v^0\left(y\right)$, i.e.,

$$R_y^0\triangleq {\displaystyle \prod _{i=1}^{N_y}{\displaystyle \underset{c_i\left({\alpha }^0\right)}{\overset{d_i\left({\alpha }^0\right)}{\int }}d{y}_i}\hspace{0.17em}G\left[v^0\left(y\right);{\alpha }^0;y\right]{\displaystyle \prod _{m=1}^M\delta \left[y_m\hspace{0.17em}-{\eta }_m\left({\alpha }^0\right)\right]}}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}M\le Z_y$$

(10)

3. Mathematical Framework of the 1st-CASAM-CP for Responses Defined at a Critical Point

The true but unknown values of the parameters will differ from their known nominal values by quantities denoted as $\delta \alpha \triangleq \left(\delta {\alpha }_1,\dots ,\delta {\alpha }_{Z_{\alpha }}\right)$, where $\delta {\alpha }_i\triangleq {\alpha }_i-{\alpha }_i^0$, $i=1,\dots ,Z_{\alpha }$. In turn, the parameter variations $\delta \alpha$ will cause variations $\delta u\left(x\right)\triangleq {\left[\delta u_1\left(x\right),\dots ,\delta u_{N_u}\left(x\right)\right]}^{†}$ and $\delta v\left(y\right)\triangleq {\left[\delta v_1\left(y\right),\dots ,\delta v_{N_v}\left(y\right)\right]}^{†}$ in the state functions, respectively. Furthermore, the variations $\delta \alpha$, $\delta u\left(x\right)$ and $\delta v\left(y\right)$ will cause variations in the response $R_x$ around the nominal response value $R_x^0$, as well as variations in the response $R_y$ around the nominal value $R_y^0$. Sensitivity analysis aims at computing the functional derivatives (called “sensitivities”) of the response under consideration to the imprecisely known parameters $\alpha$. Subsequently, these sensitivities can be used for a variety of purposes, including quantifying the uncertainties induced in responses by the uncertainties in the model and boundary parameters, combining the uncertainties in computed responses with uncertainties in measured response (“data assimilation”) to obtain more accurate predictions of responses and/or parameters (“model calibration,” “predictive modeling”), etc.

The 1st-order total sensitivity of a general vector-valued nonlinear operator $V\left[u\left(x\right),v\left(y\right);\alpha \right]$, defined on the complete domain ${\mathsf{\Omega }}_x\cup {\mathsf{\Omega }}_y$ on which the two coupled nonlinear systems are defined [details are provided in Appendix A] is defined [1,2,3,4,5] by the first-order Gateaux- (G-)variation of the respective operator. To simplify the notation in preparation for computing the G-variation of $V\left[u\left(x\right),v\left(y\right);\alpha \right]$, it is convenient to denote the functions appearing in the arguments of $V\left[u\left(x\right),v\left(y\right);\alpha \right]$ as being the components of a vector $e\triangleq {\left[u\left(x\right),v\left(y\right);\alpha \right]}^{†}$, which represents an arbitrary “point” in the combined phase-space of the state functions and parameters. The point which corresponds to the nominal values of the state functions and parameters in this phase space is denoted as $e^0\triangleq {\left[u^0\left(x\right),v^0\left(y\right);{\alpha }^0\right]}^{†}$. Analogously, it is convenient to consider the variations in the model’s state functions and parameters to be the components of a “vector of variations”, $\delta e$, defined as follows: $\delta e\triangleq {\left[\delta u\left(x\right),\delta v\left(y\right);\delta \alpha \right]}^{†}$. With these notations, the first-order G-variation of $V\left(e\right)$, denoted as $\delta V\left(e^0;\delta e\right)$, of $V\left(e\right)$ at $e^0$ in directions $\delta e$ is defined is as follows:

$$\delta V\left(e^0;\delta e\right)\triangleq {\left\{\frac{d}{d\epsilon }V\left[u^0\left(x\right)+\epsilon \delta u\left(x\right),v^0\left(y\right)+\epsilon \delta v\left(x\right);{\alpha }^0+\epsilon \delta \alpha \right]\right\}}_{\epsilon =0}$$

(11)

The G-differential $\delta V\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-differential $\delta V\left(e^0;\delta e\right)$ satisfies the relation $V\left(e^0+\epsilon \delta e\right)-V\left(e^0\right)=\delta V\left(e^0;\delta e\right)+\Delta \left(\delta e\right)\hspace{0.17em},$ with $\underset{\epsilon \to 0}{\lim }\left[\Delta \left(\epsilon \delta e\right)\right]/\epsilon =0\hspace{0.17em}.$ The existence of the G-variation $\delta V\left(e^0;\delta e\right)$ does not guarantee its numerical computability. Numerical methods most often require that $\delta V\left(e^0;\delta e\right)$ be linear in the variations $\delta e$ in a neighborhood $\left(e^0+\epsilon \delta e\right)$ around $e^0$. The necessary and sufficient conditions for the G-differential $\delta V\left(e^0;\delta e\right)$ of a nonlinear operator $V\left(e\right)$ to be linear in $\delta e$ in a neighborhood $\left(e^0+\epsilon \delta e\right)$ around $e^0$, and thus admit partial and total G-derivatives, are as follows:

(i)

For an arbitrary variation $\delta e$, the operator $V\left(e\right)$ must satisfy a weak Lipschitz condition at $e^0$, i.e.,

$$\Vert V\left(e^0+\epsilon \delta e\right)-V\left(e^0\right)\Vert \le k\Vert \epsilon e^0\Vert ,\hspace{0.17em}\hspace{0.17em}k<\infty$$

(12)

(ii)

For two arbitrary variations $\hspace{0.17em}\delta e_1$ and $\delta e_2$, the operator $V\left(e\right)$ must satisfy the following relation:

$$V\left(e^0+\epsilon \delta e_1+\epsilon \delta e_2\right)-V\left(e^0+\epsilon \delta e_1\right)-V\left(e^0+\epsilon \delta e_2\right)+V\left(e^0\right)=o\left(\epsilon \right)\hspace{0.17em}$$

(13)

In practice, however, the conditions shown in Equations (12) and (13) are seldom used directly. Rather, after the G-variation of the operator under consideration is obtained by applying the definition provided in Equation (11) and the resulting operator is subsequently verified (often by inspection) if it is linear (or not) in the variation $\delta e$. A G-variation $\delta V\left(e^0;\delta e\right)$ which is linear in $\delta e$ is called the G-differential of $V\left(e\right)$ at $e^0$ in directions $\delta e$. The G-differential of $V\left(e\right)$ is often denoted as $DV\left(e^0;\delta e\right)$ in order to emphasize its linearity in $\delta e$ and its identification with the total differential in customary calculus.

3.1. First-Order Sensitivities of the Response $R_x$ to the Imprecisely Known Model Parameters, Internal and External Boundaries

Section 3.1. presents the derivations of the first-order adjoint sensitivity system needed for the exact and efficient computation of the first-order sensitivities of the response $R_x$ with respect to the model’s uncertain parameters, including boundary and interface conditions. Section 3.2. presents the derivations of the first-order adjoint sensitivity system needed for the exact and efficient computation of the first-order sensitivities of the critical point $\xi \left(\alpha \right)\triangleq {\left[{\xi }_1\left(\alpha \right),\dots ,{\xi }_K\left(\alpha \right)\right]}^{†},\hspace{0.17em}\hspace{0.17em}K\le N_x$ at which the response $R_x$ is defined.

The total sensitivity of the response $R_x$ defined in Equation (6) with respect to variations in the parameters and state functions around their respective nominal values is provided by the 1st-order G-variation of the response $R_x$; this first-order G-variation will be denoted as $\delta R_x$. The total sensitivity $\delta R_x$ is obtained, by definition, as follows:

$$\begin{array}{l}\delta R_x\triangleq \{\frac{d}{d\epsilon }{\displaystyle \underset{a_1\left({\alpha }^0+\epsilon \delta \alpha \right)}{\overset{b_1\left({\alpha }^0+\epsilon \delta \alpha \right)}{\int }}d{x}_1\dots {\displaystyle \underset{a_{N_x}\left({\alpha }^0+\epsilon \delta \alpha \right)}{\overset{b_{N_x}\left({\alpha }^0+\epsilon \delta \alpha \right)}{\int }}d{x}_{N_x}}}F\left[u^0\left(x\right)+\epsilon \delta u\left(x\right);{\alpha }^0+\epsilon \delta \alpha ;x\right]\\ {\times {\displaystyle \prod _{k=1}^K\delta \left[x_k-{\xi }_k\left({\alpha }^0+\epsilon \delta \alpha \right)\right]}\}}_{\epsilon =0}={\displaystyle \sum _{p=1}^{Z_{\alpha }}{\left\{\frac{\partial R_x}{\partial {\alpha }_p}\right\}}_{\left(u^0,{\alpha }^0\right)}\delta {\alpha }_p}={\left\{\delta R_x\right\}}^{dir}+{\left\{\delta R_x\right\}}^{ind},\end{array}$$

(14)

where the so-called “indirect-effect term” ${\left\{\delta R_x\right\}}^{ind}$ depends on the unknown variation $\delta u\left(x\right)$ and is defined as follows:

$${\left\{\delta R_x\right\}}^{ind}\triangleq {\displaystyle \prod _{i=1}^{N_x}{\displaystyle \underset{a_i\left({\alpha }^0\right)}{\overset{b_i\left({\alpha }^0\right)}{\int }}d{x}_i}}\hspace{0.17em}{\left\{\left[\frac{\partial F\left(u;\alpha ;x\right)}{\partial u}\delta u\right]{\displaystyle \prod _{k=1}^K\delta \left[x_k-{\xi }_k\left(\alpha \right)\right]}\right\}}_{\left(u^0;{\alpha }^0\right)}\hspace{0.17em}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}K\le N_x$$

(15)

while the so-called “direct-effect term” ${\left\{\delta R_x\right\}}^{dir}$ does not depend on the unknown variation $\delta u\left(x\right)$ but only depends on the parameter variations $\delta \alpha$ and is defined as follows:

$$\begin{array}{l}{\left\{\delta R_x\right\}}^{dir}={\displaystyle \sum _{p=1}^{Z_{\alpha }}{\left\{\frac{\partial R_x^{dir}}{\partial {\alpha }_p}\right\}}_{\left(e^0\right)}\delta {\alpha }_p}=\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\displaystyle \prod _{i=1}^{N_x}{\displaystyle \underset{a_i\left({\alpha }^0\right)}{\overset{b_i\left({\alpha }^0\right)}{\int }}d{x}_i}}\hspace{0.17em}{\left\{\left[\frac{\partial F\left(u;\alpha ;x\right)}{\partial \alpha }\delta \alpha \right]{\displaystyle \prod _{k=1}^K\delta \left[x_k-{\xi }_k\left(\alpha \right)\right]}\hspace{0.17em}\right\}}_{\left(u^0;{\alpha }^0\right)}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+{\displaystyle \sum _{i=1}^{N_x}{\displaystyle \prod _{j=1,j\ne i}^{N_x}{\displaystyle \underset{a_j\left({\alpha }^0\right)}{\overset{b_j\left({\alpha }^0\right)}{\int }}d{x}_j}}}{\left\{{\left[F\left(u;\alpha ;x\right)\right]}_{x_i=b_i\left({\alpha }^0\right)}\left[\frac{\partial b_i\left(\alpha \right)}{\partial \alpha }\delta \alpha \right]{\displaystyle \prod _{k=1}^K\delta \left[x_k-{\xi }_k\left(\alpha \right)\right]}\right\}}_{\left(u^0,{\alpha }^0\right)}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-{\displaystyle \sum _{i=1}^{N_x}{\displaystyle \prod _{j=1,j\ne i}^{N_x}{\displaystyle \underset{a_j\left({\alpha }^0\right)}{\overset{b_j\left({\alpha }^0\right)}{\int }}d{x}_j}}}{\left\{{\left[F\left(u;\alpha ;x\right)\right]}_{x_i=a_i\left({\alpha }^0\right)}\left[\frac{\partial a_i\left(\alpha \right)}{\partial \alpha }\delta \alpha \right]{\displaystyle \prod _{k=1}^K\delta \left[x_k-{\xi }_k\left(\alpha \right)\right]}\right\}}_{\left(u^0,{\alpha }^0\right)}\\ -\hspace{0.17em}{\displaystyle \sum _{k=1}^K{\left\{{\displaystyle \prod _{i=1}^{N_x}{\displaystyle \underset{a_i\left({\alpha }^0\right)}{\overset{b_i\left({\alpha }^0\right)}{\int }}d{x}_i}}\hspace{0.17em}F\left(u;\alpha ;x\right){\delta }^{\prime }\left[x_k-{\xi }_k\left(\alpha \right)\right]\left(\frac{\partial {\xi }_k\left(\alpha \right)}{\partial \alpha }\delta \alpha \right){\displaystyle \prod _{m=1,m\ne k}^K\delta \left[x_m-{\xi }_m\left(\alpha \right)\right]}\right\}}_{\left(u^0;{\alpha }^0\right)}\hspace{0.17em}},\end{array}$$

(16)

The following notations have been used for the various (column) vectors appearing in Equations (15) and (16):

$$\frac{\partial F\left(u;\alpha ;x\right)}{\partial u}\triangleq {\left[\frac{\partial F}{\partial u_1},\dots ,\frac{\partial F}{\partial u_{N_u}}\right]}^{†},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\frac{\partial F\left(u;\alpha ;x\right)}{\partial u}\delta u={\displaystyle \sum _{i=1}^{N_u}\frac{\partial F\left(u;\alpha ;x\right)}{\partial u_i}\delta u_i}\hspace{0.17em},$$

(17)

$$\frac{\partial F\left(u;\alpha ;x\right)}{\partial \alpha }\triangleq {\left[\frac{\partial F}{\partial {\alpha }_1},\dots ,\frac{\partial F}{\partial {\alpha }_{Z_{\alpha }}}\right]}^{†};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\frac{\partial F\left(u;\alpha ;x\right)}{\partial \alpha }\delta \alpha ={\displaystyle \sum _{p=1}^{Z_{\alpha }}\frac{\partial F\left(u;\alpha ;x\right)}{\partial {\alpha }_p}\delta {\alpha }_p}\hspace{0.17em}.$$

(18)

The vectors $\partial {\xi }_j\left(\alpha \right)/\partial \alpha$, $\partial b_k\left(\alpha \right)/\partial \alpha$, and $\partial a_k\left(\alpha \right)/\partial \alpha$, which appear in Equation (16), each have the same structure as the vector $\partial F\left(\alpha \right)/\partial \alpha$ defined in Equation (18). Thus, the partial derivative $\partial R_x^{dir}/\partial {\alpha }_p$, which is implicitly defined in Equation (16), has the following expression:

$$\begin{array}{l}{\left\{\frac{\partial R_x^{dir}}{\partial {\alpha }_p}\right\}}_{\left(e^0\right)}=\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\displaystyle \prod _{i=1}^{N_x}{\displaystyle \underset{a_i\left({\alpha }^0\right)}{\overset{b_i\left({\alpha }^0\right)}{\int }}d{x}_i}}\hspace{0.17em}{\left\{\left[\frac{\partial F\left(u;\alpha ;x\right)}{\partial {\alpha }_p}\right]{\displaystyle \prod _{k=1}^K\delta \left[x_k-{\xi }_k\left(\alpha \right)\right]}\hspace{0.17em}\right\}}_{\left(u^0;{\alpha }^0\right)}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+{\displaystyle \sum _{i=1}^{N_x}{\displaystyle \prod _{j=1,j\ne i}^{N_x}{\displaystyle \underset{a_j\left({\alpha }^0\right)}{\overset{b_j\left({\alpha }^0\right)}{\int }}d{x}_j}}}{\left\{{\left[F\left(u;\alpha ;x\right)\right]}_{x_i=b_i\left({\alpha }^0\right)}\frac{\partial b_i\left(\alpha \right)}{\partial {\alpha }_p}{\displaystyle \prod _{k=1}^K\delta \left[x_k-{\xi }_k\left(\alpha \right)\right]}\right\}}_{\left(u^0,{\alpha }^0\right)}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-{\displaystyle \sum _{i=1}^{N_x}{\displaystyle \prod _{j=1,j\ne i}^{N_x}{\displaystyle \underset{a_j\left({\alpha }^0\right)}{\overset{b_j\left({\alpha }^0\right)}{\int }}d{x}_j}}}{\left\{{\left[F\left(u;\alpha ;x\right)\right]}_{x_i=a_i\left({\alpha }^0\right)}\frac{\partial a_i\left(\alpha \right)}{\partial {\alpha }_p}{\displaystyle \prod _{k=1}^K\delta \left[x_k-{\xi }_k\left(\alpha \right)\right]}\right\}}_{\left(u^0,{\alpha }^0\right)}\\ -\hspace{0.17em}{\displaystyle \sum _{k=1}^K{\left\{{\displaystyle \prod _{i=1}^{N_x}{\displaystyle \underset{a_i\left({\alpha }^0\right)}{\overset{b_i\left({\alpha }^0\right)}{\int }}d{x}_i}}\hspace{0.17em}F\left(u;\alpha ;x\right){\delta }^{\prime }\left[x_k-{\xi }_k\left(\alpha \right)\right]\frac{\partial {\xi }_k\left(\alpha \right)}{\partial {\alpha }_p}{\displaystyle \prod _{m=1,m\ne k}^K\delta \left[x_m-{\xi }_m\left(\alpha \right)\right]}\right\}}_{\left(u^0;{\alpha }^0\right)}\hspace{0.17em}}.\end{array}$$

(19)

As indicated by Equation (16), the “direct-effect” term ${\left\{\delta R_x\right\}}^{dir}$ depends only on the parameter variations $\delta \alpha$ and can therefore be computed immediately. On the other hand, the “indirect-effect” term ${\left\{\delta R_x\right\}}^{ind}$ depends indirectly on the parameter variations $\delta \alpha$ through the yet unknown variation $\delta u$. The variations $\delta u$ and $\delta v\left(y\right)$ are related to the parameter variations $\delta \alpha$ through the “first-level forward sensitivity system” (1st-LFSS), which is obtained by G-differentiating the original equations that underly the coupled nonlinear sub-systems I and II. Applying the definition of the first-order G-variation to Equations (1)–(3) yields the following relations:

$${\left\{\frac{d}{d\epsilon }{N}^{\left(I\right)}\left[u^0\left(x\right)+\epsilon \delta u\left(x\right);{\alpha }^0+\epsilon \delta \alpha \right]\right\}}_{\epsilon =0}=\hspace{0.17em}{\left\{\frac{d}{d\epsilon }{Q}^{\left(I\right)}\left({\alpha }^0+\epsilon \delta \alpha ;x\right)\right\}}_{\epsilon =0},x\in {\mathsf{\Omega }}_x\left({\alpha }^0\right),$$

(20)

$${\left\{\frac{d}{d\epsilon }{N}^{\left(II\right)}\left[v^0\left(y\right)+\epsilon \delta v\left(y\right);{\alpha }^0+\epsilon \delta \alpha \right]\right\}}_{\epsilon =0}=\hspace{0.17em}{\left\{\frac{d}{d\epsilon }{Q}^{\left(II\right)}\left({\alpha }^0+\epsilon \delta \alpha ;y\right)\right\}}_{\epsilon =0},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}y\in {\mathsf{\Omega }}_y\left({\alpha }^0\right),$$

(21)

$${\left\{\frac{d}{d\epsilon }B\left[u^0\left(x\right)+\epsilon \delta u\left(x\right),v^0\left(y\right)+\epsilon \delta v\left(y\right);{\alpha }^0+\epsilon \delta \alpha ;x,y\right]\right\}}_{\epsilon =0}=\hspace{0.17em}0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}x\in \partial {\mathsf{\Omega }}_x\left({\alpha }^0\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}y\in \partial {\mathsf{\Omega }}_y\left({\alpha }^0\right).$$

(22)

Performing in Equations (20)–(22) the differentiations with respect to $\epsilon$ and setting $\epsilon =0$ in the resulting expression yields the following system of equations:

$$\delta N^{\left(I\right)}\left[u^0\left(x\right),{\alpha }^0;\delta u\left(x\right),\delta \alpha \right]=\delta Q^{\left(I\right)}\left({\alpha }^0;\delta \alpha \right)\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}x\in {\mathsf{\Omega }}_x\left({\alpha }^0\right)$$

(23)

$$\delta N^{\left(II\right)}\left[v^0\left(y\right),{\alpha }^0;\delta v\left(y\right),\delta \alpha \right]=\delta Q^{\left(II\right)}\left({\alpha }^0;\delta \alpha \right)\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}y\in {\mathsf{\Omega }}_y\left({\alpha }^0\right),$$

(24)

$$\delta B\left(e^0;\delta e\right)=\hspace{0.17em}0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}x\in \partial {\mathsf{\Omega }}_x\left({\alpha }^0\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}y\in \partial {\mathsf{\Omega }}_y\left({\alpha }^0\right)$$

(25)

The existence of the G-variations in Equations (23)–(25) does not guarantee their numerical computability. Numerical methods most often require that these G-variations be linear in the variations $\delta e$ in a neighborhood $\left(e^0+\epsilon \delta e\right)$ around $e^0$, implying that the operators $N^{\left(I\right)},\hspace{0.17em}\hspace{0.17em}{Q}^{\left(I\right)},\hspace{0.17em}\hspace{0.17em}{N}^{\left(II\right)},\hspace{0.17em}{Q}^{\left(II\right)},\hspace{0.17em}\hspace{0.17em}B$ must satisfy the necessary and sufficient conditions provided in Equations (12) and (13), which will henceforth be assumed to the case, so that Equations (23)–(25) can be written in the following form:

$$\begin{array}{l}{\left\{\left[\begin{array}{cc}\frac{\partial N^{\left(I\right)}\left(u;\alpha \right)}{\partial u}& 0\\ 0& \frac{\partial N^{\left(II\right)}\left(v;\alpha \right)}{\partial v}\end{array}\right]\hspace{0.17em}\right\}}_{\left(e^0\right)}\left(\begin{array}{c}\delta u\left(x\right)\\ \delta v\left(y\right)\end{array}\right)={\left\{\left(\begin{array}{c}{Q}_1^{\left(1\right)}\left(u;\alpha ;\delta \alpha \right)\\ Q_2^{\left(1\right)}\left(v;\alpha ;\delta \alpha \right)\end{array}\hspace{0.17em}\right)\hspace{0.17em}\right\}}_{\left(e^0\right)},\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}for\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}x\in {\mathsf{\Omega }}_x\left({\alpha }^0\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}y\in {\mathsf{\Omega }}_y\left({\alpha }^0\right),\end{array}$$

(26)

$$\begin{array}{l}{\left\{\frac{\partial B\left[u\left(x\right),v\left(y\right);\alpha ;x,y\right]}{\partial u}\right\}}_{\left(e^0\right)}\delta u\left(x\right)+{\left\{\frac{\partial B\left[u\left(x\right),v\left(y\right);\alpha ;x,y\right]}{\partial v}\right\}}_{\left(e^0\right)}\delta v\left(y\right)\hspace{0.17em}\\ +{\left\{\frac{\partial B\left[u\left(x\right),v\left(y\right);\alpha ;x,y\right]}{\partial \alpha }\right\}}_{\left(e^0\right)}\delta \alpha =\hspace{0.17em}0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}x\in \partial {\mathsf{\Omega }}_x\left({\alpha }^0\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}y\in \partial {\mathsf{\Omega }}_y\left({\alpha }^0\right).\end{array}$$

(27)

where

$$\begin{array}{l}{\left\{Q_1^{\left(1\right)}\left(u;\alpha ;\delta \alpha \right)\hspace{0.17em}\right\}}_{\left(u^0;{\alpha }^0\right)}\triangleq {\left\{\frac{\partial \left[Q^{\left(I\right)}\left(\alpha ;x\right)-N^{\left(I\right)}\left[u\left(x\right);\alpha \right]\right]}{\partial \alpha }\hspace{0.17em}\right\}}_{\left(u^0;{\alpha }^0\right)}\delta \alpha \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}={\displaystyle \sum _{i=1}^{Z_{\alpha }}{\left\{\hspace{0.17em}\frac{\partial \left[Q^{\left(I\right)}\left(\alpha ;x\right)-N^{\left(I\right)}\left[u\left(x\right);\alpha \right]\right]}{\partial {\alpha }_i}\right\}}_{\left(u^0;{\alpha }^0\right)}\delta {\alpha }_i},\end{array}$$

(28)

$$\begin{array}{l}\hspace{0.17em}{\left\{Q_2^{\left(1\right)}\left(v;\alpha ;\delta \alpha \right)\hspace{0.17em}\right\}}_{\left(v^0,{\alpha }^0\right)}\triangleq {\left\{\frac{\partial \left[Q^{\left(II\right)}\left(\alpha ;y\right)-N^{\left(II\right)}\left[v\left(y\right);\alpha \right]\right]}{\partial \alpha }\hspace{0.17em}\right\}}_{\left(v^0,{\alpha }^0\right)}\delta \alpha \\ \hspace{0.17em}\hspace{0.17em}={\displaystyle \sum _{i=1}^{Z_{\alpha }}{\left\{\frac{\partial \left[Q^{\left(II\right)}\left(\alpha ;y\right)-N^{\left(II\right)}\left[v\left(y\right);\alpha \right]\right]}{\partial {\alpha }_i}\hspace{0.17em}\right\}}_{\left(v^0,{\alpha }^0\right)}\delta {\alpha }_i},\end{array}$$

(29)

The partial G-derivatives $\partial N^{\left(I\right)}\left[u\left(x\right);\alpha \right]/\partial u$, $\partial N^{\left(II\right)}\left[v\left(y\right);\alpha \right]/\partial v$, $\partial B\left[u\left(x\right),v\left(y\right);\alpha ;x,y\right]/\partial u$, $\partial B\left[u\left(x\right),v\left(y\right);\alpha ;x,y\right]/\partial v$, $\partial N^{\left(I\right)}\left[u\left(x\right);\alpha \right]/\partial \alpha$, $\partial N^{\left(II\right)}\left[v\left(y\right);\alpha \right]/\partial \alpha$, $\partial Q^{\left(I\right)}\left[u\left(x\right);\alpha \right]/\partial \alpha$, $\partial Q^{\left(II\right)}\left[u\left(x\right);\alpha \right]/\partial \alpha$ and $\partial B\left[u\left(x\right),v\left(y\right);\alpha ;x,y\right]/\partial \alpha$, which appear in Equations (26)–(29), are matrices of corresponding dimensions.

The 1st-LFSS, comprising Equations (26) and (27), could be solved to obtain the variations $\delta u$ and $\delta v$ in the state functions in terms of the parameter variations $\delta \alpha$ which appear as sources in the 1st-LFSS equations. Subsequently, the variations $\delta u$ and $\delta v$ thus obtained could be used to compute the total sensitivity $\partial R_x$ using Equation (14). However, computing repeatedly the 1st-LFSS for every possible parameter variation $\delta {\alpha }_p,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}p=1,\dots ,Z_{\alpha }$ becomes prohibitively expensive for large scale systems involving many uncertain parameters. These repeated computations can be circumvented by extending the concepts underlying the “Adjoint Sensitivity Analysis Methodology” conceived by Cacuci [2,3] to construct a “First-Level Adjoint Sensitivity System” (1st-LASS), the solution of which will be independent of the variations $\delta \alpha$, $\delta u$ and $\delta v$. Subsequently, the solution of the 1st-LASS will be used to compute the indirect-effect term ${\left\{\delta R_x\right\}}^{ind}$ by constructing an equivalent expression (for this indirect-effect term) which would not involve the unknown variations $\delta u$ and $\delta v$.

The construction of the requisite 1st-LASS is achieved by implementing the following sequence of steps:

• Consider the union of independent variables describing the combined domain ${\mathsf{\Omega }}_x\cup {\mathsf{\Omega }}_y$ to be the components of the following $N_z$-dimensional column vector $z\triangleq {\left[z_1,\dots ,z_{N_z}\right]}^{†}$.
• Introduce a Hilbert space pertaining to the domain ${\mathsf{\Omega }}_x\cup {\mathsf{\Omega }}_y$, denoted as $\mathrm{H}$, comprising square-integrable vector-valued elements of the form $f^{\left(\alpha \right)}\left(z\right)\triangleq {\left[g^{\left(\alpha \right)}\left(x\right),h^{\left(\alpha \right)}\left(y\right)\right]}^{†}\in \mathrm{H}$ and $f^{\left(\beta \right)}\left(z\right)\triangleq {\left[g^{\left(\beta \right)}\left(x\right),h^{\left(\beta \right)}\left(y\right)\right]}^{†}\in \mathrm{H}$, where $g^{\left(\alpha \right)}\left(x\right)\triangleq {\left[g_1^{\left(\alpha \right)}\left(x\right),\dots ,g_{N_u}^{\left(\alpha \right)}\left(x\right)\right]}^{†}$, $g^{\left(\beta \right)}\left(x\right)\triangleq {\left[g_1^{\left(\beta \right)}\left(x\right),\dots ,g_{N_u}^{\left(\beta \right)}\left(y\right)\right]}^{†}$, $h^{\left(\alpha \right)}\left(y\right)\triangleq {\left[h_1^{\left(\alpha \right)}\left(y\right),\dots ,h_{N_v}^{\left(\alpha \right)}\left(y\right)\right]}^{†}$, $h^{\left(\beta \right)}\left(y\right)\triangleq {\left[h_1^{\left(\beta \right)}\left(y\right),\dots ,h_{N_v}^{\left(\beta \right)}\left(y\right)\right]}^{†}$.
• Define the inner product, denoted as $\langle f^{\left(\alpha \right)}\left(z\right),f^{\left(\beta \right)}\left(z\right)\rangle$, between two elements of $\mathrm{H}$, as follows:

  $$\langle f^{\left(\alpha \right)}\left(z\right),f^{\left(\beta \right)}\left(z\right)\rangle \triangleq {\displaystyle \prod _{i=1}^{N_x}{\displaystyle \underset{a_i\left({\alpha }^0\right)}{\overset{b_i\left({\alpha }^0\right)}{\int }}d{x}_i}}\left[g^{\left(\alpha \right)}\left(x\right)·g^{\left(\beta \right)}\left(x\right)\right]+{\displaystyle \prod _{i=1}^{N_y}{\displaystyle \underset{c_i\left({\alpha }^0\right)}{\overset{d_i\left({\alpha }^0\right)}{\int }}d{y}_i}}\left[h^{\left(\alpha \right)}\left(y\right)·h^{\left(\beta \right)}\left(y\right)\right]$$

  (30)

  where

  $$g^{\left(\alpha \right)}\left(x\right)·g^{\left(\beta \right)}\left(x\right)\triangleq {\displaystyle \sum _{n=1}^{N_u}{g}_n^{\left(\alpha \right)}\left(x\right)g_n^{\left(\beta \right)}\left(x\right)}$$

  (31)

  $$h^{\left(\alpha \right)}\left(y\right)·h^{\left(\beta \right)}\left(y\right)\triangleq {\displaystyle \sum _{n=1}^{N_v}{h}_n^{\left(\alpha \right)}\left(y\right)h_n^{\left(\beta \right)}\left(y\right)}$$

  (32)
• Use the definition provided in Equation (30) to form the inner product of Equation (26) with a square-integrable vector ${\psi }^{\left(1\right)}\left(z\right)\triangleq {\left[{\psi }^{\left(I\right)}\left(x\right),{\psi }^{\left(II\right)}\left(y\right)\right]}^{†}\in \mathrm{H}$ to obtain the following relation:

  $$\begin{array}{l}\langle {\left(\begin{array}{c}{\psi }^{\left(I\right)}\left(x\right)\\ {\psi }^{\left(II\right)}\left(y\right)\end{array}\right)}^{†},\hspace{0.17em}\hspace{0.17em}{\left\{\left[\begin{array}{cc}\frac{\partial N^{\left(I\right)}\left(u;\alpha \right)}{\partial u}& 0\\ 0& \frac{\partial N^{\left(II\right)}\left(v;\alpha \right)}{\partial v}\end{array}\right]\hspace{0.17em}\right\}}_{\left(e^0\right)}\left(\begin{array}{c}\delta u\left(x\right)\\ \delta v\left(y\right)\end{array}\right)\hspace{0.17em}\rangle \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}=\langle {\left(\begin{array}{c}{\psi }^{\left(I\right)}\left(x\right)\\ {\psi }^{\left(II\right)}\left(y\right)\end{array}\right)}^{†},\hspace{0.17em}\hspace{0.17em}{\left\{\left(\begin{array}{c}{Q}_1^{\left(1\right)}\left(u;\alpha ;\delta \alpha \right)\\ Q_2^{\left(1\right)}\left(v;\alpha ;\delta \alpha \right)\end{array}\hspace{0.17em}\right)\hspace{0.17em}\right\}}_{\left(e^0\right)}\hspace{0.17em}\rangle .\end{array}$$

  (33)
• Using the definition of the adjoint operator in the Hilbert space $\mathrm{H}$, recast the left-side of Equation (33) as follows:

  $$\begin{array}{l}\langle {\left(\begin{array}{c}{\psi }^{\left(I\right)}\left(x\right)\\ {\psi }^{\left(II\right)}\left(y\right)\end{array}\right)}^{†},\hspace{0.17em}\hspace{0.17em}{\left\{\left[\begin{array}{cc}\frac{\partial N^{\left(I\right)}\left(u;\alpha \right)}{\partial u}& 0\\ 0& \frac{\partial N^{\left(II\right)}\left(v;\alpha \right)}{\partial v}\end{array}\right]\hspace{0.17em}\right\}}_{\left(e^0\right)}\left(\begin{array}{c}\delta u\left(x\right)\\ \delta v\left(y\right)\end{array}\right)\hspace{0.17em}\rangle \\ =\langle {\left(\begin{array}{c}\delta u\left(x\right)\\ \delta v\left(y\right)\end{array}\right)}^{†},\hspace{0.17em}{\left\{\left[\begin{array}{cc}{A}^{*}\left(u;\alpha \right)& 0\\ 0& B^{*}\left(v;\alpha \right)\end{array}\right]\hspace{0.17em}\right\}}_{\left(e^0\right)}\hspace{0.17em}\left(\begin{array}{c}{\psi }^{\left(I\right)}\left(x\right)\\ {\psi }^{\left(II\right)}\left(y\right)\end{array}\right)\hspace{0.17em}\rangle \\ +{\left\{B{C}^{\left(1\right)}{\left[u\left(x\right);v\left(y\right);{\psi }^{\left(I\right)}\left(x\right),{\psi }^{\left(II\right)}\left(y\right);\delta u\left(x\right),\delta v\left(y\right);\alpha ;x,y;\delta \alpha \right]}_{{{}_{\partial {\mathsf{\Omega }}_x\cup }}_{\partial {\mathsf{\Omega }}_y}}\right\}}_{\left(e^0\right)},\end{array}$$

  (34)

  where the operator $A^{*}\left(u;\alpha \right)$ denotes the formal adjoint of $\partial A\left(u;\alpha \right)/\partial u$, the operator $B^{*}\left(v;\alpha \right)$ denotes the formal adjoint of $\partial B\left(v;\alpha \right)/\partial v$, and where ${\left\{B{C}^{\left(1\right)}{\left[u\left(x\right);v\left(y\right);{\psi }^{\left(I\right)}\left(x\right),{\psi }^{\left(II\right)}\left(y\right);\delta u\left(x\right),\delta v\left(y\right);\alpha ;x,y;\delta \alpha \right]}_{{{}_{\partial {\mathsf{\Omega }}_x\cup }}_{\partial {\mathsf{\Omega }}_y}}\right\}}_{\left(e^0\right)}$ denotes the bilinear concomitant evaluated on the combined (internal and external) boundary $\delta {\mathsf{\Omega }}_x\cup \delta {\mathsf{\Omega }}_y$.
• Replace the left-side of Equation (33) with the right-side of Equation (34) and re-arrange the resulting equation to obtain the following relation:

  $$\begin{array}{l}\langle {\left(\begin{array}{c}\delta u\left(x\right)\\ \delta v\left(y\right)\end{array}\right)}^{†},\hspace{0.17em}{\left\{\left[\begin{array}{cc}{A}^{*}\left(u;\alpha \right)& 0\\ 0& B^{*}\left(v;\alpha \right)\end{array}\right]\hspace{0.17em}\right\}}_{\left(e^0\right)}\hspace{0.17em}\left(\begin{array}{c}{\psi }^{\left(I\right)}\left(x\right)\\ {\psi }^{\left(II\right)}\left(y\right)\end{array}\right)\hspace{0.17em}\rangle \\ =\langle {\left(\begin{array}{c}{\psi }^{\left(I\right)}\left(x\right)\\ {\psi }^{\left(II\right)}\left(y\right)\end{array}\right)}^{†},\hspace{0.17em}\hspace{0.17em}{\left\{\left(\begin{array}{c}{Q}_1^{\left(1\right)}\left(u;\alpha ;\delta \alpha \right)\\ Q_2^{\left(1\right)}\left(v;\alpha ;\delta \alpha \right)\end{array}\hspace{0.17em}\right)\hspace{0.17em}\right\}}_{\left(e^0\right)}\hspace{0.17em}\rangle \\ -{\left\{B{C}^{\left(1\right)}{\left[u\left(x\right);v\left(y\right);{\psi }^{\left(I\right)}\left(x\right),{\psi }^{\left(II\right)}\left(y\right);\delta u\left(x\right),\delta v\left(y\right);\alpha ;x,y;\delta \alpha \right]}_{{{}_{\partial {\mathsf{\Omega }}_x\cup }}_{\partial {\mathsf{\Omega }}_y}}\right\}}_{\left(e^0\right)}.\end{array}$$

  (35)

Require the left-side of Equation (35) to represent the indirect-effect term ${\left\{\delta R_x\right\}}^{ind}$ defined in Equation (15), which can be fulfilled by requiring the yet undetermined (adjoint) functions ${\psi }^{\left(I\right)}\left(x\right)$ and ${\psi }^{\left(II\right)}\left(y\right)$ to satisfy the following (adjoint) equations:

$${\left\{A^{*}\left(u;\alpha \right)\hspace{0.17em}\right\}}_{\left(e^0\right)}{\psi }^{\left(I\right)}\left(x\right)={\left\{\frac{\partial F\left(u;\alpha ;x\right)}{\partial u}{\displaystyle \prod _{k=1}^K\delta \left[x_k-{\xi }_k\left(\alpha \right)\right]}\right\}}_{\left(u^0;{\alpha }^0\right)},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}K\le N_x,$$

(36)

$${\left\{B^{*}\left(v;\alpha \right)\hspace{0.17em}\right\}}_{\left(e^0\right)}{\psi }^{\left(II\right)}\left(y\right)=0$$

(37)

7.

The boundary, interface, and initial/final conditions for the functions ${\psi }^{\left(I\right)}\left(x\right)$ and ${\psi }^{\left(II\right)}\left(y\right)$ are now determined as follows:

(i)

The boundary, interface and initial/final conditions given in Equation (27) are substituted into the bilinear concomitant in Equation (35).

(ii)

The remaining unknown boundary, interface and initial/final conditions involving the functions $\delta u\left(x\right)$ and $\delta v\left(y\right)$ are eliminated from ${\left\{B{C}^{\left(1\right)}{\left[u\left(x\right);v\left(y\right);{\psi }^{\left(I\right)}\left(x\right),{\psi }^{\left(II\right)}\left(y\right);\delta u\left(x\right),\delta v\left(y\right);\alpha ;x,y;\delta \alpha \right]}_{{{}_{\partial {\mathsf{\Omega }}_x\cup }}_{\partial {\mathsf{\Omega }}_y}}\right\}}_{\left(e^0\right)}$ by selecting boundary, interface and initial/final conditions for the adjoint functions ${\psi }^{\left(I\right)}\left(x\right)$ and ${\psi }^{\left(II\right)}\left(y\right)$ such that: (a) the selected conditions for these adjoint functions must be independent of unknown values of $\delta u\left(x\right)$, $\delta v\left(y\right)$ and $\delta \alpha$; and (b) Equations (36) and (37) are well posed in conjunction with the boundary conditions chosen for the adjoint functions ${\psi }^{\left(I\right)}\left(x\right)$ and ${\psi }^{\left(II\right)}\left(y\right)$.

(iii)

In operator form, the boundary, interface, and initial/final conditions thus obtained for the adjoint functions ${\psi }^{\left(I\right)}\left(x\right)$ and ${\psi }^{\left(II\right)}\left(y\right)$ can be represented as follows:

$${\left\{C_A^{\left(1\right)}\left[u\left(x\right);v\left(y\right);{\psi }^{\left(I\right)}\left(x\right),{\psi }^{\left(II\right)}\left(y\right);\alpha ;x,y\right]\right\}}_{\left(e^0\right)}=0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}x\in \partial {\mathsf{\Omega }}_x,\hspace{0.17em}\hspace{0.17em}y\in \partial {\mathsf{\Omega }}_y$$

(38)

where the subscript “A“ indicates “adjoint” and the superscript “(1)” indicates “first-level.” The system of equations represented by Equations (36) and (37) together with the boundary/interface and initial/final time conditions represented by Equation (38) constitute the 1st-Level Adjoint Sensitivity System (1st-LASS). It is important to note that the 1st-LASS is independent of the parameter variations $\delta \alpha$ and therefore needs to be solved only once to determine the (1st-level) adjoint functions ${\psi }^{\left(I\right)}\left(x\right)$ and ${\psi }^{\left(II\right)}\left(y\right)$.

8.

The selection of the boundary conditions for the adjoint functions ${\psi }^{\left(I\right)}\left(x\right)$ and ${\psi }^{\left(II\right)}\left(y\right)$ represented by Equation (38) eliminates the appearance of any unknown values of the variations $\delta u\left(x\right)$ and $\delta v\left(y\right)$ in the bilinear concomitant in Equation (35), reducing it to a residual quantity that contains boundary terms involving only known values of $\delta \alpha$, $u\left(x\right)$, $v\left(y\right)$, ${\psi }^{\left(I\right)}\left(x\right)$, ${\psi }^{\left(II\right)}\left(y\right)$, and $\alpha$. This residual bilinear concomitant will be denoted as ${\left\{{\widehat{P}}^{\left(1\right)}{\left[u\left(x\right);v\left(y\right);{\psi }^{\left(I\right)}\left(x\right),{\psi }^{\left(II\right)}\left(y\right);\alpha ;x,y;\delta \alpha \right]}_{{{}_{\delta \hspace{0.17em}{\mathsf{\Omega }}_x\cup }}_{\delta \hspace{0.17em}{\mathsf{\Omega }}_y}}\right\}}_{\left(e^0\right)}$. In general, this residual bilinear concomitant does not automatically vanish, although it may do so in particular instances. If necessary, this residual bilinear concomitant could be forced to vanish by considering extensions, in the operator sense, of the linear operators $A^{*}\left(u;\alpha \right)\hspace{0.17em}$ and/or $B^{*}\left(v;\alpha \right)$.

9.

Using Equation (35) in conjunction with Equations (36) and (37) in Equation (15) yields the following expression for the indirect-effect term ${\left\{\delta R_x\right\}}^{ind}$:

$$\begin{array}{l}{\left\{\delta R_x\right\}}^{ind}={\displaystyle \sum _{p=1}^{Z_{\alpha }}{\left\{\frac{\partial R_x^{ind}}{\partial {\alpha }_p}\right\}}_{\left(e^0\right)}\delta {\alpha }_p}=\hspace{0.17em}\langle {\left(\begin{array}{c}{\psi }^{\left(I\right)}\left(x\right)\\ {\psi }^{\left(II\right)}\left(y\right)\end{array}\right)}^{†},\hspace{0.17em}\hspace{0.17em}{\left\{\left(\begin{array}{c}{Q}_1^{\left(1\right)}\left(u;\alpha ;\delta \alpha \right)\\ Q_2^{\left(1\right)}\left(v;\alpha ;\delta \alpha \right)\end{array}\hspace{0.17em}\right)\hspace{0.17em}\right\}}_{\left(e^0\right)}\hspace{0.17em}\rangle \\ \hspace{0.17em}\hspace{0.17em}-{\left\{{\widehat{P}}^{\left(1\right)}{\left[u\left(x\right);v\left(y\right);{\psi }^{\left(I\right)}\left(x\right),{\psi }^{\left(II\right)}\left(y\right);\alpha ;x,y;\delta \alpha \right]}_{{{}_{\partial {\mathsf{\Omega }}_x\cup }}_{\partial {\mathsf{\Omega }}_y}}\right\}}_{\left(e^0\right)},\end{array}$$

(39)

where

$$\begin{array}{l}{\left\{\frac{\partial R_x^{ind}}{\partial {\alpha }_p}\right\}}_{\left(e^0\right)}={\displaystyle \prod _{i=1}^{N_x}{\displaystyle \underset{a_i\left({\alpha }^0\right)}{\overset{b_i\left({\alpha }^0\right)}{\int }}d{x}_i}}{\left\{\hspace{0.17em}{\psi }^{\left(I\right)}\left(x\right)\hspace{0.17em}·\hspace{0.17em}\frac{\partial \left[Q^{\left(I\right)}\left(\alpha ;x\right)-N^{\left(I\right)}\left[u\left(x\right);\alpha \right]\right]}{\partial {\alpha }_p}\right\}}_{\left(u^0;{\alpha }^0\right)}\\ +{\displaystyle \prod _{i=1}^{N_y}{\displaystyle \underset{c_i\left({\alpha }^0\right)}{\overset{d_i\left({\alpha }^0\right)}{\int }}d{y}_i}}{\left\{{\psi }^{\left(II\right)}\left(y\right)\hspace{0.17em}·\hspace{0.17em}\frac{\partial \left[Q^{\left(II\right)}\left(\alpha ;y\right)-N^{\left(II\right)}\left[v\left(y\right);\alpha \right]\right]}{\partial {\alpha }_p}\right\}}_{\left(v^0;{\alpha }^0\right)}\\ \hspace{0.17em}\hspace{0.17em}-{\left\{\frac{\partial }{\partial {\alpha }_p}{\widehat{P}}^{\left(1\right)}{\left[u\left(x\right);v\left(y\right);{\psi }^{\left(I\right)}\left(x\right),{\psi }^{\left(II\right)}\left(y\right);\alpha ;x,y;\delta \alpha \right]}_{{{}_{\partial {\mathsf{\Omega }}_x\cup }}_{\partial {\mathsf{\Omega }}_y}}\right\}}_{\left(e^0\right)}.\end{array}$$

(40)

As the expression in Equation (39) indicates, the desired elimination of the unknown variations $\delta u$ and $\delta v$ from the expression of ${\left\{\delta R_x\right\}}^{ind}$ has been accomplished by having used the adjoint functions ${\psi }^{\left(I\right)}\left(x\right)$ and ${\psi }^{\left(II\right)}\left(y\right)$.

The complete expression for computing exactly and efficiently the partial sensitivities $\partial R_x/\partial {\alpha }_p$ of the response $R_x$ with respect to an uncertain parameter ${\alpha }_p$ at the nominal values $\left(u^0;{\alpha }^0\right)$ is obtained by adding Equations (19) and (40) according to Equation (14), i.e.,

$${\left\{\frac{\partial R_x}{\partial {\alpha }_p}\right\}}_{\left(e^0\right)}={\left\{\frac{\partial R_x^{dir}}{\partial {\alpha }_p}+\frac{\partial R_x^{ind}}{\partial {\alpha }_p}\right\}}_{\left(e^0\right)},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}p=1,\dots ,Z_{\alpha }\hspace{0.17em}.$$

(41)

3.2. First-Order Sensitivities of the Critical Point of $R_x$ to the Imprecisely Known Model Parameters, Internal and External Boundaries

The implicit definition provided by Equation (4) for the critical point $\xi \left(\alpha \right)\triangleq {\left[{\xi }_1\left(\alpha \right),\dots ,{\xi }_K\left(\alpha \right)\right]}^{†}$ of $F\left[u\left(x\right);\alpha ;x\right]$ can be written in the following equivalent form:

$${\displaystyle \prod _{i=1}^{N_x}{\displaystyle \underset{a_i\left(\alpha \right)}{\overset{b_i\left(\alpha \right)}{\int }}d{x}_i}\frac{\partial F\left[u\left(x\right);\alpha ;x\right]}{\partial x_k}}\hspace{0.17em}{\displaystyle \prod _{j=1}^K\delta \left[x_j-{\xi }_j\left(\alpha \right)\right]}\hspace{0.17em}=0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}k=1,\dots ,K\le N_x.$$

(42)

Taking the first-order G-variation of Equation (42) at $e^0$ in directions $\delta e$ yields the following expression:

$$\begin{array}{l}\{\frac{d}{d\epsilon }{\displaystyle \prod _{i=1}^{N_x}{\displaystyle \underset{a_i\left({\alpha }^0+\epsilon \delta \alpha \right)}{\overset{b_i\left({\alpha }^0+\epsilon \delta \alpha \right)}{\int }}d{x}_i}}\frac{\partial F\left[u^0\left(x\right)+\epsilon \delta u\left(x\right);{\alpha }^0+\epsilon \delta \alpha ;x\right]}{\partial x_k}\\ \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}\hspace{0.17em}\hspace{0.17em}\times \hspace{0.17em}{{\displaystyle \prod _{m=1}^K\delta \left[x_m-{\xi }_m\left({\alpha }^0+\epsilon \delta \alpha \right)\right]}\}}_{\epsilon =0}=0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}k=1,\dots ,K\le N_x.\end{array}$$

(43)

Carrying out the differentiations indicated in Equation (43) leads to the following relation:

$$\begin{array}{l}0={\displaystyle \prod _{i=1}^{N_x}{\displaystyle \underset{a_i\left({\alpha }^0\right)}{\overset{b_i\left({\alpha }^0\right)}{\int }}d{x}_i}}{\left\{\frac{\partial }{\partial x_k}\left[\frac{\partial F\left(u;\alpha ;x\right)}{\partial \alpha }\delta \alpha +\frac{\partial F\left(u;\alpha ;x\right)}{\partial u}\delta u\right]\right\}}_{\left(u^0;{\alpha }^0\right)}{\displaystyle \prod _{m=1}^K\delta \left[x_m-{\xi }_m\left({\alpha }^0\right)\right]}\hspace{0.17em}\hspace{0.17em}\\ +{\displaystyle \sum _{i=1}^{N_x}{\displaystyle \prod _{j=1,j\ne i}^{N_x}{\displaystyle \underset{a_j\left({\alpha }^0\right)}{\overset{b_j\left({\alpha }^0\right)}{\int }}d{x}_j}}}{\left\{{\left[\frac{\partial F\left[u\left(x\right);\alpha ;x\right]}{\partial x_k}\right]}_{x_i=b_i\left({\alpha }^0\right)}\left[\frac{\partial b_i\left(\alpha \right)}{\partial \alpha }\delta \alpha \right]{\displaystyle \prod _{m=1}^K\delta \left[x_m-{\xi }_m\left(\alpha \right)\right]}\right\}}_{\left(u^0,{\alpha }^0\right)}\\ -{\displaystyle \sum _{i=1}^{N_x}{\displaystyle \prod _{j=1,j\ne i}^{N_x}{\displaystyle \underset{a_j\left({\alpha }^0\right)}{\overset{b_j\left({\alpha }^0\right)}{\int }}d{x}_j}}}{\left\{{\left[\frac{\partial F\left[u\left(x\right);\alpha ;x\right]}{\partial x_k}\right]}_{x_i=a_i\left({\alpha }^0\right)}\left[\frac{\partial a_i\left(\alpha \right)}{\partial {\alpha }_i}\delta \alpha \right]{\displaystyle \prod _{m=1}^K\delta \left[x_m-{\xi }_m\left(\alpha \right)\right]}\right\}}_{\left(u^0,{\alpha }^0\right)}\\ -\hspace{0.17em}{\displaystyle \sum _{j=1}^K{\left\{{\displaystyle \prod _{i=1}^{N_x}{\displaystyle \underset{a_i\left({\alpha }^0\right)}{\overset{b_i\left({\alpha }^0\right)}{\int }}d{x}_j\frac{\partial F\left[u\left(x\right);\alpha ;x\right]}{\partial x_k}{\delta }^{\prime }\left[x_j-{\xi }_j\left(\alpha \right)\right]\left[\frac{\partial {\xi }_j\left(\alpha \right)}{\partial \alpha }\delta \alpha \right]{\displaystyle \prod _{m=1,m\ne j}^K\delta \left[x_m-{\xi }_m\left(\alpha \right)\right]}}}\hspace{0.17em}\right\}}_{\left(u^0;{\alpha }^0\right)}\hspace{0.17em}}.\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}for\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}k=1,\dots ,K\le N_x.\end{array}$$

(44)

After having determined the variations $\delta u$, which are needed in order to evaluate the second term on the right-side of Equation (44), the $K\times K$ system of equations represented by Equation (44) can be solved, in principle, to compute each of the sensitivities $\partial {\xi }_k\left(\alpha \right)/\partial {\alpha }_p$ of the components of the location $\xi \left(\alpha \right)\triangleq {\left[{\xi }_1\left(\alpha \right),\dots ,{\xi }_K\left(\alpha \right)\right]}^{†}$ of the critical point of in phase-space. However, the determination of $\delta u$ requires solving the $\delta \alpha$-dependent 1st-LFSS, cf. Equations (26) and (27), for each parameter variation $\delta {\alpha }_p$, which is prohibitively expensive computationally for large systems with many parameter variations.

To avoid the need for having to solve repeatedly the $\delta \alpha$-dependent 1st-LFSS, cf. Equations (26) and (27), for each parameter variation $\delta {\alpha }_i$, the term containing $\delta u$ in Equation (44) can be expressed in terms of adjoint functions that will be the solutions of a $\delta \alpha$-independent first-level adjoint sensitivity system (1st-LASS), which will be constructed next by following the same procedure as the one that led to Equations (36)–(38). The following well-known relation

$$-{\displaystyle \int g\left(x,y\right){\delta }^{\prime }\left(x-\xi \right)}\hspace{0.17em}dx={\displaystyle \int \frac{\partial g\left(x,y\right)}{\partial x}\delta \left(x-\xi \right)dx}$$

(45)

is used to recast the last term on the right-side of Equation (44) as follows:

$$\begin{array}{l}-\hspace{0.17em}{\displaystyle \sum _{j=1}^K{\left\{{\displaystyle \prod _{i=1}^{N_x}{\displaystyle \underset{a_i\left({\alpha }^0\right)}{\overset{b_i\left({\alpha }^0\right)}{\int }}d{x}_i\frac{\partial F\left[u\left(x\right);\alpha ;x\right]}{\partial x_k}{\delta }^{\prime }\left[x_j-{\xi }_j\left(\alpha \right)\right]\left[\frac{\partial {\xi }_j\left(\alpha \right)}{\partial \alpha }\delta \alpha \right]{\displaystyle \prod _{m=1,m\ne j}^K\delta \left[x_m-{\xi }_m\left(\alpha \right)\right]}}}\hspace{0.17em}\right\}}_{\left(u^0;{\alpha }^0\right)}\hspace{0.17em}}\\ ={\displaystyle \sum _{j=1}^K{\displaystyle \prod _{i=1}^{N_x}{\displaystyle \underset{a_i\left({\alpha }^0\right)}{\overset{b_i\left({\alpha }^0\right)}{\int }}d{x}_i}}\hspace{0.17em}{\left\{\frac{{\partial }^{2}F\left(u;\alpha ;x\right)}{\partial x_k\partial x_j}\frac{d{\xi }_j\left(\alpha \right)}{d\alpha }\delta \alpha \right\}}_{\left(u^0;{\alpha }^0;x_j={\xi }_j\right)}{\displaystyle \prod _{m=1}^K\delta \left[x_m-{\xi }_m\left({\alpha }^0\right)\right]}}\hspace{0.17em}.\end{array}$$

(46)

The relation provided in Equation (45) is also used in order to recast the term containing $\delta u$ in Equation (44) into the following form:

$$\begin{array}{l}-{\displaystyle \prod _{i=1}^{N_x}{\displaystyle \underset{a_i\left({\alpha }^0\right)}{\overset{b_i\left({\alpha }^0\right)}{\int }}d{x}_i}}{\left\{\frac{\partial }{\partial x_k}\left[\frac{\partial F\left(u;\alpha ;x\right)}{\partial u}\delta u\right]\right\}}_{\left(u^0;{\alpha }^0\right)}{\displaystyle \prod _{m=1}^K\delta \left[x_m-{\xi }_m\left({\alpha }^0\right)\right]}\\ =\hspace{0.17em}{\displaystyle \prod _{i=1}^{N_x}{\displaystyle \underset{a_i\left({\alpha }^0\right)}{\overset{b_i\left({\alpha }^0\right)}{\int }}d{x}_i}}{\left\{\frac{\partial F\left(u;\alpha ;x\right)}{\partial u}\delta u\right\}}_{\left(u^0;{\alpha }^0\right)}{\delta }^{\prime }\left[x_k-{\xi }_k\left({\alpha }^0\right)\right]{\displaystyle \prod _{m=1,m\ne k}^K\delta \left[x_m-{\xi }_m\left({\alpha }^0\right)\right]}.\end{array}$$

(47)

The right-side of Equation (47) is similar to the expression of the indirect-effect term defined in Equation (15), with the following correspondence indicated below by the symbol “$\hspace{0.17em}\iff$”:

$$\begin{array}{l}{\left\{\frac{\partial F\left(u;\alpha ;x\right)}{\partial u}\delta u\right\}}_{\left(u^0;{\alpha }^0\right)}{\delta }^{\prime }\left[x_k-{\xi }_k\left({\alpha }^0\right)\right]{\displaystyle \prod _{m=1,m\ne k}^K\delta \left[x_m-{\xi }_m\left({\alpha }^0\right)\right]}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\iff \hspace{0.17em}\hspace{0.17em}{\left\{\frac{\partial F\left(u;\alpha ;x\right)}{\partial u}\delta u\right\}}_{\left(u^0;{\alpha }^0\right)}{\displaystyle \prod _{k=1}^K\delta \left[x_k-{\xi }_k\left({\alpha }^0\right)\right].}\end{array}$$

(48)

Based on the correspondence indicated in Equation (48) and applying the same procedure as already outlined in Section 3.1 for expressing the indirect-effect term defined in Equation (15) leads to the following expression for the functional defined in Equation (47):

$$\begin{array}{l}\hspace{0.17em}{\displaystyle \prod _{i=1}^{N_x}{\displaystyle \underset{a_i\left({\alpha }^0\right)}{\overset{b_i\left({\alpha }^0\right)}{\int }}d{x}_i}}{\left\{\frac{\partial F\left(u;\alpha ;x\right)}{\partial u}\delta u\right\}}_{\left(u^0;{\alpha }^0\right)}{\delta }^{\prime }\left[x_k-{\xi }_k\left({\alpha }^0\right)\right]{\displaystyle \prod _{m=1,m\ne k}^K\delta \left[x_m-{\xi }_m\left({\alpha }^0\right)\right]}\\ =\langle {\left(\begin{array}{c}{\phi }_k^{\left(I\right)}\left(x\right)\\ {\phi }_k^{\left(II\right)}\left(y\right)\end{array}\right)}^{†},\hspace{0.17em}\hspace{0.17em}{\left\{\left(\begin{array}{c}{Q}_1^{\left(1\right)}\left(u;\alpha ;\delta \alpha \right)\\ Q_2^{\left(1\right)}\left(v;\alpha ;\delta \alpha \right)\end{array}\hspace{0.17em}\right)\hspace{0.17em}\right\}}_{\left(e^0\right)}\hspace{0.17em}\rangle \\ \hspace{0.17em}\hspace{0.17em}-{\left\{{\widehat{P}}^{\left(1\right)}{\left[u\left(x\right);v\left(y\right);{\phi }_k^{\left(I\right)}\left(x\right),{\phi }_k^{\left(II\right)}\left(y\right);\alpha ;x,y;\delta \alpha \right]}_{{{}_{\partial {\mathsf{\Omega }}_x\cup }}_{\partial {\mathsf{\Omega }}_y}}\right\}}_{\left(e^0\right)},\hspace{0.17em}\hspace{0.17em}k=1,\dots ,K\le N_x\hspace{0.17em},\end{array}$$

(49)

where the adjoint functions ${\phi }_i^{\left(I\right)}\left(x\right)$ and ${\phi }_i^{\left(II\right)}\left(y\right)$ are the solutions of the following first-level adjoint sensitivity system (1st-LASS), for each $k=1,\dots ,K\le N_x$:

$${\left\{A^{*}\left(u;\alpha \right)\hspace{0.17em}\right\}}_{\left(e^0\right)}{\phi }_k^{\left(I\right)}\left(x\right)={\left\{\frac{\partial F\left(u;\alpha ;x\right)}{\partial u}\right\}}_{\left(u^0;{\alpha }^0\right)}{\delta }^{\prime }\left[x_k-{\xi }_k\left({\alpha }^0\right)\right]{\displaystyle \prod _{m=1,m\ne k}^K\delta \left[x_m-{\xi }_m\left({\alpha }^0\right)\right]},$$

(50)

$${\left\{B^{*}\left(v;\alpha \right)\hspace{0.17em}\right\}}_{\left(e^0\right)}{\phi }_k^{\left(II\right)}\left(y\right)=0$$

(51)

$${\left\{C_A^{\left(1\right)}\left[u\left(x\right);v\left(y\right);{\phi }_k^{\left(I\right)}\left(x\right),{\phi }_k^{\left(II\right)}\left(y\right);\alpha ;x,y\right]\right\}}_{\left(e^0\right)}=0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}x\in \partial {\mathsf{\Omega }}_x,\hspace{0.17em}\hspace{0.17em}y\in \partial {\mathsf{\Omega }}_y$$

(52)

Replacing the results obtained in Equations (46) and (49) in Equation (44) and rearranging the terms in the resulting relations yields the following equation:

$$\begin{array}{l}\langle {\left(\begin{array}{c}{\phi }_k^{\left(I\right)}\left(x\right)\\ {\phi }_k^{\left(II\right)}\left(y\right)\end{array}\right)}^{†},\hspace{0.17em}\hspace{0.17em}{\left\{\left(\begin{array}{c}{Q}_1^{\left(1\right)}\left(u;\alpha ;\delta \alpha \right)\\ Q_2^{\left(1\right)}\left(v;\alpha ;\delta \alpha \right)\end{array}\hspace{0.17em}\right)\hspace{0.17em}\right\}}_{\left(e^0\right)}\hspace{0.17em}\rangle \\ \hspace{0.17em}\hspace{0.17em}-{\left\{{\widehat{P}}^{\left(1\right)}{\left[u\left(x\right);v\left(y\right);{\phi }_k^{\left(I\right)}\left(x\right),{\phi }_k^{\left(II\right)}\left(y\right);\alpha ;x,y;\delta \alpha \right]}_{{{}_{\partial {\mathsf{\Omega }}_x\cup }}_{\partial {\mathsf{\Omega }}_y}}\right\}}_{\left(e^0\right)}\\ -{\displaystyle \prod _{i=1}^{N_x}{\displaystyle \underset{a_i\left({\alpha }^0\right)}{\overset{b_i\left({\alpha }^0\right)}{\int }}d{x}_i}}{\left\{\frac{\partial }{\partial x_k}\left[\frac{\partial F\left(u;\alpha ;x\right)}{\partial \alpha }\delta \alpha \right]\right\}}_{\left(u^0;{\alpha }^0\right)}{\displaystyle \prod _{m=1}^K\delta \left[x_m-{\xi }_m\left({\alpha }^0\right)\right]}\\ -{\displaystyle \sum _{i=1}^{N_x}{\displaystyle \prod _{j=1,j\ne i}^{N_x}{\displaystyle \underset{a_j\left({\alpha }^0\right)}{\overset{b_j\left({\alpha }^0\right)}{\int }}d{x}_j}}}{\left\{{\left[\frac{\partial F\left[u\left(x\right);\alpha ;x\right]}{\partial x_k}\right]}_{x_i=b_i\left({\alpha }^0\right)}\left[\frac{\partial b_i\left(\alpha \right)}{\partial \alpha }\delta \alpha \right]{\displaystyle \prod _{m=1}^K\delta \left[x_m-{\xi }_m\left(\alpha \right)\right]}\right\}}_{\left(u^0,{\alpha }^0\right)}\\ +{\displaystyle \sum _{i=1}^{N_x}{\displaystyle \prod _{j=1,j\ne i}^{N_x}{\displaystyle \underset{a_j\left({\alpha }^0\right)}{\overset{b_j\left({\alpha }^0\right)}{\int }}d{x}_j}}}{\left\{{\left[\frac{\partial F\left[u\left(x\right);\alpha ;x\right]}{\partial x_k}\right]}_{x_i=a_i\left({\alpha }^0\right)}\left[\frac{\partial a_i\left(\alpha \right)}{\partial \alpha }\delta \alpha \right]{\displaystyle \prod _{m=1}^K\delta \left[x_m-{\xi }_m\left(\alpha \right)\right]}\right\}}_{\left(u^0,{\alpha }^0\right)}\\ ={\displaystyle \sum _{j=1}^K{\displaystyle \prod _{i=1}^{N_x}{\displaystyle \underset{a_i\left({\alpha }^0\right)}{\overset{b_i\left({\alpha }^0\right)}{\int }}d{x}_i}}\hspace{0.17em}{\left\{\frac{{\partial }^{2}F\left(u;\alpha ;x\right)}{\partial x_k\partial x_j}\frac{\partial {\xi }_j\left(\alpha \right)}{\partial \alpha }\delta \alpha \right\}}_{\left(u^0;{\alpha }^0;x_j={\xi }_j\right)}{\displaystyle \prod _{m=1}^K\delta \left[x_m-{\xi }_m\left({\alpha }^0\right)\right]}}\hspace{0.17em}\hspace{0.17em}.\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}for\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}k=1,\dots ,K\le N_x.\end{array}$$

(53)

Replacing $\left(\partial /\partial \alpha \right)\delta \alpha ={\displaystyle \sum _{p=1}^{Z_{\alpha }}\left(\partial /\partial {\alpha }_p\right)\delta {\alpha }_p}$ in Equation (53) and collecting the expressions multiplying the parameter variations $\delta {\alpha }_p$ in the resulting relation leads to a $K\times K$ system of equations for determining the sensitivities $\left(\partial {\xi }_k/\partial {\alpha }_p\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}k=1,\dots ,K\le N_x$, for each parameter ${\alpha }_p$, $p=1,\dots ,Z_{\alpha }$. This $K\times K$ system of equations can be written in matrix form as follows:

$$\mathrm{H}{\theta }^{\left(p\right)}=s^{\left(p\right)},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{H}\triangleq \left(\hspace{0.17em}\begin{array}{ccc}{h}_{11}& ·& h_{1K}\\ ·& ·& ·\\ h_{K1}& ·& h_{KK}\end{array}\right),\hspace{0.17em}\hspace{0.17em}{\theta }^{\left(p\right)}\triangleq \left(\hspace{0.17em}\hspace{0.17em}\begin{array}{c}{\theta }_1^{\left(p\right)}\\ ·\\ {\theta }_K^{\left(p\right)}\end{array}\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{s}^{\left(p\right)}\triangleq \left(\hspace{0.17em}\hspace{0.17em}\begin{array}{c}{s}_1^{\left(p\right)}\\ ·\\ s_K^{\left(p\right)}\end{array}\right),$$

(54)

where:

(a)

the matrix $\mathrm{H}$ denotes the Hessian of the function $F\left(u;\alpha ;x\right)$ computed at the nominal value $\xi \left({\alpha }^0\right)$ of the critical point of $F\left(u;\alpha ;x\right)$; the components of $\mathrm{H}$ are defined as follows:

$$h_{k\ell }\triangleq {\displaystyle \prod _{i=1}^{N_x}{\displaystyle \underset{a_i\left({\alpha }^0\right)}{\overset{b_i\left({\alpha }^0\right)}{\int }}d{x}_i}}\hspace{0.17em}{\left\{\frac{{\partial }^{2}F\left(u;\alpha ;x\right)}{\partial x_k\partial x_{\ell }}{\displaystyle \prod _{m=1}^K\delta \left[x_m-{\xi }_m\left(\alpha \right)\right]}\right\}}_{\left(u^0;{\alpha }^0\right)};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}k,\ell =1,\dots K\le N_x;$$

(55)

(b)

the components of the vectors ${\theta }^{\left(p\right)}$ and $s^{\left(p\right)}$, respectively, are defined as follows:

$${\theta }_k^{\left(p\right)}\triangleq \hspace{0.17em}\partial {\xi }_k\left(\alpha \right)/\partial {\alpha }_p\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}k=1,\dots ,K\le N_x;\hspace{0.17em}\hspace{0.17em}p=1,\dots ,Z_{\alpha }$$

(56)

$$\begin{array}{l}{s}_k^{\left(p\right)}={\displaystyle \prod _{i=1}^{N_x}{\displaystyle \underset{a_i\left({\alpha }^0\right)}{\overset{b_i\left({\alpha }^0\right)}{\int }}d{x}_i}}{\left\{\hspace{0.17em}{\phi }_k^{\left(I\right)}\left(x\right)\hspace{0.17em}·\hspace{0.17em}\frac{\partial \left[Q^{\left(I\right)}\left(\alpha ;x\right)-N^{\left(I\right)}\left[u\left(x\right);\alpha \right]\right]}{\partial {\alpha }_p}\right\}}_{\left(u^0;{\alpha }^0\right)}\\ +{\displaystyle \prod _{i=1}^{N_y}{\displaystyle \underset{c_i\left({\alpha }^0\right)}{\overset{d_i\left({\alpha }^0\right)}{\int }}d{y}_i}}{\left\{{\phi }_k^{\left(II\right)}\left(y\right)\hspace{0.17em}·\hspace{0.17em}\frac{\partial \left[Q^{\left(II\right)}\left(\alpha ;y\right)-N^{\left(II\right)}\left[v\left(y\right);\alpha \right]\right]}{\partial {\alpha }_p}\right\}}_{\left(v^0;{\alpha }^0\right)}\\ \hspace{0.17em}-{\left\{\frac{\partial }{\partial {\alpha }_p}{\widehat{P}}^{\left(1\right)}{\left[u\left(x\right);v\left(y\right);{\phi }_k^{\left(I\right)}\left(x\right),{\phi }_k^{\left(II\right)}\left(y\right);\alpha ;x,y;\delta \alpha \right]}_{{{}_{\partial {\mathsf{\Omega }}_x\cup }}_{\partial {\mathsf{\Omega }}_y}}\right\}}_{\left(e^0\right)}\\ -{\displaystyle \prod _{i=1}^{N_x}{\displaystyle \underset{a_i\left({\alpha }^0\right)}{\overset{b_i\left({\alpha }^0\right)}{\int }}d{x}_i}}{\left\{\frac{\partial }{\partial x_k}\frac{\partial F\left(u;\alpha ;x\right)}{\partial {\alpha }_p}\right\}}_{\left(u^0;{\alpha }^0\right)}{\displaystyle \prod _{m=1}^K\delta \left[x_m-{\xi }_m\left({\alpha }^0\right)\right]}\\ -{\displaystyle \sum _{i=1}^{N_x}{\displaystyle \prod _{j=1,j\ne i}^{N_x}{\displaystyle \underset{a_j\left({\alpha }^0\right)}{\overset{b_j\left({\alpha }^0\right)}{\int }}d{x}_j}}}{\left\{{\left[\frac{\partial F\left[u\left(x\right);\alpha ;x\right]}{\partial x_k}\frac{\partial b_i\left(\alpha \right)}{\partial {\alpha }_p}\right]}_{x_i=b_i\left({\alpha }^0\right)}{\displaystyle \prod _{m=1}^K\delta \left[x_m-{\xi }_m\left(\alpha \right)\right]}\right\}}_{\left(u^0,{\alpha }^0\right)}\\ +{\displaystyle \sum _{i=1}^{N_x}{\displaystyle \prod _{j=1,j\ne i}^{N_x}{\displaystyle \underset{a_j\left({\alpha }^0\right)}{\overset{b_j\left({\alpha }^0\right)}{\int }}d{x}_j}}}{\left\{{\left[\frac{\partial F\left[u\left(x\right);\alpha ;x\right]}{\partial x_k}\frac{\partial a_i\left(\alpha \right)}{\partial {\alpha }_p}\right]}_{x_i=a_i\left({\alpha }^0\right)}{\displaystyle \prod _{m=1}^K\delta \left[x_m-{\xi }_m\left(\alpha \right)\right]}\right\}}_{\left(u^0,{\alpha }^0\right)}\\ k=1,\dots K\le N_x;\hspace{0.17em}\hspace{0.17em}p=1,\dots ,Z_{\alpha }\hspace{0.17em}.\end{array}$$

(57)

Evidently, the sensitivities ${\theta }_k^{\left(p\right)}\triangleq \hspace{0.17em}\partial {\xi }_k\left(\alpha \right)/\partial {\alpha }_p\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}k=1,\dots ,K\le N_x$ of the critical point $\xi \left(\alpha \right)\triangleq {\left[{\xi }_1\left(\alpha \right),\dots ,{\xi }_K\left(\alpha \right)\right]}^{†}$ with respect to the uncertain parameter ${\alpha }_p$ are computed at the nominal value $\xi \left({\alpha }^0\right)$ by solving Equation (54) to obtain:

$${\theta }^{\left(p\right)}={\mathrm{H}}^{-1}{s}^{\left(p\right)},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}p=1,\dots ,Z_{\alpha }$$

(58)

Examining the structure of Equation (58) reveals the following important characteristics:

• The matrix $\mathrm{H}$, which is the only matrix that needs to be inverted, does not depend on the parameters variations and has very small dimensions (its dimensions are at most equal to the number of the independent variables that characterize “subsystem I”). Therefore, $\mathrm{H}$ can be inverted once and for all, and its inverse can be easily stored for using it as many times (i.e., $Z_{\alpha }$-times) as there are system parameters.
• Before computing the right-side of Equation (58), it is necessary to determine the adjoint functions ${\phi }_k^{\left(I\right)}\left(x\right)$ and ${\phi }_k^{\left(II\right)}\left(y\right)$ by solving the 1st-LASS comprising Equations (50)–(52). Solving the 1st-LASS represents a “large-scale” computation, of the same size as solving the original system represented by Equations (1)–(3). Solving the 1st-LASS is expected to be less intensive computationally than solving the original system, since Equations (50)–(52) are linear in the dependent variables, whereas Equations (1)–(3) are nonlinear in the dependent variables. It is also important to note that the 1st-LASS is independent of parameter variations. The 1st-LASS needs to be solved once for each component of the critical point in the phase-space of independent variables, which is a relatively small number (e.g., the 1st-LASS need to be solved once for one-dimensional systems, at most twice for two-dimensional systems, at most three times for three-dimensional systems, etc.).
• Only the right-side of Equation (58) depends on derivatives of various quantities with respect to the parameters, so it would need to be computed anew for each parameter. However, these computations involve only integrals over various quantities. These integrals can be computed efficiently and inexpensively using quadrature formulas (as opposed to needing to solve large-scale differential equations).

The derivations of the first-order adjoint sensitivity system needed for the exact and efficient computation of the first-order sensitivities of the response $R_y$ and of the critical point $\xi \left(\alpha \right)\triangleq {\left[{\xi }_1\left(\alpha \right),\dots ,{\xi }_I\left(\alpha \right)\right]}^{†},\hspace{0.17em}\hspace{0.17em}I\le N_x$ at which the response $R_y$ is defined (with respect to the model’s uncertain parameters, including boundary and interface conditions) are similar to those presented in Section 3.1 and are therefore presented in Appendix B, for convenient reference.

4. Discussion

This work has presented the mathematical framework of the extended “first-order comprehensive adjoint sensitivity analysis methodology for critical points” (1st-CASAM-CP), which enables the exact and efficient computation of the first-order sensitivities of responses defined at critical points of coupled nonlinear systems characterized by imprecisely known parameters, interfaces, and boundaries. It has been shown that the computation of the first-order response sensitivities of the magnitude of a response defined at a critical point of a function in the phase-space of the systems’ independent variables requires a single large-scale computations of the 1st-LASS which corresponds to the magnitude of the critical point. Solving the 1st-LASS represents a “large-scale” computation which is not more extensive than solving the original coupled systems, since the 1st-LASS is linear in the dependent variables, whereas the original coupled systems nonlinear in the dependent variables. Furthermore, the computation of the sensitivities of the location in phase-space of each critical point requires solving one 1st-LASS for each of the components, in the phase-space of independent variables, of the respective critical point. Noteworthy, the 1st-LASS needed to solve in order to obtain the respective adjoint sensitivity functions have the same operators on their right-sides; only the sources on the left-sides of the respective 1st-LASS differ from each other. Hence, the same solver can be used to solve all of the requisite 1st-LASS (as opposed to needing to develop a different solver for each 1st-LASS) to compute all of the respective adjoint functions. Of course, the computational effort and resources needed to solve the requisite 1st-LASS are significantly less expensive than attempting the compute the sensitivities of critical points by brute-force forward re-computations using the model of the coupled nonlinear systems with altered parameter values in conjunction with finite-difference schemes, which would require at least as many large-scale computations as there are uncertain parameters in the coupled nonlinear systems. Computing higher-order sensitivities

In the accompanying work [13], the unique features of the 1st-CASAM-CP presented in this work will be illustrated by applying it to a nuclear reactor thermal-hydraulics safety benchmark [14,15,16] which has been used for the verification of the “FLUENT Adjoint Solver” [18].