Engineering multi-objective optimization design often encounters problems with constraints. While, constraints handling has not been well-addressed for the multi-fidelity multi-objective Bayesian optimization method, especially for problems with high-dimensional design variables. To that end, a constrained parallel multi-fidelity multi-objective Bayesian optimization method assisted with a two-stage optimal search approach is developed. In the first stage, a constraint lower confidence bound infill criterion is proposed to locate a point aiming to improve the accuracy of the constraint boundary and find the potential optimal solution near the boundary. In the latter stage, the general expected improvement matrix (GEIM)-based criterion for constrained problem is posed to find multiple promising candidate points. Dedicated strategy is deployed to solve the massive auxiliary sub-optimization problem resulted from the GEIM-based infill criterion. The highly competitive performance of the proposed method is demonstrated by solving six analytic problems along with the comparison with the state-of-the-art methods. To solve the 3-point optimization design problem of a 3-stage axial flow compressor with 144 design variables efficiently, two improvements are incorporated to the proposed method: an efficient modeling method for Hierarchical Kriging to relief the curse of dimensionality of the high-dimension surrogate model and an imputation strategy to deal with simulation failures. Improvements of the 3 objectives with 4 stringent constraints are achieved with 2970 high-fidelity computational fluid dynamics simulations.