Among the pillars of decarbonization of the global energy system, energy efficiency plays a key role in reducing energy consumption across end-use (industry, transport and buildings) sectors. In industrial processes, energy efficiency can be improved by exploiting heat recovery via heat exchange between process streams. This paper develops a stage-wise superstructure-based mathematical programming model for the optimization of heat exchanger networks. The model incorporates rigorous formulation to handle process streams with phase change (condensation or evaporation), and is applied to a case study of an ethylene glycol production plant in Taiwan for minimizing utility consumption. The results show a compromise between steam savings and process feasibility, as well as how the model is modified to reflect practical considerations. In the preliminary analysis, with a substantial potential steam saving of 15,476 kW (28%), the solution involves forbidden matches that pose a hazard to the process and cannot be implemented. In the further analysis without process streams that cause forbidden matches, although the space limitation in the plant renders the best solution infeasible, the compromise solution can achieve a considerable steam saving of up to 8448 kW (91%) and is being evaluated by the plant managers and operators.