The optimization of scenarios and design of real-time-control in tokamaks, especially for machines still in design phase, requires a comprehensive exploration of solutions to the Grad–Shafranov (GS) equation over a high-dimensional space of plasma and coil parameters. Emulators can bypass the numerical issues in the GS equation, if a large enough library of equilibria is available. We train an ensemble of neural networks to emulate the typical shape-control targets (separatrix at midplane, X-points, divertor strike point, flux expansion, and poloidal beta) as a function of plasma parameters and active coil currents for the range of plasma configurations relevant to spherical tokamaks with a super-X divertor, with percent-level accuracy. This allows a quick calculation of the classical-control shape matrices, potentially allowing real-time calculation at any point in a shot with submillisecond latency. We devise a hyperparameter sampler to select the optimal network architectures and quantify uncertainties on the model predictions. To generate the relevant training set, we devise a Markov-chain Monte Carlo algorithm to produce large libraries of forward Grad–Shafranov solutions without the need for user intervention. The algorithm promotes equilibria with desirable properties, while avoiding parameter combinations resulting in problematic profiles or numerical issues in the integration of the GS equation.