Robust Covariance Matrix Adaptation Evolution Strategy: Optimal Design of Magnetic Devices Considering Material Variation

Uncertainties caused by material variation can significantly impair the characteristics of devices. Therefore, it is important to design devices whose performance is not significantly damaged even when material variations occur. Robust optimization seeks for the optimal solutions that are robust to fluctuations due to uncertainties caused by material variation, geometrical variation due to assembly tolerances, and changes in physical properties over time in real-world problems. However, naive robust optimization requires iterative calculations to compute the expected values, which need a huge computational burden. This paper introduces a novel robust optimization method for magnetic devices using the covariance matrix adaptation evolution strategy (CMA-ES). In this method, called RCMA-ES (robust CMA-ES), the expected value of the objective function is evaluated using the local average of neighboring individuals without increasing the computation cost. For validation, RCM-ES and robust genetic algorithm (RGA), one of the robust optimization methods without increasing the computational load, was applied to the topology optimization of a magnetic shield and actuator, considering the uncertainty in the BH characteristics. RCM-ES was demonstrated to be particularly more effective for topology optimization with a large number of dimensions compared to RGA and provides robust optimal shapes that are insensitive to variations in BH characteristics.


I. INTRODUCTION
Uncertainties caused by material variation, geometrical variation due to assembly tolerances, and changes in physical properties over time can significantly impair the characteristics of devices and systems [1], [2]. If optimization is performed without considering such uncertainties, the solution cannot be robust to fluctuations. Robust optimization is required to obtain optimal solutions that are robust to fluctuations due to uncertainties. One naive method for robust optimization is based on the Monte-Carlo method, in which the evaluation of an objective function is repeated to obtain the expected value [3]. Although this method provides a The associate editor coordinating the review of this manuscript and approving it for publication was Yu-Da Lin . solution that is highly robust to variation, it incurs a large computational cost because the evaluation of the objective function must be repeated many times. Therefore, robust optimization methods that do not increase computational load are required. For the genetic algorithm (GA), one of the most widely used biology-based metaheuristic optimization methods [4], there is an algorithm called robust GA (RGA), which finds robust solutions without increasing the computational load [5], [6]. Meanwhile, the search performance of the covariance matrix adaptation evolution strategy (CMA-ES) [7] has been confirmed to outperform GA in real problems [8], [9]. CMA-ES, as in GA, is one of the evolutionary algorithms and is classified as a biology-based metaheuristic method [4]. Moreover, the recommended number of individuals for CMA-ES only increases proportionally to log n, where n is the number of optimization variables, whereas that for GA increases proportionally to n. Therefore, the former is more advantageous than the latter for high-dimensional optimization problems. It is, therefore, expected to be widely used in the future. However, there is no effective method for obtaining robust solutions using CMA-ES. In this study, we propose a robust CMA-ES (RCMA-ES), a method that uses CMA-ES to find robust solutions for changes in material properties due to manufacturing tolerances and environmental changes, as well as for fluctuations in material properties. A comparison of the proposed method, RCMA-ES, with other robust optimization methods is shown in Table 1.
For validation, we solved the design optimization problems by considering material variation using RCMA-ES. Because material properties often have uncertainties in realworld problems, it is important to optimize the device shape to account for material variation. Design optimization can be classified into parameter and topology optimization. Parameter optimization requires designers to set design parameters in advance, whereas topology optimization does not require setting design parameters and can produce novel shapes.
In this paper, we propose RCMA-ES, which is a robust optimization method based on CMA-ES, for the design of magnetic devices. RCMA-ES was applied to the parameter and topology optimization of the magnetic shield, and topology optimization of the actuator was performed, considering material variations. The effectiveness of RCMA-ES was verified by comparing the optimization results with those of CMA-ES, GA, and RGA.

II. OPTIMIZATION METHOD
The nomenclature is summarized in Table 2.
A. RGA RGA [5], [6] is an extension of GA to robust optimization without increasing the computational load. A flowchart of RGA is shown in Fig. 1. The Instead of varying the parameter values for each individual to obtain the expected values, RGA generates a single parameter value according to a probability distribution, which is then used to evaluate the individual. Because individuals are generated for each generation in GA, different variations are made to the parameters of these individuals for evaluation. This effectively expresses the  variability in the characteristics due to uncertainty. In other words, the Monte-Carlo method performs spatial averaging to obtain the expected characteristic value, whereas RGA performs time (generation) averaging to obtain the expected characteristic value.

B. CMA-ES
CMA-ES is an evolutionary algorithm that generates individuals from a multivariate normal distribution N m, σ 2 s C and uses the evaluated value of those individuals to update the distribution from which better individuals are generated. m and C are updated as follows [7]:  where x i:λ denotes the ith best individual out of λ individuals. w i is the weight of each individual, giving greater weight to superior individuals, such that w 1 ≥ w 2 ≥ · · · ≥ w λ . The basic procedure is as follows: i. Set hyperparameters. ii. Generate individuals from the normal distribution. iii. Evaluate the individuals. iv. Update m and C of the normal distribution. v. Repeat ii through iv. A flowchart of CMA-ES is shown in Fig. 2 (a). We considered two methods to realize robust optimization based on CMA-ES.
C. METHOD 1 In method 1, the objective function E is evaluated by adding variation to the variables when evaluating individuals, as in RGA, which perturbs an optimization parameter with a noise generated according to a probability distribution, which is then used to evaluate the individual. In the final generation, the individuals are evaluated without variation, and the individual with the best evaluation value is selected as the result. A flowchart of method 1 is shown in Fig. 2

D. RCMA-ES
Simple random variation in the variables in the evaluation of individuals in method 1 could make a negative impact on the distribution update of CMA-ES. In method 2, considering this fact, robustness is considered during the individual evaluation of CMA-ES without increasing the computation cost. It consists of the following procedure to evaluate the pseudo-expected value from neighboring individuals: i. Calculate the objective function value for individuals x (g) i , E i , by adding variation to the variables, as in method 1. ii. As shown in Fig. 3 (a), calculate the distance D j between individuals x (g) i and x (g−1) j at generations g and g − 1 from i , as shown in Fig. 3 where . v. Update the value of the objective function as follows: where α is satisfying 0 ≤ α ≤ 1.
In the final generation, the individuals are evaluated without variation, and the individual with the best evaluation    value is selected as the result. The above process is intended to evaluate the robustness from the local average of the neighbors in the previous generation. Because new function evaluations are not involved, the computing cost can be kept almost the same as that for the original CMA-ES. A flowchart of method 2 is shown in Fig. 2 (c). It would also be possible to perform a similar robustness evaluation using the individuals not in the previous generation but in the current generation while the present algorithm is simpler.

E. SIMPLE TEST PROBLEMS
Methods 1 and 2 were applied to the following two simple test problems [5] to verify their effectiveness: 1) Function f a [5]: It has one broad peak and one steep peak, and is defined by where the domain of x was set to −3 ≤ x ≤ 3. When evaluating the objective function to account for variation, function f a (x + δ) was evaluated by adding Gaussian noise δ ∼ N (0, 0.4) to x. In this study, we used the computer with the configuration shown in Table 3. The optimizations were performed 10 times each with varying random numbers for CMA-ES, RGA, and methods 1 and 2, where the same random seeds were used for each method. Table 4 summarizes the analysis condition. For the crossover and mutation rates of the RGA, the hyper parameters that lead to the best results were adopted after trials. The optimization results are shown in Fig. 4. As shown in Fig. 4 (a), all the solutions of CMA-ES are at the global optimum (1.5 ≤ x ≤ 1.7, f a = 2). From Fig. 4 (c), it can be observed that the optimization results of method 1 have solutions distributed between the globally optimal solution (1.5 ≤ x ≤ 1.7, f a = 2) and the robust solution (−1 ≤ x ≤ 1, f a = 1). The poor performance of this approach can be ascribed to the fact that simply adding variation during the evaluation of the function distorted the normal distribution for CMA-ES, resulting in poor convergence. On the other hand, all solutions of RGA and method 2 are in the robust-solution region with a wide range of peaks (−1 ≤ x ≤ 1, f a = 1). Because RGA can obtain the expected value of individuals by time (generation) averaging, and method 2 can obtain the expected value from the local average of the neighboring by previous individuals, the algorithms converged to a robust solution.
2) Function f b [5], [10]: It has five peaks in the range −1 ≤ x ≤ 1, and is defined as where the domain of x was set to 0 ≤ x ≤ 1. The variation considered in the robust methods was set to δ ∼ N (0, 0.065).
The optimizations were performed 10 times with varying random numbers. Table 5 summarizes the analysis conditions. The optimization results are shown in Fig. 5. CMA-ES has nine solutions at the first peak from the left, which is the global optimal solution, and one solution at the second peak from the left. Similarly, method 1 has nine solutions around the first peak and one solution at the second peak. Because method 1 considers variation during the evaluation, the solutions are distributed slightly off the first peak. However, method 1 has no solution in the third peak from the left, which is the robust solution. Hence, it is concluded that method 1 cannot find the robust solution. However, all the solutions of RGA and method 2 are distributed around the third peak; therefore, it is concluded that they are effective for robust optimization, at least for a simple problem. Here, method 2 uses the pseudo-expected value evaluated from the neighboring individuals in the previous generation. For this reason, the solution can slightly deviate from the peak position as shown in Fig. 5 (d).
From the above discussion, method 2 is referred to as RCMA-ES which is used in the subsequent design optimization.

III. MATERIAL MODEL A. PREISACH MODEL
The Preisach model is used to express the variation in BH characteristics. It is a hysteresis model that is widely used in finite element analysis and can express magnetization properties by adding up basic hysteresis loops. The Preisach model expresses M as follows [11], [12]: where, is a triangular domain in the Preisach plane that satisfies −H s ≤ H v ≤ H u ≤ H s .

B. CONSIDERATION OF MATERIAL VARIATION
In this study, to account for material variation, we assumed a certain function form for the Preisach distribution function K (H u , H v ) in (7) containing three parameters (A, σ 1 , σ 2 ) as follows [12], [13]: Reference initial magnetization curve (original) and 6 initial magnetization curves generated to account for variation (#1-#6). B is expressed as where M is obtained from the Preisach model based on (7). The BH curves are expressed by identifying the unknown variables (A, σ 1 , σ 2 , and µ out r ) in (8) and (9). In this study, the values (A, σ 1 , σ 2 , µ out r ) = (0.81, 23.44, 24.50, 170.21) identified in [13] were used as reference values. The BH curves were generated with variations in the material parameters (A, σ 1 , σ 2 , and µ out r ). The resulting BH curves are shown in Fig. 6. Here, the material parameters (A, σ 1 , σ 2 , and µ out r ) were generated by a Gaussian distribution with the mean as the reference value and the standard deviation as 10% of the reference value.

IV. PARAMETER OPTIMIZATION A. OPTIMIZATION PROBLEM
We considered the optimization of the two-dimensional magnetic shielding system shown in Fig. 7 (a). The aim of this optimization is to prevent the magnetic flux generated by the coil from penetrating into the target area with as little magnetic material as possible. The magnetic flux was analyzed by the finite element analysis, where the coil current and coil turns were set to 10 A and 4000 turns, respectively. The BH curve of the magnetic material of the analysis object was calculated using (8) and (9). As shown in Fig. 7 (b), six design parameters x were set in the design region, and parameter optimization was performed for the double magnetic-shield geometry, where We consider the optimization problem to minimize S mag while maintaining B T ave below B ref , which is set to 0.16 mT as an example, which corresponds to 0.145% of B T ave when the design region is air. The optimization problem is expressed as follows: Introducing the constraint as a penalty term P, we define the objective function E as follows: where W 1 and W 2 were both set to 1 in the following numerical examples.

B. OPTIMIZATION RESULTS
By solving (11) using CMA-ES, RCMA-ES, GA, and RGA, we optimized the design parameter x, where the material model presented in Section III was used. Table 6 summarizes the corresponding analysis conditions. The optimizations were performed for five times each with different random numbers, and the results with the smallest objective function E are adopted. The optimization results are shown in Fig. 8. Compared to CMA-ES, the inner shield obtained using RCMA-ES is thicker in the x-axis direction. In addition, the outer shield obtained using RGA is thicker in the x-axis direction than that of GA. We can also find that since there is no change in the base structure (double shield structure) in this parameter optimization, there are no significant changes in the magnetic flux density distribution and magnetic flux lines. The objective function E is better for the normal methods, CMA-ES and GA, than for the robust methods, RCMA-ES and RGA.
For the optimized shapes shown in Fig. 8, finite element analysis was performed 500 times by varying the material parameters, where the same random seeds were used for each method. The results are summarized in Table 7. RCMA-ES and RGA yielded robust solutions with a smaller mean and deviation than those obtained using CMA-ES and GA. The densities of E are plotted in Fig. 9 (a). It can be seen that the distributions for RCMA-ES and RGA are more sharply peaked in the small E region compared to CMA-ES and GA.   Furthermore, the densities of B T ave are plotted in Fig. 9 (b). The optimal solutions obtained using RCMA-ES and RGA have higher possibilities of satisfying the constraints than those of CMA-ES and GA. The optimized shapes for RGA and RCMA-ES are more robust owing to the thicker magnetic shields.
From the above discussion, we conclude that RGA and RCMA-ES provide robust solutions to material variability fluctuations in parameter optimization without increasing the computational load. In parameter optimization, the number of variables is relatively small, which allows sufficient generations for convergence. In contrast, topology optimization usually has a greater number of optimization variables.   The performance of these methods is discussed in the next section.

V. TOPOLOGY OPTIMIZATION A. FORMULATION
Next, we considered applying RCMA-ES to topology optimization, where we used the NGnet on-off method [14]. In this method, the physical properties of the design region are determined from the output of the following shape function: (12) Moreover, G j (p) is the Gaussian functions given by where p i is the center of the ith Gaussian function. The material attribute of element e is assumed to be magnetic (air) if y p e , x ≥ 0 (< 0), where p e is the position vector of the center of e. Topology optimization is reduced to parameter optimization with respect to x in (12). This topology-optimization method has been shown to be effective for the design of electric motors [15], [16], [17] and wireless power-transfer devices [18].

B. OPTIMIZATION RESULTS FOR MAGNETIC SHIELD
The magnetic-shield-design optimization problem discussed in Section IV was solved using topology optimization, where 96 Gaussian functions were placed in the design region, as shown in Fig. 10, where σ G was set to 7 mm. Topology optimization of the magnetic shape was performed  by solving (11) using CMA-ES, RCMA-ES, GA, and RGA to determine the value of the weighting parameter x in (12). Table 8 summarizes the analysis conditions. The optimizations were performed five times each with different random numbers, and the results with the smallest objective function E were adopted. The optimization results are shown in Fig. 11. The optimized result of CMA-ES has two separated thin shields. The optimized shape corresponding to RCMA-ES has a thicker shield than that to CMA-ES, and the inner and outer shields are connected at the bottom. This geometry is thought to mitigate magnetic saturation owing to the connection at the bottom. The objective function E is the best for CMA-ES, suggesting that the optimized shape obtained using CMA-ES is the best if the material variation is not considered. The objective-function value in Fig. 11 (a) is smaller than that in Fig. 8 (a), indicating that topology optimization allowed us to find better shapes that could not be obtained using parameter optimization. Next, the robustness of the optimized shapes against material variation was verified. For the optimized shapes shown in Fig. 11, finite element analysis was performed 500 times by varying the material parameters. The results are summarized in Table 9. As shown in Table 9, RCMA-ES is a robust solution with smaller mean and deviation values than those obtained using the other methods. Note that RGA did not  yield a robust solution because of the large values of the mean and standard deviation. The densities of E are plotted in Fig. 12 (a). It can be observed that the distribution corresponding to RCMA-ES has the sharpest peak in the small E region, whereas the distribution of CMA-ES has the lowest. This is because CMA-ES converges to the solution that has the lowest value for E but is highly sensitive to the material variation. The densities of B T ave are plotted in Fig. 12 (b). Although the peaks of all methods in Fig. 12 (b) are close to 0.16 mT, RCMA-ES satisfies the constraints in a wider domain than the other methods because the distribution of RCMA-ES has less variance and is less affected by the variations in material properties. The optimized shape corresponding to RCMA-ES is more robust to the material variation because the magnetic material area is slightly larger.
It is concluded that the proposed method, RCMA-ES, provides robust solutions for topology optimization without increasing the computational load, whereas the RGA is ineffective for this problem. Because topology optimization has a large number of variables, RGA could not have a sufficient number of generations for the number of individuals, and thus, could not obtain the expected value of individuals by time (generation) averaging. RCMA-ES could provide the pseudo-expected values evaluated from the neighboring individuals in the previous generation, and thus, could obtain the robust optimal structure.

C. OPTIMIZATION PROBLEM FOR ACTUATOR
Next, we considered the topology optimization problem for the actuator [19] shown in Fig. 13 (a), where 181 Gaussian functions were placed in the design region, as shown in Fig. 13 (b), where σ G was set to 1 mm. This problem aimed to maximize the magnetic energy in the air gap using as small amount of magnetic material as possible. We consider the optimization problem for minimizing S mag while maintaining B T ave , which is related to the magnitude of the magnetic energy in the gap, greater than B ref , where B ref is set to 80 mT which corresponds to 187 times larger than the initial-shape B T ave . The optimization problem is expressed as follows: where W 1 and W 2 were set to 50 and 1, respectively. Furthermore, the coil current and coil turns were set to 0.8 A and 200 turns, respectively.

D. OPTIMIZATION RESULTS FOR ACTUATOR
Topology optimization of the magnetic shape was performed by solving (14) using CMA-ES, RCMA-ES, GA, and RGA to determine the value of the weighting parameter x in (12). Table 10 summarizes the analysis conditions. The optimizations were performed five times each with different random numbers, and the results with the smallest objective function E were adopted. The optimization results are shown in Fig. 14. The optimized shape corresponding to CMA-ES has the magnetic material extending from the coil to the target region, forming a fan shape to cover the gap, which makes the magnetic field in the gap uniform. In addition, the optimized shape obtained using RCMA-ES is thicker than that obtained using CMA-ES. In Fig. 14 (c), the optimized shape corresponding to GA has the region of air at the bottom of the magnetic material. Because this problem contains 181 variables which is large for GA, falling into a local solution as in this case might happen by chance.   For the optimized shapes shown in Fig. 14, finite element analysis was performed 500 times with varying material parameters. The results are summarized in Table 11, and the densities of E are plotted in Fig. 15 (a). As in the case of the magnetic shield, RCMA-ES provides a robust solution with smaller mean and deviation magnitudes than those obtained using the other methods. The densities of B T ave are plotted in Fig. 15 (b). RCMA-ES satisfies the constraints in a wider domain compared with the other methods.
In summary, RCMA-ES provides a robust solution for actuator topology optimization with material variation, without increasing the computational load.

VI. CONCLUSION AND FUTURE RESEARCH
In this study, we proposed a robust design method for magnetic devices. In this method, we evaluated the expected value of the objective function without increasing the computation cost using the local average of neighboring individuals. We demonstrated that robust optimal structures for magnetic devices can be achieved using the proposed method.
We compared the performance of RCMA-ES with those of GA, RGA, and CMA-ES for the parameter and topology optimization of the magnetic shield and topology optimization of the magnetic actuator. The conclusions of this study are summarized as follows: (a) RCMA-ES is effective especially for topology optimization with a large number of variables to obtain a device shape whose performance is robust against possible changes in material property. (b) Simple extension of CMA-ES adopted in method 1, where variables are randomly variated as in RGA to consider the robustness in method 1, does not work well. (c) RCMA-ES is more effective than RGA for large optimization problems. (d) Although the robustness for the solutions obtained by RCMA-ES might be slightly compromised in comparison with those obtained by the naïve robust optimization method based on Monte-Carlo, the former works much faster than the latter. In the future, we will take into account the effects of other variations, such as geometrical variations caused in the manufacturing processes. The geometrical variation can be considered by adding noise to the design parameters for parameter optimization or to the weighting coefficients to the Gaussian functions for TO using the NGnet on-off method. One of the limitations of the proposed method is that, due to the stochastic nature of the method, there is no guarantee that the optimized device can maintain the required property under the uncertainties. This needs to be verified, e.g., by Monte Carlo simulation after optimization.