Topology Optimization for Electromagnetics: A Survey

The development of technologies for the additive manufacturing, in particular of metallic materials, is offering the possibility of producing parts with complex geometries. This opens up to the possibility of using topological optimization methods for the design of electromagnetic devices. Hence, a wide variety of approaches, originally developed for solid mechanics, have recently become attractive also in the field of electromagnetics. The general distinction between gradient-based and gradient-free methods drives the structure of the paper, with the latter becoming particularly attractive in the last years due to the concepts of artificial neural networks. The aim of this paper is twofold. On one hand, the paper aims at summarizing and describing the state-of-art on topology optimization techniques while on the other it aims at showing how the latter methodologies developed in non-electromagnetic framework (e.g., solid mechanics field) can be applied for the optimization of electromagnetic devices. Discussions and comparisons are both supported by theoretical aspects and numerical results.

ments in additive manufacturing and 3D printing technolo-23 gies [12], [13], [14]. The associate editor coordinating the review of this manuscript and approving it for publication was Debdeep Sarkar . est in applications where there is no prior knowledge 30 on the optimal arrangement or shape of the material in 31 the device. 32 General optimization problems can be tackled follow-33 ing two main approaches [15]: gradient-based and gradient-34 free techniques, each one characterized by pros and cons. 35 Gradient-based methodologies typically exhibit a faster 36 convergence speed, but suffer for local optima trapping. 37 Gradient-free methods on the other hand are more flexi-38 ble since they do not require the knowledge of the objec-39 tive function derivatives with respect to the design variables 40 (sensitivities), usually not easy to compute, sand generally 41 require a higher computational effort, especially in the case of 42 stochastic algorithm. Alternative to these methods, determin-43 istic gradient-free ones, such as the simplex gradient method 44 or the simplex-base simplicial Nelder-Mead method, can be 45 adopted [16]. 46 Denoting with ρ(x) the material distribution at any point 47 of the design domain D ⊂ R d , the TO problem can 48 When facing engineering problems, the optimal design of 100 components often requires conflicting demands to be satis-101 fied for their best functionality [24], [25], [26]. This reflects 102 into the need of solving multi-objective optimization prob-103 lems, where the scalar objective functions F i are formally 104 collected in vector form and the optimal solutions lie on the 105 so-called Pareto front. 106 When only two conflicting objectives are considered it 107 may be convenient to rewrite the vectorial multi-objective 108 optimization problem into a scalar one through the convex 109 combination where α i ∈ R are weights in [0,1]. These weights can be 112 retrieved using the Adaptive Weighted Sum (AWS) scheme, 113 which for some problems complies also with non-convex 114 Pareto fronts [27]. Since the AWS scheme at each iteration 115 treats only scalar optimization problems, this vectorial-scalar 116 transformation may be conveniently applied when using 117 ''black-box'' numerical optimization tools for which only a 118 scalar objective is required as input. For example the opti-119 mization toolbox of MATLAB R can be used for the single 120 objective optimization and the results combined following 121 AWS scheme.

122
The aim of this work is to address in detail the charac-123 teristics of the different TO methods presented so far and to 124 provide practical application examples. Moreover, it is shown 125 how numerical tools initially developed for e.g., mechanical 126 applications, can be also applied to electromagnetic prob-127 lems. In this respect, section II highlights the main peculiari-128 ties of TO approaches developed in the solid mechanics field, 129 and the ones that can be efficiently used for the electromag-130 netics problems. 131 In particular the gradient-based method described in the 132 following are: homogenization method III-C1, density-based 133 method III-C2, level-set method III-C3, phase-field III-C4 134 method, Bi-directional Evolutionary Structural Topology 135 Optimization (BESO) method III-C5, Topology Optimization 136 of Binary Structures (TOBS) method III-C6 and the Two-Step 137 method III-C7, which is a bridge between gradient-based and 138 gradient-free approaches. The methods considered in the lat-139 ter class are: boolean methods IV-A1, binary methods IV-A2, 140 Normalized Gaussian Network Methods (NGnet) IV-A3, 141 deterministic methods IV-B, deep learning methods IV-C and 142 the Proportional Topology Optimization (PTO) method IV-D. 143 The remainder of the paper is organized as follows: 144 section III describes the gradient-based methods requiring 145 the computation of objective function sensitivities. To this 146 aim, a brief introduction of the so-called adjoint variable 147 approach is given in section III-A. Next, methods that do 148 not require the evaluation of the gradient of the objec-149 tive function, thus avoiding the knowledge of classical sen-150 sitivities, are described in section IV. Section V briefly 151   The topology optimization, although it was mainly devel-163 oped in the structural mechanics framework, can be con-164 ceived in a transversal way based on physical problems of 165 different nature, for example fluid dynamics [49], heat trans-166 fer [50], photonics design [51] and so on. The investiga-167 tion of TO for electromagnetics (EM) is discussed in this 168 paper. The wide spectrum of TO applications may lead to 169 the erroneous conclusion that any method can be used ''as 170 is'' independently of the underlying physics. In fact each 171 physical problem has specific requirements to be satisfied.

172
For instance, when dealing with solid mechanics, the opti-173 mized topology has to comply with connectivity properties.

174
In this class of problems, for which the compliance is the problems. Generally, the gradient-based methods such as the 184 SIMP described in section III-C2 can be used as general 185 purpose, conversely, the level-set of section III-C3 requiring 186 the topological derivative information, is more challenging.

187
The PTO algorithm, developed for structural mechanics prob-188 lems and briefly described in section IV-D, uses deeply the 189 assumption that the objective function is an energy density 190 and cannot be considered as general purpose because if the 191 electromagnetic problem does not uses the energy density 192 (e.g., of magnetic nature) as objective function, the algorithm 193 cannot be applied ''as is ''. 194 From the computational point of view, as briefly stated in 195 the introductory section, the numerical solution of the PDE 196 governing the physical problem is a key part in the opti-197 mization process. In fact, the objective function evaluation 198 requires the solution of the discretized system many times 199 during the optimization algorithm. The FEM is extensively 200 used in solving mechanical and structural problems. The main 201 advantage of this approach is the sparsity pattern of the arising 202 system matrices, which reduces the memory requirements for 203 their storage. In most cases, this approach is used even for 204 the numerical solution of the PDE governing the electromag-205 netic problem. When modelling antenna propagation or eddy 206 currents the need of discretizing the non-conductive parts, 207 for example the air domains, can be circumvented by using 208 e.g., IEM. These methods require only the discretization 209 of the non-electromagnetic neutral domains, but the arising 210 system matrices are fully populated (dense), thus increasing 211 the computational burden for their storage up to O(N 2 ), 212 where N is the number of unknowns, and to O(N 3 ) for the 213 system solution (although some sophisticated technique is 214 used [52]). In principle the TO method can be based both 215 on FEM and IEM approaches for the solution of the EM 216 problem. However, the development of the code in the two 217 cases is different and also the general performance. As an 218 example consider a TO problem which changes the material 219 from air to magnetic or the other way round. By using FEM 220 the number of degrees of freedom of the whole discretized 221 domain (air + magnetic) remains unchanged, conversely, 222 by using e.g., an IEM when subtracting or adding magnetic 223 material, the number of unknowns is respectively reduced or 224 augmented, this because the IEM deals only with the non-air Computing the second term of (11) is a complex task since 273 the derivative of the solution array with respect to the design 274 variables is required. However, if the coefficients λ k are 275 chosen in such a way that the second term of (11) vanishes. Thus, the evaluation of the 278 objective sensitivities is a two step process, first requiring the 279 solution of the the so-called adjoint problem defined by (12), 280 from which the adjoint field λ is obtained. Then, the latter is 281 used for the evaluation of the first term of (11), which is equal 282 to the objective sensitivities dF/dρ i .

283
As an example, considering for simplicity a system whose 284 coefficients depend only on the design parameters ρ, the 285 discretized PDE can be written as The arising adjoint system (12) is written as 292 and, once the solution array λ is obtained, the sensitivity of 293 the objective function is given by The adjoint variable method can be extended to non-296 linear [55], [57] and time-domain problems such as transient 297 eddy current ones [58]. In this section, the interpolation functions mapping the design 301 variables ρ to the elemental material property (5) are defined. 302 To mitigate the oscillatory behaviour of material properties, 303 thus improving the numerical stability, filtering techniques 304 are described together with the projection schemes used to 305 reduce the ''gray scales''.  Denoting with g(x) the material property at the point 322 x ∈ , e.g., the relative permeability µ r (x) for magnetic 323 problems, the classical MIS interpolates g from the contin-324 uous density ρ(x), using a power-law [59] 325

345
In solid mechanics TO, the Rational Approximation of Mate-346 rial Properties (RAMP) proposed in [60] is usually adopted.

347
Given the parameter q ≥ 0, the RAMP interpolation is 348 expressed as In [61], a different interpolation method is proposed by D.

351
Lukàš through the scheme   A graphical comparison of the aforementioned interpo-360 lation schemes is reported in Fig. 1, while an illustrative 361 example showing the spread of material property g when 362 different values of the penalization parameter α of classical 363 MIS are used is reported in Fig. 2.

364
The advantage of using a continuous design variable 365 instead of a discrete binary one, lies in the possibility of 366 computing derivatives for sensitivity analysis [30]. Although 367 it would be better if the optimization problem was binary (i.e., 368 ρ ∈ {0, 1}), it is in any case convenient to work continuously 369 with the drawback of having to manage the gray scale issue. 370

371
Using ρ directly as input of the MIS, may produce an oscilla-372 tory behaviour of the material distribution over the finite ele-373 ment discretization [63], therefore a spatial filtering function 374 has to be adopted. The density filter for each mesh element i 375 can be written as [64]  The color bar refers to the value of µ r which can vary within [1,1000] in the design domain, whose boundary is highlighted in red.
where N i,j is the set of neighbourhoods of the ith cell within 378 the filter radius R and w(x j ) is a weight function between cells 379 i, j: The effect of filtering on the material density ρ is illustrated 391 for the case of the Helmholtz filter (24) in Fig. 3.
where the constant c is (25)

414
In this section, the gradient-based methods listed in 415 Table 1  In the homogenization method, each finite element consti-434 tuting the discretized design domain D, is composed of an 435 infinite number of microstructures usually formed by rectan-436 gular cells with rectangular holes [2]. The size of each hole 437 and its rotation angle are the design variables of the opti-438 mization. In homogenization logic, a cell becomes ''solid'' 439 if the the hole vanishes, conversely is ''void'' if the hole 440 has the same dimension of the cell. Thus, the material is 441 changed only in the microscale cell, and not in the whole 442 mesh element. In the optimization process, the material is 443 transferred between different parts of the design domain 444 giving, at the end, the optimal material distribution. This 445 approach currently seems to be abandoned in favor of the 446 others described below. 447 98598 VOLUME 10, 2022  The density method can be realized also within commercial (28) 483 Performing topology optimization with LSM, means track-484 ing the evolution of the level-set function φ(x) solving the 485 Hamilton-Jacobi partial differential equation for a fictitious 486 time t [75] where v is the velocity field When the level-set function is a signed distance function, 495 the level-set function at iteration k +1, or equivalently at time for k > 2, with τ < 1 a regularization parameter.

508
Equation (29)  is then reformulated as  To overcome this issue, the initial design domain should be 527 seeded with holes or a mechanism for hole nucleation during 528 the optimization process would need to be adopted [82]. :  A large value of η selected ∈ [0, 1] can speed up the conver-556 gence, but at the cost of increasing the probability of reaching 557 a non-optimal solution [87]. A MATLAB R numerical imple-558 mentation of the aforementioned logic for the compliance 559 minimization problem can be found at https://www. 560 topopt.mek.dtu.dk/. The phase-field method is a variation of the standard level-set 563 approach aiming at regularizing the topology optimization 564 problem. In the phase-field approach the level-set function 565 is re-defined as [89], [90] 566 These algorithms in general belong to the class of hard-586 kill methods, however as described in [41], the sensitivity 587 information ϕ i can be used to change the elemental density 588 as where ϕ th rem and ϕ th add are real numbers used as threshold 592 values to decide whether add or remove.  At the kth iteration, the problem to be solved is written as [44] 602 Minimize 2) The configuration achieved in the first step is further 632 optimized using methods involving sensitivity analysis.

633
In [48], the local search is performed evolving the mate- proposed.

639
In this section, the use of gradient-free methods avoiding devices. The section generally distinguishes between stochas-643 tic IV-A and deterministic IV-B approaches. The novel 644 deep learning methods based on neural networks are briefly 645 described in section IV-C, followed by the recently proposed 646 proportional topology optimization (PTO) method IV-D for 647 solid mechanics.

648
A list of the described methods is reported in Table 2.

686
In the class of binary methods we group all the methods  ON/OFF topology optimization problems [127]. with the latter recently applied for the topology optimization 720 of wireless power transfer (WPT) devices [132].

722
The Normalized Gaussian network (NGnet) introduced by 723 Sato et. al. in 2015 allows smooth shapes without introducing 724 additional filtering [106], [107]. Considering a finite ele-725 ment discretization of the computational domain, the material 726 property g to be assigned to each mesh element is determined 727 by the value of the shape function f (x) defined as where G i is the Gaussian function centered at the ith element 731 barycenter, b i is the ith normalized function and w i is the 732 ith weight. For the two-material case (e.g., iron and air), the 733 elemental material is expressed as [133] 734 but extensions are available also in multi-material prob-736 lems [134]. The weight vector w = [w 1 , . . . , w Ng ] is deter-737 mined by an evolutionary algorithm (e.g., µGA) in such a way 738 the optimization problem (1) is satisfied. That is, the original 739 problem becomes a parametric optimization problem [132]. 740 An illustrative example of the procedure determining the 741 shape function f (x), is shown in Fig. 8.

747
In this section, we briefly discuss the deterministic (direct search [108], [136]. An example of flowchart of hybrid 760 stochastic-deterministic algorithm can be seen in Fig. 9.

761
From the author's knowledge, its seems that electromag-  [137], [138], 776 [139], [140], [141], Kriging methods [142], [143], response 777 surface methods [144], and Space Mapping methods [145] 778 have been built, even for multi-objective problems [146], 779 [147], [148]. The Proportional Topology Optimization (PTO) is an heuris-790 tic non-sensitivity based method for solid mechanics applica-791 tions proposed by Biyikli and To [116]. Even in this method, 792 the material property is interpolated according to the MIS 793 approach (17), while the elemental densities during the itera-794 tion procedure are updated as where ρ opt e is the optimized density proportional to the objec-797 tive function, and α is a history parameter controlling the ratio 798 of dependence of elemental density to its older value from the 799 previous iteration [116]. 800 With the focus of increasing robustness and capability of 801 approaching binary distributions, some improvements have 802 been recently proposed [151].

803
A MATLAB R implementation of the PTO algorithm 804 applied to classical solid mechanics examples is available 805 at http://www.ptomethod.org/. Up to now, to the 806 author's knowledge, this method has not been used for the 807 topology optimization of electromagnetic devices.

809
Since the research on TO for electromagnetics is rapidly 810 evolving, a variety of approaches were developed in the recent 811 VOLUME 10, 2022 years, and some of them are summarized in this section.

812
In [152] the Allen-Cahn equation is used to update the design 813 variables using the phase-field method of section III-C4.

839
In this section, numerical results obtained with selected meth- example is chosen to show the capability of the proposed 867 approaches of treating 3D non-trivial geometries which will 868 be of great interest for future industrial applications. In the following sections some selected gradient-based III 880 and gradient-free IV approaches are used to solve prob-881 lem (50) and the results are compared. The objective function 882 improvement is measured with the following metric where W 0 m is the magnetic energy when the whole design 885 domain is filled with ferrite. 886

887
Due to their general interest and wide application, 888 even in electromagnetics, the density based III-C2 and 889 level-set III-C3 approaches are used to solve TO problem. 890 Moreover, the novel TOBS method III-C6 which to the 891 author's knowledge, has not yet been used in electromag-892 netics, is considered. The density-based approach is imple-893 mented in COMSOL R under the topology optimization node 894 and the classical MIS approach (17) used for the material 895 interpolation. Then, then minimization problem is solved 896 using the Method of Moving Asymptotes (MMA) [162]. 897 The developed level-set method combines COMSOL R with 898 MATLAB R environments. The physics is solved within 899 COMSOL R , while the optimization proceeds in MATLAB 900   Optimization results are reported in Table 3 and the final 911 material distribution is illustrated in Fig. 11.       The results are reported in Table 4 and the final material 942 distribution illustrated in Fig. 12 Fig. 13, the design domain D lies below the 948 Ground Assembly (GA) coil. There, the relative magnetic 949 permittivity is the subject of topology optimization, that is 950 µ r ∈ {1, µ FE }. The geometrical parameters are reported in 951 Table 5.

952
The topology optimization aims at maximising the cou-953 pling coefficient k, while keeping the amount of ferrite below 954 VOLUME 10, 2022 FIGURE 13. CAD view of WP1-Z1 3D device. Due to symmetry, only a quarter of the design is shown. The design domain D highlighted in red, correspond to the ferrite in the Ground Assembly (GA). GA and Vehicle Assembly (VA) coils are colored in cyan and green, respectively, while gray parts are aluminum shielding. half of the volume of the design domain, i.e.,  the design variable [55], thus requiring multiple constructions 988 of the system of equations. When a density method is used 989 (e.g., SIMP), the selection of MIS function penalizing the 990 material property, is extremely important to reduce numerical 991 instabilities. In addition, when the objective function exhibits 992 many local minima, these approaches may remain stuck and 993 the true global minimum may not be reached. The latter 994 is a well-known problem when dealing with gradient-based 995 optimization techniques.

996
To overcome this issues gradient-free approaches may be 997 adopted IV. As highlighted in the test-case results of Table 4, 998 one of the main disadvantages of these methods is the high 999 number of objective function evaluations, which depend on 1000 the population size. Reducing the population size in general 1001 may produce non-satisfactory results thus this approach can-1002 not be followed to improve the performances. Due to the 1003 stochastic nature of these approaches it is difficult to ensure 1004 domain connectivity or avoid checkerboards layouts, and the 1005 situation is even worse in cases of final topologies where mul-1006 tiply disconnected components are generated. As expected, 1007 the checkerboard pattern is highly evident. To overcome this issue, the user can ensure the connectivity between selected parts using specific techniques, but this approach is extremely 1011 problem dependent and unfeasible in cases when no prior 1012 information is known on the final topology. An example 1013 for this can be found in the WPT TO of [132] where the  The paper follows the general distinction between 1037 gradient-based and gradient-free methods. Although numer-1038 ical TO is highly problem dependent, so that it is difficult 1039 to assert that an approach is better or worse than another, the algorithm is coupled with commercial software that 1063 already implement the computation of sensitivity maps. 1064 If the sensitivity of the objective function is extremely 1065 costly from the computational viewpoint, a gradient-free 1066 approach may be preferred due to its simplicity, requiring 1067 only objective function evaluation. With these methods some 1068 forethought must be used to ensure domain connectivity, 1069 thus avoiding checkerboards patterns. Moreover, when using 1070 stochastic-based methods, the number of objective function 1071 evaluations is highly increased with respect to gradient-based 1072 ones. When the evaluation of objective function becomes very 1073 costly from the computational point-of-view, as is the case of 1074 large systems of equations, the novel approaches using neural 1075 networks can be adopted to reduce the computational burden 1076 through the creation of surrogate models. [20] R. Torchio, ''A volume PEEC formulation based on the cell method for electromagnetic problems from low to high frequency,'' IEEE Trans.