Many Objective Optimization of a Magnetic Micro–Electro–Mechanical (MEMS) Micromirror with Bounded MP-NSGA Algorithm

: The paper proposes the automated optimal design of a class of micro–electro–mechanical (MEMS) devices, based on a procedure of ﬁnite element analysis coupled to evolutionary optimization algorithms. A magnetic MEMS, used as an optical switch, is considered as the case study. In particular, the geometry of the device is optimized in order to maximize the actuation torque and minimize the power losses and the device volume. The optimization algorithms belong to the genetic class and, in particular, Migrated Parents - Non-Dominated Sorting Genetic Algorithm MP-NSGA, with three objective functions, is compared to NSGA-III.


Introduction
It can be stated that the miniaturization of electromechanical systems is impacting contemporary society as deeply the mass production of electronic systems in recent decades. In fact, the impressive evolution of embedded system technologies has been fostered by the availability of miniaturized devices, which can behave like sensors or actuators: in a sense, these devices are the link between a physical process and the electric circuit processing the analogue signal (sensors) or governing the movement of a mechanical component (actuator) [1]. In this respect, the interest in micro-electro-mechanical system (MEMS) devices has substantially grown after the technological birth of MEMS, dating back to 1964 with the production of the first batch device [2]. However, only recently have the analysis and design of MEMS been approached in a methodologically mature way, based on the precise formulation of forward and inverse problems, which, in turn, might give rise to non-trivial mathematical problems to solve in at least an approximate way.
In the area of MEMS, two main lines of research can be observed. The former one, theoretically oriented, is devoted to the analysis and synthesis of multi-physics models [3][4][5] of, e.g., coupled thermal-elastic systems [3,6,7], electrostatic-elastic systems [3,8,9], magnetically actuated systems, and microfluidic systems [3,[10][11][12][13][14]. In contrast, the latter is more focused on various application areas, for instance, the design and manufacturing of MEMS for biomedical systems with an emphasis on miniaturized bio-sensors and microdevices for tissue engineering [9,15].
Several research projects are characterized by advanced theoretical models which aim at obtaining classes of closed-form solutions [16,17]. Alternatively, problem formulations which prove a priori conditions for the existence, uniqueness, and regularity of the solutions without explicitly computing them, have been developed [18,19]. When an explicit solution is not obtainable, methods finding approximate solutions by means of numerical techniques have been proposed: in particular, if the conditions of existence, uniqueness, and regularity are fulfilled, then the absence of ghost solutions is ensured [6,[20][21][22].
Moreover, techniques for modeling magneto-thermo-elastic devices, considering the wave propagation in micro-domains with fixed or free boundaries, have been studied [23][24][25][26]. Obtaining the conditions of existence and uniqueness of solutions to inverse problems with moving boundaries is the core of several contributions [18,[27][28][29][30]. Due to their simple construction technology, membrane MEMS devices are another class of devices often considered in many engineering applications [5,25,[31][32][33].
There are many techniques of exciting a MEMS actuator; the ones which are commonly used can be classified into three categories, according to their relevant physical principle: (i) The electrostatic excitation, based on an electric field which causes the controlled displacement of a movable component or the deformation of an elastic membrane; (ii) The thermal excitation, which exploits the difference between the thermal expansion coefficients featuring two elastic materials when subject to a temperature gradient; (iii) The magnetic excitation, based on the Lorentz force acting on a loop of current placed in an external magnetic field.
Specifically, magnetic actuation exhibits many advantages. In fact, it has a good linearity of movement versus excitation signal for a broad interval of values of current, low voltages needed for power supply, and hence low power consumption, and finally it is simple to control by means of a sequence of pulses of current [11,34,35]. For these reasons, in the paper, the field model of a magnetic MEMS used as a micromirror is considered, and a numerical method for the optimal shape design trading off multiple design criteria is proposed. The method is a very general one because it is based on the combination of a Finite Element model and a general purpose optimization algorithm. Hence, it is suitable for a broad class of MEMS devices, provided the designer sets up a parametric model of the device and defines the design variables, the objective functions and the constraints.
The manuscript is organized as follows: in Section 2, the principle of operation of the device is described; then, direct and inverse problems are formulated, specifically, the former is solved by means of a 3D finite-element model, while the solution of the latter is based on an algorithm of evolutionary computing inspired by Pareto-like optimality. Eventually, in Section 3, the results are presented and discussed.

The Magnetic Micromirror
A micromagnetic device used as an optical switch [36] is considered. It is based on the actuation of a ferromagnetic plate under the influence of a magnetic field variation. It consists of an NdFeB magnet (neodymium, iron, and boron alloy, Label A in Figure 1), two copper conductors carrying like currents (Label B in Figure 1), and a ferromagnetic plate (Label C in Figure 1) free to rotate around its axis in the point marked with R. The cross-section of the micromirror is shown in Figure 1. The residual induction of magnet is 1.2 T and the relative magnetic permeability is equal to 1.05. The plate is made of ferromagnetic material which is assumed to be linear with a relative permeability equal to 10 3 .
The torque holding the plate at the prescribed angle, evaluated at the point R, is due to the magnetic field of the permanent magnet in the absence of current, while the actuation torque, evaluated in point R, necessary to switch the plate angle is due to the magnetic field variation caused by a current pulse in the conductors. In particular, the actuation torque is calculated as the difference between the total torque (with currents on) and the holding torque (with currents off).

Field Analysis
Given the geometry in Figure 1, the field analysis problem consists of finding the magnetic field distribution for a given plate angle and a given amplitude I of the current pulse. Since it is assumed to neglect the transient behavior, a steady-state analysis is considered. Hence, the problem to be solved, in terms of magnetic vector potential ̅ , is subject to suitable boundary conditions, where ̅ is the specific current, which is supposed to be supplied to the coil by an external source of DC current I. Moreover, µ is the magnetic permeability and is the magnetic remanence of the permanent magnet. The magnetic induction field is then calculated as The torque acting on the plate is calculated by means of the Maxwell's stress tensor [9]: the force F acting on each point Q of a closed surface surrounding the plate is calculated as where is the outward normal unit vector and is the Maxwell's magnetic stress tensor. Since the tensor is symmetric, the torque corresponding to (3) is given by The residual induction of magnet is 1.2 T and the relative magnetic permeability is equal to 1.05. The plate is made of ferromagnetic material which is assumed to be linear with a relative permeability equal to 10 3 .
Typical dimensions of the device, hereinafter called the prototype, are: x 1 = 100 µm, x 2 = 1.2 mm, x 3 = 50 µm, x 4 = 50 µm, x 5 = 200 µm, x 6 = 600 µm, x 7 = 10 A, plate length = 1 mm, and plate height = 25 µm [37,38]. The magnet exhibits a square shaped base, hence the magnet depth in the third dimension, referring to Figure 1, is equal to 2 × x 2 . The plate depth is equal to its length (1 mm) and the two conductors are as long as the magnet in the third dimension, i.e., equal to 2 × x 2 .
The torque holding the plate at the prescribed angle, evaluated at the point R, is due to the magnetic field of the permanent magnet in the absence of current, while the actuation torque, evaluated in point R, necessary to switch the plate angle is due to the magnetic field variation caused by a current pulse in the conductors. In particular, the actuation torque is calculated as the difference between the total torque (with currents on) and the holding torque (with currents off).

Field Analysis
Given the geometry in Figure 1, the field analysis problem consists of finding the magnetic field distribution for a given plate angle and a given amplitude I of the current pulse. Since it is assumed to neglect the transient behavior, a steady-state analysis is considered. Hence, the problem to be solved, in terms of magnetic vector potential A, is subject to suitable boundary conditions, where J is the specific current, which is supposed to be supplied to the coil by an external source of DC current I. Moreover, µ is the magnetic permeability and B 0 is the magnetic remanence of the permanent magnet. The magnetic induction field B is then calculated as The torque acting on the plate is calculated by means of the Maxwell's stress tensor [9]: the force F acting on each point Q of a closed surface surrounding the plate is calculated as where n is the outward normal unit vector and T is the Maxwell's magnetic stress tensor.
Since the tensor is symmetric, the torque corresponding to (3) is given by where r QR is the position vector of a given point Q on Γ with respect to the rotation axis in R, which is the center point of the plate. Based on (4), the holding torque and actuation torque were independently defined and computed. In particular, the holding torque is defined as i.e., the no-load torque due to magnet only, while the total torque is defined as i.e., the on-load torque, due to both magnet and current. The actuation torque is eventually computed as M a = M t − M h . In order to solve the field analysis problem (1)-(4) given the geometry of magnet and conductors (i.e., the forward problem), a 3D Finite Element FE model ( Figure 2) was built with MagNet by Mentor Infolytica [39]; the field simulation was obtained by means of the magnetostatic solver.
where ̅ is the position vector of a given point Q on Γ with respect to the rotation axis in R, which is the center point of the plate. Based on (4), the holding torque and actuation torque were independently defined and computed. In particular, the holding torque is defined as i.e., the no-load torque due to magnet only, while the total torque is defined as i.e., the on-load torque, due to both magnet and current. The actuation torque is eventually computed as = − . In order to solve the field analysis problem (1)-(4) given the geometry of magnet and conductors (i.e., the forward problem), a 3D Finite Element FE model ( Figure 2) was built with MagNet by Mentor Infolytica [39]; the field simulation was obtained by means of the magnetostatic solver. In the real device, the length of the copper wires is equal to the magnet depth. For modelling purposes, however, the copper conductors must start from and end at a face of the airbox; this feature is required by the FE code used, in order to have a geometrically consistent model. This is why in the model, (see Figures 2 and 3) the copper conductors are longer than the magnet.
A typical finite-element mesh, shown in Figure 3, is composed of about 350,000 tetrahedral elements; second-order polynomial Lagrangian elements are considered. The torque-angle curve has been computed based on the Maxwell's stress tensor, for both holding and actuation torque. The torques calculated for the prototype are shown in Figure 4. In the real device, the length of the copper wires is equal to the magnet depth. For modelling purposes, however, the copper conductors must start from and end at a face of the airbox; this feature is required by the FE code used, in order to have a geometrically consistent model. This is why in the model, (see Figures 2 and 3) the copper conductors are longer than the magnet.
A typical finite-element mesh, shown in Figure 3, is composed of about 350,000 tetrahedral elements; second-order polynomial Lagrangian elements are considered. The torque-angle curve has been computed based on the Maxwell's stress tensor, for both holding and actuation torque. The torques calculated for the prototype are shown in Figure 4.  The actuation torque of the prototype was computed for different current values, considering the rotation angle equal to 10 degrees. The torque versus current curve is shown in Figure 5. It can be noted that the actuation torque versus current curve is rather linear for currents up to 100 A. In fact, the linear behavior of the actuation torque is one of the advantages of magnetic MEMS.  The actuation torque of the prototype was computed for different current values, considering the rotation angle equal to 10 degrees. The torque versus current curve is shown in Figure 5. It can be noted that the actuation torque versus current curve is rather linear for currents up to 100 A. In fact, the linear behavior of the actuation torque is one of the advantages of magnetic MEMS. The actuation torque of the prototype was computed for different current values, considering the rotation angle equal to 10 degrees. The torque versus current curve is shown in Figure 5.  The actuation torque of the prototype was computed for different current values, considering the rotation angle equal to 10 degrees. The torque versus current curve is shown in Figure 5. It can be noted that the actuation torque versus current curve is rather linear for currents up to 100 A. In fact, the linear behavior of the actuation torque is one of the advantages of magnetic MEMS. It can be noted that the actuation torque versus current curve is rather linear for currents up to 100 A. In fact, the linear behavior of the actuation torque is one of the advantages of magnetic MEMS.
The values of the torques calculated with the 3D model are in accordance with those calculated in previous papers by the authors [37,38] with the 2D model. For the sake of comparison, let it be considered that the depth of the 2D model is 10 µm, while the plate in the 3D model is 1-mm long in the third dimension: for a rotation angle ϕ = 10 • the holding torque is equal to 2.53 × 10 −8 Nm for the 3D model, while it is equal to about 4.39 × 10 −11 Nm for the 2D model. Therefore, for the 2D model, the holding torque is 0.44 × 10 −8 Nm per mm of depth. The torque resulting from the 3D model (2.53 × 10 −8 Nm) is comparable with the torque of the 2D model (0.44 × 10 −8 Nm). Although a full comparison between the two models is hardly possible, because of the different depth of the various components, the two torque values are in reasonable agreement.

Optimization Problem
The design (or inverse) problem-i.e., find the shape of the MEMS device fulfilling the prescribed performance-is based on three design criteria which, in turn, depend on a set of seven degrees of freedom, here called design variables. The latter control the geometry of permanent magnet and current-carrying conductors, and are subject to a set of constraints. Design criteria are modelled via objective functions dependent on the design variables: this way, the inverse problem can be converted into a constrained optimization problem, subsequently solved by means of an algorithm of evolutionary computing (see Section 2.4).
In particular, the seven-dimensional vector of design variable vector is X = [x 1 , x 2 , x 3 , x 4 , x 5 , x 6 , x 7 ] where x 1 , . . . , x 6 are geometrical quantities (see Figure 1), while x 7 is the value of current carried by each conductor.
Moreover, the following objective functions are defined for a prescribed mirror position ϕ = 10 • : ohmic power loss in conductors, to be minimized with the copper conductivity σ = 5.77 × 10 7 Sm −1 ; Actuation torque acting on the mirror, to be maximized excitation system (i.e., conductors and magnet) volume, to be minimized All in one, meeting the design criteria would help in the design of an efficient and miniaturized device with a good torque performance. From the modelling viewpoint, it can be noted that the computation of f 1 and f 3 is inexpensive (closed-form Equations (7) and (9)), while f 2 is computationally demanding because it is a field-dependent objective (finite-element simulations).
Eventually, three behavioral constraints are prescribed, namely, lower thresholds for holding and actuation torque are set to 10 −8 and 10 −9 Nm, respectively; moreover, the pulsed current density value in the conductors is constrained to be lower than 5 × 10 9 Am −1 . In Table 1, the range of variation in the design variables and the three constraints are shown.
The optimal design problem can be formally cast as follows: starting from an initial guess, identify the family of geometries such that f 1 (X) is minimized, f 2 (X) is maximized and f 3 (X) is minimized, subject to the inequality constraints (Table 1) and to field Equations (1)- (6). The simultaneous optimization of objective vector (f 1 , f 2 , f 3 ) is intended in the Pareto sense; therefore, the family of optimal solutions called non-dominated solutions is the unknown [40]. Traditionally, methods and algorithms for solving bi-objective problems has been developed, while approaches to higher-order objective spaces have been considered more recently [40]. In fact, a tri-objective problem is proposed here.

Optimization Algorithms
As far as the selection of optimization algorithms to solve problem (7)- (9) is concerned, the following remarks can be put forward. The arguments minimizing (7) could be found analytically; the same holds for (9), while the argument maximizing (8) could be numerically approximated by resorting to a zero-order deterministic algorithm, e.g., the Nelder-Mead simplex algorithm or the Powell algorithm of conjugate directions linked to the finite-element solver; this way, the three independent solutions (i.e., two minimizers and one maximizer) would be found. In multi-objective optimization, however, the identification of the Pareto front trading off conflicting objectives is the aim: it cannot be achieved except through resorting to a strategy of evolutionary computing, because deterministic computing does not make it possible the non-dominated sorting of a set of solutions [41], and the tri-objective problem here considered is no exception.
In particular, the selected optimization algorithms belong to the class of genetic Algorithms and they are based on the classical Non-Dominated Sorting Genetic Algorithm, (NSGA-II) [41][42][43][44][45][46]. In particular, Migration -Non-Dominated Sorting Genetic Algorithm M-NSGA includes the Migration concept introduced in [47-49] with a self-adaptation mechanism that allows for the modulation of the size of the migrated population in [50][51][52]. In [47], the migrated population is introduced after the generation step and the new genetic heritage due to the migrated individuals will be active in the next generation if the individuals survive after selection. In this paper, an alternative version of the M-NSGA algorithm, named MP-NSGA, is used [53].
In the Non-Dominated-Sorting Genetic Algorithm strategy, the initial population is composed of a random set of N individuals. The new generation is obtained by applying the genetic operator to the parent population of the previous step: in general, for a pair of individuals, one or two offspring are generated by means of cross-over genetic rules [43]; in the NSGA-II strategy, the new generation is the result of the selection of the individuals that best fit with a selection criteria, i.e., the one that minimizes the objective functions.
In the NSGA-II algorithm at each generation, a population of N parents generates up to N sons (up to 2N individuals after generation step) that, after a selection step, produce a new population of N individuals [43]. In the M-NSGA algorithm, a new population [47,50,51], with a size variable in number of individuals from 5 to N, immigrates after the generation step, so improved individuals could modify the genetic heritage of the original population if they survive after selection. The flow chart is shown in Figure 6a. In the MP-NSGA algorithm, immigration occurs before generation step, so the genetic heritage of the migrated population can affect the successive generation step of the new individuals. Moreover, in this case, the size of the immigrated population is smaller, and the number of individuals is variable between 5 to N/2: then, the generation can generate up to 3N/2 new individuals. The proposed version of the M-NSGA algorithm includes the check of three problem boundaries. In case a solution does not satisfy even one of the constraints, all the objective functions are set to 10 10 . This way, these individuals are discarded during the selection step. The flow chart of the bounded version of the MP-NSGA algorithm is shown in Figure 6b. The algorithm stops after a given number of generations. NSGA-III is able to manage up to 15 objective functions and is based on NSGA-II for the individual generation and classification of Non-Dominated Sorting levels [53]. The algorithm uses reference points to preserve diversity. The initial population is formed by N individuals, with being N a multiple of four larger than H, which is the number of reference points chosen randomly in the design variable space [54]. After each generation step, the operators' non-domination sorting and selection of the reference points are applied to preserve diversity in the population. First, the current minimum of each objective function is searched for and the scaling factors for each objective function are set to the current maximum value for each objective function. Maximum values generate a hyperplane where the reference points are mapped. Each generated individual is assigned to the closer reference point mapped in the hyper-plane following the procedure in [53] related to nichepreservation. In practice, if more points are associated with one reference point, only one is chosen to belong to the new Pt+1 population. This process is repeated until N individuals are selected.

Results
The optimization problem is solved running MP-NSGA with 20 individuals for 50 generations and NSGA-III with the same number of individuals and generations. The objective space is threedimensional, however, a 3D plot of the results does not allow a good evaluation of the results. Hence, for the sake of clarity, 2D orthogonal projections of the objective space into f1-f2, f1-f3 and f2-f3 planes, respectively, are done. They are shown in Figures 7-9. NSGA-III is able to manage up to 15 objective functions and is based on NSGA-II for the individual generation and classification of Non-Dominated Sorting levels [53]. The algorithm uses reference points to preserve diversity. The initial population is formed by N individuals, with being N a multiple of four larger than H, which is the number of reference points chosen randomly in the design variable space [54]. After each generation step, the operators' non-domination sorting and selection of the reference points are applied to preserve diversity in the population. First, the current minimum of each objective function is searched for and the scaling factors for each objective function are set to the current maximum value for each objective function. Maximum values generate a hyperplane where the reference points are mapped. Each generated individual is assigned to the closer reference point mapped in the hyper-plane following the procedure in [53] related to niche-preservation. In practice, if more points are associated with one reference point, only one is chosen to belong to the new P t+1 population. This process is repeated until N individuals are selected.

Results
The optimization problem is solved running MP-NSGA with 20 individuals for 50 generations and NSGA-III with the same number of individuals and generations. The objective space is three-dimensional, however, a 3D plot of the results does not allow a good evaluation of the results. Hence, for the sake of clarity, 2D orthogonal projections of the objective space into f 1 -f 2 , f 1 -f 3       From Figures 7-9, it can be noted that MP-NSGA is able to better explore the search space than NSGA-III. This is a typical benefit of the migration operator introduced in MP-NSGA.
On the other hand, however, NSGA-III appears to better approximate the Pareto front, because there are solutions obtained by MP-NSGA, which are dominated by NSGA-III solutions.
In Table 2, the solutions obtained with the two methods are shown. From Figures 7-9, it can be noted that MP-NSGA is able to better explore the search space than NSGA-III. This is a typical benefit of the migration operator introduced in MP-NSGA.
On the other hand, however, NSGA-III appears to better approximate the Pareto front, because there are solutions obtained by MP-NSGA, which are dominated by NSGA-III solutions.
In Table 2, the solutions obtained with the two methods are shown. In the following, an in-depth analysis of the solutions obtained with MP-NSGA is made. In Figures 10-12, the starting and final solutions are shown; moreover, the prototype and other reference solutions are highlighted in each picture.    The points SPL and FPL highlighted in Figures 9-11 are the solutions with the lowest power losses in the first generation and in the last one, respectively. In Table 3, the values of the design variables and the relevant objective function values are shown. For the sake of comparison, the prototype proposed in [37,38] (P point in Figures 7-9) is also reported.  Table 3 the solutions with the lowest power losses of the first and last generations are shown, in Table 4, the solutions with the highest actuation torque of the first and last iterations are shown, along with the solutions with the lowest system volume of the first and last iterations. The points S PL and F PL highlighted in Figures 9-11 are the solutions with the lowest power losses in the first generation and in the last one, respectively. In Table 3, the values of the design variables and the relevant objective function values are shown. For the sake of comparison, the prototype proposed in [37,38] (P point in Figures 7-9) is also reported.  Table 3 the solutions with the lowest power losses of the first and last generations are shown, in Table 4, the solutions with the highest actuation torque of the first and last iterations are shown, along with the solutions with the lowest system volume of the first and last iterations. However, the solution of the first iteration with the lowest system volume, S SV, is, by chance, the same with the highest actuation torque, S AT. The geometry of the three, up to twenty final solutions, F PL , F AT and F SV, is shown in Figure 13.  The geometry of the three, up to twenty final solutions, FPL, FAT and FSV, is shown in Figure 13.

Discussion
In the past, a bi-objective formulation [37] as well as a tri-objective formulation [38] of the optimization problem, was proposed and solved by the authors. The bi-objective formulation [34] reads: having prescribed lower thresholds for holding and actuation torque, given the plate angle φ = 10°, find the geometry of magnet and conductors, as well as the amplitude of the current pulse, so that the power loss in the conductors and the magnet volume are both minimized. Moreover, there were thresholds for the holding torque and for the actuation torque, set to 1 and 0.25 nNm, respectively.
The tri-objective formulation [35] reads: having prescribed lower thresholds for holding torque (10 −8 Nm) and actuation torque (10 −9 Nm), given the mirror position φ = 10°, find the geometry of magnet and conductors, as well as the current pulse amplitude, such that power loss in conductors (f1), weighted material cost (f2), and excitation system volume (i.e., conductors and magnet, f3) are

Discussion
In the past, a bi-objective formulation [37] as well as a tri-objective formulation [38] of the optimization problem, was proposed and solved by the authors. The bi-objective formulation [34] reads: having prescribed lower thresholds for holding and actuation torque, given the plate angle ϕ = 10 • , find the geometry of magnet and conductors, as well as the amplitude of the current pulse, so that the power loss in the conductors and the magnet volume are both minimized. Moreover, there were thresholds for the holding torque and for the actuation torque, set to 1 and 0.25 nNm, respectively.
The tri-objective formulation [35] reads: having prescribed lower thresholds for holding torque (10 −8 Nm) and actuation torque (10 −9 Nm), given the mirror position ϕ = 10 • , find the geometry of magnet and conductors, as well as the current pulse amplitude, such that power loss in conductors (f 1 ), weighted material cost (f 2 ), and excitation system volume (i.e., conductors and magnet, f 3 ) are simultaneously minimized in the Pareto sense.
However, both optimization problems were based on a 2D model of the micromirror, while, in this paper, a 3D model is considered. This improvement allows good accuracy in the solution of the direct problem. In fact, the geometry of the system is clearly a 3D one due to the different depths of the device components: only a 3D model is able to accurately simulate the field distribution in the model and, hence, relevant quantities, e.g., forces and torques. A comparison of the torque calculated by means of both models (2D and 3D) is made (see the end of Section 2.2); there is an agreement between the two models, but the 3D model is more reliable because the field distribution in the third dimension is taken into account.
Moreover, the tri-objective formulation was solved by means of P-EStra, a Paretian Evolution strategy of lowest order. For each optimization run, only one solution was found. In this paper, a migration-assisted algorithm (MP-NSGA), which runs with a population of N individuals at a time, is used; the Paretian solution found is then composed of N = 20 individuals.
Looking at the results, the following comments can be put forward. From Table 3, it can be noted that the solutions with the lowest power losses are characterized by relatively thick magnets (high values of magnet thickness x 1 ), and rather thin coils (low values of x 4 ) carrying a low current (low values of x 7 ). The length of the coils and the distance between them (design variables x 5 and x 6 , respectively) have values in the middle of the feasible range. From Table 4, it can be noted that the solutions with a higher actuation torque have thick magnets (high values of magnet thickness x 1 ), small intercoil distance (low values of variable x 6 ) and high current values (variable x 7 ).
In turn, the dimensions of the coils (design variables x 4 and x 5 , respectively) have values in the middle of the feasible range.
From Table 4, it can be also noted that the solutions with lowest system volume are characterized by small magnets and thin coils located close to the magnet (low values of variables x 1 , x 2 , x 3 , x 4 ). In turn, the coils are larger and quite distant from each other; in fact, variables x 5 and x 6 do not affect the system volume calculation. The current in the coils is the maximum possible, given the current density constraint.
In general, for all the Pareto-optimal solutions identified, it can be noted that short magnets (low values for magnet length x 2 ) are preferred and the distance between the coils and the magnet (design variable x 3 ) is kept as small as possible.