In this work we develop, test, and implement a methodology that is able to perform, in an automated manner, "lumping" of a high-dimensional, discrete dynamical system onto a lower-dimensional space. Our aim is to develop an algorithm which, without any assumption about the nature of the system's slow dynamics, is able to reproduce accurately the long-time dynamics with minimal loss of information. Both the original and the lumped systems conform to master equations, related via the "lumping" analysis introduced by Wei and Kuo [Ind. Eng. Chem. Fundam. 8, 114 (1969)], and have the same limiting equilibrium probability distribution. The proposed method can be used in a variety of processes that can be modeled via a first order kinetic reaction scheme. Lumping affords great savings in the computational cost and reveals the characteristic times governing the slow dynamics of the system. Our goal is to approach the best lumping scheme with respect to three criteria, in order for the lumped system to be able to fully describe the long-time dynamics of the original system. The criteria used are: (a) the lumping error arising from the reduction process; (b) a measure of the magnitude of singular values associated with long-time evolution of the lumped system; and (c) the size of the lumped system. The search for the optimum lumping proceeds via Monte Carlo simulation based on the Wang-Landau scheme, which enables us to overcome entrapment in local minima in the above criteria and therefore increases the probability of encountering the global optimum. The developed algorithm is implemented to reproduce the long-time dynamics of a glassy binary Lennard-Jones mixture based on the idea of "inherent structures," where the rate constants for transitions between inherent structures have been evaluated via hazard plot analysis of a properly designed ensemble of molecular dynamics trajectories.