Rigorous black-box simulations are useful to describe complex systems. However, it cannot be directly integrated into mathematical programming models in some algebraic modeling environments because of the lack of symbolic formulation. In the present paper, a framework is proposed to embed the Aspen HYSYS process simulator in GAMS using kriging surrogate models to replace the simulator-dependent, black-box objective, and constraints functions. The approach is applied to the energy-efficient C3MR natural gas liquefaction process simulation optimization using multi-start nonlinear programming and the local solver CONOPT in GAMS. Results were compared with two other meta-heuristic approaches, Particle Swarm Optimization (PSO) and Genetic Algorithm (GA), and with the literature. In a small simulation evaluation budget of 20 times the number of decision variables, the proposed optimization approach resulted in 0.2538 kW of compression work per kg of natural gas and surpassed those of the PSO and GA and the previous literature from 2.45 to 15.3 %.