Determining the most stable configurations of nanoclusters, interfaces, and surfaces with a large number of constituent elements or adsorbate configurations is quite challenging due to the complexity of the potential energy surface. Additionally, reconstructing the interfaces due to the exposure to various physical conditions and the presence of different adsorbates drastically complicates the simulation processes. However, in reality, and in particular in heterogeneous catalysis, these phenomena can significantly govern the activation and/or deactivation of catalytic materials showcasing the importance of including these effects to predict the catalytic behavior more accurately. Although density functional theory (DFT) calculations can be used to estimate the stable structures at a very high computational expense, accurate prediction of most stable and meta-stable configurations requires sampling over a significant number of atomic orderings and shapes of the atomic structures. In this work, we developed a computational workflow to perform global optimizations. Specifically, a Python-based software package was developed integrating DFT, Genetic Algorithms (GA), and on-the-fly machine learning (ML). We show that incorporating the on-the-fly ML methods into the algorithms to estimate the electronic energies of the candidate structures significantly accelerated the search process reducing the dependency on DFT calculations.