Robust optimisation of computationally expensive models using adaptive multi-ﬁdelity emulation

Computationally expensive models are increasingly employed in the design process of engineering products and systems. Robust design in particular aims to obtain designs that exhibit near-optimal performance and low variability under uncertainty. Surrogate models are often employed to imitate the behaviour of expensive computational models. Surro-gates are trained from a reduced number of samples of the expensive model. A crucial component of the performance of a surrogate is the quality of the training set. Problems occur when sampling fails to obtain points located in an area of interest and/or where the computational budget only allows for a very limited number of runs of the expensive model. This paper employs a Gaussian process emulation approach to perform eﬃcient single-loop robust optimisation of expensive models. The emulator is enhanced to propagate input uncertainty to the emulator output, allowing single-loop robust optimisation. Further, the emulator is trained with multi-ﬁdelity data obtained via adaptive sampling to maximise the quality of the training set for the given computational budget. An illustrative example is presented to highlight how the method works, before it is applied to two industrial case studies.


Introduction
The chief aim of engineering design is to create systems that satisfy specific performance objectives and constraints over a period of time. Usually, there exist many feasible designs that satisfy the required objectives. For this reason, it is necessary to choose an optimal design according to some criterion. Modern engineering systems are inherently complex. This complexity means that endogenous (geometry, material properties) and exogenous (loads) information is never complete, and often varies throughout the life cycle of the system (e.g. degradation altering geometry, etc.). The objective of robust design is to determine a set of designs that exhibit high levels of performance with low variability, whilst taking uncertainties into account. The benefits of robust design include the assurance of high performance regardless of a variety of unknown factors and occurrences throughout the system's life cycle. Robust design is essentially a traditional optimisation task, but with an added constraint relating to the performance variability, or robustness, within some predefined neighbourhood of the input variables. There are various definitions of robustness. A detailed review of which is presented in Gabrel et al. [1] , leading to various methodologies for tackling the robust optimisation problem. The authors of [2] employed a reliability-based optimisation algorithm which utilised Monte Carlo integration to obtain an averaged performance value within the neighbourhood. Similarly, Ryan [3] employed a probability distribution estimation method to obtain an approximate distribution of the performance within the neighbourhood. Another approach utilised the Taylor expansion of the expectation and variance of the performance and attempted to minimise both criteria simultaneously. Alternatively, several papers chose to optimise the worst-case scenario rather than any sort of averaged performance [4,5] .
Typically, the behaviour of modern engineering systems is modelled by computationally expensive simulators, which can be seen as mappings from the input space to the output space, denoted f : x ∈ χ → y ∈ R . However, working directly with f ( x ) is often infeasible due to computational expense. A widespread approach to tackle this problem is to replace f ( x ) with a surrogate model, which has been trained using data obtained from a small amount of simulator evaluations. One option is to train a Gaussian Process Emulator (GPE), which is defined by a mean function and a covariance function respectively.
The mean function provides an inexpensive approximation to the simulator, η( x ) ≈ f ( x ) , whilst the covariance function provides a measure of output uncertainty at each set of inputs, V x [ η( x )] [6] . The result is that the robust design problem can be interpreted mathematically as Here x represents the set of input variables located within the hypercube, or neighbourhood, centered at x and bounded by x ± . Consequently, h j ( x ) and w ν ( x ) are the respective inequality and equality constraints of this neighbourhood. In this context, robust design is interpreted as a double-loop optimisation task, with the outer-loop optimising the overall performance, subject to the constraint functions, and the inner-loop optimising for robustness in neighbourhood of the input variables.
The ability of the GPE to accurately approximate f ( x ) is directly related to the quality of the training set. There are two main approaches to address this issue: adaptive sampling schemes and supplementing the training set with data from multiple levels of fidelity. Adaptive sampling approaches tend to involve a utility function which measures some form of model improvement to select additional sample points. The most popular choice is expected improvement [7] , which has been widely used in reliability [8] , optimisation [9] and robust optimisation problems [10] , amongst others. Further, the concept can be extended to multiple performance functions by considering the expected improvement of the current Pareto front via hypervolume expected improvement [11] . Other schemes include maximising the probability of improvement [12] or selecting samples with high uncertainty [13] . Multi-fidelity (MF) approaches are applicable when more than one potential simulator exists for the system under study. Lower-fidelity (LF) samples are defined by a lower-computational cost, but lower accuracy, than higher-fidelity (HF) samples. Multi-fidelity surrogate approaches exploit LF samples to gain information of the behaviour of the underlying system, and HF samples to maintain the desired accuracy. Most multi-fidelity approaches utilise LF data and adaptive sampling to attempt to sample the HF points in regions of interest and maximise the effectiveness of the surrogate [14][15][16] . Employing a surrogate model reduces the computational cost involved in robust design problems considerably. However, when there are a large number of performance functions and/or input variables, the double-loop approach becomes increasingly inefficient. A solution is to collapse the problem into a single-loop approach as done for a single-fidelity surrogate in Ryan [17] . In that paper, a GPE was enhanced to provide exact values of output uncertainty in the presence of uncertain inputs. This paper provides a framework to perform efficient robust design on computationally expensive models. The framework adapts the single-loop approach discussed above to factor in multiple levels of fidelity, and supplements it with a hybrid adaptive sampling scheme. The paper is organised as follows. Section 2 provides an overview of various forms of Gaussian process emulation. The proposed approach is introduced in Section 3 and discusses the main components. An illustrative example and two industrial CFD case studies are presented in Section 4 . The final section provides relevant conclusions and highlights future work.

Methodology overview
This section provides an overview of various forms of Gaussian process emulation. The main steps involved in the training of a single-fidelity (SF) GPE are described. Two extensions of the GPE framework are then discussed: training a GPE with MF data and factoring in input uncertainty for a SF GPE.

Single-fidelity Gaussian process emulation
Computationally expensive models are deterministic mappings from some input x ∈ R d to an output y = f ( x) . Due to the computational expense, often only a limited number of input samples can be evaluated. Under the Bayesian paradigm, f ( x) can be regarded as a random variable as the output is unknown until it is computed (and thus observed) by the modeller.
Gaussian process emulation follows a Bayesian framework to provide a statistical approximation of η( x ) ≈ f ( x ) . Initially, a Gaussian Process prior is placed on the output, in the form [18] η( (2) The first term is the mean of the emulator and provides the general trend; h ( x ) is a vector of q known regression functions, β is a vector of q unknown coefficients. The second term controls local behaviour; Z( x ) represents a Gaussian Process, with mean zero and covariance σ 2 c( x , x ; θ ) . Here σ 2 is a scalar parameter, and θ is a vector which specifies the smoothness of the inputs and ultimately dictates the behaviour of the correlation function Given a set of single-fidelity (i.e. data from only one model/simulator) training data, D = ( x 1 , y 1 ) , ( x 2 , y 2 ) , . . . , ( x n , y n ) , the emulator distribution, conditional on training data and the unknown hyperparameters β, σ , θ , is defined as where M(·) and C(·, ·) are the mean and covariance functions of the GP, respectively. Both β and σ 2 are assigned a prior and estimated via their respective maximum a posteriori estimates, ˆ β and ˆ σ 2 , while ˆ θ is often estimated via maximising a likelihood function [19] . Ultimately, the posterior distribution of at some unobserved input x * , conditional on the observed data, is given by Oakley [20] η( with posterior predictive mean function and posterior predictive covariance function Here is a matrix containing the correlation between each training point, t (x * ) is a vector containing correlation between x * and the training points, and H is a vector containing the training points evaluated at the regression function h (x ) . For the work in this paper, the Gaussian process prior is assumed to have mean zero, i.e. h (x ) T = 0 , which will be adopted from here onward.

Multi-fidelity Gaussian process emulation
Computationally expensive models are designed to capture the behaviour of an underlying physical system or product. It is often the case that more than one computational model is available, with each model corresponding to a varying degree of computational cost and accuracy. The models can usually be organised in levels of fidelity; a model with lower computational costs but less accuracy is considered to be of a lower fidelity than a more expensive and accurate model. For example, a 2-dimensional versus a 3-dimensional model, or a model with a coarse mesh versus one with a fine mesh. The concept of Gaussian Process Emulation can be extended to incorporate multiple levels of fidelity within the training process [21] , allowing for improved emulator performance or lower training costs. In the case where there are two levels of fidelity, denoted low-fidelity (LF) and high-fidelity (HF), this paper adopts the recursive multi-fidelity approach from [22] , to approximate the output from the HF model as: Here, η LF (x ) represents a GPE trained using data from the LF model, ρ LF ( x ) represents a regression function, and δ HF (x ) represents a Gaussian Process Emulator which models the discrepancy between the HF estimation of ρ LF ( x ) η LF ( x ) and the true HF simulator realisations. Both emulators are trained via the steps described in Section 2.1 . This can be generalised for t levels of fidelity in a recursive fashion:

Single-fidelity robust Gaussian process emulation
The two-looped approach to solving the robust optimisation problem (1) works by first attempting to minimise the objective functions η( x ) and V x [ η( x )] in the outer-loop. Once a potential solution is found, the inner-loop measures the robustness over the input distribution. As a result, the predictive distribution of the emulator given input uncertainty is found by marginalising over the input distribution: This marginalisation is the aforementioned inner-loop and often achieved via Monte Carlo sampling. In the case where the uncertainty within the inputs is normally distributed, i.e. for an unknown point x * ∼ N ( u , S ) , it is possible to extract the first and second moments of Eq. (10) via methods described in Quinonero-Candela et al. [23 , 24] . These moments provide analytical expressions for the mean, m ( u , S ) , and variance, v ( u , S ) , of p η( x * ) | u , S, D . Ultimately, having direct access to the mean and variance of the emulator conditional on the input uncertainty collapses the robust optimisation problem down to a single-loop: The full details and steps involved are discussed further in [17] . The resulting mean and variance functions are fundamentally Eqs. (6) and (7) corrected to factor in the input uncertainty. The added input uncertainty essentially flattens the output, with a decreased vertical amplitude and increased correlation.

Proposed approach
The goal of the proposed approach is to perform efficient robust optimisation of computationally expensive models. The method is a combination of the various forms of Gaussian process emulation discussed in the previous section, and is termed Multi-Fidelity Robust Gaussian Process Emulator (MF-RGPE). When employing a GPE for the purposes of robust optimisation, the two main considerations are the ability of the GPE to accurately portray the behaviour of the underlying expensive model, and the efficiency of the robust optimisation process. To address the former, the proposed approach utilises training data from multiple levels of fidelity obtained via an extension of the Expected Improvement (EI) criterion [7] to maximise the quality of the training set. To increase the efficiency, the proposed approach extends the robust GPE detailed in Section 2.3 to the multi-fidelity case. Further details of the steps are discussed in the following subsections.

Generating training samples
The framework begins with the design of experiment (DoE) of the LF model. Latin hypercube sampling (LHS) [25] is used as the space-filling algorithm to generate the initial samples. These samples are then evaluated on both the LF model and any relevant constraints functions, and are referred to as LF samples. To generate the initial HF samples, the LF samples are first sorted according to their objective and constraint values. A proportion of the top performing samples are selected to be part of the initial HF samples. The remaining initial HF samples are selected by filling the remaining space using a spacefilling algorithm. This is done to encourage sampling of high-interest areas, whilst not neglecting the general performance of the GPE elsewhere. The proportion used in this work was 20% of initial samples from the top performing LF samples and 80% resulting from the space-filling algorithm. The HF samples were then evaluated on both the HF model and any relevant constraint functions.

Constructing the MF RGPE
The MF RGPE provides an approximation of the HF output whilst considering input uncertainty, and is constructed in a similar fashion to the standard MF GPE described in Eq. (8) : Here, η RLF ( x ) represents a robust GPE trained using data from the LF model via the steps described in Section 2.3 and ρ RLF ( x ) represents a regression function. The last term, δ HF ( x ) , represents a Gaussian Process Emulator which models the discrepancy between the estimation of the output of the HF training data, without accounting for input uncertainty, i.e. ρ LF ( x ) η LF ( x ) , and the actual HF simulator output. In an industrial context, there will usually be a predetermined computational budget and the stopping criterion will be met once this budget is exhausted. Other stopping criterion may include reaching a certain threshold of performance, such as obtaining a suitable design or reducing overall GPE uncertainty below some required value.

Adaptive sampling
For a SF GPE, the expected improvement (EI) [7] at some point x is defined as Here y min represents the current best performing objective value amongst the training data, s denotes the standard deviation of the GPE and (·) and φ(·) represent the cumulative and probability density functions of a standard Gaussian random variable, respectively. EI attempts to locate samples that offer improved nominal performance against the current best sample. The method balances higher probability of a relatively small improvement (exploitation) versus a lower probability of high improvement (exploration). The concept of EI can also be applied to cases with more than one objective function by considering a hypervolume of improvement. Following the steps described in Li et al. [26] , given some reference point r , the HVEI at some point x against a set of solutions U is defined as where N U is the number of points in the set U and N Ob j is the number of objective functions. Both EI and HVEI are designed for optimising nominal performance. The proposed approach extends them for the purposes of robust design by employing the mean and standard deviation GPE output from the SF RGPE and MF RGPE within the EI process. Consequently, the robust EI and robust HVEI can therefore be defined as where y R min is the current best performing objective value amongst the training data whilst also taking input uncertainty into account. Fig. 1 illustrates the concept of robust HVEI in the case of two objective functions. The set of solutions, P , represents the robust Pareto solutions taken from the current training data. Utilising these solutions and some reference point r , a set of local upper bounds U can be constructed such that U i lies at the intercept of P i and P i +1 . The robust HVEI is thus the summation of the robust EI against each local upper bound, to give an overall value of improvement. To obtain promising adaptive samples, an algorithm known as subset simulation is used to explore the input space and locate samples with high values of robust HVEI. Subset simulation [ 27 ] is an efficient Monte Carlo technique that employs Markov Chain Monte Carlo (MCMC) and a simple evolutionary strategy to converge to and sample from extremely small, or rare, regions of the input domain. Samples within these rare regions are associated with superior performance according to some user-defined criteria of interest. Originally developed and applied to reliability problems with great success [28,29] , the authors of [30][31][32] proposed the analogy between sampling from a small failure region and sampling from a small region exhibiting high performance. Consequently, the algorithm was adapted for the purposes of optimisation, allowing it to be used to locate promising samples even in the case where the areas offering improvement are extremely rare. Moreover, it is suitable when dealing with high-dimensional problems. To increase efficiency, several samples are adaptively sampled in one optimisation iteration. An influence function [33] , denoted τ ( x ) , is employed to discourage the adaptive samples clustering in one area, by scaling the robust EI values after each new adaptive sample is taken, i.e.
where x AS represents the latest adaptive sample.

Robust design
Once the computational budget is exhausted and the final batch of adaptive sampling completed, the final MF RGPEs can be utilised for robust design. Subset simulation is employed to locate the input regions corresponding to samples with high performance according to the MF RGPEs. These samples should be insensitive to perturbation in the values of the input variables, and given the computational resources, be validated on the HF model. The steps involved the the proposed method are outlined in a flowchart provided in Fig. 2 . Steps 1-4 involve generating the LF and initial HF training data, and are described in Section 3.1 . This provides the foundation for the construction of the initial MF RGPE in step 5, which is detailed in Section 3.2 . Provided the stopping criterion (step 6) has not been met, this MF RGPE is then used as a tool to attempt to locate samples with improved performance in step 7, using the adaptive sampling process from Section 3.3 . This procedure is repeated on a loop, with an improved MF RGPE constructed at each generation until the stopping criterion is met. Optimisation of the MF RGPE(s) takes place in step 8.

Numerical examples
This section provides three examples showcasing the MF RGPE approach discussed in the previous section. A synthetic example is first presented to showcase the concept of the approach before it is applied to two industrially-relevant test cases. In all examples, the regression function ρ RLF from Eq. (12) is set to one, as in each example there is no assumed prior knowledge regarding the relationship between the LF simulator output and HF simulator output.

Synthetic example
The motivation behind this synthetic example was to illustrate the main concepts of the proposed approach. The HF and LF functions are defined as The functions were constructed such that the LF function exhibited similar behaviour to the HF function, and as such could be used to infer regions of high interest. Additionally, both were designed to possess two maxima; a global maximum that was more sensitive to input uncertainty, and a local more robust maximum. The goal is to maximise f HF in the face of some input uncertainty, with the intention to favour the more robust local maximum. An initial batch of 50 LF samples were selected via LHS. The 4 samples with the highest objective values were then selected, alongside 16 further samples from LHS, to populate the HF training set. A MF RGPE was then constructed with training data normalised between 0 and 1, and input uncertainty for an unknown point x * is defined by the probability distribution x * ∼ N ( u , diag[0 . 01 , 0 . 01]) . Here u is the mean approximation of x * while diag[0 . 01 , 0 . 01] is a diagonal matrix containing the variance with respect to each input variable. A further 3 samples were obtained via the robust EI adaptive sampling algorithm, and the retrained MF RGPE employed for robust optimisation. Finally, the inputs were transformed back to their original domains.  2 ) , ( 3 π 2 ) and a local, more robust, optimum around x ≈ ( π 2 ) , ( π 2 ) in the bottom left. The adaptive samples all lie within these regions of interest, with a preference to the local, more robust optimum. The local optimum in the bottom left is considered more robust as it has a wider base, meaning there is a lower drop in performance given any perturbation in the inputs. Note that several of the initial batch of HF samples were already in proximity to the two optima, highlighting the importance of utilising the best performing LF samples. Further, the LF data provided valuable information in the regions where HF samples were sparse (e.g. top left), saving an adaptive sample being wasted in an area of low interest. The illustrative example was repeated 10,0 0 0 times, and the normalised error from the true robust optimum presented in Fig. 4 . The error was normalised to illustrate the discrepancy between the true robust optimum and the actual values more clearly. The goal of the study was to showcase the individual steps described in Fig. 2 and illustrate the merits of the approach. Overall, the majority of cases were within 1% of the true robust optimal input values.

Industrial examples
Design engineers often utilise computationally expensive models in their design process. It is often desirable to factor input uncertainty into this process. The proposed approach has been designed to assist design engineers in this task, within a reasonable computational budget. Computational Fluid Dynamics (CFD) [34] models are a common tool in engineering design. They are usually computationally expensive, which limits their ability to be used directly in practical applications, but makes them a prime candidate for the MF RGPE approach.

Turbulated duct case study
A frequent feature in turbine blades is the presence of turbulated internal cooling ducts. The presence of rib turbulators repeatedly perturb the boundary layer, which can result in significant heat transfer by promoting convective mixing with the core cooling flow. A downside is that this heat transfer comes at the cost of higher-pressure drop [35] . However, due to manufacturing constraints and degradation during the life cycle, the duct will likely diverge from the initial design at some point. The challenge is therefore to select a design that maximises heat flow, in this case Nusselt number, whilst minimising pressure drop in the face of input uncertainty. To address this challenge, a model of the turbulated duct was constructed using ANSYS software according to four geometric parameters that control the cross-sectional profile and angle of the turbulators, as shown in Fig. 5 . The range of parameter values are shown in Table 1 in the appendix. Within the ANSYS software, each combination of these four parameters would result in a unique turbulated duct geometry. This geometry was then meshed using an unstructured tetrahedral grid and solved using Reynolds-averaged NavierStokes (RANS) [34] to output the Nusselt number and pressure coefficient for that particular design.
For the MF RGPE approach, the overall computational budget assigned was equivalent to 44 HF samples. The LF model consisted of a mesh of approximately one million elements and solved using k − ω SST RANS on ANSYS, whilst the HF model consisted of a mesh of approximately five million elements and solved using k − ω SST RANS. The approximate computational cost was 20 LF samples ≈ 1 HF sample. Consequently, two separate MF RGPEs were trained for the Nusselt number and pressure coefficient respectively, with the initial MF RGPEs trained using 80 LF samples and 20 HF samples. A further 20 HF samples were adaptively selected in two batches of 10 samples to supplement the training set. The training data was normalised between 0 and 1, with input uncertainty represented via a diagonal matrix S = diag [0 . 025 , 0 . 025 , 0 . 025 , 0 . 025]) . The final MF RGPEs were then optimised using a multi-objective subset simulation algorithm. For comparison, the case study was repeated with the same computational budget, but using only HF samples to construct two SF RGPEs. The initial SF RGPEs were trained using 40 HF samples. A further 4 HF samples were adaptively selected in four batches of a single sample to supplement the training set. Fig. 6 demonstrates the adaptive sampling process and the final Pareto front for the MF RGPE approach. On inspection, several of the initial HF samples (green dots, top left plot) were located in close proximity to the eventual Pareto front, again highlighting the advantages of incorporating the LF training data to locate regions of interest. Indeed, the general performance of the adaptive points (blue stars) is significantly better than that of the randomly sampled points, showcasing the benefits of adaptive sampling. Furthermore, by comparing the two batches of adaptive sampling, it is clear how the adaptive sampling process attempts to converge towards the true Pareto front. It should be noted that the adaptive sampling procedure was assisted by the LF data to discard areas of low interest. The optimisation process placed constraints on the output variance of the respective GPEs to ensure a certain level of performance. This is highlighted in the close proximity between the Pareto front and the best performing training samples, placing further emphasis on the importance of quality training data. A single training point lies above the Pareto front, however this point was deemed to lack the necessary robustness according to the final MF RGPEs. Fig. 7 contains the Pareto fronts from the MF RGPE approach (red stars) and the SF RGPE approach (blue dots). In general, the two Pareto fronts possess similar behaviour, although the MF RGPE Pareto solutions exhibit superior performance than the SF RGPE Pareto solutions. A potential contributor to this discrepancy is the fact that the SF approach was made up  of a higher proportion of randomly sampled training data. However, it is standard practice to use a budget of at least ten training samples per input dimension [36] in order to have sufficient confidence in the output of the underlying surrogate. The MF RGPE approach circumvents this issue by utilising LF data to make up for any loss of information. As such, it is reasonable that increasing the proportion of adaptively sampled data in the SF RGPE case would not necessarily improve the performance, due to surrogate inaccuracy and added uncertainty disrupting the sampling process. A second contributor is the fact that the MF RGPE approach was able to infer regions of high interest from the LF training data to aid in the adaptive sampling procedure. Overall, the MF RGPE offered superior performance than the SF alternative for the same computational budget.

Aerofoil case study
The aerofoil test case involved obtaining a set of aerofoil solutions that maximise lift-to-drag ratio whilst minimising maximum blade thickness of a turbine blade in the face of potential perturbation of input values caused by uncertainty. A prospective aerofoil geometry was defined using the Class-Shape Transformation (CST) method [37] . In particular the Au and Al parameters are the weighting coefficients that help prescribe the thickness/shape at various locations along the upper and lower surfaces respectively. The parameters and their respective ranges are shown in Table 2 in the appendix. The LF model consisted of the aerofoil being solved over a range of angles of attack in XFOIL software, which performed a potential flow calculation without taking into account viscosity or a boundary layer. The HF model consisted of the aerofoil being solved via k-ω RANS in ANSYS. Unlike the turbulated duct case study, where the level of fidelity was solely due to mesh resolution, the fidelity in this case is dictated by two separate methods of varying accuracy and cost. It should be noted that the the definition of varying levels of fidelity is problem specific, with the only requirement that they exhibit similar behaviour in attempting to model the same underlying phenomena. The computational budget for the test case was approximately 240 HF samples. The comparative computational costs were approximately 20 LF samples per HF sample. Two separate MF RGPEs were trained for the lift-to-drag ratio and maximum thickness respectively, with the initial MF RGPEs trained using 600 LF samples and 120 HF samples. A further 80 HF samples were adaptively sampled in four batches of 20 to supplement the training set. The training data was normalised between 0 and 1, with input uncertainty represented via a 20 ×20 diagonal matrix S with each of the entries equal to 0 . 025 . The final MF RGPEs were then optimised using a multi-objective subset simulation algorithm. As in the previous example, the case study was repeated using the same computational budget comprising of only HF samples. The initial SF RGPEs were trained using 200 HF samples, with a further four batches of a 10 samples added via adaptive sampling.  Fig. 9 presents the adaptive sampling process and the final MF RGPE Pareto front. As in the turbulated duct study, several of the initial HF samples (green dots, top left plot) exhibited high performance and there was a clear convergence towards the suspected true Pareto front as the number of adaptive samples increased. The Pareto front closely followed the path of the best performing training samples. The training sample with the lowest maximum thickness was omitted from the Pareto front, as the performance of this point was particularly sensitive to input perturbations. Fig. 10 contains the Pareto fronts from the MF RGPE approach (red stars) and the SF RGPE approach (blue dots). Both Pareto fronts initially rise relatively sharply before reaching a plateau with respect to the lift-to-drag coefficient. However, there is a significant discrepancy between the respective performance of the two Pareto fronts, with that of the MF RGPE completely dominating the SF RGPE counterpart. Whilst the SF RGPE approach wasted a number of samples searching in uncertain but ultimately low interest areas, the MF RGPE approach was able to discard these areas and target more promising locations due to the information provided by the LF data. Fig. 11 displays 20 validation samples for 4 designs taken from the MF RGPE Pareto front. As in the turbulated duct case study, the number of validation samples were limited due to the computational costs involved. The validation samples were selected to verify the performance of the MF RGPE approach across the Pareto front. It should be noted that there was zero discrepancy between the LF simulator and HF simulator output for the maximum thickness. As a result, there was significantly less GPE uncertainty for this objective, and the majority of the uncertainty bounds with respect to the maximum thickness is due to input uncertainty. Each validation sample was within the 2 σ uncertainty bounds of the original Pareto solution. Moreover, the aerofoils corresponding to the validation points for the third (in ascending order of L/D ratio) Pareto solution were plotted against the original for a visual depiction of input uncertainty.

Conclusion
A Gaussian process emulation approach to perform efficient single-loop robust optimisation of expensive models, denoted MF RGPE, was presented. The approach combines various enhancements of the Gaussian process emulation approach, utilising MF training data and factoring input uncertainty into the output of the emulator. MF RGPE addresses the two main issues found when employing emulation-based approaches for robust optimisation, namely the quality of the emulator training set, and the efficiency of the robust optimisation process itself. Provided lower-fidelity simulators exhibit similar behaviour to their higher-fidelity counterparts (as expected), the approach offers improvement over single-fidelity methods. This is due to the increased information available in both the adaptive sampling phase and in estimating the performance in areas without HF training data. In the situation where the LF simulator exhibits dissimilar, or even misleading behaviour, the approach essentially reverts to a SF GPE, albeit with a slightly reduced training set. Compared to other MF methods, MF RGPE offers increased efficiency in performing robust optimisation, by collapsing the problem down to a single-loop. An illustrative example highlighted some of the key concepts of the approach before two industrial test cases demonstrated its ability to outperform produce quality results. Future work involves augmenting the adaptive sampling regime with the ability to choose both the location and fidelity of an adaptive sample, as well as applying the approach to cases with more than two levels of fidelity.