Abstract: Phase change material (PCM)-based latent heat storage (LHS) systems help address the mismatch between renewable energy supply and thermal demand. However, their practical implementation is constrained by the strongly nonlinear and multiphysics nature of phase change, which makes high-fidelity simulations and real-time applications computationally expensive. This review examines mathematical reduced-order modeling (ROM) as an effective strategy to overcome this limitation by combining physics-based simplifications, projection methods, interpolation techniques, and data-driven models for PCM-based LHS systems. While physical simplifications (such as dimensional reduction and effective property approximations) represent an important first layer of model reduction, the primary focus of this work is on the mathematical ROM methodologies that operate on the governing equations after such physical simplifications have been applied. The review covers approaches including two-temperature non-equilibrium and analytical thermal-resistance models, Proper Orthogonal Decomposition (POD), CFD-derived look-up tables, kriging and ε-NTU grey/black-box metamodels, and machine-learning methods such as artificial neural networks and gradient-boosted regressors trained from CFD data. These ROM techniques have been applied to packed beds, PCM-integrated heat exchangers, finned enclosures, triplex-tube systems, and solar thermal components, achieving speed-ups from tens to over 80,000 times faster than full CFD simulations while maintaining prediction errors typically below 5% or within sub-Kelvin temperature deviations. A critical comparative analysis exposes the fundamental trade-off between interpretability, data dependence, and computational efficiency, leading to a practical decision-making framework that guides method selection for specific applications such as design optimization, real-time control, and system-level simulation. Remaining challenges—including accurate representation of phase change nonlinearity, moving phase boundaries, multi-timescale dynamics, generalization across geometries, experimental validation, and integration into industrial workflows—motivate a structured roadmap for future hybrid physics–machine learning developments, standardized validation protocols, and pathways toward industrial deployment.