<p>Solving engineering design problems involving uncertainties is often intractable for high-fidelity models due to the computational cost of sampling expensive simulations to estimate output statistics. We propose a surrogate-based method for efficiently solving constrained, robust design optimization problems involving integrated statistical moments of high-fidelity simulation outputs. Similar to multi-fidelity deterministic optimization techniques, the method progressively increases the accuracy of the statistical moments by increasing sampling rate based on the optimizer’s proximity to convergence. Both the optimization and statistical integration leverage a partition-of-unity, gradient-enhanced surrogate model constructed over both the design variable and uncertain parameter spaces. We also develop an error estimate for optimality and a corresponding adaptive sampling algorithm, which is used at the end of each design cycle to improve the accuracy of the surrogate model and, consequently, the statistical moments. We apply the method to various analytical benchmark optimization-under-uncertainty problems, demonstrating increased efficiency compared to surrogate-free, single-fidelity approaches to reach local minimizers. In addition, we apply the method to a robust optimization involving shock–boundary layer interaction for which each sample requires the solution of a fully-coupled aerostructural analysis.</p>