Reduced-order models (ROMs) are today an indispensable element of contemporary science and engineering. ROMs surpass the capabilities of traditional high-fidelity simulations by striking a unique, optimal balance between cost efficiency and accuracy, opening the door to fast and effective numerical analyses. This transforms the applicability of modeling with regard to both scientific discovery and engineering developments. The paradigm of physics-informed/physics-inspired machine learning (ML) serves to further transcend the ROMs in terms of their interpretability and generalizability. Predictive, interpretable, and generalizable ML-enabled ROMs facilitate scientific research into yet-unresolved complex phenomena. On the applied side, one of the strong motivators for such ROMs is the transformative technology of digital twins. The Electric Propulsion (EP) sector and the broader field of plasma physics can benefit significantly from advances in ML-enabled physics-inspired ROMs that feature low-dimensionality and sparsity. This paper aims to cast light on reduced-order modeling within the context of plasma systems relevant to EP. Following on introductory discussions on the practical relevance of ROMs, particularly in the framework of envisioned EP digital twins, we describe the two central elements of ROM development, i.e., finding a proper reduced coordinate system, and describing the system’s dynamics within that reduced coordinate. We emphasize how ML and data-driven (DD) approaches can majorly aid each of these elements. Taking the example of a plasma discharge that represents a radial-azimuthal Hall thruster configuration, we demonstrate how autoencoder networks can find a highly low-dimensional space that fully represents the high-dimensional plasma state. We discuss multiple possibilities for ML/DD learning of the time dynamics in the low-dimensional latent space of the autoencoder. We finally present dynamics forecasting results from a fast and predictive ROM obtained using the optimized dynamic mode decomposition method, which amounts to a special case of combining a linear shallow autoencoder with a linear operator for describing the time dynamics.