The development of artificial intelligence (AI) greatly boosts scientific and engineering innovation. As one of the promising candidates for transiting the carbon intensive economy to zero emission future, proton exchange membrane (PEM) fuel cells has aroused extensive attentions. The gas diffusion layer (GDL) strongly affects the water and heat management during PEM fuel cells operation, therefore multi-variable optimization, including thickness, porosity, conductivity, channel/rib widths and compression ratio, is essential for the improved cell performance. However, traditional experiment-based optimization is time consuming and economically unaffordable. To break down the obstacles to rapidly optimize GDLs, physics-based simulation and machine-learning-based surrogate modelling are integrated to build a sophisticated M5 model, in which multi-physics and multi-phase flow simulation, machine-learning-based surrogate modelling, multi-variable and multi-objects optimization are included. Two machine learning methodologies, namely response surface methodology (RSM) and artificial neural network (ANN) are compared. The M5 model is proved to be effective and efficient for GDL optimization. After optimization, the current density and standard deviation of oxygen distribution at 0.4 V are improved by 20.8% and 74.6%, respectively. Pareto front is obtained to trade off the cell performance and homogeneity of oxygen distribution, e.g., 20.5% higher current density is achieved when sacrificing the standard deviation of oxygen distribution by 26.0%.