This paper presents an interpretable machine-learning framework integrating RK4IP simulation, surrogate modeling, SHAP analysis, and genetic algorithm (GA) optimization to address the challenge of optimizing nonlinear optical loop mirror (NOLM) based passively mode-locked fiber lasers in their high-dimensional coupled parameter space. The Ginzburg-Landau equation is solved via the fourth-order Runge-Kutta in the interaction picture (RK4IP) method to simulate pulse evolution in the laser cavity. A comprehensive dataset of 500 samples is generated by uniformly varying six key parameters: second-order dispersion (β₂, 11.2-29.2 ps²/km), nonlinear refractive index (n₂, 2.3-2.9×10^-16 cm²/W), small-signal gain (g_ss, 36-39 dB), gain saturation power (P_sat, 22-25 dBm), gain fiber length (L_AF, 0.7-1.3 m), and gain bandwidth (Δλ, 50~110 nm). Four surrogate models—Random Forest (RF), XGBoost (XGB), backpropagation neural network (BPNN), and one-dimensional convolutional neural network (1D CNN) — are systematically compared. BPNN and 1D CNN demonstrate superior generalization, achieving test-set R² values of 0.951 for peak power, 0.981 for pulse energy, and 0.996 for pulse duration, with ratio of performance to deviation (RPD) reaching 4.56, 7.31, and 15.94, respectively. RF and XGB exhibit certain overfitting tendencies in specific tasks. SHAP analysis quantifies feature contributions and uncovers the physical mechanisms governing output characteristics. For peak power, n₂ and g_ss are the most influential parameters, with SHAP dependence plots revealing that n₂ exerts a negative effect while g_ss shows a strong positive correlation. For pulse energy and pulse duration, g_ss, P_sat, 22-25 dBm, and L_AF collectively dominate model predictions, while β₂, n₂, and Δλ play marginal roles. Based on these interpretability results, a two-stage GA optimization is implemented. Baseline optima are first established through full six-dimensional search for peak power, pulse energy, and pulse duration. Subsequently, parameters with consistently low SHAP importance across models are fixed at nominal values, and GA is re-executed in the reduced subspace. The optimized outputs from the reduced search closely match full-space results: RF and BPNN retain over 99% of baseline peak power performance; BPNN and 1D CNN preserve 86.5% of pulse energy; pulse duration retention exceeds 94% for RF, BPNN, and 1D CNN. Meanwhile, computational time is reduced by 60%-70% and convergence accelerates by 30%-60%, demonstrating substantial efficiency gains without sacrificing solution quality. This work establishes an efficient, interpretable framework for NOLM laser optimization, combining simulation, surrogate modeling, SHAP interpretation, and GA search to enable accurate prediction, global optimization, and physical insight.