The DFT+U method is frequently employed to improve the first-principles description of strongly correlated materials. However, it is prone to deliver metastable electronic minima as in some cases different self-consistent electronic structure solutions are obtained depending on the initial occupation choice in the strongly correlated orbital space. While these local minima of the DFT+U method are often considered as computational artifacts, their physical meaning and relationship to true excited states remains unclear. We investigated the possibility of theoretically modeling transformations in the solid state which require thermal or optical excitations of electrons by taking into account the metastable states of the computationally undemanding DFT+U formalism on the example of the VO2 metal-insulator transition. Vanadium dioxide is well-suited as model system due to its distinct change in V d orbital occupations during its metal-insulator transition which can be used to manipulate the resulting electronic structure. We find the total DFT+U energy of VO2 to have a pronounced multi-minima character with strongly varying electronic structure characteristics and relative phase stabilities between the occurring local minima which can be sampled by choosing different initial V d orbital occupations. These local minima lie on different electronic potential energy surfaces and ionic relaxation results in different structures which can be assigned to VO2 high- and low-temperature phases, as well as other metastable phases, matching experimental literature. The identified metastable electronic states can hence be used to model the collapse of the VO2 band gap at elevated temperatures, upon photoexcitation, or during other phase transformations between semiconducting VO2 phases. Our results suggest that local DFT+U minima can indeed carry physical meaning while they remain underreported in theoretical literature on transition metal oxides like VO2. Tuning the U parameter to obtain multiple minima with the DFT+U approach can also potentially allow to model excitations in other suitable solid-state systems at much lower computational cost compared to e.g. time-dependent DFT, while providing decent agreement with experiment.