In this paper, we address aerodynamic shape optimization problems including uncertain operating conditions. After a review of the possible approaches to take into account uncertainty, we propose to use meta-modeling techniques in order to develop a two-level modeling procedure for statistics estimation. Radial basis functions are employed to approximate the aerodynamic coefficients as operating conditions vary. Then, a Monte-Carlo method is employed to estimate statistics using the approximate model. The proposed approach is applied to the robust optimization of the wing shape of a business aircraft, by minimizing the mean and the variance of the drag coefficient with uncertain free-stream Mach number.