# Multi-Fidelity Gradient-Based Strategy for Robust Optimization in Computational Fluid Dynamics
## 3 Uncertainty Quantification Methods
### Bayesian Kriging and Gradient-Enhanced Kriging
#### 3.1.1 Bayesian Kriging  
A surface-response method based on BK was implemented to achieve a non-intrusive estimate of the cost function statistics. The QoI \(\mathbf{J}\) is modeled as a regression function of the form:  
\[\mathbf{J}(\mathbf{\xi})=\sum_{i}f_{i}(\mathbf{\xi})\beta_{i}+Z(\mathbf{\xi})=m(\mathbf{\xi})+Z( \mathbf{\xi}) \tag{5}\]  
where we dropped the dependency on the deterministic quantities for simplicity, so that \(\mathbf{\xi}\) defines a generic point in the \(n_{unc}\)-dimensional uncertain space, \(f_{i}\) are basis functions of the regression model, \(\beta_{i}\) are regression coefficients to be determined, and \(Z(\mathbf{\xi})\) is a zero mean Gaussian process (GP) [60], modeling the deviation between the regression function and the data. Choosing a zeroth-order polynomial as the regression model leads to the so-called ordinary kriging. The deviation \(Z\) is represented as a multi-variate normal distribution \(Z\sim N(0,R)\), with zero mean and covariance matrix \(\mathbf{R}\). The kriging approximation can be formulated in a Bayesian statistical framework. In this case, the BK surrogate predicts a set of \(M\) values of the QoI \(\mathbf{J}(\mathbf{\xi})\), conditional on \(N\) observed data \(\mathbf{J}^{\star}\), where \(\mathbf{J}^{\star}\) is selected from the full set of \(J\) through the observation matrix \(\mathbf{H}\) of size \(N\times M\): \(\mathbf{J}^{\star}=\mathbf{H}\mathbf{J}\) with:  
\[H_{ij}= \begin{cases}1,&\text{if}\qquad i=j\\ 0,&\text{otherwise}\end{cases}  \tag{6}\]  
for \(i=1,...,M\) and \(j=1,...,N\). The vector of the unknown QoI \(\mathbf{J}(\mathbf{\xi})\) is supposed to follow a normal prior distribution \(p(J)\), \(\mathbf{J}{\sim}N(0,\mathbf{P})\), where the covariance matrix \(P\) has to be estimated. Besides, the conditional probability to observe the \(N\) data \(\mathbf{J}^{\star}\) given the unknowns \(\mathbf{J}\) (known as the likelihood function) is also supposed to be a normal distribution of the form:  
\[\mathbf{J}^{\star}|\mathbf{J}\sim N(\mathbf{H}\mathbf{J},\mathbf{R}) \tag{7}\]  
The observation error is modeled as uniform and uncorrelated, resulting in a covariance matrix of the form \(\mathbf{R}=\epsilon^{2}\mathbf{I}\), where \(I\) is the unit matrix and \(\epsilon\) a pre-specified error of the observed variable values. The posterior distribution of the unknowns, conditional to the observed data, can be inferred by means of the Bayes theorem, such that:  
\[p(\mathbf{J}|\mathbf{J}^{\star})\sim p(\mathbf{J}^{\star}|\mathbf{J})p(\mathbf{J}) \tag{8}\]The posterior being the product of two normal distributions results in being a normal distribution as well, whose parameters depend on the expression of the prior covariance \(\mathbf{P}\):  
\[p(\mathbf{I}|\mathbf{J}^{\star})\sim N(E\left[\mathbf{I}|\mathbf{J}^{\star}\right],var\left[\mathbf{ J}|\mathbf{J}^{\star}\right]) \tag{9}\]  
Assuming that \(\mathbf{J}\) and \(\mathbf{J}^{\star}\) are normalized such that \(\mathbf{J}^{\star}\) has zero mean and unit variance, the covariance matrix \(\mathbf{P}\) of the prior distribution \(p(\mathbf{J})\) is defined such that its elements are generated by the squared exponential kernel (reported below in the 1D case, for simplicity):  
\[\mathbf{P}=cov(h_{ij},\theta)=\left[exp\left(-\frac{h_{ij}^{2}}{2\theta^{2}}\right)\right] \tag{10}\]  
where \(h_{ij}=\zeta_{i}-\zeta_{j}\) is the correlation range, \(\zeta_{i}\) the generic coordinate of the observation point \(\mathbf{J}^{\star}\), and \(\theta\) a hyperparameter that has to be estimated (for multi-dimensional problems, it is a vector whose dimension is the cardinality of the parameter space).