Next Article in Journal
Assessment of Overvoltage and Insulation Coordination in Mixed HVDC Transmission Lines Exposed to Lightning Strikes
Next Article in Special Issue
Comprehensive Second-Order Adjoint Sensitivity Analysis Methodology (2nd-ASAM) Applied to a Subcritical Experimental Reactor Physics Benchmark: I. Effects of Imprecisely Known Microscopic Total and Capture Cross Sections
Previous Article in Journal
Experimental Investigation of Displacer Seal Geometry Effects in Stirling Cycle Machines
Previous Article in Special Issue
Comprehensive Second-Order Adjoint Sensitivity Analysis Methodology (2nd-ASAM) Applied to a Subcritical Experimental Reactor Physics Benchmark: II. Effects of Imprecisely Known Microscopic Scattering Cross Sections
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Towards Overcoming the Curse of Dimensionality: The Third-Order Adjoint Method for Sensitivity Analysis of Response-Coupled Linear Forward/Adjoint Systems, with Applications to Uncertainty Quantification and Predictive Modeling

Department of Mechanical Engineering, University of South Carolina, Columbia, SC 29208, USA
Energies 2019, 12(21), 4216; https://doi.org/10.3390/en12214216
Submission received: 29 September 2019 / Revised: 27 October 2019 / Accepted: 30 October 2019 / Published: 5 November 2019

Abstract

:
This work presents the Third-Order Adjoint Sensitivity Analysis Methodology (3rd-ASAM) for response-coupled forward and adjoint linear systems. The 3rd-ASAM enables the efficient computation of the exact expressions of the 3rd-order functional derivatives (“sensitivities”) of a general system response, which depends on both the forward and adjoint state functions, with respect to all of the parameters underlying the respective forward and adjoint systems. Such responses are often encountered when representing mathematically detector responses and reaction rates in reactor physics problems. The 3rd-ASAM extends the 2nd-ASAM in the quest to overcome the “curse of dimensionality” in sensitivity analysis, uncertainty quantification and predictive modeling. This work also presents new formulas that incorporate the contributions of the 3rd-order sensitivities into the expressions of the first four cumulants of the response distribution in the phase-space of model parameters. Using these newly developed formulas, this work also presents a new mathematical formalism, called the 2nd/3rd-BERRU-PM “Second/Third-Order Best-Estimated Results with Reduced Uncertainties Predictive Modeling”) formalism, which combines experimental and computational information in the joint phase-space of responses and model parameters, including not only the 1st-order response sensitivities, but also the complete hessian matrix of 2nd-order second-sensitivities and also the 3rd-order sensitivities, all computed using the 3rd-ASAM. The 2nd/3rd-BERRU-PM uses the maximum entropy principle to eliminate the need for introducing and “minimizing” a user-chosen “cost functional quantifying the discrepancies between measurements and computations,” thus yielding results that are free of subjective user-interferences while generalizing and significantly extending the 4D-VAR data assimilation procedures. Incorporating correlations, including those between the imprecisely known model parameters and computed model responses, the 2nd/3rd-BERRU-PM also provides a quantitative metric, constructed from sensitivity and covariance matrices, for determining the degree of agreement among the various computational and experimental data while eliminating discrepant information. The mathematical framework of the 2nd/3rd-BERRU-PM formalism requires the inversion of a single matrix of size Nr × Nr, where Nr denotes the number of considered responses. In the overwhelming majority of practical situations, the number of responses is much less than the number of model parameters. Thus, the 2nd-BERRU-PM methodology overcomes the curse of dimensionality which affects the inversion of hessian matrices in the parameter space.

1. Introduction

The functional derivatives (also called “sensitivities”) of results (also called “responses”) 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. As is well known, even the approximate determination of the first-order sensitivities R / α i , i = 1, …, Nα of a model response R to Nα parameters αi using conventional finite-difference methods would already require Nα large-scale computations with altered parameter values, which is unfeasible for large-scale models comprising many parameters. The computation of higher-order sensitivities by conventional methods is limited in practice by the so-called “curse of dimensionality,” since the number of such large-scale computations increases exponentially with the order of the response sensitivities. For the exact computation of the first- and second-order response sensitivities to parameters, the “curse of dimensionality” has been overcome by the Second-Order Adjoint Sensitivity Analysis Methodology (2nd-ASAM) conceived and developed by Cacuci [1,2,3]. The unique capability of the 2nd-ASAM to compute comprehensively and efficiently the exact first-and second-order sensitivities of a response to parameters in a large-scale physical system has been demonstrated [4,5,6] by an application to a reactor physics system which comprises 21,976 first-order sensitivities and 482,944,576 second-order sensitivities.
Section 2 and Section 3 of this work present the Third-Order Adjoint Sensitivity Analysis Methodology (3rd-ASAM) for coupled forward and adjoint linear systems, which evidently extends and generalizes the 2nd-ASAM. The 3rd-ASAM aims at the efficient computation of the exact expressions of the 3rd-order functional derivatives (“sensitivities”) of a system response that depends on both the forward and adjoint state functions with respect to all of the parameter underlying the respective forward and adjoint systems. Such responses are often encountered when representing mathematically detector responses and reaction rates in reactor physics problems. The 3rd-ASAM will be applied to the reactor physics system analyzed in [4,5,6,7] to compute the exact magnitude of the 3rd-order sensitivities to the model parameters that were found in [4,5,6,7] to have unexpectedly large 2nd-order sensitivities. Furthermore, the 3rd-order sensitivities computed using the 3rd-ASAM are incorporated in the new formulas, presented in Section 4, for computing to 3rd-order the first four cumulants of the response distribution in the phase-space of model parameters. Section 5 presents a new mathematical formalism, which will be called the “Second/Third-Order Best-Estimated Results with Reduced Uncertainties Predictive Modeling (2nd/3rd-BERRU-PM).” Set in the joint phase-space of responses and parameters, the 2nd/3rd-BERRU-PM incorporates experimental and computational information, including the complete (as opposed to partial) vector of 1st-order response sensitivities, the complete hessian matrix of 2nd-order second-sensitivities and also the 3rd-order sensitivities, all computed using the 3rd-ASAM presented in Section 3. Thus, the 2nd/3rd-BERRU-PM extends the “BERRU Predictive Modeling” [7], thereby generalizing and including, as particular cases, similar formulas used in other fields e.g., [8,9,10]. The 2nd/3rd-BERRU-PM uses the maximum entropy (MaxEnt) principle [11] to eliminate the need for introducing and “minimizing” a user-chosen “cost functional quantifying the discrepancies between measurements and computations.” Incorporating correlations, including those between the imprecisely known model parameters and computed model responses, the 2nd/3rd-BERRU-PM also provides a quantitative metric, constructed from sensitivity and covariance matrices, for determining the degree of agreement among the various computational and experimental data and helping eliminate discrepant information. Conclusions regarding the significance of this work’s novel results in the quest to overcome the curse of dimensionality in sensitivity analysis, uncertainty quantification and predictive modeling are presented in Section 6.

2. Mathematical Description of the Physical System

A linear physical system is generally represented by means of Nu coupled operator equations of the form j = 1 N u L i j ( α ) φ j ( x ) = Q i ( α ) ,   i = 1 , , N u ,   x Ω x , in which the operators L i j ( α ) act linearly on the state functions φ j ( x ) . The above system of equations can be written in matrix form as follows:
L [ α ( x ) ] φ ( x ) = Q [ α ( x ) ]   ,   x Ω x
Matrices and vectors will be denoted using bold letters. Since the right-side of Equation (1) may contain distributions, the equality in this equation is considered to hold in the weak (“distributional”) sense. Similarly, all of the equalities that involve differential equations in this work will be considered to hold in the weak/distributional sense. All vectors in this work are considered to be column vectors, and transposition will be indicated by a dagger ( ) superscript. The vectors, matrices and operators appearing in Equation (1) are defined as follows:
  • α ( α 1 , , α N α ) denotes a N α -dimensional column vector whose components are the physical system’s imprecisely known parameters, which are subject to uncertainties; α E α N α , where E α denotes a subset of the N α -dimensional real vector space N α . The symbol “ ” will be used to denote “is defined as” or “is by definition equal to.” The vector α E α N α is considered to include any imprecisely known model parameters that may enter into defining the system’s boundary in the phase space of independent variables.
  • x ( x 1 , , x N x ) N x denotes the   N x -dimensional phase-space position vector, defined on a phase-space domain denoted as Ω x .
  • φ ( x ) [ φ 1 ( x ) , , φ N u ( x ) ] denotes a N u -dimensional column vector whose components represent the system’s dependent variables (also called “state functions”). In virtually all of the physical systems represented by Equation (1), the components φ i ( x ) ,   i = 1 , , N u are square-integrable functions and φ ( x ) H φ , where H φ is a Hilbert space endowed with an inner product that will be denoted as φ ( x ) , ψ ( x ) φ , ψ ( x ) H φ , and which is defined as follows:
    φ ( x ) , ψ ( x ) φ Ω x φ ( x ) · ψ ( x ) d x i = 1 N u Ω x φ i ( x ) ψ i ( x ) d x .
  • Q ( α ) [ Q 1 ( α ) , , Q N u ( α ) ] denotes a N u -dimensional column vector. The components of Q are operators acting (in general nonlinearly) on α and x ;
  • L ( α ) [ L 1 ( α ) , , L N u ( α ) ] denotes a N u -component column vector. The components of L ( α ) are operators acting linearly on φ and nonlinearly on α . When L ( α ) contains differential operators, a set of boundary and/or initial conditions which define the domain of L ( α ) must also be given. Since L ( α ) is considered to act linearly on φ ( x ) , the accompanying boundary and/or initial conditions must also be linear in φ ( x ) . Such linear boundary and/or initial conditions are represented in the following operator form:
    B F ( α ) φ ( x ) C F ( α ) = 0 ,   x Ω x .
In Equation (3), the operator B F ( α ) [ B i j ( α ) ] ;   i = 1 , , N B ;   j = 1 , , N u is a matrix comprising, as components, operators that act linearly on φ ( x ) and nonlinearly on α ; the quantity N B denotes the total number of boundary and initial conditions. The operator C F ( α ) [ C 1 ( α ) , , C N B ( α ) ] is a N B -dimensional vector comprising components that are operators acting, in general, nonlinearly on α . The subscript “F” in Equation (3) indicates boundary conditions associated with the “forward” system of equations.
In most practical situations the Hilbert space H φ is self-dual. The operator L ( α ) admits an adjoint (operator), which will be denoted as L + ( α ) , and which is defined through the following relation for an arbitrary vector ψ ( x ) H φ :
ψ ( x ) ,   L ( α )   φ ( x ) φ = L + ( α )   ψ ( x ) ,   φ ( x ) φ
In Equation (4), the formal adjoint operator L + ( α ) is the N u × N u matrix:
L + ( α ) [ L j i + ( α ) ] ,   i , j = 1 , , N u   ,
comprising elements L j i + ( α ) obtained by transposing the formal adjoints of the operators L i j ( α ) . Thus, the system adjoint to Equations (1) and (3) has the following general representation in operator form:
L + ( α ) ψ ( x ) = Q + ( α )   ,   x Ω x ,
B A ( α ) ψ ( x ) C A ( α ) = 0 ,   x Ω x .
The domain of L + ( α ) is determined by selecting the adjoint boundary and/or initial conditions represented in operator form in Equation (7), where the letter “A” indicates “adjoint” and the letter “B” indicates “boundary and/or initial conditions.” These adjoint boundary and/or initial conditions are selected so as to ensure that the boundary terms that arise when going from the left-side to the right side of Equation (4) vanish, in conjunction with the forward boundary conditions given in Equation (3). The source term Q + [ α ( x ) ] in Equation (6) is associated with the system’s response which, in this work, is considered to be a scalar-valued nonlinear functional of the adjoint and forward fluxes, which will be denoted as R [ φ ( x ) , ψ ( x ) ; α ] . Such responses are often encountered when representing mathematically detector responses and reaction rates in reactor physics problems, and can be generally represented in the following inner-product form:
R [ φ ( x ) ; ψ ( x ) ; α ] S ( φ ; ψ ; α ) , 1 φ
where S ( φ ; ψ ; α ) denotes a suitably differentiable function of its arguments. The nominal solution of Equations (1) and (3) is denoted as φ 0 ( x ) , and is obtained by solving these equations at the nominal parameter values α 0 . The superscript “zero” will henceforth be used to denote “nominal” or “expected” or “mean” values. Thus, the vectors φ 0 ( x ) and α 0 satisfy the following equations:
L ( α 0 ) φ 0 ( x ) = Q [ α 0 ( x ) ]   ,   x Ω x ,
B F ( α 0 ) φ 0 ( x ) C F ( α 0 ) = 0 ,   x Ω x .
Equations (9) and (10) represent the “base-case” or nominal state of the forward physical system. Similarly, the “base-case” or nominal state of the adjoint physical system is given by the following equations:
L + ( α 0 ) ψ 0 ( x ) = Q + ( α 0 )   ,   x Ω x ,
B A ( α 0 ) ψ 0 ( x ) C A ( α 0 ) = 0 ,   x Ω x .
The nominal value of the response, R [ φ 0 ( x ) , ψ 0 ( x ) ; α 0 ] , is determined by using the nominal parameter values α 0 , the nominal value of the forward function φ 0 ( x ) obtained by solving Equations (9) and (10), and the nominal value of the adjoint function obtained by solving Equations (11) and (12).

3. The Third-Order Adjoint Sensitivity Analysis Methodology (3rd-ASAM) for Coupled Linear Forward and Adjoint Systems: Another Step Towards Overcoming the Curse of Dimensionality in the Exact Computation of High-Order Response Sensitivities

The model parameters α i are imprecisely known quantities, so their actual values may differ from their nominal values by quantities denoted as δ α i α i α i 0 ,   i = 1 , , N α . Since the model parameters α and the state functions are related to each other through the forward and adjoint systems, it follows that variations δ α ( δ α 1 , , δ α N α ) in the model parameters will cause corresponding variations δ φ ( δ φ 1 , , δ φ N u ) ,   δ φ i φ i φ i 0 ,   i = 1 , , N u and δ ψ ( δ ψ 1 , , δ ψ N u ) ,   δ ψ i ψ i ψ i 0 ,   i = 1 , , N u in the forward and adjoint state functions. In turn, the variations δ α , δ φ , and δ ψ cause a response variation R ( φ 0 + δ φ ; ψ 0 + δ ψ ; α 0 + δ α ) around the nominal response value R [ φ 0 ( x ) , ψ 0 ( x ) ; α 0 ] .

3.1. The First-Level Adjoint System (1st-LASS) for Computing Exactly and Efficiently the First-Order Model Response Sensitivities to Parameters

The total first-order sensitivity of the response R [ φ ( x ) , ψ ( x ) ; α ] to variations δ α [ δ α 1 , , δ α N α ] in the model parameters is given by the definition of the Gateaux (G-) differential, denoted as δ R ( φ 0 ; ψ 0 ; α 0 ; δ φ ; δ ψ ; δ α ) of R [ φ ( x ) , ψ ( x ) ; α ] around the nominal values ( φ 0 ; ψ 0 ; α 0 ) . By definition, this G-differential is:
δ R ( φ 0 ; ψ 0 ; α 0 ; δ φ ; δ ψ ; δ α ) d d ε { R ( φ 0 + ε δ φ ; ψ 0 + ε δ ψ ; α 0 + ε δ α ) } ε = 0 = { δ R   ( φ 0 ; ψ 0 ; α 0 ;   δ α ) } d i r + { δ R   ( φ 0 ; ψ 0 ; α 0 ;   δ φ ; δ ψ ) } i n d ,
where the “direct-effect term” { δ R   ( φ 0 ; ψ 0 ; α 0 ;   δ α ) } d i r depends solely on the parameter variations δ α and is generally defined as follows:
{ δ R   ( φ ; ψ ; α ;   δ α ) } d i r [ S ( φ ; ψ ; α ) / α ]   δ α , 1 φ
while the “indirect-effect term” { δ R   ( φ 0 ; ψ 0 ; α 0 ;   δ φ ; δ ψ ) } i n d depends solely on the variations   δ φ and δ ψ in the forward and, respectively, adjoint functions, and is generally defined as follows:
{ δ R   ( φ ; ψ ; α ;   δ φ ; δ ψ ) } i n d [ S ( φ ; ψ ; α ) / φ ]   δ φ , 1 φ + [ S ( φ ; ψ ; α ) / ψ ]   δ ψ , 1 φ .
Since the nominal values of the forward and adjoint functions are known after having solved Equations (9) through (12), it follows that the direct-effect term { δ R   ( φ 0 ; ψ 0 ; α 0 ;   δ α ) } d i r can already be computed at this stage. The indirect-effect term { δ R   ( φ 0 ; ψ 0 ; α 0 ;   δ φ ; δ ψ ) } i n d can be computed only after having determined the functions   δ φ and δ ψ . These functions are obtained by solving the 1st-Level Forward Sensitivity System (1st-LFSS), which is obtained by G-differentiating the original forward and adjoint transport equations and respective boundary conditions given in Equations (1), (3), (6) and (7). Performing these differentiations yields the following 1st-LFSS:
{ ( L ( α ) 0 0 L + ( α ) ) } ( α 0 ) ( δ φ ( x ) δ ψ ( x ) ) = ( Q 1 ( 1 ) ( φ ; α ; δ α ) Q 2 ( 1 ) ( ψ ; α ; δ α ) ) ( φ 0 ; ψ 0 ; α 0 ) ,
together with the following boundary conditions:
{ B F ( α ) δ φ ( x ) } ( α 0 ) = { [ C F ( α ) B F ( α ) φ ( x ) ] / α } ( φ 0 ; α 0 ) δ α ,   x Ω x , { B A ( α ) δ ψ ( x ) } ( α 0 ) = { [ C A ( α ) B A ( α ) ψ ( x ) ] / α } ( ψ 0 ; α 0 ) δ α ,   x Ω x .
The source-terms Q 1 ( 1 ) ( φ ; α ; δ α ) and Q 2 ( 1 ) ( ψ ; α ; δ α ) in Equation (16) are defined as follows:
Q 1 ( 1 ) ( φ ; α ; δ α ) { [ Q [ α ( x ) ] L ( α ) φ ( x ) ] / α } δ α Q 2 ( 1 ) ( ψ ; α ; δ α ) { [ Q + [ α ( x ) ] L + ( α ) ψ ( x ) ] / α } δ α ,
Solving the 1st-LFSS defined by Equations (16) and (17) is computationally expensive, since the 1st-LFSS would need to be solved anew for every variation δ α i ,   i = 1 , , N α in the model parameters, as each such variation would affect the source term on the right-side of Equation (16). The computationally expensive evaluation of the indirect-effect term by using the 1st-LFSS can be avoided by expressing this indirect-effect term { δ R   ( φ 0 ; ψ 0 ; α 0 ;   δ φ ; δ ψ ) } i n d in terms of the solution of the 1st-Level Adjoint Sensitivity System (1st-LASS), which is constructed by implementing the following sequence of steps:
(i)
Consider two vector-valued functions u ( 1 ) ( x ) [ u 1 ( 1 ) ( x ) , u 2 ( 1 ) ( x ) ] and v ( 1 ) ( x ) [ v 1 ( 1 ) ( x ) , v 2 ( 1 ) ( x ) ] , each having two N u -dimensional vector-components defined as follows: u i ( 1 ) ( x ) [ u i , 1 ( 1 ) ( x ) , , u i , j ( 1 ) ( x ) , , u i , N u ( 1 ) ( x ) ] ,   i = 1 , 2 , and v i ( 1 ) ( x ) [ v i , 1 ( 1 ) ( x ) , , v i , j ( 1 ) ( x ) , , v i , N u ( 1 ) ( x ) ] ,   i = 1 , 2 , . The components of these vectors are assumed to be square-integrable functions.
(ii)
Introduce a Hilbert space, denoted as H ( 1 ) , endowed with the following inner product, denoted as u ( 1 ) ( x ) , v ( 1 ) ( x ) ( 1 ) , u ( 1 ) ( x ) H ( 1 ) , v ( 1 ) ( x ) H ( 1 ) , between the two functions defined in item (i), above:
u ( 1 ) ( x ) , v ( 1 ) ( x ) ( 1 ) i = 1 2 u i ( 1 ) ( x ) , v i ( 1 ) ( x ) φ = i = 1 2 j = 1 N u Ω x [ u i , j ( 1 ) ( x ) ] [ v i , j ( 1 ) ( x ) ] d x
(iii)
In the Hilbert space H ( 1 ) , form the inner product of Equation (16) with a yet undefined vector-valued function ψ ( 1 ) ( x ) [ ψ 1 ( 1 ) ( x ) , ψ 2 ( 1 ) ( x ) ] H ( 1 ) to obtain the following relation, evaluated at ( φ 0 ; ψ 0 ; α 0 ) , in which the superscript “zero” is omitted to simplify the notation:
( ψ 1 ( 1 ) ψ 2 ( 1 ) ) ,   ( L ( α ) 0 0 L + ( α ) ) ( δ φ δ ψ ) ( 1 ) = ( ψ 1 ( 1 ) ψ 2 ( 1 ) ) ,   ( Q 1 ( 1 ) ( φ ; α ; δ α ) Q 2 ( 1 ) ( ψ ; α ; δ α ) ) ( 1 )
(iv)
Use the definition of the adjoint operator in the Hilbert space H ( 1 ) to recast the left-side of Equation (20) as follows:
( ψ 1 ( 1 ) ψ 2 ( 1 ) ) ,   ( L ( α ) 0 0 L + ( α ) ) ( δ φ δ ψ ) ( 1 ) = ( δ φ δ ψ ) ,   ( L + ( α ) 0 0 L ( α ) ) ( ψ 1 ( 1 ) ψ 2 ( 1 ) ) ( 1 ) + P ( 1 ) [ δ φ ; δ ψ ; ψ 1 ( 1 ) , ψ 2 ( 1 ) ]
where the bilinear concomitant P ( 1 ) [ δ φ ; δ ψ ; ψ 1 ( 1 ) , ψ 2 ( 1 ) ] is defined on the phase-space boundary x Ω x . The superscript “zero” denoting nominal values for the quantities ( φ 0 ; ψ 0 ; α 0 ) was also omitted in Equations (20) and (21), in order to simplify the notation. Omitting henceforth the superscript “zero” denoting nominal values for the quantities ( φ 0 ; ψ 0 ; α 0 ) should not cause any loss of clarity, since all quantities are to be evaluated/computed using the respective nominal values of the model parameters and using the nominal values of the forward and adjoint functions evaluated at preceding stages/steps at nominal parameter values.
(v)
Identify the term on the left-side of Equation (21) with the indirect effect term defined in Equation (15), i.e., require that
( ψ 1 ( 1 ) ψ 2 ( 1 ) ) ,   ( L ( α ) 0 0 L + ( α ) ) ( δ φ δ ψ ) ( 1 ) = [ S ( φ ; ψ ; α ) / φ ]   δ φ , 1 φ + [ S ( φ ; ψ ; α ) / ψ ]   δ ψ , 1 φ = ( ψ 1 ( 1 ) ψ 2 ( 1 ) ) ,   ( Q 1 ( 1 ) ( φ ; α ; δ α ) Q 2 ( 1 ) ( ψ ; α ; δ α ) ) ( 1 )
and use Equation (21) in conjunction with the boundary conditions given in Equation (17) on the above relation to construct the following 1st-Level Adjoint Sensitivity System (1st-LASS):
( L + ( α ) 0 0 L ( α ) ) ( ψ 1 ( 1 ) ( x ) ψ 2 ( 1 ) ( x ) ) = ( [ S ( φ ; ψ ; α ) / φ ] [ S ( φ ; ψ ; α ) / ψ ] ) .
(vi)
The boundary conditions given in Equation (17) are now implemented in Equation (21), thereby reducing by half the number of unknown boundary-values in the bilinear concomitant P ( 1 ) [ δ φ ; δ ψ ; ψ 1 ( 1 ) , ψ 2 ( 1 ) ] . The boundary conditions for the adjoint functions ψ 1 ( 1 ) ( x ) and ψ 2 ( 1 ) ( x ) are chosen next so as to eliminate the remaining unknown boundary-values of the functions δ φ and δ ψ while ensuring that Equation (22) is well posed. The boundary conditions thus chosen for the adjoint functions ψ 1 ( 1 ) ( x ) and ψ 2 ( 1 ) ( x ) can be represented in operator form as follows:
B A ( 1 ) [ φ ( x ) ; ψ ( x ) ; ψ 1 ( 1 ) ( x ) , ψ 2 ( 1 ) ( x ) ; α ] = 0 ,   x Ω x .
(vii)
In most cases, the above choice of boundary conditions for the 1st-level adjoint function ψ ( 1 ) ( x ) will cause the bilinear concomitant P ( 1 ) [ δ φ , δ ψ , ψ 1 ( 1 ) , ψ 2 ( 1 ) ] in Equation (21) to vanish. When the boundary conditions for the original system are non-homogeneous, however, the bilinear concomitant P ( 1 ) [ δ φ , δ ψ , ψ 1 ( 1 ) , ψ 2 ( 1 ) ] may not vanish. Even when it does not vanish, however, this bilinear concomitant will be reduced to a quantity, denoted here as P ^ ( 1 ) [ φ ; ψ ; ψ 1 ( 1 ) , ψ 2 ( 1 ) ; α ; δ α ] , which will contain only known values of its arguments.
(viii)
Use the 1st-LASS defined by Equations (22) and (23) together with Equations (20) and (21) to obtain the following expression for the indirect-effect term defined in Equation (15), in terms of the adjoint functions ψ 1 ( 1 ) ( x ) and ψ 2 ( 1 ) ( x ) :
{ δ R   ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ; δ α ) } i n d = ( ψ 1 ( 1 ) ψ 2 ( 1 ) ) ,   ( Q 1 ( 1 ) ( φ ; α ; δ α ) Q 2 ( 1 ) ( ψ ; α ; δ α ) ) ( 1 ) P ^ ( 1 ) [ φ ; ψ ; ψ 1 ( 1 ) , ψ 2 ( 1 ) ; α ; δ α ] .
As indicated in Equation (22), the function ψ 2 ( 1 ) ( x ) is computed by solving the original forward equation with the source S ( φ ; ψ ; α ) / ψ while the function ψ 1 ( 1 ) ( x ) is obtained by solving the original adjoint equation with the source [ S ( φ ; ψ ; α ) / φ ] . Thus, after the 1st-LASS is solved to determine these two adjoint functions, the “indirect-effect term” is computed efficiently and exactly by simply performing the integrations (“quadratures”) indicated by the inner-product on the right-side of Equation (24).
Replacing Equations (24) and (14) in Equation (13) eliminates the appearance of the functions δ φ ( x ) , δ ψ ( x ) in the resulting expression. Consequently, the total 1st-order response sensitivity can be expressed in terms of the adjoint functions ψ 1 ( 1 ) ( x ) and ψ 2 ( 1 ) ( x ) as follows:
δ R   ( φ ; ψ ; ψ ( 1 ) ; α ; δ α ) = S ( φ ; ψ ; α ) α   δ α ,   1 φ + ψ 1 ( 1 ) ,   Q 1 ( 1 ) ( φ ; α ; δ α ) φ + ψ 2 ( 1 ) ,   Q 2 ( 1 ) ( ψ ; α ; δ α ) φ P ^ ( 1 ) [ φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ; δ α ] i = 1 N α { R ( φ ; ψ ;   ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) α i } δ α i .
All of the quantities shown in Equation (25) are to be evaluated at the nominal values ( φ 0 ; ψ 0 ;   ψ 1 ( 1 ) , 0 ; ψ 2 ( 1 ) , 0 ; α 0 ) . The partial 1st-order response sensitivities, denoted in Equation (25) as R ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) / α i ,   i = 1 , , N α , of the response R ( φ ; ψ ; α ) to a generic parameter α i are obtained by identifying the quantities that multiply the various parameter variations δ α i in Equation (25).

3.2. The Second-Level Adjoint System (2nd-LASS) for Computing Exactly and Efficiently the Second-Order Model Response Sensitivities to Parameters

The second order sensitivities of the response R ( φ ; ψ ; α ) with respect to the parameters α i ,   i = 1 , , N α , are obtained by determining the first-order G-differentials of the 1st-order sensitivities R ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) / α i . For this purpose, it is convenient to use the notation R i ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) R ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) / α i . The first-order G-differential of R i ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) is obtained, by definition, as follows:
{ δ R i ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ; δ φ ; δ ψ ; δ ψ 1 ( 1 ) ; δ ψ 2 ( 1 ) ; δ α ) } ( φ 0 ; ψ 0 ; ψ 1 ( 1 ) , 0 ; ψ 2 ( 1 ) , 0 ; α 0 ) d d ε { R i ( 1 ) ( φ 0 + ε δ φ ; ψ 0 + ε δ ψ ; ψ 1 ( 1 ) , 0 + ε δ ψ 1 ( 1 ) ; ψ 2 ( 1 ) , 0 + ε δ ψ 2 ( 1 ) ; α 0 + ε δ α ) } ε = 0 = { R i ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) α } ( φ 0 ; ψ 0 ; ψ ( 1 ) , 0 ; α 0 ) δ α + { δ R i ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ; δ φ ; δ ψ ; δ ψ 1 ( 1 ) ; δ ψ 2 ( 1 ) ) } i n d ,
where the “indirect-effect term” { δ R i ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ; δ φ ; δ ψ ; δ ψ 1 ( 1 ) ; δ ψ 2 ( 1 ) ) } i n d depends solely on the variations   δ φ , δ ψ , δ ψ 1 ( 1 ) and δ ψ 2 ( 1 ) in the forward and, respectively, adjoint functions and is defined as follows:
{ δ R i ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ; δ φ ; δ ψ ; δ ψ 1 ( 1 ) ; δ ψ 2 ( 1 ) ) } i n d { [ R i ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) φ ] ,   δ φ φ + [ R i ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) ψ ] ,   δ ψ φ + [ R i ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) ψ 1 ( 1 ) ] ,   δ ψ 1 ( 1 ) φ + [ R i ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) ψ 2 ( 1 ) ] ,   δ ψ 2 ( 1 ) φ } ( φ 0 ; ψ 0 ; ψ 1 ( 1 ) , 0 ; ψ 2 ( 1 ) , 0 ; α 0 ) .
The indirect-effect term defined in Equation (27) can be computed only after having obtained both the solution the functions   δ φ , δ ψ , δ ψ 1 ( 1 ) and δ ψ 2 ( 1 ) . The functions   δ φ and δ ψ are the solutions of the 1st-LFSS defined by Equations (16). Furthermore, the functions δ ψ 1 ( 1 ) and δ ψ 2 ( 1 ) are the solutions of the following system of equations obtained by G-differentiating the 1st-LASS:
L + ( α ) δ ψ 1 ( 1 ) ( x ) 2 S ( φ ; ψ ; α ) φ φ δ φ ( x ) 2 S ( φ ; ψ ; α ) φ   ψ δ ψ ( x ) = Q 1 ( 2 ) ( φ ; ψ ; ψ 1 ( 1 ) ; α ; δ α ) L ( α ) δ ψ 2 ( 1 ) ( x ) 2 S ( φ ; ψ ; α ) ψ φ δ φ ( x ) 2 S ( φ ; ψ ; α ) ψ ψ δ ψ ( x ) = Q 2 ( 2 ) ( φ ; ψ ; ψ 2 ( 1 ) ; α ; δ α ) }
δ B A ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) B A ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) φ δ φ + B A ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) ψ δ ψ + B A ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) ψ 1 ( 1 ) δ ψ 1 ( 1 ) + B A ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) ψ 2 ( 1 ) δ ψ 2 ( 1 ) + B A ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) α δ α = 0 ,   x Ω x .
where:
Q 1 ( 2 ) ( φ ; ψ ; ψ 1 ( 1 ) ; α ; δ α ) { 2 S ( φ ; ψ ; α ) φ α [ L + ( α ) ψ 1 ( 1 ) ( x ) ] α } δ α ,
Q 2 ( 2 ) ( φ ; ψ ; ψ 2 ( 1 ) ; α ; δ α ) { 2 S ( φ ; ψ ; α ) ψ α [ L ( α ) ψ 2 ( 1 ) ( x ) ] α } δ α .
The system comprising Equations (16) and (28), together with the boundary conditions provided in Equations (17) and (29) constitutes the 2nd-Level Forward Sensitivity System (2nd-LFSS). Since the source-terms of the 2nd-LFSS depend on the parameters variations δ α i , it follows that the determination of the functions δ ψ 1 ( 1 ) and δ ψ 2 ( 1 ) is at least as expensive computationally as determining the functions   δ φ and δ ψ by solving the 1st-LFSS. To avoid the need for solving the 2nd-LFSS, the indirect-effect term defined in Equation (27) will be expressed in terms of a 2nd-Level Adjoint Sensitivity System (2nd-LASS), which will be constructed by following the general principles introduced by Cacuci [1,2,3], comprising the following sequence of steps:
(i) Define a Hilbert space, denoted as H ( 2 ) , having vector-valued elements of the form the u ( 2 ) ( x ) [ u 1 ( 2 ) ( x ) , u 2 ( 2 ) ( x ) , u 3 ( 2 ) ( x ) , u 4 ( 2 ) ( x ) ] H ( 2 ) , with components that are N u -dimensional vectors of the form u i ( 2 ) ( x ) [ u i , 1 ( 2 ) ( x ) , , u i , j ( 2 ) ( x ) , , u i , N u ( 2 ) ( x ) ] ,   i = 1 , , 4 , with square-integrable components u i , j ( 2 ) ( x ) . In H ( 2 ) , define the inner-product, denoted as u ( 2 ) ( x ) , v ( 2 ) ( x ) ( 2 ) , of two functions u ( 2 ) ( x ) H ( 2 ) and v ( 2 ) ( x ) H ( 2 ) as follows:
u ( 2 ) ( x ) , v ( 2 ) ( x ) ( 2 ) i = 1 4 u i ( 2 ) ( x ) , v i ( 2 ) ( x ) φ = i = 1 4 j = 1 N u Ω x u i , j ( 2 ) ( x ) v i , j ( 2 ) ( x )   d x .
Using the definition provided in Equation (32), construct the inner product of a vector ψ i ( 2 ) ( x ) [ ψ 1 , i ( 2 ) ( x ) , ψ 2 , i ( 2 ) ( x ) , ψ 3 , i ( 2 ) ( x ) , ψ 4 , i ( 2 ) ( x ) ] H ( 2 ) with Equations (16) and (28) to obtain the following relation:
( ψ 1 , i ( 2 ) ( x ) ψ 2 , i ( 2 ) ( x ) ψ 3 , i ( 2 ) ( x ) ψ 4 , i ( 2 ) ( x ) ) ,   ( L ( α ) 0 0 0 0 L + ( α ) 0 0 2 S ( φ ; ψ ; α ) φ φ 2 S ( φ ; ψ ; α ) ψ φ L + ( α ) 0 2 S ( φ ; ψ ; α ) φ   ψ 2 S ( φ ; ψ ; α ) ψ ψ 0 L ( α ) ) ( δ φ ( x ) δ ψ ( x ) δ ψ 1 ( 1 ) ( x ) δ ψ 2 ( 1 ) ( x ) ) ( 2 ) = = ( ψ 1 , i ( 2 ) ( x ) ψ 2 , i ( 2 ) ( x ) ψ 3 , i ( 2 ) ( x ) ψ 4 , i ( 2 ) ( x ) ) ,   ( Q 1 ( 1 ) ( φ ; α ; δ α ) Q 2 ( 1 ) ( ψ ; α ; δ α ) Q 1 ( 2 ) ( φ ; ψ ; ψ 1 ( 1 ) ; α ; δ α ) Q 2 ( 2 ) ( φ ; ψ ; ψ 2 ( 1 ) ; α ; δ α ) ) ( 2 ) .
(ii)
Use the definition of the adjoint operator in the Hilbert space H ( 2 ) to recast the left-side of Equation (33) as follows:
( ψ 1 , i ( 2 ) ( x ) ψ 2 , i ( 2 ) ( x ) ψ 3 , i ( 2 ) ( x ) ψ 4 , i ( 2 ) ( x ) ) ,   ( L ( α ) 0 0 0 0 L + ( α ) 0 0 2 S ( φ ; ψ ; α ) φ φ 2 S ( φ ; ψ ; α ) ψ φ L + ( α ) 0 2 S ( φ ; ψ ; α ) φ   ψ 2 S ( φ ; ψ ; α ) ψ ψ 0 L ( α ) ) ( δ φ ( x ) δ ψ ( x ) δ ψ 1 ( 1 ) ( x ) δ ψ 2 ( 1 ) ( x ) ) ( 2 ) = = ( δ φ ( x ) δ ψ ( x ) δ ψ 1 ( 1 ) ( x ) δ ψ 2 ( 1 ) ( x ) ) ,   ( L + ( α ) 0 2 S ( φ ; ψ ; α ) φ φ 2 S ( φ ; ψ ; α ) φ   ψ 0 L ( α ) 2 S ( φ ; ψ ; α ) ψ φ 2 S ( φ ; ψ ; α ) ψ ψ 0 0 L ( α ) 0 0 0 0 L + ( α ) ) ( ψ 1 , i ( 2 ) ( x ) ψ 2 , i ( 2 ) ( x ) ψ 3 , i ( 2 ) ( x ) ψ 4 , i ( 2 ) ( x ) ) ( 2 ) + P ( 2 ) [ δ φ ; δ ψ ; δ ψ 1 ( 1 ) , δ ψ 2 ( 1 ) ; ψ 1 , i ( 2 ) ( x ) , ψ 2 , i ( 2 ) ( x ) , ψ 3 , i ( 2 ) ( x ) , ψ 4 , i ( 2 ) ( x ) ]   ,
where the bilinear concomitant P ( 2 ) [ δ φ ; δ ψ ; δ ψ 1 ( 1 ) , δ ψ 2 ( 1 ) ; ψ 1 , i ( 2 ) ( x ) , ψ 2 , i ( 2 ) ( x ) , ψ 3 , i ( 2 ) ( x ) , ψ 4 , i ( 2 ) ( x ) ] is defined on the phase-space boundary x Ω x .
(iii)
Identify the first term on the right-side of Equation (34) with the indirect-effect term defined in Equation (27) by requiring that the following system of equations be satisfied for i = 1 , N α :
  L + ( α ) ψ 1 , i ( 2 ) ( x ) 2 S ( φ ; ψ ; α ) φ φ ψ 3 , i ( 2 ) ( x ) 2 S ( φ ; ψ ; α ) φ   ψ ψ 4 , i ( 2 ) ( x ) = [ R i ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) / φ ] L ( α ) ψ 2 , i ( 2 ) ( x ) 2 S ( φ ; ψ ; α ) ψ φ ψ 3 , i ( 2 ) ( x ) 2 S ( φ ; ψ ; α ) ψ ψ ψ 4 , i ( 2 ) ( x ) = [ R i ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) / ψ ] L ( α ) ψ 3 , i ( 2 ) ( x ) = [ R i ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) / ψ 1 ( 1 ) ] L + ( α ) ψ 4 , i ( 2 ) ( x ) = [ R i ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) / ψ 2 ( 1 ) ]
(iv)
The boundary conditions given in Equations (17) and (29) are now implemented in Equation (34), thereby reducing by half the number of unknown boundary-values in the bilinear concomitant P ( 2 ) [ δ φ ; δ ψ ; δ ψ 1 ( 1 ) , δ ψ 2 ( 1 ) ; ψ 1 , i ( 2 ) ( x ) , ψ 2 , i ( 2 ) ( x ) , ψ 3 , i ( 2 ) ( x ) , ψ 4 , i ( 2 ) ( x ) ] . The boundary conditions for the 2nd-level adjoint functions ψ 1 , i ( 2 ) ( x ) , ψ 2 , i ( 2 ) ( x ) , ψ 3 , i ( 2 ) ( x ) , ψ 4 , i ( 2 ) ( x ) are now chosen so as to eliminate the remaining unknown boundary-values of the functions δ φ , δ ψ , δ ψ 1 ( 1 ) and δ ψ 2 ( 1 ) while ensuring that Equation (35) is well posed. The boundary conditions thus chosen for the adjoint functions ψ 1 ( 1 ) ( x ) and ψ 2 ( 1 ) ( x ) can be represented in operator form as follows:
B A ( 2 ) [ φ ( x ) ; ψ ( x ) ; ψ 1 ( 1 ) ( x ) ; ψ 2 ( 1 ) ( x ) ; ψ 1 , i ( 2 ) ( x ) ; ψ 2 , i ( 2 ) ( x ) ; ψ 3 , i ( 2 ) ( x ) ; ψ 4 , i ( 2 ) ( x ) ; α ] = 0 ,   x Ω x .
In most cases, the above choice of boundary conditions for the 1st-level adjoint function ψ ( 1 ) ( x ) will cause the bilinear concomitant P ( 2 ) [ δ φ ; δ ψ ; δ ψ 1 ( 1 ) , δ ψ 2 ( 1 ) ; ψ 1 , i ( 2 ) ( x ) , ψ 2 , i ( 2 ) ( x ) , ψ 3 , i ( 2 ) ( x ) , ψ 4 , i ( 2 ) ( x ) ] in Equation (34) to vanish. Even when it does not vanish, however, this bilinear concomitant will be reduced to a quantity, denoted here as P ^ ( 2 ) [ φ ; ψ ; ψ 1 ( 1 ) , ψ 2 ( 1 ) ; ψ 1 , i ( 2 ) ( x ) , ψ 2 , i ( 2 ) ( x ) , ψ 3 , i ( 2 ) ( x ) , ψ 4 , i ( 2 ) ( x ) ; α ; δ α ] , which will contain only known values of its arguments.
The system of equations comprising Equations (35) and (36) will be called the 2nd-Level Adjoint Sensitivity System (2nd-LASS) for the 2nd-level adjoint functions ψ 1 , j ( 2 ) ( x ) , ψ 2 , j ( 2 ) ( x ) , ψ 3 , j ( 2 ) ( x ) , ψ 4 , j ( 2 ) ( x ) , i = 1 , N α . These 2nd-level adjoint functions are obtained by solving the 2nd-LASS successively by using two “forward” and two “adjoint” computations, for each of the imprecisely known scalar model parameters.
(v)  Use the 2nd-LASS defined by Equations (35) and (36) together with Equations (33)–(35) to obtain the following expression for the indirect-effect term defined in Equation (27), in terms of the 2nd-level adjoint functions ψ 1 , j ( 2 ) ( x ) , ψ 2 , j ( 2 ) ( x ) , ψ 3 , j ( 2 ) ( x ) , ψ 4 , j ( 2 ) ( x )
{ δ R i ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) , ψ 2 ( 1 ) ; ψ 1 , i ( 2 ) ( x ) , ψ 2 , i ( 2 ) ( x ) , ψ 3 , i ( 2 ) ( x ) , ψ 4 , i ( 2 ) ( x ) ; α ; δ α ) } i n d = ( ψ 1 , i ( 2 ) ( x ) ψ 2 , i ( 2 ) ( x ) ψ 3 , i ( 2 ) ( x ) ψ 4 , i ( 2 ) ( x ) ) ,   ( Q 1 ( 1 ) ( φ ; α ; δ α ) Q 2 ( 1 ) ( ψ ; α ; δ α ) Q 1 ( 2 ) ( φ ; ψ ; ψ 1 ( 1 ) ; α ; δ α ) Q 2 ( 2 ) ( φ ; ψ ; ψ 2 ( 1 ) ; α ; δ α ) ) ( 2 ) P ^ ( 2 ) [ φ ; ψ ; ψ 1 ( 1 ) , ψ 2 ( 1 ) ; ψ 1 , i ( 2 ) ( x ) , ψ 2 , i ( 2 ) ( x ) , ψ 3 , i ( 2 ) ( x ) , ψ 4 , i ( 2 ) ( x ) ; α ; δ α ] .
As Equation (37) indicates, the indirect-effect term can be computed speedily by quadratures once the 2nd-level adjoint functions ψ 1 , j ( 2 ) ( x ) , ψ 2 , j ( 2 ) ( x ) , ψ 3 , j ( 2 ) ( x ) , ψ 4 , j ( 2 ) ( x ) become available.
(vi)  Replace Equation (37) in Equation (26) to obtain the following expression for the total 2nd-order response sensitivity to model parameters:
δ R i ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ; δ φ ; δ ψ ; δ ψ 1 ( 1 ) ; δ ψ 2 ( 1 ) ; δ α ) = j = 1 N α R i j ( 2 ) [ φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; ψ 1 , i ( 2 ) ( x ) , ψ 2 , i ( 2 ) ( x ) , ψ 3 , i ( 2 ) ( x ) , ψ 4 , i ( 2 ) ( x ) ; α ] δ α j   ,   i = 1 , , N α .
where R i j ( 2 ) [ φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; ψ 1 , i ( 2 ) ( x ) , ψ 2 , i ( 2 ) ( x ) , ψ 3 , i ( 2 ) ( x ) , ψ 4 , i ( 2 ) ( x ) ; α ] 2 R / α i α j ;   i , j = 1 , , N α denotes the 2nd-order partial sensitivity of the response to the model parameters and is defined as follows:
R i j ( 2 ) [ φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; ψ 1 , i ( 2 ) ( x ) , ψ 2 , i ( 2 ) ( x ) , ψ 3 , i ( 2 ) ( x ) , ψ 4 , i ( 2 ) ( x ) ; α ] R i ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) α j + ψ 1 , i ( 2 ) ( x ) ,   [ Q [ α ( x ) ] L ( α ) φ ( x ) ] / α j φ + ψ 2 , i ( 2 ) ( x ) ,   [ Q + [ α ( x ) ] L + ( α ) ψ ( x ) ] / α j φ + ψ 3 , i ( 2 ) ( x ) ,   2 S ( φ ; ψ ; α ) φ α j [ L + ( α ) ψ 1 ( 1 ) ( x ) ] α j φ + ψ 4 , i ( 2 ) ( x ) ,   2 S ( φ ; ψ ; α ) ψ α j [ L ( α ) ψ 2 ( 1 ) ( x ) ] α j φ P ^ ( 2 ) [ φ ; ψ ; ψ 1 ( 1 ) , ψ 2 ( 1 ) ; ψ 1 , i ( 2 ) ( x ) , ψ 2 , i ( 2 ) ( x ) , ψ 3 , i ( 2 ) ( x ) , ψ 4 , i ( 2 ) ( x ) ; α ] α j ,   i = 1 , , N α ;   j = 1 , , i .
Note that the 2nd-LASS is independent of parameter variations δ α . Thus, the exact computation of all of the partial second-order sensitivities, R i j ( 2 ) [ φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; ψ 1 , i ( 2 ) ( x ) , ψ 2 , i ( 2 ) ( x ) , ψ 3 , i ( 2 ) ( x ) , ψ 4 , i ( 2 ) ( x ) ; α ] ,   i ,   j = 1 , , N α , requires at most N α large-scale (adjoint) computations using the 2nd-LASS, rather than O ( N α 2 ) large-scale computations as would be required by forward methods. It is also important to note that by solving the 2nd-LASS N α -times, the “off-diagonal” 2nd-order mixed sensitivities R i j ( 2 ) [ φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; ψ 1 , i ( 2 ) ( x ) , ψ 2 , i ( 2 ) ( x ) , ψ 3 , i ( 2 ) ( x ) , ψ 4 , i ( 2 ) ( x ) ; α ] will be computed twice, in two different ways (using distinct 2nd-level adjoint functions), thereby providing an independent intrinsic (numerical) verification that the 1st- and 2nd-order sensitivities are computed accurately. In practice, it is useful to prioritize the computation of the 2nd-order sensitivities by using the rankings of the relative magnitudes of the 1st-order sensitivities as a “priority indicator”: the larger the magnitude of the relative 1st-order sensitivity, the higher the priority for computing the corresponding 2nd-order sensitivities. Also, since vanishing 1st-order sensitivities may indicate critical points of the response in the phase-space of model parameters, it is also of interest to compute the 2nd-order sensitivities that correspond to vanishing 1st-order sensitivities. Thus, only the 2nd-order partial sensitivities of the response R [ φ ( x ) , ψ ( x ) ; α ] which are deemed important will need to be computed. Information provided by the 1st-order sensitivities might indicate which 2nd-order sensitivities could be neglected.

3.3. The Third-Level Adjoint System (3rd-LASS) for Computing Exactly and Efficiently the Third-Order Model Response Sensitivities to Parameters

The third-order sensitivities of the response R ( φ ; ψ ; α ) with respect to the model parameters α i ,   i = 1 , , N α , are obtained by determining the first-order G-differential, δ R i j ( 2 ) [ φ 0 ; ψ 0 ; ψ 1 ( 1 ) , 0 ; ψ 2 ( 1 ) , 0 ; ψ 1 , i ( 2 ) , 0 ( x ) , ψ 2 , i ( 2 ) , 0 ( x ) , ψ 3 , i ( 2 ) , 0 ( x ) , ψ 4 , i ( 2 ) , 0 ( x ) ; α 0 ] , of the 2nd-order sensitivities R i j ( 2 ) [ φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; ψ 1 , i ( 2 ) ( x ) , ψ 2 , i ( 2 ) ( x ) , ψ 3 , i ( 2 ) ( x ) , ψ 4 , i ( 2 ) ( x ) ; α ] , computed in Section 3.2, which is given by the following expression:
δ R i j ( 2 ) [ φ 0 ( x ) ; ψ 0 ( x ) ; ψ 1 ( 1 ) , 0 ( x ) ; ψ 2 ( 1 ) , 0 ( x ) ; ψ 1 , i ( 2 ) , 0 ( x ) , ψ 2 , i ( 2 ) , 0 ( x ) , ψ 3 , i ( 2 ) , 0 ( x ) , ψ 4 , i ( 2 ) , 0 ( x ) ; α 0 ; δ φ ( x ) ; δ ψ ( x ) ; δ ψ 1 ( 1 ) ( x ) ; δ ψ 2 ( 1 ) ( x ) ; δ ψ 1 , i ( 2 ) ( x ) , δ ψ 2 , i ( 1 ) , δ ψ 3 , i ( 1 ) , 0 , δ ψ 4 , i ( 1 ) , 0 ; δ α ] d d ε { R i j ( 2 ) [ φ 0 ( x ) + ε δ φ ( x ) ; ψ 0 ( x ) + ε δ ψ ( x ) ; ψ 1 ( 1 ) , 0 ( x ) + ε δ ψ 1 ( 1 ) ( x ) ; ψ 2 ( 1 ) , 0 ( x ) + ε δ ψ 2 ( 1 ) ( x ) ; ψ 1 , i ( 2 ) ( x ) + ε δ ψ 1 , i ( 2 ) ( x ) , ψ 2 , i ( 1 ) , 0 ( x ) + ε δ ψ 2 , i ( 1 ) , ψ 3 , i ( 2 ) ( x ) + ε δ ψ 3 , i ( 1 ) , 0 , ψ 4 , i ( 2 ) ( x ) + ε δ ψ 4 , i ( 1 ) , 0 ; α 0 + ε δ α ] } ε = 0 = { R i j ( 2 ) ( φ ; ψ ; ; α ) α } ( φ 0 ; ψ 0 ; ; α 0 ) δ α + { δ R i j ( 2 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; ψ 1 , i ( 2 ) , , ψ 4 , i ( 2 ) ; α ; δ φ ; δ ψ ; δ ψ 1 ( 1 ) ; δ ψ 2 ( 1 ) ; δ ψ 1 , i ( 2 ) , , δ ψ 4 , i ( 2 ) ) } i n d ,   i = 1 , , N α ;   j = 1 , , i .
The quantity { δ R i j ( 2 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; ψ 1 , i ( 2 ) , , ψ 4 , i ( 2 ) ; α ; δ φ ; δ ψ ; δ ψ 1 ( 1 ) ; δ ψ 2 ( 1 ) ; δ ψ 1 , i ( 2 ) , , δ ψ 4 , i ( 2 ) ) } i n d denotes the “indirect-effect term”, which depends solely on the variations ( δ φ ; δ ψ ; δ ψ 1 ( 1 ) ; δ ψ 2 ( 1 ) ; δ ψ 1 , i ( 2 ) , , δ ψ 4 , i ( 2 ) ) in the forward and adjoint functions, respectively, and is defined as follows:
{ δ R i j ( 2 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; ψ 1 , i ( 2 ) , , ψ 4 , i ( 2 ) ; α ; δ φ ; δ ψ ; δ ψ 1 ( 1 ) ; δ ψ 2 ( 1 ) ; δ ψ 1 , i ( 2 ) , , δ ψ 4 , i ( 2 ) ) } i n d { [ R i , j ( 2 ) ( φ α ) φ ] ,   δ φ φ + [ R i , j ( 2 ) ( φ ; ψ α ) ψ ] ,   δ ψ φ + [ R i , j ( 2 ) ( ; ψ 1 ( 1 ) α ) ψ 1 ( 1 ) ] ,   δ ψ 1 ( 1 ) φ + [ R i , j ( 2 ) ( ; ψ 2 ( 1 ) ) ψ 2 ( 1 ) ] ,   δ ψ 2 ( 1 ) φ + [ R i , j ( 2 ) ( ; ψ 1 , i ( 2 ) ) ψ 1 , i ( 2 ) ] ,   δ ψ 1 , i ( 2 ) φ + [ R i , j ( 2 ) ( ; ψ 2 , i ( 2 ) ) ψ 1 , i ( 2 ) ] ,   δ ψ 1 , i ( 2 ) φ + [ R i , j ( 2 ) ( ; ψ 3 , i ( 2 ) ) ψ 3 , i ( 2 ) ] ,   δ ψ 3 , i ( 2 ) φ + [ R i , j ( 2 ) ( ; ψ 4 , i ( 2 ) ) ψ 4 , i ( 2 ) ] ,   δ ψ 4 , i ( 2 ) φ } ( φ 0 ; ; α 0 ) ,   i = 1 , , N α ;   j = 1 , , i .
The indirect-effect term defined in Equation (41) can be computed only after having computed the functions δ ψ 1 , i ( 2 ) ,   δ ψ 2 , i ( 2 ) ,   δ ψ 3 , i ( 2 ) ,   δ ψ 4 , i ( 2 ) , in addition to the functions   δ φ , δ ψ , δ ψ 1 ( 1 ) and δ ψ 2 ( 1 ) . Altogether, the functions   δ φ , δ ψ , δ ψ 1 ( 1 ) , δ ψ 2 ( 1 ) δ ψ 1 , i ( 2 ) ,   δ ψ 2 , i ( 2 ) ,   δ ψ 3 , i ( 2 ) ,   δ ψ 4 , i ( 2 ) , are the solutions of Equation (28), augmented by the solutions of the system of equations obtained by G-differentiating the 2nd-LASS, which can be written in matrix form as follows:
( L ( α ) 0 0 0 0 0 0 0 0 L + ( α ) 0 0 0 0 0 0 F 31 ( 2 ) F 32 ( 2 ) L + ( α ) 0 0 0 0 0 F 41 ( 2 ) F 42 ( 2 ) 0 L ( α ) 0 0 0 0 F 51 ( 3 ) F 52 ( 3 ) F 53 ( 3 ) F 54 ( 3 ) L + ( α ) 0 F 57 ( 3 ) F 58 ( 3 ) F 61 ( 3 ) F 62 ( 3 ) F 63 ( 3 ) F 64 ( 3 ) 0 L ( α ) F 67 ( 3 ) F 68 ( 3 ) F 71 ( 3 ) F 72 ( 3 ) F 73 ( 3 ) F 74 ( 3 ) 0 0 L ( α ) 0 F 81 ( 3 ) F 82 ( 3 )   F 83 ( 3 ) F 84 ( 3 ) 0 0 0 L + ( α ) ) ( δ φ ( x ) δ ψ ( x ) δ ψ 1 ( 1 ) ( x ) δ ψ 2 ( 1 ) ( x ) δ ψ 1 , i ( 2 ) ( x ) δ ψ 2 , i ( 2 ) ( x ) δ ψ 3 , i ( 2 ) ( x ) δ ψ 4 , i ( 2 ) ( x ) ) = ( Q 1 ( 1 ) ( φ ; α ; δ α ) Q 2 ( 1 ) ( ψ ; α ; δ α ) Q 1 ( 2 ) ( φ ; ; δ α ) Q 2 ( 2 ) ( φ ; ; δ α ) Q 1 , i ( 3 ) ( φ ; ; δ α ) Q 2 , i ( 3 ) ( φ ; ; δ α ) Q 3 , i ( 3 ) ( φ ; ; δ α ) Q 4 , i ( 3 ) ( φ ; ; δ α ) ) ,
where:
F 31 ( 2 ) 2 S ( φ ; ψ ; α ) φ φ ,   F 32 ( 2 ) 2 S ( φ ; ψ ; α ) ψ φ ;   F 41 ( 2 ) 2 S ( φ ; ψ ; α ) φ   ψ ;   F 42 ( 2 ) 2 S ( φ ; ψ ; α ) ψ ψ ;
F 51 ( 3 ) 3 S ( φ ; ψ ; α ) φ φ φ ψ 3 , i ( 2 ) ( x ) 3 S ( φ ; ψ ; α ) φ φ   ψ ψ 4 , i ( 2 ) ( x ) φ [ R i ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) / φ ] ;
F 52 ( 3 ) 3 S ( φ ; ψ ; α ) ψ φ φ ψ 3 , i ( 2 ) ( x ) 3 S ( φ ; ψ ; α ) φ ψ   ψ ψ 4 , i ( 2 ) ( x ) [ R i ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) / φ ] ψ ;
F 53 ( 3 ) ψ 1 ( 1 ) [ R i ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) / φ ] ;   F 54 ( 3 ) ψ 2 ( 1 ) [ R i ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) / φ ] ;
F 55 ( 3 ) L + ( α ) ;   F 56 ( 3 ) 0 ;   F 57 ( 3 ) 2 S ( φ ; ψ ; α ) φ φ ;   F 58 ( 3 ) 2 S ( φ ; ψ ; α ) φ   ψ ;
F 61 ( 3 ) 3 S ( φ ; ψ ; α ) ψ φ φ ψ 3 , i ( 2 ) ( x ) 3 S ( φ ; ψ ; α ) φ ψ ψ ψ 4 , i ( 2 ) ( x ) φ [ R i ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) / ψ ] ;
F 62 ( 3 ) 3 S ( φ ; ψ ; α ) ψ ψ φ ψ 3 , i ( 2 ) ( x ) 3 S ( φ ; ψ ; α ) ψ ψ ψ ψ 4 , i ( 2 ) ( x ) ψ [ R i ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) / ψ ] ;
F 63 ( 3 ) ψ 1 ( 1 ) [ R i ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) / ψ ] ;   F 64 ( 3 ) ψ 2 ( 1 ) [ R i ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) / ψ ] ;
F 65 ( 3 ) 0 ;   F 66 ( 3 ) L ( α ) ;   F 67 ( 3 ) 2 S ( φ ; ψ ; α ) ψ φ ;   F 68 ( 3 ) 2 S ( φ ; ψ ; α ) ψ ψ ;
F 71 ( 3 ) φ [ R i ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) / ψ 1 ( 1 ) ] ;   F 72 ( 3 ) ψ [ R i ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) / ψ 1 ( 1 ) ] ;
F 73 ( 3 ) ψ 1 ( 1 ) [ R i ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) / ψ 1 ( 1 ) ] ;   F 73 ( 3 ) ψ 2 ( 1 ) [ R i ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) / ψ 1 ( 1 ) ] ;
F 75 ( 3 ) 0 ;   F 76 ( 3 ) 0 ;   F 77 ( 3 ) L ( α ) ;   F 78 ( 3 ) 0 ;
F 81 ( 3 ) φ [ R i ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) / ψ 2 ( 1 ) ] ;   F 82 ( 3 ) ψ [ R i ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) / ψ 2 ( 1 ) ] ;
F 83 ( 3 ) ψ 1 ( 1 ) [ R i ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) / ψ 2 ( 1 ) ] ;   F 84 ( 3 ) ψ 2 ( 1 ) [ R i ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) / ψ 2 ( 1 ) ] ;
F 85 ( 3 ) 0 ;   F 66 ( 3 ) 0 ;   F 87 ( 3 ) 0 ;   F 88 ( 3 ) L + ( α ) ;
Q 1 , i ( 3 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; ψ 1 , i ( 2 ) ; ; ψ 4 , i ( 2 ) ; α ; δ α ) α { [ R i ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) / φ ] L + ( α ) ψ 1 , i ( 2 ) ( x ) 2 S ( φ ; ψ ; α ) φ φ ψ 3 , i ( 2 ) ( x ) 2 S ( φ ; ψ ; α ) φ   ψ ψ 4 , i ( 2 ) ( x ) } δ α ,
Q 2 , i ( 3 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; ψ 1 , i ( 2 ) ; ; ψ 4 , i ( 2 ) ; α ; δ α ) α { [ R i ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) / ψ ]   L ( α ) ψ 2 , i ( 2 ) ( x ) 2 S ( φ ; ψ ; α ) ψ φ ψ 3 , i ( 2 ) ( x ) 2 S ( φ ; ψ ; α ) ψ ψ ψ 4 , i ( 2 ) ( x ) } δ α ,
Q 3 , i ( 3 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; ψ 1 , i ( 2 ) ; ; ψ 4 , i ( 2 ) ; α ; δ α ) α { [ R i ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) / ψ 1 ( 1 ) ] L ( α ) ψ 3 , i ( 2 ) ( x ) } δ α ,
Q 4 , i ( 3 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; ψ 1 , i ( 2 ) ; ; ψ 4 , i ( 2 ) ; α ; δ α ) α { [ R i ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) / ψ 2 ( 1 ) ] L + ( α ) ψ 4 , i ( 2 ) ( x ) } δ α   .
The boundary conditions for the functions   δ φ , δ ψ , δ ψ 1 ( 1 ) , δ ψ 2 ( 1 ) δ ψ 1 , i ( 2 ) ,   δ ψ 2 , i ( 2 ) ,   δ ψ 3 , i ( 2 ) ,   δ ψ 4 , i ( 2 ) , are those provided in Equations (17) and (29), augmented by the boundary conditions obtained by G-differentiating Equation (36), i.e.:
δ B A ( 2 ) ( φ α ) B A ( 2 ) ( φ ) φ δ φ + B A ( 2 ) ( ψ ) ψ δ ψ + B A ( 2 ) ( ψ 1 ( 1 ) ) ψ 1 ( 1 ) δ ψ 1 ( 1 ) + B A ( 2 ) ( ψ 2 ( 1 ) ) ψ 2 ( 1 ) δ ψ 2 ( 1 ) + B A ( 2 ) ( ψ 1 , i ( 2 ) ) ψ 1 , i ( 2 ) δ ψ 1 , i ( 2 ) + B A ( 2 ) ( ψ 2 , i ( 2 ) ) ψ 2 , i ( 2 ) δ ψ 2 , i ( 2 ) + B A ( 2 ) ( ψ 3 , i ( 2 ) ) ψ 3 , i ( 2 ) δ ψ 3 , i ( 2 ) + B A ( 2 ) ( ψ 4 , i ( 2 ) ) ψ 4 , i ( 2 ) δ ψ 4 , i ( 2 ) + B A ( 1 ) ( α ) α δ α = 0 ,   x Ω x .
The system comprising the 2nd-LFSS together with Equation (42) and the boundary conditions provided in Equation (62) is called the 3rd-Level Forward Sensitivity System (3rd-LFSS). Since the source-terms of the 3rd-LFSS depend on the parameters variations δ α i , it follows that solving this system of equations is prohibitively expensive computationally. To avoid the need for solving the 3rd-LFSS, the indirect-effect term defined in Equation (41) will be expressed in terms of a 3rd-Level Adjoint Sensitivity System (3rd-LASS), which will be constructed by implementing the following sequence of steps:
(i) Define a Hilbert space, denoted as H ( 3 ) , having vector-valued elements of the form the u ( 3 ) ( x ) [ u 1 ( 3 ) ( x ) , , u 8 ( 3 ) ( x ) ] H ( 3 ) , with components that are N α -dimensional vectors of the form u i ( 3 ) ( x ) [ u i , 1 ( 3 ) ( x ) , , u i , j ( 3 ) ( x ) , , u i , N u ( 3 ) ( x ) ] ,   i = 1 , , 8 , with square-integrable components u i , j ( 3 ) ( x ) . In H ( 3 ) , define the inner-product, denoted as u ( 3 ) ( x ) , v ( 3 ) ( x ) ( 3 ) , of two functions u ( 3 ) ( x ) H ( 3 ) and v ( 3 ) ( x ) H ( 3 ) as follows:
u ( 3 ) ( x ) , v ( 3 ) ( x ) ( 3 ) i = 1 8 u i ( 3 ) ( x ) , v i ( 3 ) ( x ) φ = i = 1 8 j = 1 N u Ω x u i , j ( 3 ) ( x ) v i , j ( 3 ) ( x )   d x .
Using the definition provided in Equation (63), construct the inner product of a vector ψ i j ( 3 ) ( x ) [ ψ 1 , i j ( 3 ) ( x ) , , ψ 8 , i j ( 3 ) ( x ) ] H ( 3 ) with Equations (16) and (28) to obtain the following relation:
( ψ 1 , i j ( 3 ) ψ 2 , i j ( 3 ) ψ 3 , i j ( 3 ) ψ 4 , i j ( 3 ) ψ 5 , i j ( 3 ) ψ 6 , i j ( 3 ) ψ 7 , i j ( 3 ) ψ 8 , i j ( 3 ) ) ,   ( L ( α ) 0 0 0 0 0 0 0 0 L + ( α ) 0 0 0 0 0 0 F 31 ( 2 ) F 32 ( 2 ) L + ( α ) 0 0 0 0 0 F 41 ( 2 ) F 42 ( 2 ) 0 L ( α ) 0 0 0 0 F 51 ( 3 ) F 52 ( 3 ) F 53 ( 3 ) F 54 ( 3 ) L + ( α ) 0 F 57 ( 3 ) F 58 ( 3 ) F 61 ( 3 ) F 62 ( 3 ) F 63 ( 3 ) F 64 ( 3 ) 0 L ( α ) F 67 ( 3 ) F 68 ( 3 ) F 71 ( 3 ) F 72 ( 3 ) F 73 ( 3 ) F 74 ( 3 ) 0 0 L ( α ) 0 F 81 ( 3 ) F 82 ( 3 )   F 83 ( 3 ) F 84 ( 3 ) 0 0 0 L + ( α ) ) ( δ φ δ ψ δ ψ 1 ( 1 ) δ ψ 2 ( 1 ) δ ψ 1 , i ( 2 ) δ ψ 2 , i ( 2 ) δ ψ 3 , i ( 2 ) δ ψ 4 , i ( 2 ) ) ( 3 ) = = ( ψ 1 , i j ( 3 ) ( x ) ψ 2 , i j ( 3 ) ( x ) ψ 3 , i j ( 3 ) ( x ) ψ 4 , i j ( 3 ) ( x ) ψ 5 , i j ( 3 ) ( x ) ψ 6 , i j ( 3 ) ( x ) ψ 7 , i j ( 3 ) ( x ) ψ 8 , i j ( 3 ) ( x ) ) ,   ( Q 1 ( 1 ) ( φ ; α ; δ α ) Q 2 ( 1 ) ( ψ ; α ; δ α ) Q 1 ( 2 ) ( φ ; ; δ α ) Q 2 ( 2 ) ( φ ; ; δ α ) Q 1 , i ( 3 ) ( φ ; ; δ α ) Q 2 , i ( 3 ) ( φ ; ; δ α ) Q 3 , i ( 3 ) ( φ ; ; δ α ) Q 4 , i ( 3 ) ( φ ; ; δ α ) ) ( 3 ) .
(ii)
Use the definition of the adjoint operator in the Hilbert space H ( 3 ) to recast the left-side of Equation (64) as follows:
( ψ 1 , i j ( 3 ) ψ 2 , i j ( 3 ) ψ 3 , i j ( 3 ) ψ 4 , i j ( 3 ) ψ 5 , i j ( 3 ) ψ 6 , i j ( 3 ) ψ 7 , i j ( 3 ) ψ 8 , i j ( 3 ) ) ,   ( L ( α ) 0 0 0 0 0 0 0 0 L + ( α ) 0 0 0 0 0 0 F 31 ( 2 ) F 32 ( 2 ) L + ( α ) 0 0 0 0 0 F 41 ( 2 ) F 42 ( 2 ) 0 L ( α ) 0 0 0 0 F 51 ( 3 ) F 52 ( 3 ) F 53 ( 3 ) F 54 ( 3 ) L + ( α ) 0 F 57 ( 3 ) F 58 ( 3 ) F 61 ( 3 ) F 62 ( 3 ) F 63 ( 3 ) F 64 ( 3 ) 0 L ( α ) F 67 ( 3 ) F 68 ( 3 ) F 71 ( 3 ) F 72 ( 3 ) F 73 ( 3 ) F 74 ( 3 ) 0 0 L ( α ) 0 F 81 ( 3 ) F 82 ( 3 )   F 83 ( 3 ) F 84 ( 3 ) 0 0 0 L + ( α ) ) ( δ φ δ ψ δ ψ 1 ( 1 ) δ ψ 2 ( 1 ) δ ψ 1 , i ( 2 ) δ ψ 2 , i ( 2 ) δ ψ 3 , i ( 2 ) δ ψ 4 , i ( 2 ) ) ( 3 ) = ( δ φ δ ψ δ ψ 1 ( 1 ) δ ψ 2 ( 1 ) δ ψ 1 , i ( 2 ) δ ψ 2 , i ( 2 ) δ ψ 3 , i ( 2 ) δ ψ 4 , i ( 2 ) ) ,   ( L + ( α ) 0 F 31 ( 2 ) F 41 ( 2 ) F 51 ( 3 ) F 61 ( 3 ) F 71 ( 3 ) F 81 ( 3 ) 0 L ( α ) F 32 ( 2 ) F 42 ( 2 ) F 52 ( 3 ) F 62 ( 3 ) F 72 ( 3 ) F 82 ( 3 ) 0 0 L ( α ) 0 F 53 ( 3 ) F 63 ( 3 ) F 73 ( 3 ) F 83 ( 3 ) 0 0 0 L + ( α ) F 54 ( 3 ) F 64 ( 3 ) F 74 ( 3 ) F 84 ( 3 ) 0 0 0 0 L ( α ) 0 0 0 0 0 0 0 0 L + ( α ) 0 0 0 0 0 0 F 57 ( 3 ) F 67 ( 3 ) L + ( α ) 0 0 0 0 0 F 58 ( 3 ) F 68 ( 3 ) 0 L ( α ) ) ( ψ 1 , i j ( 3 ) ψ 2 , i j ( 3 ) ψ 3 , i j ( 3 ) ψ 4 , i j ( 3 ) ψ 5 , i j ( 3 ) ψ 6 , i j ( 3 ) ψ 7 , i j ( 3 ) ψ 8 , i j ( 3 ) ) ( 3 ) + P ( 3 ) [ δ φ ; δ ψ ; δ ψ 1 ( 1 ) , δ ψ 2 ( 1 ) ; δ ψ 1 , i ( 2 ) , , δ ψ 4 , i ( 2 ) ; ψ 1 , i j ( 3 ) , , ψ 8 , i j ( 3 ) ]   ,
where the bilinear concomitant P ( 3 ) [ δ φ ; δ ψ ; δ ψ 1 ( 1 ) , δ ψ 2 ( 1 ) ; δ ψ 1 , i ( 2 ) , , δ ψ 4 , i ( 2 ) ; ψ 1 , i j ( 3 ) , , ψ 8 , i j ( 3 ) ] is defined on the phase-space boundary x Ω x .
(iii)
Identify the first term on the right-side of Equation (65) with the indirect-effect term defined in Equation (41) by requiring that:
( L + ( α ) 0 F 31 ( 2 ) F 41 ( 2 ) F 51 ( 3 ) F 61 ( 3 ) F 71 ( 3 ) F 81 ( 3 ) 0 L ( α ) F 32 ( 2 ) F 42 ( 2 ) F 52 ( 3 ) F 62 ( 3 ) F 72 ( 3 ) F 82 ( 3 ) 0 0 L ( α ) 0 F 53 ( 3 ) F 63 ( 3 ) F 73 ( 3 ) F 83 ( 3 ) 0 0 0 L + ( α ) F 54 ( 3 ) F 64 ( 3 ) F 74 ( 3 ) F 84 ( 3 ) 0 0 0 0 L ( α ) 0 0 0 0 0 0 0 0 L + ( α ) 0 0 0 0 0 0 F 57 ( 3 ) F 67 ( 3 ) L + ( α ) 0 0 0 0 0 F 58 ( 3 ) F 68 ( 3 ) 0 L ( α ) ) ( ψ 1 , i j ( 3 ) ψ 2 , i j ( 3 ) ψ 3 , i j ( 3 ) ψ 4 , i j ( 3 ) ψ 5 , i j ( 3 ) ψ 6 , i j ( 3 ) ψ 7 , i j ( 3 ) ψ 8 , i j ( 3 ) ) = ( [ R i j ( 2 ) / φ ] [ R i j ( 2 ) / ψ ] [ R i j ( 2 ) / ψ 1 ( 1 ) ] [ R i j ( 2 ) / ψ 2 ( 1 ) ] [ R i j ( 2 ) / ψ 1 , i ( 2 ) ] [ R i j ( 2 ) / ψ 2 , i ( 2 ) ] [ R i j ( 2 ) / ψ 3 , i ( 2 ) ] [ R i j ( 2 ) / ψ 4 , i ( 2 ) ] ) ,   f o r   i = 1 , N α ;   j = 1 , , i .
(iv)
The boundary conditions given for the 2nd-LFSS and those given in Equation (62) are now implemented in Equation (65), thereby reducing by half the number of unknown boundary-values in the bilinear concomitant P ( 3 ) [ δ φ ; δ ψ ; δ ψ 1 ( 1 ) , δ ψ 2 ( 1 ) ; δ ψ 1 , i ( 2 ) , , δ ψ 4 , i ( 2 ) ; ψ 1 , i j ( 3 ) , , ψ 8 , i j ( 3 ) ] . The boundary conditions for the 3rd-level adjoint functions ψ 1 , i j ( 3 ) , , ψ 8 , i j ( 3 ) are chosen next so as to eliminate the remaining unknown boundary-values of the functions δ φ ; δ ψ ; δ ψ 1 ( 1 ) , δ ψ 2 ( 1 ) ; δ ψ 1 , i ( 2 ) , , δ ψ 4 , i ( 2 ) while ensuring that Equation (66) is well posed. The boundary conditions thus chosen for the adjoint functions ψ 1 , i j ( 3 ) , , ψ 8 , i j ( 3 ) can be represented in operator form as follows:
B A ( 3 ) [ φ ( x ) ; ψ ( x ) ; ψ 1 ( 1 ) ( x ) ; ψ 2 ( 1 ) ( x ) ; ψ 1 , i ( 2 ) ( x ) , , ψ 4 , i ( 2 ) ( x ) ; ψ 1 , i j ( 3 ) , , ψ 8 , i j ( 3 ) ; α ] = 0 ,   x Ω x .
In most cases, the above choice of boundary conditions for the 3rd-level adjoint functions ψ 1 , i j ( 3 ) , , ψ 8 , i j ( 3 ) will cause the bilinear concomitant P ( 3 ) [ δ φ ; δ ψ ; δ ψ 1 ( 1 ) , δ ψ 2 ( 1 ) ; δ ψ 1 , i ( 2 ) , , δ ψ 4 , i ( 2 ) ; ψ 1 , i j ( 3 ) , , ψ 8 , i j ( 3 ) ] in Equation (65) to vanish. Even when it does not vanish, however, this bilinear concomitant will be reduced to a quantity, denoted here as P ^ ( 3 ) [ φ ; ψ ; ψ 1 ( 1 ) , ψ 2 ( 1 ) ; ψ 1 , i ( 2 ) , , ψ 4 , i ( 2 ) ; ψ 1 , i j ( 3 ) , , ψ 8 , i j ( 3 ) ; α ; δ α ] , which will contain only known values of its arguments. The system of equations comprising Equations (66) and (67) will be called the 3rd-Level Adjoint Sensitivity System (3rd-LASS) for the 3rd-level adjoint functions ψ 1 , i j ( 3 ) , , ψ 8 , i j ( 3 ) , i = 1 , N α ;   j = 1 , , i .
(v)  Use the 3rd-LASS together with Equations (65) and (64) to obtain the following expression for the indirect-effect term defined in Equation (41), in terms of the 3rd-level adjoint functions ψ 1 , i j ( 3 ) , , ψ 8 , i j ( 3 ) :
{ δ R i j ( 2 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; ψ 1 , i ( 2 ) , , ψ 4 , i ( 2 ) ; ψ 1 , i j ( 3 ) , , ψ 8 , i j ( 3 ) ; α ; δ α ) } i n d = ψ 1 , i j ( 3 ) ( x ) ,   Q 1 ( 1 ) ( φ ; α ; δ α ) φ + ψ 2 , i j ( 3 ) ( x ) ,   Q 2 ( 1 ) ( ψ ; α ; δ α ) φ + ψ 3 , i j ( 3 ) ( x ) ,   Q 1 ( 2 ) ( φ ; ; δ α ) φ   + ψ 4 , i j ( 3 ) ( x ) ,   Q 2 ( 2 ) ( φ ; ; δ α ) φ + ψ 5 , i j ( 3 ) ( x ) ,   Q 1 , i ( 3 ) ( φ ; ; δ α ) φ + ψ 6 , i j ( 3 ) ( x ) ,   Q 2 , i ( 3 ) ( φ ; ; δ α ) φ + ψ 7 , i j ( 3 ) ( x ) ,   Q 3 , i ( 3 ) ( φ ; ; δ α ) φ + ψ 8 , i j ( 3 ) ( x ) ,   Q 4 , i ( 3 ) ( φ ; ; δ α ) φ P ^ ( 3 ) [ φ ; ψ ; ψ 1 ( 1 ) , ψ 2 ( 1 ) ; ψ 1 , i ( 2 ) , , ψ 4 , i ( 2 ) ; ψ 1 , i j ( 3 ) , , ψ 8 , i j ( 3 ) ; α ; δ α ] ,   i = 1 , , N α ;   j = 1 , , i .
As Equation (68) indicates, the indirect-effect term can be computed speedily by quadratures once the 3rd-level adjoint functions ψ 1 , i j ( 3 ) , , ψ 8 , i j ( 3 ) become available.
(vi)
Replace Equation (68) in Equation (40) to obtain the following expression for the total 2nd-order response sensitivity to model parameters:
δ R i j ( 2 ) [ φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; ψ 1 , i ( 2 ) , , ψ 4 , i ( 2 ) ; ψ 1 , i j ( 3 ) , , ψ 8 , i j ( 3 ) ; α ; δ α ] k = 1 N α R i j k ( 3 ) [ φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; ψ 1 , i ( 2 ) , , ψ 4 , i ( 2 ) ; ψ 1 , i j ( 3 ) , , ψ 8 , i j ( 3 ) ; α ] δ α k ,   i = 1 , , N α ;   j = 1 , , i ,
where the quantity R i j k ( 3 ) [ φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; ψ 1 , i ( 2 ) , , ψ 4 , i ( 2 ) ; ψ 1 , i j ( 3 ) , , ψ 8 , i j ( 3 ) ; α ] denotes the 3rd-order partial sensitivity of the response to the model parameters and is defined as follows:
R i j k ( 3 ) [ φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; ψ 1 , i ( 2 ) , , ψ 4 , i ( 2 ) ; ψ 1 , i j ( 3 ) , , ψ 8 , i j ( 3 ) ; α ] R i j ( 2 ) ( φ ; ψ ; ; α ) α k ψ 1 , i j ( 3 ) ( x ) ,   [ Q [ α ( x ) ] L ( α ) φ ( x ) ] / α k φ + ψ 2 , i j ( 3 ) ( x ) ,   [ Q + [ α ( x ) ] L + ( α ) ψ ( x ) ] / α k φ + ψ 3 , i j ( 3 ) ( x ) ,   2 S ( φ ; ψ ; α ) φ α k [ L + ( α ) ψ 1 ( 1 ) ( x ) ] α k φ + ψ 4 , i j ( 3 ) ( x ) ,   2 S ( φ ; ψ ; α ) ψ α k [ L ( α ) ψ 2 ( 1 ) ( x ) ] α k φ + ψ 5 , i j ( 3 ) ( x ) ,   α k { [ R i ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) / φ ] L + ( α ) ψ 1 , i ( 2 ) ( x )      2 S ( φ ; ψ ; α ) φ φ ψ 3 , i ( 2 ) ( x ) 2 S ( φ ; ψ ; α ) φ   ψ ψ 4 , i ( 2 ) ( x ) } φ + ψ 6 , i j ( 3 ) ( x ) ,   α k { [ R i ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) / ψ ] L ( α ) ψ 2 , i ( 2 ) ( x )      2 S ( φ ; ψ ; α ) ψ φ ψ 3 , i ( 2 ) ( x ) 2 S ( φ ; ψ ; α ) ψ ψ ψ 4 , i ( 2 ) ( x ) } φ + ψ 7 , i j ( 3 ) ( x ) ,   α k { [ R i ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) / ψ 1 ( 1 ) ] L ( α ) ψ 3 , i ( 2 ) ( x ) } φ + ψ 8 , i j ( 3 ) ( x ) ,   α k { [ R i ( 1 ) ( φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; α ) / ψ 2 ( 1 ) ] L + ( α ) ψ 4 , i ( 2 ) ( x ) } φ P ^ ( 3 ) [ φ ; ψ ; ψ 1 ( 1 ) , ψ 2 ( 1 ) ; ψ 1 , i ( 2 ) , , ψ 4 , i ( 2 ) ; ψ 1 , i j ( 3 ) , , ψ 8 , i j ( 3 ) ; α ] α k ,   i = 1 , , N α ;   j = 1 , , i ;   k = 1 , , j .
(vii)
Note that the 2nd-LASS is independent of parameter variations δ α . Thus, the exact computation of all of the partial third-order sensitivities, R i j k ( 3 ) [ φ ; ψ ; ψ 1 ( 1 ) ; ψ 2 ( 1 ) ; ψ 1 , i ( 2 ) , , ψ 4 , i ( 2 ) ; ψ 1 , i j ( 3 ) , , ψ 8 , i j ( 3 ) ; α ] , i = 1 , , N α ;   j = 1 , , i ;   k = 1 , , j , requires at most N α ( N α + 1 ) / 2 large-scale (adjoint) computations using the 3rd-LASS, rather than N α ( N α + 1 ) ( N α + 2 ) / 6 large-scale computations as would be required by forward methods. In order to implement the practical computation of the 3rd-level adjoint functions, ψ 1 , i j ( 3 ) , , ψ 8 , i j ( 3 ) it is important to note that, in component form, Equation (66) has the following structure, for each i = 1 , N α ;   j = 1 , , i :
L ( α ) ψ 5 , i j ( 3 ) = [ R i j ( 2 ) / ψ 1 , i ( 2 ) ] ,
L + ( α ) ψ 6 , i j ( 3 ) = [ R i j ( 2 ) / ψ 2 , i ( 2 ) ] ,
L + ( α ) ψ 7 , i j ( 3 ) = [ R i j ( 2 ) / ψ 3 , i ( 2 ) ] F 57 ( 3 ) ψ 5 , i j ( 3 ) F 67 ( 3 ) ψ 6 , i j ( 3 )   ,
L ( α ) ψ 8 , i j ( 3 ) = [ R i j ( 2 ) / ψ 4 , i ( 2 ) ] F 58 ( 3 ) ψ 5 , i j ( 3 ) F 68 ( 3 ) ψ 6 , i j ( 3 )   ,
L + ( α ) ψ 4 , i j ( 3 ) = [ R i j ( 2 ) / ψ 2 ( 1 ) ] F 54 ( 3 ) ψ 5 , i j ( 3 ) F 64 ( 3 ) ψ 6 , i j ( 3 ) F 74 ( 3 ) ψ 7 , i j ( 3 ) F 84 ( 3 ) ψ 8 , i j ( 3 )   ,
L ( α ) ψ 3 , i j ( 3 ) = [ R i j ( 2 ) / ψ 1 ( 1 ) ] F 53 ( 3 ) ψ 5 , i j ( 3 ) F 63 ( 3 ) ψ 6 , i j ( 3 ) F 73 ( 3 ) ψ 7 , i j ( 3 ) F 83 ( 3 ) ψ 8 , i j ( 3 )   ,
L ( α ) ψ 2 , i j ( 3 ) = [ R i j ( 2 ) / ψ ] F 32 ( 2 ) ψ 3 , i j ( 3 ) F 42 ( 2 ) ψ 4 , i j ( 3 ) F 52 ( 3 ) ψ 5 , i j ( 3 ) F 62 ( 3 ) ψ 6 , i j ( 3 ) F 72 ( 3 ) ψ 7 , i j ( 3 ) F 82 ( 3 ) ψ 8 , i j ( 3 )   ,
L + ( α ) ψ 1 , i j ( 3 ) = [ R i j ( 2 ) / φ ] F 31 ( 2 ) ψ 3 , i j ( 3 ) F 41 ( 2 ) ψ 4 , i j ( 3 ) F 51 ( 3 ) ψ 5 , i j ( 3 ) F 61 ( 3 ) ψ 6 , i j ( 3 ) F 71 ( 3 ) ψ 7 , i j ( 3 ) F 81 ( 3 ) ψ 8 , i j ( 3 )   .
Thus, the 3rd-LASS it is to be solved successively, by first using Equations (71) and (72) to compute the 3rd-level adjoint functions ψ 5 , i j ( 3 ) and ψ 6 , i j ( 3 ) . Note that solving using Equations (71) and (72) would be performed by using the same forward and adjoint solvers (i.e., computer codes) as used for solving the original forward and adjoint systems, namely Equations (1) and (6) subject to the corresponding boundary conditions, except that the right-sides of the respective solvers would have as “sources” the terms R i j ( 2 ) / ψ 1 , i ( 2 ) and R i j ( 2 ) / ψ 2 , i ( 2 ) , respectively. After having obtained the 3rd-level adjoint functions ψ 5 , i j ( 3 ) and ψ 6 , i j ( 3 ) , the next round of computations would be to solve Equations (73) and (74) in order to determine the 3rd-level adjoint functions ψ 7 , i j ( 3 ) and ψ 8 , i j ( 3 ) , respectively. As before, solving using Equations (73) and (74) would be performed by using the same forward and adjoint solvers (i.e., computer codes) as used for solving the original forward and adjoint systems, namely Equations (1) and (6) subject to the corresponding boundary conditions, except that the right-sides of the respective solvers would have different “sources.” The next round of computations would require the use of the forward and adjoint solvers to solve Equations (75) and (76) in order to determine the 3rd-level adjoint functions ψ 4 , i j ( 3 ) and ψ 3 , i j ( 3 ) , respectively. The final set of computations would require the use of the forward and adjoint solvers to solve Equations (77) and (78) in order to determine the 3rd-level adjoint functions ψ 2 , i j ( 3 ) and ψ 1 , i j ( 3 ) , respectively. Thus, solving the 3rd-ASAM in order to determine the 3rd-level adjoint functions does not require any significant “code development,” since the original forward and adjoint solvers (codes) do not need to be modified; only the right-sides (i.e., “sources”) for these solvers/codes would need to be programmed accordingly.
Using the 3rd-LASS enables the specific computation of the 3rd-order sensitivities in the priority order set by the user, so that only the important 3rd-order partial sensitivities R [ φ ( x ) , ψ ( x ) ; α ] would be computed. Information provided by the first- and second-order sensitivities might indicate which 3rd-order sensitivities could be neglected.

4. Third-Order Expressions for the Cumulants of the Response Distribution in Parameter Space

The 3rd-ASAM presented in Section 3, above, provides the most efficient way for computing exactly the first-, second- and third-order sensitivities of a response that couples the forward and adjoint systems that describe physical problems which are linear in the state-functions. The availability of these sensitivities enables the use of a third-order multivariate Taylor series-expansion (of the response around the known nominal parameter values) for quantifying the cumulants of the response distribution in the phase-space of model parameters. For a model’s computed-response, denoted as r i 1 c ( α ) , where the superscript “c” denotes “computed” and the subscript i 1 = 1 , , N r denotes one of a total of N r responses that would be of interest, the third-order Taylor-series of r a c ( α ) around the model’s parameters’ mean (expected) values α 0 ( α 1 0 , , α N α 0 ) is:
r i 1 c ( α ) = r i 1 c ( α 0 ) + i = 1 N α { r i 1 c ( α ) α i } α 0 ( α i α i 0 ) + 1 2 i , j = 1 N α { 2 r i 1 c ( α ) α i α j } α 0 ( α i α i 0 ) ( α j α j 0 ) + 1 6 i , j , k = 1 N α { 3 r i 1 c ( α ) α i α j α k } α 0 ( α i α i 0 ) ( α j α j 0 ) ( α k α k 0 ) + ;   a = 1 , , N r .
where r i 1 c ( α 0 ) denotes the nominal value of the response computed at the nominal (mean) parameter values α 0 ( α 1 0 , , α N α 0 ) .
In practice, the values of the parameters α n are determined experimentally. Therefore, these parameters can be considered to be variates that behave stochastically, obeying a multivariate probability distribution function, denoted as p α ( α ) , which is seldom known in practice, particularly for large-scale systems involving many parameters. Considering that the multivariate distribution p α ( α ) is formally defined on a domain D α , the various moments (e.g., mean values, covariance and variances, etc.) of p α ( α ) can be defined in a standard manner by using the following notation:
u ( α ) α D α u ( α ) p α ( α ) d α
where u ( α ) is a continuous function of the parameters α . Using the notation defined in Equation (80), the expected (or mean) value of a model parameter α i , denoted as α i 0 , is defined as follows:
α i 0 α i α D α α i p α ( α ) d α ,   i = 1 , , N α
The covariance, c o v ( α i , α j ) , of two parameters, α i and α j , is defined as follows:
μ 2 i j ( α ) c o v ( α i , α j ) ( α i α i 0 ) ( α j α j 0 ) α ,   i , j = 1 , , N α
The variance, v a r ( α i ) , of a parameter α i , is defined as follows:
v a r ( α i ) ( α i α i 0 ) 2 α ,   i = 1 , , N α
The standard deviation, σ i , of α i , is defined as follows: σ i v a r ( α i ) ;
The correlation, ρ i j , between two parameters α i and α j , is defined as follows:
ρ i j c o v ( α i , α j ) / ( σ i σ j ) ;   i , j = 1 , , N α
The 3rd-order moment, μ 3 i j k , of the multivariate parameter distribution function p ( α ) , and the 3rd-order parameter correlation, t i j k , respectively, are defined as follows:
μ 3 i j k ( α ) D α ( α i α i 0 ) ( α j α j 0 ) ( α k α k 0 ) p ( α ) d α t i j k σ i σ j σ k ,   i , j , k = 1 , , N α ;
The 4th-order moment, μ 4 i j k l , of the multivariate parameter distribution function p ( α ) , and the 4th-order parameter correlation, q i j k l , respectively, are defined as follows:
μ 4 i j k l ( α ) D α ( α i α i 0 ) ( α j α j 0 ) ( α k α k 0 ) ( α l α l 0 ) p ( α ) d α q i j k l σ i σ j σ k σ l ;   i , j , k , l = 1 , , N α .
Using Equations (79) together with Equations (81) through (85) yields the following expression for the expected (mean) value, denoted as E [ r i 1 c ( α ) ] , of a response r i 1 c ( α ) :
E [ r i 1 c ( α ) ] = r i 1 c ( α 0 ) + 1 2 i , j = 1 N α { 2 r i 1 c ( α ) α i α j } α 0 ρ i j σ i σ j + 1 6 i , j , k = 1 N α { 3 r i 1 c ( α ) α i α j α k } α 0 t i j k σ i σ j σ k ;   i 1 = 1 , , N r .
Using Equations (79) together with Equations (81) through (86) yields the following expression for the covariance, denoted as   c o v ( r i 1 c ,   r i 2 c ) , of two responses, r i 1 c ( α ) and r i 2 c ( α ) f o r   i 1 ,   i 2 = 1 , , N r :
c o v ( r i 1 c ,   r i 2 c ) = i = 1 N α j = 1 N α ( r i 1 c α i r i 2 c α j ) ρ i j σ i σ j + 1 2 i = 1 N α j = 1 N α μ = 1 N α ( 2 r i 1 c α i α j r i 2 c α μ + r i 1 c α i 2 r i 2 c α j α μ ) t i j μ σ i σ j σ μ + 1 4 i = 1 N α j = 1 N α μ = 1 N α ν = 1 N α ( 2 r i 1 c α i α j ) ( 2 r i 2 c α μ α ν ) ( q i j μ ν ρ i j ρ μ ν ) σ i σ j σ μ σ ν + 1 6 i = 1 N α j = 1 N α μ = 1 N α ν = 1 N α ( r i 1 c α i 3 r i 2 c α j α μ α ν + r i 2 c α i 3 r i 1 c α j α μ α ν )   q i j μ ν σ i σ j σ μ σ ν   .
In particular, the variance of a response r i 1 c ( α ) is obtained as by setting i 1 = i 2 in Equation (88).
The covariance of a response, r i 1 c ( α ) and a parameter α , i 1 = 1 , , N r and = 1 , , N α , which is denoted as   c o v ( r i 1 c ,   α ) and is given by the following expression:
c o v ( r i 1 c ,   α ) = i = 1 N α { r i 1 c ( α ) α i } α 0 c o v ( α i , α ) + 1 2 i , j = 1 N α { 2 r i 1 c ( α ) α i α j } α 0 t i j σ i σ j σ   + 1 6 i , j , k = 1 N α { 3 r i 1 c ( α ) α i α j α k } α 0 q i j μ ν σ i σ j σ k σ ;   i 1 = 1 , , N r ;   = 1 , , N α .
The third-order cumulant for three responses, r i 1 c ( α ) , r i 2 c ( α ) and r i 3 c ( α ) , which is denoted below as μ 3 ( r i 1 c , r i 2 c , r i 3 c ) , f o r   i 1 , i 2 , i 3 = 1 , , N r , is obtained similarly by using Equations (79) together with Equations (81) through (86), and has the following expression:
μ 3 ( r i 1 c , r i 2 c , r i 3 c ) [ r i 1 c E ( r i 1 c ) ] [ r i 2 c E ( r i 2 c ) ] [ r i 3 c E ( r i 3 c ) ] α = i = 1 N α j = 1 N α μ = 1 N α r i 1 c α i r i 2 c α j r i 3 c α μ t i j μ σ i σ j σ μ + 1 2 i = 1 N α j = 1 N α μ = 1 N α ν = 1 N α r i 1 c k α i r i 2 c α j 2 r i 3 c α μ α ν ( q i j μ ν ρ i j ρ μ ν ) σ i σ j σ μ σ ν + 1 2 i = 1 N α j = 1 N α μ = 1 N α ν = 1 N α r i 1 c α i 2 r i 2 c α j α μ r i 3 c α ν ( q i j μ ν ρ i ν ρ j μ ) σ i σ j σ μ σ ν + 1 2 i = 1 N α j = 1 N α μ = 1 N α ν = 1 N α 2 r i 1 c α i α j r i 2 c α μ r i 3 c α ν ( q i j μ ν ρ i j ρ μ ν ) σ i σ j σ μ σ ν   .
In particular, the skewness of a single response is customarily denoted as γ 1 ( R ) , and is defined as follows: γ 1 ( R ) = μ 3 ( R ) [ v a r ( R ) ] 3 2 , where μ 3 ( R ) denotes the third central moment of the response distribution. For normally-distributed uncorrelated parameters and a single response, the expression in Equation (90) simplifies to yield μ 3 ( R ) = 3 i = 1 N α ( R α i ) 2 2 R α i 2 σ i 4 . As is well-known, the skewness provides a quantitative measure of the asymmetries in the respective distribution.
The first-order sensitivities contribute the leading terms to the second-, third-, and fourth-order moments of the response distribution, thus providing the leading contributions to the responses variance/covariances, skewness, and kurtosis. Obtaining the exact and complete set of first-order sensitivities of responses to model parameters is of paramount importance for any analysis of a computational model.
The second-order sensitivities contribute the leading correction terms to the response’s expected value (causing it to differ from the response’s computed value). The second-order sensitivities also contribute to the response variances and covariances. If the parameters follow a normal (Gaussian) multivariate distribution, the second-order sensitivities contribute the leading terms to the response’s third-order moment. Thus, neglecting the second-order response sensitivities to normally distributed parameters would nullify the third-order response correlations and hence would nullify the skewness of a response.
The fourth-order cumulant, denoted as μ 4 ( r i 1 c , r i 2 c , r i 3 c , r i 4 c ) f o r   i 1 , i 2 , i 3 , i 4 = 1 , , N r , among four responses r i 1 c ( α ) , r i 2 c ( α ) , r i 3 c ( α ) and r i 4 c ( α ) , is also obtained by using Equations (79). Up to fourth-order correlations in the model parameters, μ 4 ( r i 1 c , r i 2 c , r i 3 c , r i 4 c ) has the following expression:
μ 4 ( r i 1 c , r i 2 c , r i 3 c , r i 4 c ) [ r i 1 c E ( r i 1 c ) ] [ r i 2 c E ( r i 2 c ) ] [ r i 3 c E ( r i 3 c ) ] [ r i 4 c E ( r i 4 c ) ] α = i = 1 N α j = 1 N α μ = 1 N α ν = 1 N α r i 1 c k α i r i 2 c α j r i 3 c α μ r i 4 c α ν q i j μ ν σ i σ j σ μ σ ν   .

5. 2nd/3rd-Order Best-Estimated Results with Reduced Uncertainties Predictive Modeling (2nd/3rd-BERRU-PM) in the Joint Phase-Space of Responses and Parameters

Cacuci [7] has summarized the scope of “BERRU-PM” as follows: “BERRU-PM commences by identifying and characterizing the uncertainties involved in every step in the sequence of the numerical simulation processes that ultimately lead to a prediction. This includes characterizing: (a) errors and uncertainties in the data used in the simulation (e.g., input data, model parameters, initial conditions, boundary conditions, sources and forcing functions), (b) numerical discretization errors, and (c) uncertainties in (e.g., lack of knowledge of) the processes being modeled. Under ideal circumstances, the result of is a probabilistic description of possible future outcomes based on all recognized errors and uncertainties.”
Consider a vector-valued variate x ( x 1 , , x N ) , the components of which obey an unknown multivariate distribution p ( x ) . Of course, the probability distribution function p ( x ) would need to be properly normalized, i.e., it must also satisfy the constraint:
p ( x ) d x   = 1 .
Consider further that the moments of several known functions F k ( x ) over the unknown distribution p ( x ) , denoted as F k , and defined as:
F k = p ( x ) F k ( x )   d x   ,   k = 1 ,   2 ,   ,   K
are also known. The problem of reconstructing a function from a finite number of its moments has been investigated for many decades in the mathematical and physical sciences. For the purposes of predictive modeling, the main goal is to determine a probability density function p ( x ) which is consistent with the knowledge expressed by Equation (93) and introduces no unwarranted information. Such a probability density function, p ( x ) , can be constructed using the method of maximum entropy (“MaxEnt”), which generates the most conservative estimate of a probability distribution with the given information and the most non-committal one with regard to missing information [11]. According to the MaxEnt principle, the unknown probability density function p ( x ) must satisfy the constraints expressed by Equation (93) while having its Boltzmann–Shannon–Gibbs (BSG) entropy (also referred to as the “information entropy”) as large as possible. For a continuous distribution having a probability density function p ( x ) , the expression for its information/BSG-entropy is:
H ( p ) = d x   p ( x ) l n p ( x ) m ( x )
where m ( x ) is a prior density function that ensures form invariance under change of variable. Intuitively, in a bounded domain, the most conservative distribution, i.e., the distribution of maximum entropy, is the one that assigns equal probability to all the accessible states. Hence, the method of maximum entropy can be thought of as choosing the most “uniform” distribution p ( x ) that satisfies the given moment constraints expressed by Equation (93) and introduces no unwarranted information. Any probability density function p ( x ) satisfying the constraints which has smaller entropy will contain more information (less uncertainty), and thus would predict something stronger than warranted by our knowledge and/or assumptions. The probability density function with maximum entropy, satisfying the constraints imposed, is the one which should be least surprising in terms of the predictions it makes.
Selecting the unique p ( x ) which would maximize the information entropy defined by Equation (94) while simultaneously satisfying the known constraints given in Equations (93) and (92) is a variational problem that can be solved by the well-known method of using Lagrange multipliers, λ k ,   k = 1 ,   2 ,   ,   K , to construct the following Lagrangian functional:
L [ p ( x ) ] d x   p ( x ) l n p ( x ) m ( x ) k = 1 K [ λ k p ( x ) F k ( x )   d x F k ] λ 0 [ p ( x ) d x   1 ] .
The critical point of L [ p ( x ) ] is obtained by solving the equation that results from setting the first Gateaux-differential of L [ p ( x ) ] to zero, namely:
δ L [ p ; δ p ] { d d ε L [ p ( x ) + ε   δ p ( x ) ] } ε = 0 = 0 = l n p ( x ) m ( x )   δ p ( x ) d x δ p ( x )   d x k = 1 K λ k δ p ( x ) F k ( x )   d x λ 0 δ p ( x ) d x .
It follows from Equation (96) that:
p ( x ) = m ( x ) e x p [ λ 0 1 k = 1 K λ k F k ( x ) ] .
Replacing the results obtained in Equation (97) into Equation (92) and eliminating the Lagrange multiplier λ 0 from the resulting expression leads to the following expression for the probability density function p ( x ) :
p ( x ) = 1 Z m ( x ) e x p [ k = 1 K λ k F k ( x ) ]
where the normalization constant Z in Equation (98) is defined as follows:
Z d x   m ( x ) e x p [ k = 1 K λ k F k ( x ) ]
In statistical mechanics, the normalization constant Z is called the partition function (or sum over states), and carries all of the information available about the possible states of the system.
The expected integral data is obtained by differentiating Z with respect to the Lagrange multiplier λ k , to obtain the following relationships:
F k = λ k l n Z ,   k = 1 ,   2 ,   ,   K
When the integral data F k are not yet known, the uniform distribution m ( x ) = 1 is the most appropriate to consider. In this case, the maximum entropy (MaxEnt) algorithm yields the uniform distribution as would be required by the principle of insufficient reason. Thus, the MaxEnt principle generalizes the principle of insufficient reason. The MaxEnt can be applied to both discrete and to continuous distributions. The MaxEnt method has been shown [12] to be equivalent to constrained variational inference, thus establishing the link between MaxEnt and Bayesian approximations. The MaxEnt method has been used in many fields; enumerating these fields is beyond the scope of this work, which is limited to nuclear engineering applications. The pioneering application of the MaxEnt method to time-independent nuclear reactor physics problems was initiated in the 1970s [13,14,15,16,17]. The first application of the MaxEnt method to a time-dependent nuclear energy system was by Barhen et al. [18]. This work was subsequently extended by Cacuci and Ionescu-Bujor [19], which presented analytical formulas for predicted mean values and covariance matrices, for both predicted model parameters and responses, which have generalized the previous results presented in data assimilation procedures for geophysical sciences and linear Bayesian models. Other applications to nuclear energy systems are presented in [20,21,22,23,24,25,26,27].
To the author’s knowledge, none of the analytical results published thus far include comprehensively all of the second- and third-order response sensitivities to all of the model’s parameters. The end-results to be presented in this section will extend the results published thus far in the open literature by presenting analytical formulas for both predicted model responses and parameters by including all of the second-and third-order sensitivities of computed model responses to the model’s parameters.

5.1. 2nd/3rd-BERRU-PM: A Priori Information

This Subsection will present the mathematical form of the information that will be ultimately used for predicting the optimal, best-estimate mean values for both the model responses and model parameters, with reduced predicted uncertainties, in the combined parameter-response phase space.

5.1.1. Expected Values and Covariances of Measured Responses

Consider that N r quantities of interest, henceforth called “system responses” and denoted as r i m , i = 1 , , N r , have been experimentally measured, yielding their expected values, as well as their corresponding covariances (i.e., standard deviations and correlations). For the mathematical derivations to follow, it is convenient to consider the responses r i m to constitute the components of the N r -dimensional column vector r m defined as follows:
r m ( r 1 m , , r N r m )
Since the experimentally measured responses cannot be measured exactly, they are usually considered to be variates that follow an unknown multivariate distribution function of the observations, denoted as p r ( r ) , which is formally defined on a domain D r . Methods for finding estimates of a measured quantity and indicators of the quality of the estimates depend on the assumed form of the unknown distribution function p r ( r ) of the observations. The moments of the distribution of measured responses can be conveniently denoted by introducing the following notation for the expectation (or mean value) of a function w ( r m ) :
E [ w ( r m ) ] D r w ( r m ) p r ( r ) d r
When the distribution is discrete, the integral in Equation (102) denotes a sum over the respective discrete probabilities. Using the notation introduced in Equation (102), the expectation (or mean value) of the experimentally measured responses r i m is formally defined as follows:
E ( r i m ) D r r i m   p r ( r ) d r   ,   i = 1 , , N r ,
The expected values of the measured responses will be considered to constitute the components of the vector r m defined as:
E ( r m )   [ E ( r 1 m ) , , E ( r N r m ) ]
The covariance matrix of measured responses will be denoted as C m , and is defined as:
C m   [ r m E ( r m ) ] [ r m E ( r m ) ] r = [ c o v ( r i m , r j m ) ] N r × N r   , [ c o v ( r i m , r j m ) ] N r × N r [ r i m E ( r i m ) ] [ r j m E ( r j m ) ] r ;   i , j = 1 , , N r   .

5.1.2. Expectations and Covariances of Computed Responses Including Second- and Third-Order Sensitivities to Model Parameters

For subsequent computations, it is convenient to consider that the expected values α i 0 of the components α i of the N α -dimensional vector of model parameters α ( α 1 , , α N α ) are components of the following vector of mean (expected) values:
α 0 ( α 1 0 , , α N α 0 )
The variances and covariances defined in Equations (82) and (83) are considered to constitute the elements of a symmetric, positive-definite parameter covariance matrix of dimension N α × N α , denoted as C α and defined as follows:
C α ( α α 0 ) ( α α 0 ) α .
Consider that the values of the N r experimentally measured responses r i m can be computed using a multi-physics model that comprises N α model parameters α n ,   n = 1 , , N α which are related to the model’s independent and dependent variables through the model’s underlying equations, correlations, tables, etc. Of course, the computed response values will not coincide with the measured ones, because, just like the experimentally measured responses, the model’s parameters and numerical solution of the underlying equations, and consequently the computed response values, are also subject to uncertainties. The computed responses r k c ( α ) ,   k = 1 , , N r , are considered to be elements of an N r -dimensional vector r c ( α ) , defined as follows:
r c ( α ) [ r 1 c ( α ) , , r N r c ( α ) ]
The expectation values E [ r k c ( α ) ] given by Equation (87) are considered to be the components of the following vector of “expected values of the computed response”:
E [ r c ( α ) ]   [ E ( r 1 c ) , , E ( r N r c ) ]
The response covariances defined in Equation (88) are considered to be the components of a ( N r × N r ) -dimensional matrix denoted as C r and defined as follows:
C r [ r c E ( r c ) ] [ r c E ( r c ) ] α .
The covariances between the computed responses and the model –parameters defined in Equation (89) are considered to be the components of an ( N r × N α ) -dimensional matrix denoted as C r α and defined as follows:
C r α [ r c E ( r c ) ] ( α α 0 ) α = C α r
The joint covariance matrix, denoted as C M , of the model parameters and model-computed responses is defined as follows:
C M ( C α C α r C r α C r ) ,

5.2. 2nd/3rd-BERRU-PM: Analytical Expressions for Best-Estimate Results with Reduced Uncertainties for Responses and Parameters in the Joint Phase-Space of Responses and Parameters

Consider the joint probability function p ( α , r ) of the multi-variates α and r m , which is defined on the domain D D α × D r and is properly normalized such that:
D p ( α , r )   d α   d r = 1
The exact form of p ( α , r ) is unknown, of course. Since the (multi)variates α and r m are statistically independent of each other, it follows that p ( α , r ) = p r ( r ) p α ( α ) . Therefore, the expected value of a function w ( r m ) satisfies the following relations:
w ( r m ) D w ( r m ) p ( α , r )   d α   d r = D w ( r m ) p r ( r ) p α ( α )   d α   d r = [ D α p α ( α ) d α ] [ D r w ( r m ) p r ( r ) d r ] = w ( r m ) r
while the expected value of a function u ( α ) satisfies the following relations:
u ( α ) D u ( α ) p ( α , r )   d α   d r = D u ( α ) p r ( r ) p α ( α )   d α   d r = [ D r p r ( r ) d r ] [ D α u ( α ) p α ( α ) d α ] = u ( α ) α
Therefore, the a priori information about the model parameters, computed and measured responses can be conveniently summarized by considering that the physical system under consideration is described mathematically by a multivariate vector:
x [ α , r c ( α ) , r m ]
obeying an unknown joint multivariate distribution function:
p ( α , r ) = p r ( r ) p α ( α )
but having a known vector of expected values denoted as:
x 0 [ α 0 , E ( r c ) , E ( r m ) ]
and a covariance matrix denoted as:
C ( C α C α r 0 C r α C r 0 0 0 C m )
Applying the MaxEnt principle as described in the Appendix A to the information provided in Equations (116) through (119) indicates that the MaxEnt form, denoted as p 2 ( x | x 0 ,   C ) , of the unknown distribution p ( α , r ) will have the following multivariate Gaussian form:
p 2 ( x | x 0 ,   C )   d x   =   e x p [ 1 2 ( x x 0 ) C 1 ( x x 0 ) ] d x , d e t ( 2 π C ) ,   < x j < .
The MaxEnt-Gaussian form shown in Equation (120) can also be written in the equivalent form:
p 2 ( x | x 0 ,   C ) =   e x p [ Q ( α , r c , r m ) ] d e t ( 2 π C m ) d e t [ 2 π ( C α C α r C r α C r ) ] , Q ( α , r c , r m ) 1 2 ( α α 0 r c E ( r c ) ) ( C α C α r C r α C r ) 1 ( α α 0 r c E ( r c ) ) 1 2 [ r m E ( r m ) ] C m 1 [ r m E ( r m ) ] .
The expression provided in Equation (121) highlights the “Bayesian” construction underlying the expression of the a posteriori joint maximum-entropy provability distribution function p 2 ( x | x 0 ,   C ) of the physical system’s responses and parameters. In order for the measured and computed responses to represent the same physical quantity, it is necessary that:
r m = r c ( α ) = r ,
where the vector r represents both the computed and measured responses. Determining the moments of p 2 ( x | x 0 ,   C ) for subsequent predictions will require evaluations of integrals of the following form:
E [ g ( α , r ) ] = D g ( α , r ) e x p [ h ( α , r ) ]   d α   d r D e x p [ h ( α , r ) ]   d α   d r
Expressions such as shown in Equation (123) can be evaluated to a high degree of accuracy (a priory controlled) by using the saddle-point method (also called Laplace approximation, or steepest descent method), which relies on evaluating the respective integrals at the so-called “saddle point(s).” For the integral in the denominator of Equation (123), the saddle point, denoted as x D ( α D , r D ) , is the point at which the gradient of the function h ( x ) vanishes, i.e., x h ( x D ) = 0 , so that h ( x ) can be expanded in the following Taylor series:
h ( x ) = h ( x D ) + S O T ( x ) ,
where the quantity " S O T ( x ) " denotes terms of second- and higher-order in x . When the function g ( x ) varies slowly, it is simply evaluated at the respective saddle point, and the resulting expression of E [ g ( α , r ) ] in Equation (123) becomes:
E [ g ( α , r ) ] = g ( α D , r D ) e x p [ h ( α D , r D ) ] D e x p [ S O T ( α , r ) ]   d α   d r e x p [ h ( α D , r D ) ] D e x p [ S O T ( α , r ) ]   d α   d r g ( α D , r D ) ,
where the saddle-point ( α D , r D ) is defined as the point in phase-space where the gradients of h ( α , r ) vanish, i.e.,
{ α h ( α , r ) } ( α D , r D ) = 0 ,   { r h ( α , r ) } ( α D , r D ) = 0
When the function h ( x ) has a Taylor-series containing powers higher than second-order in x , the gradient x h ( x ) of the function h ( x ) may vanish at multiple saddle points, in which case the contributions from all of these saddle points would need to be accounted for when evaluating the integrals in Equation (123).

5.2.1. Predicted Best-Estimate Expected Values for the Responses and Parameters in the Joint Phase-Space of Responses and Parameters

The MaxEnt Gaussian p 2 ( x | x 0 ,   C ) has a quadratic form and hence possesses a single saddle point, which is determined by requiring the first-variation δ Q ( α D   ,   r D ; δ α , δ r ) of the exponential term in Equation (121) to vanish at the saddle point x D ( α D , r D ) , namely:
δ Q ( α D   ,   r D ; δ α , δ r ) { d d ε Q ( α D + ε δ α , r D + ε δ r ) } ε = 0 = Q α δ α + Q r δ r = 0   .
Imposing the requirement indicated in Equation (127) while taking Equation (122) into account yields the following system of equations:
α [ 1 2 ( α α 0 r E ( r c ) ) ( C α C α r C r α C r ) 1 ( α α 0 r E ( r c ) ) 1 2 [ r E ( r m ) ] C m 1 [ r E ( r m ) ] ] = 0 ,
r [ 1 2 ( α α 0 r E ( r c ) ) ( C α C α r C r α C r ) 1 ( α α 0 r E ( r c ) ) 1 2 [ r E ( r m ) ] C m 1 [ r E ( r m ) ] ] =   0 .
Recalling that:
( x A x ) x = x ( A + A )
and applying the result shown in Equation (130) to evaluate the expressions provided in Equations (128) and (129), respectively, yields the following system of equations at the saddle point, denoted as ( α D , r D ) :
( C α C α r C r α C r ) 1 ( α D   α 0 r D E ( r c ) ) + ( 0 0 0 C m 1 ) ( α D   α 0 r D E ( r c ) ) +   ( 0 C m 1 [ E ( r c ) E ( r m ) ] ) =   ( 0 0 ) .
Solving Equation (131) yields the following intermediate expressions:
( α D α 0 + C α r C m 1 [ r D E ( r m ) ] r D E ( r c ) + C r C m 1 [ r D E ( r m ) ] ) =   ( 0 0 ) .
For computational purposes, the results obtained in Equation (132) are recast into the following forms:
r D = E ( r m ) + C m ( C m + C r ) 1 [ E ( r c ) E ( r m ) ] ,
α D = α 0 C α r ( C m + C r ) 1 [ E ( r c ) E ( r m ) ] .
In view of Equation (125), it follows that “best-estimate predicted” results for the posterior expectation values for the responses and parameters, which will be denoted as ( r b e , α b e ) , will have the same expressions as shown in Equations (133) and (134), namely:
r b e D r p 2 ( x | x 0 ,   C )   d α   d r = r D = r m + C m ( C m + C r ) 1 [ E ( r c ) E ( r m ) ] ,
α b e D α p 2 ( x | x 0 ,   C )   d α   d r = α D = α 0 C α r ( C m + C r ) 1 [ E ( r c ) E ( r m ) ] .
Since the components of the vector E ( r c ) , and the components of the matrices C r and C α r contain 2nd-order and 3rd-order sensitivities, the formulas presented in Equations (135) and (136) generalize all of the previous formulas of this type found in data assimilation/assimilation procedures published to date (which contain at most first-order sensitivities).

5.2.2. Predicted Best-Estimate Covariances for the Responses and Parameters in the Joint Phase-Space of Responses and Parameters

The second-order moments of the posterior distribution p ( α , r ) = p r ( r ) p α ( α ) comprise the covariances between the best estimated response, which are denoted as C r b e , the covariances between the best-estimate parameters, which are denoted as C α b e , and the covariances between the best-estimate parameters and responses., which are denoted as C α r b e . The expression of the “best-estimate” posterior parameter covariance matrix C r b e for the best-estimate responses r b e is derived by using the results given in Equations (133) and (135) to obtain:
C r b e D ( r r b e ) ( r r b e ) p ( α , r )   d α   d r = C m C m ( C m + C r ) 1 C m   = C m [ I ( C m + C r ) 1 C m ] .
As indicated in Equation (137), the initial covariance matrix C m is multiplied by the matrix [ I ( C m + C r ) 1 C m ] , which means that the variances contained on the diagonal of the best-estimate matrix C r b e will be smaller than the experimentally measured variances contained in C m . Hence, the addition of new experimental information has reduced the predicted best-estimate response variances in C r b e by comparison to the measured variances contained a priori in C m . Since the components of the matrix C r contain 2nd-order and 3rd-order sensitivities, the formula presented in Equation (137) generalizes all of the previous formulas of this type found in data assimilation/assimilation procedures published to date (which contain at most first-order sensitivities).
The expression of the “best-estimate” posterior parameter covariance matrix C α b e for the best-estimate parameters α b e is derived by using the result given in Equation (136) to obtain:
C α b e D ( α α b e ) ( α α b e ) p ( α , r )   d α   d r = C α C α r ( C m + C r ) 1 C r α .
Both matrices C α and C α r ( C m + C r ) 1 C r α are symmetric and positive definite. Therefore, the subtraction indicated in Equation (138) implies that the components of the main diagonal of C α b e must have smaller values than the corresponding elements of the main diagonal of C α . In this sense, the introduction of new computational and experimental information has reduced the best-estimate parameter variances on the diagonal of C α b e . Since the components of the matrices C α , C α r , and C r contain 2nd-order and 3rd-order sensitivities, the formula presented in Equation (138) generalizes all of the previous formulas of this type found in data assimilation/assimilation procedures published to date (which contain at most first-order sensitivities).
The expression of the “best-estimate” posterior parameter covariance matrix C α b e for the best-estimate parameters α b e and best-estimate responses r b e is derived by using the results given in Equations (135) and (136) to obtain:
C α r b e D ( α α b e ) ( r r b e ) p ( α , r )   d α   d r = C α r ( C m + C r ) 1 C m   .
As before, since the components of the matrices C α r , and C r contain 2nd-order and 3rd-order sensitivities, the formula presented in Equation (139) generalizes all of the previous formulas of this type found in data assimilation/assimilation procedures published to date (which contain at most first-order sensitivities).
The expression of the best-estimate covariance matrix C r α b e is derived by performing a sequence of operations similar to that shown in Equation (139) to obtain:
C r α b e D ( r r b e ) ( α α b e ) p ( α , r )   d α   d r = C m ( C m + C r ) 1 C r α = ( C α r b e ) .
It is important to note from the results shown in Equations (135) through (140) that the computation of the best estimate parameter and response values, together with their corresponding best-estimate covariance matrices, only requires the computation of ( C m + C r ) 1 , which entails the inversion of a matrix of size N r × N r . This is computationally very advantageous, since N r N α (i.e., the number of responses is much less than the number of model parameters) in the overwhelming majority of practical situations.

5.2.3. Data Consistency Indicator

As will be shown in the following, the minimum value, Q min , of the functional Q ( α b e , r b e ) in the exponential in Equation (121), provides a “consistency indicator” which quantifies the mutual and joint consistency of the information available for model calibration. The minimum value, Q m i n , of the functional Q ( α b e , r b e ) in the exponential in Equation (121) has the following expression:
Q min Q ( α b e , r b e ) = 1 2 ( α b e α 0 r b e E ( r c ) ) ( C α C α r C r α C r ) 1 ( α b e α 0 r b e E ( r c ) )   1 2 [ r b e E ( r m ) ] C m 1 [ r b e E ( r m ) ] .
For subsequent matrix algebra, it is convenient to use the following matrix:
( B 11 B 12 B 21 B 22 ) ( C α C α r C r α C r ) 1 .
Using Equations (142) and (132) in Equation (141) yields the following result:
Q m i n = 1 2 [ E ( r c ) E ( r m ) ] ( C m + C r ) 1 × [ C r α ( B 11 C α r + B 12 C r ) + C r ( B 21 C α r + B 22 C r ) + C m ] ( C m + C r ) 1 [ E ( r c ) E ( r m ) ]   = 1 2 d ( C m + C r ) 1 d .
where
d [ E ( r c ) E ( r m ) ]
As the expression obtained in Equation (143) indicates, the quantity Q m i n represents the square of the length of the vector d [ E ( r c ) E ( r m ) ] , measuring (in the corresponding metric) the deviations between the experimental and nominally computed responses. The quantity Q m i n can be evaluated directly from the given data (i.e., given parameters and responses, together with their original uncertainties) after having inverted the covariance matrix ( C m + C r ) . It is also important to note that Q m i n is independent of calibrating (or adjusting) the original data. As the dimension of [ E ( r c ) E ( r m ) ] indicates, the number of degrees of freedom characteristic of the calibration under consideration is equal to the number N r of experimental responses. In the extreme case of absence of experimental responses, no actual calibration takes place. An actual calibration (adjustment) occurs only when including at least one experimental response.
It is noteworthy that the expressions presented in Equations (135) through (139) and (143) look superficially similar to the expressions presented in [7,19]. This superficial similarity stems from the fact that the 2nd/3rd-BERRU-PM formalism only includes means and covariances of the model parameters and computed and measured responses. However, contrary to previously published works, the 2nd/3rd-BERRU-PM formalism presented in in Equations (135) through (139) and (143) comprise all of the 2nd-order and 3rd-order sensitivities of the computed responses with respect to the model parameters. None of these 2nd- and 3rd-order sensitivities appear in the works previous to this work. Consequently, the 2nd/3rd-BERRU-PM formalism comprises, as particular cases, the previously published formulas used in Kalman filters, Bayesian linear statistics [10], data adjustment [13,14,15,16,17], data assimilation [8,9] and predictive modeling [7,18,19,20,21,22,23,24,25,26,27].

6. Conclusions

This work has presented the Third-Order Adjoint Sensitivity Analysis Methodology (3rd-ASAM), which enables the efficient computation of the exact expressions of the 3rd-order functional derivatives (“sensitivities”) of a general system response that depends on both the forward and adjoint state functions, with respect to all of the parameters underlying the respective forward and adjoint systems. Such responses are often encountered when representing mathematically detector responses and reaction rates in reactor physics problems. The 3rd-ASAM extends the 2nd-ASAM in the quest to overcome the “curse of dimensionality” in sensitivity analysis, uncertainty quantification and predictive modeling.
Very importantly, the computation of the 2nd-level adjoint functions ψ 1 , j ( 2 ) ( x ) , ψ 2 , j ( 2 ) ( x ) , ψ 3 , j ( 2 ) ( x ) , ψ 4 , j ( 2 ) ( x ) , and of the 3rd-level adjoint functions, ψ 1 , i j ( 3 ) , , ψ 8 , i j ( 3 ) , is performed by using the same forward and adjoint solvers (i.e., computer codes) as used for solving the original forward and adjoint systems, namely Equations (1) and (6) subject to the corresponding boundary conditions. Thus, solving the 2nd-LASS and 3rd-ASAM would not require any significant “code development,” since the original forward and adjoint solvers (codes) would not need to be modified; only the right-sides (i.e., “sources”) for these solvers/codes would need to be programmed accordingly. Of course, if the response depends only on the original forward or original adjoint function, than only half of the equations underlying the 2nd-ASAM and, correspondingly, the 3rd-ASAM will need to be solved.
This work also presents new formulas that incorporate the contributions of the 3rd-order sensitivities into the expressions of the first four cumulants of the response distribution in the phase-space of model parameters. Using these newly developed formulas, this work also presents a new mathematical formalism, called the 2nd/3rd-BERRU-PM (“Second/Third-Order Best-Estimated Results with Reduced Uncertainties Predictive Modeling”), which combines experimental and computational information in the joint phase-space of responses and model parameters, including not only the 1st-order response sensitivities, but also the complete hessian matrix of 2nd-order second-sensitivities and also the 3rd-order sensitivities, all computed using the 3rd-ASAM. The 2nd/3rd-BERRU-PM formalism uses the maximum entropy principle to eliminate the need for introducing and “minimizing” a user-chosen “cost functional quantifying the discrepancies between measurements and computations,” thus yielding results that are free of subjective user-interferences and generalizing and significantly extending the 4D-VAR data assimilation procedures. Incorporating correlations, including those between the imprecisely known model parameters and computed model responses, the 2nd/3rd-BERRU-PM also provides a quantitative metric, constructed from sensitivity and covariance matrices, for determining the degree of agreement among the various computational and experimental data while eliminating discrepant information. The mathematical framework of the 2nd/3rd-BERRU-PM formalism requires the inversion of a single matrix of size N r × N r , where N r denotes the number of considered responses. In the overwhelming majority of practical situations, the number of responses is much less than the number of model parameters. Thus, the 2nd-BERRU-PM methodology overcomes the curse of dimensionality which affects the inversion of hessian matrices in the parameter space.

Funding

This work was partially funded by the United States National Nuclear Security Administration’s Office of Defense Nuclear Nonproliferation Research & Development.

Conflicts of Interest

The author declares no conflict of interest.

Appendix A

For convenient referencing, this Appendix summarizes the derivation of the MaxEnt approximation of distribution, p ( x ) , when only the first- and second-order moments of this distribution are known. These moments are generally defined as follows:
(i)
The vector x 0 ( x 1 0 , , x N 0 ) of mean values (first-order moments), denoted as x i 0 , of x i , and defined as x 1 0 x i p ( x ) d x ;
(ii)
The second-order moment or covariance, μ 2 i j ( x ) c o v ( x i , x j ) , of two parameters, x i and x j , defined as μ 2 i j c o v ( x i , x j ) ( x i x i 0 ) ( x j x j 0 ) p ( x ) d x ,   i , j = 1 , , N . The covariances μ 2 i j ( x ) constitute the elements of a symmetric, positive-definite parameter covariance matrix of dimension N × N , denoted as C [ μ 2 11 . μ 2 1 N . μ 2 i j . μ 2 N 1 . μ 2 N N ]
Taking m ( x ) = 1 and particularizing the generic constraints represented by Equation (98) yields the following particular probability density function, denoted as p 2 ( x ) , where the subscript “2” indicates “up to second-order moments:” p 2 ( x ) = 1 Z 2 e x p [ i = 1 N a i x i 1 2 i = 1 N j = 1 N b i j ( x i x i 0 ) ( x j x j 0 ) ] , where the normalization constant Z 2 is defined as Z 2 d x   e x p [ i = 1 N a i x i 1 2 i = 1 N j = 1 N b i j ( x i x i 0 ) ( x j x j 0 ) ] . The normalization constant Z 2 can be evaluated in closed form as follows: introduce the N-dimensional vectors z ( z 1 , , z N ) x x 0 ( x 1 x 1 0 , , x N x N 0 ) , a ( a 1 , , a N ) and the N × N -dimensional symmetric matrix B ( b i j ) to transform the expression in the argument of the exponential as follows:
i = 1 N a i x i 1 2 i = 1 N j = 1 N b i j ( x i x i 0 ) ( x j x j 0 ) = a x 0 a z 1 2 z B z = 1 2 ( z + B 1 a ) B ( z + B 1 a ) + a x 0 + 1 2 a B 1 a   .
Since d x e x p [ 1 2 ( z + B 1 a ) B ( z + B 1 a ) ] = ( 2 π ) N d e t ( B ) , it follows that Z 2 = ( 2 π ) N d e t ( B ) e x p ( a x 0 1 2 a B 1 a ) . The Lagrange multipliers a i and b i j are now obtained as follows:
a i l n Z 2 = x i 0 + ( B 1 a ) i = x i 0 ;     a i = 0 ,   i = 1 ,   2 ,   ,   N . b i j l n Z 2 = 1 2 [ d e t ( B ) ] 1 [ d e t ( B ) ] b i j = c o f a c t o r ( b i j ) d e t ( B ) = [ B 1 ] i j = μ 2 i j ,   i , j = 1 ,   2 ,   ,   N ;     B 1 = C .
Inserting the above results into the expression of p 2 ( x ) yields:
p 2 ( x | x ,   C )   d x   =   e x p [ 1 2 ( x x 0 ) C 1 ( x x 0 ) ] d x , ( 2 π ) N d e t ( C ) ,   < x j < .
The above derivation demonstrates that, when only means and covariances are known, the maximum entropy algorithm yields the Gaussian probability distribution p 2 ( x | x ,   C )   as the most objective probability distribution, where: x is the data vector with coordinates x j , C is the covariance matrix with elements μ 2 i j ( x ) c o v ( x i , x j ) , and d x j d x j is the volume element in the data space. It often occurs in practice that the variances μ 2 i i ( x ) are known but the covariances μ 2 i j ( x ) ,   i j , are not known. In this case, the covariance matrix C would a priori be diagonal. Consequently, only the Lagrange parameters b i i would be non-zero, so that the matrix B ( b i j ) would also be a priori diagonal. In other words, in the absence of information about correlations, the maximum entropy algorithm indicates that unknown covariances/correlations can be taken to be zero.

References

  1. Cacuci, D.G. Second-Order Adjoint Sensitivity Analysis Methodology (2nd-ASAM) for Computing Exactly and Efficiently First- and Second-Order Sensitivities in Large-Scale Linear Systems: I. Computational Methodology. J. Comput. Phys. 2015, 284, 687–699. [Google Scholar] [CrossRef]
  2. Cacuci, D.G. Second-Order Adjoint Sensitivity Analysis Methodology (2nd-ASAM) for Large-Scale Nonlinear Systems: I. Theory. Nucl. Sci. Eng. 2016, 184, 16–30. [Google Scholar] [CrossRef]
  3. Cacuci, D.G. The Second-Order Adjoint Sensitivity Analysis Methodology; CRC Press, Taylor & Francis Group: Boca Raton, FL, USA, 2018. [Google Scholar]
  4. Cacuci, D.G.; Fang, R.; Favorite, J.A. Comprehensive Second-Order Adjoint Sensitivity Analysis Methodology (2nd-ASAM) Applied to a Subcritical Experimental Reactor Physics Benchmark: I. Effects of Imprecisely Known Microscopic Total and Capture Cross Sections. Energies 2019, 12, 4100. [Google Scholar] [CrossRef]
  5. Fang, R.; Cacuci, D.G. Comprehensive Second-Order Adjoint Sensitivity Analysis Methodology (2nd-ASAM) Applied to a Subcritical Experimental Reactor Physics Benchmark: II. Effects of Imprecisely Known Microscopic Scattering Cross Sections. Energies 2019, 12, 4114. [Google Scholar] [CrossRef]
  6. Cacuci, D.G.; Fang, R.; Favorite, J.A.; Badea, M.C.; Di Rocco, F. Comprehensive Second-Order Adjoint Sensitivity Analysis Methodology (2nd-ASAM) Applied to a Subcritical Experimental Reactor Physics Benchmark: III. Effects of Imprecisely Known Microscopic Fission Cross Sections and Average Number of Neutrons per Fission. Energies 2019, 12, 4100. [Google Scholar] [CrossRef]
  7. Cacuci, D.G. BERRU Predictive Modeling: Best Estimate Results with Reduced Uncertainties; Springer: Berlin/Heidelberg, Germany; New York, NY, USA, 2018. [Google Scholar]
  8. Lewis, J.M.; Lakshmivarahan, S.; Dhall, S.K. Dynamic Data Assimilation: A Least Square Approach; Cambridge University Press: Cambridge, UK, 2006. [Google Scholar]
  9. Lahoz, W.; Khattatov, B.; Ménard, R. Data Assimilation: Making Sense of Observations; Springer: New York, NY, USA, 2010. [Google Scholar]
  10. Goldstein, M.; Woof, D. Bayes Linear Statistics: Theory and Methods; John Wiley & Sons: Chichester, UK, 2007. [Google Scholar]
  11. Jaynes, E.T. Information Theory and Statistical Mechanics. Phys. Rev. 1957, 106, 620–630. [Google Scholar] [CrossRef]
  12. Granziol, D.; Ru, B.; Zohren, S.; Dong, X.; Osborne, M.; Roberts, S. MEMe: An Accurate Maximum Entropy Method for Efficient Approximations in Large-Scale Machine Learning. Entropy 2019, 21, 551. [Google Scholar] [CrossRef]
  13. Rowlands, J.L.; Dean, C.J.; MacDougall, J.D.; Smith, R.W. The Production and Performance of the Adjusted Cross-Section Set FGL5. In Proceedings of the International Symposium “Physics of Fast Reactors”, Tokio, Japan, 16–23 October 1973. [Google Scholar]
  14. Gandini, A.; Petilli, M. AMARA: A Code Using the Lagrange Multipliers Method for Nuclear Data Adjustment; CNEN-RI/FI(73)39; Comitato Nazionale Energia Nucleare: Rome, Italy, 1973.
  15. Kuroi, H.; Mitani, H. Adjustment to Cross-Section Data to Fit Integral Experiments by Least Squares Method. J. Nucl. Sci. Technol. 1975, 12, 663. [Google Scholar] [CrossRef]
  16. Dragt, J.B.; Dekker, J.W.M.; Gruppelaar, H.; Janssen, A.J. Methods of Adjustment and Error Evaluation of Neutron Capture Cross Sections. Nucl. Sci. Eng. 1977, 62, 11. [Google Scholar] [CrossRef]
  17. Weisbin, C.R.; Oblow, E.M.; Marable, J.H.; Peelle, R.W.; Lucius, J.L. Application of Sensitivity and Uncertainty Methodology to Fast Reactor Integral Experiment Analysis. Nucl. Sci. Eng. 1978, 66, 307. [Google Scholar] [CrossRef]
  18. Barhen, J.; Cacuci, D.G.; Wagschal, J.J.; Bjerke, M.A.; Mullins, C.B. Uncertainty Analysis of Time-Dependent Nonlinear Systems: Theory and Application to Transient Thermal Hydraulics. Nucl. Sci. Eng. 1982, 81, 23–44. [Google Scholar] [CrossRef]
  19. Cacuci, D.G.; Ionescu-Bujor, M. Sensitivity and Uncertainty Analysis, Data Assimilation and Predictive Best-Estimate Model Calibration. In Handbook of Nuclear Engineering; Cacuci, D.G., Ed.; Springer: New York, NY, USA; Berlin/Heidelberg, Germnay, 2010; Volume 3, Chapter 17; pp. 1913–2051. ISBN 978-0-387-98150-5. [Google Scholar]
  20. Badea, M.C.; Cacuci, D.G.; Badea, A.F. Best-Estimate Predictions and Model Calibration for Reactor Thermal-Hydraulics. Nucl. Sci. Eng. 2012, 172, 1–19. [Google Scholar] [CrossRef]
  21. Arslan, E.; Cacuci, D.G. Predictive Modeling of Liquid-Sodium Thermal-Hydraulics Experiments and Computations. Ann. Nucl. Energy 2014, 63, 355–370. [Google Scholar] [CrossRef]
  22. Cacuci, D.G.; Arslan, E. Reducing Uncertainties via Predictive Modeling: FLICA4 Calibration Using BFBT Benchmarks. Nucl. Sci. Eng. 2014, 176, 339–349. [Google Scholar] [CrossRef]
  23. Peltz, J.J.; Cacuci, D.G.; Badea, A.F.; Badea, M.C. Predictive Modeling Applied to a Spent Fuel Dissolver Model: II. Uncertainty Quantification and Reduction. Nucl. Sci. Eng. 2016, 183, 332–346. [Google Scholar] [CrossRef]
  24. Badea, A.F.; Cacuci, D.G. Predictive Uncertainty Reduction in Coupled Neutron-Kinetics/Thermal Hydraulics Modeling of the BWR-TT2 Benchmark. Nucl. Eng. Des. 2017, 313, 330–344. [Google Scholar] [CrossRef]
  25. Peltz, J.J.; Cacuci, D.G. Inverse Predictive Modeling of a Spent Fuel Dissolver Model. Nucl. Sci. Eng. 2016, 184, 1–15. [Google Scholar] [CrossRef]
  26. Cacuci, D.G. Inverse predictive modeling of radiation transport through optically thick media in the presence of counting uncertainties. Nucl. Sci. Eng. 2017, 186, 199–223. [Google Scholar] [CrossRef]
  27. Cacuci, D.G.; Fang, R.; Badea, M.C. MULTI-PRED: A Software Module for Predictive Modeling of Coupled Multi-Physics Systems. Nucl. Sci. Eng. 2018, 191, 187–202. [Google Scholar] [CrossRef]

Share and Cite

MDPI and ACS Style

Cacuci, D.G. Towards Overcoming the Curse of Dimensionality: The Third-Order Adjoint Method for Sensitivity Analysis of Response-Coupled Linear Forward/Adjoint Systems, with Applications to Uncertainty Quantification and Predictive Modeling. Energies 2019, 12, 4216. https://doi.org/10.3390/en12214216

AMA Style

Cacuci DG. Towards Overcoming the Curse of Dimensionality: The Third-Order Adjoint Method for Sensitivity Analysis of Response-Coupled Linear Forward/Adjoint Systems, with Applications to Uncertainty Quantification and Predictive Modeling. Energies. 2019; 12(21):4216. https://doi.org/10.3390/en12214216

Chicago/Turabian Style

Cacuci, Dan Gabriel. 2019. "Towards Overcoming the Curse of Dimensionality: The Third-Order Adjoint Method for Sensitivity Analysis of Response-Coupled Linear Forward/Adjoint Systems, with Applications to Uncertainty Quantification and Predictive Modeling" Energies 12, no. 21: 4216. https://doi.org/10.3390/en12214216

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop