Symmetries in an open quantum system lead to degenerated Liouvillians that physically imply the existence of multiple steady states. In such cases, obtaining the initial condition independent steady states is highly nontrivial since any linear combination of the true asymptotic states, which may not necessarily be a density matrix, is also a valid asymptote for the Liouvillian. Thus, in this work, we consider different approaches to obtain the true steady states of a degenerated Liouvillian. In the ideal scenario, when the open system symmetry operators are known, we show how these can be used to obtain the invariant subspaces of the Liouvillian and hence the steady states. We then discuss two other approaches that do not require any knowledge of the symmetry operators. These could be powerful numerical tools to deal with quantum many-body complex open systems. The first approach that is based on Gram-Schmidt orthonormalization of density matrices allows us to obtain all the steady states, whereas the second one based on large deviations allows us to obtain the non-degenerated maximum and minimum current carrying states. We discuss the symmetry-decomposition and the orthonormalization methods with the help of an open para-benzene ring and examine interesting scenarios such as the dynamical restoration of Hamiltonian symmetries in the long-time limit and apply the method to study the eigenspacing statistics of the nonequilibrium steady state.