Classical integer-order models provide useful insights into alcoholism dynamics but cannot adequately capture the memory and hereditary effects associated with addiction, relapse, and recovery. In this work, the Caputo fractional derivative is used to study an alcoholism model describing susceptible to alcohol consumption ( \({\textbf {P}}\) ), moderate drinkers ( \({\textbf {M}}\) ), heavy drinkers ( \({\textbf {H}}\) ), affluent heavy drinkers receiving care in private treatment facilities ( \({\textbf {T}}^{\textsf {r}}\) ), disadvantaged heavy drinkers attending public treatment centers ( \({\textbf {T}}^{\textsf {p}}\) ), and those who have stopped drinking ( \({\textbf {Q}}\) ). Using the concepts of the fixed point theory and nonlinear analysis, the existence, uniqueness, and stability of solutions of the considered fractional alcoholism model are rigorously proved. Using the Newton interpolation polynomial, the numerical analysis and simulations of the fractional alcoholism model are carried out to visualize the behavior of the considered system and to examine the effects of various crucial parameters on its dynamics. Optimal control analysis for the fractional order model is demonstrated to show its dynamics with and without control. A deep neural network surrogate is also developed to approximate the numerical trajectories generated by the fractional-order model.