Multi-objective optimization of coaxial magnetic gears

Abstract

In this paper multi-objective optimization techniques are applied to the design of magnetic gears in order to maximize the performances and to obtain a competitive device. A parallelized stochastic algorithm is implemented and an analytical tool based on magnetic fields computation and losses estimation is used to compute the multi-objective function. This FEM-free computational model can be used in the initial stages of product development.

Keywords

Multi-objective optimization magnetic gears

1. Introduction

In all applications where high torques are required, a common solution is the use of mechanical gearboxes. These devices are characterized by low complexity and high torque to weight ratios but suffer from very strict maintenance requirements. Mechanical gears may also suffer from noise, vibrations and depending on the design parameters a complex lubrication system may be required. Magnetic gears are their electro-mechanical counterparts. Since in these devices the high speed rotor and the low speed rotor have no mechanical contact, the transmission efficiency is high, bulky lubrication systems are not necessary and a maintenance free operation may be guaranteed [1]. The main drawbacks are the greater design complexity and the lower torque to weight ratio that can be achieved unless proper design optimization techniques are used, as in this work. Magnetic gears differ from traditional electrical machines and the design procedures have to be reviewed, but the targets are similar: a good design has to maximize the torque per unit weight and the conversion efficiency while minimizing torque ripple and costs. The design problem is therefore a multi-objective optimization problem, where classical stochastic optimization strategies may be applied if a computationally sufficiently fast and accurate model of the magnetic gear can be provided.

2. Gear fundamental equations

If $p$ is the number of poles of one permanent magnet rotor, $q$ is the number of iron poles and $p_{m,k}$ the number of poles of the other rotor, a magnetic gear has to satisfy Eq. (1):

$\displaystyle\begin{cases}p_{m,k}=|mp+kn_{s}|\\ m=1,3,5,...\infty\\ k=0,\pm\;1,\pm\;2,\pm\;3,...,\pm\;\infty\end{cases}$ (1)

If the iron poles are kept fixed, the relation between the rotor speeds are:

$\displaystyle\omega_{m,k}=\frac{m\cdot p}{m\cdot p+k\cdot n_{s}}\cdot\omega_{r}$ (2)

where $\omega_{m,k}$ is the rotational speed of the rotor with number of poles $p_{m,k}$ and $\omega_{r}$ is the rotational speed of the rotor with number of poles $p$ . According to [1] the couple that results in the maximum interaction is $m=1$ ; $k=-1$ . In this hypothesis magnetic gears are machines capable of transmitting torque with a reduction ratio equal to the ratio between outer and inner rotor poles.

3. Magnetic gear model

Two different approaches are adopted to solve the magnetic problem inside the low and high permeability regions.

3.1 Low permeability regions modelling

Figure 1.

Magnetic gear structure and subdomains. According to [2], the regions are: Inner magnets $\rightarrow\ I$ ; Inner airgap $\rightarrow\ \textit{II}$ ; Outer airgap $\rightarrow\ \textit{III}$ ; Outer magnets $\rightarrow\ \textit{IV}$ .

The formulation adopted for low permeability regions is the one proposed in [2]. The iron permeability is assumed to be infinite and the eddy currents induced in the magnets are neglected. From the 2D Gauss law, the magnetic flux density is $\mathbf{B}=\nabla\times(0,0,A_{z})$ where $A_{z}$ is the $z$ component of the magnetic vector potential. In the magnet regions the material model is:

$\displaystyle\mathbf{B}=\mu_{0}\cdot\mu_{r}\cdot\mathbf{H}+\mu_{0}\mathbf{M}$ (3)

where $\mathbf{M}$ is the magnetization, $\mu_{0}$ and $\mu_{r}$ are the vacuum and relative permeability. Ampere’s law can then be written as:

$\displaystyle\nabla\times\mathbf{H}=\nabla\times\Big{(}\frac{1}{\mu_{0}\mu_{r}% }\nabla\times(0,0,A_{z})-\frac{1}{\mu_{r}}\mathbf{M}\Big{)}=0$ (4)

and therefore in cylindrical coordinates:

$\displaystyle-\nabla^{2}A_{z}=\mu_{0}\nabla\times\mathbf{M}\cdot\vec{e_{z}}=-% \mu_{0}\frac{\partial M_{r}}{r\partial\theta}$ (5)

In the air regions this reduces to the Laplace equation:

$\displaystyle\nabla^{2}A_{z}=0$ (6)

At the inner and outer boundaries the flux lines have to be normal to the iron regions in the hypothesis of infinite permeability:

$\displaystyle\begin{cases}\frac{\partial A_{z}^{I}(R1,\theta)}{\partial r}=0\\ \frac{\partial A_{z}^{I}V(R6,\theta)}{\partial r}=0\\ \end{cases}$ (7)

At the interface between the inner magnets and inner airgap, the magnetic field in cylindrical coordinates is $\mathbf{H}=\frac{1}{\mu}\nabla\times(0,0,A_{z})$ and its tangential component $H_{\theta}$ is conserved:

$\displaystyle\begin{cases}A_{z}^{I}(R2,\theta)=A_{z}^{\textit{II}}(R2,\theta)% \\ \frac{\partial A_{z}^{\textit{II}}(R2,\theta)}{\partial r}=\frac{\partial A_{z% }^{I}(R2,\theta)}{\partial r}\\ \end{cases}$ (8)

At the interface between inner airgap and slots, the continuity of the potential and of the tangential component of the magnetic field have to be guaranteed:

$\displaystyle\begin{cases}A_{z}^{\textit{slot}}(R3,\theta)=A_{z}^{\textit{II}}% (R3,\theta)\\ \frac{\partial A_{z}^{\textit{II}}(R3,\theta)}{\partial r}=\frac{\partial A_{z% }^{\textit{slot}}(R3,\theta)}{\partial r}\\ \end{cases}$ (9)

while at the interface between inner airgap and iron poles the boundary condition is an homogeneous Neumann constraint since the flux lines are normal to the iron regions with infinite permeability:

$\displaystyle\frac{\partial A_{z}^{\textit{II}}(R3,\theta)}{\partial r}=0$ (10)

The same condition of normal flux lines is imposed at the interfaces between iron poles and slots:

$\displaystyle\begin{cases}\frac{\partial A_{z}^{\textit{slot}}(r,\alpha_{1})}{% \partial\theta}=0\\ \frac{\partial A_{z}^{\textit{slot}}(r,\alpha_{2})}{\partial\theta}=0\\ \end{cases}$ (11)

Equations (8) and (9) hold also for the interface between outer magnets and airgap. The PDEs that have to be solved are:

$\displaystyle\begin{cases}\frac{\partial A_{z}}{\partial r^{2}}+\frac{1}{r}% \frac{\partial A_{z}}{\partial r}+\frac{1}{r^{2}}\frac{\partial^{2}A_{z}}{% \partial\theta^{2}}=0\\ \frac{\partial A_{z}}{\partial r^{2}}+\frac{1}{r}\frac{\partial A_{z}}{% \partial r}+\frac{1}{r^{2}}\frac{\partial^{2}A_{z}}{\partial\theta^{2}}=\frac{% \mu_{0}}{r}\frac{\partial M_{r}}{\partial\theta}\\ \end{cases}$ (12)

Using the separation of variables, the general solution inside the air regions is:

$\displaystyle A_{z}(r,\theta)=\sum_{n=1}^{\infty}\Big{(}(E_{n}r^{n}+F_{n}r^{-n% })\cos(n\theta)+(G_{n}r^{n}+H_{n}\cdot r^{-n})\sin(n\theta)\Big{)}+E_{0}\ln r+% F_{0}$ (13)

where all the coefficients have to be calculated for each harmonic $n$ and depend on the geometrical parameters of the gear. An expression similar to Eq. (13) is obtained for the slot regions; in fact the summation is limited to $K$ terms inside the slots and $N$ elsewhere. Substituting Eq. (13) into Eq. (12) and applying the boundary conditions of each subdomain, the differential operators are discretized into a matrix $A$ while the source terms are discretized into a vector $b$ . The solution of the problem $A\cdot x=b$ represents the set of Fourier coefficient for each region of the magnetic gear.

3.2 High permeability regions modelling

For high permeability regions a finite difference approach has been adopted. The magnetostatic equation in terms of vector potential reads:

$\displaystyle\nabla\times\frac{1}{\mu}\nabla\times\mathbf{A}=0$ (14)

which for the 2D case in cylindrical coordinates becomes:

$\displaystyle\frac{\partial^{2}Az}{\partial r^{2}}+\frac{1}{r}\frac{\partial Az% }{\partial r}+\frac{1}{r^{2}}\frac{\partial^{2}Az}{\partial\theta^{2}}=0$ (15)

which is solved using a finite difference method applied to a regular grid in cylindrical coordinates (Fig. 2). In Eq. (15) the first and second derivatives have to be computed. The first and second derivatives with respect to an arbitrary coordinate $\xi$ are expressed as:

$\displaystyle\frac{\partial U}{\partial\xi}=\frac{U_{i+1,j}-U_{i-1,j}}{2\cdot% \Delta\xi}$ (16) $\displaystyle\frac{\partial^{2}U}{\partial^{2}\xi}=\frac{1}{\Delta\xi}\Bigg{(}% \frac{U_{i+1,j}-U_{i,j}}{\Delta\xi}-\frac{U_{i,j}-U_{i-1,j}}{\Delta\xi}\Bigg{)% }=\frac{U_{i+1,j}+U_{i-1,j}-2\cdot U_{i,j}}{\Delta\xi^{2}}$ (17)

Figure 2.

Magnetic vector potential boundary condition for iron regions. $A z$ is computed along the red line and applied on the black one as boundary condition.

Once a regular grid is built inside the iron regions, the previously computed boundary potential $A z$ is applied as Dirichlet condition. When the vector potential is known in the iron regions, the expressions of the magnetic flux densities can be easily computed as $B_{r}=\frac{1}{r}\frac{\partial A_{z}}{\partial\theta}$ and $B_{\theta}=-\frac{\partial A_{z}}{\partial r}$ . Using Steinmetz based approaches [3] the iron losses can be computed, giving an estimation of the gear efficiency. Multiple magneto-static simulations sweeping the inner and outer rotor positions are performed; the magnetic flux densities in each point of the iron regions are stored and consequently iron losses are computed. Eddy current losses inside the magnets are estimated using the approach proposed in [4].

3.3 FEM validation

Table 1 shows the validation of the analytical model with results obtained with a commercial FEM software. Iron losses have been computed using the same Steinmetz based approaches while magnets losses are computed as eddy current losses. The analytical results are in good agreement with the FEM ones.

4. Optimization objectives

In order to obtain a competitive device, the magnetic gear has to be optimized, as discussed in [5]. Depending on the application requirements, the objectives can be different e.g. overall weight, dimensions, torque density, torque ripple, magnet quantity, efficiency. The optimization variables are all the geometrical dimensions, the number of inner or outer poles and the materials used to build the magnetic gear.

The problem to be solved is heavily constrained. The constraints can be divided in two categories: the explicit ones depend directly on variables or functions of variables and their evaluation does not require any simulation results while the implicit ones involve the objectives or other quantities computed through the simulation results. While the former can be checked directly before evaluating the objective function, the latter can be verified only after the objective function has been calculated and thus computational resources have been spent. In the case of real constrained problems, such as the present one, the probability of unfeasible population members is generally high due to both implicit and explicit constraints. The detection of constraint violations at the beginning of the parallel function evaluation reduces memory and processor usage significantly. Depending on the adopted parallel structure, the population size can be chosen to be bigger than the number of cores relying on the fact that in each iteration some individuals will be unfeasible. Depending on the model accuracy, several constraints can be included. The ones that have been implemented are:

$\displaystyle\begin{cases}\textit{Inner\ temperature\ preliminary\ check;}\\ \textit{Outer\ temperature\ preliminary\ check;}\\ \textit{Poles\ stresses\ exceeding\ limit,\ preliminary\ check;}\\ \textit{Radial\ displacement\ preliminary\ check;}\\ \textit{Excessive\ inner\ speed\ preliminary\ check;}\\ \textit{Pole\ stresses\ exceeding\ limit,\ post\ processing;}\\ \textit{Fatigue\ limit\ reached;}\\ \textit{Efficiency\ low;}\\ \textit{Saturation\ high;}\\ \textit{High\ torque\ ripple;}\\ \textit{Inner\ shaft\ resonance;}\\ \textit{Outer\ shaft\ resonance;}\\ \textit{Post\ processing\ outer\ limit\ temperature\ reached;}\\ \textit{Post\ processing\ inner\ limit\ temperature\ reached;}\\ \textit{Post\ processing\ Curie\ temperature\ reached;}\\ \textit{High\ conditioning\ number\ of\ matrix\ A\ in\ Section 3.1}\\ \end{cases}$ (18)

where the preliminary checks refer to the evaluations of constraints in the expected worst cases of rotors alignment. These checks are run at the beginning of the objective evaluation and their violation terminates immediately the objective function calculation.

Table 1

Optimal magnetic gears: semi-analytical (S.A.) validation through a FEM code

Pi	Po	ratio	L (m)	D (m)	Mag segm	$\textit{rpm}_{\textit{in}}$	Mass (kg)	Tool	mat	(Nm/kg)	$T_{o}(\textit{Nm})$		$\eta$
$5$	$42$	$8.4\phantom{0}$	$0.1$	$0.4$	$15\times 2,4\times 2$	$1400$	$43.46$	S.A	Hiperco50	$50.08$	$2177$		$0.996$
								FEM lin.	lossless	$47.92$	$2082$	$.7$	$0.997$
								FEM	Hiperco50	$47.87$	$2080$	$.5$	$0.996$
$4$	$34$	$8.5\phantom{0}$	$0.1$	$0.4$	$15\times 2,4\times 2$	$1400$	$36.38$	S.A	Hiperco50	$45.89$	$1669$		$0.996$
								FEM lin.	lossless	$44.03$	$1601$	$.8$	$0.997$
								FEM	Hiperco50	$41.4\phantom{0}$	$1506$		$0.995$
$3$	$25$	$8.33$	$0.1$	$0.4$	$15\times 2,4\times 2$	$1400$	$27.64$	S.A	Hiperco50	$35.82$	$990$		$0.996$
								FEM lin.	lossless	$34.95$	$966$		$0.997$
								FEM	Hiperco50	$32.54$	$899$	$.3$	$0.996$
$2$	$17$	$8.5\phantom{0}$	$0.1$	$0.4$	$15\times 2,4\times 2$	$1400$	$26.93$	S.A	Hiperco50	$15.61$	$420$	$.4$	$0.995$
								FEM lin.	lossless	$14.81$	$399$		$0.997$
								FEM	Hiperco50	$14.56$	$392$	$.1$	$0.996$

For the sake of brevity only the preliminary temperature constraints and the poles displacements constraints are discussed. The temperature in the inner and outer magnets is limited because of demagnetization issues. Considering the $B H$ curve as a function of the temperature Eq. (19), the maximum working temperature can be found once the minimum magnetic flux density $B_{\textit{lim}}$ in the permanent magnets is known Eq. (20).

$\displaystyle H_{ci}|^{T}=H_{ci}|^{T_{0}}\cdot(1+\beta(T-T_{0}))$ (19) $\displaystyle T_{\textit{max}}=T_{0}+\frac{B_{\textit{lim}}-\mu_{0}\cdot\mu_{r% }\cdot H_{ci}-Br}{\beta\cdot\mu_{0}\cdot\mu_{r}\cdot H_{ci}}\leqslant T_{% \textit{design}}$ (20)

In Eqs (19) and (20) $\beta$ is a material parameter, $H_{ci}$ is the intrinsic coercive force, $T_{0}$ is the reference temperature, $\mu_{0}$ and $\mu_{r}$ are the vacuum and relative magnetic permeability, $B r$ is the remanence flux density and $T_{\textit{design}}$ is the design working temperature. Equation (20) is implemented for both inner and outer permanent magnets.

The mechanical constraints are related to the forces acting on the gear domains which are computed through the Maxwell stress tensor integration. In Eq. (21) the Maxwell tensor $\mathbf{R_{r\theta z}}$ is integrated along the iron poles surface and the value obtained has an upper bound $\mathbf{F_{\text{MAX}}}$ due to the stresses that the material can withstand.

$\displaystyle\mathbf{F_{r\theta z}}=\frac{1}{\mu_{0}}\cdot\int_{S}\mathbf{R_{r% \theta z}}\cdot\mathbf{u_{n}}\cdot dS\leqslant\mathbf{F_{\text{MAX}}}$ (21)

In the simple hypothesis of beams fixed at both extremities, the maximum displacement $\Delta$ occurs in the middle of the iron pole axial length:

$\displaystyle\Delta_{x=\frac{L}{2}}=\frac{F_{r}\cdot L^{3}}{384\cdot E\cdot I}% \leqslant\Delta_{\textit{max}}$ (22)

where $E$ is the equivalent Young Modulus (since the iron poles are laminated along the $z$ axis, $E$ has to be computed for glued iron sheets), $I$ is the moment of inertia, $L$ is the pole axial length, $F_{r}$ is the radial component of the force in a cylindrical coordinate system and $\Delta_{\textit{max}}$ is the maximum allowed displacement in the radial direction.

Further thermal and mechanical models are not discussed here for the sake of brevity but are included in the optimization loop. The algorithm that has been chosen for solving the magnetic gear optimization problem is a multi-objective differential evolution (DE) described in [10]. According to [6], constraint handling in an evolutionary algorithm (EAs) is not straightforward because the mutation and recombination operators are “blind” to constraints. Effectively there is no guarantee that if the parents satisfy the constraints the offspring will satisfy them as well. EAs therefore perform an unconstrained search and are intrinsically not very well suited to handle constrained problems unless such problem is properly addressed.

Constraints can be handled directly or indirectly: in the former approach the constraints are directly applied to the EA while in the latter the distance from the feasible area is added as an additional objective. Using the direct approach unfeasible points can be repaired, eliminated or preserved through special operators [7]. Elimination is the simplest approach but can easily result in a very inefficient procedure. A more sophisticated procedure is the indirect implementation: since the distance from the feasible region is an objective, the algorithm tends to minimize this distance therefore the population drifts towards more feasible solutions. This approach allows unfeasible elements to survive and to explore in a more efficient way the search space.

5. Multi-objective differential evolution

Mutation and crossover are fundamental concepts of the Differential Evolution method [8]. Starting from an initial random set, mutation is applied according to:

$\displaystyle y^{t+1}=x_{a1}^{t}+F\cdot(x_{a2}^{t}-x_{a3}^{t})$ (23)

where $x_{a1}^{t}$ , $x_{a2}^{t}$ , $x_{a3}^{t}$ are the variable vectors of $3$ randomly chosen population members $a1$ , $a2$ , $a3$ taken at the generation $t$ ; the result $y_{i}^{t+1}$ is is the ith of $m$ elements. The crossover operator allows to combine a previous element $x_{i}^{t}$ with a mutated element $y_{i}^{t+1}$ through a random procedure:

$\displaystyle u_{i,j}^{t+1}=\begin{cases}y_{i,j}^{t+1}\ \textit{if rand}_{i,j}% \leqslant CR\\ x_{i,j}^{t}\ \textit{otherwise}\\ \end{cases}$ (24)

where $t+1$ indicates the next generation and the elements $u^{t+1}$ are the result of mutation and crossover operation.

One of the main problems in evolution strategies is the choice of the control parameters: in the case of DE mutation and crossover have to be correctly set in order to prevent premature convergence or stagnation and to ensure acceptable convergence ratio towards the global optimum. Since in DE the perturbation is based on differences between population elements, when the population looses diversity the perturbation vanishes and therefore the algorithm stagnates. In order to avoid premature convergence it is important to calculate the population diversity e.g. by measuring the variance.

According to [9] if $x=(x_{1}...x_{m})$ is the current population, $y=(y_{1}...y_{m})$ is the intermediate population after applying mutation and $u=(u_{1}...u_{m})$ is the population obtained crossing $x$ and $y$ , the expected values $E()$ are:

$\displaystyle E(\textit{Var}(y))=\Bigg{(}2F^{2}+\frac{m-1}{m}\Bigg{)}\textit{% Var}(x)$ (25) $\displaystyle E(\textit{Var}(u))=\Bigg{(}2F^{2}CR-\frac{2CR}{m}+\frac{CR^{2}}{% m}+1\Bigg{)}\textit{Var}(x)$ (26)

where $C R$ is the crossover probability, $m$ is the population size, $F$ is the mutation coefficient and Var is the variance $\textit{Var}(x)=E[(x-\overline{x})^{2}]$ . If $2F^{2}CR-\frac{2CR}{m}+\frac{CR^{2}}{m}<0$ the population variance tends to decrease after mutation and crossover otherwise the variance tends to increase; this equation has been implemented in order to guarantee the variance increase.

The implemented version of multiobjective DE [10] improves the distribution of the population along the Pareto front by means of ideas originating from the NSGA-II algorithm: in each step the current individuals and the ones stored inside an archive are ranked in different fronts and the elements that belongs to the same front are ordered according to their crowding distance.

Figure 3.

Direct implementation: minimum IGD Pareto front of the CF6 function [11].

Figure 4.

Direct implementation with reseeding scheme: minimum IGD Pareto front of the CF6 function [11].

Figure 5.

Indirect implementation with reseeding scheme: minimum IGD Pareto front of the CF6 function [11].

Figure 6.

Direct implementation: minimum IGD Pareto front of the CF7 function [11].

Figure 7.

Direct implementation with reseeding scheme: minimum IGD Pareto front of the CF7 function [11].

Figure 8.

Indirect implementation with reseeding: minimum IGD Pareto front of the CF7 function [11].

Table 2

Direct and indirect optimization comparison on CEC-09 benchmarks [11]. Maximum iterations $n=300\ k$ , $F=0.5$ , $CR=0.55$ , $\textit{Pop. size}=100$ . In the “reseed” versions population reseeding is applied when the algorithm stagnates

Problem	Algorithm	% unfeasible	IGD
			Mean	Min	Max
CF-06	DE-direct	$13.69$	$0.1266\pm 0.047\phantom{0}$	$0.0406$	$0.1984$
CF-06	DE-direct-reseed	$16.8\phantom{0}$	$0.0599\pm 0.029\phantom{0}$	$0.0283$	$0.1182$
CF-06	DE-indirect-reseed	$27.42$	$0.0234\pm 0.0027$	$0.0186$	$0.0299$
CF-07	DE-direct	$12.11$	$0.3255\pm 0.1266$	$0.161\phantom{0}$	$0.6982$
CF-07	DE-direct-reseed	$16.35$	$0.0881\pm 0.0335$	$0.0435$	$0.1583$
CF-07	DE-indirect-reseed	$29.81$	$0.0597\pm 0.0135$	$0.044\phantom{0}$	$0.0858$

Figure 9.

Pareto front. Red points are the current population while the blue ones are previous points stored inside the archive.

Figure 10.

Pareto set. The variables are respectively: Inner radius, Inner iron thickness, Inner magnets thickness, Inner airgap thickness, Iron poles thickness, Outer gap thickness, Outer magnets thickness, Outer iron thickness, Number of inner poles.

Figure 11.

Population storage, violated constraints flag. Yellow points correspond to high saturation value check failure.

Figure 12.

Example of the geometry of a point on the Pareto Front of Fig. 9 with $Pi=5$ .

5.1 Constraints handling

Since multiple constraints handling procedures have been introduced in the algorithm, the benchmark problems $CF6$ and $CF7$ [11] have been used in order to test the performances. The benchmark variable space has the same dimension of the magnetic gear case and the objectives are non-convex subjected to non-convex constraints. The maximum number of iterations in the CEC contest was fixed to $n_{\textit{max}}=300\ k$ .

Three different algorithms have been tested and the results are shown in Figs 3–8 and in Table 2. DE-direct is a simple DE algorithm with the NSGA-II concepts of multiple front ranking and crowding distance. $\textit{DE-direct}+\textit{reseed}$ is similar to the DE-direct but additionally a population reseed scheme has been implemented: when the algorithm stagnates, the population is restarted keeping the previous Pareto front in an archive. $\textit{DE-indirect}+\textit{reseed}$ scheme handles the constraints in an indirect way and reseeding is applied as in the previous algorithm.

Figures 3–5 show the results for $CF6$ . The first part of the Pareto front is well approximated by all the algorithms; the constraints limit the second and third portion of the front. Clearly the simple direct implementation is really inefficient when constraints are present. Adding the reseeding scheme the performance improves but the standard deviation of the results is high. The indirect scheme guarantees better results with a far better standard deviation.

Figures 6–8 show the results relatively to the $CF7$ benchmark. In this case the constraints are much more severe and the first approach can find only few points of the true front. The reseeding scheme gives a much better performance but the optimization results are affected by a high standard deviation. Again, the indirect procedure allows to find a better Pareto front approximation with lower standard deviations.

6. Multiobjective magnetic gear optimization results

In the case studied in this paper, the objectives are: torque density, efficiency and number of magnets of the inner rotor. The geometrical dimensions can vary inside a given range, the number of inner rotor poles is variable between $2$ and $6$ , the gear ratio is fixed around $k\simeq 8.5$ , the iron adopted is 2VPermendur and the magnets are made of high-grade NdFeB with remanent flux density $Br=1.4\ T$ . The inner number of pole pairs is not fixed: indeed once the gear ratio is fixed, there are many combinations of inner and outer pole pairs that result in this gear ratio. The choice is affected by other factors, such as torque ripple or magnets minimum size.

All the constraints in Eq. (18) have been implemented: the design temperature for the initial checks is $T=$ 80 ${}^{\circ}$ C, the maximum stresses in the material have to be lower than $\sigma=200\ \textit{MPa}$ , the maximum displacements are $\Delta_{\textit{max}}\leqslant 0.3\cdot L$ , the efficiency has to be $\eta\geqslant 0.9$ and the torque ripple $T_{\textit{ripple}}\leqslant 15\\$ . The constraints have been handled directly, giving a good Pareto front with $6000$ objective function evaluations.

6.1 Analytic optimal solution

Figure 9 shows the obtained Pareto front. The front appears discontinuous because the inner number of poles has to be an integer number. Increasing the number of magnets the torque density $(\textit{Nm/kg})$ becomes higher. This happens because higher number of magnets implies lower pole pitch therefore the iron yokes become less saturated.

Figure 10 shows the variables set. Red points correspond to the current population while blue ones refers to the stored positions. The picture shows that the optimal gears are not concentrated around a single positions but are rather distributed. From Fig. 10 it appears that the outer airgap size is crucial for the gear performance, while the inner airgap seems to have a lower effect on the selected objectives.

Figure 11 shows the cause of constraint violation. When approaching high torque to weight ratios, the main cause of failure is the over-saturation flag that deeply affects the optimal solution at low number of inner poles.

Figure 12 shows one geometry of the Pareto front. Choosing the torque to weight ratio, the optimal machines tend to be hollow, therefore an electrical machine can be set up in the cavity giving an optimal Vernier machine.

7. Conclusion

In this paper an optimization tool for for coaxial magnetic gears has been introduced and discussed. The optimization algorithm is based on a differential evolution scheme, where special care has been taken in the handling of constraints. The objective function algorithm have been developed in Matlab and validated through a FEM software. This semi-analytic tool is fully parallel and can avoid long-time consuming FEM simulations.

References

Atallah

and Howe

, A novel high-performance magnetic gear, IEEE Transactions on Magnetics 37 (2001), 2844–2846.

Lubin

Mezani

and Rezzoug

, Analytical computation of the magnetic field distribution in a magnetic gear, IEEE Transactions on Magnetics 46 (2010), 2611–2621.

Krings

and Soulard

, Overview and comparison of iron loss models for electrical machines, in: International Conference and Exhibition on Ecological Vehicles and Renewable Energies, Monte-Carlo, Monaco, MAR 25-28, 2010.

Atallah

Howe

Mellor

P.H.

and Stone

D.A.

, Rotor loss in permanent-magnet brushless AC machines, IEEE Transactions on Industry Applications 36 (2000), 1612–1618.

Jungmayr

et al., Magnetic Gear: Radial Force, Cogging Torque, Skewing, and Optimization, IEEE Transaction on Industry Applications 52 (2016), 3822–3830.

Craenen

B.G.W.

Eiben

A.E.

and Marchiori

, How to Handle Constraints with Evolutionary Algorithms, Practical Handbook of Genetic Algorithms: Applications, 2001, 341–361.

Mezura-Montes

Velazquez-Reyes

and Coello Coello

C.A.

, Modified Differential Evolution for Constrained Optimization, in: IEEE International Conference on Evolutionary Computation, Vancouver (BC), 2006, pp. 25–32.

Storn

and Price

, Differential Evolution A Simple and Efficient Heuristic for global Optimization over Continuous Spaces, Journal of Global Optimization 11 (1997), 341–359.

Zahaire

, Critical Values for the Control Parameters of Differential Evolution Algorithms, Proceedings of MENDEL, vol 2002, 2002/6/5.

10.

Alotto

, A hybrid multiobjective differential evolution method for electromagnetic device optimization, COMPEL 30 (2011), 1815–1828.

11.

Zhang

et al., Multiobjective optimization Test Instances for the CEC 2009 Special Session and Competition, Working Report, CES-887, School of Computer Science and Electrical Engineering, University of Essex, 2008.