Fifth-Order Comprehensive Adjoint Sensitivity Analysis Methodology for Nonlinear Systems (5th-CASAM-N): I. Mathematical Framework

Dan Gabriel Cacuci

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

DOI: 10.4236/ajcm.2022.121005                

Abstract

This work presents the mathematical framework of the “Fifth-Order Comprehensive Adjoint Sensitivity Analysis Methodology for Nonlinear Systems (5^th-CASAM-N),” which generalizes and extends all of the previous works performed to date on this subject. The 5^th-CASAM-N enables the exact and efficient computation of all sensitivities, up to and including fifth-order, of model responses to uncertain model parameters and uncertain boundaries of the system’s domain of definition, thus enabling, inter alia, the quantification of uncertainties stemming from manufacturing tolerances. The 5^th-CASAM-N provides a fundamental step towards overcoming the curse of dimensionality in sensitivity and uncertainty analysis.

Keywords

Nonlinear Models of Physical Systems, High-Order Sensitivity Analysis, 1st-Order Sensitivities, 2nd-Order Sensitivities, 3rd-Order Sensitivities, 4th-Order Sensitivities, 5th-Order Sensitivities

1. Introduction

This work presents the mathematical framework of the “Fifth-Order Comprehensive Adjoint Sensitivity Analysis Methodology for Nonlinear Systems” abbreviated as “5^th-CASAM-N.” The 5^th-CASAM-N generalizes the previous mathematical works on this topic, which stems from the framework of the first-order adjoint sensitivity analysis methodology for generic nonlinear systems established in [1]. Numerous specific applications using adjoint functions for the deterministic computation of first-order sensitivities of scalar-valued model responses to model parameters have been published over the years, as discussed and referenced in the book by Cacuci [2]. The history regarding the deterministic computation of second-order sensitivities of model responses to model parameters reveals that although many particular applications which used specifically-computed 2^nd-order response sensitivities have been published over the years, the generic mathematical framework of the 2^nd-order adjoint sensitivity analysis methodology was presented in [3] for linear systems and in [4] nonlinear systems. The mathematical frameworks and notable applications of the 2^nd-order adjoint sensitivity analysis methodology for both linear and nonlinear systems are presented and referenced in the book by Cacuci [5].

The largest application to date of the 2^nd-order adjoint sensitivity analysis methodology for linear systems to date has been presented in [6] - [11] for the polyethylene-reflected plutonium (acronym: PERP) OECD/NEA reactor physics benchmark [12]. The numerical model [6] - [11] of the PERP benchmark comprises 21,976 uncertain parameters; 7477 of these uncertain model parameters have nonzero nominal values, as follows: 180 microscopic total cross sections; 7101 microscopic scattering sections; 60 microscopic fission cross sections; 60 parameters that characterize the average number of neutron per fission; 60 parameters that characterize the fission spectrum; 10 parameters that characterize the fission source; and 6 parameters that characterize the isotope number densities. Applying the second-order adjoint sensitivity analysis methodology made it possible to compute exactly the 7477 non-zero first-order sensitivities and the (7477)² second-order sensitivities of the PERP benchmark’s leakage response with respect to the benchmark’s imprecisely known parameters. The results of these computations have indicated that 13 first-order sensitivities attain values between 1.0 and 10.0, while 126 second-order relative sensitivities have values greater than 10.0, and 1853 second-order relative sensitivities have values between 1.0 and 10.0. These results were contrary to the previously held belief that 2^nd-order relative sensitivities are smaller than 1^st-order relative sensitivities for neutron transport models such as the PERP benchmark’s model.

The finding that many 2^nd-order sensitivities were significantly larger than the 1^st-order ones has motivated the subsequent computation of the 3^rd-order sensitivities of the leakage response with respect to the PERP benchmark’s total cross sections. The largest 3^rd-order sensitivities were computed in [13] [14] [15] by applying the 3^rd-order adjoint sensitivity analysis methodology. It has been found that the number of 3^rd-order mixed relative sensitivities that have large values (and are therefore important) is significantly greater than the number of important 2^nd- and 1^st-order sensitivities. For example, the first-order relative sensitivity of the benchmark’s leakage response with respect to the total microscopic cross section of hydrogen in the lowest energy group, denoted as $S^{\left(1\right)}\left({\sigma }_{t,6}^{30}\right)=-9.366$, was the largest of the 7477 first-order sensitivities. But the unmixed second-order and third-order relative sensitivities of the benchmark’s leakage response with respect to the same parameter had the following values: $S^{\left(2\right)}\left({\sigma }_{t,6}^{30},{\sigma }_{t,6}^{30}\right)=429.6$ and $S^{\left(3\right)}\left({\sigma }_{t,1}^{g=30},{\sigma }_{t,6}^{g^{\prime }=30},{\sigma }_{t,6}^{g^{″}=30}\right)=-1.88\times {10}^5$, respectively. The results obtained in [13] [14] [15] have motivated the development in [16] of the 4^th-Order Comprehensive Adjoint Sensitivity Analysis Methodology for Linear Systems (4^th-CASAM-L), which was applied in [17] [18] [19] [20] to the PERP benchmark for computing exactly and efficiently the most important 4^th-order sensitivities of the benchmark’s total leakage response with respect to the benchmark’s 180 microscopic total cross sections, which include 180 4^th-order unmixed sensitivities and 360 4^th-order mixed sensitivities corresponding to the largest 3^rd-order ones. The numerical results obtained in [17] [18] [19] [20] indicated, in particular, that the largest overall 4^th-order relative sensitivity was the 4^th-order relative sensitivity of the benchmark’s leakage response with respect to the total microscopic cross section of hydrogen in the lowest energy group, namely $S^{\left(4\right)}\left({\sigma }_{t,6}^{g=30},{\sigma }_{t,6}^{g=30},{\sigma }_{t,6}^{g=30},{\sigma }_{t,6}^{g=30}\right)=2.720\times {10}^6$. This value is around 291,000 times larger than the 1^st-order relative sensitivity $S^{\left(1\right)}\left({\sigma }_{t,6}^{30}\right)=-9.366$, 6350 times larger than the 2^nd-order relative sensitivity $S^{\left(2\right)}\left({\sigma }_{t,6}^{30},{\sigma }_{t,6}^{30}\right)=429.6$, and 90 times larger than the 3^rd-order relative sensitivity $S^{\left(3\right)}\left({\sigma }_{t,1}^{g=30},{\sigma }_{t,6}^{g^{\prime }=30},{\sigma }_{t,6}^{g^{″}=30}\right)=-1.88\times {10}^5$. The results obtained in [6] - [19] have indicated that higher-order sensitivities cannot be simply ignored out of hand but must be computed and their impact (e.g., on uncertainty analysis) must be evaluated in the context of the application under consideration.

The results obtained in [6] - [20] have also motivated the development of the general mathematical framework for computing exactly and efficiently arbitrarily-high-order sensitivities of model responses to model parameters. Since only linear systems admit bona-fide adjoint operators (in contradistinction to nonlinear operators, which do not admit adjoint operators), Cacuci has developed [21] “The n^th-Order Comprehensive Adjoint Sensitivity Analysis Methodology for Response-Coupled Forward/Adjoint Linear Systems” (n^th-CASAM-L), which enables the exact and efficient computation of sensitivities, of any order, of model responses to model parameters, including imprecisely known domain boundaries, thus enabling the quantification of uncertainties stemming, among other factors, from manufacturing tolerances. The n^th-CASAM-L overcomes the “curse of dimensionality” [22] in sensitivity and uncertainty analysis, as detailed in [23].

For nonlinear systems, the first general methodology which also enabled the exact and efficient computation of model response sensitivities to uncertain domains of definition of the model’s independent variables (in addition to response sensitivities to model parameters) were the works [24] [25] [26] [27] [28] on the “1^st-Order Comprehensive Adjoint Sensitivity Analysis Methodology for Nonlinear Systems.” The “Fifth-Order Comprehensive Adjoint Sensitivity Analysis Methodology for Nonlinear Systems (5^th-CASAM-N)” to be presented in this work generalizes and extends the mathematical framework presented in [24] [25] [26] in order to enable the computation of 5^th-order sensitivities. This work is structured as follows: Section 2 presents the mathematical framework of the 5^th-CASAM-N, which builds on the lower-order adjoints sensitivity analysis methodologies [24] [25] [26] [27] [28]. The significance of the potential applications of the innovative 5^th-CASAM-N are discussed in Section 3. A paradigm illustrative application to a Bernoulli model with uncertain parameters and boundaries will be presented in an accompanying work [29].

2. The Fifth-Order Comprehensive Adjoint Sensitivity Analysis Methodology for Nonlinear Systems (5^th-CASAM-N) Methodology

The mathematical framework of the “Fifth-Order Comprehensive Adjoint Sensitivity Analysis Methodology for Nonlinear Systems (5^th-CASAM-N)” builds upon the framework of the 4^th-CASAM-N, which is presented in [27]. It is therefore necessary to recall the mathematical framework underlying the 4^th-CASAM-N, which will be summarized in this Section. Matrices will be denoted using capital bold letters while vectors will be denoted using either capital or lower-case bold letters. The symbol “ $\triangleq$ ” will be used to denote “is defined as” or “is by definition equal to.” Transposition will be indicated by a dagger $(†)$ superscript.

The computational model of a physical system comprises equations that relate the system’s state variables to system’s independent variables and parameters, which are considered to be afflicted by uncertainties. The information customarily available about the model parameters comprises their nominal (expected/mean) values and, possibly, higher-order moments or cumulants (i.e., variance/covariances, skewness, kurtosis), which are usually determined from evaluation processes external to the physical system under consideration. Occasionally, only lower and/or upper bounds may be known for some model parameters. The parameters will be denoted as ${\alpha }_1,\cdots ,{\alpha }_{TP}$, where the subscript TP indicates “total number of imprecisely known parameters.” Without loss of generality, the imprecisely known model parameters can be considered to be real-valued scalars which are considered to be the components of a “vector of parameters” denoted as $\alpha \triangleq {\left({\alpha }_1,\cdots ,{\alpha }_{TP}\right)}^{†}\in {ℝ}^{TP}$, where ${ℝ}^{TP}$ denotes the TP-dimensional subset of the set of real scalars. The components of $\alpha \in {ℝ}^{TP}$ are considered to include imprecisely known geometrical parameters that characterize the physical system’s boundaries in the phase-space of the model’s independent variables. The nominal parameter values will be denoted as ${\alpha }^0\triangleq {\left[{\alpha }_1^0,\cdots ,{\alpha }_i^0,\cdots ,{\alpha }_{TP}^0\right]}^{†}$ ; the superscript “0” will be used throughout this work to denote “nominal values.”

The generic nonlinear model is considered to comprise TI independent variables which will be denoted as $x_i,i=1,\cdots ,TI$, where the sub/superscript “TI” denotes the “total number of independent variables.” The independent variables are considered to be components of a TI-dimensional column vector denoted as $x\triangleq {\left(x_1,\cdots ,x_{TI}\right)}^{†}\in {ℝ}^{TI}$. The vector $x\in {ℝ}^{TI}$ is considered to be defined on a phase-space domain, denoted as $\Omega \left(\alpha \right)$ and defined as follows: $\Omega \left(\alpha \right)\triangleq \left\{-\infty \le {\lambda }_i\left(\alpha \right)\le x_i\le {\omega }_i\left(\alpha \right)\le \infty ;i=1,\cdots ,TI\right\}$. The lower boundary-point of an independent variable is denoted as ${\lambda }_i\left(\alpha \right)$ and the corresponding upper boundary-point is denoted as ${\omega }_i\left(\alpha \right)$. The boundary of $\Omega \left(\alpha \right)$, which will be denoted as $\partial \Omega \left(\alpha \right)$, comprises the set of all of the endpoints ${\lambda }_i\left(\alpha \right),{\omega }_i\left(\alpha \right),i=1,\cdots ,TI$ of the respective intervals on which the components of $x$ are defined, i.e., $\partial \Omega \left(\alpha \right)\triangleq \left\{{\lambda }_i\left(\alpha \right)\cup {\omega }_i\left(\alpha \right),i=1,\cdots ,TI\right\}$. The boundary $\Omega \left(\alpha \right)$ is also considered to be imprecisely known since it may depend on both geometrical parameters and material properties. For example, the “extrapolated boundary” in models based on diffusion theory depends both on the imprecisely known physical dimensions of the problem’s domain and also on the medium’s properties (atomic number densities, microscopic transport cross sections, etc.).

The model of a nonlinear physical system comprises coupled equations which can be represented in operator form as follows:

$N\left[u\left(x\right),\alpha \right]=Q\left(x,\alpha \right),x\in {\Omega }_x\left(\alpha \right).$ (1)

The quantities which appear in Equation (1) are defined as follows: 1) $u\left(x\right)\triangleq {\left[u_1\left(x\right),\cdots ,u_{TD}\left(x\right)\right]}^{†}$ is a TD-dimensional column vector of dependent variables (also called “state functions”), where “TD” denotes “total number of dependent variables;” 2) $N\left[u\left(x\right);\alpha \right]\triangleq {\left[N_1\left(u;\alpha \right),\cdots ,N_{TD}\left(u;\alpha \right)\right]}^{†}$ denotes a TD-dimensional column vector, having components $N_i\left(u;\alpha \right),i=1,\cdots ,TD$, which are operators that act on the dependent variables $u\left(x\right)$, the independent variables $x$ and the model parameters $\alpha$ ; 3) $Q\left(x,\alpha \right)\triangleq {\left[q_1\left(x;\alpha \right),\cdots ,q_{TD}\left(x;\alpha \right)\right]}^{†}$ is a TD-dimensional column vector which represents inhomogeneous source terms, which usually depend nonlinearly on the uncertain parameters $\alpha$ ; 4) since the right-side of Equation (1) may contain “generalized functions/functionals” (e.g., Dirac-distributions and derivatives thereof), the equalities in this work are considered to hold in the distributional (“weak”) sense.

When differential operators appear in Equation (1), their domains of definition must be specified by providing boundary and/or initial conditions. Mathematically, these boundaries and/or initial conditions can be represented in operator form as follows:

$B\left[u\left(x\right);\alpha ;x\right]-C\left(x,\alpha \right)=0,x\in \partial {\Omega }_x\left(\alpha \right).$ (2)

where the column vector $0$ has TD components, all of which are zero. The components $B_i\left(u;\alpha \right),i=1,\cdots ,TD$ of $B\left(u;\alpha \right)\triangleq {\left[B_1\left(u;\alpha \right),\cdots ,B_{TD}\left(u;\alpha \right)\right]}^{†}$ are nonlinear operators in $u\left(x\right)$ and $\alpha$, which are defined on the boundary $\partial {\Omega }_x\left(\alpha \right)$ of the model’s domain ${\Omega }_x\left(\alpha \right)$. The components $C_i\left(x;\alpha \right),i=1,\cdots ,TD$ of $C\left(x;\alpha \right)\triangleq {\left[C_1\left(x;\alpha \right),\cdots ,C_{TD}\left(x;\alpha \right)\right]}^{†}$ comprise inhomogeneous boundary sources which are nonlinear functions of $\alpha$.

The model’s nominal solution, denoted as $u^0\left(x\right)$, is obtained by solving Equations (1) and (2) at the nominal parameter values, namely:

$N\left[u^0\left(x\right);{\alpha }^0\right]=Q\left(x,{\alpha }^0\right),x\in {\Omega }_x$, (3)

$B\left[u^0\left(x\right);{\alpha }^0;x\right]-C\left(x,{\alpha }^0\right)=0,x\in \partial {\Omega }_x\left({\alpha }^0\right).$ (4)

The model response considered in this work is a nonlinear functional of the model’s state functions and parameters which can be generically represented as follows:

$R\left[u\left(x\right);\alpha \right]\triangleq {\displaystyle \underset{{\lambda }_1\left(\alpha \right)}{\overset{{\omega }_1\left(\alpha \right)}{\int }}\cdots {\displaystyle \underset{{\lambda }_{TI}\left(\alpha \right)}{\overset{{\omega }_{TI}\left(\alpha \right)}{\int }}S\left[u\left(x\right);\alpha ;x\right]\text{d}{x}_1\cdots \text{d}{x}_{TI}}}$, (5)

where $S\left[u\left(x\right);\alpha \right]$ is suitably differentiable nonlinear function of $u\left(x\right)$ and of $\alpha$. Noteworthy, the components of $\alpha$ also include parameters that may occur just in the definition of the response under consideration, in addition to the parameters that appear in Equations (1) and (2). Since the system domain’s boundary, $\partial \Omega \left(\alpha \right)$, is considered to be subject to uncertainties (e.g., stemming from manufacturing uncertainties), the model response $R\left[u\left(x\right);\alpha \right]$ will also be affected by the uncertainties that affect the endpoints ${\lambda }_i\left(\alpha \right),{\omega }_i\left(\alpha \right),i=1,\cdots ,TI$, of $\partial \Omega \left(\alpha \right)$.

2.1. The First-Order Comprehensive Adjoint Sensitivity Analysis Methodology for Nonlinear Systems (1^st-CASAM-N)

The model and boundary parameters $\alpha$ are considered to be uncertain quantities, having unknown true values. The nominal (or mean) parameter vales ${\alpha }^0$ are considered to be known, and these will differ from the true values by quantities denoted as $\delta \alpha \triangleq \left(\delta {\alpha }_1,\cdots ,\delta {\alpha }_{TP}\right)$, where $\delta {\alpha }_i\triangleq {\alpha }_i-{\alpha }_i^0$. Since the forward state functions $u\left(x\right)$ are related to the model and boundary parameters $\alpha$ through Equations (1) and (2), it follows that the variations $\delta \alpha$ in the model and boundary parameters will cause corresponding variations $v^{\left(1\right)}\left(x\right)\triangleq {\left[\delta u_1\left(x\right),\cdots ,\delta u_{TD}\left(x\right)\right]}^{†}$ around the nominal solution $u^0\left(x\right)$ in the forward state functions. In turn, the variations $\delta \alpha$ and $v^{\left(1\right)}\left(x\right)$ will induce variations in the system’s response.

As shown in [1], the 1^st-order sensitivities of a model response $R\left(e\right)$, where $e\triangleq \left(\alpha ,u\right)\in E$, with respect to variations $h\triangleq \left(\delta \alpha ,v^{\left(1\right)}\right)$ in the model parameters and state functions in a neighborhood around the nominal functions and parameter values $e^0\triangleq \left({\alpha }^0,u^0\right)\in E$, are obtained by determining the 1^st-order Gateaux- (G-) variation of the response. The 1^st-order Gateaux- (G-) variation, denoted as $\delta R\left(e^0;h\right)$, of the response will exist and will be linear in $h\triangleq \left(\delta \alpha ,v^{\left(1\right)}\right)$ if and only if the following two conditions are satisfied by $R\left(e\right)$ :

1) $R\left(e\right)$ satisfies a weak Lipschitz condition at $e^0$, which has the following form:

$\Vert R\left(e^0+\epsilon h\right)-R\left(e^0\right)\Vert \le k\Vert \epsilon \left(e^0\right)\Vert ,k<\infty$, (6)

2) $R\left(e\right)$ satisfies the following condition:

$\begin{array}{l}R\left(e^0+\epsilon h_1+\epsilon h_2\right)-R\left(e^0+\epsilon h_1\right)\\ -R\left(e^0+\epsilon h_2\right)+R\left(e^0\right)=o\left(\epsilon \right);h_1\text{,}{h}_2\in E\text{;}\epsilon \in F.\end{array}$ (7)

In Equation (7), the symbol $F$ denotes the underlying field of scalars. Numerical methods (e.g., Newton’s method and variants thereof) for solving Equations (1) and (2) also require the existence of the first-order G-derivatives of original model equations. Therefore, the conditions provided in Equations (6) and (7) are henceforth considered to be satisfied by the model responses and also by the operators underlying the physical system modeled by Equations (1) and (2). When the response $R\left(e\right)$ satisfies the conditions provided in Equations (6) and (7), the 1^st-order G-differential $\delta R\left(e^0;h\right)$ can be written as follows:

$\begin{array}{l}\delta R\left(e^0;h\right)\equiv {\left\{\delta R\left[u\left(x\right);\alpha ;v^{\left(1\right)}\left(x\right);\delta \alpha \right]\right\}}_{{\alpha }^0}\\ \triangleq {\left\{\delta R\left[u\left(x\right);\alpha ;\delta \alpha \right]\right\}}_{dir}+{\left\{\delta R\left[u\left(x\right);\alpha ;v^{\left(1\right)}\left(x\right)\right]\right\}}_{ind}.\end{array}$ (8)

In Equation (8), the “direct-effect” term ${\left\{\delta R\left[u\left(x\right);\alpha ;\delta \alpha \right]\right\}}_{dir}$ comprises only dependencies on $\delta \alpha$ and is defined as follows:

${\left\{\delta R\left[u\left(x\right);\alpha ;\delta \alpha \right]\right\}}_{dir}\triangleq {\left\{\frac{\partial R\left(u;\alpha \right)}{\partial \alpha }\right\}}_{{\alpha }^0}\delta \alpha \triangleq {\displaystyle \underset{j_1=1}{\overset{TP}{\sum }}{\left\{R^{\left(1\right)}\left[j_1;u\left(x\right);\alpha \right]\right\}}_{dir}\delta {\alpha }_{j_1}},$ (9)

where $\partial R\left(u;\alpha \right)/\partial \alpha$ denotes the partial G-derivatives of $R\left(e\right)$ with respect to $\alpha$, evaluated at the nominal parameter values, and where the following definitions were used:

$\frac{\partial \left[\right]}{\partial \alpha }\delta \alpha \triangleq {\displaystyle \underset{i=1}{\overset{TP}{\sum }}\frac{\partial \left[\right]}{\partial {\alpha }_i}\delta {\alpha }_i}$. (10)

$\begin{array}{l}{\left\{R^{\left(1\right)}\left[j_1;u\left(x\right);\alpha \right]\right\}}_{dir}\triangleq {\left\{{\displaystyle \underset{{\lambda }_1\left(\alpha \right)}{\overset{{\omega }_1\left(\alpha \right)}{\int }}\cdots {\displaystyle \underset{{\lambda }_{TI}\left(\alpha \right)}{\overset{{\omega }_{TI}\left(\alpha \right)}{\int }}\frac{\partial S\left(u;\alpha ;\alpha \right)}{\partial {\alpha }_{j_1}}\text{d}{x}_1\cdots \text{d}{x}_{TI}}}\right\}}_{{\alpha }^0}\\ +{\displaystyle \underset{j=1}{\overset{TI}{\sum }}{\left\{{\displaystyle \underset{{\lambda }_1\left(\alpha \right)}{\overset{{\omega }_1\left(\alpha \right)}{\int }}\cdots {\displaystyle \underset{{\lambda }_{j-1}\left(\alpha \right)}{\overset{{\omega }_{j-1}\left(\alpha \right)}{\int }}{\displaystyle \underset{{\lambda }_{j+1}\left(\alpha \right)}{\overset{{\omega }_{j+1}\left(\alpha \right)}{\int }}\cdots {\displaystyle \underset{{\lambda }_{TI}\left(\alpha \right)}{\overset{{\omega }_{TI}\left(\alpha \right)}{\int }}S\left[u\left(x_1,.,{\omega }_j\left(\alpha \right),.,x_{N_x}\right);\alpha \right]\frac{\partial {\omega }_j\left(\alpha \right)}{\partial {\alpha }_{j_1}}\text{d}{x}_1\cdots \text{d}{x}_{TI}}}}}\right\}}_{{\alpha }^0}}\\ -{\displaystyle \underset{j=1}{\overset{TI}{\sum }}{\left\{{\displaystyle \underset{{\lambda }_1\left(\alpha \right)}{\overset{{\omega }_1\left(\alpha \right)}{\int }}\cdots {\displaystyle \underset{{\lambda }_{j-1}\left(\alpha \right)}{\overset{{\omega }_{j-1}\left(\alpha \right)}{\int }}{\displaystyle \underset{{\lambda }_{j+1}\left(\alpha \right)}{\overset{{\omega }_{j+1}\left(\alpha \right)}{\int }}\cdots {\displaystyle \underset{{\lambda }_{TI}\left(\alpha \right)}{\overset{{\omega }_{TI}\left(\alpha \right)}{\int }}S\left[u\left(x_1,.,{\lambda }_j\left(\alpha \right),.,x_{N_x}\right);\alpha \right]\frac{\partial {\lambda }_j\left(\alpha \right)}{\partial {\alpha }_{j_1}}\text{d}{x}_1\cdots \text{d}{x}_{TI}}}}}\right\}}_{{\alpha }^0}}.\end{array}$ (11)

The direct-effect term can be computed once the nominal values $e^0=\left(u^0,{\alpha }^0\right)$ are available. The notation ${\left\{\right\}}_{{\alpha }^0}$ will be used in this work to indicate that the quantity enclosed within the bracket is to be evaluated at the respective nominal parameter and state functions values.

On the other hand, the quantity ${\left\{\delta R\left[u\left(x\right);\alpha ;v^{\left(1\right)}\left(x\right)\right]\right\}}_{ind}$ in Equation (8) comprises only variations in the state functions and is therefore called the “indirect-effect term,” having the following expression:

${\left\{\delta R\left[u\left(x\right);\alpha ;v^{\left(1\right)}\left(x\right)\right]\right\}}_{ind}\triangleq {\displaystyle \underset{{\lambda }_1\left({\alpha }^0\right)}{\overset{{\omega }_1\left({\alpha }^0\right)}{\int }}\mathrm{...}{\displaystyle \underset{{\lambda }_{TI}\left({\alpha }^0\right)}{\overset{{\omega }_{TI}\left({\alpha }^0\right)}{\int }}{\left\{\frac{\partial S\left(u;\alpha ;x\right)}{\partial u}\right\}}_{{\alpha }^0}{v}^{\left(1\right)}\left(x\right)\text{d}{x}_1\cdots \text{d}{x}_{TI}}}$, (12)

where

$\frac{\partial \left[\right]}{\partial u}{v}^{\left(1\right)}\left(x\right)\triangleq {\displaystyle \underset{i=1}{\overset{TD}{\sum }}\frac{\partial \left[\right]}{\partial u_i\left(x\right)}\delta u_i\left(x\right)}$. (13)

The “indirect-effect” term induces variations in the response through the variations in the state functions, which are, in turn, caused by the parameter variations through the equations underlying the model. Evidently, the indirect-effect term can be quantified only after having determined the variations $v^{\left(1\right)}\left(x\right)$ in terms of the variations $\delta \alpha$. The first-order relationship between the vectors $v^{\left(1\right)}\left(x\right)$ and $\delta \alpha$ is determined by solving the equations obtained by applying the definition of the G-differential to Equations (1) and (2), which yields the following equations:

${\left\{V^{\left(1\right)}\left(u;\alpha \right)v^{\left(1\right)}\left(x\right)\right\}}_{{\alpha }^0}={\left\{q_V^{\left(1\right)}\left(u;\alpha ;\delta \alpha \right)\right\}}_{{\alpha }^0},x\in {\Omega }_x,$ (14)

${\left\{b_V^{\left(1\right)}\left(u;\alpha ;v^{\left(1\right)};\delta \alpha \right)\right\}}_{{\alpha }^0}=0,x\in \partial {\Omega }_x\left({\alpha }^0\right).$ (15)

In Equations (14) and (15), the superscript “(1)” indicates “1^st-Level” and the various quantities which appear in these equations are defined as follows:

$V^{\left(1\right)}\left(u;\alpha \right)\triangleq \left\{\frac{\partial N\left(u;\alpha \right)}{\partial u}\right\}\triangleq \left(\begin{array}{ccc}\frac{\partial N_1}{\partial u_1}& \dots & \frac{\partial N_1}{\partial u_{TD}}\\ ⋮& \ddots & ⋮\\ \frac{\partial N_{TD}}{\partial u_1}& \cdots & \frac{\partial N_{TD}}{\partial u_{TD}}\end{array}\right);$ (16)

$q_V^{\left(1\right)}\left(u;\alpha ;\delta \alpha \right)\triangleq \frac{\partial \left[Q\left(\alpha \right)-N\left(u;\alpha \right)\right]}{\partial \alpha }\delta \alpha \triangleq {\displaystyle \underset{j_1=1}{\overset{TP}{\sum }}{s}_V^{\left(1\right)}\left(j_1;u;\alpha \right)\delta {\alpha }_{j_1}};$ (17)

$s_V^{\left(1\right)}\left(j_1;u;\alpha \right)\triangleq \frac{\partial \left[Q\left(\alpha \right)-N\left(u;\alpha \right)\right]}{\partial {\alpha }_{j_1}};$ (18)

${\left\{b_V^{\left(1\right)}\left(u;\alpha ;v^{\left(1\right)};\delta \alpha \right)\right\}}_{{\alpha }^0}\triangleq {\left\{\frac{\partial B\left(u;\alpha \right)}{\partial u}\right\}}_{{\alpha }^0}{v}^{\left(1\right)}+{\left\{\frac{\partial \left[B\left(u;\alpha \right)-C\left(\alpha \right)\right]}{\partial \alpha }\right\}}_{{\alpha }^0}\delta \alpha$. (19)

The system comprising Equations (14) and (15) is called the “1^st-Level Variational Sensitivity System” (1^st-LVSS). In order to determine the solutions of the 1^st-LVSS that would correspond to every parameter variation $\delta {\alpha }_{j_1}$, $j_1=1,\cdots ,TP$, the 1^st-LVSS would need to be solved TP times, with distinct right-sides for each $\delta {\alpha }_{j_1}$, thus requiring TP large-scale computations. In other words, the actual form of the 1^st-LVSS that would need to be solved in practice is as follows:

${\left\{V^{\left(1\right)}\left(u;\alpha \right)v^{\left(1\right)}\left(j_1;x\right)\right\}}_{{\alpha }^0}={\left\{s_V^{\left(1\right)}\left(j_1;u;\alpha \right)\right\}}_{{\alpha }^0},j_1=1,\cdots ,TP;x\in {\Omega }_x,$ (20)

${\left\{b_V^{\left(1\right)}\left[u;\alpha ;v^{\left(1\right)}\left(j_1;x\right)\right]\right\}}_{{\alpha }^0}=0;j_1=1,\cdots ,TP;x\in \partial {\Omega }_x\left({\alpha }^0\right).$ (21)

Evidently, Equations (14), (15), (20) and (21) indicate that $v^{\left(1\right)}\left(x\right)={\displaystyle \underset{j_1=1}{\overset{TP}{\sum }}{v}^{\left(1\right)}\left(j_1;x\right)\delta {\alpha }_{j_1}}$. In most practical situations, the number of model parameters significantly exceeds the number of functional responses of interest, i.e., $TR\ll TP$, so it would be advantageous to perform just TR (rather than TP) computations. The goal of the “1^st-Order Comprehensive Adjoint Sensitivity Analysis Methodology for Nonlinear Systems (1^st-CASAM-N)” is to compute exactly and efficiently the “indirect effect term” defined in Equation (12) without needing to compute explicitly the vectors $v^{\left(1\right)}\left(j_1;x\right)$, $j_1=1,\cdots ,TP$.

As has been originally shown by Cacuci [1], the need for computing the vectors $v^{\left(1\right)}\left(j_1;x\right)$, $j_1=1,\cdots ,TP$ is eliminated by expressing the indirect-effect term defined in Equation (12) in terms of the solutions of the “1^st-Level Adjoint Sensitivity System” (1^st-LASS), the construction of which requires the introduction of adjoint operators. This is accomplished by introducing a (real) Hilbert space denoted as $H_1\left({\Omega }_x\right)$, endowed with an inner product of two vectors $u^{\left(a\right)}\left(x\right)\in H_1$ and $u^{\left(b\right)}\left(x\right)\in H_1$ denoted as ${\langle u^{\left(a\right)},u^{\left(b\right)}\rangle }_1$ and defined as follows:

${\langle u^{\left(a\right)},u^{\left(b\right)}\rangle }_1\triangleq {\left\{{\displaystyle \underset{{\lambda }_1\left(\alpha \right)}{\overset{{\omega }_1\left(\alpha \right)}{\int }}\cdots {\displaystyle \underset{{\lambda }_{TI}\left(\alpha \right)}{\overset{{\omega }_{TI}\left(\alpha \right)}{\int }}\left[u^{\left(a\right)}\left(x\right)\cdot u^{\left(b\right)}\left(x\right)\right]\text{d}{x}_1\cdots \text{d}{x}_{TI}}}\right\}}_{{\alpha }^0}$, (22)

where the dot indicates the scalar product $u^{\left(a\right)}\left(x\right)\cdot u^{\left(b\right)}\left(x\right)\triangleq {\displaystyle \underset{i=1}{\overset{TD}{\sum }}{u}_i^{\left(a\right)}\left(x\right)u_i^{\left(b\right)}}\left(x\right)$.

Using the inner product defined in Equation (22), construct the inner product of Equation (14) with a vector $a^{\left(1\right)}\left(x\right)$ to obtain the following relation:

${\left\{{\left\{{\langle a^{\left(1\right)},V^{\left(1\right)}\left(u;\alpha \right)v^{\left(1\right)}\rangle }_1\right\}}_{{\alpha }^0}\right\}}_{{\alpha }^0}={\left\{{\langle a^{\left(1\right)},q_V^{\left(1\right)}\left(u;\alpha ;\delta \alpha \right)\rangle }_1\right\}}_{{\alpha }^0},x\in {\Omega }_x.$ (23)

Using the definition of the adjoint operator in $H_1\left({\Omega }_x\right)$, the left-side of Equation (23) is transformed as follows:

$\begin{array}{l}{\left\{{\langle a^{\left(1\right)},V^{\left(1\right)}\left(u;\alpha \right)v^{\left(1\right)}\rangle }_1\right\}}_{{\alpha }^0}\\ ={\left\{{\langle A^{\left(1\right)}\left(u;\alpha \right)a^{\left(1\right)},v^{\left(1\right)}\rangle }_1\right\}}_{{\alpha }^0}+{\left\{{\left[P^{\left(1\right)}\left(u;\alpha ;a^{\left(1\right)};v^{\left(1\right)}\right)\right]}_{\partial {\Omega }_x}\right\}}_{{\alpha }^0},\end{array}$ (24)

where $A^{\left(1\right)}\left(u;\alpha \right)$ is the operator adjoint to $V^{\left(1\right)}\left(u;\alpha \right)$, i.e., $A^{\left(1\right)}\left(u;\alpha \right)\triangleq {\left[V^{\left(1\right)}\left(u;\alpha \right)\right]}^{*}$, and where ${\left[P^{\left(1\right)}\left(u;\alpha ;a^{(1)};v^{\left(1\right)}\right)\right]}_{\partial {\Omega }_x}$ denotes the associated bilinear concomitant evaluated on the space/time domain’s boundary $\partial {\Omega }_x\left({\alpha }^0\right)$. The symbol ${\left[\right]}^{\ast }$ is used in this work to indicate “adjoint” operator. In certain situations, it might be computationally advantageous to include certain boundary components of ${\left[P^{\left(1\right)}\left(u;\alpha ;a^{\left(1\right)};v^{\left(1\right)}\right)\right]}_{\partial {\Omega }_x}$ into the components of $A^{\left(1\right)}\left(u;\alpha \right)$.

The first term on the right-side of Equation (24) is required to represent the indirect-effect term defined in Equation (12) by imposing the following relationship:

${\left\{A^{\left(1\right)}\left(u;\alpha \right)a^{\left(1\right)}\left(x\right)\right\}}_{{\alpha }^0}={\left\{\partial S\left(u;\alpha \right)/\partial u\right\}}_{{\alpha }^0}\triangleq q_A^{\left(1\right)}\left[u\left(x\right);\alpha \right],x\in {\Omega }_x,$ (25)

The domain of $A^{\left(1\right)}\left(u;\alpha \right)$ is determined by selecting appropriate adjoint boundary and/or initial conditions, which will be denoted in operator form as:

${\left\{b_A^{\left(1\right)}\left(u;a^{\left(1\right)};\alpha \right)\right\}}_{{\alpha }^0}=0,x\in \partial {\Omega }_x\left({\alpha }^0\right).$ (26)

The above boundary conditions for $A^{\left(1\right)}\left(u;\alpha \right)$ are usually inhomogeneous, i.e., $b_A^{\left(1\right)}\left(0;0;\alpha \right)\ne 0$, and are obtained by imposing the following requirements: 1) they must be independent of unknown values of $v^{\left(1\right)}\left(x\right)$ and $\delta \alpha$ ; 2) the substitution of the boundary and/or initial conditions represented by Equations (15) and (26) into the expression of ${\left\{{\left[P^{\left(1\right)}\left(u;\alpha ;a^{\left(1\right)};v^{\left(1\right)}\right)\right]}_{\partial {\Omega }_x}\right\}}_{{\alpha }^0}$ must cause all terms containing unknown values of $v^{\left(1\right)}\left(x\right)$ to vanish. Constructing the adjoint initial and/or boundary conditions for $A^{\left(1\right)}\left(u;\alpha \right)$ as described above and implementing them together with the variational boundary and initial conditions represented by Equations (15) into Equation (24) reduces the bilinear concomitant

${\left\{{\left[P^{\left(1\right)}\left(u;\alpha ;a^{\left(1\right)};v^{\left(1\right)}\right)\right]}_{\partial {\Omega }_x}\right\}}_{{\alpha }^0}$ to a quantity denoted as ${\left\{{\left[{\hat{P}}^{\left(1\right)}\left(u;\alpha ;a^{\left(1\right)};\delta \alpha \right)\right]}_{\partial {\Omega }_x}\right\}}_{{\alpha }^0}$, which will contain boundary terms involving only known values of $\delta \alpha$, ${\alpha }^0$, $u^0$, and ${\psi }^{\left(1\right)}$ Since ${\left\{{\left[{\hat{P}}^{\left(1\right)}\left(u;\alpha ;a^{\left(1\right)};\delta \alpha \right)\right]}_{\partial {\Omega }_x}\right\}}_{{\alpha }^0}$ is linear in $\delta \alpha$, it can be expressed in the following form: ${\left[{\hat{P}}^{\left(1\right)}\left(u;\alpha ;a^{\left(1\right)};\delta \alpha \right)\right]}_{\partial {\Omega }_x}={\displaystyle \underset{j_1=1}{\overset{TP}{\sum }}\left[\partial {\hat{P}}^{\left(1\right)}\left(u;\alpha ;a^{\left(1\right)}\right)/\partial {\alpha }_{j_1}\right]\delta {\alpha }_{j_1}}$.

The results obtained in Equations (24) and (25) are now replaced in Equation (12) to obtain the following expression of the indirect-effect term as a function of $a^{\left(1\right)}\left(x\right)$ :

$\begin{array}{l}{\left\{\delta R\left[u\left(x\right);\alpha ;v^{\left(1\right)}\left(x\right)\right]\right\}}_{ind}\\ ={\left\{{\langle a^{\left(1\right)},q_V^{\left(1\right)}\left(u;\alpha ;\delta \alpha \right)\rangle }_1\right\}}_{{\alpha }^0}-{\left\{{\left[{\hat{P}}^{\left(1\right)}\left(u;\alpha ;a^{\left(1\right)};\delta \alpha \right)\right]}_{\partial {\Omega }_x}\right\}}_{{\alpha }^0},\end{array}$ (27)

Replacing in Equation (8) the result obtained in Equation (27) together with the expression for the direct-effect term provided in Equation (9) yields the following expression for the first G-differential of the response $R\left[u\left(x\right);\alpha \right]$ :

$\begin{array}{l}{\left\{\delta R\left[u\left(x\right);\alpha ;v^{\left(1\right)}\left(x\right);\delta \alpha \right]\right\}}_{{\alpha }^0}\\ ={\left\{\delta R\left[u\left(x\right);\alpha ;\delta \alpha \right]\right\}}_{dir}+{\left\{{\langle a^{\left(1\right)},q_V^{\left(1\right)}\left(u;\alpha ;\delta \alpha \right)\rangle }_1\right\}}_{{\alpha }^0}\\ -{\left\{{\left[{\hat{P}}^{\left(1\right)}\left(u;\alpha ;a^{\left(1\right)};\delta \alpha \right)\right]}_{\partial {\Omega }_x}\right\}}_{{\alpha }^0}\\ \triangleq {\displaystyle \underset{j_1=1}{\overset{TP}{\sum }}{\left\{R^{\left(1\right)}\left[j_1;u\left(x\right);a^{\left(1\right)}\left(x\right);\alpha \right]\right\}}_{{\alpha }^0}\delta {\alpha }_{j_1}},\end{array}$ (28)

where, for each $j_1=1,\cdots ,TP$, the quantity $R^{\left(1\right)}\left[j_1;u\left(x\right);a^{\left(1\right)}\left(x\right);\alpha \right]$ denotes the 1^st-order sensitivities of the response $R\left[u\left(x\right);\alpha \right]$ with respect to the model parameters ${\alpha }_{j_1}$ and has the following expression:

$\begin{array}{l}{R}^{\left(1\right)}\left[j_1;u\left(x\right);a^{\left(1\right)}\left(x\right);\alpha \right]\\ ={\displaystyle \underset{{\lambda }_1\left(\alpha \right)}{\overset{{\omega }_1\left(\alpha \right)}{\int }}\cdots {\displaystyle \underset{{\lambda }_{TI}\left(\alpha \right)}{\overset{{\omega }_{TI}\left(\alpha \right)}{\int }}{a}^{\left(1\right)}\left(x\right)\cdot \frac{\partial \left[Q\left(\alpha \right)-N\left(u;\alpha \right)\right]}{\partial {\alpha }_{j_1}}\text{d}{x}_1\cdots \text{d}{x}_{TI}}}\\ -\frac{\partial {\hat{P}}^{\left(1\right)}\left(u;\alpha ;a^{\left(1\right)}\right)}{\partial {\alpha }_{j_1}}+{\displaystyle \underset{{\lambda }_1\left(\alpha \right)}{\overset{{\omega }_1\left(\alpha \right)}{\int }}\cdots {\displaystyle \underset{{\lambda }_{TI}\left(\alpha \right)}{\overset{{\omega }_{TI}\left(\alpha \right)}{\int }}\frac{\partial S\left(u;\alpha ;\alpha \right)}{\partial {\alpha }_{j_1}}\text{d}{x}_1\cdots \text{d}{x}_{TI}}}\end{array}$

$\begin{array}{l}+{\displaystyle \underset{j=1}{\overset{TI}{\sum }}{\left\{{\displaystyle \underset{{\lambda }_1\left(\alpha \right)}{\overset{{\omega }_1\left(\alpha \right)}{\int }}\cdots {\displaystyle \underset{{\lambda }_{j-1}\left(\alpha \right)}{\overset{{\omega }_{j-1}\left(\alpha \right)}{\int }}{\displaystyle \underset{{\lambda }_{j+1}\left(\alpha \right)}{\overset{{\omega }_{j+1}\left(\alpha \right)}{\int }}\cdots {\displaystyle \underset{{\lambda }_{TI}\left(\alpha \right)}{\overset{{\omega }_{TI}\left(\alpha \right)}{\int }}S\left[u\left(x_1,.,{\omega }_j\left(\alpha \right),.,x_{N_x}\right);\alpha \right]\frac{\partial {\omega }_j\left(\alpha \right)}{\partial {\alpha }_{j_1}}\text{d}{x}_1\cdots \text{d}{x}_{TI}}}}}\right\}}_{{\alpha }^0}}\\ -{\displaystyle \underset{j=1}{\overset{TI}{\sum }}{\left\{{\displaystyle \underset{{\lambda }_1\left(\alpha \right)}{\overset{{\omega }_1\left(\alpha \right)}{\int }}\cdots {\displaystyle \underset{{\lambda }_{j-1}\left(\alpha \right)}{\overset{{\omega }_{j-1}\left(\alpha \right)}{\int }}{\displaystyle \underset{{\lambda }_{j+1}\left(\alpha \right)}{\overset{{\omega }_{j+1}\left(\alpha \right)}{\int }}\cdots {\displaystyle \underset{{\lambda }_{TI}\left(\alpha \right)}{\overset{{\omega }_{TI}\left(\alpha \right)}{\int }}S\left[u\left(x_1,.,{\lambda }_j\left(\alpha \right),.,x_{N_x}\right);\alpha \right]\frac{\partial {\lambda }_j\left(\alpha \right)}{\partial {\alpha }_{j_1}}\text{d}{x}_1\cdots \text{d}{x}_{TI}}}}}\right\}}_{{\alpha }^0}}.\end{array}$ (29)

As indicated by Equation (29), each of the 1^st-order sensitivities $R^{\left(1\right)}\left[j_1;u\left(x\right);a^{\left(1\right)}\left(x\right);\alpha \right]$ of the response $R\left[u\left(x\right);\alpha \right]$ with respect to the model parameters ${\alpha }_{j_1}$ (including boundary and initial conditions) can be computed inexpensively after having obtained the function $a^{\left(1\right)}\left(x\right)\in H_1$, using quadrature formulas to evaluate the various inner products involving $a^{\left(1\right)}\left(x\right)\in H_1$. The function $a^{\left(1\right)}\left(x\right)\in H_1$ is obtained by solving numerically Equations (25) and (26), which is the only large-scale computation needed for obtaining all of the first-order sensitivities. Equations (26) and (25) are called the 1^st-Level Adjoint Sensitivity System (1^st-LASS), and its solution, $a^{\left(1\right)}\left(x\right)\in H_1\left({\Omega }_x\right)$, is called the 1^st-level adjoint function. It is very important to note that the 1^st-LASS is independent of parameter variation $\delta {\alpha }_{j_1}$, $j_1=1,\cdots ,TP$, and therefore needs to be solved only once, regardless of the number of model parameters under consideration. Furthermore, since Equation (25) is linear in $a^{\left(1\right)}\left(x\right){\psi }_{1,i_1}^{\left(2\right)}\left(x\right)$, solving it requires less computational effort than solving the original Equation (1), which is nonlinear in $u\left(x\right)$.

2.2. The Second-Order Comprehensive Adjoint Sensitivity Analysis Methodology for Nonlinear Systems (2^nd-CASAM-N)

The 2^nd-CASAM-N relies on the same fundamental concepts as introduced in [4], but in addition to the capabilities described in [4], the 2^nd-CASAM-N also enables the computation of response sensitivities with respect to imprecisely known domain boundaries, thus including all possible types of uncertain parameters. Fundamentally, the 2^nd-order sensitivities are defined as the “1^st-order sensitivities of the 1^st-order sensitivities.” This definition stems from the inductive definition of the 2^nd-order total G-differential of correspondingly differentiable function, which is also defined inductively as “the total 1^st-order differential of the 1^st-order total differential” of a function. The 1^st-order sensitivities $R^{\left(1\right)}\left[j_1;u\left(x\right);a^{\left(1\right)}\left(x\right);\alpha \right]$ are assumed to satisfy the conditions stated in Equations (6) and (7), for each $j_1=1,\cdots ,TP$, which ensures the existence of the 2^nd-order sensitivities. The G-variation ${\left\{\delta R^{\left(1\right)}\left[j_1;u\left(x\right);a^{\left(1\right)}\left(x\right);\alpha ;v^{\left(1\right)}\left(x\right);\delta a^{\left(1\right)}\left(x\right);\delta \alpha \right]\right\}}_{{\alpha }^0}$ of a 1^st-order sensitivity $R^{\left(1\right)}\left[j_1;u\left(x\right);a^{\left(1\right)}\left(x\right);\alpha \right]$ has the following expression:

$\begin{array}{l}{\left\{\delta R^{\left(1\right)}\left[j_1;u\left(x\right);a^{\left(1\right)}\left(x\right);\alpha ;v^{\left(1\right)}\left(x\right);\delta a^{\left(1\right)}\left(x\right);\delta \alpha \right]\right\}}_{{\alpha }^0}\\ ={\left\{\delta R^{\left(1\right)}\left[j_1;u\left(x\right);a^{\left(1\right)}\left(x\right);\alpha ;\delta \alpha \right]\right\}}_{dir}\\ +{\left\{\delta R^{\left(1\right)}\left[j_1;u\left(x\right);a^{\left(1\right)}\left(x\right);\alpha ;v^{\left(1\right)}\left(x\right);\delta a^{\left(1\right)}\left(x\right)\right]\right\}}_{ind}.\end{array}$ (30)

In Equation (30), the quantity ${\left\{\delta R^{\left(1\right)}\left[j_1;u\left(x\right);a^{\left(1\right)}\left(x\right);\alpha ;\delta \alpha \right]\right\}}_{dir}$ denotes the direct-effect term, which comprises all dependencies on the vector $\delta \alpha$ of parameter variations, and is defined as follows:

${\left\{\delta R^{\left(1\right)}\left[j_1;u\left(x\right);a^{\left(1\right)}\left(x\right);\alpha ;\delta \alpha \right]\right\}}_{dir}\triangleq {\left\{\frac{\partial R^{\left(1\right)}\left[j_1;u\left(x\right);a^{\left(1\right)}\left(x\right);\alpha \right]}{\partial \alpha }\delta \alpha \right\}}_{{\alpha }^0}$. (31)

Also in Equation (30), the indirect-effect term ${\left\{\delta R^{\left(1\right)}\left[j_1;u\left(x\right);a^{\left(1\right)}\left(x\right);\alpha ;v^{\left(1\right)}\left(x\right);\delta a^{\left(1\right)}\left(x\right)\right]\right\}}_{ind}$ comprises all dependencies on the vectors $v^{\left(1\right)}\left(x\right)$ and $\delta a^{\left(1\right)}\left(x\right)$ of variations in the state functions $u\left(x\right)$ and $a^{\left(1\right)}\left(x\right)$, respectively, and is defined as follows:

$\begin{array}{l}{\left\{\delta R^{\left(1\right)}\left[j_1;u\left(x\right);a^{\left(1\right)}\left(x\right);\alpha ;v^{\left(1\right)}\left(x\right);\delta a^{\left(1\right)}\left(x\right)\right]\right\}}_{ind}\\ \triangleq {\left\{\partial R^{\left(1\right)}\left[j_1;\cdots ;\alpha \right]/\partial u\right\}}_{{\alpha }^0}{v}^{\left(1\right)}\left(x\right)+{\left\{\partial R^{\left(1\right)}\left[j_1;\cdots ;\alpha \right]/\partial a^{\left(1\right)}\right\}}_{{\alpha }^0}\delta a^{\left(1\right)}\left(x\right).\end{array}$ (32)

The functions $v^{\left(1\right)}\left(x\right)$ and $\delta a^{\left(1\right)}\left(x\right)$ are obtained by solving the following 2^nd-Level Variational Sensitivity System (2^nd-LVSS):

$\begin{array}{l}{\left\{V{M}^{\left(2\right)}\left[2\times 2;U^{\left(2\right)}\left(2;x\right);\alpha \right]V^{\left(2\right)}\left(2;x\right)\right\}}_{{\alpha }^0}\\ ={\left\{Q_V^{\left(2\right)}\left[2;U^{\left(2\right)}\left(2;x\right);\alpha ;\delta \alpha \right]\right\}}_{{\alpha }^0},x\in {\Omega }_x,\end{array}$ (33)

${\left\{B_V^{\left(2\right)}\left[2;U^{\left(2\right)}\left(2;x\right);V^{\left(2\right)}\left(2;x\right);\alpha ;\delta \alpha \right]\right\}}_{{\alpha }^0}=0\left[2\right],0\left[2\right]\triangleq {\left[0,0\right]}^{†},x\in \partial {\Omega }_x\left({\alpha }^0\right)$.(34)

The argument “2” which appears in the list of arguments of the vector $U^{\left(2\right)}\left(2;x\right)$ and the “variational vector” $V^{\left(2\right)}\left(2;x\right)$ in Equation (33) indicates that each of these vectors is a 2-block column vector (each block comprising a column-vector of dimension TD), defined as follows:

$U^{\left(2\right)}\left(2;x\right)\triangleq \left(\begin{array}{c}{u}^{\left(1\right)}\left(x\right)\\ a^{\left(1\right)}\left(x\right)\end{array}\right);V^{\left(2\right)}\left(2;x\right)\triangleq \delta U^{\left(2\right)}\left(2;x\right)\triangleq \left(\begin{array}{c}{v}^{\left(2\right)}\left(1;x\right)\\ v^{\left(2\right)}\left(2;x\right)\end{array}\right)\triangleq \left(\begin{array}{c}{v}^{\left(1\right)}\left(x\right)\\ \delta a^{\left(1\right)}\left(x\right)\end{array}\right).$ (35)

To distinguish block-vectors from block matrices, two capital bold letter have been used (and will henceforth be used) to denote block matrices, as in the case of “the second-level variational matrix” $V{M}^{\left(2\right)}\left[2\times 2;u^{\left(2\right)}\left(x\right);\alpha \right]$. The “2^nd-level” is indicated by the superscript “(2)”. Subsequently in this work, levels higher than second will also be indicated by a corresponding superscript attached to the appropriate block-vectors and/or block-matrices. The argument “ $2\times 2$ ”, which appears in the list of arguments of $V{M}^{\left(2\right)}\left[2\times 2;u^{\left(2\right)}\left(x\right);\alpha \right]$, indicates that this matrix is a $2\times 2$ -dimensional block-matrix comprising four matrices, each of dimensions $TD\times TD$, having the following structure:

$V{M}^{\left(2\right)}\left[2\times 2;U^{\left(2\right)}\left(2;x\right);\alpha \right]\triangleq \left(\begin{array}{cc}{V}^{\left(1\right)}& 0\\ V_{21}^{\left(2\right)}& V_{22}^{\left(2\right)}\end{array}\right).$ (36)

The other quantities which appear in Equations (33) and (34) are 2-block vectors having the same structure as $V^{\left(2\right)}\left(2;x\right)$, and are defined as follows:

$Q_V^{\left(2\right)}\left[2;U^{\left(2\right)}\left(2;x\right);\alpha ;\delta \alpha \right]\triangleq \left(\begin{array}{c}{q}_V^{\left(2\right)}\left(1;U^{\left(2\right)}\left(2;x\right);\alpha ;\delta \alpha \right)\\ q_V^{\left(2\right)}\left(2;U^{\left(2\right)}\left(2;x\right);\alpha ;\delta \alpha \right)\end{array}\right)\triangleq \left(\begin{array}{c}{q}_V^{\left(1\right)}\left(u;\alpha ;\delta \alpha \right)\\ q_2^{\left(2\right)}\left(u;a^{\left(1\right)};\alpha ;\delta \alpha \right)\end{array}\right);$ (37)

$q_2^{\left(2\right)}\left(u;\alpha ;a^{\left(1\right)};\delta \alpha \right)\triangleq \frac{\partial q_A^{\left(1\right)}\left[u\left(x\right);\alpha \right]}{\partial \alpha }\delta \alpha -\frac{\partial \left[A^{\left(1\right)}\left(u;\alpha \right)a^{\left(1\right)}\left(x\right)\right]}{\partial \alpha }\delta \alpha$ ; (38)

$\begin{array}{l}{B}_V^{\left(2\right)}\left[2;U^{\left(2\right)}\left(2;x\right);V^{\left(2\right)}\left(2;x\right);\alpha ;\delta \alpha \right]\triangleq \left(\begin{array}{c}{b}_V^{\left(2\right)}\left[1;U^{\left(2\right)}\left(2;x\right);V^{\left(2\right)}\left(2;x\right);\alpha ;\delta \alpha \right]\\ b_V^{\left(2\right)}\left[2;U^{\left(2\right)}\left(2;x\right);V^{\left(2\right)}\left(2;x\right);\alpha ;\delta \alpha \right]\end{array}\right)\\ \triangleq \left(\begin{array}{c}{b}_V^{\left(1\right)}\left(u^{\left(1\right)};\alpha ;\delta u^{\left(1\right)};\delta \alpha \right)\\ \delta b_A^{\left(1\right)}\left[U^{\left(2\right)}\left(2;x\right);V^{\left(2\right)}\left(2;x\right);\alpha ;\delta \alpha \right]\end{array}\right).\end{array}$ (39)

$V_{21}^{\left(2\right)}\left(u;a^{\left(1\right)};\alpha \right)\triangleq \frac{\partial \left[A^{\left(1\right)}\left(u;\alpha \right)a^{\left(1\right)}\right]}{\partial u}-\frac{q_A^{\left(1\right)}\left[u\left(x\right);\alpha \right]}{\partial u};$ (40)

$V_{22}^{\left(2\right)}\left(u;\alpha \right)\triangleq A^{\left(1\right)}\left(u;\alpha \right)$ ; (41)

$\begin{array}{c}\delta b_A^{\left(1\right)}\left(u;a^{\left(1\right)};\alpha \right)\triangleq \frac{\partial b_A^{\left(1\right)}\left(u;a^{\left(1\right)};\alpha \right)}{\partial u}{v}^{\left(1\right)}\left(x\right)+\frac{\partial b_A^{\left(1\right)}\left(u;a^{\left(1\right)};\alpha \right)}{\partial a^{\left(1\right)}}\delta a^{\left(1\right)}\left(x\right)\\ +\frac{\partial b_A^{\left(1\right)}\left(u;a^{\left(1\right)};\alpha \right)}{\partial \alpha }\delta \alpha .\end{array}$ (42)

The structure of the second component of the source-term $Q_V^{\left(2\right)}\left[2;U^{\left(2\right)}\left(2;x\right);\alpha ;\delta \alpha \right]$ defined in Equation (37) is as follows:

$q_V^{\left(2\right)}\left[2;U^{\left(2\right)}\left(2;x\right);\alpha ;\delta \alpha \right]\triangleq {\displaystyle \underset{j_2=1}{\overset{TP}{\sum }}{s}_V^{\left(2\right)}\left[2;j_2;U^{\left(2\right)}\left(2;x\right);\alpha \right]\delta {\alpha }_{j_2},}$ (43)

where

$s_V^{\left(2\right)}\left[2;j_2;U^{\left(2\right)}\left(2;x\right);\alpha \right]\triangleq \frac{\partial q_A^{\left(1\right)}\left[u\left(x\right);\alpha \right]}{\partial {\alpha }_{j_2}}-\frac{\partial \left[A^{\left(1\right)}\left(u;\alpha \right)a^{\left(1\right)}\left(x\right)\right]}{\partial {\alpha }_{j_2}}$. (44)

Taking into account the expressions in Equations (43) and (44) while recalling the expressions in Equations (17) and (18) indicates the actual form of 2^nd-LVSS to be solved (if one would wish to solve it) would be as follows:

$\begin{array}{l}{\left\{V{M}^{\left(2\right)}\left[2\times 2;U^{\left(2\right)}\left(2;x\right);\alpha \right]V^{\left(2\right)}\left(2;j_2;j_1;x\right)\right\}}_{{\alpha }^0}\\ ={\left\{Q_V^{\left(2\right)}\left[2;U^{\left(2\right)}\left(2;j_2;j_1;x\right);\alpha ;\delta \alpha \right]\right\}}_{{\alpha }^0},j_1=1,\cdots ,TP;j_2=1,\cdots ,TP;x\in {\Omega }_x,\end{array}$ (45)

${\left\{B_V^{\left(2\right)}\left[2;j_2;j_1;U^{\left(2\right)}\left(2;x\right);V^{\left(2\right)}\left(2;x\right);\alpha ;\delta \alpha \right]\right\}}_{{\alpha }^0}=0\left[2\right],x\in \partial {\Omega }_x\left({\alpha }^0\right)$. (46)

Thus, there would be $T{P}^2$ “variational vectors” $V^{\left(2\right)}\left(2;j_2;j_1;x\right)$ to be computed. The need to avoid such impractical, extremely expensive, computations provides the fundamental motivation underlying the development of the adjoint sensitivity analysis methodologies for computing sensitivities (of all orders) of model responses with respect to the model’s parameters. Thus, since “variational sensitivity systems” will never need to be actually solved if the 5^th-CASAM-N methodology developed and presented in this work is utilized, the dependence on the indices $j_1$, $j_2$ of the variation sensitivity systems will be suppressed in this work, in order to simplify the mathematical notation. On the other hand, since the solutions of the adjoint sensitivity systems of various levels will actually be computed in practice, the dependence on the indices $j_1,j_1=1,\cdots ,TP$ will be displayed explicitly.

The need for solving the 2^nd-LVSS is circumvented by deriving an alternative expression for the indirect-effect term defined in Equation (32), in which the function $V^{\left(2\right)}\left(2;x\right)$ is replaced by a 2^nd-level adjoint function which is independent of variations in the model parameter and state functions, and is the solution of a 2^nd-Level Adjoint Sensitivity System (2^nd-LASS) which is constructed by using the 2^nd-LVSS as starting point and following the same principles as outlined in Section 2.1. The 2^nd-LASS is constructed in a Hilbert space, denoted as $H_2\left({\Omega }_x\right)$, which comprises as elements block-vectors of the same form as $V^{\left(2\right)}\left(2;x\right)$. The inner product of two vectors ${\Psi }^{\left(2\right)}\left(2;x\right)\triangleq {\left[{\psi }^{\left(2\right)}\left(1;x\right),{\psi }^{\left(2\right)}\left(2;x\right)\right]}^{†}\in H_2\left({\Omega }_x\right)$ and ${\Phi }^{\left(2\right)}\left(2;x\right)\triangleq {\left[{\phi }^{\left(2\right)}\left(1;x\right),{\phi }^{\left(2\right)}\left(2;x\right)\right]}^{†}\in H_2\left({\Omega }_x\right)$ in the Hilbert space $H_2\left({\Omega }_x\right)$ will be denoted as ${\langle {\Psi }^{\left(2\right)}\left(2;x\right),{\Phi }^{\left(2\right)}\left(2;x\right)\rangle }_2$ and defined as follows:

${\langle {\Psi }^{\left(2\right)}\left(2;x\right),{\Phi }^{\left(2\right)}\left(2;x\right)\rangle }_2\triangleq {\displaystyle \underset{i=1}{\overset{2}{\sum }}{\langle {\psi }^{\left(2\right)}\left(i;x\right),{\phi }^{\left(2\right)}\left(i;x\right)\rangle }_1}$. (47)

Following the same principles as outlined in Section 2.1, the inner product defined in Equation (47) is used to construct the following 2^nd-Level Adjoint Sensitivity System (2^nd-LASS) for the 2^nd-level adjoint function $A^{\left(2\right)}\left(2;j_1;x\right)\triangleq {\left[a^{\left(2\right)}\left(1;j_1;x\right),a^{\left(2\right)}\left(2;j_1;x\right)\right]}^{†}\in H_2\left({\Omega }_x\right)$, for each $j_1=1,\cdots ,TP$ :

$\begin{array}{l}{\left\{A{M}^{\left(2\right)}\left[2\times 2;U^{\left(2\right)}\left(2;x\right);\alpha \right]A^{\left(2\right)}\left(2;j_1;x\right)\right\}}_{{\alpha }^0}\\ ={\left\{Q_A^{\left(2\right)}\left[2;j_1;U^{\left(2\right)}\left(2;x\right);\alpha \right]\right\}}_{{\alpha }^0},j_1=1,\cdots ,TP;x\in {\Omega }_x,\end{array}$ (48)

subject to boundary conditions represented as follows:

${\left\{B_A^{\left(2\right)}\left[2;U^{\left(2\right)}\left(2;x\right);A^{\left(2\right)}\left(2;j_1;x\right);\alpha \right]\right\}}_{{\alpha }^0}=0\left[2\right];j_1=1,\cdots ,TP;x\in \partial {\Omega }_x\left({\alpha }^0\right)$.(49)

where:

$\begin{array}{l}{Q}_A^{\left(2\right)}\left[2;j_1;U^{\left(2\right)}\left(2;x\right);\alpha \right]\triangleq \left(\begin{array}{c}{q}_A^{\left(2\right)}\left(1;j_1;U^{\left(2\right)};\alpha \right)\\ q_A^{\left(2\right)}\left(2;j_1;U^{\left(2\right)};\alpha \right)\end{array}\right)\\ \triangleq \left(\begin{array}{c}\partial R^{\left(1\right)}\left[j_1;u\left(x\right);a^{\left(1\right)}\left(x\right);\alpha ;v^{\left(1\right)}\left(x\right)\right]/\partial u\\ \partial R^{\left(1\right)}\left[j_1;u\left(x\right);a^{\left(1\right)}\left(x\right);\alpha ;v^{\left(1\right)}\left(x\right)\right]/\partial a^{\left(1\right)}\end{array}\right),j_1=1,\cdots ,TP.\end{array}$ (50)

$\begin{array}{l}A{M}^{\left(2\right)}\left[2\times 2;U^{\left(2\right)}\left(2;x\right);\alpha \right]\\ \triangleq {\left[V{M}^{\left(2\right)}\left(2\times 2;U^{\left(2\right)}\left(2;x\right);\alpha \right)\right]}^{*}=\left(\begin{array}{cc}{\left[V^{\left(1\right)*}\right]}^{†}& {\left[V_{21}^{\left(2\right)*}\right]}^{†}\\ 0& {\left[V_{22}^{\left(2\right)*}\right]}^{†}\end{array}\right).\end{array}$ (51)

The matrix $A{M}^{\left(2\right)}\left[2\times 2;u^{\left(2\right)}\left(x\right);\alpha \right]$ comprises $\left(2\times 2\right)$ block-matrices, each of dimensions $T{D}^2$, thus comprising a total of $\left(2\times 2\right)T{D}^2$ components (or elements) and is obtained from the following relation:

$\begin{array}{l}{\left\{{\langle A^{\left(2\right)}\left(2;x\right),V{M}^{\left(2\right)}{V}^{\left(2\right)}\left(2;x\right)\rangle }_2\right\}}_{{\alpha }^0}\\ ={\left\{{\left[P^{\left(2\right)}\left(U^{\left(2\right)};A^{\left(2\right)};V^{\left(2\right)};\alpha \right)\right]}_{\partial {\Omega }_x}\right\}}_{{\alpha }^0}\\ +{\left\{{\langle V^{\left(2\right)}\left(2;x\right),A{M}^{\left(2\right)}\left[2\times 2;U^{\left(2\right)}\left(2;x\right);\alpha \right]A^{\left(2\right)}\left(2;x\right)\rangle }_2\right\}}_{{\alpha }^0},\end{array}$ (52)

where the quantity ${\left\{{\left[P^{\left(2\right)}\left(U^{\left(2\right)};A^{\left(2\right)};V^{\left(2\right)};\alpha \right)\right]}_{\partial {\Omega }_x}\right\}}_{{\alpha }^0}$ denotes the corresponding bilinear concomitant on the domain’s boundary, evaluated at the nominal values for the parameters and respective state functions. The 2^nd-level adjoint boundary/initial conditions represented by Equation (49) are determined by requiring that: 1) they must be independent of unknown values of $V^{\left(2\right)}\left(2;x\right)$ ; 2) the substitution of the boundary and/or initial conditions represented by Equations (49)and (34) into the expression of ${\left\{{\left[P^{\left(2\right)}\left(U^{\left(2\right)};A^{\left(2\right)};V^{\left(2\right)};\alpha \right)\right]}_{\partial {\Omega }_x}\right\}}_{{\alpha }^0}$ must cause all terms containing unknown values of $V^{\left(2\right)}\left(2;x\right)$ to vanish.

Using the 2^nd-LASS to obtain the alternative expression for the indirect-effect term in terms of $A^{\left(2\right)}\left(2;j_1;x\right)\triangleq {\left[a^{\left(2\right)}\left(1;j_1;x\right),a^{\left(2\right)}\left(2;j_1;x\right)\right]}^{†}$ and the expression for the direct-effect term provided in Equation (31) yields the following expression for the total differential defined by Equation (30):

$\begin{array}{l}{\left\{\delta R^{\left(1\right)}\left[j_1;U^{\left(2\right)}\left(2;x\right);A^{(2)}\left(2;j_1;x\right);\alpha ;\delta \alpha \right]\right\}}_{{\alpha }^0}\\ ={\left\{\frac{\partial R^{\left(1\right)}\left[j_1;u\left(x\right);a^{\left(1\right)}\left(x\right);\alpha \right]}{\partial \alpha }\delta \alpha \right\}}_{{\alpha }^0}\\ +{\left\{{\langle A^{(2)}\left(2;j_1;x\right),Q_V^{\left(2\right)}\left[2;U^{\left(2\right)}\left(2;x\right);\alpha ;\delta \alpha \right]\rangle }_2\right\}}_{{\alpha }^0}\\ -{\left\{{\left[{\hat{P}}^{\left(2\right)}\left(U^{\left(2\right)};A^{\left(2\right)};\alpha ;\delta \alpha \right)\right]}_{\partial {\Omega }_x}\right\}}_{{\alpha }^0}\\ ={\displaystyle \underset{j_2=1}{\overset{TP}{\sum }}{\left\{R^{\left(2\right)}\left[j_2;j_1;U^{\left(2\right)}\left(2;x\right);A^{\left(2\right)}\left(2;j_1;x\right);\alpha \right]\right\}}_{{\alpha }^0}}\delta {\alpha }_{j_2},j_1=1,\cdots ,TP.\end{array}$ (53)

where ${\left\{{\left[{\hat{P}}^{\left(2\right)}\left(U^{\left(2\right)};A^{(2)};\alpha ;\delta \alpha \right)\right]}_{\partial {\Omega }_x}\right\}}_{{\alpha }^0}$ denotes residual boundary terms which may not have vanished after having used the boundary and/or initial conditions represented by Equations (34) and (49). The detailed operations leading to the expression given in Equation (53) are provided in [27]. The quantity $R^{\left(2\right)}\left[j_2;j_1;U^{\left(2\right)}\left(2;x\right);A^{\left(2\right)}\left(2;j_1;x\right);\alpha \right]$ denotes the second-order sensitivity of the generic scalar-valued response $R\left[u\left(x\right);\alpha \right]$ with respect to the parameters ${\alpha }_{j_1}$ and ${\alpha }_{j_2}$ computed at the nominal values of the parameters and respective state functions, and has the following expression:

$\begin{array}{l}\text{For}{j}_1,j_2=1,\cdots ,TP:R^{\left(2\right)}\left[j_2;j_1;U^{\left(2\right)}\left(2;x\right);A^{\left(2\right)}\left(2;j_1;x\right);\alpha \right]\\ \triangleq \frac{\partial R^{\left(1\right)}\left[j_1;u\left(x\right);a^{\left(1\right)}\left(x\right);\alpha \right]}{\partial {\alpha }_{j_2}}-\frac{{\left\{\partial {\hat{P}}^{\left(2\right)}\left[U^{\left(2\right)}\left(2;x\right);A^{\left(2\right)}\left(2;j_1;x\right);\alpha \right]\right\}}_{\partial {\Omega }_x}}{\partial {\alpha }_{j_2}}\\ +{\displaystyle \underset{i=1}{\overset{2}{\sum }}{\langle a^{\left(2\right)}\left(i;j_1;x\right),s_V^{\left(2\right)}\left[i;j_2;U^{\left(2\right)}\left(2;x\right);\alpha \right]\rangle }_1}\triangleq \frac{{\partial }^{2}R\left[u\left(x\right);\alpha \right]}{\partial {\alpha }_{j_2}\partial {\alpha }_{j_1}}.\end{array}$ (54)

If the 2^nd-LASS is solved TP-times, the 2^nd-order mixed sensitivities $R^{\left(2\right)}\left[j_2;j_1;U^{\left(2\right)}\left(2;x\right);A^{\left(2\right)}\left(2;j_1;x\right);\alpha \right]\equiv {\partial }^{2}R/\partial {\alpha }_{j_2}\partial {\alpha }_{j_1}$ will be computed twice, in two different ways, in terms of two distinct 2^nd-level adjoint functions. Consequently, the symmetry property ${\partial }^{2}R\left[u\left(x\right);\alpha \right]/\partial {\alpha }_{j_2}\partial {\alpha }_{j_1}={\partial }^{2}R\left[u\left(x\right);\alpha \right]/\partial {\alpha }_{j_1}\partial {\alpha }_{j_2}$ enjoyed by the second-order sensitivities provides an intrinsic (numerical) verification that the components of the 2^nd-level adjoint function $A^{\left(2\right)}\left(2;j_1;x\right)$ and the 1^st-level adjoint function $a^{\left(1\right)}\left(x\right)$ are computed accurately.

2.3. The Third-Order Comprehensive Adjoint Sensitivity Analysis Methodology for Nonlinear Systems (3^rd-CASAM-N)

The 3^rd-order sensitivities will be computed by considering them to be the “sensitivities of a 2^nd-order sensitivity.” Thus, each of the 2^nd-order sensitivities $R^{\left(2\right)}\left[j_2;j_1;U^{\left(2\right)}\left(2;x\right);A^{\left(2\right)}\left(2;j_1;x\right);\alpha \right]\equiv {\partial }^{2}R/\partial {\alpha }_{j_2}\partial {\alpha }_{j_1}$ will be considered to be a “model response” which is assumed to satisfy the conditions stated in Equations (6) and (7) for each $j_1,j_2=1,\cdots ,TP$, so that the 1^st-order total G-differential of $R^{\left(2\right)}\left[j_2;j_1;U^{\left(2\right)}\left(2;x\right);A^{\left(2\right)}\left(2;j_1;x\right);\alpha \right]$ will exist and will be linear in the variations $V^{\left(2\right)}\left(2;x\right)$ and $\delta A^{\left(2\right)}\left(2;j_1;x\right)$ in a neighborhood around the nominal values of the parameters and the respective state functions. By definition, the 1^st-order total G-differential of $R^{\left(2\right)}\left[j_2;j_1;U^{\left(2\right)}\left(2;x\right);A^{\left(2\right)}\left(2;j_1;x\right);\alpha \right]$, which will be denoted as ${\left\{\delta R^{\left(2\right)}\left[j_2;j_1;U^{\left(2\right)}\left(2;x\right);A^{\left(2\right)}\left(2;j_1;x\right);\alpha ;V^{\left(2\right)}\left(2;x\right);\delta A^{\left(2\right)}\left(2;j_1;x\right);\delta \alpha \right]\right\}}_{{\alpha }^0}$, is given by the following expression:

$\begin{array}{l}{\left\{\delta R^{\left(2\right)}\left[j_2;j_1;U^{\left(2\right)}\left(2;x\right);A^{\left(2\right)}\left(2;j_1;x\right);\alpha ;V^{\left(2\right)}\left(2;x\right);\delta A^{\left(2\right)}\left(2;j_1;x\right);\delta \alpha \right]\right\}}_{{\alpha }^0}\\ \triangleq {\left\{\frac{\partial R^{\left(2\right)}\left[j_2;j_1;U^{\left(2\right)}\left(2;x\right);A^{\left(2\right)}\left(2;j_1;x\right);\alpha \right]}{\partial \alpha }\delta \alpha \right\}}_{{\alpha }^0}\\ +{\left\{\delta R^{\left(2\right)}\left[j_2;j_1;U^{\left(2\right)}\left(2;x\right);A^{\left(2\right)}\left(2;j_1;x\right);\alpha ;V^{\left(2\right)}\left(2;x\right);\delta A^{\left(2\right)}\left(2;j_1;x\right)\right]\right\}}_{ind},\end{array}$ (55)

where:

$\begin{array}{l}{\left\{\delta R^{\left(2\right)}\left[j_2;j_1;U^{\left(2\right)}\left(2;x\right);A^{\left(2\right)}\left(2;j_1;x\right);\alpha ;V^{\left(2\right)}\left(2;x\right);\delta A^{\left(2\right)}\left(2;j_1;x\right)\right]\right\}}_{ind}\\ \triangleq {\left\{\frac{\partial R^{\left(2\right)}\left[j_2;j_1;U^{\left(2\right)};A^{\left(2\right)};\alpha \right]}{\partial U^{\left(2\right)}\left(2;x\right)}\right\}}_{{\alpha }^0}{V}^{\left(2\right)}\left(2;x\right)\\ \text{}\text{}\text{}+{\left\{\frac{\partial R^{\left(2\right)}\left[j_2;j_1;U^{\left(2\right)};A^{\left(2\right)};\alpha \right]}{\partial A^{\left(2\right)}\left(2;j_1;x\right)}\right\}}_{{\alpha }^0}\delta A^{\left(2\right)}\left(2;j_1;x\right).\end{array}$ (56)

The indirect-effect term ${\left\{\delta R^{\left(2\right)}\left[j_2;j_1;U^{\left(2\right)}\left(2;x\right);A^{\left(2\right)}\left(2;j_1;x\right);\alpha ;\delta \alpha \right]\right\}}_{dir}$ can be computed after having determined the vectors $V^{\left(2\right)}\left(2;x\right)$ and $\delta A^{\left(2\right)}\left(2;j_1;x\right)$, which are the solutions of the following 3^rd-Level Variational Sensitivity System (3^rd-LVSS):

$\begin{array}{l}{\left\{V{M}^{\left(3\right)}\left[4\times 4;U^{\left(3\right)}\left(4;x\right);\alpha \right]V^{\left(3\right)}\left(4;x\right)\right\}}_{{\alpha }^0}\\ ={\left\{Q_V^{\left(3\right)}\left[4;U^{\left(3\right)}\left(4;x\right);\alpha ;\delta \alpha \right]\right\}}_{{\alpha }^0},x\in {\Omega }_x,\end{array}$ (57)

${\left\{B_V^{\left(3\right)}\left[4;U^{\left(3\right)}\left(4;x\right);V^{\left(3\right)}\left(4;x\right);\alpha ;\delta \alpha \right]\right\}}_{{\alpha }^0}=0\left[4\right];x\in \partial {\Omega }_x\left({\alpha }^0\right)$ (58)

where $0\left[4\right]\triangleq {\left[0,0,0,0\right]}^{†}$ and where:

$V{M}^{\left(3\right)}\left(4\times 4;U^{\left(3\right)};\alpha \right)\triangleq \left(\begin{array}{cc}V{M}^{\left(2\right)}\left(2\times 2\right)& 0\left[2\times 2\right]\\ V{M}_{21}^{\left(3\right)}\left(2\times 2\right)& V{M}_{22}^{\left(3\right)}\left(2\times 2\right)\end{array}\right);$ (59)

$U^{\left(3\right)}\left(4;x\right)\triangleq \left(\begin{array}{c}{U}^{\left(2\right)}\left(2;x\right)\\ A^{\left(2\right)}\left(2;j_1;x\right)\end{array}\right);\text{}\text{}\text{}\text{}\text{}{V}^{\left(3\right)}\left(4;x\right)\triangleq \delta U^{\left(3\right)}\left(4;x\right)=\left(\begin{array}{c}{V}^{\left(2\right)}\left(2;x\right)\\ \delta A^{\left(2\right)}\left(2;j_1;x\right)\end{array}\right);$ (60)

$V{M}_{21}^{\left(3\right)}\left(2\times 2;x\right)\triangleq \frac{\partial \left\{A{M}^{\left(2\right)}\left[2\times 2;U^{\left(2\right)};\alpha \right]A^{\left(2\right)}\left(2;x\right)\right\}}{\partial U^{\left(2\right)}\left(2;x\right)}-\frac{\partial Q_A^{\left(2\right)}\left[2;u^{\left(2\right)}\left(x\right);\alpha \right]}{\partial U^{\left(2\right)}\left(2;x\right)};$ (61)

$V{M}_{22}^{\left(3\right)}\left(2\times 2;x\right)\triangleq A{M}^{\left(2\right)}\left[2\times 2;U^{\left(2\right)};\alpha \right];0\left[2\times 2\right]\triangleq \left(\begin{array}{cc}0& 0\\ 0& 0\end{array}\right);$ (62)

$\begin{array}{l}{Q}_V^{\left(3\right)}\left[4;U^{\left(3\right)}\left(4;x\right);\alpha ;\delta \alpha \right]\triangleq \left(\begin{array}{c}{Q}_V^{\left(2\right)}\left[2;U^{\left(2\right)}\left(2;x\right);\alpha ;\delta \alpha \right]\\ Q_2^{\left(3\right)}\left[2;U^{\left(3\right)}\left(4;x\right);\alpha ;\delta \alpha \right]\end{array}\right)\\ \triangleq {\left\{q_V^{\left(3\right)}\left[1;U^{\left(3\right)}\left(4;x\right);\alpha ;\delta \alpha \right],\cdots ,q_V^{\left(3\right)}\left[4;U^{\left(3\right)}\left(4;x\right);\alpha ;\delta \alpha \right]\right\}}^{†};\end{array}$ (63)

$q_V^{\left(3\right)}\left[i;U^{\left(3\right)}\left(4;x\right);\alpha ;\delta \alpha \right]\equiv {\displaystyle \underset{j_3=1}{\overset{TP}{\sum }}{s}_V^{\left(3\right)}\left[i;j_3;U^{\left(3\right)}\left(4;x\right);\alpha \right]\delta {\alpha }_{j_3}};i=1,2,3,4;$ (64)

$\begin{array}{l}{Q}_2^{\left(3\right)}\left[2;U^{\left(3\right)}\left(4;x\right);\alpha ;\delta \alpha \right]\\ \triangleq \frac{\partial Q_A^{\left(2\right)}\left[2;u^{\left(2\right)}\left(x\right);\alpha \right]}{\partial \alpha }\partial \alpha -\frac{\partial \left\{A{M}^{\left(2\right)}\left[2\times 2;U^{\left(2\right)}\left(2;x\right);\alpha \right]A^{\left(2\right)}\left(2;j_1;x\right)\right\}}{\partial \alpha }\partial \alpha ;\end{array}$ (65)

$B_V^{\left(3\right)}\left[4;U^{\left(3\right)}\left(4;x\right);V^{\left(3\right)}\left(4;x\right);\alpha ;\delta \alpha \right]\triangleq \left(\begin{array}{c}{B}_V^{\left(2\right)}\left[2;U^{\left(2\right)}\left(2;x\right);V^{\left(2\right)}\left(2;x\right);\alpha ;\delta \alpha \right]\\ \delta B_A^{\left(2\right)}\left[2;U^{\left(3\right)}\left(4;x\right);V^{\left(3\right)}\left(4;x\right);\alpha ;\delta \alpha \right]\end{array}\right).$ (66)

The right-side of the 3^rd-LVSS actually depends on the indices $j_1,j_2,j_3=1,\cdots ,TP$, so the 3^rd-LVSS would need to be solved $T{P}^3$ times to obtain each of the variational functions $V^{\left(3\right)}\left(4;j_1,j_2,j_3;x\right)$. Thus, solving the 3^rd-LVSS would require $T{P}^3$ large-scale computations, which is unrealistic for large-scale systems comprising many parameters. Since the 3^rd-LVSS is never actually solved but is only used to construct the corresponding adjoint sensitivity system, the actual dependence of the 3^rd-LVSS on the indices $j_1,j_2,j_3=1,\cdots ,TP$ has been suppressed.

The 3^rd-CASAM-N circumvents the need for solving the 3^rd-LVSS by deriving an alternative expression for the indirect-effect term defined in Equation (56), in which the function $V^{\left(3\right)}\left(4;x\right)$ is replaced by a 3^rd-level adjoint function which is independent of parameter variations. This 3^rd-level adjoint function is the solution of a 3^rd-Level Adjoint Sensitivity System (3^rd-LASS) which is constructed by applying the same principles as those used for constructing the 1^st-LASS and the 2^nd-LASS. The Hilbert space appropriate for constructing the 3^rd-LASS, denoted as $H_3\left({\Omega }_x\right)$, comprises as elements block-vectors of the same form as $V^{\left(3\right)}\left(4;x\right)$. Thus, a generic block-vector in $H_3\left({\Omega }_x\right)$, denoted as ${\Psi }^{\left(3\right)}\left(4;x\right)\triangleq {\left[{\psi }^{\left(3\right)}\left(1;x\right),{\psi }^{\left(3\right)}\left(2;x\right),{\psi }^{\left(3\right)}\left(3;x\right),{\psi }^{\left(3\right)}\left(4;x\right)\right]}^{†}\in H_3\left({\Omega }_x\right)$, comprises four TD-dimensional vector-components of the form ${\psi }^{\left(3\right)}\left(i;x\right)\triangleq {\left[{\psi }_1^{\left(3\right)}\left(i;x\right),\cdots ,{\psi }_{TD}^{\left(3\right)}\left(i;x\right)\right]}^{†}\in H_1\left({\Omega }_x\right)$, $i=1,2,3,4$, where each of these four components is a TD-dimensional column vector. The inner product of two vectors ${\Psi }^{\left(3\right)}\left(4;x\right)\in H_3\left({\Omega }_x\right)$ and ${\Phi }^{\left(3\right)}\left(4;x\right)\in H_3\left({\Omega }_x\right)$ in the Hilbert space $H_3\left({\Omega }_x\right)$ will be denoted as ${\langle {\Psi }^{\left(3\right)}\left(4;x\right),{\Phi }^{\left(3\right)}\left(4;x\right)\rangle }_3$ and defined as follows:

${\langle {\Psi }^{\left(3\right)}\left(4;x\right),{\Phi }^{\left(3\right)}\left(4;x\right)\rangle }_3\triangleq {\displaystyle \underset{i=1}{\overset{4}{\sum }}{\langle {\psi }^{\left(3\right)}\left(i;x\right),{\phi }^{\left(3\right)}\left(i;x\right)\rangle }_1}$. (67)

The steps for constructing the 3^rd-LASS are conceptually similar to those described in Sections 2.1 and 2.2 and are detailed in [27]. The final expressions for the 3^rd-order sensitivities are as follows:

$\begin{array}{l}{\left\{\delta R^{\left(2\right)}\left[j_2;j_1;U^{\left(3\right)}\left(4;x\right);A^{\left(3\right)}\left(4;j_2;j_1;x\right);\alpha ;\delta \alpha \right]\right\}}_{{\alpha }^0}\\ ={\left\{\frac{\partial R^{\left(2\right)}\left[j_2;j_1;U^{\left(3\right)};\alpha \right]}{\partial \alpha }\delta \alpha \right\}}_{{\alpha }^0}-{\left\{{\left[{\hat{P}}^{\left(3\right)}\left(U^{\left(3\right)};A^{\left(3\right)};\delta \alpha \right)\right]}_{\partial {\Omega }_x}\right\}}_{{\alpha }^0}\\ +{\left\{{\langle A^{\left(3\right)}\left(4;j_2;j_1;x\right),Q_V^{\left(3\right)}\left[4;U^{\left(3\right)};\alpha ;\delta \alpha \right]\rangle }_3\right\}}_{{\alpha }^0},\end{array}$ (68)

where ${\left\{{\left[{\hat{P}}^{\left(3\right)}\left(U^{\left(3\right)};A^{\left(3\right)};\delta \alpha \right)\right]}_{\partial {\Omega }_x}\right\}}_{{\alpha }^0}$ denotes residual boundary terms which may have not vanished automatically, and where the 3^rd-level adjoint function $A^{\left(3\right)}\left(4;x\right)\triangleq {\left[a^{\left(3\right)}\left(1;x\right),a^{\left(3\right)}\left(2;x\right),a^{\left(3\right)}\left(3;x\right),a^{\left(3\right)}\left(4;x\right)\right]}^{†}\in H_3\left({\Omega }_x\right)$ is the solution of the following 3^rd-LASS:

$\begin{array}{l}{\left\{A{M}^{\left(3\right)}\left[4\times 4;U^{\left(3\right)}\left(4;x\right);\alpha \right]A^{\left(3\right)}\left(4;j_2;j_1;x\right)\right\}}_{{\alpha }^0}\\ ={\left\{Q_A^{\left(3\right)}\left[4;j_2;j_1;U^{\left(3\right)}\left(4;x\right);\alpha \right]\right\}}_{{\alpha }^0},j_1=1,\cdots ,TP;j_2=1,\cdots ,j_1;\end{array}$ (69)

where:

$Q_A^{\left(3\right)}\left[4;j_2;j_1;U^{\left(3\right)}\left(4;x\right);\alpha \right]\triangleq \left[q_A^{\left(3\right)}\left(1;j_2;j_1;U^{\left(3\right)};\alpha \right),\cdots ,q_A^{\left(3\right)}\left(4;j_2;j_1;U^{\left(3\right)};\alpha \right)\right];$ (70)

$q_A^{\left(3\right)}\left(1;j_2;j_1;U^{\left(3\right)};\alpha \right)\triangleq \partial R^{\left(2\right)}\left[j_2;j_1;u^{\left(2\right)};a^{\left(2\right)};\alpha \right]/\partial u^{\left(1\right)};$ (71)

$q_A^{\left(3\right)}\left(2;j_2;j_1;U^{\left(3\right)};\alpha \right)\triangleq \partial R^{\left(2\right)}\left[j_2;j_1;u^{\left(2\right)};a^{\left(2\right)};\alpha \right]/\partial a^{\left(1\right)};$ (72)

$q_A^{\left(3\right)}\left(3;j_2;j_1;U^{\left(3\right)};\alpha \right)\triangleq \partial R^{\left(2\right)}\left[j_2;j_1;u^{\left(2\right)};a^{\left(2\right)};\alpha \right]/\partial a^{\left(2\right)}\left(1;j_1;x\right);$ (73)

$q_A^{\left(3\right)}\left(3;j_2;j_1;U^{\left(3\right)};\alpha \right)\triangleq \partial R^{\left(2\right)}\left[j_2;j_1;u^{\left(2\right)};a^{\left(2\right)};\alpha \right]/\partial a^{\left(2\right)}\left(2;j_1;x\right).$ (74)

$\begin{array}{l}A{M}^{\left(3\right)}\left[4\times 4;U^{\left(3\right)}\left(4;x\right);\alpha \right]\triangleq {\left[V{M}^{\left(3\right)}\left(4\times 4;U^{\left(3\right)};\alpha \right)\right]}^{*}\\ =\left(\begin{array}{cc}{\left\{{\left[V{M}^{\left(2\right)}\left(2\times 2\right)\right]}^{*}\right\}}^{†}& {\left\{{\left[V{M}_{21}^{\left(3\right)}\left(2\times 2\right)\right]}^{*}\right\}}^{†}\\ 0\left[2\times 2\right]& {\left\{{\left[V{M}_{22}^{\left(3\right)}\left(2\times 2\right)\right]}^{*}\right\}}^{†}\end{array}\right),\end{array}$ (75)

The boundary conditions to be satisfied by each of the 3^rd-level adjoint functions $\begin{array}{l}{A}^{\left(3\right)}\left(4;j_2;j_1;x\right)\\ \triangleq {\left[a^{\left(3\right)}\left(1;j_2;j_1;x\right),a^{\left(3\right)}\left(2;j_2;j_1;x\right),a^{\left(3\right)}\left(3;j_2;j_1;x\right),a^{\left(3\right)}\left(4;j_2;j_1;x\right)\right]}^{†}\end{array}$ can be represented in operator form as follows:

$\begin{array}{l}{\left\{B_A^{\left(3\right)}\left[4;U^{\left(3\right)}\left(4;x\right);A^{\left(3\right)}\left(4;j_2;j_1;x\right);\alpha \right]\right\}}_{{\alpha }^0}=0\left[4\right];\\ \text{for}{j}_1=1,\cdots ,TP;j_2=1,\cdots ,j_1;x\in \partial {\Omega }_x\left({\alpha }^0\right).\end{array}$ (76)

In component form, the total differential expressed by Equation (68) has the following expression:

$\begin{array}{l}{\left\{\delta R^{\left(2\right)}\left[j_2;j_1;U^{\left(3\right)}\left(4;x\right);A^{\left(3\right)}\left(4;j_2;j_1;x\right);\alpha ;\delta \alpha \right]\right\}}_{{\alpha }^0}\\ ={\displaystyle \underset{j_3=1}{\overset{TP}{\sum }}{\left\{R^{\left(3\right)}\left[j_3;j_2;j_1;U^{\left(3\right)}\left(4;x\right);A^{\left(3\right)}\left(4;j_2;j_1;x\right);\alpha \right]\right\}}_{{\alpha }^0}}\delta {\alpha }_{j_3},j_1;j_2=1,\cdots ,TP,\end{array}$ (77)

where the quantity $R^{\left(3\right)}\left[j_3;j_2;j_1;U^{\left(3\right)}\left(4;x\right);A^{\left(3\right)}\left(4;j_2;j_1;x\right);\alpha \right]$ denotes the third-order sensitivity of the generic scalar-valued response $R\left[u\left(x\right);\alpha \right]$ with respect to any three model parameters ${\alpha }_{j_1}$, ${\alpha }_{j_2}$, ${\alpha }_{j_3}$, and has the following expression:

$\begin{array}{l}{R}^{\left(3\right)}\left[j_3;j_2;j_1;U^{\left(3\right)}\left(4;x\right);A^{\left(3\right)}\left(4;j_2;j_1;x\right);\alpha \right]\\ \triangleq \frac{\partial R^{\left(2\right)}\left[j_2;j_1;U^{\left(3\right)}\left(4;j_1;x\right);\alpha \right]}{\partial {\alpha }_{j_3}}-\frac{{\left[\partial {\hat{P}}^{\left(3\right)}\left(U^{\left(3\right)};A^{\left(3\right)};\delta \alpha \right)\right]}_{\partial {\Omega }_x}}{\partial {\alpha }_{j_3}}\\ +{\displaystyle \underset{i=1}{\overset{4}{\sum }}{\langle a^{\left(3\right)}\left(i;j_2;j_1;x\right),s_V^{\left(3\right)}\left[i;j_3;j_1;U^{\left(3\right)}\left(4;x\right);\alpha \right]\rangle }_1}\\ \triangleq \frac{{\partial }^{3}R\left[u\left(x\right);\alpha \right]}{\partial {\alpha }_{j_3}\partial {\alpha }_{j_2}\partial {\alpha }_{j_1}},\text{for}{j}_1,j_2,j_3=1,\cdots ,TP.\end{array}$ (78)

2.4. The Fourth-Order Comprehensive Adjoint Sensitivity Analysis Methodology for Nonlinear Systems (4^th-CASAM-N)

Assuming that the 3^rd-order sensitivities $R^{\left(3\right)}\left[j_3;j_2;j_1;U^{\left(3\right)};A^{\left(3\right)};\alpha \right]\equiv {\partial }^{3}R\left[u\left(x\right);\alpha \right]/\partial {\alpha }_{j_3}\partial {\alpha }_{j_2}\partial {\alpha }_{j_1}$ satisfy the conditions stated in Equations (6) and (7) for each $j_1,j_2,j_3=1,\cdots ,TP$, the 1^st-order total G-differential of $R^{\left(3\right)}\left[j_3;j_2;j_1;U^{\left(3\right)};A^{\left(3\right)};\alpha \right]$ will exist and will be linear in the variations $V^{\left(3\right)}\left(4;x\right)\triangleq \delta U^{\left(3\right)}\left(4;x\right)$ and $\delta A^{\left(3\right)}\left(4;x\right)$ in a neighborhood around the nominal values of the parameters and the respective state functions. By definition, the 1^st-order total G-differential of $R^{\left(3\right)}\left[j_3;j_2;j_1;U^{\left(3\right)}\left(4;j_1;x\right);A^{\left(3\right)}\left(4;j_2;j_1;x\right);\alpha \right]$, which will be denoted as ${\left\{\delta R^{\left(3\right)}\left[j_3;j_2;j_1;U^{\left(3\right)}\left(4;j_1;x\right);A^{\left(3\right)}\left(4;j_2;j_1;x\right);\alpha ;V^{\left(3\right)}\left(4;j_1;x\right);\delta A^{\left(3\right)}\left(4;j_2;j_1;x\right);\delta \alpha \right]\right\}}_{{\alpha }^0}$,

is given by the following expression:

$\begin{array}{l}{\left\{\delta R^{\left(3\right)}\left[j_3;j_2;j_1;U^{\left(3\right)}\left(4;j_1;x\right);A^{\left(3\right)}\left(4;j_2;j_1;x\right);\alpha ;V^{\left(3\right)}\left(4;j_1;x\right);\delta A^{\left(3\right)}\left(4;j_2;j_1;x\right);\delta \alpha \right]\right\}}_{{\alpha }^0}\\ \triangleq {\left\{\delta R^{\left(3\right)}\left[j_3;j_2;j_1;U^{\left(3\right)}\left(4;j_1;x\right);A^{\left(3\right)}\left(4;j_2;j_1;x\right);\alpha ;\delta \alpha \right]\right\}}_{dir}\\ +{\left\{\delta R^{\left(3\right)}\left[j_3;j_2;j_1;U^{\left(3\right)}\left(4;j_1;x\right);A^{\left(3\right)}\left(4;j_2;j_1;x\right);\alpha ;V^{\left(3\right)}\left(4;j_1;x\right);\delta A^{\left(3\right)}\left(4;j_2;j_1;x\right)\right]\right\}}_{ind}.\end{array}$ (79)

where:

$\begin{array}{l}{\left\{\delta R^{\left(3\right)}\left[j_3;j_2;j_1;U^{\left(3\right)}\left(4;j_1;x\right);A^{\left(3\right)}\left(4;j_2;j_1;x\right);\alpha ;\delta \alpha \right]\right\}}_{dir}\\ \triangleq {\left\{\frac{\partial R^{\left(3\right)}\left[j_3;j_2;j_1;U^{\left(3\right)}\left(4;j_1;x\right);A^{\left(3\right)}\left(4;j_2;j_1;x\right);\alpha \right]}{\partial \alpha }\delta \alpha \right\}}_{{\alpha }^0},\end{array}$ (80)

$\begin{array}{l}{\left\{\delta R^{\left(3\right)}\left[j_3;j_2;j_1;U^{\left(3\right)}\left(4;j_1;x\right);A^{\left(3\right)}\left(4;j_2;j_1;x\right);\alpha ;V^{\left(3\right)}\left(4;j_1;x\right);\delta A^{\left(3\right)}\left(4;j_2;j_1;x\right)\right]\right\}}_{ind}\\ \triangleq {\left\{\frac{\partial R^{\left(3\right)}\left[j_3;j_2;j_1;U^{\left(3\right)};A^{\left(3\right)};\alpha \right]}{\partial U^{\left(3\right)}\left(4;j_1;x\right)}\right\}}_{{\alpha }^0}{V}^{\left(3\right)}\left(4;j_1;x\right)\\ +{\left\{\frac{\partial R^{\left(3\right)}\left[j_3;j_2;j_1;U^{\left(3\right)};A^{\left(3\right)};\alpha \right]}{\partial A^{\left(3\right)}\left(4;j_2;j_1;x\right)}\right\}}_{{\alpha }^0}\delta A^{\left(3\right)}\left(4;j_2;j_1;x\right),\end{array}$ (81)

The vectors $V^{\left(3\right)}\left(4;x\right)$ and $\delta A^{\left(3\right)}\left(4;x\right)$ are the solutions of the following “4^th-order variational sensitivity system” (4^th-LVSS), which is derived in detail in [27]:

${\left\{V{M}^{\left(4\right)}\left[8\times 8;U^{\left(4\right)};\alpha \right]V^{\left(4\right)}\left(8;x\right)\right\}}_{{\alpha }^0}={\left\{Q_V^{\left(4\right)}\left[8;U^{\left(4\right)}\left(8;x\right);\alpha ;\delta \alpha \right]\right\}}_{{\alpha }^0},x\in {\Omega }_x,$ (82)

$\begin{array}{l}{\left\{B_V^{\left(4\right)}\left[8;U^{\left(4\right)}\left(8;x\right);V^{\left(4\right)}\left(8;x\right);\alpha ;\delta \alpha \right]\right\}}_{{\alpha }^0}\\ =0\left[8\right]\triangleq {\left[0,0,0,0,0,0,0,0\right]}^{†};x\in \partial {\Omega }_x\left({\alpha }^0\right),\end{array}$ (83)

where

$V{M}^{\left(4\right)}\left(8\times 8\right)\triangleq \left(\begin{array}{cc}V{M}^{\left(3\right)}\left(4\times 4\right)& 0\left[4\times 4\right]\\ V{M}_{21}^{\left(4\right)}\left(4\times 4\right)& V{M}_{22}^{\left(4\right)}\left(4\times 4\right)\end{array}\right);$ (84)

$U^{\left(4\right)}\left(8;j_2;j_1;x\right)\triangleq \left(\begin{array}{c}{U}^{\left(3\right)}\left(4;j_1;x\right)\\ A^{\left(3\right)}\left(4;j_2;j_1;x\right)\end{array}\right);$ (85)

$V^{\left(4\right)}\left(8;j_2;j_1;x\right)\triangleq \delta U^{\left(4\right)}\left(8;j_2;j_1;x\right)=\left(\begin{array}{c}{V}^{\left(3\right)}\left(4;j_1;x\right)\\ \delta A^{\left(3\right)}\left(4;j_2;j_1;x\right)\end{array}\right)$

$\begin{array}{l}=[v^{\left(1\right)}\left(x\right),\delta a^{\left(1\right)}\left(x\right),\delta a^{\left(2\right)}\left(1;j_1;x\right),\delta a^{\left(2\right)}\left(2;j_1;x\right),\delta a^{\left(3\right)}\left(1;j_2;j_1;x\right),\\ {\cdots ,\delta a^{\left(3\right)}\left(4;j_2;j_1;x\right)]}^{†};\end{array}$ (86)

$V{M}_{21}^{\left(4\right)}\left(4\times 4;x\right)\triangleq \frac{\partial \left\{A{M}^{\left(3\right)}\left(4\times 4;U^{\left(3\right)};\alpha \right)A^{\left(3\right)}\left(4;x\right)\right\}}{\partial U^{\left(3\right)}\left(4;j_1;x\right)}-\frac{\partial Q_A^{\left(3\right)}\left[4;U^{\left(3\right)}\left(4;x\right);\alpha \right]}{\partial U^{\left(3\right)}\left(4;j_1;x\right)};$ (87)

$V{M}_{22}^{\left(4\right)}\left(4\times 4;j_1;x\right)\triangleq A{M}^{\left(3\right)}\left[4\times 4;U^{\left(3\right)}\left(4;j_1;x\right);\alpha \right];$ (88)

$\begin{array}{l}{Q}_V^{\left(4\right)}\left[8;U^{\left(4\right)}\left(4;x\right);\alpha ;\delta \alpha \right]\triangleq \left(\begin{array}{c}{Q}_V^{\left(3\right)}\left[4;U^{\left(3\right)}\left(4;x\right);\alpha ;\delta \alpha \right]\\ Q_2^{\left(4\right)}\left[4;U^{\left(4\right)}\left(4;x\right);\alpha ;\delta \alpha \right]\end{array}\right)\\ \triangleq {\left\{q_V^{\left(4\right)}\left[1;U^{\left(4\right)}\left(4;x\right);\alpha ;\delta \alpha \right],\cdots ,q_V^{\left(4\right)}\left[8;U^{\left(4\right)}\left(4;x\right);\alpha ;\delta \alpha \right]\right\}}^{†};\end{array}$ (89)

$\begin{array}{l}{q}_V^{\left(4\right)}\left[i;U^{\left(4\right)}\left(4;j_2;j_1;x\right);\alpha ;\delta \alpha \right]\\ \equiv {\displaystyle \underset{j_4=1}{\overset{TP}{\sum }}{s}_V^{\left(4\right)}\left[i;j_4;U^{\left(4\right)}\left(4;j_2;j_1;x\right);\alpha \right]\delta {\alpha }_{j_4}};i=1,\cdots ,8;\end{array}$ (90)

$\begin{array}{l}{Q}_2^{\left(4\right)}\left[4;U^{\left(4\right)}\left(4;j_2;j_1;x\right);\alpha ;\delta \alpha \right]\triangleq \frac{\partial Q_A^{\left(3\right)}\left[4;j_2;j_1;U^{\left(3\right)}\left(4;j_1;x\right);\alpha \right]}{\partial \alpha }\partial \alpha \\ -\frac{\partial \left\{A{M}^{\left(3\right)}\left[4\times 4;U^{\left(3\right)}\left(4;j_1;x\right);\alpha \right]A^{\left(3\right)}\left(4;j_2;j_1;x\right)\right\}}{\partial \alpha }\partial \alpha ;\end{array}$ (91)

$\begin{array}{l}{B}_V^{\left(4\right)}\left[8;U^{\left(4\right)}\left(8;j_2;j_1;x\right);V^{\left(4\right)}\left(8;j_2;j_1;x\right);\alpha ;\delta \alpha \right]\\ \triangleq \left(\begin{array}{c}{B}_V^{\left(3\right)}\left[4;U^{\left(3\right)}\left(4;j_1;x\right);V^{\left(3\right)}\left(4;j_1;x\right);\alpha ;\delta \alpha \right]\\ \delta B_A^{\left(3\right)}\left[4;U^{\left(4\right)}\left(8;j_2;j_1;x\right);V^{\left(4\right)}\left(8;j_2;j_1;x\right);\alpha ;\delta \alpha \right]\end{array}\right).\end{array}$ (92)

The right-side of the 4^th-LVSS actually depends on the indices $j_1,j_2,j_3,j_4=1,\cdots ,TP$, so the 4^th-LVSS would need to be solved $T{P}^4$ times in order to obtain each of the variational functions $V^{\left(4\right)}\left(4;j_1,j_2,j_3,j_4;x\right)$, which is unrealistic for large-scale systems comprising many parameters. Since the 4^th-LVSS is never actually solved but is only used to construct the corresponding adjoint sensitivity system, the actual dependence of the 4^th-LVSS on the indices $j_1,j_2,j_3=1,\cdots ,TP$ has been suppressed.

The 4^th-CASAM-N circumvents the need for solving the 4^th-LVSS by deriving an alternative expression for the indirect-effect term defined in Equation (81), in which the function $V^{\left(4\right)}\left(8;x\right)$ is replaced by a 4^th-level adjoint function which is independent of parameter variations. This 4^th-level adjoint function will be the solution of a 4^th-Level Adjoint Sensitivity System (4^th-LASS) which will be constructed by applying the same principles as those used for constructing the 1^st-LASS, the 2^nd-LASS and the 3^rd-LASS. The Hilbert space appropriate for constructing the 4^th-LASS will be denoted as $H_4\left({\Omega }_x\right)$ and comprises as elements block-vectors of the same form as $V^{\left(4\right)}\left(8;j_2;j_1;x\right)$. Thus, a generic block-vector in $H_4\left({\Omega }_x\right)$ will have the structure ${\Psi }^{\left(4\right)}\left(8;x\right)\triangleq {\left[{\psi }^{\left(4\right)}\left(1;x\right),\cdots ,{\psi }^{\left(4\right)}\left(8;x\right)\right]}^{†}\in H_4\left({\Omega }_x\right)$, comprising 8 TD-dimensional vectors of the form ${\psi }^{\left(4\right)}\left(i;x\right)\triangleq {\left[{\psi }_1^{\left(4\right)}\left(i;x\right),\cdots ,{\psi }_{TD}^{\left(4\right)}\left(i;x\right)\right]}^{†}\in H_1\left({\Omega }_x\right)$, $i=1,\cdots ,8$. The inner product of two vectors ${\Psi }^{\left(4\right)}\left(8;x\right)\in H_4\left({\Omega }_x\right)$ and ${\Phi }^{\left(4\right)}\left(8;x\right)\in H_4\left({\Omega }_x\right)$ in the Hilbert space $H_4\left({\Omega }_x\right)$ will be denoted as ${\langle {\Psi }^{\left(4\right)}\left(8;x\right),{\Phi }^{\left(4\right)}\left(8;x\right)\rangle }_4$ and defined as follows:

${\langle {\Psi }^{\left(4\right)}\left(8;x\right),{\Phi }^{\left(4\right)}\left(8;x\right)\rangle }_4\triangleq {\displaystyle \underset{i=1}{\overset{8}{\sum }}{\langle {\psi }^{\left(4\right)}\left(i;x\right),{\phi }^{\left(4\right)}\left(i;x\right)\rangle }_1}$. (93)

The steps for constructing the 4^th-LASS are conceptually similar to those described in Sections 2.1 and 2.2 and are detailed in [27]. The final expressions for the 4^th-order sensitivities are as follows:

$\begin{array}{l}{\left\{\delta R^{\left(3\right)}\left[j_3;j_2;j_1;U^{\left(4\right)}\left(8;j_2;j_1;x\right);A^{\left(4\right)}\left(8;j_3;j_2;j_1;x\right);\delta \alpha \right]\right\}}_{{\alpha }^0}\\ ={\left\{\frac{\partial R^{\left(3\right)}\left[j_3;j_2;j_1;U^{\left(3\right)}\left(4;j_1;x\right);A^{\left(3\right)}\left(4;j_2;j_1;x\right);\alpha \right]}{\partial \alpha }\delta \alpha \right\}}_{{\alpha }^0}\\ -{\left\{{\left[{\hat{P}}^{\left(4\right)}\left(U^{\left(4\right)};A^{\left(4\right)};\delta \alpha \right)\right]}_{\partial {\Omega }_x}\right\}}_{{\alpha }^0}\\ +{\left\{{\langle A^{\left(4\right)}\left(8;j_3;j_2;j_1;x\right),Q_V^{\left(4\right)}\left[8;j_1;U^{\left(4\right)}\left(8;j_2;j_1;x\right);\alpha ;\delta \alpha \right]\rangle }_3\right\}}_{{\alpha }^0}.\end{array}$ (94)

In component form, the total differential expressed by Equation (94) can be written in the following form, for each $j_1=1,\cdots ,TP$ ; $j_2=1,\cdots ,j_1$ ; $j_3=1,\cdots ,j_2$ :

$\begin{array}{l}{\left\{\delta R^{\left(3\right)}\left[j_3;j_2;j_1;U^{\left(4\right)}\left(8;j_2;j_1;x\right);A^{\left(4\right)}\left(8;j_3;j_2;j_1;x\right);\delta \alpha \right]\right\}}_{{\alpha }^0}\\ ={\displaystyle \underset{j_4=1}{\overset{TP}{\sum }}{\left\{R^{\left(4\right)}\left[j_4;j_3;j_2;j_1;U^{\left(4\right)}\left(8;j_2;j_1;x\right);A^{\left(4\right)}\left(8;j_3;j_2;j_1;x\right);\alpha \right]\right\}}_{{\alpha }^0}}\delta {\alpha }_{j_4},\end{array}$ (95)

where the quantity $R^{\left(4\right)}\left[j_4;j_3;j_2;j_1;U^{\left(4\right)}\left(8;j_2;j_1;x\right);A^{\left(4\right)}\left(8;j_3;j_2;j_1;x\right);\alpha \right]$ denotes the fourth-order sensitivity of the generic scalar-valued response $R\left[u\left(x\right);\alpha \right]$ with respect to any four model parameters ${\alpha }_{j_1}$, ${\alpha }_{j_2}$, ${\alpha }_{j_3}$, ${\alpha }_{j_4}$, and has the following expression:

$\begin{array}{l}{R}^{\left(4\right)}\left[j_4;j_3;j_2;j_1;U^{\left(4\right)}\left(8;j_2;j_1;x\right);A^{\left(4\right)}\left(8;j_3;j_2;j_1;x\right);\alpha \right]\\ \triangleq \frac{\partial R^{\left(3\right)}\left[j_3;j_2;j_1;U^{\left(3\right)}\left(4;j_1;x\right);A^{\left(3\right)}\left(4;j_2;j_1;x\right);\alpha \right]}{\partial {\alpha }_{j_4}}-\frac{{\left[\partial {\hat{P}}^{\left(4\right)}\left(U^{\left(4\right)};A^{\left(4\right)};\delta \alpha \right)\right]}_{\partial {\Omega }_x}}{\partial {\alpha }_{j_4}}\\ +{\displaystyle \underset{i=1}{\overset{8}{\sum }}{\langle a^{\left(4\right)}\left(i;j_3;j_2;j_1;x\right),s_V^{\left(4\right)}\left[i;j_4;j_1;U^{\left(4\right)}\left(4;j_2;j_1;x\right);\alpha \right]\rangle }_1}\\ \equiv {\partial }^{4}R\left[u\left(x\right);\alpha \right]/\partial {\alpha }_{j_{}}\end{array}$