Atomic vacancies and vacancy aggregates control the thermodynamic stability and the functional response of graphene, yet the configurational space spanned by many vacancies at variable concentration and separation is too large to be mapped exhaustively by first-principles methods. Here, we map and rationalize this stability landscape by combining semiempirical atomistic thermodynamics, interpretable machine learning, and symbolic regression. Several hundred defective supercells, built from a 72-atom cell by varying the vacancy concentration and the inter-vacancy distance up to the fourth neighbor, were relaxed with the PM7 Hamiltonian in MOPAC, and the heat of formation was adopted as the stability metric. Each structure was encoded with the Dynamic Collision Fingerprint, a translationally and rotationally invariant descriptor that maps the local topology onto transport-like statistics of virtual probe particles. A gradient boosted decision tree model, optimized by Bayesian hyperparameter search, reproduces the heat of formation of an independent test set with a root mean squared error of approximately 23.5 kcal/mol and no evidence of overfitting, and a SHAP analysis identifies the defect concentration and the inter-vacancy distance as the two variables that dominate the stability. Symbolic regression then condenses the learned mapping into a compact closed-form expression that reproduces the heat of formation with a coefficient of determination of R² = 0.9966, a mean absolute error of 15.51 kcal/mol, and a root mean squared error of 24.18 kcal/mol. The workflow turns a high-dimensional structure–stability problem into an interpretable analytical law, providing a transferable route to rational defect engineering in two-dimensional materials.