This paper presents a novel method for computing three-dimensional magnetohydrodynamic (MHD) equilibrium states using artificial neural networks. Compared to equilibrium states computed using conventional solvers, we minimize the entire nonlinear global force residual in real space using a first-order optimization technique. We achieve the same minimum residual as those computed using conventional codes at a competitive computational cost. At higher computational cost, we achieve a lower residual minimum using neural networks, establishing a new lower bound on the force residual. We expect significant improvements in computing neural network models valid not only for single equilibrium states but also for continuous distributions of equilibrium states.