We present a parameter set for obtaining the maximum number of atoms in a grating magneto-optical trap (gMOT) by employing a machine learning algorithm. In the multi-dimensional parameter space, which imposes a challenge for global optimization, the atom number is efficiently modeled via Bayesian optimization with the evaluation of the trap performance given by a Monte-Carlo simulation. Modeling gMOTs for six representative atomic species - 7Li, 23Na, 87Rb, 88Sr, 133Cs, 174Yb - allows us to discover that the optimal grating reflectivity is consistently higher than a simple estimation based on balanced optical molasses. Our algorithm also yields the optimal diffraction angle which is independent of the beam waist. The validity of the optimal parameter set for the case of 87Rb is experimentally verified using a set of grating chips with different reflectivities and diffraction angles.