Abstract: An airfoil optimization method based on a MOGP (Muliple Output Gaussian Process) model and the simulated annealing algorithm is developed. The Latin hypercube method is used to construct a series of sample points in the design space, where the CST parameterization method is used to define the airfoil geometry. The initial MOGP model is established based on the response values obtained by CFD (Computational Fluid Dynamics) calculation. The proposed numerical test is to minimize the drag coefficient with consideration of some constrained conditions such as the area of airfoil, the lift coefficient and so on. Experimental results show that the drag coefficient is significantly reduced and the aera of airfoil and the lift coefficient remains, which is exactly meet the expectation of the optimization. Furthermore, the proposed method can be extended to 3D cases straightforwardly and is applicable to practical engineering problems.