Multi-Objective Optimization of a Permanent Magnet Actuator for High Voltage Vacuum Circuit Breaker Based on Adaptive Surrogate Modeling Technique

A novel mono-stable permanent magnet actuator (PMA) for high voltage vacuum circuit breaker (VCB) and its optimal design method are proposed in this paper. The proposed PMA is featured with a structure of separated magnetic circuits, which makes the holding part and closing driving part work independently without interference. The application of an auxiliary breaking coil decreases the response time in the initial phase of opening operation, and an external disc spring is adopted to accelerate the opening movement, which makes the PMA meet the fast-breaking requirement of high voltage VCB. As calculating the characteristics of the PMA accurately through numerical simulation is a time-consuming process, a multi-objective optimization (MOO) algorithm based on surrogate modeling technique and adaptive samples adding strategy are proposed to reduce the workload of numerical simulations during optimization. Firstly, initial surrogate models are constructed and evaluated, and then iteratively updated to improve their global approximating abilities. Secondly, according to the approximate MOO results obtained by the global surrogate models, additional samples are added to constantly update the surrogate models to gradually improve the models’ local accuracies in optimal solution regions and finally guide the algorithm to the true Pareto front. The efficiency and accuracy of the proposed algorithm are verified by test functions. By applying the optimization strategy to the design of the proposed PMA, a set of satisfying Pareto optimal solutions, which improve the overall performance of the PMA obviously, can be derived at a reasonable computation cost.


Introduction
Vacuum circuit breakers (VCBs) have been widely used in the power system of middle/low voltage levels (3.6~40.5 kV) and are now being developed to high voltage levels (72.5~252 kV) due to their advantages of being environmentally friendly, maintenance-free, and having excellent insulating and breaking performances [1][2][3]. As a key component of a VCB, the operating mechanism drives a movable contact in the vacuum interrupter to switch on or off the circuit through an insulating pole. The permanent magnet actuator (PMA) is a type of operating mechanism and has received much attention because of its advantages such as high reliability, high controllability, and excellent consistency in operating time, which make it extensively used in medium-voltage VCBs [4,5]. However, the VCB applied in the high voltage field has a much longer contact stroke and a higher velocity requirement for operating mechanism, which limits the application of PMAs in high voltage VCBs. In order to extend the PMAs' application to high voltage power systems, much research work has been method and trust region updating strategy. Yin et al. [34] presented a novel metamodel-based MOO algorithm combining adaptive radial basis functions with the NSGA-II algorithm and verified its superiority over the traditional metamodel-based MOO method. However, the above mentioned two methods have not carried out the global accuracy evaluation or global accuracy improvement for the constructed initial surrogate models, which makes the subsequent optimizations easily fall into local optimal solutions. This paper proposes a MOO method, which takes into account both the global accuracy of surrogate models and the local accuracy near the Pareto frontier.
This paper focuses on the design and MOO of a novel PMA for high voltage VCBs. Section 2 introduces the topology structure and working principle of PMA. The PMA with a mono-stable and axisymmetric structure separates the holding part from the closing driving part to eliminate the interference between their magnetic fields. The combination of breaking spring and auxiliary breaking coil can make it possible to rapidly decrease the holding force in the initial phase of the breaking operation, and then realize the fast-breaking movement driven by the breaking spring. Section 3 describes the proposed algorithm and its properties and conducts numerical tests. In order to reduce the optimization time, a MOO algorithm based on adaptive surrogate modeling techniques is proposed. In the first stage, the algorithm emphasizes enhancing the global accuracies of surrogate models through a samples adding strategy, thus as to make the surrogate models able to describe the overall features of real systems as accurately as possible. In the second stage, the algorithm pays attention to adding samples around the Pareto front to improve the local accuracies of surrogate models in optimal solution regions and finally guides the searching process into fast convergence. The proposed algorithm was tested by mathematical functions to verify its accuracy and efficiency. Section 4 applies the algorithm to the optimization of the PMA to validate its practicability. Section 5 gives conclusions and future work.

Topology Structure
A single-phase VCB for 126 kV high voltage level is sketched in Figure 1. It is a direct drive type of VCB, and its vacuum interrupter and operating mechanism are connected by the insulating pole and contact spring. The vacuum interrupter plays a role in extinguishing the arc and breaking the circuit, while the operating mechanism functions as a driving source for the movable contact.
Energies 2019, 12, x FOR PEER REVIEW 3 of 20 much better Pareto frontier in case of a limited number of function evaluations [32]. Chen et al. [33] proposed a sequential approximate MOO method that employs the adaptive surrogate modeling method and trust region updating strategy. Yin et al. [34] presented a novel metamodel-based MOO algorithm combining adaptive radial basis functions with the NSGA-II algorithm and verified its superiority over the traditional metamodel-based MOO method. However, the above mentioned two methods have not carried out the global accuracy evaluation or global accuracy improvement for the constructed initial surrogate models, which makes the subsequent optimizations easily fall into local optimal solutions. This paper proposes a MOO method, which takes into account both the global accuracy of surrogate models and the local accuracy near the Pareto frontier. This paper focuses on the design and MOO of a novel PMA for high voltage VCBs. Section 2 introduces the topology structure and working principle of PMA. The PMA with a mono-stable and axisymmetric structure separates the holding part from the closing driving part to eliminate the interference between their magnetic fields. The combination of breaking spring and auxiliary breaking coil can make it possible to rapidly decrease the holding force in the initial phase of the breaking operation, and then realize the fast-breaking movement driven by the breaking spring. Section 3 describes the proposed algorithm and its properties and conducts numerical tests. In order to reduce the optimization time, a MOO algorithm based on adaptive surrogate modeling techniques is proposed. In the first stage, the algorithm emphasizes enhancing the global accuracies of surrogate models through a samples adding strategy, thus as to make the surrogate models able to describe the overall features of real systems as accurately as possible. In the second stage, the algorithm pays attention to adding samples around the Pareto front to improve the local accuracies of surrogate models in optimal solution regions and finally guides the searching process into fast convergence. The proposed algorithm was tested by mathematical functions to verify its accuracy and efficiency. Section 4 applies the algorithm to the optimization of the PMA to validate its practicability. Section 5 gives conclusions and future work.

Topology Structure
A single-phase VCB for 126 kV high voltage level is sketched in Figure 1. It is a direct drive type of VCB, and its vacuum interrupter and operating mechanism are connected by the insulating pole and contact spring. The vacuum interrupter plays a role in extinguishing the arc and breaking the circuit, while the operating mechanism functions as a driving source for the movable contact. A novel PMA is proposed, as shown in Figure 2. As a type of operating mechanism, the PMA consists of a driving part in the lower half and a holding part in the upper half, which are A novel PMA is proposed, as shown in Figure 2. As a type of operating mechanism, the PMA consists of a driving part in the lower half and a holding part in the upper half, which are encapsulated in a cylindrical non-magnetic shell and separated by a non-magnetic spacer. The driving part includes a static core and a movable core for driving, a closing coil, and a disc spring that is installed between the shell and the movable core for holding. The holding part includes a PM and a static iron core for holding, a movable core for holding, and an auxiliary breaking coil. The two movable cores are fixed on a shaft, through which the electromagnetic force can be transferred to the movable contact of VCB.
Energies 2019, 12, x FOR PEER REVIEW 4 of 20 encapsulated in a cylindrical non-magnetic shell and separated by a non-magnetic spacer. The driving part includes a static core and a movable core for driving, a closing coil, and a disc spring that is installed between the shell and the movable core for holding. The holding part includes a PM and a static iron core for holding, a movable core for holding, and an auxiliary breaking coil . The two  movable cores are fixed on a shaft, through which the electromagnetic force can be transferred to the  movable contact of VCB.  The main parameters of a 126 kV vacuum interrupter and the requirements for the actuator are  listed in Table 1.

Working Principle
As shown in Figure 3a, the studied PMA was assumed to be initially maintained in the open state by the compression force of the spring exerted on the movable core for holding. Once a closing command comes, the closing coil is energized, and there arises a magnetic field around the coil accompanied by the generation of current in the coil. Once the current in the coil increases to a certain extent, the electromagnetic force exerted on the movable core in driving part exceeds the anti-force, and the movable parts begin to move towards the closed position, as shown in Figure 3b. Ultimately, the closing coil is de-energized after the movable parts arrive at the closed position, and the PMA is held in the closed state under the action of enough force generated by the PM, as shown in Figure 3c.
When it is time to break, the auxiliary breaking coil is excited, and the magnetic field induced by the current in the coil counteracts the PM field, which leads to a rapid decrease of holding force, as shown in Figure 3d. Once the current in the coil increases to a certain value, the movable core for holding starts to move down driven by the breaking spring. Finally, the auxiliary breaking coil is deenergized after the moving parts arrive at the open position and the PMA is held in the open state.
The auxiliary breaking coil helps reduce the response time in the initial phase of the breaking operation by rapidly decreasing the holding force. With the disc spring as the driving source in the breaking operation, it is able to satisfy the requirement for high breaking velocity. To sum up, the proposed PMA exhibits good performance in closing and breaking operations.  The main parameters of a 126 kV vacuum interrupter and the requirements for the actuator are listed in Table 1.

Working Principle
As shown in Figure 3a, the studied PMA was assumed to be initially maintained in the open state by the compression force of the spring exerted on the movable core for holding. Once a closing command comes, the closing coil is energized, and there arises a magnetic field around the coil accompanied by the generation of current in the coil. Once the current in the coil increases to a certain extent, the electromagnetic force exerted on the movable core in driving part exceeds the anti-force, and the movable parts begin to move towards the closed position, as shown in Figure 3b. Ultimately, the closing coil is de-energized after the movable parts arrive at the closed position, and the PMA is held in the closed state under the action of enough force generated by the PM, as shown in Figure 3c.
When it is time to break, the auxiliary breaking coil is excited, and the magnetic field induced by the current in the coil counteracts the PM field, which leads to a rapid decrease of holding force, as shown in Figure 3d. Once the current in the coil increases to a certain value, the movable core for holding starts to move down driven by the breaking spring. Finally, the auxiliary breaking coil is de-energized after the moving parts arrive at the open position and the PMA is held in the open state.
The auxiliary breaking coil helps reduce the response time in the initial phase of the breaking operation by rapidly decreasing the holding force. With the disc spring as the driving source in the breaking operation, it is able to satisfy the requirement for high breaking velocity. To sum up, the proposed PMA exhibits good performance in closing and breaking operations.

Analysis Method
The characteristics of the PMA were obtained by solving a multi-physics coupling problem, including electric circuit, magnetic field, and mechanical motion. The magnetic field can be regarded as a magneto-static field by ignoring eddy currents. The magneto-static field equation is described as: where μ is the permeability of the material, A is the magnetic vector potential, Jc is the current density vector.
In addition to the magnetic field equations, the PMA also follows the voltage balance equation and D'Alembert Equation. The state equations for the PMA are described as , , where ψ is the flux linkage across the coil; Uc and C represent the voltage and capacitance of the capacitor, respectively; i and R represent the coil current and coil resistance, respectively; Fmag is the magnetic force exerted on the moving part; Ff is the reaction force; m, v, and x denote the mass, velocity, and displacement of the moving part, respectively. When the PMA is held in the open state or closed state as shown in Figure 3a or Figure 3c, the static characteristics of the PMA involving PM holding force and distribution of magnetic field can be achieved by solving the magnetic field Equation (1) using the finite-element method. When the coils in Figure 2 are energized, the magnetic field begins to change. The dynamic characteristics involving coil current, electromagnetic force exerted on the movable cores, displacement, and velocity of the movable cores in Figure 3b,d can be acquired by solving the state Equations (2) using the time-differential method.

General Overview of the Algorithm
In order to reduce the amount of computation while ensuring the accuracy of the solution, a new method based on the adaptive surrogate modeling technique was utilized to effectively solve the MOO problems. The flowchart of the proposed algorithm is shown in Figure 4. The method divided

Analysis Method
The characteristics of the PMA were obtained by solving a multi-physics coupling problem, including electric circuit, magnetic field, and mechanical motion. The magnetic field can be regarded as a magneto-static field by ignoring eddy currents. The magneto-static field equation is described as: where µ is the permeability of the material, A is the magnetic vector potential, J c is the current density vector. In addition to the magnetic field equations, the PMA also follows the voltage balance equation and D'Alembert Equation. The state equations for the PMA are described as where ψ is the flux linkage across the coil; U c and C represent the voltage and capacitance of the capacitor, respectively; i and R represent the coil current and coil resistance, respectively; F mag is the magnetic force exerted on the moving part; F f is the reaction force; m, v, and x denote the mass, velocity, and displacement of the moving part, respectively. When the PMA is held in the open state or closed state as shown in Figure 3a or Figure 3c, the static characteristics of the PMA involving PM holding force and distribution of magnetic field can be achieved by solving the magnetic field Equation (1) using the finite-element method. When the coils in Figure 2 are energized, the magnetic field begins to change. The dynamic characteristics involving coil current, electromagnetic force exerted on the movable cores, displacement, and velocity of the movable cores in Figure 3b,d can be acquired by solving the state Equations (2) using the time-differential method.

General Overview of the Algorithm
In order to reduce the amount of computation while ensuring the accuracy of the solution, a new method based on the adaptive surrogate modeling technique was utilized to effectively solve the MOO problems. The flowchart of the proposed algorithm is shown in Figure 4. The method divided the optimization of complex problems into 2 stages: The first stage focused on the global accuracy improvement of surrogate models, based on which the optimization algorithm was able to acquire a rough Pareto front. This front may not be precise enough, but it contained all the potential optimal solutions. The second stage placed additional samples on the regions around the rough Pareto front in order to gradually improve the accuracy of the front, and finally obtained the desired Pareto optimal solutions set. The detailed steps are described as follows:

Design of Experiments
The design of experiment (DOE) method was used to generate the initial samples distribution in the design space. Based on the idea of space-filling, the modern DOE methods generated limited samples spread out over the entire design space as evenly as possible. In this study, the translational propagation Latin hypercube design (TPLHD) [35] was used to obtain optimal or near-optimal LHD designs with minimal computational effort.

Construction of Surrogated Models
The RBF model was a popular surrogate modeling technique and had been successfully developed for approximating the non-linear and complex systems in practical engineering [36]. It was constructed by a linear combination of basis functions, each of which was a function of the Euclidean distance between unknown points and existed sample points. With a set of scattered sample points xi (i = 1, 2, …, n) the mathematical form of RBF model can be given by: are the weight coefficients, n is the number of samples. As an interpolation method, the RBF model should satisfy: (5) where yi is the real response value of xi. It can be expressed in a matrix form as:  Step 1: Define the objectives and constraints of the problem, initialize the number and bounds of variables.
• STAGE I Step 2: Generate initial samples by carrying out the design of experiments (DOE).
Step 3: Calculate the objective and constraint function values of the initial samples. In practical engineering problems, the time-consuming numerical simulation is used to evaluate the samples.
Step 4: Adopt the radial basis function (RBF) model to construct the approximate models of objectives and constraints with the input and output data of the initial samples. The parameters of the RBF models are also tuned to reach the optimal fitting level.
Step 5: Assess the accuracies of the surrogate models obtained in Step 4. If the accuracies satisfy the preset requirement, the obtained optimal surrogate models are saved, and the process moves on to STAGE II. Otherwise, the process goes to Step 6.
Step 6: Use an adaptive sampling strategy to select several extra samples according to the global exploration and local exploitation criteria. After obtaining the corresponding objective and constraint values of extra samples by numerical simulations, the process goes to Step 4 to update the surrogate models. •

STAGE II
Step 7: Implement the MOO based on the surrogate models to obtain the approximate Pareto front set.
Step 8: Select several additional samples, which are uniformly distributed along the approximate Pareto front, and calculate the objective and constraint function values using numerical simulations.
Step 9: Check convergence by calculating the following two error criteria: where f i denotes the real response value calculated by numerical simulation, ∼ f i denotes the approximate response value calculated by surrogate models, − f i denotes the mean value of real response values on samples of the validation set, n denotes the number of validation samples. R 2 assesses the model accuracy on the Pareto front from the perspective of global characteristics, while relative maximum absolute error (RMAE) from the perspective of local characteristics. The closer R 2 is to 1 and the closer RMAE to 0 means the higher model accuracy on the Pareto front and the better approximation of the real Pareto front.
The entire optimization process will terminate if one of the following two conditions is satisfied: (a) The error criteria R 2 of all the response models (objective and constraint models) on the approximate Pareto front are greater than 0.9, and the RMAEs of all response models are less than 0.1. (b) The total number of numerical simulations on samples is greater than a given value.
Step 10: If the termination condition is not satisfied, add the evaluated new samples to the sample set and update the surrogate models of the objectives and constraints, then go to Step 7. Otherwise, rank all the solutions, which are selected from the appropriate Pareto front in every iteration, using a non-dominated sorting method and get a set of non-dominated solutions as the ultimate optimal solutions set.

Design of Experiments
The design of experiment (DOE) method was used to generate the initial samples distribution in the design space. Based on the idea of space-filling, the modern DOE methods generated limited samples spread out over the entire design space as evenly as possible. In this study, the translational propagation Latin hypercube design (TPLHD) [35] was used to obtain optimal or near-optimal LHD designs with minimal computational effort.

Construction of Surrogated Models
The RBF model was a popular surrogate modeling technique and had been successfully developed for approximating the non-linear and complex systems in practical engineering [36]. It was constructed by a linear combination of basis functions, each of which was a function of the Euclidean distance between unknown points and existed sample points. With a set of scattered sample points x i (i = 1, 2, . . . , n) the mathematical form of RBF model can be given by: . . , n) are the weight coefficients, n is the number of samples. As an interpolation method, the RBF model should satisfy: where y i is the real response value of x i . It can be expressed in a matrix form as: where and β = (β 1 , β 2 , . . . , β n ) T is a vector of weight coefficients, Υ = (y 1 , y 2 , . . . , y n ) T is a vector of real response values. Therefore, β can be determined by: The common basis functions are listed in Table 2, where r = x − x i , and c is the model parameter. In this study, the genetic algorithm was used to find the optimal model parameters, in which case the RBF surrogate models possessed the best accuracies. Table 2. Commonly used basis functions.

Basic Function Formula Basic Function Formula
Linear

Model Evaluation Indicator and Modeling Termination Judgment
Once the construction process was finished, the obtained surrogate models have to be assessed for their accuracies, thus as to ensure their effectiveness for subsequent applications. The cross-validation was a simple and classic method for creating an efficient estimation of model accuracy. It divided the sample set into 2 parts: One for the construction of the surrogate model, the other for the accurate evaluation of the constructed model. In this section, we chose the 5-fold cross-validation, and the root mean square error (RMSE) to calculate the model accuracy. The RMSE is defined as: where f is the real objective response, ∼ f is the surrogate model, x i is the sample point in the test set. If the surrogate model accuracy satisfies the predetermined requirement, the process of surrogate modeling will be terminated. Actually, the accuracies of the surrogate models for some complex real responses will not reach the preset values with a relatively small number of real function calls. Given the expensive computation cost, the modeling process will also be stopped if the total number of real function calls is greater than a specified value.

Sequential Sampling Strategy
The traditional surrogate modeling technique was based on data points sampled all at once, which lead to a potential waste of simulations. The sequential sampling strategy gets samples step by step based on previous samples and modeling results. It uses the information from the last simulations to select new samples as optimally as possible.
In order to obtain accurate global surrogate models, we proposed a sequential sampling strategy combining 2 metrics: One for exploration and another for exploitation. The Monte Carlo Voronoi approximation was chosen to explore new regions because of its simplicity and computational efficiency [37]. The cross-validation approach was used to estimate the local error behavior for exploitation. This method helped realize the trade-off between global exploration and local exploitation and made it possible to achieve more accurate global surrogate models with fewer samples.
The Voronoi diagram was an efficient way to identify the distribution density of samples in the design space, depending on which we can place additional samples in the region with low density to improve the model accuracy. Assume there already exists a set of sample points P = {p 1 , p 2 , . . . , p m } ∈ R d , the Voronoi diagram method partitions the whole design space into a set of Voronoi cells C = {C 1 , Energies 2019, 12, 4695 9 of 19 C 2 , . . . , C m }. C i denotes the region surrounding the point p i , in which all points are closer to p i than any other points in P. Figure 5 shows the partition result of a 2D design space.  (9) where dom(pi, pj) represents the closed half plane bounded by the perpendicular bisector of pi and pj.
The bisector separates the points closer to pi from those closer to pj.
Once the partition was accomplished, it was possible to obtain the density distribution of samples according to the volume of Voronoi cells. Actually, it was difficult to exactly calculate the Voronoi diagram of a given set of samples and design space. The Monte Carlo approach was adopted to realize the approximate Voronoi partition and estimate the volume of Voronoi cells. Firstly, a large number of randomly and uniformly distributed points were generated in the design space. Then for each of these random points, the distances between it and all the existing sample points were calculated. Each random point will be assigned to a sample point that is closest to it. Ultimately, this algorithm generates an approximate description of a Voronoi diagram and a size estimation of each Voronoi cell.
Prediction errors on samples can reflect the accuracy of the surrogate model on sample points and in the local regions around samples. The leave-one-out (LOO) cross-validation approach is adopted to calculate the prediction errors of existed samples. The LOO error of sample pi is defined as: (10) where ( ) By ranking samples according to the sizes of their Voronoi cells, we can find low-density regions and add extra samples in them to help improve the global accuracy of the surrogate model. By ranking samples according to their LOO errors, the high nonlinear and multi-modal regions can be found and adding more samples in them helps to improve the local accuracy of the model. By combining these two methods, a new sequential sampling strategy was obtained, which improved model quality both in global accuracy and local accuracy.
Before combining these 2 methods, the 2 metrics including the sizes of the Voronoi cells and LOO errors should be normalized. The Voronoi cell size was already in the range [0, 1], while the LOO error of each sample had to be scaled to [0, 1], thus the hybrid metric for a sample pi ∈ P is defined as The cell C i is defined as where dom(p i , p j ) represents the closed half plane bounded by the perpendicular bisector of p i and p j . The bisector separates the points closer to p i from those closer to p j . Once the partition was accomplished, it was possible to obtain the density distribution of samples according to the volume of Voronoi cells. Actually, it was difficult to exactly calculate the Voronoi diagram of a given set of samples and design space. The Monte Carlo approach was adopted to realize the approximate Voronoi partition and estimate the volume of Voronoi cells. Firstly, a large number of randomly and uniformly distributed points were generated in the design space. Then for each of these random points, the distances between it and all the existing sample points were calculated. Each random point will be assigned to a sample point that is closest to it. Ultimately, this algorithm generates an approximate description of a Voronoi diagram and a size estimation of each Voronoi cell.
Prediction errors on samples can reflect the accuracy of the surrogate model on sample points and in the local regions around samples. The leave-one-out (LOO) cross-validation approach is adopted to calculate the prediction errors of existed samples. The LOO error of sample p i is defined as: where f (p i ) represents the actual response value of sample p i , ∼ f P\p i (p i ) represents the predicted value of p i calculated by the surrogate model, which is constructed on the current sample set without p i .
By ranking samples according to the sizes of their Voronoi cells, we can find low-density regions and add extra samples in them to help improve the global accuracy of the surrogate model. By ranking samples according to their LOO errors, the high nonlinear and multi-modal regions can be found and adding more samples in them helps to improve the local accuracy of the model. By combining these two methods, a new sequential sampling strategy was obtained, which improved model quality both in global accuracy and local accuracy.
Before combining these 2 methods, the 2 metrics including the sizes of the Voronoi cells and LOO errors should be normalized. The Voronoi cell size was already in the range [0, 1], while the LOO error of each sample had to be scaled to [0, 1], thus the hybrid metric for a sample p i ∈ P is defined as Energies 2019, 12,4695 where V(p i ) denotes the size of the Voronoi cell of p i , and the hybrid metric values for all samples are then used to rank the existing samples. Finally, n new samples were selected from n Voronoi cells, which were ranked the highest. For each of n Voronoi cells, the selection was realized by generating several random points and picking the one farthest away from the sample in the Voronoi cell.
When there were several objectives to be modeled simultaneously, the ratio of the number of newly added points for every objective was allocated by model accuracy during each iteration.

Multi-Objective Genetic Algorithm
Based on the above mentioned adaptive surrogate modeling section, the surrogate models, which were considered to be reasonably accurate to replace the real system, were obtained. A popular multi-objective genetic algorithm called NSGA-II [21] was employed to optimize the surrogate-based system to get the approximate Pareto front set. The NSGA-II adopted a fast non-dominated sorting approach and an elitism mechanism to reduce the computational complexity and speed up the convergence to the true Pareto front. The crowded-comparison operator proposed in the algorithm guided the selection process toward a uniformly spread-out Pareto front.

Adaptive Samples Adding Strategy
Actually, the optimal solutions to engineering MOO problems existed in the local regions of design space. Based on the approximate solutions set obtained by using the NSGA-II algorithm, an adaptive samples adding strategy was utilized to place more samples on the regions around the approximate Pareto front, through which the optimization result could approach the true Pareto front more exactly. The selection of additional samples included 2 steps. As shown in Figure 6, the real samples on the Pareto front were firstly identified through a non-dominated sorting of all the existing samples and the approximate Pareto front points. Secondly, a specified number of samples were selected among approximate Pareto front points to achieve the effect that the selected ones and the real samples on the Pareto front were widely and uniformly distributed on the Pareto front. The selection task was realized by implementing a simulated annealing (SA) algorithm. where V(pi) denotes the size of the Voronoi cell of pi, and the hybrid metric values for all samples are then used to rank the existing samples. Finally, n new samples were selected from n Voronoi cells, which were ranked the highest. For each of n Voronoi cells, the selection was realized by generating several random points and picking the one farthest away from the sample in the Voronoi cell. When there were several objectives to be modeled simultaneously, the ratio of the number of newly added points for every objective was allocated by model accuracy during each iteration.

Multi-Objective Genetic Algorithm
Based on the above mentioned adaptive surrogate modeling section, the surrogate models, which were considered to be reasonably accurate to replace the real system, were obtained. A popular multi-objective genetic algorithm called NSGA-II [21] was employed to optimize the surrogate-based system to get the approximate Pareto front set. The NSGA-II adopted a fast non-dominated sorting approach and an elitism mechanism to reduce the computational complexity and speed up the convergence to the true Pareto front. The crowded-comparison operator proposed in the algorithm guided the selection process toward a uniformly spread-out Pareto front.

Adaptive Samples Adding Strategy
Actually, the optimal solutions to engineering MOO problems existed in the local regions of design space. Based on the approximate solutions set obtained by using the NSGA-II algorithm, an adaptive samples adding strategy was utilized to place more samples on the regions around the approximate Pareto front, through which the optimization result could approach the true Pareto front more exactly. The selection of additional samples included 2 steps. As shown in Figure 6, the real samples on the Pareto front were firstly identified through a non-dominated sorting of all the existing samples and the approximate Pareto front points. Secondly, a specified number of samples were selected among approximate Pareto front points to achieve the effect that the selected ones and the real samples on the Pareto front were widely and uniformly distributed on the Pareto front. The selection task was realized by implementing a simulated annealing (SA) algorithm. There were two criteria to evaluate the space-filling qualities of the set of samples. The first one was the intersite distance, which was the smallest distance between 2 samples of the design.  There were two criteria to evaluate the space-filling qualities of the set of samples. The first one was the intersite distance, which was the smallest distance between 2 samples of the design. and q i = q 1 i , q 2 i , . . . , q d i are normalized to [0, 1] d . The intersite distance criterion of all the samples on the Pareto front corresponding to the selected P is defined as: The second criterion was the projected distance, which reflects the projective property of a set of samples. The property means all the samples are distributed as evenly as possible in any one dimension of the design space. The projected distance criterion corresponding to the selected P is defined as A design with good space-filling quality should satisfy the two criteria simultaneously. By integrating the two criteria, a synthesized indicator can be obtained as: The larger the value of D(P) is, the better the spatial distribution characteristic will be. Thus, the set of samples P selected from the approximate Pareto front should make the value of D(P) as large as possible. The selection of the optimal samples set P can be regarded as a combination optimization problem, and the SA algorithm is suitable for solving the problem.

Numerical Test and Results
In this section, the efficiency and accuracy of the proposed MOO algorithm were evaluated by using 3 test functions [38,39]. The test functions and their characteristics are listed in Table 3. TEST2 was an optimization problem with multiple constraints. KUR was a test function which contains high nonlinear functions especially with the sharp oscillation term sin(x 3 ), the Pareto front of which was nonconvex and disconnected. ZDT3 was a test function which contains a high nonlinear function with five variables, and its Pareto front consisted of five discrete parts. All the objective functions were to be minimized.
The distributions of the modeling samples and optimizing samples in the objective space for the 3 test functions are shown in Figures 7a, 8a and 9a, respectively. The modeling samples were the samples used to build surrogate models in STAGE I. The optimizing samples were the samples selected for finding the Pareto front accurately in STAGE II. The numbers of the 2 types of samples were 20 and 25 in TEST2, 130 and 115 in KUR, and 122 and 78 in ZDT3, respectively. As can be seen from the 3 figures, the proposed algorithm can efficiently add samples around optimal regions and finally approach the true Pareto front.
The optimization result of the proposed algorithm was compared to that of the popular algorithm NSGA-II with the same number of functions calls in Figures 7b, 8b and 9b. Figures 7c, 8c and 9c shows the comparison results between the proposed algorithm and the NSGA-II, which calls enough functions. As the numerical simulations on samples in practical engineering account for most of the time of the entire optimization process, the number of function calls was chosen as the efficiency indicator to evaluate the algorithm. Besides, there were 2 performance metrics to evaluate the accuracy of the algorithm: One metric Υ for measuring the extent of convergence to the true Pareto-optimal solutions, and another metric ∆ for measuring the extent of spread over the true Pareto front and the uniformity of distribution. They are defined as: (15) where N is the number of the obtained Pareto-optimal solutions, d i is the minimum Euclidean distance between the ith solution and the true Pareto front, d f and d l are the distances between the extreme ones of the obtained Pareto solutions and the boundary solutions of the true Pareto front, d i,i+1 is the distance between two adjacent solutions in the obtained Pareto solutions, which are sorted by the values of one objective, d is the average value of these (N−1) distances. A smaller value of Υ means a better convergence, and a smaller value of ∆ means the obtained set of solutions cover a larger area of the true Pareto front and have better uniformity of distribution. The quantitative comparison results of different algorithms for the 3 test functions are listed in Table 4. The optimization result of the proposed algorithm was compared to that of the popular algorithm NSGA-II with the same number of functions calls in Figures 7b, 8b, 9b. Figures 7c, 8c, 9c shows the comparison results between the proposed algorithm and the NSGA-II, which calls enough functions. As the numerical simulations on samples in practical engineering account for most of the time of the entire optimization process, the number of function calls was chosen as the efficiency indicator to evaluate the algorithm. Besides, there were 2 performance metrics to evaluate the accuracy of the algorithm: One metric  for measuring the extent of convergence to the true Paretooptimal solutions, and another metric  for measuring the extent of spread over the true Pareto front and the uniformity of distribution. They are defined as: where N is the number of the obtained Pareto-optimal solutions, di is the minimum Euclidean distance between the ith solution and the true Pareto front, df and dl are the distances between the extreme ones of the obtained Pareto solutions and the boundary solutions of the true Pareto front, di,i+1 is the distance between two adjacent solutions in the obtained Pareto solutions, which are sorted by the values of one objective, d is the average value of these (N−1) distances. A smaller value of  means a better convergence, and a smaller value of  means the obtained set of solutions cover a larger area of the true Pareto front and have better uniformity of distribution. The quantitative comparison results of different algorithms for the 3 test functions are listed in Table 4.  It can be seen from Figure 7 to Figure 9 that the proposed algorithm performs much better both in the convergence and the distribution of solutions than NSGA-II when they were implemented with the same number of function calls. This was also verified by the performance indicators  and  in Table 4. When the number of function calls increased to a greater extent in NSGA-II, the 2 It can be seen from Figure 7 to Figure 9 that the proposed algorithm performs much better both in the convergence and the distribution of solutions than NSGA-II when they were implemented with the same number of function calls. This was also verified by the performance indicators  and  in Table 4. When the number of function calls increased to a greater extent in NSGA-II, the 2  It can be seen from Figure 7 to Figure 9 that the proposed algorithm performs much better both in the convergence and the distribution of solutions than NSGA-II when they were implemented with the same number of function calls. This was also verified by the performance indicators Υ and ∆ in Table 4. When the number of function calls increased to a greater extent in NSGA-II, the 2 algorithms performed almost equally well in accuracy. However, the computation cost of NSGA-II was much higher due to the increase in function calls.
The PSP method was also a kind of MOO method based on an adaptive surrogate modeling technique and had been proven that it can effectively solve MOO problems involving expensive black-box functions. In order to further verify the superiority of the proposed algorithm, the optimization results of the proposed algorithm for the 3 test functions were compared with those solved by the PSP method under the same number of function calls, as shown in Figures 7d, 8d and 9d. Both of the Pareto fronts generated by the two methods for TEST2 were accurate enough. As for KUR and ZDT3 who have more complex characteristics, the optimal solutions obtained by the proposed algorithm can approach more closely to the true Pareto front and have better distribution characteristics compared with those obtained by the PSP method, which can be shown from corresponding figures and performance indicators in Table 4. The comparisons above show that the proposed algorithm can acquire more accurate Pareto fronts compared with the PSP method under the same computation cost.
Through the comparisons with the multi-objective genetic algorithm NSGA-II and the PSP method, we can conclude that the proposed algorithm has satisfactory accuracy and efficiency and is suitable for solving practical engineering problems. In the next section, the algorithm is used to carry out optimization for the studied PMA.

Problem Description
The optimization of the PMA aimed to improve the characteristics of the holding force, breaking operation, and closing operation, which corresponded to the optimizations of the PM holding part, breaking spring, and closing driving part. The optimizations of these 3 parts have different design variables and objective functions. A multi-step optimization method was proposed in reference [40] to optimize the PMA, which divided the whole optimization into 3 sub-optimization modules, including the optimizations of the breaking spring mechanism, PM holding mechanism, and closing driving mechanism. These 3 sub-modules were optimized sequentially to get an optimal design for the PMA. However, the optimization of the closing driving mechanism in the last sub-module focused on only one objective about the terminal velocity of the closing operation, while neglecting the improvement of other performances in the closing process. Besides, the optimizations about the dimension parameters and the excitation circuit parameters were carried out separately, and the surrogate models were constructed based on a one-shot sampling approach. All these mean that the obtained results were not optimal. The proposed algorithm in this paper was utilized to realize the MOO of the closing driving mechanism and attempted to further improve the performances of the closing operation.

Optimization Variables
The closing driving mechanism consisted of a mechanical part and an excitation circuit part, as shown in Figure 10. The main dimension parameters of the mechanism are marked in Figure 10a. The movable core thickness L 1 and coil frame height h 2 were chosen as a part of optimization variables since they were identified as the most influential factors to the closing operation among the 5 dimension variables by sensitivity analysis. Besides, the excitation circuit parameters, including initial voltage of capacitor U 0 , capacitance C and coil turns N were also selected as optimization variables. All the optimization variables are summarized as follows: dimension parameters and the excitation circuit parameters were carried out separately, and the surrogate models were constructed based on a one-shot sampling approach. All these mean that the obtained results were not optimal. The proposed algorithm in this paper was utilized to realize the MOO of the closing driving mechanism and attempted to further improve the performances of the closing operation.

Optimization Variables
The closing driving mechanism consisted of a mechanical part and an excitation circuit part, as shown in Figure 10. The main dimension parameters of the mechanism are marked in Figure 10a. The movable core thickness L1 and coil frame height h2 were chosen as a part of optimization variables since they were identified as the most influential factors to the closing operation among the 5 dimension variables by sensitivity analysis. Besides, the excitation circuit parameters, including initial voltage of capacitor U0, capacitance C and coil turns N were also selected as optimization variables. All the optimization variables are summarized as follows: where Uend is the voltage of the capacitor when the movable core reaches the closed position, rshaft is the radius of driving shaft, Lstroke is the stroke of the mechanism.

Constraint Functions
The closing driving mechanism should firstly ensure that the movable core driven by the electromagnetic force can arrive at the closed position to finish the closing operation successfully. Secondly, the average velocity in 30 mm stroke before the contacts close vav plays a key role in

Objective Functions
Firstly, the terminal velocity of the closing operation v end should be as low as possible in order to reduce the bouncing time and the wear of contacts. Besides, the energy consumption of the whole closing operation E con and the volume of the closing driving mechanism V should also be as low as possible. Hence, all the objective functions are described as: where U end is the voltage of the capacitor when the movable core reaches the closed position, r shaft is the radius of driving shaft, L stroke is the stroke of the mechanism.

Constraint Functions
The closing driving mechanism should firstly ensure that the movable core driven by the electromagnetic force can arrive at the closed position to finish the closing operation successfully. Secondly, the average velocity in 30 mm stroke before the contacts close v av plays a key role in successful closing, and it should be restricted within a specified range. Thus, the constraint functions can be described as: S unfin (x) = 0 (18)  (19) where S unfin denotes the unfinished stroke in the closing operation.

Optimization Results and Analysis
The above mentioned objective functions and constraint functions were regarded as black-box functions except for the calculation of the mechanism volume. They were substituted by the adaptive surrogate models constructed based on the numerical simulation results in the optimization process.
The distribution of all the 398 samples calculated by the numerical simulations in the optimization process in the objective space is shown in Figure 11. It can be seen that the 308 samples used to build the surrogate models in STAGE I of the algorithm were distributed in a large area of objective space, while the 90 samples obtained through 15 times of adaptive samples addition in STAGE II were distributed in a small area where all the three objective values were relatively small. This indicated that the added samples in STAGE II converged to the real Pareto front continually and ultimately a set of satisfactory solutions could be obtained.

Optimization Results and Analysis
The above mentioned objective functions and constraint functions were regarded as black-box functions except for the calculation of the mechanism volume. They were substituted by the adaptive surrogate models constructed based on the numerical simulation results in the optimization process.
The distribution of all the 398 samples calculated by the numerical simulations in the optimization process in the objective space is shown in Figure 11. It can be seen that the 308 samples used to build the surrogate models in STAGE I of the algorithm were distributed in a large area of objective space, while the 90 samples obtained through 15 times of adaptive samples addition in STAGE II were distributed in a small area where all the three objective values were relatively small. This indicated that the added samples in STAGE II converged to the real Pareto front continually and ultimately a set of satisfactory solutions could be obtained. Figure 12 shows that the final Pareto front obtained consisted of 23 non-dominated solutions screened out from all the samples via non-dominated sorting. The final scheme can be chosen from it according to the requirement or preferences of the designers. By considering the three objectives comprehensively, we chose three representative Pareto optimal solutions A, B, and C, as marked in Figure 12, to conduct comparative analysis.    Figure 12 shows that the final Pareto front obtained consisted of 23 non-dominated solutions screened out from all the samples via non-dominated sorting. The final scheme can be chosen from it according to the requirement or preferences of the designers. By considering the three objectives comprehensively, we chose three representative Pareto optimal solutions A, B, and C, as marked in Figure 12, to conduct comparative analysis. successful closing, and it should be restricted within a specified range. Thus, the constraint functions can be described as: where Sunfin denotes the unfinished stroke in the closing operation.

Optimization Results and Analysis
The above mentioned objective functions and constraint functions were regarded as black-box functions except for the calculation of the mechanism volume. They were substituted by the adaptive surrogate models constructed based on the numerical simulation results in the optimization process.
The distribution of all the 398 samples calculated by the numerical simulations in the optimization process in the objective space is shown in Figure 11. It can be seen that the 308 samples used to build the surrogate models in STAGE I of the algorithm were distributed in a large area of objective space, while the 90 samples obtained through 15 times of adaptive samples addition in STAGE II were distributed in a small area where all the three objective values were relatively small. This indicated that the added samples in STAGE II converged to the real Pareto front continually and ultimately a set of satisfactory solutions could be obtained. Figure 12 shows that the final Pareto front obtained consisted of 23 non-dominated solutions screened out from all the samples via non-dominated sorting. The final scheme can be chosen from it according to the requirement or preferences of the designers. By considering the three objectives comprehensively, we chose three representative Pareto optimal solutions A, B, and C, as marked in Figure 12, to conduct comparative analysis.   The comparison of the MOO results with the single objective optimization result [40] and the initial design is listed in Table 5. The MOO of the closing driving mechanism was carried out based on the optimal design obtained by single objective optimization. The PM holding characteristics and breaking characteristics of the multi-objective optimal designs A, B, and C were almost the same as those in single objective optimal design. The key differences existed in the characteristics of the closing operations. It can be seen that the average velocities of closing of A, B, and C satisfied the constraint, and the terminal velocities of the closing and the volumes of the driving mechanism of A, B, and C were further decreased compared to those in single objective optimal design. Only the energy consumptions of closing in A and B were slightly higher than that in the single objective optimal design. The comparison results validated the effectiveness of the proposed MOO algorithm. Besides, each of the three schemes A, B, and C have their own merits, which offers designers the flexibility of choice.

Conclusions
This paper firstly introduces a new mono-stable PM actuator with separated magnetic circuits for application in high voltage VCBs. The PMA can satisfy the fast-breaking requirement for operating mechanisms in high voltage circuit breakers by means of an auxiliary breaking coil and a breaking spring. A MOO algorithm based on the adaptive surrogate modeling technique is subsequently proposed to reduce the number of numerical simulations needed in the optimization and then reduce the optimization time. Test results of mathematical functions show that the algorithm outperforms the multi-objective genetic algorithm in computation cost and is competitive with other efficient surrogate-based MOO methods in accuracy and efficiency. The proposed optimization method is finally applied to the optimal design of the closing driving mechanism of the PMA. By comparing this with the results of the initial design and single objective optimization, the MOO results of the proposed optimization strategy have a notable improvement in the terminal velocity of closing and the volume of the closing driving mechanism. Furthermore, the MOO strategy enhances the performances of the closing driving mechanism more comprehensively and can generate a set of Pareto optimal solutions for selection. The successful application of the proposed MOO algorithm to the PMA shows that the algorithm has great potential in solving practical MOO problems involving time-consuming simulations.
In this research, the RBF model is chosen to construct surrogate models. Other kinds of surrogate models, such as Kriging, support vector regression (SVR), and artificial neural network (ANN), and their applications in the proposed algorithm will be studied in the future. Examining the applicability of the proposed algorithm to high-dimensional problems with many constraints is also an important task.