Pillar stability is an important condition for safe work in room-and-pillar mines. The instability of pillars will lead to large-scale collapse hazards, and the accurate estimation of induced stresses at different positions in the pillar is helpful for pillar design and guaranteeing pillar stability. There are many modeling methods to design pillars and evaluate their stability, including empirical and numerical method. However, empirical methods are difficult to be applied to places other than the original environmental characteristics, and numerical methods often simplify the boundary conditions and material properties, which cannot guarantee the stability of the design. Currently, machine learning (ML) algorithms have been successfully applied to pillar stability assessment with higher accuracy. Thus, the study adopted a back-propagation neural network (BPNN) and five elements including the sparrow search algorithm (SSA), gray wolf optimizer (GWO), butterfly optimization algorithm (BOA), tunicate swarm algorithm (TSA), and multi-verse optimizer (MVO). Combining metaheuristic algorithms, five hybrid models were developed to predict the induced stress within the pillar. The weight and threshold of the BPNN model are optimized by metaheuristic algorithms, in which the mean absolute error (MAE) is utilized as the fitness function. A database containing 149 data samples was established, where the input variables were the angle of goafline (A), depth of the working coal seam (H), specific gravity (G), distance of the point from the center of the pillar (C), and distance of the point from goafline (D), and the output variable was the induced stress. Furthermore, the predictive performance of the proposed model is evaluated by five metrics, namely coefficient of determination (R2), root mean squared error (RMSE), variance accounted for (VAF), mean absolute error (MAE), and mean absolute percentage error (MAPE). The results showed that the five hybrid models developed have good prediction performance, especially the GWO-BPNN model performed the best (Training set: R2 = 0.9991, RMSE = 0.1535, VAF = 99.91, MAE = 0.0884, MAPE = 0.6107; Test set: R2 = 0.9983, RMSE = 0.1783, VAF = 99.83, MAE = 0.1230, MAPE = 0.9253).