Drag reduction by laminarization, an innovative technology in aerodynamics design, provides the possibility to improve significantly the aerodynamic performance of aircraft and remains to be one of the most promising and effective technologies. In this paper, the design method for transonic natural laminar flow (NLF) wing at high Reynolds number is investigated in details and a new optimization model is established to solve conflicts between the increase of a laminar flow region and the shock wave strength at the trailing edge. In order to numerically handle this competitive situation efficiently, a two-objective hierarchical variable fidelity Pareto evolutionary optimization model based on search space contraction is proposed and implemented. Numerical experiments show that it can capture simultaneously a Pareto front of the two-objectives: a wave drag minimization and a laminar flow region maximization. The results also show that both wave drag and skin friction drag performances of non-dominated Pareto members can be significantly improved via the optimal laminar wing shape and shock wave control device. Nearly 40% of the laminar flow area on the wing surface is achieved at the Reynolds number of 107. It is concluded that the present hierarchical variable-fidelity multi-objective evolutionary optimization method is computationally efficient, and illustrate its potential for solving complex engineering design problems.