Abstract: GH4065A, a wrought nickel-based superalloy commonly used for aero-engine turbine disks, is susceptible to fatigue crack initiation under elevated-temperature, centrifugal, and cyclic loading conditions. Shot peening can improve fatigue performance by introducing residual compressive stress and modifying the near-surface microstructure. A random multi-shot peening model was developed using ABAQUS and Python scripting, and a simulation scheme based on the Box- Behnken response surface methodology was employed to systematically investigate the effects of shot material (ceramic, cast steel, and glass), shot diameter, impact velocity, and coverage on surface roughness, surface residual compressive stress, and the overall residual compressive stress layer depth. Correlation analysis and analysis of variance were conducted to examine the effects and interactions among the process parameters. An XGBoost surrogate models were constructed, and its prediction performance was evaluated using ten repetitions of five-fold grouped cross-validation. Analysis of feature importance was performed, and multi objective Pareto optimization was further implemented. The results show that, within the investigated parameter range, shot diameter and impact velocity are the primary factors affecting surface roughness and residual compressive stress layer depth, whereas coverage has the greatest influence on surface residual compressive stress. Ceramic shot exhibits favorable residual compressive stress uniformity and overall strengthening performance. A total of 313 Pareto-nondominated candidations are identified. The compromise candidate employs ceramic shot with a diameter of 0.42 mm, an impact velocity of 80 m/s, and a coverage of 200%. The predicted surface roughness, surface residual compressive stress, and overall residual compressive stress layer depth are 22.351 μm, -981.457 MPa, and 0.3309 mm, respectively. The relative errors between the independent finite element back-analysis results and the surrogate-model predictions are 10.71%, 2.38%, and 4.38%, respectively. Moreover, the absolute errors of all three responses are lower than the corresponding root mean square errors(RMSE) obtained from repeated grouped cross-validation, which demonstrates reasonable agreement between the surrogate-model predictions and the finite element results. These findings provide a reference for optimizing shot peening process parameters for GH4065A superalloy.