Refrigeration systems are complex, non-linear, multi-modal, and multi-dimensional. However, traditional methods are based on a trial and error process to optimize these systems, and a global optimum operating point cannot be guaranteed. Therefore, this work aims to study a two-stage vapor compression refrigeration system (VCRS) through a novel and robust hybrid multi-objective grey wolf optimizer (HMOGWO) algorithm. The system is modeled using response surface methods (RSM) to investigate the impacts of design variables on the set responses. Firstly, the interaction between the system components and their cycle behavior is analyzed by building four surrogate models using RSM. The model fit statistics indicate that they are statistically significant and agree with the design data. Three conflicting scenarios in bi-objective optimization are built focusing on the overall system following the Technique for Order of Preference by Similarity to Ideal Solution (TOPSIS) and Linear Programming Technique for Multidimensional Analysis of Preference (LINMAP) decision-making methods. The optimal solutions indicate that for the first to third scenarios, the exergetic efficiency (EE) and capital expenditure (CAPEX) are optimized by 33.4 and 7.5, and the EE and operational expenditure (OPEX) are improved by 27.4 and 19.0. The EE and global warming potential (GWP) are also optimized by 27.2 and 19.1, where the proposed HMOGWO outperforms the MOGWO and NSGA-II. Finally, the K-means clustering technique is applied for Pareto characterization. Based on the research outcomes, the combined RSM and HMOGWO techniques have proved an excellent solution to simulate and optimize two-stage VCRS. © 2023 Ahmed et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.