Introducing the n^th-Order Features Adjoint Sensitivity Analysis Methodology for Nonlinear Systems (n^th-FASAM-N): I. Mathematical Framework

Dan Gabriel Cacuci

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

DOI: 10.4236/ajcm.2024.141002                

Abstract

This work presents the “n^th-Order Feature Adjoint Sensitivity Analysis Methodology for Nonlinear Systems” (abbreviated as “n^th-FASAM-N”), which will be shown to be the most efficient methodology for computing exact expressions of sensitivities, of any order, of model responses with respect to features of model parameters and, subsequently, with respect to the model’s uncertain parameters, boundaries, and internal interfaces. The unparalleled efficiency and accuracy of the n^th-FASAM-N methodology stems from the maximal reduction of the number of adjoint computations (which are considered to be “large-scale” computations) for computing high-order sensitivities. When applying the n^th-FASAM-N methodology to compute the second- and higher-order sensitivities, the number of large-scale computations is proportional to the number of “model features” as opposed to being proportional to the number of model parameters (which are considerably more than the number of features).When a model has no “feature” functions of parameters, but only comprises primary parameters, the n^th-FASAM-N methodology becomes identical to the extant n^th CASAM-N (“n^th-Order Comprehensive Adjoint Sensitivity Analysis Methodology for Nonlinear Systems”) methodology. Both the n^th-FASAM-N and the n^th-CASAM-N methodologies are formulated in linearly increasing higher-dimensional Hilbert spaces as opposed to exponentially increasing parameter-dimensional spaces thus overcoming the curse of dimensionality in sensitivity analysis of nonlinear systems. Both the n^th-FASAM-N and the n^th-CASAM-N are incomparably more efficient and more accurate than any other methods (statistical, finite differences, etc.) for computing exact expressions of response sensitivities of any order with respect to the model’s features and/or primary uncertain parameters, boundaries, and internal interfaces.

Keywords

Computation of High-Order Sensitivities, Sensitivities to Features of Model Parameters, Sensitivities to Domain Boundaries, Adjoint Sensitivity Systems

1. Introduction

The computational model of a physical system comprises the following fundamental conceptual components: (a) a well-posed system of equations that relate the system’s independent variables and parameters to the system’s state (i.e., dependent) variables; (b) nominal or mean values for the system parameters, along with information regarding the probability distributions, moments thereof, inequality and/or equality constraints that define the range of variations of the system’s parameters; (c) one or several quantities, customarily referred to as system responses (or objective functions, or indices of performance), which are computed using the mathematical model. Since the physical processes themselves are seldom known precisely and since most of the model’s parameters stem from experimental procedures that are also subject to imprecisions and/or uncertainties, the results predicted by these models are also imprecise, being affected by the uncertainties underlying the respective model. The analysis of the accuracy of responses computed by models relies fundamentally on the functional derivatives (also called “sensitivities”) of the respective responses with respect to the imprecisely known parameters underlying the computational model. Such sensitivities are needed for many purposes, including: (i) understanding the model by ranking the importance of the various parameters; (ii) performing “reduced-order modeling” by eliminating unimportant parameters and/or processes; (iii) quantifying the uncertainties induced in a model response due to model parameter uncertainties; (iv) performing “model validation,” by comparing computations to experiments to address the question “does the model represent reality?” (v) prioritizing improvements in the model; (vi) performing data assimilation and model calibration as part of forward “predictive modeling” to obtain best-estimate predicted results with reduced predicted uncertainties; (vii) performing inverse “predictive modeling”; (viii) designing and optimizing the system.

Response sensitivities are computed by using either deterministic or statistical methods. The “statistical methods” construct an approximate response distribution (often called “response surface”) in the parameters space by performing many “forward” computations using the model with altered parameter values, and subsequently use scatter plots, regression, rank transformation, correlations, and/or so-called “partial correlation analysis,” in order to identify approximate expectation values, variances and covariances for the responses. These statistical quantities are subsequently used to construct quantities that play the role of (approximate) first-order response sensitivities. Thus, statistical methods commence with “uncertainty analysis” and subsequently attempt an approximate “sensitivity analysis” of the approximately computed model response (called a “response surface”) in the phase-space of the parameters under consideration. The currently popular statistical methods for uncertainty and sensitivity analysis are broadly categorized as variance-based methods (see, e.g., [1] [2] [3] ), sampling-based methods (see, e.g., [4] [5] ), or Bayesian methods (see, e.g., [6] ). Various variants of the statistical methods for uncertainty and sensitivity analysis are reviewed in the book edited by Saltarelli et al. [7] . The conclusions that emerge from examining these statistical methods are as follows:

1) The main advantage of using statistical methods for uncertainty and sensitivity analysis is that they are conceptually easy to implement.

2) Even first-order sensitivities cannot be computed exactly.

3) Statistical methods are subject to the curse of dimensionality and have not been developed for producing higher-order sensitivities.

4) Since the response sensitivities and parameter uncertainties are amalgamated, inherently and inseparably, within the results produced by statistical methods, improvements in parameter uncertainties cannot be directly propagated to improve response uncertainties; rather, the entire set of simulations and statistical post-processing must be repeated anew.

5) A “fool-proof” statistical method for analyzing correctly models involving highly correlated parameters does not seem to exist currently so particular care must be used when interpreting regression results obtained using such models.

The simplest deterministic method for computing response sensitivities is to use finite-difference schemes in conjunction with re-computations using the model with “judiciously chosen” altered parameter values. Evidently, such methods can at best compute approximate values of a very limited number of sensitivities. The earliest deterministic methods that could compute more exactly the values of first-order sensitivities include the “Green’s function method” [8] , the “forward sensitivity analysis methodology” [9] , and the “direct method” [10] , which rely on analytical or numerical differentiation of the computational model under investigation to compute local response sensitivities exactly. However, for a computational model comprising many parameters, the conventional deterministic methods become impractical for computing sensitivities higher than first-order because they are subject to the “curse of dimensionality,” a term coined by Belmann [11] to describe phenomena in which the number of computations increases exponentially in the respective phase-space. In the particular case of sensitivity analysis using conventional deterministic methods, the number of large-scale computations increases exponentially in the phase-space of the model parameter as the order of sensitivities increases.

It is known that the “adjoint method of sensitivity analysis” has been the most efficient method for computing exactly first-order sensitivities, since it requires a single large-scale (adjoint) computation for computing all of the first-order sensitivities, regardless of the number of model parameters. The idea underlying the computation of response sensitivities with respect to model parameters using adjoint operators was first used by Wigner [12] to analyze first-order perturbations in nuclear reactor physics and shielding models based on the linear neutron transport (or diffusion) equation, as subsequently described in textbooks on these subjects (see, e.g., [13] [14] [15] [16] [17] ). Cacuci [18] is given credit (see, e.g., [19] [20] ) for having conceived the rigorous mathematical framework of the “1st-order adjoint sensitivity analysis methodology” for generic large-scale nonlinear (as opposed to linearized) systems involving generic operator responses, and for having introduced these principles to the earth, atmospheric and other sciences.

Cacuci [21] [22] has extended his 1st-order adjoint sensitivity analysis methodology to enable the comprehensive computation of 2nd-order sensitivities of model responses to model parameters (including imprecisely known domain boundaries and interfaces) for large-scale linear and nonlinear systems. The unparalleled efficiency of the 2nd-order adjoint sensitivity analysis methodology for linear systems [21] was demonstrated by Cacuci and Fang [23] by applying this methodology to compute exactly the 21,976 first-order sensitivities and 482,944,576 second-order sensitivities (of which 241,483,276 are distinct from each other) for an OECD/NEA reactor physics benchmark [24] . This benchmark is modeled by the neutron transport equation involving 21,976 uncertain parameters, the solving of which is representative of “large-scale computations.” The neutron transport equation was solved using the software package PARTISN [25] in conjunction with the MENDF71X cross section library [26] , which comprises 618-group cross sections based on ENDF/B-VII.1 nuclear data [27] . The spontaneous fission source has been computed using the code SOURCES4C [28] .

Contrary to the widely held belief that second- and higher-order sensitivities are negligible for reactor physics systems, Cacuci and Fang [23] found that many 2nd-order sensitivities of the OECD benchmark’s response to the benchmark’s uncertain parameters were much larger than the largest 1st-order ones, which motivated the investigation of the largest 3rd-order sensitivities, many of which were found to be even larger than the 2nd-order ones. This finding has motivated the development of the mathematical framework for determining and computing the 4th-order sensitivities, many of which were found to be larger than the 3rd-order ones. This sequence of findings has motivated the development by Cacuci [29] of the “nth-Order Comprehensive Adjoint Sensitivity Analysis Methodology for Response-Coupled Forward/Adjoint Linear Systems” (which is abbreviated as “nth-CASAM-L”). The “nth-CASAM-L” mathematical framework was developed specifically for linear systems because important model responses produced by such systems are various Lagrangian functionals such as eigenvalues or “CONTRIBUTONS” [30] , which depend simultaneously on both the forward and adjoint state functions governing the respective linear system. Such responses can occur only for linear systems because responses in nonlinear systems can only depend on the system’s forward state functions, since nonlinear operators do not admit adjoint operators. The nth-CASAM-L overcomes the curse of dimensionality in sensitivity analysis of linear systems, enabling the efficient computation of exactly-determined expressions of arbitrarily high-order sensitivities of a generic system response (that can depend on both the forward and adjoint state functions) with respect to all of the parameters that characterize the physical system.

In parallel with developing the nth-CASAM-L, Cacuci [31] has also developed the nth-Order Comprehensive Adjoint Sensitivity Analysis Methodology for Nonlinear Systems (nth-CASAM-N). Just like the nth-CASAM-L, the nth-CASAM-N is also formulated in linearly increasing higher-dimensional Hilbert spaces (as opposed to exponentially increasing parameter-dimensional spaces), thus overcoming the curse of dimensionality in sensitivity analysis of nonlinear systems, enabling the most efficient computation of exactly-determined expressions of arbitrarily high-order sensitivities of generic nonlinear system responses with respect to model parameters, uncertain boundaries and internal interfaces in the model’s phase-space.

Recently, Cacuci [32] has introduced the “Second-Order Function/Feature Adjoint Sensitivity Analysis Methodology for Nonlinear Systems” (2nd-FASAM-N), which enables a considerable reduction (by comparison to the 2nd-CASAM-N) of the number of large-scale computations needed to compute the second-order sensitivities of a model response with respect to the model parameters. The construction of the 2nd-FASAM-N is based on the same principles as those underlying the construction of the 2nd-CASAM-N. This fact indicates that the principles which enabled the construction of the nth-CASAM-N methodology [31] could also be used to generalize the 2nd-FASAM-N to enable the most efficient computation of the exact expressions of arbitrarily-high (“nth”) order sensitivities of model responses with respect to features/functions of model parameters. This is indeed the case, as will be shown in Section 2 of this work, which presents the construction of the general mathematical framework underlying the nth-FASAM-N. This construction uses “mathematical induction” mirroring the construction of the nth-CASAM-N methodology [31] . Section 3 concludes this work by discussing the significance of the nth-CASAM-N and preparing the groundwork for a paradigm illustrative application to the Nordheim-Fuchs reactor dynamics/safety model [33] [34] to be presented in the accompanying “Part 2” [35] .

2. The nth-Order Feature/Functions Adjoint Sensitivity Analysis Methodology for Nonlinear Systems (nth-FASAM-N) Methodology

The generic mathematical model of a nonlinear system which comprises uncertain (i.e., imprecisely known) model parameters, domain boundaries and interfaces, is the same as was introduced by Cacuci [31] and is reproduced in Appendix A, for convenient referencing.

The validity of the mathematical methodology underlying the nth-FASAM-N methodology will be established in this Section by using the “proof by mathematical induction” comprising the usual steps, as follows:

1) Surmise the general pattern underlying the nth-FASAM-N methodology for an arbitrarily high-order denoted as “n.” It is expected that the general pattern underlying the nth-FASAM-N can be surmised based on the expected similarities with the general pattern underlying the nth-CASAM-N [31] .

2) Prove that the general pattern underlying the nth-FASAM-N is valid for the lowest values of n, i.e., n = 1, which should reduce to the mathematical framework underlying the 1st-FASAM-N, which is reproduced from Cacuci [32] in Appendix B.

3) Assuming that the pattern underlying the nth-FASAM-N is valid for an arbitrarily high-order, n, prove that this pattern is also valid for n+1, i.e., for the (n+1)th-FASAM-N.

2.1. Establishing the Mathematical Framework of the nth-Order Feature Adjoint Sensitivity Analysis Methodology for Nonlinear Systems (nth-FASAM-N)

Comparing the mathematical framework of the 1^st-FASAM-N [32] to the framework of the 1^st-CASAM-N [31] suggests that the component “features” $f_i\left(\alpha \right),i=1,\cdots ,TF$ , of the vector-valued “feature function” $f\left(\alpha \right)\triangleq {\left[f_1\left(\alpha \right),\cdots ,f_{TF}\left(\alpha \right)\right]}^{†}$ , where TF denotes the total number of the components, play within the 1^st-FASAM-N the same role as played by the components ${\alpha }_j,j=1,\cdots ,TP$ , of the “vector of primary model parameters” $\alpha \triangleq {\left({\alpha }_1,\cdots ,{\alpha }_{TP}\right)}^{†}$ within the framework of the 1^st-CASAM-N. It is important to underscore at the outset that the total number of model parameters is by definition much larger than the total number of components of the feature function $f\left(\alpha \right)$ , i.e., $TP\gg TF$ ; TF coincides with TP if the mathematical/computational model under consideration possesses no “features.”

An examination of the pattern underlying the n^th-CASAM-N methodology [31] indicates that the n^th-order sensitivity of the model’s response $R\left[u\left(x\right);f\left(\alpha \right)\right]$ with respect to the components $f_{j_1},\cdots ,f_{j_n}$ of the “feature” function $f\left(\alpha \right)\triangleq {\left[f_{j_1}\left(\alpha \right),\cdots ,f_{j_n}\left(\alpha \right)\right]}^{†}$ , for each $j_1,\cdots ,j_n=1,\cdots ,TF$ is expected to have the functional form $R^{\left(n\right)}\left[j_n;\cdots ;j_1;U^{\left(n\right)};A^{\left(n\right)};f\left(\alpha \right)\right]\triangleq {\partial }^{n}R\left[u\left(x\right);f\left(\alpha \right)\right]/\partial f_{j_n}\cdots \partial f_{j_1}$ . It can also be surmised that the n^th-order sensitivity $R^{\left(n\right)}\left[j_n;\cdots ;j_1;U^{\left(n\right)};A^{(n)};f\left(\alpha \right)\right]$ stems from the total first-order G-differential of each of the (n-1)^th-order sensitivities, each of which is expected to have the following expression for $j_1,\cdots ,j_{n-1}=1,\cdots ,TF$ :

$\begin{array}{l}{R}^{\left(n-1\right)}\left[j_{n-1};\cdots ;j_1;U^{\left(n-1\right)};A^{\left(n-1\right)};f\left(\alpha \right)\right]\equiv {\partial }^{n-1}R\left[u\left(x\right);f\left(\alpha \right)\right]/\partial {\alpha }_{j_{n-1}}\cdots \partial {\alpha }_{j_1}\\ \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(n-1\right)}\left[j_{n-1};\cdots ;j_1;U^{\left(n-1\right)};A^{\left(n-1\right)};f\left(\alpha \right)\right]\text{d}{x}_1\cdots \text{d}{x}_{TI}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}}}.\end{array}$ (1)

Analogous with the pattern underlying the n^th-CASAM-Nmethodology [31] , it is surmised that the total first-order G-differential of the (n-1)^th-order sensitivity $R^{\left(n-1\right)}\left[j_{n-1};\cdots ;j_1;U^{\left(n-1\right)};A^{\left(n-1\right)};f\left(\alpha \right)\right]$ with respect to the feature-function $f\left(\alpha \right)$ is expected to have the following expression:

$\begin{array}{l}{\left\{\delta R^{\left(n-1\right)}\left[j_{n-1};\cdots ;j_1;U^{\left(n-1\right)};A^{(n-1)};f\right]\right\}}_{{\alpha }^0}\\ \triangleq {\left\{\frac{\text{d}}{\text{d}\epsilon }\delta R^{\left(n-1\right)}\left[j_{n-1};\cdots ;j_1;U^{\left(n-1\right)}+\epsilon V^{\left(n-1\right)};A^{\left(n-1\right)}+\epsilon \delta A^{\left(n-1\right)};f+\epsilon \delta f\right]\right\}}_{\epsilon =0}\\ ={\displaystyle \underset{j_n=1}{\overset{TF}{\sum }}{\left\{\frac{\partial R^{\left(n-1\right)}\left[j_{n-1};\cdots ;j_1;U^{\left(n-1\right)};A^{\left(n-1\right)};f\right]}{\partial f_{j_n}}\right\}}_{{\alpha }^0}}\delta f_{j_n}\\ +{\left\{\delta R^{\left(n-1\right)}\left[j_{n-1};\cdots ;j_1;U^{\left(n-1\right)};A^{\left(n-1\right)};f;V^{\left(n-1\right)};\delta A^{\left(n-1\right)}\right]\right\}}_{ind},\end{array}$ (2)

where the quantity ${\left\{\delta R^{\left(n-1\right)}\left[j_{n-1};\cdots ;j_1;U^{\left(n-1\right)};A^{\left(n-1\right)};f;V^{\left(n-1\right)};\delta A^{\left(n-1\right)}\right]\right\}}_{ind}$ denotes

the so-called “indirect-effect term” which is defined as follows:

$\begin{array}{l}{\left\{\delta R^{\left(n-1\right)}\left[j_{n-1};\cdots ;j_1;U^{\left(n-1\right)};A^{\left(n-1\right)};f;V^{\left(n-1\right)};\delta A^{\left(n-1\right)}\right]\right\}}_{ind}\\ \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 }}{\left\{\frac{\partial S^{\left(n-1\right)}}{\partial U^{\left(n-1\right)}\left(x\right)}{V}^{\left(n-1\right)}\left(2^{n-2};x\right)+\frac{\partial S^{\left(n-1\right)}}{\partial A^{(n-1)}\left(x\right)}\delta A^{\left(n-1\right)}\left(2^{n-2};x\right)\right\}}_{{\alpha }^0}\text{d}{x}_1\cdots \text{d}{x}_{TI}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}}}.\end{array}$

(3)

The vectors $V^{\left(n-1\right)}\left(2^{n-2};j_{n-3};\cdots ;j_1;x\right)\triangleq \delta U^{\left(n-1\right)}\left(2^{n-2};j_{n-3};\cdots ;j_1;x\right)$ and $\delta A^{\left(n-1\right)}\left(2^{n-2};j_{n-2};\cdots ;j_1;x\right)$ are the solution of the following n^th-Level Variational Sensitivity System (n^th-LVSS), which is obtained by concatenating the (n-1)^th-LVSS together with the G-differentiated (n-1)^th-Level Adjoint Sensitivity System, for $j_1,\cdots ,j_{n-1}=1,\cdots ,TF$ :

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

${\left\{B_V^{\left(n\right)}\left[2^{n-1};U^{\left(n\right)}\left(2^{n-1};x\right);V^{\left(n\right)}\left(2^{n-1};x\right);f;\delta f\right]\right\}}_{{\alpha }^0}=0\left[2^{n-1}\right];x\in \partial {\Omega }_x\left({\alpha }^0\right),$ (5)

The various matrices and vectors appearing in Equations (4) and (5) are defined as follows:

1) The variational matrix $V{M}^{\left(n\right)}\left[2^{n-1}\times 2^{n-1};U^{\left(n\right)}\left(2^{n-1};j_{n-2};\cdots ;j_1;x\right);f\right]$ comprises $\left(2^{n-1}\times 2^{n-1}\right)$ block-matrices, each comprising TD² components/elements, defined as follows:

$V{M}^{\left(n\right)}\left[2^{n-1}\times 2^{n-1};x;f\right]\triangleq \left(\begin{array}{cc}V{M}^{\left(n-1\right)}\left[2^{n-2}\times 2^{n-2};x\right]& 0\left[2^{n-2}\times 2^{n-2}\right]\\ V{M}_{21}^{\left(n\right)}\left(2^{n-2}\times 2^{n-2};x\right)& V{M}_{22}^{\left(n\right)}\left(2^{n-2}\times 2^{n-2};x\right)\end{array}\right).$ (6)

The matrix $V{M}^{\left(n\right)}\left[2^{n-1}\times 2^{n-1};U^{\left(n\right)}\left(2^{n-1};j_{n-2};\cdots ;j_1;x\right);f\right]$ comprises a total of $\left(2^{n-1}\times 2^{n-1}\right)T{D}^2$ components/elements, as indicated by its first argument. The submatrix $V{M}^{\left(n-1\right)}\left[2^{n-2}\times 2^{n-2};x\right]$ is defined recursively, while the remaining submatrices are defined as follows:

$\begin{array}{l}V{M}_{21}^{\left(n\right)}\left(2^{n-2}\times 2^{n-2};j_{n-2};\cdots ;j_1;x\right)\\ \triangleq -\frac{\partial Q_A^{\left(n-1\right)}\left[2^{n-2};j_{n-2};\cdots ;j_1;U^{\left(n-2\right)}\left(2^{n-2};j_{n-3};\cdots ;j_1;x\right);f\right]}{\partial U^{\left(n-1\right)}\left(2^{n-2};j_{n-3};\cdots ;j_1;x\right)}\\ +\frac{\partial \left\{A{M}^{\left(n-1\right)}\left[2^{n-2}\times 2^{n-2};U^{\left(n-1\right)}\left(2^{n-2};x\right);f\right]A^{(n-1)}\left(2^{n-2};j_{n-2};\cdots ;j_1;x\right)\right\}}{\partial U^{\left(n-1\right)}\left(2^{n-2};j_{n-3};\cdots ;j_1;x\right)};\end{array}$ (7)

$V{M}_{22}^{\left(n\right)}\left(2^{n-2}\times 2^{n-2};x\right)\triangleq A{M}^{\left(n-1\right)}\left[2^{n-2}\times 2^{n-2};U^{\left(n-1\right)}\left(2^{n-2};j_{n-3};\cdots ;j_1;x\right);f\right]$ (8)

2) The vector $V^{\left(n\right)}\left(2^{n-1};j_{n-2};\cdots ;j_1;x\right)$ is defined recursively, below, and comprises $2^{n-1}$ blocks of TD-dimensional vectors [thus comprising a total of $\left(2^{n-1}\times TD\right)$ components/elements]; this fact is indicated by the first argument of each of these quantities.

$\begin{array}{l}{V}^{\left(n\right)}\left(2^{n-1};j_{n-2};\cdots ;j_1;x\right)\\ \triangleq \delta U^{\left(n\right)}\left(2^{n-1};j_{n-2};\cdots ;j_1;x\right)=\left(\begin{array}{c}{V}^{\left(n-1\right)}\left(2^{n-2};j_{n-3};\cdots ;j_1;x\right)\\ \delta A^{\left(n-1\right)}\left(2^{n-2};j_{n-2};\cdots ;j_1;x\right)\end{array}\right)\\ =[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),\cdots ,\\ {\delta a^{\left(n-1\right)}\left(1;j_{n-2};\cdots ;j_1;x\right),\cdots ,\delta a^{\left(n-1\right)}\left(2^{n-2};j_{n-2};\cdots ;j_1;x\right)]}^{†};\end{array}$ (9)

3) The vector $U^{\left(n\right)}\left(2^{n-1};j_{n-2};\cdots ;j_1;x\right)$ has the same structure as the vector $V^{\left(n\right)}\left(2^{n-1};j_{n-2};\cdots ;j_1;x\right)$ , comprising $2^{n-1}$ blocks of TD-dimensional vectors, as defined below.2

$U^{\left(n\right)}\left(2^{n-1};j_{n-2};\cdots ;j_1;x\right)\triangleq \left(\begin{array}{c}{U}^{\left(n-1\right)}\left(2^{n-2};j_{n-3};\cdots ;j_1;x\right)\\ A^{\left(n-1\right)}\left(2^{n-2};j_{n-2};\cdots ;j_1;x\right)\end{array}\right);$ (10)

4) The vector $Q_V^{\left(n\right)}\left[2^{n-1};U^{\left(n\right)}\left(2^{n-1};j_{n-2};\cdots ;j_1;x\right);f;\delta f\right]$ also comprises $2^{n-1}$ blocks of TD-dimensional vectors, as defined recursively below.

$\begin{array}{l}{Q}_V^{\left(n\right)}\left[2^{n-1};j_1;U^{\left(n\right)}\left(2^{n-1};j_{n-2};\cdots ;j_1;x\right);f;\delta f\right]\triangleq \left(\begin{array}{c}{Q}_V^{\left(n-2\right)}\left[2^{n-2};U^{\left(n-1\right)}\left(2^{n-2};x\right);f;\delta f\right]\\ Q_2^{\left(n\right)}\left[2^{n-2};U^{\left(n\right)}\left(2^{n-1};x\right);f;\delta f\right]\end{array}\right)\\ \triangleq {\left\{q_V^{\left(n\right)}\left[1;U^{\left(n\right)}\left(2^{n-1};j_{n-2};\cdots ;j_1;x\right);f;\delta f\right],\cdots ,q_V^{\left(n\right)}\left[2^{n-1};U^{\left(n\right)}\left(2^{n-1};j_{n-2};\cdots ;j_1;x\right);f;\delta f\right]\right\}}^{†};\end{array}$

(11)

$\begin{array}{l}{q}_V^{\left(n\right)}\left[i;U^{\left(n\right)}\left(x\right);f;\delta f\right]\equiv {\displaystyle \underset{j_n=1}{\overset{TP}{\sum }}{s}_V^{\left(n\right)}\left[i;j_{n-2};\cdots ;j_1;U^{\left(n\right)}\left(x\right);f\right]\delta f_{j_n}};\\ i=1,\cdots ;2^{n-1};\end{array}$ (12)

$\begin{array}{l}{Q}_2^{\left(n\right)}\left[2^{n-2};U^{\left(n\right)}\left(2^{n-1};j_{n-2};\cdots ;j_1;x\right);f;\delta f\right]\\ \triangleq \frac{\partial Q_A^{\left(n-1\right)}\left[2^{n-2};j_1;U^{\left(n-1\right)}\left(x\right);f\right]}{\partial f}\partial f\\ -\frac{\partial \left\{A{M}^{\left(n-1\right)}\left[2^{n-2}\times 2^{n-2};U^{\left(n-1\right)}\left(2^{n-2};x\right);f\right]A^{\left(n-1\right)}\left(2^{n-2};x\right)\right\}}{\partial f}\partial f;\end{array}$ (13)

5) The vector block-vector $0\left[2^{n-1}\right]$ comprises $2^{n-1}$ components, each component being a TD-dimensional vector having identically zero components.

6) The boundary terms are represented by the block-vector $B_V^{\left(n\right)}\left[2^{n-1};U^{\left(n\right)}\left(2^{n-1};j_{n-2};\cdots ;j_1;x\right);V^{\left(n\right)}\left(2^{n-1};j_{n-2};\cdots ;j_1;x\right);f;\delta f\right]$ , which is defined recursively as shown below.

$\begin{array}{l}{B}_V^{\left(n\right)}\left[2^{n-1};U^{\left(n\right)}\left(2^{n-1};j_{n-2};\cdots ;j_1;x\right);V^{\left(n\right)}\left(2^{n-1};j_{n-2};\cdots ;j_1;x\right);f;\delta f\right]\\ \triangleq \left(\begin{array}{c}{B}_V^{\left(n-1\right)}\left[2^{n-2};U^{\left(n-1\right)}\left(2^{n-2};j_{n-3};\cdots ;j_1;x\right);V^{\left(n-1\right)}\left(2^{n-2};j_{n-3};\cdots ;j_1;x\right);f;\delta f\right]\\ \delta B_A^{\left(n-1\right)}\left[2^{n-2};U^{\left(n\right)}\left(2^{n-1};j_{n-2};\cdots ;j_1;x\right);V^{\left(n\right)}\left(2^{n-1};j_{n-2};\cdots ;j_1;x\right);f;\delta f\right]\end{array}\right).\end{array}$ (14)

The need for solving the following n^th-Level Variational Sensitivity System (n^th-LVSS) defined by Equations (4) and (5) is circumvented by deriving an alternative expression for the indirect-effect term defined in Equation (3), in which the function $V^{\left(n\right)}\left(2^{n-1};j_{n-2};\cdots ;j_1;x\right)$ is replaced by a n^th-level adjoint function which will be independent of parameter variations, and which will be denoted as $A^{\left(n\right)}\left(2^{n-1};j_{n-1};\cdots ;j_1;x\right)\triangleq {\left[a^{\left(n\right)}\left(1;j_{n-1};\cdots ;j_1;x\right),\cdots ,a^{\left(n\right)}\left(2^{n-1};j_{n-1};\cdots ;j_1;x\right)\right]}^{†}\in H_n\left({\Omega }_x\right)$ . The elements of the Hilbert space $H_n\left({\Omega }_x\right)$ are surmised to be block-vectors of the form ${\Psi }^{\left(n\right)}\left(2^{n-1};j_{n-1};\cdots ;j_1;x\right)\triangleq {\left[{\psi }^{\left(n\right)}\left(1;j_{n-1};\cdots ;j_1;x\right),\cdots ,{\psi }^{(n)}\left(2^{n-1};j_{n-1};\cdots ;j_1;x\right)\right]}^{†}$ , comprising as elements $2^{n-1}$ TD-dimensional vectors of the form ${\psi }^{\left(n\right)}\left(i;x\right)\triangleq {\left[{\psi }_1^{\left(n\right)}\left(i;x\right),\cdots ,{\psi }_{TD}^{\left(n\right)}\left(i;x\right)\right]}^{†}\in H_1\left({\Omega }_x\right)$ , $i=1,\cdots ,2^n$ . The inner product

of two vectors ${\Psi }^{\left(n\right)}\left(2^{n-1};j_{n-1};\cdots ;j_1;x\right)$ and ${\Phi }^{\left(n\right)}\left(2^{n-1};j_{n-1};\cdots ;j_1;x\right)$ in the

Hilbert space $H_n\left({\Omega }_x\right)$ is denoted as ${\langle {\Psi }^{\left(n\right)}\left(2^{n-1};x\right),{\Phi }^{\left(n\right)}\left(2^{n-1};x\right)\rangle }_n$ and defined as follows:

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

The n^th-Level Adjoint Sensitivity System (n^th-LASS) for the n^th-level adjoint function $A^{\left(n\right)}\left(2^{n-1};j_{n-1};\cdots ;j_1;x\right)$ is obtained by using the inner product defined in Equation (15), as follows:

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

where the quantity ${\left\{{\left[P^{\left(n\right)}\left(U^{\left(n\right)};A^{\left(n\right)};V^{\left(n\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.

In terms of the n^th-Level adjoint function $A^{\left(n\right)}\left(2^{n-1};j_{n-1};\cdots ;j_1;x\right)$ , the indirect-effect term defined by Equation (3) will have the following expression:

$\begin{array}{l}{\left\{\delta R^{\left(n-1\right)}\left[j_{n-1};\cdots ;j_1;U^{\left(n-1\right)};A^{\left(n-1\right)};\alpha ;V^{\left(n-1\right)};\delta A^{\left(n-1\right)}\right]\right\}}_{ind}\\ =-{\left\{{\left[{\hat{P}}^{\left(n\right)}\left(U^{\left(n\right)};A^{\left(n\right)};\delta \alpha \right)\right]}_{\partial {\Omega }_x}\right\}}_{{\alpha }^0}\\ +{\left\{{\langle A^{\left(n\right)}\left(2^{n-1};j_{n-1};\cdots ;j_1;x\right),Q_V^{\left(n\right)}\left[2^{n-1};U^{\left(n\right)}\left(2^{n-1};j_{n-2};\cdots ;j_1;x\right);\alpha ;\delta \alpha \right]\rangle }_n\right\}}_{{\alpha }^0}.\end{array}$ (17)

where $A^{\left(n\right)}\left(2^{n-1};j_{n-1};\cdots ;j_1;x\right)$ is the solution of the following n^th-Level Adjoint Sensitivity System (n^th-LASS):

$\begin{array}{l}A{M}^{\left(n\right)}\left[2^{n-1}\times 2^{n-1};U^{\left(n\right)}\left(2^{n-1};x\right);\alpha \right]A^{\left(n\right)}\left(2^{n-1};j_{n-1};\cdots ;j_1;x\right)\\ =Q_A^{\left(n\right)}\left[2^{n-1};j_{n-1};\cdots ;j_1;U^{\left(n\right)}\left(2^{n-1};j_{n-2};\cdots ;j_1;x\right);\alpha \right],\end{array}$ (18)

subject to boundary conditions represented in operator form as follows:

$\begin{array}{l}{\left\{B_A^{\left(n\right)}\left[2^{n-1};U^{\left(n\right)}\left(2^{n-1};j_{n-2};\cdots ;j_1;x\right);A^{\left(n\right)}\left(2^{n-1};j_{n-1};\cdots ;j_1;x\right);\alpha \right]\right\}}_{{\alpha }^0}=0\left[2^{n-1}\right],\\ x\in \partial {\Omega }_x\left({\alpha }^0\right);j_1=1,\cdots ,TP;j_2=1,\cdots ,j_1;\cdots ;j_{n-1}=1,\cdots ,j_{n-2}.\end{array}$ (19)

In Equation (17), the quantity ${\left\{{\left[{\hat{P}}^{\left(n\right)}\left(U^{\left(n\right)};A^{\left(n\right)};\delta f\right)\right]}_{\partial {\Omega }_x}\right\}}_{{\alpha }^0}$ denotes residual

boundary terms which may have not vanished automatically after having used the boundary conditions provided in Equations (5) and (19) to eliminate from Equation (17) all unknown values of the nth-level variational function $V^{\left(n\right)}\left(2^{n-1};j_{n-2};\cdots ;j_1;x\right)$ .

The quantities which appear in the definition of the n^th-LASS represented by Equations (18) and (19), are defined as follows:

$\begin{array}{l}A{M}^{\left(n\right)}\left[2^{n-1}\times 2^{n-1};U^{\left(n\right)}\left(2^{n-1};j_{n-2};\cdots ;j_1;x\right);f\right]\\ \triangleq {\left[V{M}^{\left(n\right)}\left(2^{n-1}\times 2^{n-1};U^{\left(n\right)};f\right)\right]}^{*}\\ =\left(\begin{array}{cc}{\left\{{\left[V{M}^{\left(n-1\right)}\left(2^{n-2}\times 2^{n-2}\right)\right]}^{*}\right\}}^{†}& {\left\{{\left[V{M}_{21}^{\left(n\right)}\left(2^{n-2}\times 2^{n-2}\right)\right]}^{*}\right\}}^{†}\\ 0\left[2^{n-2}\times 2^{n-2}\right]& {\left\{{\left[V{M}_{22}^{\left(n\right)}\left(2^{n-2}\times 2^{n-2}\right)\right]}^{*}\right\}}^{†}\end{array}\right),\end{array}$ (20)

$\begin{array}{l}{Q}_A^{\left(n\right)}\left[2^{n-1};j_{n-1};\cdots ;j_1;U^{\left(n\right)}\left(2^{n-1};x\right);f\right]\triangleq \{q_A^{\left(n\right)}\left[1;j_{n-1};\mathrm{...};j_1;U^{\left(n\right)}\left(2^{n-1};x\right);f\right],\\ \cdots ,{q_A^{\left(n\right)}\left[2^{n-1};j_{n-1};\cdots ;j_1;U^{\left(n\right)}\left(2^{n-1};x\right);f\right]\}}^{†},j_1=1,\cdots ,TF;\cdots ;j_{n-1}=1,\cdots ,j_{n-2};\end{array}$ (21)

$q_A^{\left(n\right)}\left[1;j_{n-1};\cdots ;j_1;U^{\left(n\right)}\left(x\right);f\right]\triangleq \partial S^{\left(n-1\right)}\left[j_{n-1};\cdots ;j_1;U^{\left(n\right)};f\right]/\partial u\left(x\right);$ (22)

$q_A^{\left(n\right)}\left[2;j_{n-1};\cdots ;j_1;U^{\left(n\right)}\left(x\right);f\right]\triangleq \partial S^{\left(n-1\right)}\left[j_{n-1};\cdots ;j_1;U^{\left(n\right)};f\right]/\partial a^{\left(1\right)}\left(x\right);$ (23)

$\begin{array}{l}\text{for}n\ge 3:q_A^{\left(n\right)}\left[2^k+i;j_{n-1};\cdots ;j_1;U^{\left(n\right)}\left(x\right);f\right]\\ \triangleq \frac{\partial S^{\left(n-1\right)}\left[j_{n-1};\cdots ;j_1;U^{\left(n\right)};f\right]}{\partial a^{\left(k+1\right)}\left(i;j_k;\cdots ;j_1;x\right)};k=1,\cdots ,n-2;i=1,\cdots ,2^k.\end{array}$ (24)

The final expression of the total differential expressed by Equation (17) is obtained by inserting in Equation (2) the expression for the indirect-effect term obtained in Equation (17), which yields the following expression for the n^th-Order sensitivities $R^{\left(n\right)}\left[j_n;\cdots ;j_1;U^{\left(n\right)};A^{(n)};f\right]\equiv {\partial }^{n}R\left[u\left(x\right);f\right]/\partial f_{j_n}\cdots \partial f_{j_1}$ of the response $R\left[u\left(x\right);\alpha \right]$ with respect to the parameters ${\alpha }_{j_1},\cdots ,{\alpha }_{j_n}$ , for $j_1=1,\cdots ,TF;\cdots ;j_{n-1}=1,\cdots ,j_{n-2}$ :

$\begin{array}{l}{R}^{\left(n\right)}\left[j_n;\cdots ;j_1;U^{\left(n\right)};A^{(n)};f\right]\equiv {\partial }^{n}R\left[u\left(x\right);f\right]/\partial f_{j_n}\cdots \partial f_{j_1}\\ =\frac{\partial R^{\left(n-1\right)}\left[j_{n-1};\cdots ;j_1;U^{\left(n-1\right)}\left(x\right);A^{\left(n-1\right)}\left(x\right);f\right]}{\partial f_{j_n}}-\frac{{\left[\partial {\hat{P}}^{\left(n\right)}\left(U^{\left(n\right)};A^{\left(n\right)};\delta f\right)\right]}_{\partial {\Omega }_x}}{\partial f_{j_n}}\\ +{\displaystyle \underset{i=1}{\overset{2^{n-1}}{\sum }}{\langle a^{\left(n\right)}\left(i;j_{n-1};\cdots ;j_1;x\right),s_V^{\left(n\right)}\left[i;j_n;\cdots ;j_1;U^{\left(n\right)}\left(x\right);f\right]\rangle }_1}.\\ \triangleq {\displaystyle \underset{{\lambda }_1\left(\alpha \right)}{\overset{{\omega }_1\left(\alpha \right)}{\int }}\mathrm{...}{\displaystyle \underset{{\lambda }_{TI}\left(\alpha \right)}{\overset{{\omega }_{TI}\left(\alpha \right)}{\int }}{S}^{\left(n\right)}\left[j_n;\cdots ;j_1;U^{\left(n\right)}\left(2^{n-1};x\right);A^{\left(n\right)}\left(2^{n-1};x\right);f\right]\text{d}{x}_1\cdots \text{d}{x}_{TI}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}}}.\end{array}$ (25)

2.2. The Particular Form of the n^th-FASAM-N Framework for n = 1

Setting n = 1 into the mathematical framework of the n^th-FASAM-N conjectured in Section 2.1 above yields the following expressions for the particular case n = 1:

i) Expression of the model response:

$\begin{array}{l}{R}^{\left(0\right)}\left[U^{\left(0\right)};f\right]\equiv R\left[u\left(x\right);f\left(\alpha \right)\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);g\left(\alpha \right);h\left(\alpha \right);x\right]\text{d}{x}_1\cdots \text{d}{x}_{TI}}}.\end{array}$ (26)

ii) Expression of the 1^st-order response sensitivities to the components of the feature function of model parameters, for $j_1=1,\cdots ,TF$ :

$\begin{array}{l}{R}^{\left(1\right)}\left[j_1;U^{\left(1\right)}\left(2^0;x\right);A^{\left(1\right)}\left(2^0;x\right);f\left(\alpha \right)\right]=\frac{\partial R\left[u\left(x\right);f\left(\alpha \right)\right]}{\partial f_{j_1}}\\ =-\frac{{\left[\partial {\hat{P}}^{\left(1\right)}\left(U^{\left(1\right)};A^{\left(1\right)};\delta f\right)\right]}_{\partial {\Omega }_x}}{\partial f_{j_1}}+{\langle a^{\left(1\right)}\left(x\right),s_V^{\left(1\right)}\left[U^{\left(1\right)}\left(x\right);\alpha \right]\rangle }_1.\end{array}$ (27)

The various quantities which appear in Equation (27) take on the following particular forms for n = 1:

$\begin{array}{l}{U}^{\left(1\right)}\left(2^0;x\right)\equiv u\left(x\right);A^{\left(1\right)}\left(2^0;x\right)\triangleq a^{\left(1\right)}\left(x\right);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 f_{j_1}}\end{array}$ (28)

The 1^st-level adjoint sensitivity function $A^{\left(1\right)}\left(2^0;x\right)\equiv a^{\left(1\right)}\left(x\right)$ is the solution of the following 1^st-LASS obtained by setting n=1 in Equations (18) and (19):

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

$\begin{array}{l}{\left\{B_A^{\left(1\right)}\left[2^0;U^{\left(1\right)}\left(2^0;x\right);V^{\left(1\right)}\left(2^0;x\right);\alpha ;\delta f\right]\right\}}_{{\alpha }^0}=0,\\ x\in \partial {\Omega }_x\left({\alpha }^0\right).\end{array}$ (30)

where:

$A{M}^{\left(1\right)}\left[2^0\times 2^0;U^{\left(1\right)}\left(2^0;x\right);f\right]\equiv {\left[N^{\left(1\right)}\left(u;f\right)\right]}^{*};$ (31)

$Q_A^{\left(1\right)}\left[2^0;U^{\left(1\right)}\left(2^0;x\right);\alpha ;\delta f\right]\equiv q_A^{\left(1\right)}\left[u\left(x\right);f\right]\triangleq {\left\{\partial S\left(u;f\right)/\partial u\right\}}_{{\alpha }^0};$ (32)

$B_V^{\left(1\right)}\left[2^0;U^{\left(1\right)}\left(2^0;x\right);V^{\left(1\right)}\left(2^0;x\right);f;\delta f\right]\equiv b_A^{\left(1\right)}\left(u;a^{\left(1\right)};f\right).$ (33)

Comparing the expressions obtained above in Equations (26)-(33) with the expressions obtained independently, from first principles, in Appendix B reveals that the corresponding expressions are identical to each other, which proves the correctness of the conjectured general expressions obtained within the n^th-FASAM-N methodology for the particular case n = 1.

2.3. Derivation of the Mathematical Framework of the (n+1)^th-FASAM-N

It will be assumed that, for each $j_1,\cdots ,j_n=1,\cdots ,TF$ , the 1^st-order total G-differential of the n^th-order sensitivities $R^{\left(n\right)}\left[j_n;\cdots ;j_1;U^{\left(n\right)};A^{\left(n\right)};f\right]$ will exist and will be linear in the variations

$V^{\left(n\right)}\left(2^{n-1};j_{n-2};\cdots ;j_1;x\right)\triangleq \delta U^{\left(n-1\right)}\left(2^{n-1};j_{n-2};\cdots ;j_1;x\right)$ and

$\delta A^{\left(n\right)}\left(2^{n-1};j_{n-1};\cdots ;j_1;x\right)$ in a neighborhood around the nominal values of the “feature functions,” the parameters, and the respective state functions. By definition, the 1^st-order total G-differential of $R^{\left(n\right)}\left[j_n;\cdots ;j_1;U^{\left(n\right)};A^{\left(n\right)};f\right]$ is given by the following expression:

$\begin{array}{l}{\left\{\delta R^{\left(n\right)}\left[j_n;\cdots ;j_1;U^{\left(n\right)};A^{\left(n\right)};f\right]\right\}}_{{\alpha }^0}\\ \triangleq \frac{\text{d}}{\text{d}\epsilon }\{R^{\left(n\right)}[j_n;\cdots ;j_1;U^{\left(n\right)}+\epsilon V^{\left(n\right)};{A^{\left(n\right)}+\epsilon \delta A^{\left(n\right)};f+\epsilon \delta f]\}}_{\epsilon =0}\\ \triangleq {\displaystyle \underset{j_{n+1}=1}{\overset{TP}{\sum }}{\left\{\frac{\partial R^{\left(n\right)}\left[\cdots ;U^{\left(n\right)};A^{\left(n\right)};f\right]}{\partial f_{j_{n+1}}}\right\}}_{{\alpha }^0}}\delta f_{j_{n+1}}\\ +{\left\{\delta R^{\left(n\right)}\left[j_n;\cdots ;j_1;U^{\left(n\right)};A^{\left(n\right)};f;V^{\left(n\right)}\left(x\right);\delta A^{\left(n\right)}\left(x\right)\right]\right\}}_{ind},\end{array}$ (34)

where the quantity ${\left\{\delta R^{\left(n\right)}\left[j_n;\cdots ;j_1;U^{\left(n\right)};A^{\left(n\right)};f;V^{\left(n\right)};\delta A^{\left(n\right)}\right]\right\}}_{ind}$ denotes the “indirect-effect term” which is defined as follows:

$\begin{array}{l}{\left\{\delta R^{\left(n\right)}\left[j_n;\cdots ;j_1;U^{\left(n\right)};A^{\left(n\right)};f;V^{\left(n\right)};\delta A^{\left(n\right)}\right]\right\}}_{ind}\\ \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 }}{\left\{\frac{\partial S^{\left(n\right)}}{\partial U^{\left(n\right)}\left(x\right)}{V}^{\left(n\right)}\left(2^{n-1};x\right)+\frac{\partial S^{\left(n\right)}}{\partial A^{\left(n\right)}\left(x\right)}\delta A^{\left(n\right)}\left(2^{n-1};x\right)\right\}}_{{\alpha }^0}\text{d}{x}_1\cdots \text{d}{x}_{TI}\text{}}}\end{array}$ (35)

The vectors $V^{\left(n\right)}\left(2^{n-1};j_{n-2};\cdots ;j_1;x\right)\triangleq \delta U^{\left(n\right)}\left(2^{n-1};j_{n-2};\cdots ;j_1;x\right)$ and

$\delta A^{\left(n\right)}\left(2^{n-1};j_{n-1};\cdots ;j_1;x\right)$ , which are needed in order to evaluate the indirect-

effect term ${\left\{\delta R^{\left(n\right)}\left[j_n;\cdots ;j_1;U^{\left(n\right)};A^{\left(n\right)};f;V^{\left(n\right)};\delta A^{\left(n\right)}\right]\right\}}_{ind}$ , are the solutions of the

following (n+1)^th-Level Variational Sensitivity System, which is obtained by concatenating the nth-LVSS defined by Equations (4) and (5) together with the G-differentiated nth-LASS, for $j_1,\cdots ,j_n=1,\cdots ,TP$ :

$\begin{array}{l}{\left\{V{M}^{\left(n+1\right)}\left[2^n\times 2^n;U^{\left(n+1\right)}\left(2^n;j_{n-2};\cdots ;j_1;x\right);f\right]V^{\left(n+1\right)}\left(2^n;j_{n-2};\cdots ;j_1;x\right)\right\}}_{{\alpha }^0}\\ ={\left\{Q_V^{\left(n+1\right)}\left[2^n;U^{\left(n\right)}\left(2^n;j_{n-2};\cdots ;j_1;x\right);f;\delta f\right]\right\}}_{{\alpha }^0},x\in {\Omega }_x,\end{array}$ (36)

${\left\{B_V^{\left(n+1\right)}\left[2^n;U^{\left(n+1\right)}\left(2^n;x\right);V^{\left(n+1\right)}\left(2^n;x\right);f;\delta f\right]\right\}}_{{\alpha }^0}=0\left[2^n\right];x\in \partial {\Omega }_x\left({\alpha }^0\right)$ (37)

where:

$U^{\left(n+1\right)}\left(2^n;j_{n-1};\cdots ;j_1;x\right)\triangleq \left(\begin{array}{c}{U}^{\left(n\right)}\left(2^{n-1};j_{n-2};\cdots ;j_1;x\right)\\ A^{(n)}\left(2^{n-1};j_{n-1};\cdots ;j_1;x\right)\end{array}\right);$ (38)

$\begin{array}{l}{V}^{\left(n+1\right)}\left(2^n;j_{n-1};\cdots ;j_1;x\right)\\ \triangleq \delta U^{\left(n+1\right)}\left(2^n;j_{n-1};\cdots ;j_1;x\right)=\left(\begin{array}{c}{V}^{\left(n\right)}\left(2^{n-1};j_{n-2};\cdots ;j_1;x\right)\\ \delta A^{(n)}\left(2^{n-1};j_{n-1};\cdots ;j_1;x\right)\end{array}\right)\\ =[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),\cdots ,\\ {\delta a^{\left(n\right)}\left(1;j_{n-1};\cdots ;j_1;x\right),\cdots ,\delta a^{(n)}\left(2^{n-1};j_{n-1};\cdots ;j_1;x\right)]}^{†};\end{array}$ (39)

$\begin{array}{l}V{M}^{\left(n+1\right)}\left[2^n\times 2^n;U^{\left(n+1\right)}\left(2^n;j_{n-1};\cdots ;j_1;x\right);f\right]\\ \triangleq \left(\begin{array}{cc}V{M}^{\left(n\right)}\left[2^{n-1}\times 2^{n-1};x;f\right]& 0\left[2^{n-1}\times 2^{n-1}\right]\\ V{M}_{21}^{\left(n+1\right)}\left(2^{n-1}\times 2^{n-1};x;f\right)& V{M}_{22}^{\left(n+1\right)}\left(2^{n-1}\times 2^{n-1};x;f\right)\end{array}\right);\end{array}$ (40)

$\begin{array}{l}V{M}_{21}^{\left(n+1\right)}\left(2^{n-1}\times 2^{n-1};j_{n-1};\mathrm{...};j_1;x;f\right)\\ \triangleq -\frac{\partial Q_A^{\left(n\right)}\left[2^{n-1};j_{n-1};\cdots ;j_1;U^{\left(n\right)}\left(2^{n-1};j_{n-2};\cdots ;j_1;x\right);f\right]}{\partial U^{\left(n\right)}\left(2^{n-1};j_{n-2};\cdots ;j_1;x\right)}\\ +\frac{\partial \left\{A{M}^{\left(n\right)}\left[2^{n-1}\times 2^{n-1};U^{\left(n\right)}\left(2^{n-1};x\right);f\right]A^{\left(n\right)}\left(2^{n-1};j_{n-1};\cdots ;j_1;x\right)\right\}}{\partial U^{\left(n\right)}\left(2^{n-1};j_{n-2};\cdots ;j_1;x\right)};\end{array}$ (41)

$V{M}_{22}^{\left(n+1\right)}\left(2^{n-1}\times 2^{n-1};x;f\right)\triangleq A{M}^{\left(n\right)}\left[2^{n-1}\times 2^{n-1};U^{\left(n\right)}\left(2^{n-1};j_{n-2};\cdots ;j_1;x\right);f\right];$ (42)

$\begin{array}{l}{Q}_V^{\left(n+1\right)}\left[2^n;j_1;U^{\left(n+1\right)}\left(2^n;j_{n-1};\cdots ;j_1;x\right);f;\delta f\right]\\ \triangleq \left(\begin{array}{c}{Q}_V^{\left(n\right)}\left[2^{n-1};U^{\left(n\right)}\left(2^{n-1};x\right);f;\delta f\right]\\ Q_2^{\left(n+1\right)}\left[2^{n-1};U^{\left(n+1\right)}\left(2^n;x\right);f;\delta f\right]\end{array}\right)\\ \triangleq \{q_V^{\left(n+1\right)}\left[1;U^{\left(n+1\right)}\left(2^n;j_{n-1};\cdots ;j_1;x\right);f;\delta f\right],\\ {\cdots ,q_V^{\left(n+1\right)}\left[2^n;U^{\left(n+1\right)}\left(2^n;j_{n-1};\cdots ;j_1;x\right);f;\delta f\right]\}}^{†};\end{array}$ (43)

$\begin{array}{l}{q}_V^{\left(n\right)}\left[i;U^{\left(n\right)}\left(x\right);\alpha ;\delta f\right]\equiv {\displaystyle \underset{j_n=1}{\overset{TP}{\sum }}{s}_V^{\left(n\right)}\left[i;j_{n-2};\cdots ;j_1;U^{\left(n\right)}\left(x\right);f\right]\delta f_{j_n}};\\ i=1,\cdots ;2^{n-1};\end{array}$ (44)

$\begin{array}{l}{Q}_2^{\left(n+1\right)}\left[2^{n-1};U^{\left(n+1\right)}\left(2^n;j_{n-1};\cdots ;j_1;x\right);f;\delta f\right]\\ \triangleq \frac{\partial Q_A^{\left(n\right)}\left[2^{n-1};j_{n-1};\cdots ;j_1;U^{\left(n\right)}\left(2^{n-1};j_{n-2};\cdots ;j_1;x\right);f\right]}{\partial f}\partial f\\ -\frac{\partial \left\{A{M}^{\left(n\right)}\left[2^{n-1}\times 2^{n-1};U^{\left(n\right)}\left(2^{n-1};x\right);f\right]A^{(n)}\left(2^{n-1};j_{n-1};\cdots ;j_1;x\right)\right\}}{\partial f}\partial f;\end{array}$ (45)

$\begin{array}{l}{B}_V^{\left(n+1\right)}\left[2^n;U^{\left(n+1\right)}\left(2^n;j_{n-1};j_1;x\right);V^{\left(n+1\right)}\left(2^n;j_{n-1};\cdots ;j_1;x\right);f;\delta f\right]\\ \triangleq \left(\begin{array}{c}{B}_V^{\left(n\right)}\left[2^{n-1};U^{\left(n\right)}\left(2^{n-1};j_{n-2};\cdots ;j_1;x\right);V^{\left(n\right)}\left(2^{n-1};j_{n-2};\cdots ;j_1;x\right);f;\delta f\right]\\ \delta B_A^{\left(n\right)}\left[2^{n-1};U^{\left(n\right)}\left(2^{n-1};j_{n-2};\cdots ;j_1;x\right);A^{(n)}\left(2^{n-1};j_{n-1};\cdots ;j_1;x\right);f\right]\end{array}\right).\end{array}$ (46)

Solving the (n+1)^th-LVSS would require $O\left(T{P}^{n+1}\right)$ large-scale computations, which is unrealistic for large-scale systems comprising many parameters. The (n+1)^th-FASAM-N circumvents the need for solving the (n+1)^th-LVSS by deriving an alternative expression for the indirect-effect term defined in Equation , in which the function $V^{\left(n+1\right)}\left(2^n;j_{n-1};\cdots ;j_1;x\right)$ is replaced by a (n+1)^th-level adjoint function, denoted as $A^{\left(n+1\right)}\left(2^n;j_n;\cdots ;j_1;x\right)\triangleq {\left[a^{\left(n+1\right)}\left(1;j_n;\mathrm{...};j_1;x\right),\cdots ,a^{\left(n+1\right)}\left(2^n;j_n;\cdots ;j_1;x\right)\right]}^{†}\in H_{n+1}\left({\Omega }_x\right)$ , which is independent of parameter variations. The elements of the Hilbert space $H_{n+1}\left({\Omega }_x\right)$ are block-vectors of the form

${\Psi }^{\left(n+1\right)}\left(2^n;j_n;\cdots ;j_1;x\right)\triangleq {\left[{\psi }^{\left(n+1\right)}\left(1;j_n;\cdots ;j_1;x\right),\cdots ,{\psi }^{\left(n+1\right)}\left(2^n;j_n;\cdots ;j_1;x\right)\right]}^{†}$ , comprising $2^n$ TD-dimensional block-vectors of the form

${\psi }^{\left(n+1\right)}\left(i;x\right)\triangleq {\left[{\psi }_1^{\left(n+1\right)}\left(i;x\right),\cdots ,{\psi }_{TD}^{\left(n+1\right)}\left(i;x\right)\right]}^{†}\in H_1\left({\Omega }_x\right)$ , $i=1,\cdots ,2^{n+1}$ . The inner

product of two vectors ${\Psi }^{\left(n+1\right)}\left(2^n;j_n;\cdots ;j_1;x\right)$ and ${\Phi }^{\left(n+1\right)}\left(2^n;j_n;\cdots ;j_1;x\right)$ in the Hilbert space $H_{n+1}\left({\Omega }_x\right)$ is denoted as ${\langle {\Psi }^{\left(n+1\right)}\left(2^n;x\right),{\Phi }^{\left(n+1\right)}\left(2^n;x\right)\rangle }_{n+1}$

and defined as follows:

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

The (n+1)^th-Level Adjoint Sensitivity System [abbreviated as: (n+1)^th-LASS ] for the (n+1)^th-level adjoint function $A^{\left(n+1\right)}\left(2^n;j_n;\cdots ;j_1;x\right)$ is obtained by using the inner product defined in Equation (47), as follows:

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

where the quantity ${\left\{{\left[P^{\left(n+1\right)}\left(U^{\left(n+1\right)};A^{\left(n+1\right)};V^{\left(n+1\right)};f\right)\right]}_{\partial {\Omega }_x}\right\}}_{{\alpha }^0}$ denotes the corres-

ponding bilinear concomitant on the domain’s boundary, evaluated at the nominal values for the parameters and respective state functions.

In terms of the (n+1)^th-level adjoint function $A^{\left(n+1\right)}\left(2^n;j_n;\cdots ;j_1;x\right)$ , the indirect-effect term defined by Equation (35) will have the following expression:

$\begin{array}{l}{\left\{\delta R^{\left(n\right)}\left[j_n;\cdots ;j_1;U^{\left(n\right)};A^{(n)};\alpha ;V^{\left(n\right)};\delta A^{(n)}\right]\right\}}_{ind}\\ =-{\left\{{\left[{\hat{P}}^{\left(n+1\right)}\left(U^{\left(n+1\right)};A^{\left(n+1\right)};\delta f\right)\right]}_{\partial {\Omega }_x}\right\}}_{{\alpha }^0}\\ +{\left\{{\langle A^{\left(n+1\right)}\left(2^n;j_n;\cdots ;j_1;x\right),Q_V^{\left(n+1\right)}\left[2^n;U^{\left(n+1\right)}\left(2^n;j_{n-1};\cdots ;j_1;x\right);f;\delta f\right]\rangle }_n\right\}}_{{\alpha }^0}.\end{array}$ (49)

where $A^{\left(n+1\right)}\left(2^n;j_n;\cdots ;j_1;x\right)$ is the solution of the following (n+1)^th-LASS:

$\begin{array}{l}A{M}^{\left(n+1\right)}\left[2^n\times 2^n;U^{\left(n+1\right)}\left(2^n;x\right);f\right]A^{\left(n+1\right)}\left(2^n;j_n;\cdots ;j_1;x\right)\\ =Q_A^{\left(n+1\right)}\left[2^n;j_n;\cdots ;j_1;U^{\left(n+1\right)}\left(2^n;j_{n-1};\cdots ;j_1;x\right);f\right],\end{array}$ (50)

subject to boundary conditions represented in operator form as follows:

$\begin{array}{l}{\left\{B_A^{\left(n+1\right)}\left[2^n;U^{\left(n\right)}\left(2^n;j_{n-1};\cdots ;j_1;x\right);A^{\left(n\right)}\left(2^n;j_n;\cdots ;j_1;x\right);f\right]\right\}}_{{\alpha }^0}=0\left[2^n\right],\\ x\in \partial {\Omega }_x\left({\alpha }^0\right);j_1=1,\cdots ,TF;j_2=1,\cdots ,j_1;\cdots ;j_n=1,\cdots ,j_{n-1}.\end{array}$ (51)

In Equation (49), the quantity ${\left\{{\left[{\hat{P}}^{\left(n+1\right)}\left(U^{\left(n+1\right)};A^{\left(n+1\right)};\delta f\right)\right]}_{\partial {\Omega }_x}\right\}}_{{\alpha }^0}$ denotes

residual boundary terms which may have not vanished automatically after having used the boundary conditions provided in Equations (37) and (51) to eliminate from Equation (48) all unknown values of the (n+1)^th-level variational function $V^{\left(n+1\right)}\left(2^n;j_{n-1};\cdots ;j_1;x\right)$ .

The quantities that appear in the definition of the (n+1)^th-LASS represented by Equations (50) and (51), are defined as follows:

$\begin{array}{l}A{M}^{\left(n+1\right)}\left[2^n\times 2^n;U^{\left(n+1\right)}\left(2^n;j_{n-1};\cdots ;j_1;x\right);f\right]\\ \triangleq {\left[V{M}^{\left(n+1\right)}\left(2^n\times 2^n;U^{\left(n+1\right)};f\right)\right]}^{*}\\ =\left(\begin{array}{cc}{\left\{{\left[V{M}^{\left(n\right)}\left(2^{n-1}\times 2^{n-1}\right)\right]}^{*}\right\}}^{†}& {\left\{{\left[V{M}_{21}^{\left(n+1\right)}\left(2^{n-1}\times 2^{n-1}\right)\right]}^{*}\right\}}^{†}\\ 0\left[2^{n-1}\times 2^{n-1}\right]& {\left\{{\left[V{M}_{22}^{\left(n+1\right)}\left(2^{n-1}\times 2^{n-1}\right)\right]}^{*}\right\}}^{†}\end{array}\right),\end{array}$ (52)

$\begin{array}{l}{Q}_A^{\left(n+1\right)}\left[2^n;j_n;\cdots ;j_1;U^{\left(n+1\right)}\left(2^n;x\right);f\right]\triangleq \{q_A^{\left(n+1\right)}\left[1;j_n;\cdots ;j_1;U^{\left(n+1\right)}\left(2^n;x\right);f\right],\\ \cdots ,{q_A^{\left(n+1\right)}\left[2^n;j_n;\cdots ;j_1;U^{\left(n+1\right)}\left(2^n;x\right);f\right]\}}^{†},j_1=1,\cdots ,TF;\cdots ;j_n=1,\cdots ,j_{n-1};\end{array}$ (53)

$q_A^{\left(n+1\right)}\left[1;j_n;\cdots ;j_1;U^{\left(n+1\right)}\left(x\right);f\right]\triangleq \partial S^{\left(n\right)}\left[j_n;\cdots ;j_1;U^{\left(n+1\right)};f\right]/\partial u\left(x\right);$ (54)

$q_A^{\left(n+1\right)}\left[2;j_n;\cdots ;j_1;U^{\left(n+1\right)}\left(x\right);\alpha \right]\triangleq \partial S^{\left(n\right)}\left[j_n;\cdots ;j_1;U^{\left(n+1\right)};\alpha \right]/\partial a^{\left(1\right)}\left(x\right);$ (55)

$\begin{array}{l}\text{for}n\ge 3:q_A^{\left(n+1\right)}\left[2^k+i;j_n;\cdots ;j_1;U^{\left(n+1\right)}\left(x\right);f\right]\\ \triangleq \frac{\partial S^{\left(n\right)}\left[j_n;\cdots ;j_1;U^{\left(n+1\right)};f\right]}{\partial a^{\left(k+1\right)}\left(i;j_k;\cdots ;j_1;x\right)};k=1,\cdots ,n-1;i=1,\cdots ,2^k.\end{array}$ (56)

The final expression of the total differential expressed by Equation (17) is obtained by inserting in Equation (34) the expression for the indirect-effect term obtained in Equation (49), which yields the following expression for the n^th-order sensitivities $R^{\left(n+1\right)}\left[j_{n+1};\cdots ;j_1;U^{\left(n+1\right)};A^{\left(n+1\right)};f\right]\equiv {\partial }^{n+1}R\left[u\left(x\right);f\left(\alpha \right)\right]/\partial f_{j_{n+1}}\cdots \partial f_{j_1}$ of the response $R\left[u\left(x\right);f\left(\alpha \right)\right]$ with respect to the components $f_{j_1},\cdots ,f_{j_{n+1}}$ of the “features function” $f\left(\alpha \right)$ , for $j_1=1,\cdots ,TF;\cdots ;j_n=1,\cdots ,j_{n-1}$ :

$\begin{array}{l}{R}^{\left(n+1\right)}\left[j_{n+1};\cdots ;j_1;U^{\left(n+1\right)};A^{\left(n+1\right)};f\left(\alpha \right)\right]\equiv {\partial }^{n+1}R\left[u\left(x\right);\alpha \right]/\partial f_{j_{n+1}}\cdots \partial f_{j_1}\\ =\frac{\partial R^{\left(n\right)}\left[j_n;\cdots ;j_1;U^{\left(n\right)}\left(x\right);A^{\left(n\right)}\left(x\right);f\left(\alpha \right)\right]}{\partial f_{j_{n+1}}}-\frac{{\left[\partial {\hat{P}}^{\left(n+1\right)}\left(U^{\left(n+1\right)};A^{\left(n+1\right)};\delta f\right)\right]}_{\partial {\Omega }_x}}{\partial f_{j_{n+1}}}\\ +{\displaystyle \underset{i=1}{\overset{2^n}{\sum }}{\langle a^{\left(n+1\right)}\left(i;j_n;\cdots ;j_1;x\right),s_V^{\left(n+1\right)}\left[i;j_{n+1};\cdots ;j_1;U^{\left(n+1\right)}\left(x\right);f\right]\rangle }_1}.\\ \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(n+1\right)}\left[j_{n+1};\cdots ;j_1;U^{\left(n+1\right)}\left(2^n;x\right);A^{\left(n+1\right)}\left(2^n;x\right);f\right]\text{d}{x}_1\cdots \text{d}{x}_{TI}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}}}.\end{array}$ (57)

The expression obtained in Equation (57) is identical with the expression that would be obtained by replacing the index n with (n+1) in the expression obtained in Equation (25), thus completing the proof, by mathematical induction, of the validity/correctness of the conjectured general expressions underlying the n^th-FASAM-N.

3. Conclusions

This work has presented the “n^th-Order Feature Adjoint Sensitivity Analysis Methodology for Nonlinear Systems” (abbreviated as “nt^h-FASAM-N”), which is the most efficient methodology for computing exact expressions of sensitivities of model responses with respect to features of model parameters and, subsequently, with respect to the model parameters themselves. This efficiency stems from the reduction of the number of adjoint computations (which are considered to be “large-scale” computations), by comparison to the extant [31] high-order adjoint sensitivity analysis methodology n^th-CASAM-N (the “n^th-Order Comprehensive Adjoint Sensitivity Analysis Methodology for Nonlinear Systems”), as follows:

1) Comparing the mathematical framework of the n^th-FASAM-N methodology presented in this work to the framework of the n^th-CASAM-N methodology [31] indicates that the components $f_i\left(\alpha \right),i=1,\cdots ,TF$ , of the “feature function” $f\left(\alpha \right)\triangleq {\left[f_1\left(\alpha \right),\cdots ,f_{TF}\left(\alpha \right)\right]}^{†}$ play within the n^th-FASAM-N the same role as played by the components ${\alpha }_j,j=1,\cdots ,TP$ , of the “vector of primary model parameters” $\alpha \triangleq {\left({\alpha }_1,\cdots ,{\alpha }_{TP}\right)}^{†}$ within the framework of the n^th-CASAM-N [31] . It is paramount to underscore, at the outset, that the total number of model parameters is always larger (usually by a wide margin) than the total number of components of the feature function $f\left(\alpha \right)$ , i.e., $TP\gg TF$ .

2) The 1^st-FASAM-N and the 1^st-CASAM-N methodologies require a single large-scale “adjoint” computations for solving the 1^st-LASS (1^st-Level Adjoint Sensitivity System), so they are equally efficient for computing the exact expressions of the first-order sensitivities of a model response to the model’s uncertain parameters, boundaries, and internal interfaces.

3) For computing the exact expressions of the second-order response sensitivities with respect to the primary model’s parameters, the 2^nd-FASAM-N methodology [32] requires as many large-scale “adjoint” computations as there are “feature functions of parameters” $f_i\left(\alpha \right),i=1,\cdots ,TF$ (where TF denotes the total number of feature functions) for solving the left-side of the 2^nd-LASS with TF distinct sources on its right-side. By comparison, the 2^nd-CASAM-N methodology [31] requires TP (where TP denotes the total number of model parameters) large-scale computations for solving the same left-side of the 2^nd-LASS but with TP distinct sources. Since $TF\ll TP$ , the 2^nd-FASAM-N methodology is considerably more efficient than the 2^nd-CASAM-N methodology for computing the exact expressions of the second-order sensitivities of a model response to the model’s uncertain parameters, boundaries, and internal interfaces.

4) For computing the exact expressions of the third-order response sensitivities with respect to the primary model’s parameters, it can be deduced from the general theory underlying the 3^rd-FASAM-N presented in this work that, for $n=3$ , the 3^rd-FASAM-N requires at most $TF\left(TF+1\right)/2$ large-scale “adjoint” computations while the 3^rd-CASAM-N methodology [31] requires at most. The same computational-count of ”large-scale computations” caries over when computing the higher-order sensitivities, i.e., the formula for calculating the “number of large-scale adjoint computations” is formally the same for both the n^th-FASAM-N and the n^th-CASAM-N methodologies, but the “variable” in the formula for determining the number of adjoint computations for the n^th-FASAM-N methodology is TF (i.e., the total number of feature functions) while the counterpart for the formula for determining the number of adjoint computations for the n^th-CASAM-N is methodology is TP (i.e., the total number of model parameters). Since $TF\ll TP$ , it follows that the higher the order of computed sensitivities, the more efficient the n^th-FASAM-N methodology becomes by comparison to the n^th-CASAM-N methodology.

5) When a model has no “feature” functions of parameters, but only comprises primary parameters, the n^th-FASAM-N methodology becomes identical to the n^th-CASAM-N methodology [31] .

Both the n^th-FASAM-N and the n^th-CASAM-N methodologies are formulated in linearly increasing higher-dimensional Hilbert spaces as opposed to exponentially increasing parameter-dimensional spaces thus overcoming the curse of dimensionality in sensitivity analysis of nonlinear systems. Both the n^th-FASAM-N and the n^th-CASAM-N methodologies are incomparably more efficient and more accurate than any other methods (statistical, finite differences, etc.) for computing exact expressions of response sensitivities (of any order) with respect to the model’s uncertain parameters, boundaries, and internal interfaces.

The question of “when to stop computing progressively higher-order sensitivities?” has been addressed by Cacuci [29] [31] in conjunction with the question of convergence of the Taylor-series expansion of the response in terms of the uncertain model parameters [see Equations (65)‒(67) in Appendix A]. This Taylor-series expansion is the fundamental premise for obtaining meaningful (convergent) expressions provided by the “propagation of errors” methodology, as originally proposed heuristically by Tukey [36] , for the cumulants of the model response distribution in the phase-space of model parameters. The convergence of this Taylor-series, which depends on both the response sensitivities to parameters and the uncertainties associated with the parameter distribution, must be ensured. This can be done by ensuring that the combination of parameter uncertainties and response sensitivities are sufficiently small to fall inside the radius of convergence of this Taylor-series expansion.

The companion work “Part 2” [35] will present an illustrative paradigm application of the nth-FASAM-N methodology to the Nordheim-Fuchs reactor dynamics/safety model [33] [34] .

Appendix A: Mathematical Modeling of a Generic Nonlinear System Comprising Functions (“Features”) of Uncertain Parameters and Boundaries

The mathematical model that underlies the numerical evaluation of a process and/or state of a physical system comprises equations that relate the system’s independent variables and parameters to the system’s state/dependent variables. These coupled equations, which are in general nonlinear, can be represented generically in operator form as follows:

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

$\begin{array}{l}B\left[u\left(x\right);g\left(\alpha \right);\lambda \left(\alpha \right);\omega \left(\alpha \right)\right]=C\left[g\left(\alpha \right);\lambda \left(\alpha \right);\omega \left(\alpha \right)\right],\\ x\in \partial {\Omega }_x\left[\lambda \left(\alpha \right);\omega \left(\alpha \right)\right]\end{array}$ . (59)

The results computed using a mathematical model are customarily called “model responses” (or “system responses” or “objective functions” or “indices of performance”). As has been discussed by Cacuci (2022, 2023a, 2023b), all responses can be fundamentally analyzed in terms of the following generic integral representation:

$R\left[u\left(x\right);f\left(\alpha \right)\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);g\left(\alpha \right);h\left(\alpha \right);x\right]\text{d}{x}_1\cdots \text{d}{x}_{TI}}}$ , (60)

where $S\left[u\left(x\right);g\left(\alpha \right);h\left(\alpha \right);x\right]$ is a suitably differentiable nonlinear function of $u\left(x\right)$ and of $\alpha$ .

Without loss of generality, the quantities which appear in Equations (58)-(60) can be considered to be real-valued, having the following meanings:

1) Matrices are 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 equalities in this work are considered to hold in the weak (“distributional”) sense. Both sides of Equations (58)-(60) may contain “generalized functions/functionals”, particularly Dirac-distributions and derivatives thereof.

2) The TP-dimensional column-vector $\alpha \triangleq {\left({\alpha }_1,\cdots ,{\alpha }_{TP}\right)}^{†}\in {ℝ}^{TP}$ represents the “vector of primary model parameters” and has components denoted as ${\alpha }_1,\cdots ,{\alpha }_{TP}$ , where TP denotes the “total number of parameters” involved in the model under consideration. These model parameters usually stem from processes that are external to the system under consideration and afflicted by uncertainties (i.e., they are not known precisely). The known characteristics of the model parameters usually includes their nominal (expected/mean) values and, possibly, higher-order moments or cumulants (i.e., variance/covariances, skewness, kurtosis) of their unknown distribution; these characteristics are usually determined from experimental data and/or processes external to the physical system under consideration. Occasionally, just the lower and the upper bounds may be known for some model parameters, expressed by inequality and/or equality constraints that delimit the ranges of the system’s parameters are known. Without loss of generality, the imprecisely known model parameters can be considered to be real-valued scalar quantities. It is important to note that the components of the vector $\alpha$ include not only parameters that appear in in Equations (58) and (59), but also include parameters that may specifically occur only in the definition of the model’s response provided in Equation (60). 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 monograph to denote “nominal” or “mean” values.

3) The TI-dimensional column vector $x\triangleq {\left(x_1,\cdots ,x_{TI}\right)}^{†}\in {ℝ}^{TI}$ comprises the model’s independent variables, denoted as $x_i,i=1,\cdots ,TI$ , where the sub/superscript “TI” denotes the “total number of independent variables.” The vector $x\in {ℝ}^{TI}$ is considered to be defined on a phase-space domain denoted as ${\Omega }_x\left(x\right)\triangleq \left\{{\lambda }_i\left(\alpha \right)\le x_i\le {\omega }_i\left(\alpha \right);i=1,\cdots ,TI\right\}$ , including the particular cases when ${\lambda }_i\left(\alpha \right)=-\infty$ , ${\omega }_i\left(\alpha \right)=\infty$ for some independent variables $x_i,i=1,\cdots ,TI$ . The domain boundary $\partial {\Omega }_x\left[\lambda \left(\alpha \right);\omega \left(\alpha \right)\right]\triangleq \left\{{\lambda }_i\left(\alpha \right)\cup {\omega }_i\left(\alpha \right),i=1,\cdots ,TI\right\}$ of ${\Omega }_x\left[f\left(\alpha \right)\right]$ is defined to comprise the set of all of the endpoints ${\lambda }_i\left(\alpha \right),{\omega }_i\left(\alpha \right),i=1,\cdots ,TI$ . For subsequent mathematical developments, it is convenient to consider that the endpoints ${\lambda }_i\left(\alpha \right),{\omega }_i\left(\alpha \right),i=1,\cdots ,TI$ are components of column-vectors $\lambda \left(\alpha \right)\triangleq {\left[{\lambda }_1\left(\alpha \right),\cdots ,{\lambda }_{TI}\left(\alpha \right)\right]}^{†}$ and $\omega \left(\alpha \right)\triangleq {\left[{\omega }_1\left(\alpha \right),\cdots ,{\omega }_{TI}\left(\alpha \right)\right]}^{†}$ , respectively These endpoints depend on the physical system’s geometrical dimensions, which may be imprecisely known because of manufacturing tolerances, and are considered therefore to be components of the vector $\alpha \triangleq {\left({\alpha }_1,\cdots ,{\alpha }_{TP}\right)}^{†}\in {ℝ}^{TP}$ of primary model parameters. Furthermore, the boundary-endpoints ${\lambda }_i\left(\alpha \right),{\omega }_i\left(\alpha \right),i=1,\cdots ,TI$ , may also depend on the parameters that define the material properties of the respective medium. For example, in models based on diffusion theory, the boundary conditions for materials facing air/vacuum are imposed on a physics-based mathematical construct called the “extrapolated boundary” of the respective spatial domain. The “extrapolated boundary” depends both on the imprecisely known physical dimensions of the system’s materials and also on the material’s properties, such as atomic number densities and microscopic transport cross sections. Therefore, the boundary end-points can be considered, in general, to be functions of (some of) the primary model parameters.

4) The components $u_i\left(x\right),i=1,\cdots ,TD$ of the TD-dimensional column vector $u\left(x\right)\triangleq {\left[u_1\left(x\right),\cdots ,u_{TD}\left(x\right)\right]}^{†}$ represent the model’s dependent variables (also called “state functions”); the abbreviation “TD” denotes “total number of dependent variables.”

5) The vector $g\left(\alpha \right)\triangleq \left[g_1\left(\alpha \right),\cdots ,g_{TG}\left(\alpha \right)\right]$ is a TG-dimensional vector having components $g_i\left(\alpha \right),i=1,\cdots ,TG$ , which are real-valued functions of (some of) the primary model parameters $\alpha \in {ℝ}^{TP}$ . The quantity TG denotes the total number of such functions which appear exclusively in the definition of the model’s underlying equations. Such functions customarily appear in models in the form of correlations that describe “features” of the system under consideration, such as material properties, flow regimes. etc. Usually, the number of functions $g_i\left(\alpha \right)$ is considerably smaller than the total number of model parameters, i.e., $TG\ll TP$ . For example, the numerical model (Cacuci and Fang, 2023) of the OECD/NEA reactor physics benchmark (Valentine 2006) comprises 21,976 uncertain primary model parameters (including microscopic cross sections and isotopic number densities) but the neutron transport equation, which is solved to determine the neutron flux distribution within the benchmark, does not use these primary parameters directly but instead uses several hundreds of “group-averaged macroscopic cross sections” which are functions/features of the microscopic cross sections and isotopic number densities (which in turn are uncertain quantities that would be components of the vector of primary model parameters). In particular, a component $g_j\left(\alpha \right)$ may simply be one of the primary model parameters ${\alpha }_j$ , i.e., $g_j\left(\alpha \right)\equiv {\alpha }_j$ .

6) The TD-dimensional column vector $N\left[u\left(x\right);x;g\left(\alpha \right)\right]\triangleq {\left(N_1,\cdots ,N_{TD}\right)}^{†}$ comprises components $N_i\left[u\left(x\right);x;g\left(\alpha \right)\right],i=1,\cdots ,TD$ , which are operators (including differential, difference, integral, distributions, and/or finite or infinite matrices) acting nonlinearly (in general) on the dependent variables $u\left(x\right)$ , the independent variables $x$ and on the functions $g\left(\alpha \right)$ of model parameters $\alpha$ .

7) The TD-dimensional column vector $Q\left[u\left(x\right);x;g\left(\alpha \right)\right]\triangleq {\left(q_1,\cdots ,q_{TD}\right)}^{†}$ , having components $q_i\left[u\left(x\right);x;g\left(\alpha \right)\right],i=1,\cdots ,TD$ , denotes inhomogeneous source terms, which usually depend nonlinearly on uncertain parameters $\alpha$ .

8) The components of $B\left[u\left(x\right);g\left(\alpha \right);x\right]$ are nonlinear operators, while the components of $C\left[x,g\left(\alpha \right)\right]$ represent inhomogeneous boundary sources, all defined on the boundary $\partial {\Omega }_x$ .

9) The integral representation of the response provided in Equation (60) can represent “averaged” and/or “point-valued” quantities in the phase-space of independent variables. For example, if $R\left[u\left(x\right);f\left(\alpha \right)\right]$ represents the computation or the measurement (which would be a “detector-response”) of a quantity of interest at a point $x_d$ in the phase-space of independent variables, then $S\left[u\left(x\right);g\left(\alpha \right);h\left(\alpha \right);x\right]$ would contain a Dirac-delta functional of the form $\delta \left(x-x_d\right)$ . Responses that represent “differentials/derivatives of quantities” would contain derivatives of Dirac-delta functionals in the definition of $S\left[u\left(x\right);g\left(\alpha \right);h\left(\alpha \right);x\right]$ . The vector $h\left(\alpha \right)\triangleq \left[h_1\left(\alpha \right),\cdots ,h_{TH}\left(\alpha \right)\right]$ , having components $h_i\left(\alpha \right),i=1,\cdots ,TH$ , which appears among the arguments of the function $S\left[u\left(x\right);g\left(\alpha \right);h\left(\alpha \right);x\right]$ , represents functions of primary parameters that often appear solely in the definition of the response but do not appear in the mathematical definition of the model, i.e., in Equations (58) and (59). The quantity TH denotes the total number of such functions which appear exclusively in the definition of the model’s response. Evidently, the response will depend directly and/or indirectly (through the “feature”-functions) on all of the primary model parameters. This fact has been indicated in Equation (60) by using the vector-valued function $f\left(\alpha \right)$ as an argument in the definition of the response $R\left[u\left(x\right);f\left(\alpha \right)\right]$ to represent the concatenation of all of the “features” of the model and response under consideration. The vector $f\left(\alpha \right)$ of “model features” is thus defined as follows:

$\begin{array}{c}f\left(\alpha \right)\triangleq {\left[g\left(\alpha \right);h\left(\alpha \right);\lambda \left(\alpha \right);\omega \left(\alpha \right)\right]}^{†}\\ \triangleq {\left[f_1\left(\alpha \right),\cdots ,f_{TF}\left(\alpha \right)\right]}^{†};\text{}\text{}TF\triangleq TG+TH+2TI.\end{array}$ (61)

As defined in Equation (61), the quantity TF denotes the total number of “feature functions of the model’s parameters” which appear in the definition of the nonlinear model’s response and underlying equations.

Solving Equations (58) and (59) at the nominal parameter values, ${\alpha }^0\triangleq {\left[{\alpha }_1^0,\cdots ,{\alpha }_i^0,\cdots ,{\alpha }_{TP}^0\right]}^{†}$ , provides the “nominal solution” $u^0\left(x\right)$ ; which means that the vectors $u^0\left(x\right)$ and ${\alpha }^0$ satisfy the following equations:

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

$\begin{array}{l}B\left[u^0\left(x\right);g\left({\alpha }^0\right);\lambda \left({\alpha }^0\right);\omega \left({\alpha }^0\right)\right]\\ =C\left[g\left({\alpha }^0\right);\lambda \left({\alpha }^0\right);\omega \left({\alpha }^0\right)\right],x\in \partial {\Omega }_x\left[\lambda \left({\alpha }^0\right);\omega \left({\alpha }^0\right)\right].\end{array}$ (63)

Using the nominal parameter values ${\alpha }^0$ together with the “nominal solution” $u^0\left(x\right)$ in Equation (60) yields the nominal value of the response, namely:

$R^0\triangleq R\left[u^0\left(x\right);f\left({\alpha }^0\right)\right]\triangleq {\displaystyle \underset{{\lambda }_1\left({\alpha }^0\right)}{\overset{{\omega }_1\left({\alpha }^0\right)}{\int }}\cdots {\displaystyle \underset{{\lambda }_{TI}\left({\alpha }^0\right)}{\overset{{\omega }_{TI}\left({\alpha }^0\right)}{\int }}S\left[u^0\left(x\right);g\left({\alpha }^0\right);h\left({\alpha }^0\right);x\right]\text{d}{x}_1\cdots \text{d}{x}_{TI}}}$ .(64)

In view of Equation (64), each model response $R_k\left[f\left(\alpha \right)\right],\text{}\text{}\text{}k=1,\cdots ,TR$ , where TR denotes the “total number of responses,” can be considered to depend directly on the feature functions $f\left(\alpha \right)$ , and would therefore admit a Taylor-series expansion around the nominal value $f^0\triangleq f\left({\alpha }^0\right)$ , having the following form:

$\begin{array}{c}{R}_k\left[f\left(\alpha \right)\right]=R_k\left(f^0\right)+{\displaystyle \underset{j_1=1}{\overset{TF}{\sum }}{\left\{\frac{\partial R_k\left(f\right)}{\partial f_{j_1}}\right\}}_{f^0}\delta f_{j_1}}\\ +\frac{1}{2}{\displaystyle \underset{j_1=1}{\overset{TF}{\sum }}{\displaystyle \underset{j_2=1}{\overset{TF}{\sum }}{\left\{\frac{{\partial }^2{R}_k\left(f\right)}{\partial f_{j_1}\partial f_{j_2}}\right\}}_{f^0}}}\delta f_{j_1}\delta f_{j_2}+\cdots \end{array}$ (65)

where $\delta f_j\triangleq \left[f_j\left(\alpha \right)-f_j^0\right];f_j^0\triangleq f_j\left({\alpha }^0\right);j=1,\cdots ,TF$ . The “sensitivities of the model response with respect to the (feature) functions” are naturally defined as being the functional derivatives of $R_k\left[f\left(\alpha \right)\right]$ with respect to the components (“features”) $f_j\left(\alpha \right)$ of $f\left(\alpha \right)$ . The notation ${\{·\}}_{f^0}$ indicates that the quantity enclosed within the braces is to be evaluated at the nominal values $f^0\triangleq f\left({\alpha }^0\right)$ . Since $TF\ll TP$ , the computations of the functional derivatives of $R_k\left[f\left(\alpha \right)\right]$ with respect to the functions $f_j\left(\alpha \right)$ , which appear in Equation (65), will be considerably less expensive computationally than the computation of the functional derivatives involved in the Taylor-series of the response with respect to the model parameters. The functional derivatives of the response with respect to the parameters can be obtained from the functional derivatives of the response with respect to the “feature” functions $f_j\left(\alpha \right)$ by simply using the chain rule, i.e.:

$\begin{array}{l}{\left\{\frac{\partial R_k\left(\alpha \right)}{\partial {\alpha }_{j_1}}\right\}}_{{\alpha }^0}={\displaystyle \underset{i_1=1}{\overset{TF}{\sum }}{\left\{\frac{\partial R_k\left(f\right)}{\partial f_{i_1}}\frac{\partial f_{i_1}\left(\alpha \right)}{\partial {\alpha }_{j_1}}\right\}}_{{\alpha }^0}};\\ {\left\{\frac{{\partial }^2{R}_k\left(\alpha \right)}{\partial {\alpha }_{j_1}\partial {\alpha }_{j_2}}\right\}}_{{\alpha }^0}=\frac{\partial }{\partial {\alpha }_{j_2}}{\displaystyle \underset{i_1=1}{\overset{TF}{\sum }}{\left\{\frac{\partial R_k\left(f\right)}{\partial f_{i_1}}\frac{\partial f_{i_1}\left(\alpha \right)}{\partial {\alpha }_{j_1}}\right\}}_{{\alpha }^0}};\end{array}$ (66)

and so on. The evaluation/computation of the functional derivatives $\partial f_{i_1}\left(\alpha \right)/\partial {\alpha }_{j_1}$ , ${\partial }^2{f}_{i_1}\left(\alpha \right)/\partial {\alpha }_{j_1}\partial {\alpha }_{j_2}$ , etc., does not require computations involving the model, and is therefore computationally trivial by comparison to the evaluation of the functional derivatives (“sensitivities”) of the response with respect to either the functions (“features”) $f_j\left(\alpha \right)$ or the model parameters ${\alpha }_i,i=1,\cdots ,TP$ .

The range of validity of the Taylor-series shown in Equation (65) is defined by its radius of convergence. The accuracy—as opposed to the “validity”—of the Taylor-series in predicting the value of the response at an arbitrary point in the phase-space of model parameters depends on the order of sensitivities retained in the Taylor-expansion: the higher the respective order, the more accurate the respective response value predicted by the Taylor-series. In the particular cases when the response happens to be a polynomial function of the “feature” functions $f_j\left(\alpha \right)$ , the Taylor series is actually exact.

In turn, the functions $f_i\left(\alpha \right)$ can also be formally expanded in a multivariate Taylor-series around the nominal (mean) parameter values ${\alpha }^0$ , namely:

$\begin{array}{l}{f}_i\left(\alpha \right)=f_i\left({\alpha }^0\right)+{\displaystyle \underset{j_1=1}{\overset{TP}{\sum }}{\left\{\frac{\partial f_i\left(\alpha \right)}{\partial {\alpha }_{j_1}}\right\}}_{{\alpha }^0}\delta {\alpha }_{j_1}}+\frac{1}{2}{\displaystyle \underset{j_1=1}{\overset{TP}{\sum }}{\displaystyle \underset{j_2=1}{\overset{TP}{\sum }}{\left\{\frac{{\partial }^2{f}_i\left(\alpha \right)}{\partial {\alpha }_{j_1}\partial {\alpha }_{j_2}}\right\}}_{{\alpha }^0}}}\delta {\alpha }_{j_1}\delta {\alpha }_{j_2}\\ +\frac{1}{3!}{\displaystyle \underset{j_1=1}{\overset{TP}{\sum }}{\displaystyle \underset{j_2=1}{\overset{TP}{\sum }}{\displaystyle \underset{j_3=1}{\overset{TP}{\sum }}{\left\{\frac{{\partial }^3{f}_i\left(\alpha \right)}{\partial {\alpha }_{j_1}\partial {\alpha }_{j_2}\partial {\alpha }_{j_3}}\right\}}_{{\alpha }^0}}}}\delta {\alpha }_{j_1}\delta {\alpha }_{j_2}\delta {\alpha }_{j_3}+\cdots ,\end{array}$ (67)

The domain of validity of the Taylor-series in Equation (67) is defined by its own radius of convergence.

Appendix B: The 1^st-FASAM-N: First-Order Function/Feature Adjoint Sensitivity Analysis Methodology for Nonlinear Systems

The 1^st-order G-differential, denoted as $\delta R\left[u^0\left(x\right);f^0\left(\alpha \right);v^{\left(1\right)}\left(x\right);\delta f\left(\alpha \right)\right]$ , of the response $R\left[u\left(x\right);f\left(\alpha \right)\right]$ at a point $\left(u^0;f^0\right)$ in the phase-space of dependent variables $u\left(x\right)$ and feature-functions $f\left(\alpha \right)$ is defined as follows:

$\begin{array}{l}\delta R\left[u^0\left(x\right);f^0\left(\alpha \right);v^{\left(1\right)}\left(x\right);\delta f\left(\alpha \right)\right]\\ \triangleq \{\frac{\text{d}}{\text{d}\epsilon }{\displaystyle \underset{{\lambda }_1\left({\alpha }^0\right)+\epsilon \left(\delta {\lambda }_1\right)}{\overset{{\omega }_1\left({\alpha }^0\right)+\epsilon \left(\delta {\omega }_1\right)}{\int }}\cdots }{{\displaystyle \underset{{\lambda }_{TI}\left({\alpha }^0\right)+\epsilon \left(\delta {\lambda }_{TI}\right)}{\overset{{\omega }_{TI}\left({\alpha }^0\right)+\epsilon \left(\delta {\omega }_{TI}\right)}{\int }}S\left[u^0+\epsilon v^{\left(1\right)},g^0+\epsilon \left(\delta g\right);h^0+\epsilon \left(\delta h\right);x\right]\text{d}{x}_1\cdots \text{d}{x}_{TI}}\}}_{\epsilon =0}.\end{array}$

(68)

The necessary and sufficient conditions for the 1^st-order G-differential $\delta R\left[u^0\left(x\right);f^0\left(\alpha \right);v^{\left(1\right)}\left(x\right);\delta f\left(\alpha \right)\right]$ to be linear in the variations $v^{\left(1\right)}\triangleq \delta u$ and $\delta f$ (and hence admit partial G-derivatives with respect to $u$ and $f$ ) at a point $\left(u^0;f^0\right)$ in the phase-space of dependent variables and feature-functions are as follows:

i) $R\left[u\left(x\right);f\left(\alpha \right)\right]$ satisfies the following inequality:

$\Vert R\left[u^0+\epsilon v^{\left(1\right)};f^0+\epsilon \left(\delta f\right)\right]-R\left(u^0;f^0\right)\Vert \le k\Vert \epsilon \left(u^0;f^0\right)\Vert ,k<\infty$ . (69)

ii) $R\left[u\left(x\right);f\left(\alpha \right)\right]$ satisfies the following condition for a scalar $\epsilon$ and vectors $v_2\text{,}{f}_2$ :

$\begin{array}{l}R\left[u^0+\epsilon v^{\left(1\right)}+\epsilon v_2;f^0+\epsilon \left(\delta f\right)+\epsilon f_2\right]-R\left[u^0+\epsilon v^{\left(1\right)};f^0+\epsilon \left(\delta f\right)\right]\\ -R\left[u^0+\epsilon v_2;f^0+\epsilon f_2\right]+R\left(u^0;f^0\right)=o\left(\epsilon \right).\end{array}$ (70)

In practice, the relations provided in Equations (69) and (70) are seldom used directly since the computation of the expression on the right-side of Equation (68) reveals immediately if the respective expression is linear (or not) in the vectors $v^{\left(1\right)}\left(x\right)$ and/or $\delta f\left(\alpha \right)$ .

Numerical methods (e.g., Newton’s method and variants thereof) for solving Equations (58) and (59) also require the existence of the first-order G-derivatives of the original model equations. Therefore, the conditions provided in Equations (69) and (70) are henceforth considered to be satisfied by the model responses and also by the operators underlying the physical system modeled by Equations (58)-(60), which implies that all of the operators/functions considered in this work admit G-derivatives.

When the 1^st-order G-variation $\delta R\left[u^0\left(x\right);f^0\left(\alpha \right);v^{\left(1\right)}\left(x\right);\delta f\left(\alpha \right)\right]$ satisfies the conditions provided in (69) and (70), it can be written as follows:

$\begin{array}{l}\delta R\left[u^0\left(x\right);f^0\left(\alpha \right);v^{\left(1\right)}\left(x\right);\delta f\left(\alpha \right)\right]\\ \triangleq {\left\{\text{d}R\left[u\left(x\right);f\left(\alpha \right);\delta f\right]\right\}}_{dir}+{\left\{\text{d}R\left[u\left(x\right);f\left(\alpha \right);v^{\left(1\right)}\left(x\right)\right]\right\}}_{ind}.\end{array}$ (71)

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

$\begin{array}{l}{\left\{\text{d}R\left[u\left(x\right);f;\delta f\right]\right\}}_{dir}\triangleq {\left\{\frac{\partial R\left(u;f\right)}{\partial f}\right\}}_{{\alpha }^0}\delta f\\ \triangleq {\displaystyle \underset{j_1=1}{\overset{TF}{\sum }}{\left\{R^{\left(1\right)}\left[j_1;u\left(x\right);f\left(\alpha \right)\right]\right\}}_{dir}\delta f_{j_1}},\end{array}$ (72)

where:

$\frac{\partial \left[\right]}{\partial f}\delta f\triangleq {\displaystyle \underset{i=1}{\overset{TF}{\sum }}\frac{\partial \left[\right]}{\partial f_i}\delta f_i}={\displaystyle \underset{i=1}{\overset{TG}{\sum }}\frac{\partial \left[\right]}{\partial g_i}\delta g_i}+{\displaystyle \underset{i=1}{\overset{TH}{\sum }}\frac{\partial \left[\right]}{\partial h_i}\delta h_i}+{\displaystyle \underset{i=1}{\overset{TI}{\sum }}\frac{\partial \left[\right]}{\partial {\lambda }_i}\delta {\lambda }_i}+{\displaystyle \underset{i=1}{\overset{TI}{\sum }}\frac{\partial \left[\right]}{\partial {\omega }_i}\delta {\omega }_i},$ (73)

and where:

i) For $j_1=1,\cdots ,TG$ :

${\left\{R^{\left(1\right)}\left[j_1;u\left(x\right);f\left(\alpha \right)\right]\right\}}_{dir}\delta f_{j_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 }}\frac{\partial S\left(u;g;h\right)}{\partial g_i}\text{d}{x}_1\cdots \text{d}{x}_{TI}}}\right\}}_{{\alpha }^0}\delta g_i;$ (74)

ii) For $j_1=TG+1,\cdots ,TG+TH$ :

$\begin{array}{l}{\left\{R^{\left(1\right)}\left[j_1;u\left(x\right);f\left(\alpha \right)\right]\right\}}_{dir}\delta f_{j_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 }}\frac{\partial S\left(u;h\right)}{\partial h_i}\text{d}{x}_1\cdots \text{d}{x}_{TI}}}\right\}}_{{\alpha }^0}\delta h_i;i=j_1-TG;\end{array}$ (75)

iii) For $j_1=TG+TH+1,\cdots ,TG+TH+TI$ :

$\begin{array}{l}{\left\{R^{\left(1\right)}\left[j_1;u\left(x\right);f\left(\alpha \right)\right]\right\}}_{dir}\delta f_{j_1}\\ \triangleq {\left\{{\displaystyle \underset{{\lambda }_1}{\overset{{\omega }_1}{\int }}\text{d}{x}_1\cdots {\displaystyle \underset{{\lambda }_{i-1}}{\overset{{\omega }_{i-1}}{\int }}\text{d}{x}_{i-1}{\displaystyle \underset{{\lambda }_{1+1}}{\overset{{\omega }_{i+1}}{\int }}\text{d}{x}_{i+1}\cdots }}{\displaystyle \underset{{\lambda }_{TI}}{\overset{{\omega }_{TI}}{\int }}\text{d}{x}_{TI}S\left[u\left(x_1,.,{\omega }_i\left(\alpha \right),.,x_{N_x}\right);g;h\right]}}\right\}}_{{\alpha }^0}\delta {\omega }_i\left(\alpha \right),\\ \text{where}i=j_1-TG-TH;\end{array}$ (76)

iv) For $j_1=TG+TH+TI+1,\cdots ,TG+TH+2TI$ :

$\begin{array}{l}{\left\{R^{\left(1\right)}\left[j_1;u\left(x\right);f\left(\alpha \right)\right]\right\}}_{dir}\delta f_{j_1}\\ \triangleq -{\left\{{\displaystyle \underset{{\lambda }_1}{\overset{{\omega }_1}{\int }}\text{d}{x}_1\cdots {\displaystyle \underset{{\lambda }_{i-1}}{\overset{{\omega }_{i-1}}{\int }}\text{d}{x}_{i-1}{\displaystyle \underset{{\lambda }_{1+1}}{\overset{{\omega }_{i+1}}{\int }}\text{d}{x}_{i+1}}}\cdots {\displaystyle \underset{{\lambda }_{TI}}{\overset{{\omega }_{TI}}{\int }}\text{d}{x}_{TI}S\left[u\left(x_1,.,{\lambda }_i\left(\alpha \right),.,x_{N_x}\right);g;h\right]}}\right\}}_{{\alpha }^0}\delta {\lambda }_i,\\ \text{where}i=j_1-TG-TH-TI.\end{array}$ (77)

The notation on the left-side of Equation (73) represents the inner product between two vectors (comprising an implied multiplication of a vector with a transposed vector), but the symbol “( $†$ )” which indicates “transposition” has been omitted in order to keep the notation as simple as possible. “Daggers” indicating transposition will also be omitted in other inner products, whenever possible, while avoiding ambiguities.

The direct-effect term can be computed after having solved Equations (58) and (59) to obtain the nominal values, $u^0\left(x\right)$ , of the dependent variables. On the other hand, the quantity ${\left\{\text{d}R\left[u\left(x\right);f\left(\alpha \right);v^{\left(1\right)}\left(x\right)\right]\right\}}_{ind}$ defined in Equation (71) comprises only variations in the state functions and is therefore called the “indirect-effect term” being defined as follows:

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

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)}$ . (79)

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. 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 g$ , $\delta \lambda$ , and $\delta \omega$ .

The first-order relationship between the vectors $v^{\left(1\right)}\left(x\right)$ and the variations $\delta g$ , $\delta \lambda$ , and $\delta \omega$ is determined by solving the equations obtained by applying the definition of the G-differential to Equations (58) and (59), which yields the following equations:

${\left\{\frac{\text{d}}{\text{d}\epsilon }N\left[u^0+\epsilon v^{\left(1\right)}\left(x\right);g\left({\alpha }^0\right)+\epsilon \left(\delta g\right)\right]\right\}}_{\epsilon =0}={\left\{\frac{\text{d}}{\text{d}\epsilon }Q\left[x,g\left({\alpha }^0\right)+\epsilon \left(\delta g\right)\right]\right\}}_{\epsilon =0}$ , (80)

$\begin{array}{l}{\left\{\frac{\text{d}}{\text{d}\epsilon }B\left[u^0+\epsilon v^{\left(1\right)}\left(x\right);g\left({\alpha }^0\right)+\epsilon \left(\delta g\right);\lambda \left({\alpha }^0\right)+\epsilon \left(\delta \lambda \right);\omega \left({\alpha }^0\right)+\epsilon \left(\delta \omega \right)\right]\right\}}_{\epsilon =0}\\ -{\left\{\frac{\text{d}}{\text{d}\epsilon }C\left[g\left({\alpha }^0\right)+\epsilon \left(\delta g\right);\lambda \left({\alpha }^0\right)+\epsilon \left(\delta \lambda \right);\omega \left({\alpha }^0\right)+\epsilon \left(\delta \omega \right)\right]\right\}}_{\epsilon =0}=0,\end{array}$ (81)

Carrying out the differentiations with respect to $\epsilon$ in Equations (80) and (81), and setting $\epsilon =0$ in the resulting expressions yields the following equations:

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

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

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

$N^{\left(1\right)}\left(u;g\right)\triangleq \left\{\frac{\partial N\left(u;g\right)}{\partial u}\right\}\triangleq {\left\{\frac{\partial N_i}{\partial u_j}\right\}}_{TD\times TD};$ (84)

$q_V^{\left(1\right)}\left(u;g;\delta g\right)\triangleq \frac{\partial \left[Q\left(g\right)-N\left(u;g\right)\right]}{\partial g}\delta g\triangleq {\displaystyle \underset{j_1=1}{\overset{TG}{\sum }}{s}_V^{\left(1\right)}\left(j_1;u;g\right)\delta g_{j_1}}$ ; (85)

$\begin{array}{l}{\left\{b_V^{\left(1\right)}\left(u;g;\lambda ;\omega ;v^{\left(1\right)};\delta g;\delta \lambda ;\delta \omega \right)\right\}}_{{\alpha }^0}\\ \triangleq {\left\{\frac{\partial B\left(u;g;\lambda ;\omega \right)}{\partial u}\right\}}_{{\alpha }^0}{v}^{\left(1\right)}+\{\frac{\partial \left[B\left(u;g;\lambda ;\omega \right)-C\left(g;\lambda ;\omega \right)\right]}{\partial g}\delta g\\ +\frac{\partial \left[B\left(u;g;\lambda ;\omega \right)-C\left(g;\lambda ;\omega \right)\right]}{\partial \lambda }\delta \lambda +{\frac{\partial \left[B\left(u;g;\lambda ;\omega \right)-C\left(g;\lambda ;\omega \right)\right]}{\partial \omega }\delta \omega \}}_{{\alpha }^0}.\end{array}$ (86)

The system of equations comprising Equations (82) and (83) is called the “1^st-Level Variational Sensitivity System” (1^st-LVSS). The solution, $v^{\left(1\right)}\left(x\right)$ , of the 1^st-LVSS will be a function of the variations $\delta g$ , $\delta \lambda$ , and $\delta \omega$ . Hence, when the function $v^{\left(1\right)}\left(x\right)$ is introduced into the expression of the indirect-effect term defined in Equation (78), it will introduce dependencies of the response sensitivities on $\delta g$ , $\delta \lambda$ , and $\delta \omega$ , which will be in addition to the dependencies displayed by the direct-effect term defined in Equation (72). As Equations (82) and (83) indicate, every parameter variation $\delta {\alpha }_{j_1}$ , $j_1=1,\cdots ,TP$ , will induce, directly or indirectly, a variations in the model’s state variables. In principle, therefore, for every parameter variation $\delta {\alpha }_{j_1}$ , $j_1=1,\cdots ,TP$ , there would correspond a solution $v^{\left(1\right)}\left(j_1;x\right)$ , $j_1=1,\cdots ,TP$ , of the 1^st-LVSS. Thus, if the effect of every parameter variation were of interest, then the 1^st-LVSS would need to be solved TP times, with distinct right-sides and boundary conditions for each parameter variation $\delta {\alpha }_{j_1}$ , which would require at least TP large-scale computations. If only variations induced by the functions $\delta g$ , $\delta \lambda$ , and/or $\delta \omega$ were of interest, then fewer than TP large-scale computations would be required.

However, solving the 1^st-LVSS can be avoided altogether by using the ideas underlying the “adjoint sensitivity analysis methodology” as originally conceived by Cacuci (1981) and subsequently generalized by Cacuci (2022, 2023a) to enable the computation of arbitrarily high-order response sensitivities to model parameters. Thus, 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 in terms of the solutions of the “1^st-Level Adjoint Sensitivity System” (1^st-LASS), which is constructed by introducing a (real) Hilbert space denoted as $H_1\left({\Omega }_x\right)$ , endowed with an inner product of two vectors $w^{\left(1\right)}\left(x\right)\in H_1$ and $w^{\left(2\right)}\left(x\right)\in H_1$ which is denoted as ${\langle w^{\left(1\right)},w^{\left(2\right)}\rangle }_1$ and defined as follows:

${\langle w^{\left(1\right)},w^{\left(2\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[w^{\left(1\right)}\left(x\right)\cdot w^{\left(2\right)}\left(x\right)\right]\text{d}{x}_1\cdots \text{d}{x}_{TI}}}\right\}}_{{\alpha }^0}$ . (87)

In Equation (87), the “dagger” ( $†$ ), which indicates “transposition,” has been omitted to simplify the notation for the scalar product $w^{\left(1\right)}\left(x\right)\cdot w^{\left(2\right)}\left(x\right)\triangleq {\displaystyle \underset{i=1}{\overset{TD}{\sum }}{w}_i^{\left(1\right)}\left(x\right)w_i^{\left(2\right)}}\left(x\right)$ .

The 1^st-LASS is now constructed by considering a vector $a^{\left(1\right)}\left(x\right)\in H_1$ , which is an element in $H_1\left({\Omega }_x\right)$ but is otherwise arbitrary at this stage, and by using Equation (87) to form the inner product of $a^{\left(1\right)}\left(x\right)\in H_1$ with the relation provided in Equation (82) to obtain:

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

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

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

where ${\left[P^{\left(1\right)}\left(u;g;\lambda ;\omega ;a^{\left(1\right)};v^{\left(1\right)}\right)\right]}_{\partial {\Omega }_x}$ denotes the associated bilinear concomitant evaluated on the domain’s boundary $\partial {\Omega }_x$ , and where $A^{\left(1\right)}\left(u;g\right)\triangleq {\left[N^{\left(1\right)}\left(u;g\right)\right]}^{*}$ denotes the operator formally adjoint to $N^{\left(1\right)}\left(u;g\right)$ .

The symbol ${\left[\right]}^{\ast }$ indicates “formal adjoint” operator.

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

${\{A^{\left(1\right)}(}_{}$

Conflicts of Interest

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

References

[1]	Cukier, R.I., Levine, H.B. and Shuler, K.E. (1978) Nonlinear Sensitivity Analysis of Multiparameter Model Systems. Journal of Computational Physics, 26, 1-42. https://doi.org/10.1016/0021-9991(78)90097-9
[2]	Hora, S.C. and Iman, R.L. (1986) A Comparison of Maximum/Bounding and Bayesian/Monte Carlo for Fault Tree Uncertainty Analysis. Technical Report SAND85-2839, Sandia National Laboratories, Albuquerque. https://doi.org/10.2172/5824798
[3]	Sobol, I.M. (1993) Sensitivity Analysis for Nonlinear Mathematical Models. Ath. Model. Comput. Exp, 1, 407-414.
[4]	Iman, R.L., Helton, J.C. and Campbell, J.E. (1981) An Approach to Sensitivity Analysis of Computer Models, Part 1. Introduction, Input Variable Selection and Preliminary Variable Assessment. Journal of Quality Technology, 13, 174-183. https://doi.org/10.1080/00224065.1981.11978748
[5]	Iman, R.L., Helton, J.C. and Campbell, J.E. (1981) An Approach to Sensitivity Analysis of Computer Models, Part 2. Ranking of Input Variables, Response Surface Validation, Distribution Effect and Technique Synopsis. Journal of Quality Technology, 13, 232-240. https://doi.org/10.1080/00224065.1981.11978763
[6]	Rios Insua, D. (1990) Sensitivity Analysis in Multiobjective Decision Making. In: Rios Insua, D., Eds., Sensitivity Analysis in Multi-Objective Decision Making, Springer, Berlin, 74-126. https://doi.org/10.1007/978-3-642-51656-6_3
[7]	Saltarelli, A., Chan, K. and Scott, E.M. (2000) Sensitivity Analysis. John Wiley & Sons, Chichester.
[8]	Kramer, M.A., Calo, J.M. and Rabitz, H. (1981) An Improved Computational Method for Sensitivity Analysis: Green’s Function Method with “AIM”. Applied Mathematical Modelling, 5, 432-442 https://doi.org/10.1016/S0307-904X(81)80027-3
[9]	Cacuci, D.G. (1981) Sensitivity Theory for Nonlinear Systems: I. Nonlinear Functional Analysis Approach. Journal of Mathematical Physics, 22, 2794-2802. https://doi.org/10.1063/1.525186
[10]	Dunker, A.M. (1984) The Decoupled Direct Method for Calculating Sensitivity Coefficients in Chemical Kinetics. The Journal of Chemical Physics, 81, 2385-2393. https://doi.org/10.1063/1.447938
[11]	Bellman, R.E. (1957) Dynamic Programming. Rand Corporation, Princeton University Press, Princeton.
[12]	Wigner, E.P. (1945) Effect of Small Perturbations on Pile Period. Chicago Report CP-G-3048, Chicago.
[13]	Weiberg, A.M. and Wigner, E.P. (1958) The Physical Theory of Neutron Chain Reactors. University of Chicago Press, Chicago.
[14]	Weisbin, C.R., et al. (1978) Application of Sensitivity and Uncertainty Methodology to Fast Reactor Integral Experiment Analysis. Nuclear Science and Engineering, 66, 307-333. https://doi.org/10.13182/NSE78-3
[15]	Williams, M.L. (1986) Perturbation Theory for Nuclear Reactor Analysis. In: Ronen, Y., Ed., Handbook of Nuclear Reactor Calculations, CRC Press, Boca Raton, 63-188.
[16]	Shultis, J.K. and Faw, R.E. (2000) Radiation Shielding. American Nuclear Society, Illinois.
[17]	Stacey, W.M. (2001) Nuclear Reactor Physics. John Wiley & Sons, New York.
[18]	Cacuci, D.G. (1981) Sensitivity Theory for Nonlinear Systems: II. Extensions to Additional Classes of Responses. Journal of Mathematical Physics, 22, 2794-2802. https://doi.org/10.1063/1.525186
[19]	Práger, T. and Kelemen, F.D. (2014) Adjoint Methods and Their Application in Earth Sciences. In: Faragó, I., Havasi, Á. and Zlatev, Z., Eds., Advanced Numerical Methods for Complex Environmental Models: Needs and Availability, Bentham Science Publishers, Oak Park, 203-275. https://doi.org/10.2174/9781608057788113010011
[20]	Luo, Z., Wang, X. and Liu, D. (2020) Prediction on the Static Response of Structures with Large-Scale Uncertain-But-Bounded Parameters Based on the Adjoint Sensitivity Analysis. Structural and Multidisciplinary Optimization, 61, 123-139. https://doi.org/10.1007/s00158-019-02349-w
[21]	Cacuci, D.G. (2015) Second-Order Adjoint Sensitivity Analysis Methodology for Computing Exactly and Efficiently First- and Second-Order Sensitivities in Large-Scale Linear Systems: I. Computational Methodology. Journal of Computational Physics, 284, 687-699. https://doi.org/10.1016/j.jcp.2014.12.042
[22]	Cacuci, D.G. (2016) Second-Order Adjoint Sensitivity Analysis Methodology for Large-Scale Nonlinear Systems: I. Theory. Nuclear Science and Engineering, 184, 16-30. https://doi.org/10.13182/NSE16-16
[23]	Cacuci, D.G. and Fang, R. (2023) The nth-Order Comprehensive Adjoint Sensitivity Analysis Methodology (nth-CASAM): Overcoming the Curse of Dimensionality in Sensitivity and Uncertainty Analysis, Volume II: Application to a Large-Scale System. Springer Nature Switzerland, Cham.
[24]	Valentine, T.E. (2006) Polyethylene-Reflected Plutonium Metal Sphere Subcritical Noise Measurements, SUB-PU-METMIXED-001. International Handbook of Evaluated Criticality Safety Benchmark Experiments, NEA/NSC/DOC(95)03/I-IX, Organization for Economic Co-Operationand Development, Nuclear Energy Agency, Paris.
[25]	Alcouffe, R.E., Baker, R.S., Dahl, J.A., Turner, S.A. and Ward, R. (2008) PARTISN: A Time-Dependent, Parallel Neutral Particle Transport Code System. LA-UR-08-07258, Los Alamos National Laboratory, Los Alamos.
[26]	Conlin, J.L., Parsons, D.K., Gardiner, S.J., Gray, M., Lee, M.B. and White, M.C. (2013) MENDF71X: Multigroup Neutron Cross-Section Data Tables Based upon ENDF/B-VII.1X. Los Alamos National Laboratory Report LA-UR-15-29571, Los Alamos National Laboratory, Los Alamos. https://doi.org/10.2172/1063914
[27]	Chadwick, M.B., Herman, M., Obložinský, P., Dunn, M.E., Danon, Y., Kahler, A.C., Smith, D.L., Pritychenko, B., Arbanas, G., Brewer, R., et al. (2011) ENDF/B-VII.1: Nuclear Data for Science and Technology: Cross Sections, Covariances, Fission Product Yields and Decay Data. Nuclear Data Sheets, 112, 2887-2996. https://doi.org/10.1016/j.nds.2011.11.002
[28]	Wilson, W.B., Perry, R.T., Shores, E.F., Charlton, W.S., Parish, T.A., Estes, G.P., Brown, T.H., Arthur, E.D., Bozoian, M., England, T.R., et al. (2002) SOURCES4C: A Code for Calculating (α, n), Spontaneous Fission, and Delayed Neutron Sources and Spectra. Proceedings of the American Nuclear Society/Radiation Protection and Shielding Division 12th Biennial Topical Meeting, Santa Fe, 14-18 April 2002.
[29]	Cacuci, D.G. (2022) The nth-Order Comprehensive Adjoint Sensitivity Analysis Methodology (nth-CASAM): Overcoming the Curse of Dimensionality in Sensitivity and Uncertainty Analysis, Volume I: Linear Systems. Springer, New York.
[30]	Williams, M.L. and Engle, W.W. (1977) The Concept of Spatial Channel Theory Applied to Reactor Shielding Analysis. Nuclear Science and Engineering, 62, 92-104. https://doi.org/10.13182/NSE77-A26941
[31]	Cacuci, D.G. (2023) The nth-Order Comprehensive Adjoint Sensitivity Analysis Methodology (nth-CASAM): Overcoming the Curse of Dimensionality in Sensitivity and Uncertainty Analysis, Volume III: Nonlinear Systems. Springer Nature Switzerland, Cham.
[32]	Cacuci, D.G. (2023) Computation of High-Order Sensitivities of Model Responses to Model Parameters. II: Introducing the Second-Order Adjoint Sensitivity Analysis Methodology for Computing Response Sensitivities to Functions/Features of Parameters. Energies, 16, Article 6356. https://doi.org/10.3390/en16176356
[33]	Hetrick, D.L. (1993) Dynamics of Nuclear Reactors. American Nuclear Society, Inc., La Grange Park, 164-174.
[34]	Lamarsh, J.R. (1966) Introduction to Nuclear Reactor Theory. Adison-Wesley Publishing Co., Reading, 491-492.
[35]	Cacuci, D.G. (2024) Introducing the nth-Order Features Adjoint Sensitivity Analysis Methodology for Nonlinear Systems (nth-FASAM-N): II. Illustrative Application.
[36]	Tukey, J.W. (1957) The Propagation of Errors, Fluctuations and Tolerances. Technical Reports No. 10-12, Princeton University, Princeton.