We present a novel framework, enabled by machine learning (ML) trainable models, for collisional plasma flow simulations governed by the Rosenbluth-Fokker-Planck (RFP) equation. The maximum entropy objective is applied to close the particle velocity distribution. The closure offers a flexible tool to reconstruct a large variety of distributions, using a few low-order statistical moments. The methodology is extended by introducing an artificial neural network (ANN) to couple the moments and the transport coefficients. In order to ensure the model's training and input data quality, we generate diverse training data sets using the high-fidelity Direct Simulation Monte-Carlo (DSMC) method. Moreover, energy constraints are integrated in the model training, thereby improving prediction accuracy. As a final step, the ANN is combined with a stochastic particle system, where the accuracy of the whole machinery for non-equilibrium thermal relaxation scenarios is evaluated. Simulation results indicate that our approach reproduces the correct relaxation behavior at lower computational cost compared to standard methods.