An approach to an optimal design of permanent magnet synchronous machines for battery electric vehicles

The described design of a permanent magnet synchronous motor for battery electric vehicles is based on a drive cycle. For a given vehicle a set of operating points of the electric machine can be derived from the cycle. The next step in the design of the traction motor is to find a suitable lamination layout. Common motor design procedures take only a few representative operating points into account. This paper describes an approach which includes a large number of operating points at early stages of the motor design process. It is based on a specialized motor model and an optimization strategy which are described in this paper. As an example a traction motor for the VW Golf CityStromer is designed by using the proposed method.


Introduction
In this paper an approach to the design of the electric traction motor of a battery electric vehicle is described. As an example the design process is demonstrated for a permanent magnet synchronous machine (PMSM) used in a VW Golf CityStromer [1]. Input of the design is a drive cycle measured during a trip from Berlin-Centre to a suburban location. A car having a 75 kW combustion engine was used to measure the cycle. The design goal is to find an optimal lamination layout for the given drive cycle. An example for a lamination layout is shown in fig. 5. Since the space inside the CityStromer is limited, the new motor has to have the same size as the original motor. Especially the outer diameter is set constant in this paper. The electric machine is operated at different torque and rotational speeds. Each operating point of the electric machine is defined by a certain torque T op and speed n op . These points depend on the drive cycle and the battery electric vehicle itself [8]. The distribution of the operating points in the torque speed plane is shown in fig. 1. Each point is marked by a black cross. The red line indicates the desired maximum torque of the new motor. An electric machine designer has to define the parameters of a lamination layout as it is shown in figure 5.
If the electric machine is optimized for a single point of operation this design may not be suitable for the whole set of operating points. For example figure 2 shows a comparison of two lamination designs which are determined by the optimization algorithm described in section 4. Each optimization has been carried out for a single operating point. The first operating point is given by T op = 52 Nm and n op = 3800 rpm and the second by T op = 20 Nm and n op = 5400 rpm. Normally there are many operating points in the torque speed plane. State of the art design tools, such as the finite element method, need computational effort. Including all points of operation would end in a very time consuming process. There are methods which allow to reduce the number of points by choosing a limited number of representative points [7]. These methods are typically based on the energy consumption. The torque speed plane is roughly divided into a small number of areas and for each area an representative operating point is calculated similar to the procedure described in section 2. The number of regions is limited due to the computational effort. In this paper an approach is presented which allows taking a large number of operating points into account. The computational effort is reduced by using a specialized motor model in combination with an optimization algorithm. This can be used to divide the torque speed plane into a larger number of areas.

Relative energy throughput
In a first step the motoring area of the torque speed plane is divided into 15 x 15 tiles.
The relative energy throughput of the tiles is shown in fig. 3. Each tile is represented by an operating point which is located in the centre of a tile.

Modelling the PMSM
On the one hand the model of the permanent magnet machine must be developed such that it needs little computational effort. On the other hand it must include all physical effects which are playing an important role in the machine design. The model of the PMSM is based on the fundamental wave phasor diagram expressed in the direct and quadrature axis coordinate system (see fig. 4). The parameters are the phase resistance R 1 , the direct and quadrature axis inductances L d and L q , and the magnet flux linkage Ψ pm . The components of the current vector are the direct-and quadrature-axis current.
The components of the voltage vector are the direct-and quadrature-axis voltages.
Accounting for the optimization process, the model is chosen such that the dependency of the torque T shaft , the phase voltage U and the total losses P loss on the current vector I is expressed by quadratic functions. The matrices Q loss , Q t and Q u are 2 x 2 matrices while g loss , g t and g u are 2 x 1 vectors. The matrix Q loss is positive definite. The elements of Q loss , Q t and Q u as well as g loss , g t and g u are functions of the parameters R 1 , L d , L q and Ψ pm (see section 8.1). The model is extended by a simple iron loss model which is based on the Steinmetz equation.
The parameter σ hyst is the hysteresis loss coefficient and σ ed is the eddy current loss coefficient. The frequency f is the electric frequency and B is the local flux density. Normally the dependency of the hysteresis losses on the flux density is not quadratic. Typical exponents of B in the hysteresis loss term are in the range from 1.5 to 1.8. Since the model should have a quadratic behaviour an exponent of 2 is chosen. The iron losses in the stator tooth P iron,t and rotor yoke P iron,y are treated separately and are derived from (13).

Basic layout
To determine the connection between the parameters of the phasor diagram and the geometric parameters of a lamination an initial lamination layout must be chosen by the designer. The chosen lamination layout is shown in fig. 5. This lamination allows an analytical calculation of the parameters R 1 , L d , L q and Ψ pm , which is important for the processing time. This choice does not limit the results of the optimization to the chosen lamination as described in section 6. Typically the magnetic characteristic of a lamination does not depend on absolute values but on a set of ratios of geometric parameters. Therefore, the following ratios are defined: The parameter τ s is the slot pitch. The lamination parameters are collected in the parameter vector Γ : The dependencies of R 1 , L d , L q and Ψ pm on Γ are expressed by analytic functions as far as possible. This allows calculating the gradients of these functions. They are used in the optimization process.

Permanent magnet flux linkage
For an analytical calculation of the permanent magnet flux linkage a simplified model is used.
The relative permeability of all iron parts is set to infinity. The effective magnetic air gap ′ δ is set to a constant value. Considering these simplifications, the permanent magnet flux linkage can be calculated in different ways. All presented possibilities are based on the radial component of the magnet flux density in the air gap. The fundamental wave form of the air gap flux density is (17) The angular position on the rotor is represented by α. According to the coil pitch the air gap flux density is integrated along the circumference to evaluate the permanent magnet flux linkage. For this purpose the peak value B δ must be known. The easiest way to calculate the air gap flux density is to assume a rectangular dependency of the flux density on the rotor position. The flux density reaches its maximum value in the magnet area and is set to zero in the magnet gap. In this case the peak value of the fundamental flux density distribution caused by the magnets can be calculated by using the load line.
The magnet is specified by its remanent flux density B r and the recoil permeability µ m . The second possibility is solving the Laplace equation. The solution is derived similar to [11]. First the magnetization of the permanent magnets is described by a Fourier series.
The fundamental component of the flux density in the middle of the air gap is given by The number of pole pairs p must exceed one. The dependency of the factor A l on Γ is given in the appendix.
Since the magnetic air gap is not constant, assuming a constant air gap leads to an error in the fundamental waveform of the magnet flux linkage. Therefore the peak value of the air gap flux density for a non constant air gap is compared to ,Laplacê B δ and ,load linê B δ . The reference air gap flux density is calculated by the finite element method. A slotless stator is assumed and the relative permeability in the soft iron parts is set to infinity. Figure 6 shows the comparison. Since the fundamental waveform of the solution of the Laplace equation is closer to the finite element waveform this model is chosen.

Inductances
The synchronous inductance is the ratio of the phase flux linkage divided by the phase current when all phases are powered by a symmetric three phase current. The field of the magnets is removed by setting its coercive field strength to zero.
To find an analytical model for the direct-axis inductance a reference air gap flux density is computed by using the finite element method. The finite element model consists of a slotless stator and a sine distributed winding. The relative permeability of all iron parts is set to infinity. The winding is powered by a direct-axis current. Figure 7 shows the air gap flux density in the middle of the air gap (black line).
In the magnet gap regions the magnetic effective air gap is smaller than in the magnet region.
2 q δ = δ + δ (24) By using the rotating field theory [9] the fundamental air gap flux density for each air gap can be computed. It is assumed, that the air gap flux density for δ 1 is valid in the magnet region and the air gap flux density for δ 2 is valid in the magnet gap region. The blue line in fig. 7 shows the resulting air gap flux density. This leads to the definition of two inductances: The air gap direct-axis inductance is given by a weighted sum of L h1 and L h2 .
( ) The stack length is referred to as l fe and w 1 is the number of turns per phase. The parameter ξ p is the fundamental winding factor. The air gap quadrature-axis inductance is calculated in a similar manner ( ) 1 l i α = − α (33) Besides the air gap direct-and quadrature-axis inductances the slot leakage inductance L σs and the differential inductance due to a non sinusoidal winding distribution are also included in the model. The differential leakage inductance L σd is modelled by a leakage coefficient σ d which is derived from the Görges diagram [10].

Saturation Model
Here only the saturation due to the permanent magnet field is taken into account.

Model validation
The results of the analytical model are compared to a static finite element analysis (FEA). The comparison is done for the operating point T op = 52 Nm and n op = 3800 rpm and the optimized lamination layout given in table 2. The nonlinear behaviour of the electric steel is computed by using (40) and included into the finite element analysis. For calculating the inductances the fixed permeability method described in [2] and [5] is used. The comparison is shown in table 1. The values given for the finite element calculation are mean values.

Reduced parameters and per unit system
The goal of the presented design process is finding suitable values for the parameter vector Γ which minimizes the total losses for a given drive cycle. By Γ a lamination design is defined. A lamination is two dimensional and different windings can be inserted. The number of turns per phase or the stack length does not influence the flux densities in the magnetic circuit of the electric machine. The characteristic of the magnetic circuit is controlled by the lamination design. To account for this fact a reduced parameter set is introduced: In a second step the model is normalized by defining a reference current I b , a reference flux linkage Ψ b , a reference radius r B and a reference iron loss parameter k ub . The iron loss parameter describes the dependency of the iron losses on the frequency and the parameter vector Γ (see section 8.1). The resulting normalized parameters are marked by a Π.
The optimization is carried out for the normalized system. All operating points have also to be normalized.

Optimization process
The optimization process consists of an inner and outer optimization step.

Inner optimization step
The inner optimization is carried out for each operating point. Here the parameter vector Γ is constant. It implements a minimum loss strategy. The optimization variables are the components of the current vector I .   (10), (11) and (12)).

( )
The Jacobi matrix J t , i is calculated numerically. Equation (58) offers the possibility to measure the error of the numerical calculation. The gradient of the shaft torque is calculated by (59). If the absolute value of the result is larger than a small threshold value the results are assumed to be incorrect.
Since the outer optimization is based on the results of the inner optimization it is important that the inner optimization is reliable and fast. For this reason the model of the permanent magnet synchronous machine is limited to quadratic models with respect to the current vector I .

Outer optimization step
The optimization variables of the outer optimization steps are the components of the parameter vector Γ . The target function is the weighted sum of losses of the operating points. The gradient of the target function with respect to Γ is calculated by using (57). In this paper the operating points used in the optimization correspond to the tiles centres as described in section 2. The weights are the relative energy throughput of each tile. The formulation of the outer optimization problem is: Now the current limit is added to the optimization process. The voltage inequality constrained appears in the inner and outer optimization loop. Additional constraints can be added to ensure that the elements of the parameter vector Γ stay inside a reasonable range.

Optimization algorithm
The inner and outer optimization uses the same nonlinear optimization algorithm. The algorithm was described in [6]. It is based on an exact penalty function and belongs to the class of trust region optimization algorithms [4]. By using an exact penalty function the constrained optimization is transferred into a non constraint optimization. Therefore the objective and the constraints are combined into a new objective. The optimization algorithm uses gradients of the objective and constraints. The gradients are used to build a quadratic model of the exact penalty function around the current iteration. It is assumed that this model is valid inside a given trust region radius. This radius is controlled by the algorithm. The model is used to forecast a reduction in the penalty function when the parameter vector Γ is updated to a new value. If the predicted reduction is not close to the decrease of the penalty function the trust region radius is reduced. If the model and the penalty function agree well, the trust region radius is enlarged. More detailed explanations on trust region optimization algorithms can be found in [4]. Using an exact penalty function allows initial values of the parameter vector Γ which do not satisfy the constraints. For example the initial parameter vector could be chosen such that the voltage constraint (64) is not achieved for all operating points at the beginning of the optimization. On the other hand it must be ensured that the values of the parameter vector Γ stay inside a meaningful range and the model can be validated.

Results
The optimization has been carried out for 80 operating points. Each operating point corresponds to a tile centre as described in section 2. Their position in the torque speed plane is marked by a red cross in fig. 10. The points are chosen according to their energy throughput. By using 80 points 85% of the total energy throughput in motoring operation is covered. The average torque is 39.8 Nm and the average rotational speed is 3935 rpm. The average values are calculated by using the relative energy throughputs of each tile (see equation (65)). The initial lamination and the final layout after the optimization are shown in fig. 11. For both layouts the weighted sum of copper losses P cu,cycle and iron losses P fe,cycle are calculated (see table 2).
The minimization of the normalized losses during the optimization process is shown in fig. 12. The overall cycle losses are reduced by 11.8% of the initial losses. The average torque of the cycle is approximately half of the maximum desired torque. Hence the copper losses are the dominant losses. To reduce these losses the optimization algorithm increases the slot area. This results in a small stator tooth width and high flux densities in the stator tooth. Because of this it is important to include saturation into the model of the permanent magnet Iteration k normalized cycle loss Figure 12: Reduction of the normalized total cycle losses synchronous machine. A lower limit of the stator tooth is given by saturation effects. Since the flux density in the stator tooth is limited by the magnetization curve of the electric steel a very small tooth width would lead to a reduced permanent magnet flux and this would increase the copper losses.
As it can be seen from the load line equation (18) and (14) a small value of α δhm leads to high air gap flux densities and to an increased permanent magnet flux linkage. This also decreases the copper losses. A small value of α δhm results in a large magnet height h m , especially if there is a lower limit for the air gap. An important design goal is to minimize the used magnet volume. Therefore an lower limit α dhm,min = 0.3 is included in the optimization. This limit is chosen such that there is no danger of demagnetization of the magnets. Due to manufacturing tolerances a lower air gap limit is also introduced. Additionally for surface mounted magnets a rotor tape may be used which requires a minimum air gap.

Further steps
In the previous sections it has been described how the distribution of operating points in the torque speed plane can be matched to the parameters of a lamination. Due to processing time a specialized model was used which is based on the lamination shown in fig. 5. This lamination layout may not be the desired final layout. This leads to the question how the results can be transferred to a different lamination. Most magnetic circuits of permanent magnet synchronous machines can be described by the parameters L d , L q and Ψ pm . Suitable values for these parameters have been found by the optimization described in the previous section. These values are referred to as L d1 , L q1 and Ψ pm1 . Their values are directly related to the given drive cycle. The permanent magnet flux linkage and the inductances depend on the parameters ′ Γ of the desired lamination and the current vector I . For the ith operating point three error functions could be defined: This allows defining an objective for a second optimization.

Conclusion
In the previous sections an optimization process is described which allows to find optimal parameter values of a lamination layout for a given drive cycle. In contrast to the common methods of machine design for variable loads a large number of operating points can be included into the optimization. This is achieved by using a specialized model of the permanent magnet synchronous machine in combination with an optimization strategy. In the example given in this paper the overall cycle losses are reduced by 11.8 % in comparison to the initial losses.