Probability-space surrogate modeling for fast multidisciplinary optimization under uncertainty

Journal Pre-proof

Probability-space surrogate modeling for fast multidisciplinary optimization under uncertainty Saideep Nannapaneni , Sankaran Mahadevan PII: DOI: Reference:

S0951-8320(19)30321-7 https://doi.org/10.1016/j.ress.2020.106896 RESS 106896

To appear in:

Reliability Engineering and System Safety

Please cite this article as: Saideep Nannapaneni , Sankaran Mahadevan , Probability-space surrogate modeling for fast multidisciplinary optimization under uncertainty, Reliability Engineering and System Safety (2020), doi: https://doi.org/10.1016/j.ress.2020.106896

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Published by Elsevier Ltd.

Highlights     

Probability-space surrogate models with analytical solutions Surrogate training using single iterations of multidisciplinary analysis Simultaneous uncertainty and multidisciplinary analyses using Bayesian strategy Computationally efficient reliability-based and robustness-based optimization Illustration using a three-discipline aircraft design example

Probability-space surrogate modeling for fast multidisciplinary optimization under uncertainty Saideep Nannapaneni1 *, Sankaran Mahadevan2 1

Department of Industrial, Systems, and Manufacturing Engineering, Wichita State University, Wichita, KS, 67260, USA 2

Department of Civil & Environmental Engineering, Vanderbilt University, Nashville, TN, 37235, USA

Abstract This paper proposes a probability-space surrogate modeling approach for computationally efficient multidisciplinary design optimization under uncertainty. This paper uses a probabilityspace surrogate as opposed to an algebraic surrogate so that the probability distributions of the required outputs at a given design input can naturally be obtained without repeated Monte Carlo runs of an algebraic surrogate at different realizations of the uncertain variables. We consider three probability-space surrogates with analytical solutions for prediction and inference - Multivariate Gaussian, Gaussian Copula, and Gaussian Mixture Model, and investigate their applicability to perform multidisciplinary design optimization under uncertainty. All the input design and random variables, coupling variables, objective and constraint functions are incorporated within the probability-space surrogate, which helps analytically obtain the distributions of coupling variables, objective and constraint functions at desired design inputs while enforcing multidisciplinary compatibility. The training points for the probability-space surrogates are obtained by performing one-pass analysis through the disciplinary models at different realizations of the input variables. The proposed methodology is demonstrated for reliability-based design optimization (RBDO) and

*

Corresponding author. Email: [email protected], Tel.: +1 316-978-6240

2

reliability-based robust design optimization (RBRDO) in an aircraft design example. The performance of the probability-space surrogates is compared against a Kriging algebraic surrogate and fully coupled Monte Carlo analysis. Keywords: Multidisciplinary, Multi-physics, Optimization, Probability-space surrogate, Gaussian mixture, Gaussian copula, Coupling, Uncertainty, Bayesian network Notation FEA

Finite element analysis

CFD

Computational fluid dynamics

RBDO

Reliability-based design optimization

RDO

Robust design optimization

RBRDO

Reliability-based robust design optimization

MDOUU

Multidisciplinary design optimization under uncertainty

BN

Bayesian network

MCMC

Markov Chain Monte Carlo

MG

Multivariate Gaussian

GC

Gaussian Copula

GMM

Gaussian Mixture Model

GP

Gaussian Process

FORM

First order reliability method

FOSM

First order second moment

MDO

Multidisciplinary design optimization

LAMDA

Likelihood-based approach for multidisciplinary analysis

3

DIRECT

DIviding RECTangles (optimization algorithm)

AR

Wing Aspect Ratio

SW

Wing area

SWEEP

Wing sweep

TOGW

Gross Weight

FPRD

Fan Pressure Ratio

OPRD

Overall Pressure Ratio

FRWI

Wing Weight Factor

FRFU

Fuselage Weight Factor

FCDI

Induced drag factor

FCDO

Profile drag factor

𝜂𝐹𝐴𝑁

Fan design efficiency curve factor

𝜂𝐻𝑃𝐶

High pressure compressor design efficiency curve factor

TOFL

Takeoff field length

𝑽𝒂𝒑𝒑

Approach speed

DNAC

Nacelle Max Diameter

XNAC

Nacelle Length

𝐶𝑟𝑠𝑊𝑡

Cruise weight

𝐶𝑟𝑠𝐷𝑟𝑎𝑔

Cruise drag

𝐶𝑟𝑠𝐿𝑜𝐷

Cruise lift over drag

𝑊𝑒𝑛𝑔

Weight of engine

BF

Block Fuel

OEW

Operating empty weight

4

𝐶𝑟𝑠𝑆𝐹𝐶

Cruise specific fuel consumption

𝑿𝑫

A vector of deterministic design inputs

𝑿𝑼

A vector of inputs associated with aleatory uncertainty

𝑿𝑫,𝑼

A vector of design variables associated with aleatory uncertainty

𝐴𝑖

𝑖 𝑡ℎ disciplinary system in the multi-physics system

𝒇

A vector of system-level objective functions

𝒄

A vector of system-level constraints

𝒈𝒊

Output from the 𝑖 𝑡ℎ disciplinary system analysis

𝒖𝒊𝒋

A vector of coupling variables from 𝑖 𝑡ℎ to 𝑗𝑡ℎ disciplinary systems

𝒖𝒊𝒋,𝒊𝒏

A vector of coupling variables that are inputs to the 𝑗𝑡ℎ disciplinary system

𝒖𝒊𝒋,𝒐𝒖𝒕

A vector of coupling variables that are outputs from the 𝑖 𝑡ℎ disciplinary system

𝑫𝒊𝒋

The difference between 𝒖𝒊𝒋,𝒐𝒖𝒕 and 𝒖𝒊𝒋.𝒊𝒏

𝑿𝒐𝒃𝒔

A vector of observable variables

𝑿𝒖𝒏𝒐𝒃𝒔

A vector of unobservable variables

𝒙𝑫 𝒐𝒃𝒔

Data available on observable variables

𝝁

Expectation vector

𝚺

Covariance matrix

I. Introduction Multidisciplinary systems, which are characterized by coupled interactions between multiple physics disciplines, are often encountered in engineering applications such as hypersonic aircraft [1] (aero, structural and thermal interactions), and long-span bridges [2] (aero-structure

5

interaction). Design optimization of such coupled multi-physics systems is computationally expensive, as it requires a nested three-loop analysis (Fig. 1), with the multidisciplinary convergence analysis in the innermost loop, uncertainty and/or reliability analysis in the middle loop, and optimization analysis in the outermost loop. Both the multidisciplinary analysis and optimization require iterative analyses, and uncertainty or reliability analysis requires many samples; thus, computational effort is a prominent challenge for multidisciplinary optimization under uncertainty.

Fig. 1. Nested three-loop analysis for multidisciplinary design under uncertainty

One approach to reduce the large computational expense is to replace the expensive disciplinary computational models with inexpensive surrogate models. A variety of surrogate modeling techniques have been studied, such as polynomial chaos models, Kriging (or Gaussian Process), and neural networks [3–7]. In this paper, these types of surrogates are termed as physicsspace algebraic surrogates as they seek to provide output value corresponding to an individual realization of the design inputs and uncertain variables. Another approach to reduce the computational effort is to decouple the bi-directional coupling between individual disciplines in multidisciplinary analysis [8–11]. In the decoupled approaches, the coupling variables have been estimated using first-order methods such as First-Order Reliability Methods (FORM) and First-

6

Order Second Moment (FOSM) [12], and considered as additional inputs in the uni-directional analysis. A brief review of the decoupled approaches is presented in Section II. B. The existing methods for multidisciplinary optimization under uncertainty (MDOUU) have several shortcomings. Repeated runs of the surrogate model at the same design values with multiple realizations of the uncertainty variables are required in order to compute the reliability or robustness constraints (or probabilistic objectives). A drawback in the existing decoupled approaches is that they lose the statistical relationship between the design inputs and the coupling variables, i.e., the distribution of the coupling variable remains the same irrespective of the values of the design input variables. The distributions of the coupling variables are dependent on the design and uncertain inputs. The distributions of the individual disciplinary outputs are dependent on the values of the design inputs, the distributions of the uncertain inputs, and the distributions of the coupling variables. Furthermore, the overall distributions of the performance functions (objective and constraint functions) are dependent on the distributions of individual disciplinary outputs. Incorrect estimation of the joint distribution of the coupling variables (due to losing their statistical dependence) can lead to incorrect distributions of the performance functions, which will impact the optimization results. Thus, there exists a need for the development of an MDOUU approach that is computationally efficient yet overcomes the approximations introduced by first-order estimates of coupling variable distributions and reliability/robustness, and is scalable to coupled problems with multiple disciplines. This paper also uses a surrogate-based methodology but a key distinction with the existing approaches is the use of a probability-space surrogate (or distribution surrogate [13,14]) as opposed to an algebraic surrogate in the physical space. An algebraic surrogate seeks to provide the value of the output for a given value of input, whereas a probability-space surrogate (such as a

7

Bayesian network) provides a distribution of the output for a given probability distribution of the inputs. In other words, the probability-space surrogate is constructed in the probability space whereas the algebraic surrogate is constructed in the variable space. A Bayesian network (BN) is a probability-space surrogate where the joint probability distribution is specified as a combination of several marginal and conditional probability distributions. When developing algebraic surrogates, the inputs for such surrogates are the design variables and uncertain variables. To estimate the uncertainty associated with the output prediction at a given input, the algebraic surrogate needs to be evaluated multiple times at the same design input for multiple realizations of the uncertain variables. In contrast, a single evaluation of the probability-space surrogate provides the entire output distribution considering all the uncertain variables at a given value of the input variable [14]. In a generic BN, prediction (i.e., forward propagation) and inference (i.e., backward or inverse propagation) do not have analytical solutions; instead, sampling-based techniques such as Markov Chain Monte Carlo (MCMC) methods [15] need to be used, which are computationally expensive and approximate. This paper investigates three types of multivariate distribution models (that have analytical solutions for the forward and inverse problems) as approximations to a BN. The three models are: (1) Multivariate Gaussian (MG), (2) Gaussian Copula (GC), and (3) Gaussian Mixture Model (GMM). Whenever appropriate, the original BN may be replaced by one of these three probability-space surrogates for computationally efficient MDOUU. Copulas have been used in reliability analysis and RBDO for correlated random inputs [16]. Liang and Mahadevan [17] used a Gaussian copula approach for single objective reliability-based MDO and multi-objective RBDO [18]. A mixture model represents the joint distribution of variables through a weighted combination of individual multivariate distributions; this is

8

particularly useful for multi-modal distributions. When the individual components are modeled using Gaussian distributions, then the mixture model is referred to as a Gaussian Mixture Model. Mixture models have been studied in reliability analysis [19], traffic flow forecasting [20], process monitoring [21], uncertainty quantification in dynamic systems [22] and probabilistic community discovery [23]. In this paper, we exploit the availability of analytical solutions in MG, GC and GMM models for efficient multidisciplinary optimization. We also develop an efficient one-pass analysis approach to generate training points for constructing the probability-space surrogates and perform MDOUU. We demonstrate the one-pass analysis approach for a conceptual two-discipline system, three-discipline system, and a three-discipline system with both bi-directional and unidirectional coupling between individual disciplines. The first contribution of the paper is to investigate a probability-space surrogate modeling approach to carry out uncertainty and reliability analysis in multi-physics systems. The second contribution is to use the probability-space surrogate for multidisciplinary design optimization under uncertainty. In addition, this paper also develops efficient strategies for the generation of training points to construct the probability-space surrogates, particularly for coupled multi-physics systems. A simplified three-discipline aircraft design problem is used to demonstrate the proposed approach; however, the approach can be used for any multi-physics system such as long-span bridges with aero-structure interactions, or electronic packaging problems with thermal-electrical interactions. The results from the three probability-space surrogates are compared against a Gaussian process algebraic surrogate and fully coupled Monte Carlo analysis. The remainder of the paper is organized as follows. Section II briefly introduces multidisciplinary optimization under uncertainty, literature review on computationally efficient approaches for MDOUU, and probability-space surrogate modeling. Section III discusses the

9

proposed methodology for probability-space surrogates, multidisciplinary optimization, and training point generation for surrogate construction using the one-pass analysis approach. Section IV demonstrates the proposed methodology for reliability-based design optimization (RBDO) and reliability-based robust design optimization (RBRDO) of a three-way coupled aircraft design example. Concluding remarks are provided in Section V. II. Background A. Multidisciplinary design optimization under uncertainty Consider a conceptual two-way coupled multi-physics system, as shown in Fig. 2. Let 𝑿𝑫 represent the set of deterministic design variables, 𝑿𝑼 the uncertain but non-design variables, and 𝑿𝑫,𝑼 the design variables associated with uncertainty. Let 𝒈𝟏 and 𝒈𝟐 represent the outputs of the coupled systems, which are propagated through a third discipline 𝐴3 to obtain the objective functions 𝒇 and constraint functions 𝒄. The coupling variables between the coupled disciplines 𝐴1 and 𝐴2 are represented as 𝒖𝟏𝟐 and 𝒖𝟐𝟏 respectively. An RBDO optimization, where the mean of the objective function is optimized can be written as 𝑀𝑖𝑛 {𝜇[𝑓𝑖 (𝑿𝑫 , 𝑿𝑼 , 𝑿𝑫,𝑼 )]} 𝑖 = 1,2, … 𝑘 𝑠. 𝑡

𝑃𝑟(𝑐𝑗 (𝑿𝑫 , 𝑿𝑼 , 𝑿𝑫,𝑼 , 𝒖𝟏𝟐 , 𝒖𝟐𝟏 ) < 0) > 𝜸𝑗 𝑗 = 1,2, … 𝑚 𝒍𝒃𝒅 ≤ 𝑿𝑫 ≤ 𝒖𝒃𝒅

(1)

Pr(𝑿𝑫,𝑼 ≥ 𝒍𝒃𝑿𝑫,𝑼 ) ≥ 𝒑𝒍𝒃 Pr(𝑿𝑫,𝑼 ≤ 𝒖𝒃𝑿𝑫,𝑼 ) ≥ 𝒑𝒖𝒃 In the above formulation, 𝑓𝑖 ( 𝑖 = 1 𝑡𝑜 𝑘) and 𝑐𝑗 (𝑗 = 1 𝑡𝑜 𝑚) represent the objective and constraint functions respectively. 𝑐𝑗 < 0 represents the safe region and 𝛾𝑗 represents the reliability threshold for the 𝑗𝑡ℎ constraint. Note that for design variables associated with known uncertainty

10

(𝑿𝑫,𝑼 ), we optimize their mean values (𝝁𝑿𝑫,𝑼 ). For deterministic design variables (𝑿𝑫 ), we have strict lower and upper bounds (𝒍𝒃𝒅 and 𝒖𝒃𝒅), and design variables with associated uncertainty (𝑿𝑫,𝑼 ) have probabilistic bounds. 𝒑𝒍𝒃 and 𝒑𝒖𝒃 represent the probability thresholds for 𝑿𝑫,𝑼 that they are in between their lower and upper bounds (𝒍𝒃𝑿𝑫,𝑼 and 𝒖𝒃𝑿𝑫,𝑼 ). Similar formulations for RDO and RBRDO are available in Zaman and Mahadevan [24] and Youn and Xi [25] respectively.

Fig. 2. A conceptual two-way coupled multi-physics system

B. Review of computationally efficient approaches for MDOUU In this section, we review some of the computationally efficient approaches for MDOUU. Mahadevan and Smith [9] proposed a decoupled framework, where the coupling variables between the individual disciplinary analyses are estimated through First-Order Second Moment (FOSM) approximations [12] and later used in reliability analysis. Du et al [10] proposed an efficient RBMDO (Reliability-based Multidisciplinary Design optimization) under uncertainty framework by decoupling it into a deterministic MDO problem and Reliability analysis; these analyses are repeated in cycles, where the results from the reliability analysis in one cycle are used to inform the deterministic MDO analysis in the following cycle. Sankararaman and Mahadevan [11] developed a likelihood-based approach (LAMDA) to obtain the distributions of the coupling

11

variables by using the First-Order Reliability Methods [12] (FORM). Du and Chen [8] proposed a collaborative reliability analysis framework, where the multidisciplinary compatibility conditions are considered as optimization constraints when performing reliability analysis using First-Order Reliability Methods (FORM) [12]. Liang and Mahadevan [26] used a Gaussian copula to connect all the design, coupling variables, objective, constraint functions to perform multidisciplinary and uncertainty analyses simultaneously. The decoupled approaches have primarily used first-order methods (such as FORM and FOSM) to estimate the coupling variables, and considered the bidirectional coupled analysis as uni-directional analysis. As mentioned in Section I, a drawback in the above techniques is that the statistical relationship between the design inputs and coupling variables is disregarded.

C. Probability-space surrogate A probability-space surrogate approximates the joint probability distribution of inputs and outputs, and provides the conditional distribution of the output when the inputs are conditioned at particular values. Some examples of a probability-space surrogate are: BN, MG, GC and GMM. The benefits of a probability-space surrogate are illustrated below using a BN. The other types of probability-space surrogates are discussed later in Section III. A Bayesian network is a directed acyclic graphical model that represents a joint distribution of a set of random variables (both inputs and outputs) as a product of marginal PDFs and conditional PDFs as 𝑛

𝑃(𝑿) = ∏ 𝑃(𝑋𝑖 |Π𝑋𝑖 ) 𝑖=1

12

(2)

where 𝑿 represents a set of random variables, Π𝑋𝑖 represents the set of parent nodes of a random variable 𝑋𝑖 and 𝑃(𝑋𝑖 |Π𝑋𝑖 ) represents the conditional PDF of 𝑋𝑖 , given its parent nodes. When 𝑋𝑖 has no parent nodes (also referred to as a root node), then 𝑃(𝑋𝑖 |Π𝑋𝑖 ) represents its marginal PDF, i.e., 𝑃(𝑋𝑖 ). Fig. 3(a) represents a conceptual system with two inputs 𝑋1 and 𝑋2 , and two outputs 𝑋3 and 𝑋4 . An associated BN is represented as shown in Fig. 3(b). The nodes represent random variables (𝑋1 , 𝑋2 , 𝑋3 and 𝑋4 ), and edges represent the dependence between them through conditional probability distributions.

(a)

(b)

Fig. 3. (a) A system with two inputs and two outputs, and (b) the corresponding BN

A BN can be used for prediction i.e., estimation of outputs 𝑋3 , 𝑋4 for given values of inputs 𝑋1 , 𝑋2 , i.e., 𝑓(𝑋1 , 𝑋2 |𝑋3 = 𝑥3 , 𝑋4 = 𝑥4 ) or model inference 𝑓(𝑋1 , 𝑋2 |𝑋3 = 𝑥3 , 𝑋4 = 𝑥4 ), for inferring the possible values of the inputs that may have resulted in given values of outputs. Here, 𝑓(. ) represents a conditional probability distribution, and 𝑥1 , 𝑥2 , 𝑥3 and 𝑥4 represent the values at which the variables 𝑋1 , 𝑋2 , 𝑋3 and 𝑋4 are conditioned upon. Forward prediction can be evaluated using Monte Carlo sampling while inference analysis can be performed using Markov Chain Monte Carlo (MCMC) sampling [27].

13

III. Proposed Methodology Fig. 4 presents a flow chart of various steps in the proposed methodology for MDOUU using a probability-space surrogate. Each of the steps is discussed in detail later in this section. First, we describe several probability-space surrogates in Section III. A, constructed using the data on various design input, uncertain and coupling variables, and objective and constraint functions. Then, we evaluate the system quantities of interest (QoIs) such as objective and constraint functions, and coupling variables, at multidisciplinary compatibility in Section III. B, using the constructed probability-space surrogate. We detail the one-pass analysis approach in Section III. C for generating training data for constructing the probability-space surrogate, and then the multidisciplinary optimization analysis is discussed in Section III. D.

Fig. 4: Flow chart for MDOUU using a probability-space surrogate 14

A. Construction of a probability-space surrogate As discussed in Section I, there exist two types of options for constructing a probability-space surrogate: (1) A Bayesian network, where the dependences between variables are accurately quantified, but the inference is approximate due to the use of computationally expensive Markov Chain Monte Carlo (MCMC) methods, or (2) an approximate probability-space surrogate, where the dependences are not accurately quantified (in comparison to a Bayesian network), but where the inference is exact (i.e., analytical, not numerical). The computational expense for inference analysis using MCMC methods (first option) multiplies with the number of design iterations, whereas the computational expense in the second option is minimal (compared to the first option).

Fig. 5. Bayesian network of a two-way coupled multi-physics system

Therefore, this paper investigates the efficacy of the second option in performing MDOUU, and compares its performance against a fixed point iteration analysis using Monte Carlo simulations. However, a brief description of the first option is detailed below. Fig. 5 shows the dependence relationships between several variables for the system in Fig. 2. In the system described in Fig. 2, when 𝒖𝟏𝟐 is considered as an input for one-pass analysis, then we have 𝒖𝟏𝟐,𝒊𝒏

15

and 𝒖𝟏𝟐,𝒐𝒖𝒕 , which represent the values of 𝒖𝟏𝟐 before and after the one pass-analysis. 𝒖𝟏𝟐,𝒊𝒏 is input to discipline 𝐴2 from which 𝒖𝟐𝟏 is computed, and then input to discipline 𝐴1 . 𝒖𝟏𝟐,𝒐𝒖𝒕 is the output from 𝐴1 . 𝑫𝒖𝟏𝟐 represents the difference between the ‘in’ and ‘out’ values of 𝒖𝟏𝟐 in a onepass analysis. The outputs 𝒈𝟏 and 𝒈𝟐 from the disciplines 𝐴1 and 𝐴2 , which are dependent on the inputs (design and uncertain - 𝑿𝑫 , 𝑿𝑼 and 𝑿𝑫,𝑼 ) and the coupling variables, are then used to calculate the system-level objective and constraint functions. This dependence information can now be used to obtain the interactions between variables shown in Fig. 4. Note that 𝒈𝟏 and 𝒈𝟐 are not shown as they are intermediate outputs. Note that the values of 𝒖𝟏𝟐,𝒐𝒖𝒕 in one iteration become the values of 𝒖𝟏𝟐,𝒊𝒏 in the following iteration. 𝑫𝒖𝟏𝟐 is not a stochastic variable, but a deterministic variable as it is the difference between the same set of coupling variables in successive iterations (𝑫𝒖𝟏𝟐 = 𝒖𝟏𝟐,𝒐𝒖𝒕 − 𝒖𝟏𝟐,𝒊𝒏 ). As 𝑫𝒖𝟏𝟐 are deterministic nodes, the arcs between 𝒖𝟏𝟐,𝒐𝒖𝒕 and 𝑫𝒖𝟏𝟐 can be reversed such that 𝑫𝒖𝟏𝟐 and 𝒖𝟏𝟐,𝒐𝒖𝒕 become stochastic and deterministic nodes respectively. To achieve multidisciplinary compatibility, we conditionalize the difference variable at zero (𝑫𝒖𝟏𝟐 = 0). The deterministic design variables (𝑿𝑫 ) and design variables with uncertainty (𝑿𝑫,𝑼 ) are treated differently when constructing the probability-space surrogates. For design inputs with aleatory uncertainty (𝑿𝑫,𝑼 ) we typically find its mean value (𝝁𝑿𝑫,𝑼 ), if the aleatory uncertainty in modeled using a Gaussian distribution. For each setting of the deterministic design variable (𝑿𝑫 ), the output distribution will also incorporate the aleatory uncertainty associated with 𝑿𝑫,𝑼 . The proposed probability-space surrogate approach is comprehensive incorporating all the available uncertainty sources, and thus, it can be used for complex design problems with high reliability requirement. The transformed BN surrogate is provided in Fig. 6. Since inference in a Bayesian network using MCMC methods is computationally expensive and approximate, we discuss below three probability-space surrogates 16

with analytical inference: (1) Multivariate Gaussian, (2) Gaussian Copula, and (3) Gaussian Mixture Model. An introduction to the three models and their application to multidisciplinary analysis are discussed below.

Fig. 6. Transformed Bayesian network of a two-way coupled multi-physics system

1. Multivariate Gaussian distribution If 𝑿 = {𝑋1 … 𝑋𝑝 } represent 𝑝 random variables, then their joint probability distribution modeled using a multivariate Gaussian distribution can be represented as 𝑓(𝑿) = Φ(𝑿|𝝁, 𝚺) =

1 𝑝 (2𝜋) 2

1

1 |𝚺|2

exp(− 2 (𝑿 − 𝝁)𝑇 𝚺 −1 (𝑿 − 𝝁))

(3)

where 𝝁 is a 𝑝-dimensional expectation vector and 𝚺 is the 𝑝 × 𝑝 covariance matrix. Given the training data, the mean vector 𝝁 and covariance matrix 𝚺 can be calculated from data as 𝜇𝑖 = 𝐸[𝑋𝑖 ] 𝑇

and 𝛴𝑖,𝑗 = 𝐸 [(𝑋𝑖 − 𝜇𝑖 )(𝑋𝑗 − 𝜇𝑗 ) ] where E[.] is the expectation operator, 𝑋𝑖 represents the training data on the 𝑖 𝑡ℎ variable, and 𝛴𝑖,𝑗 is the 𝑖 𝑡ℎ row and 𝑗𝑡ℎ column element of 𝚺. |𝚺| represents

17

the determinant of the covariance matrix. As discussed in Section I, an MG assumes that every random variable follows a univariate Gaussian distribution. Let 𝑿 represent the set of random variables modeled using the MG, 𝑿𝑜𝑏𝑠 ⊂ 𝑿, and 𝑿𝑢𝑛𝑜𝑏𝑠 represent the observable and unobservable variables. Let 𝒙𝐷 𝑜𝑏𝑠 represent the data available on 𝑿𝑜𝑏𝑠 . In the inference process, the goal is to infer 𝑿𝑢𝑛𝑜𝑏𝑠 using 𝒙𝐷 𝑜𝑏𝑠 on 𝑿𝑜𝑏𝑠 . The posterior distributions of the unobserved variables can be obtained as

𝑓 (𝑿𝑢𝑛𝑜𝑏𝑠 |𝑿𝑜𝑏𝑠 =

=

=

𝒙𝐷 𝑜𝑏𝑠

𝑓(𝑿𝑢𝑛𝑜𝑏𝑠 , 𝒙𝑫 𝑜𝑏𝑠 ) )= 𝑓 (𝑿𝑜𝑏𝑠 = 𝒙𝐷 𝑜𝑏𝑠 )

𝛷(𝑿𝑢𝑛𝑜𝑏𝑠 , 𝑿𝑜𝑏𝑠 = 𝒙𝐷 𝑜𝑏𝑠 |𝝁, 𝜮) 𝐷 ∫ 𝛷(𝑿𝑢𝑛𝑜𝑏𝑠 , 𝑿𝑜𝑏𝑠 = 𝒙𝑜𝑏𝑠 |𝝁, 𝜮) 𝑑𝑿𝑢𝑛𝑜𝑏𝑠

(4)

𝐷 𝛷(𝑿𝑢𝑛𝑜𝑏𝑠 |𝑿𝑜𝑏𝑠 = 𝒙𝐷 𝑜𝑏𝑠 , 𝝁, 𝜮) × 𝛷(𝑿𝑜𝑏𝑠 = 𝒙𝑜𝑏𝑠 |𝝁, 𝜮) ∫ 𝛷(𝑿𝑢𝑛𝑜𝑏𝑠 , 𝑿𝑜𝑏𝑠 = 𝒙𝐷 𝑜𝑏𝑠 |𝝁, 𝜮) 𝑑𝑿𝑢𝑛𝑜𝑏𝑠

where 𝛷(. ) represents a Gaussian distribution, and 𝝁, 𝜮 represent the mean vector and covariance matrix respectively.

2. Gaussian Copula If 𝐹1 (𝑥1 ) … 𝐹𝑝 (𝑥𝑝 ) represent the marginal CDFs of 𝑝 random variables 𝑋1 … 𝑋𝑝 , then a copula function to represent the joint CDF of the random variables can be represented as 𝐶(𝑢1 … 𝑢𝑝 ) = 𝑃[𝐹1 (𝑥1 ) ≤ 𝑢1 , … , 𝐹𝑝 (𝑥𝑝 ) ≤ 𝑢𝑝 ]

(5)

where 𝐶(. ) is the copula function and 𝐶(𝑢1 … 𝑢𝑝 ) represents the joint CDF of marginal CDFs associated with 𝑋1 … 𝑋𝑝 . Several types of copula functions are available to model various dependence relationships between random variables such as Gaussian, Gumbel, Clayton and

18

Independence copulas [28]. Among them, only the Gaussian copula has analytical inference whereas other copulas require sampling-based inference. Following the notation in Section III. A, the dependence between variables (both unobserved and observed) when a Gaussian copula is used can be represented as 𝐶𝐺 (𝒖𝑿𝒐𝒃𝒔 , 𝒖𝑿𝒖𝒏𝒐𝒃𝒔 ) = Φ(Φ −1 (𝒖𝑿𝒐𝒃𝒔 ), Φ−1 (𝒖𝑿𝒖𝒏𝒐𝒃𝒔 ) )

(6)

where 𝐶𝐺 (𝒖𝑿𝒐𝒃𝒔 , 𝒖𝑿𝒖𝒏𝒐𝒃𝒔 ) represents the Gaussian copula, Φ−1 (. ) represents the inverse CDF of a standard Gaussian and Φ(. ) represents the multivariate Gaussian joint CDF. The joint PDF 𝑐𝐺 (𝒖𝑿𝒐𝒃𝒔 , 𝒖𝑿𝒖𝒏𝒐𝒃𝒔 ) can be obtained from the joint CDF as 𝑇

𝛷 −1 (𝒖𝑿𝒐𝒃𝒔 ) 1 𝛷 −1 (𝒖𝑿𝒐𝒃𝒔 ) −𝟏 ( ) (− 𝑐𝐺 (𝒖𝑿𝒐𝒃𝒔 , 𝒖𝑿𝒖𝒏𝒐𝒃𝒔 ) = 𝑒𝑥𝑝 ( ) × 𝑹 − 𝐼 × ( )) 1 2 𝛷 −1 (𝒖𝑿𝒖𝒏𝒐𝒃𝒔 ) 𝛷 −1 (𝒖𝑿𝒖𝒏𝒐𝒃𝒔 ) | 𝑹| 2 1

(7)

where 𝑹 is the correlation matrix of the multivariate standard Gaussian distribution, and |𝑹| represents the determinant of 𝑹. The training data on several variables are transformed into the standard Gaussian space using the inverse CDF technique, which are then used to calculate the correlation matrix. We now discuss inference in a Gaussian copula used to estimate the objective and constraint functions, and coupling variables when the multidisciplinary compatibility is enforced, i.e., 𝑫𝒖𝟏𝟐 = 0. Conditional sampling with the Gaussian copula assumption can be implemented as the joint PDF in Equation (7) can be converted to conditional PDF analytically. All the variables, both unobserved and observed are transformed into an equivalent standard normal space, where the conditional distribution of the unobserved variables is obtained using Equation (7) in the standard normal space. Several samples of the unobserved variables are generated in the standard normal space, and they are transformed back to the original space using the inverse CDF technique. Let the data on observed variables be represented as 𝑿𝒐𝒃𝒔 = 𝒙𝑫 𝒐𝒃𝒔 . The 19

′ equivalent standard normals corresponding to 𝑿𝒐𝒃𝒔 = 𝒙𝑫 𝒐𝒃𝒔 are first calculated as 𝒙𝒐𝒃𝒔 =

𝚽 −𝟏 (𝑭𝑿𝒐𝒃𝒔 (𝒙𝑫 𝒐𝒃𝒔 )) , where 𝑭𝑿𝒐𝒃𝒔 (. ) is the CDF of 𝒙𝒐𝒃𝒔 . Let the associated standard normal values of 𝑿𝒖𝒏𝒐𝒃𝒔 be represented as 𝑿′𝒖𝒏𝒐𝒃𝒔 . The conditional distribution 𝑓 (𝑿′𝑢𝑛𝑜𝑏𝑠 |𝑿𝑜𝑏𝑠 = 𝒙′𝑜𝑏𝑠 ) is computed using Equation (4), which is transformed back to the original space to obtain 𝑓 (𝑿𝑢𝑛𝑜𝑏𝑠 |𝑿𝑜𝑏𝑠 = 𝒙𝐷 𝑜𝑏𝑠 ) using inverse CDF technique. 3. Gaussian Mixture Model As discussed in Section I, a Gaussian mixture model represents the joint probability distribution of 𝑝 random variables {𝑋1 , 𝑋2 … 𝑋𝑝 }

as a weighted sum of several Gaussian

distributions. The joint PDF 𝑓(𝑿) can be represented using an 𝑁-component GMM as 𝑁

𝑓 (𝑿) = ∑ 𝑤𝑖 × 𝛷𝑖 (𝑿|𝝁𝒊 , 𝚺𝒊 )

(8)

𝑖=1

where 𝛷𝑖 (𝑿|𝝁𝑖 , 𝚺𝑖 ) and 𝑤𝑖 represent the 𝑖 𝑡ℎ Gaussian Mixture component and its corresponding weight respectively. The sum of the weights of all Gaussian components is equal to 1, i.e., ∑𝑁 𝑖=1 𝑤𝑖 = 1. The first step in approximating a joint PDF through an N-component GMM is the selection of the value of 𝑁. The optimal number of components can be estimated by maximizing a model selection score such as the Akaike Information Criterion (AIC), Bayesian Information Criterion (BIC) or Mutual Information (MI) [29]. For illustration, this paper uses the BIC score which is the defined as 𝐵𝐼𝐶 = 𝑘 × ln(𝑀) − 2 ln(𝐿)

(9)

where 𝐿, 𝑘 and 𝑀 refer to the likelihood of observing the data given a GMM, the number of parameters to be estimated in the GMM and the amount of available data respectively. The

20

parameters include the weights (𝑤𝑖 ), the expectation and variance matrices of all N-components in the GMM. Parameter estimation in a GMM is typically carried out through the ExpectationMaximization (EM) [30] algorithm. The EM algorithm performs the parameter estimation in two steps: the E step and the M step. The E step estimates the weights of components using Maximum Likelihood Estimation (MLE) conditioned on the mean vectors and covariance matrices of the Gaussian Mixture components, and the M step estimates the mean vector and covariance matrices of the components conditioned on the weights. The E and M steps are repeated several times until the values of weights, mean vectors and covariance matrices of several components converge. The posterior distributions of the unobserved variables can be obtained as

𝑓(𝑿𝑢𝑛𝑜𝑏𝑠 |𝑿𝑜𝑏𝑠 = 𝒙𝐷 𝑜𝑏𝑠 ) =

𝑓(𝑿𝑢𝑛𝑜𝑏𝑠 , 𝒙𝑫 𝑜𝑏𝑠 ) 𝐷 𝑓(𝑿𝑜𝑏𝑠 = 𝒙𝑜𝑏𝑠 )

𝐷 ∑𝑁 𝑖=1 𝑤𝑖 × 𝛷𝑖 (𝑿𝑢𝑛𝑜𝑏𝑠 , 𝑿𝑜𝑏𝑠 = 𝒙𝑜𝑏𝑠 |𝝁𝒊 , 𝚺𝒊 ) = 𝐷 ∫ ∑𝑁 𝑖=1 𝑤𝑖 × 𝛷𝑖 (𝑿𝑢𝑛𝑜𝑏𝑠 , 𝑿𝑜𝑏𝑠 = 𝒙𝑜𝑏𝑠 |𝝁𝒊 , 𝚺𝒊 ) 𝑑𝑿𝑢𝑛𝑜𝑏𝑠

=

(10)

𝐷 𝐷 ∑𝑁 𝑖=1 𝑤𝑖 × 𝛷𝑖 (𝑿𝑢𝑛𝑜𝑏𝑠 |𝑿𝑜𝑏𝑠 = 𝒙𝑜𝑏𝑠 , 𝝁𝒊 , 𝚺𝒊 ) × 𝛷𝑖 (𝑿𝑜𝑏𝑠 = 𝒙𝑜𝑏𝑠 |𝝁𝒊 , 𝚺𝒊 ) 𝐷 ∫ ∑𝑁 𝑖=1 𝑤𝑖 × 𝛷𝑖 (𝑿𝑢𝑛𝑜𝑏𝑠 , 𝑿𝑜𝑏𝑠 = 𝒙𝑜𝑏𝑠 |𝝁𝒊 , 𝚺𝒊 ) 𝑑𝑿𝑢𝑛𝑜𝑏𝑠

𝑁 𝐷 where 𝑤𝑖 , 𝛷𝑖 (𝑿𝑜𝑏𝑠 = 𝒙𝐷 𝑜𝑏𝑠 ) and ∫ ∑𝑖=1 𝑤𝑖 × 𝛷𝑖 (𝑿𝑢𝑛𝑜𝑏𝑠 , 𝑿𝑜𝑏𝑠 = 𝒙𝑜𝑏𝑠 |𝝁𝒊 , 𝚺𝒊 ) 𝑑𝑿𝑢𝑛𝑜𝑏𝑠 are

constants, and therefore Equation (10) can be simplified as 𝑁 𝐷 𝑓 (𝑿𝑢𝑛𝑜𝑏𝑠 |𝑿𝑜𝑏𝑠 = 𝒙𝐷 𝑜𝑏𝑠 ) = ∑ 𝛾𝑖 × 𝛷𝑖 (𝑿𝑢𝑛𝑜𝑏𝑠 | 𝑿𝑜𝑏𝑠 = 𝒙𝑜𝑏𝑠 )

(11)

𝑖=1

where 𝛾𝑖 is a new constant combining all the three constants described above. 𝑁𝑖 (𝑿𝑢𝑛𝑜𝑏𝑠 | 𝑿𝑜𝑏𝑠 = 𝑡ℎ 𝒙𝐷 Gaussian mixture component, which is also a 𝑜𝑏𝑠 ) represents the posterior distribution of the 𝑖

multivariate

Gaussian

distribution.

Therefore,

the

overall

posterior

distribution,

𝑓 (𝑿𝑢𝑛𝑜𝑏𝑠 |𝑿𝑜𝑏𝑠 = 𝒙𝐷 𝑜𝑏𝑠 ), represents another GMM. Thus, a GMM allows for analytical inference.

21

4. Model Characteristics Here, we discuss the underlying assumptions in each of these distribution models. The Multivariate Gaussian assumes that each individual random variable follows a Gaussian distribution, and can be used to model a uni-modal joint distribution. A less restrictive model compared to Multivariate Gaussian is the Gaussian Copula, which does not assume that the individual random variables follow Gaussian distributions. The Gaussian copula can include nonGaussian individual distributions, but it still only models a uni-modal joint distribution. A Gaussian Mixture Model is the least restrictive of all the three models, as it does not assume the individual variables to follow Gaussian distributions and can be used to model multi-modal joint distributions. By changing the number of Gaussian components in the GMM, any order of modality can be represented; however, it should be noted that a multivariate Gaussian distribution is used in each of the individual GMM components. We need only one probability-space surrogate to perform the optimization. We can use model selection techniques such as the BIC score (Equation 9) to identify the best probability-space surrogate among MG, GC and GMM, and then use it for optimization analysis.

B. Evaluation of system QoIs at multidisciplinary compatibility In order to obtain the values of QoIs (objective and constraint functions, coupling variables) at convergence between individual disciplines, the multidisciplinary compatibility constraints are imposed by conditioning the difference variables 𝑫𝒖 at zero, and the design variables are conditioned at some chosen design values. Use the above conditionalization; conditional sampling is carried out to obtain the conditional PDFs of the objective and constraint functions. Let 𝑥𝐷∗ and

22

∗ 𝜇𝐷,𝑈 refer to the design values of 𝑋𝐷 and 𝜇𝑋𝐷,𝑈 at which the objective (𝒇), constraint functions (𝒄)

and the coupling variables (𝒖𝟏𝟐 and 𝒖𝟐𝟏 ) need to be evaluated. The Gaussian copula model between the design and coupling variables, objective and constraint functions is given in Equation (12), where 𝛴 refers to the covariance matrix. Given observed data, i.e., 𝑿𝑫 = 𝒙∗𝑫 , 𝝁𝑫,𝑼 = 𝝁∗𝑫,𝑼 and 𝑫𝒖 = 𝟎 (multidisciplinary compatibility), the conditional joint distribution of the objectives, constraints and coupling variables can be evaluated following the procedure laid down in Section III.A.2. 𝑐𝐺 (𝒖𝑿𝑫 , 𝒖𝝁𝑿𝑫,𝑼 , 𝒖𝑿𝑼 , 𝒖𝒖𝟏𝟐 , 𝒖𝒖𝟐𝟏 , 𝒖𝑫𝒖 , 𝒖𝒇 , 𝒖𝒄 ) Φ−1 (𝒖𝑿𝑫 )

1

1 = exp − 2 √det(𝛴 )

(

If 𝑿 = {𝑿𝑫 , 𝒖𝝁𝑿

𝑫,𝑼

𝑇

Φ−1 (𝒖𝑿𝑫 )

Φ−1 (𝒖𝝁𝑿𝑫,𝑼 )

Φ −1 (𝒖𝝁𝑿𝑫,𝑼 )

Φ−1 (𝒖𝑿𝑼 )

Φ−1 (𝒖𝑿𝑼 )

Φ −1 ( 𝒖𝒖𝟏𝟐 ) Φ−1 (𝒖𝒖𝟐𝟏 )

× (𝛴 −1 − 𝐼 ) ×

(12)

Φ−1 ( 𝒖𝒖𝟏𝟐 ) Φ−1 (𝒖𝒖𝟐𝟏 )

Φ−1 (𝒖𝑫𝒖 )

Φ −1 (𝒖𝑫𝒖 )

Φ−1 (𝒖𝒇 ) ( Φ−1 (𝒖𝒄 ) )

Φ−1 (𝒖𝒇 ) ( Φ−1 (𝒖𝒄 ) )

)

, 𝑿𝑼 , 𝒖𝟏𝟐,𝒊𝒏 , 𝒖𝟐𝟏 , 𝑫𝒖𝟏𝟐 , 𝒇, 𝒄}, then their joint PDF 𝑓(𝑿) when represented

as a GMM is given in Equation (8) and the inference can be carried out using Equation (10). Here, the observed variables (𝑿𝑜𝑏𝑠 ) include the design variables and the multidisciplinary compatibility, and unobserved variables (𝑿𝑢𝑛𝑜𝑏𝑠 ) include the design variables with uncertainty, the coupling variables, the objective and constraint functions. Therefore, 𝑿𝒐𝒃𝒔 = {𝑿𝑫 , 𝝁𝑫,𝑼 , 𝑫𝒖 } and 𝒙𝐷 𝑜𝑏𝑠 = {𝒙∗𝑫 , 𝝁∗𝑫,𝑼 , 𝟎} and 𝑿𝑢𝑛𝑜𝑏𝑠 = {𝑿𝑼 , 𝒖𝟏𝟐,𝒊𝒏 , 𝒖𝟐𝟏 , 𝒇, 𝒄}. From the joint posterior distribution, the marginal posterior distributions of the objectives and constraints can be obtained for further 23

analysis such as reliability calculations for the constraints. The use of a probability-space surrogate enables us to perform the uncertainty analysis and multidisciplinary convergence analysis simultaneously. This offers a tremendous computational advantage over an algebraic surrogate, since the algebraic surrogate needs to be trained with samples of fully converged multi-disciplinary analyses, whereas the probability-space surrogate is trained with only one-pass samples of the multi-disciplinary analysis. (In the probability-space surrogate, the multi-disciplinary compatibility is cleverly enforced through Bayesian updating, for which analytical solutions are available for the proposed probability-space surrogates). Thus, the nested three-loop design optimization (shown in Fig. 1) that is required when performing fully coupled analysis using physics models or algebraic surrogates can now be collapsed into a double loop analysis as shown in Fig. 7. This advantage becomes important when the disciplinary physics models are expensive.

Optimization analysis Uncertainty analysis and Multidisciplinary analysis

Fig. 7. Nested three-loop analysis (in Fig. 1) simplified into a two-loop analysis using a probability-space surrogate

C. Generation of training points using one-pass analysis Section III.A discussed the construction of a probability-space surrogate using available data on several design input, uncertain and coupling variables, and objective and constraint functions. Here, we detail the one-pass analysis approach for the generation of those training points in multiphysics systems with more than two coupled disciplines. Performing one-pass analysis requires

24

the identification of a feasible one-pass path through the individual disciplinary models; this path determines the sequence in which the individual disciplinary analyses need to be carried out to generate the training points. We detail an algorithm to identify a one-pass analysis path in a generic multi-physics system with bi-directional interactions and a combination of uni-directional and bidirectional interactions. We demonstrate the algorithm for (1) a three-way coupled system with bidirectional interactions, and (2) a three-way coupled system with both bi-directional and unidirectional interactions. 1. Algorithm to generate a one-pass analysis path Let 𝑛 represent the number of disciplinary systems; this includes the analysis to compute the system-level objectives and constraints. For example, 𝑛 = 3 for the system in Fig. 2. Let us represent the interactions between 𝑛 disciplinary analysis using an 𝑛 × 𝑛 interaction matrix with elements 𝐼(𝑖, 𝑗) where 𝐼(𝑖, 𝑗) = 1 if there a dependence of 𝑖 𝑡ℎ discipline on the 𝑗𝑡ℎ discipline, i.e., there exists an arrow from 𝐴𝑖 to 𝐴𝑗 . Since dependence to oneself are not valid, the diagonal values are zero, i.e., 𝐼(𝑖, 𝑖) = 0. Note that in the presence of two-way interactions between two disciplines 𝑖 and 𝑗, then 𝐼(𝑖, 𝑗) = 1 and 𝐼(𝑗, 𝑖) = 1. Given the interaction matrix, the steps below need to be followed for generating a one-pass analysis path. (1) Create a list of ordered pairs (𝑖, 𝑗) of disciplines with two way interactions between them. Let 𝑘 represent the number of such ordered pairs. When disciplines are numbered, then it is not required that 𝑖 < 𝑗 (2) For index ‘𝑚’ going from 1 to 𝑘 – 1, 𝐼(𝑖𝑚 , 𝑗𝑚 ) = 0. (3) With regard to 𝑘 𝑡ℎ pair, either 𝐼(𝑖𝑘 , 𝑗𝑘 ) = 0 or 𝐼(𝑗𝑘 , 𝑖𝑘 ) = 0; the substitution that results in a column of all zeros is the chosen.

25

The resultant interaction matrix provides a one-pass analysis path. Two lists are created, namely, Analyzed() and NotAnalyzed(). These lists specify the order in which various disciplines are to be analyzed. Analyzed() and NotAnalyzed() represent the lists of disciplines that are analyzed and to be analyzed. Initially, Analyzed() is an empty list while NotAnalyzed() represents all the disciplines. First, we start with the one-pass analysis with analyzing the discipline ℎ for which all the associated column elements are zero, i.e., 𝐼(: , ℎ) = 0. If there are many such disciplines, then all of them are analyzed sequentially in any random order. All the analyzed disciplines are then added to Analyzed() and removed from NotAnalyzed(). Next, we analyze the disciplines which only have interactions starting from the disciplines in Analyzed(). Such disciplines are then removed from NotAnalyzed() and added to Analyzed(). This process is repeated until all the disciplines are analyzed and NotAnalyzed() becomes an empty list. We demonstrate this algorithm for two-way and three-way coupled systems below.

Fig. 8. Conceptual three-way coupled multi-physics system

2. Conceptual three-way coupled system

26

Consider a conceptual coupled three-disciplined system as shown in Fig. 8, where 𝐴1 , 𝐴2 and 𝐴3 represent the three coupled disciplines, and 𝐴4 represents a fourth discipline that does not have coupling with 𝐴1 , 𝐴2 or 𝐴3 . Let 𝑢𝑖𝑗 represent the coupling variables between any two disciplines, which are output from the 𝑖𝑡ℎ discipline and input to the 𝑗𝑡ℎ discipline. 𝑔1 , 𝑔2 and 𝑔3 represent the outputs from 𝐴1 , 𝐴2 and 𝐴3 respectively that are input to 𝐴4 , the outputs of which are the values of quantities of interest. The interaction matrix for the system in Fig. 8 is given in Fig. 9a. The ordered pairs that correspond to two-way interactions between disciplines are (1, 2), (1, 3) and (2, 3). Here, k = 3. Following step 2, we choose 𝐼 (1, 2) = 0 and 𝐼 (1, 3) = 0. For the 3rd pair, either 𝐼 (2, 3) = 0 or 𝐼(3, 2) = 0 will result in a column with all zeros. Let us choose 𝐼 (2, 3) = 0. The resultant interaction matrix is given in Fig. 9b. Let us create two lists – Analyzed() and NotAnalyzed(). Initially, NotAnalyzed() = [𝐴1 , 𝐴2 , 𝐴3 , 𝐴4 ] and Analyzed() = []. We start the onepass analysis with 𝐴3 , since 𝐼 (: ,3) = 0. Therefore, Analyzed() = [𝐴3 ] and NotAnalyzed() = [𝐴1 , 𝐴2 , 𝐴4 ]. The next analyzed discipline is 𝐴2 as it has an interaction from 𝐴3 , which is present in Analyzed(). Now, Analyzed() = [𝐴3 , 𝐴2 ] and NotAnalyzed() = [𝐴1 , 𝐴4 ]. Following this approach, the next analyzed discipline is 𝐴1 and then 𝐴4 . Therefore, Analyzed() = [𝐴3 , 𝐴2 , 𝐴1 , 𝐴4 ], which provides the order of analysis of various disciplines. The resultant interaction matrix for one-pass analysis path is given in Fig. 9b and the system is shown in Fig. 10.

(a) (b) Fig. 9. Interaction matrix in (a) original system, and (b) one-pass analysis path

27

Fig. 10. A one-pass analysis path for a three-discipline system

Fig. 11. Three-discipline coupled system with one-directional coupling between two disciplines

3. Three-way coupled system with one-directional interaction between two disciplines In some problems, there might be one-directional interaction between some of the individual disciplines as shown in Fig. 11. In Fig. 11, one-directional coupling exists between 𝐴1 and 𝐴2 whereas bi-directional couplings exist between 𝐴1 and 𝐴3 , and 𝐴2 and 𝐴3 . For the three-discipline system with both one-way and two-way interactions, the interaction matrix is given in Fig. 12a. The ordered pairs with two-way interactions are (1, 3) and (2, 3). Following step 2, we have

28

𝐼 (1,3) = 0. For the last pair, either 𝐼 (2, 3) = 0 or 𝐼 (3, 2) = 0 will result in a column with all zeros. Let us choose 𝐼 (3, 2) = 0. The original and resultant interaction matrices is given in Fig. 12a and 12b respectively. In this case, the order of analysis of various disciplines is [𝐴2 , 𝐴3 , 𝐴1 , 𝐴4 ]. The system with the one-pass analysis path is shown in Fig. 13.

(a) (b) Fig. 12. Interaction matrix in (a) original system, and (b) one-pass analysis path

Fig. 13. One-pass analysis path of a three-discipline coupled system with one-directional coupling between two disciplines, resulting in one-directional coupling between all disciplines

D. Optimization under uncertainty In this section, we discuss an optimization algorithm that can be used for design optimization analysis and a procedure to handle multi-objective optimization formulations. 1. Optimization algorithm

29

Optimization algorithms are broadly classified into two categories, based on the usage of gradients, as - (1) gradient-based, and (2) gradient-free [31]. Gradient-based algorithms use the gradients of objective functions while gradient-free algorithms do not use gradients but only the function evaluations to reach the optimum solution. Examples of gradient-based algorithms include Newton-Raphson method, Steepest-descent method and Conjugate-gradient method [32]. Some examples of gradient-free algorithms include DIviding RECTangles (DIRECT) [33,34], genetic algorithm, particle swarm optimization and simulated annealing [32]. This paper uses the DIRECT algorithm following our previous work [35]. For the sake of completeness, a brief description of the algorithm is given below. The idea behind the DIRECT algorithm is the systematic adaptive meshing of the design space until a global optimum is reached. The original design space is transformed into a unit hypercube. First, the center of the hypercube is evaluated and then the hypercube is trisected along all dimensions resulting in 2𝐻 + 1 hyper-rectangles; 𝐻 represents the dimension of the design input space From the available hyper-rectangles, a set of potentially optimal hyper-rectangles are selected based on two factors: (1) function value, and (2) distance of a hyper-rectangle from the center of the hypercube. The selected hyper-rectangles are further trisected, and from them, a set of hyper-rectangles are again selected. This process of optimal hyper-rectangle selection and adaptive meshing is carried out until the global optimum is reached. 2. Converting a multi-objective formulation to several single objective formulations In this paper, we adopt a widely-used technique whereby a multi-objective formulation is converted to several single optimization problems through a weighted-sum of the individual objective functions. The new single objective function is defined as

30

𝐹 = ∑ 𝛽𝑖

𝐻[𝑓𝑖 ] , ∑ 𝛽𝑖 = 1 𝑓𝑖,𝑡ℎ

(13)

where 𝑓𝑖 , 𝑖 = 1,2 … 𝑘 refer to 𝑘 competing objective functions. 𝑓𝑖,𝑡ℎ and 𝛽𝑖 correspond to the normalization values and weights associated with each of the objective functions. 𝐻[. ] is an operator to compute either the expectation, standard deviation or variance of an objective function. Different values of weight coefficients (𝛽𝑖 ) result in different single optimization problems, and for each of these problems, multidisciplinary optimization is carried out and an optimum design point is obtained; thereby, resulting in several optimum design points, one for each value of 𝛽𝑖 . The optimal design points can be visualized through a Pareto surface, which can be used for further decision-making such as obtaining the best design configuration among all the optimal design configurations considering trade-offs between individual objective functions. The number of weight coefficient combinations depends on the desired accuracy of the Pareto surface.

Fig. 14. Individual systems and the coupling variables in a multidisciplinary aircraft system

31

IV. Illustration Example: Aircraft Design Problem Description: The three disciplinary analyses (Weights, Engine and Aerodynamics) and their coupling are shown in Fig. 14. The system in Fig. 14 is a three-discipline system with bidirectional couplings between Weights analysis and Aerodynamics, and Engine analysis and Aerodynamics, but only a uni-directional interaction between Engine and Weight analyses. The model for the Weights analysis was obtained from empirical data, while the models for Engine and Aerodynamics analyses are obtained from a combination of physics principles and empirical data. 𝐷𝑉 and 𝑈𝑉 represent the vectors of design and uncertain variables. In this example, there are six design and six uncertain variables (mentioned in Table 1 and Table 2). The bounds of design variables are given in Table 1 and the distribution parameters of uncertain variables are given in Table 2. The performance of the aircraft is evaluated using Take-off Field Length (𝑇𝑂𝐹𝐿), Range, Approach speed (𝑉𝑎𝑝𝑝 ), Nacelle Max Diameter (𝐷𝑁𝐴𝐶) and Block Fuel (𝐵𝐹). Using the above system description, we consider two different optimization formulations (RBDO and RBRDO) and compare the results from probability-space surrogates against an algebraic surrogate and fully coupled analysis using physics-based models and Monte Carlo simulations. Table 1. Design variables and their bounds Parameter Wing Aspect Ratio (𝐴𝑅) Wing area (𝑆𝑊) (𝑓𝑡 2 ) Wing Sweep (𝑆𝑊𝐸𝐸𝑃) Gross Weight (𝑇𝑂𝐺𝑊) (𝑙𝑏𝑠) Fan Pressure Ratio (𝐹𝑃𝑅𝐷) Overall Pressure Ratio (𝑂𝑃𝑅𝐷)

Minimum 9 1118.4 20 150000 1.35 45

32

Maximum 12 1677.6 30 188000 1.55 55

Table 2. Parameters associated with uncertainty and their distributions Parameter

Distribution

Minimum

Mode

Maximum

Wing weight factor (𝐹𝑅𝑊𝐼)

Triangular

0.9

1.05

1.1

Fuselage weight factor (𝐹𝑅𝐹𝑈)

Triangular

0.9

1.05

1.1

Induced drag factor (𝐹𝐶𝐷𝐼)

Triangular

0.9

1.02

1.1

Profile drag factor (𝐹𝐶𝐷𝑂)

Triangular

0.8955

1.0149

1.0945

Fan design efficiency curve factor (𝜂𝐹𝐴𝑁 )

Triangular

-0.052

-0.027

0.048

High pressure compressor design efficiency curve factor (𝜂𝐻𝑃𝐶 )

Triangular

-0.05234

-0.02734

0.04766

A. Reliability-based design optimization (RBDO) An RBDO problem is formulated to minimize the Block Fuel consumption and maximize the Range while satisfying a reliability constraint defined with respect to 𝐷𝑁𝐴𝐶, 𝑉𝑎𝑝𝑝 and 𝑇𝑂𝐹𝐿 as 𝑀𝑖𝑛𝑖𝑚𝑖𝑧𝑒 𝐵𝐹 𝑎𝑛𝑑 𝑀𝑎𝑥𝑖𝑚𝑖𝑧𝑒 𝑅𝑎𝑛𝑔𝑒 (14) 𝑠𝑢𝑐ℎ 𝑡ℎ𝑎𝑡 Pr(𝑉𝑎𝑝𝑝 < 𝑉𝑎𝑝𝑝,𝑡ℎ ∩ 𝐷𝑁𝐴𝐶 < 𝐷𝑁𝐴𝐶𝑡ℎ ∩ 𝑇𝑂𝐹𝐿 < 𝑇𝑂𝐹𝐿𝑡ℎ ) > 𝛾 where 𝑉𝑎𝑝𝑝,𝑡ℎ = 145 𝑘𝑛𝑜𝑡𝑠 (𝑘𝑡𝑠), 𝐷𝑁𝐴𝐶𝑡ℎ = 7.9 𝑓𝑡 and 𝑇𝑂𝐹𝐿𝑡ℎ = 8350 𝑙𝑏𝑠 are the threshold values of 𝑉𝑎𝑝𝑝 , 𝐷𝑁𝐴𝐶 and 𝑇𝑂𝐹𝐿 respectively. 𝛾 represents the reliability threshold and is equal to 0.9. It should be noted that 𝐷𝑁𝐴𝐶 is a performance indicator and also happens to be a coupling variable between the Aerodynamics and the Engine disciplines. Since it is a coupling variable, it is not shown as an output from the Performance discipline, whereas the other performance indicators such as the objective function variables (Block Fuel and Range), and other constraint function variables (in addition to 𝐷𝑁𝐴𝐶) such as 𝑉𝑎𝑝𝑝 and 𝑇𝑂𝐹𝐿 are the outputs from the Performance discipline shown in Fig. 14. A reliability threshold of 0.9 implies a failure probability of 0.1. Since the number of Monte Carlo samples required for failure probability estimation (for a 33

given tolerance) is inversely related to the square root of the failure probability [12], we chose a failure probability of 0.1 in order to reduce the effort in computing our Monte Carlo results using fully coupled analysis. The value is chosen only for illustration of the proposed methods and is not meant to be realistic. Realistic applications will require a much smaller failure probability. The multiple objectives are combined into a single objective using the weighted-sum approach as

𝑀𝑖𝑛𝑖𝑚𝑖𝑧𝑒 𝑂 = 𝛽 ×

𝐸[𝐵𝐹] 𝐸[𝑅𝑎𝑛𝑔𝑒] − (1 − 𝛽 ) × 𝐵𝐹𝑏𝑙 𝑅𝑎𝑛𝑔𝑒𝑏𝑙

(15)

where 𝑂 represents the new single objective function, 𝐵𝐹𝑏𝑙 = 36734 and 𝑅𝑎𝑛𝑔𝑒𝑏𝑙 = 2952 represent the baseline values of BF and Range respectively. In Equation (15), 𝐸[. ] represents the expectation function; 𝛽 and (1 − 𝛽) represent the weights of BF and Range used to obtain a weighted sum of the objective functions. The range of 𝛽, between 0 and 1, is uniformly discretized in intervals of 0.05 totaling to 21 values of 𝛽 including 0 and 1.

(a) (b) Fig. 15. Interaction matrix between various disciplines in the multidisciplinary aircraft design example: (a) original and (b) one-pass analysis path

One-pass analysis path to generate training points: The next step is the organization of onepass analysis of the multidisciplinary system for the generation of training points to build the probability-space surrogate. Using the algorithm in Section III.C, the interaction matrix for the original system is given in Fig. 15a, and the list of ordered pairs of bi-directional interactions are

34

(W,A) and (E,A). In this example, 𝐶𝑟𝑠𝐷𝑟𝑎𝑔 is a common interaction variable between W and A, and between E and A. When 𝐶𝑟𝑠𝐷𝑟𝑎𝑔 is considered as input, the link from A to E, and from A to W can be removed. The resultant interaction matrix and the corresponding one-pass analysis path are shown in Fig. 15b and 16 respectively. For this one-pass analysis path, only 𝐶𝑟𝑠𝐷𝑟𝑎𝑔 is considered as input, whose range is given as [6800, 10500].

Fig. 16. A one-pass analysis path for the multidisciplinary aircraft system

Table 3. Variables in the probability-space surrogate Type of variable 𝐷𝑉 𝑈𝑉 ‘In’ coupling variables

Quantity 6 6 1

Variables 𝐴𝑅, 𝑆𝑊, 𝑆𝑊𝐸𝐸𝑃, 𝑇𝑂𝐺𝑊, 𝐹𝑃𝑅𝐷, 𝑂𝑃𝑅𝐷 𝐹𝑅𝑊𝐼, 𝐹𝑅𝐹𝑈, 𝐹𝐶𝐷𝐼, 𝐹𝐶𝐷𝑂, 𝜂𝐹𝐴𝑁 , 𝜂𝐻𝑃𝐶 𝐶𝑟𝑠𝐷𝑟𝑎𝑔,𝑖𝑛 ,

‘Difference’ coupling variables Coupling variables not considered as inputs Intermediate variables Objective functions Constraint functions

1 4

𝐶𝑟𝑠𝐷𝑟𝑎𝑔,𝑑𝑖𝑓𝑓 𝐶𝑟𝑠𝑊𝑡 , 𝑊𝑒𝑛𝑔 , 𝐷𝑁𝐴𝐶, 𝑋𝑁𝐴𝐶

𝐶𝑟𝑠𝐿𝑜𝐷 , 𝑂𝐸𝑊, 𝐶𝑟𝑠𝑆𝐹𝐶 𝐵𝐹, Range 𝑉𝑎𝑝𝑝 , 𝑇𝑂𝐹𝐿 ** DNAC is both a coupling variable and a performance indicator in this case study 3 2 2

35

For the one-pass analysis, several values of design variables, uncertain variables, and ‘in’ coupling variables are generated to obtain the corresponding ‘out’ coupling variables, objectives and constraint values. We used second iteration analysis data since the ‘in’ and ‘out’ values of the coupling variables in the first one-pass analysis were found to be not stable; this instability was also previously observed in [14]. For the second pass analysis, the values of the ‘out’ coupling variables in the first pass are used as values of the ‘in’ coupling variables in the second pass, and the values of the ‘out’ coupling variables, objective and constraint functions are calculated. After performing the first iteration of the one-pass analysis, the difference values are computed between the ‘in’ and ‘out’ of the coupling variables that are considered as inputs. The ideal outcome is a sufficient number of positive and negative difference values; this is important when performing Bayesian inference by conditionalizing the difference variables at zero value. If most (or all) of the difference values are either positive or negative, then the value of zero will be at the tail end or even outside the fitted probability distribution of the difference values. Performing inference in such cases may provide inaccurate and biased results. Therefore, it is necessary to check the distribution of the difference variables at the end of the one-pass analysis. This is what we mean by stability analysis here. If the first iteration of multi-disciplinary analysis does not provide a balanced proportion of positive and negative values, then we use the second iteration results of the multi-disciplinary analysis to construct the probability-space surrogate. After generating the training points, we fit the three probability-space surrogates (Multivariate Gaussian, Gaussian Copula and Gaussian Mixture Models) discussed in Section III and compare their results with the results from fully coupled multidisciplinary analysis and Monte Carlo simulations. In total, there are 25 variables, and their details are given in Table 3. The intermediate

36

variables refer to variables that are not involved in coupling between individual disciplines but are used in the calculation of objective and constraint functions. To identify the optimal number of training points, we compared the performance of the probability-space surrogate against the fully coupled analysis with varying training sample size. The prediction accuracy of the probability-space surrogates with respect to Monte Carlo results with varying number of training points was calculated, and the results are shown in Fig. 17. The set of training points chosen for the study are 100, 250, 500, 1000, 1500, 2000, 2500, 3000, 3500, 4000, 4500, 5000, 5500, 6000, 6500, and 7000. For a given number of training points, the best probability-space surrogate is considered based on the BIC score, and its performance is compared against Monte Carlo results (MG was found to be the best surrogate for 100 and 250 training points while GMM is found to be the best surrogate at all the other number of training points, from 500 to 7000). 25 validation inputs obtained using Latin Hypercube sampling were chosen. At each validation input, the prediction of the mean and standard deviation of the Block Fuel was considered as it was later used in the objective functions of both RBDO and RBRDO formulations. The predicted mean and standard deviation values using the probability-space surrogate were then compared against those computed from the fully coupled analysis. The performance of several surrogates built with varied sample sizes are compared at the same 25 validation points.

(a)

(b)

37

Fig. 17. Average of the magnitude percentage difference in the prediction of (a) mean, and (b) standard deviation of Block Fuel using the probability-space surrogate when compared against the Monte Carlo result at 25 validation points Figs. 17a and 17b show the performance of the probability-space surrogate in the prediction of mean and standard deviation. The average of the magnitude percentage difference in the prediction from the probability-space surrogate when compared against the fully coupled analysis was used as the performance measure, which is calculated as

𝑛

|𝑣𝑠,𝑖 − 𝑣𝑓𝑐,𝑖 | 1 𝜖 (%) = ∑ × 100 𝑛 𝑣𝑓𝑐,𝑖

(16)

𝑖=1

In Eq. (16), 𝑛 represents the number of validation points (here, 𝑛 = 25), 𝑣𝑠,𝑖 and 𝑣𝑓𝑐,𝑖 represents the estimates from the probability-space surrogate and fully coupled analysis respectively at the 𝑖 𝑡ℎ validation point, and |. | is the absolute value operator. It can be observed from Fig. 17a that the average magnitude of the percentage difference in the mean prediction is less than 1% for all the probability-space surrogates built with different numbers of training points. However, there is a higher percentage difference in the standard deviation prediction (Fig. 17b), about 10% or lower for 1000 or more training points.

Fig. 18. Average of the magnitude percentage difference in the joint reliability prediction used in the RBRDO analysis using the probability-space surrogate when compared against the Monte Carlo result at 25 validation points

38

In addition to mean and standard deviation, we also investigated the performance of the probability-space surrogate in predicting the reliability over the performance variables later used in the optimization analysis. We illustrated the joint reliability prediction used in the RBRDO analysis (Section IV. B), i.e., Pr(𝑉𝑎𝑝𝑝 < 145 ∩ 𝐷𝑁𝐴𝐶 < 7.9 ∩ 𝑇𝑂𝐹𝐿 < 8350 ∩ 𝑅𝑎𝑛𝑔𝑒 > 2952) as it is a more stricter constraint when compared to the reliability constraint in RBDO analysis. Fig. 18, which provides the variation in the average of the magnitude percentage difference (computed using Eq. (16)) in the reliability prediction of the probability-space surrogate when compared against the fully coupled Monte Carlo analysis. We considered 25 points that satisfied the RBRDO reliability constraint, i.e., points at which the Pr(𝑉𝑎𝑝𝑝 < 145 ∩ 𝐷𝑁𝐴𝐶 < 7.9 ∩ 𝑇𝑂𝐹𝐿 < 8350 ∩ 𝑅𝑎𝑛𝑔𝑒 > 2952) is greater than 0.9, and compared the probability-space surrogate reliability prediction against that from the fully coupled analysis. It can be observed that the average of the magnitude percentage difference in the reliability prediction is less than 1.5% when the number of training points is more than 250. This figure illustrates the ability of the probability-space surrogate to predict the reliability values required to carry out further optimization analysis. Since the generation of training points is associated with computational expense, we chose 1500 training points (associated with the first locally minimum value in Fig. 17b) in further analysis. The optimum number of components in the GMM are obtained by minimizing the BIC score, which was found to be 5. The BIC scores corresponding to MG, GC and GMM trained with 1500 training points are 1.415 × 105 , 1.416 × 105 and 1.294 × 105 respectively. Among the three surrogates, the GMM model provides the best fit to the training data (lowest BIC score). Since GMM has the lowest BIC score, it can be used for optimization analysis. However, we performed optimization analysis with all the three surrogates and compared their performance against the

39

fully coupled analysis. Multiple single optimization problems are solved using the DIRECT [34] global optimizer by varying the values of 𝛽 to obtain a Pareto surface. The performance of various probability-space surrogates are also compared against a Gaussian Process (GP) (also known as Kriging) algebraic surrogate and the fully coupled Monte Carlo analysis. A GP surrogate is used for its ability to capture non-linear relationships between variables, and for its ability to provide an uncertainty estimate of the prediction in addition to a mean prediction; this uncertainty is also considered for a comprehensive uncertainty analysis [36]. Also, GP models have been previously used for multidisciplinary optimization and optimization under uncertainty [6,37]. Note that the training points for an algebraic surrogate (GP) should be the results of fully converged multi-disciplinary analysis. For a fair comparison of the performance of various probability-space surrogates and the GP, we considered the same number of function evaluations in generating training points. Since we considered the second iteration values of the multidisciplinary analysis as training points for the probability-space surrogates, the total number of function evaluations are equal to 3000 (1500 training points in the second iteration). On an average, we observed that the multidisciplinary analysis reaches convergence in about 6.4 iterations. Therefore, we considered 469 (3000/6.4) points for training the GP algebraic surrogate. The Pareto surfaces obtained using the three probability-space surrogates, GP algebraic surrogate, and through the fully coupled analysis are shown in Fig. 19. It can be observed that the Pareto surfaces using all the three probability-space surrogates are comparable to the Pareto surfaces using the GP algebraic surrogate and fully coupled analysis; these results can be used to infer that all the three probability-space surrogates and the GP surrogate are able to capture the expectation values of the objective functions very well. Next,

40

we consider the same case study and solve an RBRDO problem, where both the expectation and standard deviation of the objective function are considered.

Fig. 19. Pareto surfaces using various probability-space surrogates, Gaussian process algebraic surrogate, and fully coupled analysis for RBDO with 1500 training points B. Reliability-based robust design optimization (RBRDO) Similar to the RBDO problem, an RBRDO problem is formulated to minimize the expectation and standard deviation of block fuel subject to a joint reliability constraint involving 𝑉𝑎𝑝𝑝 , 𝑇𝑂𝐹𝐿, 𝐷𝑁𝐴𝐶 and 𝑅𝑎𝑛𝑔𝑒 defined as 𝑀𝑖𝑛𝑖𝑚𝑖𝑧𝑒 𝐸 [𝐵𝐹 ]𝑎𝑛𝑑 𝑀𝑖𝑛𝑖𝑚𝑖𝑧𝑒 𝑆𝑡𝑑[𝐵𝐹] (17) 𝑠. 𝑡 Pr(𝑉𝑎𝑝𝑝 < 𝑉𝑎𝑝𝑝,𝑡ℎ ∩ 𝐷𝑁𝐴𝐶 < 𝐷𝑁𝐴𝐶𝑡ℎ ∩ 𝑇𝑂𝐹𝐿 < 𝑇𝑂𝐹𝐿𝑡ℎ ∩ 𝑅𝑎𝑛𝑔𝑒 > 𝑅𝑎𝑛𝑔𝑒𝑡ℎ ) > 𝛾

Here, the threshold values of 𝑉𝑎𝑝𝑝 , 𝑇𝑂𝐹𝐿, 𝐷𝑁𝐴𝐶 and 𝑅𝑎𝑛𝑔𝑒 are assumed to be equal to 145 𝑘𝑡𝑠, 8350 𝑙𝑏𝑠, 7.9 𝑓𝑡 and 2952 𝑚𝑖𝑙𝑒𝑠 respectively, and 𝛾 is equal to 0.9. Note that for the RBDO formulation, 𝑅𝑎𝑛𝑔𝑒 variable is considered in the objective function while in RBRDO formulation,

41

it is considered as a constraint. The multi-objective function is converted to a single-objective function as shown in Section IV.A as

𝑀𝑖𝑛𝑖𝑚𝑖𝑧𝑒 𝑂 = 𝛽 ×

𝐸 [𝐵𝐹 ] 𝑆𝑡𝑑[𝐵𝐹] 𝑚 + (1 − 𝛽 ) × 𝑠 𝐵𝐹𝑏𝑙 𝐵𝐹𝑏𝑙

(18)

where 𝛽 and 1 − 𝛽 represent the weights corresponding to the expectation and standard deviation 𝑚 𝑠 of block fuel. 𝑆𝑡𝑑[. ] represents the standard deviation operator. 𝐵𝐹𝑏𝑙 and 𝐵𝐹𝑏𝑙 represent the

normalization factors and are assumed to be equal to 10000 and 100 respectively. These values of normalization factors are chosen to make the mean and standard deviation be in the same order of magnitude for considering weighted combination in the objective function. Typically, the normalization factors are chosen by the modeler. The design solutions using various probabilityspace surrogates, GP algebraic surrogate, and fully coupled analysis are shown in Fig. 20. From Fig. 20, it can be observed that the Pareto surface corresponding to the GMM is closer to that of the fully coupled analysis. The Pareto surfaces corresponding to the MG and GC show higher standard deviations (Y-axis values in Fig. 20) when compared against the fully coupled analysis. Similar observations have been made in our earlier work [38], where an RBRDO analysis was carried out for an aircraft wing airfoil design. Thus, it can be concluded from this illustrative example that the GMM performs better when compared against the Multivariate Gaussian and Gaussian Copula when the standard deviation is also considered in the optimization formulation such as in RBRDO. Also, we can observe that the Pareto surface associated with the GMM is much closer to that of the fully coupled Monte Carlo analysis when compared against the Pareto surface associated with the GP algebraic surrogate. Therefore, we can infer from this illustrative example that the GMM is better able to capture the standard deviation of the objective function

42

(Block Fuel) when compared to the GP algebraic surrogate built with the same number of function evaluations.

Fig. 20. Pareto surfaces using various probability-space surrogates, Gaussian process algebraic surrogate, and fully coupled analysis for RBRDO with 1500 training points

The above RBDO and RBRDO analyses were carried out for a target reliability of 0.9. To study the performance of the probability-space surrogates for higher target reliabilities, both the RBDO and RBRDO analyses were repeated for target reliabilities of 0.95 and 0.99 (𝛾 in Equation 17). Figs. 21 and 22 provide the Pareto surfaces of the RBDO and RBRDO analyses using various probability-space surrogates, GP algebraic surrogate, and fully coupled analysis.

(a)

(b)

43

Fig. 21 Pareto surfaces using various probability-space surrogates, GP surrogate, and fully coupled analysis for RBDO using 1500 training points and with targeted reliabilities of (a) 0.95 and (b) 0.99

(a) (b) Fig. 22 Pareto surfaces using various probability-space surrogates, GP surrogate, and fully coupled analysis for RBRDO analysis using 1500 training points and with targeted reliabilities of (a) 0.95 and (b) 0.99

From Fig. 21, it was observed that various surrogates (probability-space and algebraic) were able to capture the Pareto surface of the fully coupled analysis well when the target reliability is 0.95 (Fig. 21 (a)); however, when the target reliability is 0.99, the end of the Pareto surface (when the expected value of the Block Fuel is about 5.5 × 104 in Fig. 21(b)) was not well captured by the surrogates. This can be attributed to the high target reliability requirement in the optimization formulation. From Fig. 22, it can be observed that the accuracy of all the surrogates to capture the Pareto surface of the fully coupled analysis has reduced when compared against Fig. 20 (target reliability of 0.9) due to the high target reliabilities of 0.95 and 0.99 when compared against the target reliability of 0.9 in Fig. 20. The Pareto surface of the GMM is closer to that of the fully coupled analysis when compared against the GP surrogate. Therefore, it can be concluded from this study that the Gaussian Mixture Model can be used to carry out both RBDO and RBRDO analyses in the presence of high target reliabilities, and its performance is better when compared to the GP surrogate and other probability-space surrogates. 44

Moreover, to compare the performance of the probability-space surrogates with that of the GP algebraic surrogate with respect to the number of training points, we repeated the RBDO and RBRDO analyses with 500, 1000, and 2000 training points, and also investigated the variation of reliability prediction performance of the probability-space and GP algebraic surrogates with the number of training points.

Fig. 23. Variation in the reliability prediction performance using various probability-space surrogates and GP algebraic surrogate with the number of training points

Fig. 23 provides the variation of the average of the magnitude percentage difference in the reliability prediction (computed using Eq. (16)) when compared against fully coupled Monte Carlo analysis using the three probability-space surrogates and also the GP surrogate at the same 25 validation points used in Fig. 18. It can concluded from this figure that the reliability prediction performance of the GMM is higher than the Multivariate Gaussian, Gaussian Copula, and also the GP surrogate. The average of the magnitude percentage difference of the GP surrogate has reduced the most among all the surrogates (about 3% error at 500 training points to less than 2% error at 2000 training points); however, its values are more than those of the GMM. The magnitude percentage difference of the Multivariate Gaussian model has reduced with increased training points, and this could be due to its inability to capture the all non-linear relationships between variables using the mean vector and covariance matrix of the Multivariate Gaussian distribution.

45

(a)

(b)

(c) Fig. 24. Pareto surfaces using various probability-space surrogates, GP algebraic surrogate, and fully coupled analysis for RBDO with (a) 500, (b) 1000, and (c) 2000 training points

(a)

(b)

(c)

46

Fig. 25. Pareto surfaces using various probability-space surrogates, GP algebraic surrogate, and fully coupled analysis for RBRDO with (a) 500 (b) 1000 and (c) 2000 training points

Figs. 24 and 25 provide the Pareto surfaces of the RBDO and RBRDO analyses using various probability-space surrogates, GP surrogate and fully coupled analysis. From Fig. 24, it can be observed that various surrogates (probability-space and GP algebraic) were able to capture the Pareto surface obtained from fully coupled analysis. In Fig. 25, when the number of training points was 500, both the probability-space surrogates and GP surrogate were unable to capture the Pareto surface obtained from fully coupled analysis. As the number of training points increased, the performance of the surrogates also generally improved; however, the rate of improvement in the GMM performance was higher than that of the GP surrogate since the GMM Pareto surface became closer to that obtained from fully coupled analysis faster than the GP Pareto surface. At 1000 training points, the performance of the GMM has significantly improved when compared to that of the GP surrogate, and this can be concluded from their corresponding Pareto surfaces. Similar improvements in performance were observed with 1500 training points and also with 2000 training points. From this study, we can conclude that the GMM probability-space surrogate was able to achieve better optimization results when compared to that of the GP algebraic surrogate with the same number of training points. Discussion: When the training samples cannot be fit with any of the proposed probabilityspace surrogates with desired accuracy, then other types of probability-space surrogates with computationally expensive sampling-based inference methods can be used, such as Bayesian networks, other types of copulas (elliptical and Archimedean) [39] and Vine Copulas[40]. In some cases, the uncertainty and reliability analysis results can be significantly impacted by the presence of tail dependence between random variables. Before using the proposed surrogates, 47

their ability to capture any tail dependence (if available) needs to be investigated. Since the proposed surrogates (Multivariate Gaussian, Gaussian Copula and Gaussian Mixture Models) are learnt from the training samples, and if there is an insufficient number of samples in the tails of the random variables, then it may not be possible to accurately capture the tail dependence. If there is a sufficient number of training samples in the tail regions, then the Gaussian Mixture Model is best suited to capture such tail dependence (compared to Multivariate Gaussian and Gaussian Copula) since there exists flexibility in choosing the number of components in the Gaussian mixture. Additional Gaussian components may be chosen in the tail regions of the random variables to capture the tail dependence accurately. If adding more Gaussian components still cannot capture the tail dependence with desired accuracy, then other types of probability-space surrogates with computationally intensive sampling-based inference can be used, such as Vine Copulas [39,40] or Bayesian networks. In high-dimensional systems where the probability-space surrogate cannot accurately capture the nonlinear relationships between the variables, several dimension reduction techniques can be applied first before constructing the probability-space surrogate. These techniques include principal component analysis of the output to map the high-dimensional output to a lower dimensional space [26], variance-based global sensitivity analysis to reduce the dimension of the uncertain variables [41], and the active subspace technique to map the high-dimensional input to a lower dimensional space [42]. 5. Conclusion This paper developed a probability-space surrogate modeling approach to perform efficient multidisciplinary design optimization under uncertainty. Multidisciplinary optimization under uncertainty, in general, is computationally expensive as it involves a nested three-loop process 48

where the multidisciplinary analysis is performed in the inner loop, uncertainty and reliability analysis in the middle loop and optimization in the outer loop. In contrast to the commonly used algebraic surrogates in the physical space, this paper considers a probability-space surrogate, one run of which provides a distribution prediction considering the uncertain variables. Two options exist for building a probability-space surrogate: (1) build an accurate surrogate such as a Bayesian network, which has an approximate and computationally expensive sampling-based inference through Markov Chain Monte Carlo methods (MCMC), or (2) build an approximate surrogate but with analytical inference (fast and exact) such as a multivariate Gaussian, Gaussian copula or Gaussian mixture; the second option was investigated in this paper. The training points for building a probability-space surrogate are obtained from a one-pass analysis through the individual disciplines in a multidisciplinary system. The probability-space surrogate is then trained, and used to infer the distributions of the coupling variables, the objectives, and the constraints by imposing the multidisciplinary compatibility condition. By using the probability-space surrogate, the three-loop optimization analysis is collapsed into a two-loop analysis as the uncertainty analysis and multidisciplinary analysis are performed simultaneously. The proposed methodologies are demonstrated for both reliability-based optimization (RBDO) and reliability-based robust design optimization (RBRDO) of a simplified three-discipline aircraft design example. Further work is needed on scaling up the proposed surrogate modeling strategy to realistic problems with large numbers of coupled disciplines, high dimensionality of design and uncertain variables, and consideration of complicated uncertainty sources, including spatiotemporal variability [43] and model errors [44,45].

49

Acknowledgement The research reported in this paper was partly supported by funds from the Air Force Office of Scientific Research (Grant No. FA9550-15-1-0018, Program Officers: David Stargel and Jaimie Tiley). The support is gratefully acknowledged. The authors also thank Dr. Jimmy C. Tai and Dr. Dimitri Mavris at the Georgia Institute of Technology for providing the aircraft design problem statement and associated models, Dr. Thomas A. Zang from Zang M&S Consulting for discussions regarding the aircraft design problem, and Dr. Roger Cooke from Resources for the Future (RFF) for discussions regarding copula models, during an earlier research project funded by Airbus (Contact: Dr. Sanjiv Sharma), and the National Institute of Aerospace (Contact: Dr. David Peake). Credit Author Statement Saideep Nannapaneni: Methodology, Data Curation, Formal Analysis, Software, Visualization, Writing – Original draft preparation Sankaran Mahadevan: Conceptualization, Research supervision, Resource acquisition and deployment, Writing - Review & Editing, Project administration

Declaration of interests The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

References [1]

Löhner, R., Yang, C., Cebral, J., Baum, J., Luo, H., Pelessone, D., and Charman, C., 1998, “Fluid-Structure-Thermal Interaction Using a Loose Coupling Algorithm and Adaptive Unstructured Grids,” Proceedings of the 29th AIAA Fluid Dynamics Conference.

[2]

Morgenthal, G., and McRobie, A., 2002, “A Comparative Study of Numerical Methods for Fluid Structure Interaction Analysis in Long-Span Bridge Design,” Wind & Structures, 5(2), pp. 101–114.

[3]

Jin, R., Du, X., and Chen, W., 2003, “The Use of Metamodeling Techniques for Optimization under Uncertainty,” Structural and Multidisciplinary Optimization, 25(2), pp. 99–116.

50

[4]

Koch, P., Wujek, B., Golovidov, O., and Simpson, T., 2002, “Facilitating Probabilistic Multidisciplinary Design Optimization Using Kriging Approximation Models,” Proceedings of the 9th AIAA/ISSMO Symposium on Multidisciplinary Analysis and Optimization.

[5]

Paiva, R. M., D. Carvalho, A. R., Crawford, C., and Suleman, A., 2010, “Comparison of Surrogate Models in a Multidisciplinary Optimization Framework for Wing Design,” AIAA Journal, 48(5), pp. 995–1006.

[6]

Simpson, T. W., Mauery, T. M., Korte, J., and Mistree, F., 2001, “Kriging Models for Global Approximation in Simulation-Based Multidisciplinary Design Optimization,” AIAA Journal, 19(12), pp. 2233–2241.

[7]

Ghosh, S., and Mavris, D. N., 2016, “A Methodology for Probabilistic Analysis of Distributed Multidisciplinary Architecture (PADMA),” Proceedings of the 17th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference, p. 3210.

[8]

Du, X., and Chen, W., 2005, “Collaborative Reliability Analysis under the Framework of Multidisciplinary Systems Design,” Optimization and Engineering, 6, pp. 63–84.

[9]

Mahadevan, S., and Smith, N., 2006, “Efficient First-Order Reliability Analysis of Multidisciplinary Systems,” International Journal of Reliability and Safety, 1, pp. 137– 154.

[10] Du, X., Guo, J., and Beeram, H., 2008, “Sequential Optimization and Reliability Assessment for Multidisciplinary Systems Design,” Structural and Multidisciplinary Optimization, 35(2), pp. 117–130. [11] Sankararaman, S., and Mahadevan, S., 2012, “Likelihood-Based Approach to Multidisciplinary Analysis Under Uncertainty,” Journal of Mechanical Design, 134(3), p. 031008. [12] Haldar, A., and Mahadevan, S., 2000, Probability, Reliability and Statistical Methods in Engineering Design, John Wiley & Sons Inc, New York, NY. [13] Hu, Z., and Mahadevan, S., 2019, “Global Sensitivity Analysis Using Efficient Distribution Surrogates,” AIAA Scitech 2019 Forum, p. 0440. [14] Liang, C., 2016, “Multidisciplinary Analysis and Optimization under Uncertainty,” PhD Dissertation, Vanderbilt University. [15] Brooks, S., 1998, “Markov Chain Monte Carlo Method and Its Application,” Journal of the Royal Statistical Society: Series D (the Statistician), 47(1), pp. 69–100. [16] Noh, Y., Choi, K., and Du, L., 2009, “Reliability-Based Design Optimization of Problems with Correlated Input Variables Using a Gaussian Copula,” Structural and Multidisciplinary Optimization, 38, pp. 1–16. [17] Liang, C., and Mahadevan, S., 2016, “Multidisciplinary Optimization under Uncertainty Using Bayesian Network,” SAE International Journal of Materials and Manufacturing, 9, pp. 419–429.

51

[18] Liang, C., and Mahadevan, S., 2015, “Reliability-Based Multi-Objective Optimization Under Uncertainty,” Proceedings of the 16th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference. [19] Zuo, M., Jiang, R., and Yam, R., 1999, “Approaches for Reliability Modeling of Continuous-State Devices,” IEEE Transactions on Reliability, 48(1), pp. 9–18. [20] Sun, S., Zhang, C., and Yu, G., 2006, “A Bayesian Network Approach to Traffic Flow Forecasting,” IEEE Transactions on Intelligent Transportation Systems, 7(1), pp. 124– 132. [21] Yu, J., and Qin, S., 2008, “Multimode Process Monitoring with Bayesian Inference‐based Finite Gaussian Mixture Models,” AIChE Journal, 54(7), pp. 1811–1829. [22] Terejanu, G., Singla, P., Singh, T., and Scott, P., 2008, “Uncertainty Propagation for Nonlinear Dynamic Systems Using Gaussian Mixture Models,” Journal of Guidance, Control, and Dynamics, 31(6), pp. 1623–1633. [23] Zhang, H., Giles, C., Foley, H., and Yen, J., 2007, “Probabilistic Community Discovery Using Hierarchical Latent Gaussian Mixture Model,” AAAI, 7, pp. 663–668. [24] Zaman, K., and Mahadevan, S., 2013, “Robustness-Based Design Optimization of Multidisciplinary System Under Epistemic Uncertainty,” AIAA Journal, 51(5), pp. 1021– 1031. [25] Youn, B., and Xi, Z., 2009, “Reliability-Based Robust Design Optimization Using the Eigenvector Dimension Reduction (EDR) Method,” Structural and Multidisciplinary Optimization, 37(5), pp. 475–492. [26] Liang, C., and Mahadevan, S., 2016, “Stochastic Multidisciplinary Analysis with High Dimensional Coupling,” AIAA Journal, 54(4), pp. 1209–1219. [27] Nannapaneni, S., Mahadevan, S., and Rachuri, S., 2016, “Performance Evaluation of a Manufacturing Process under Uncertainty Using Bayesian Networks,” Journal of Cleaner Production, 113, pp. 947–959. [28] Bedford, T., and Cooke, R. M., 2001, “Probability Density Decomposition for Conditionally Dependent Random Variables Modeled by Vines,” Annals of Mathematics and Artificial Intelligence, 32(1–4), pp. 245–268. [29] Posada, D., and Buckley, T. R., 2004, “Model Selection and Model Averaging in Phylogenetics: Advantages of Akaike Information Criterion and Bayesian Approaches over Likelihood Ratio Tests,” Systematic Biology, 53(5), pp. 793–808. [30] Mahmoud, M. S., and Khalid, H. M., 2013, “Expectation Maximization Approach to Data-Based Fault Diagnostics,” Information Sciences, 235, pp. 80–96. [31] Yao, W., Chen, X., Luo, W., Van Tooren, M., and Guo, J., 2011, “Review of UncertaintyBased Multidisciplinary Design Optimization Methods for Aerospace Vehicles,” Progress in Aerospace Sciences, 47(6), pp. 450–479. [32] Koziel, S., and Yang, X., 2011, Computational Optimization, Methods and Algorithms,

52

Springer. [33] Jones, D. R., Perttunen, C. D., and Stuckman, B. E., 1993, “Lipschitzian Optimization without the Lipschitz Constant,” Journal of Optimization Theory and Applications, 79(1), pp. 157–181. [34] Finkel, D., 2005, “Global Optimization with the DIRECT Algorithm,” PhD Dissertation, North Carolina State University. [35] Liang, C., and Mahadevan, S., 2017, “Pareto Surface Construction for Multi-Objective Optimization under Uncertainty,” Structural and Multidisciplinary Optimization, 55(5), pp. 1865–1882. [36] Nannapaneni, S., Hu, Z., and Mahadevan, S., 2016, “Uncertainty Quantification in Reliability Estimation with Limit State Surrogates,” Structural and Multidisciplinary Optimization, 54(6), pp. 1509–1526. [37] Dubourg, V., Sudret, B., and Bourinet, J. M., 2011, “Reliability-Based Design Optimization Using Kriging Surrogates and Subset Simulation,” Structural and Multidisciplinary Optimization, 44(5), pp. 673–690. [38] Nannapaneni, S., Liang, C., and Mahadevan, S., 2017, “Bayesian Network Approach to Multidisciplinary, Multi-Objective Design Optimization under Uncertainty,” 18th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference, p. 3825. [39] Kurowicka, D., and Cooke, R., 2006, Uncertainty Analysis with High Dimensional Dependence Modelling, John Wiley & Sons. [40] Joe, H., Li, H., and Nikoloulopoulos, A. K., 2010, “Tail Dependence Functions and Vine Copulas,” Journal of Multivariate Analysis, 101(1), pp. 252–270. [41] Hu, Z., and Mahadevan, S., 2019, “Probability Models for Data-Driven Global Sensitivity Analysis,” Reliability Engineering & System Safety, 187, pp. 40–57. [42] Constantine, P. G., 2015, Active Subspaces: Emerging Ideas for Dimension Reduction in Parameter Studies, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, Pennsylvania, USA. [43] Hu, Z., and Mahadevan, S., 2017, “A Surrogate Modeling Approach for Reliability Analysis of a Multidisciplinary System with Spatio-Temporal Output,” Structural and Multidisciplinary Optimization, 56(3), pp. 553–569. [44] Mahadevan, S., and Rebba, R., 2006, “Inclusion of Model Errors in Reliability-Based Optimization,” Journal of Mechanical Design, 128(4), p. 936. [45] Nannapaneni, S., and Mahadevan, S., 2016, “Reliability Analysis under Epistemic Uncertainty,” Reliability Engineering & System Safety, 155, pp. 9–20.

53