Multi-objective optimization problems of electromagnetic devices generally involve finding the optimal solution set by iteratively computing objective functions corresponding to design variables, which incurs high computational costs. This issue is particularly pronounced when the optimization process relies on high-fidelity finite element analysis (FEA). To address this issue, a Kriging surrogate model-assisted multi-objective optimization method is proposed in this paper. In order to balance model training time with the exploration of information in the design space, representative and promising samples are selected from the non-dominated solution set and the dominated solution set to form the model’s training dataset. A two-stage hybrid infill strategy is developed to enhance the model’s capability to explore true non-dominated solutions and improve its global prediction accuracy in the meantime. The effectiveness of the proposed method is validated by performance test of test functions. Besides, the proposed method is applied to the optimal design of the TEAM Problem 22-Superconducting Magnetic Energy Storage (SMES) device. Numerical results demonstrate that the proposed method can significantly reduce the simulation number of calling FEA, thereby speeding up the optimal design process.