In this paper, a methodology is presented to compute the plasma properties (e.g.,, density and temperature) accounting for the dynamics of the excited states. The proposed strategy applies to both zero-dimensional (0D) models and multidimensional fluid and hybrid codes handling low-pressure (<50 mTorr) plasma discharges filled with argon, neon, krypton, and xenon gases. The paper focuses on two main aspects: (i) a lumping methodology is proposed to reduce the number of reactions and species considered in order to keep at bay the computational cost without a major loss of accuracy; (ii) the influence that different datasets of cross sections have on the results has been assessed. First, the lumping methodology has been implemented in a 0D model accounting for singly charged ions, neutrals, along with 1s and 2p excited states (Paschen notation). Metastable and resonant are treated as two separate species within the 1s energy level ($1s_M$ and $1s_R$ , respectively). The results have been benchmarked against those obtained treating each energy level of the excited states as an individual species. Differences lower than 1% have been obtained. Second, the results of the 0D model have been compared against measurements of electron density and temperature performed on an inductively coupled plasma. Numerical predictions and experiments present a disagreement up to 20%–30%, which is comparable to the uncertainty band of the measurements. Finally, the lumping strategy has been implemented in a 2D fluid code to assess its computational affordability, and the results have been compared against the experiments as well. A variance up to 30% in electron density and temperature is registered adopting different datasets of cross sections.