Minimum Set of Rotor Parameters for Synchronous Reluctance Machine and Improved Optimization Convergence via Forced Rotor Barrier Feasibility †

for actuation of Electric Multi-purpose Take-Off” International Conference Electrical Machines (ICEM). Abstract: Although rare earth materials are the critical component in high torque density permanent magnet machines, their use has historically been a commercial risk. The alternatives that have been in the recent industry focus are synchronous reluctance machines (SyRM). They have lower torque density but also relatively low material cost and higher overload capability. Multi-layer IPM and SyRM machines have signiﬁcant geometric complexity, resulting in a high number of parameters. Considering that modern machine design requires the use of optimization algorithms with computational load proportional to the number of parameters, the whole design process can take several days. This paper presents novel SyRM parameterization with reduced number of parameters. Furthermore, the paper introduces the novel forced feasibility concept, applied on rotor barrier parameters, resulting in improved optimization convergence with overall optimization time reduced by 12.3%. Proposed approaches were demonstrated using optimization procedure based on the existing differential evolution algorithm (DE) framework.


Introduction
To reduce environmental impact, global legislation is pushing to increase electric vehicle (EV) production [1]. In addition to regulatory requirements, consumers are demanding cleaner and safer vehicles. This has led to a tectonic shift in the automotive industry, both in terms of knowledge and production, forcing the industry to evolve rapidly.
At the moment, passenger vehicles are leading the market development (Tesla, Toyota, BMW), mainly because of lighter vehicles that require smaller traction batteries. On the other hand, commercial vehicles are much heavier and require large battery capacity, resulting in significant production costs. Therefore, the commercial vehicle industry is forced to move into niche markets, such as medium-duty, short-haul, and last mile applications. Long-haul vehicle development is likely to take more time and have even more sensitive financing.
Examples of commercial vehicles suitable for SyRM adoption include electric multipurpose vehicles (eMPVs), such as refuse trucks, hook loader trucks, or vacuum trucks [2]. eMPVs must actuate additional body systems (usually through some type of hydraulic pump) in addition to electric propulsion. Traditionally, this actuation is done by a diesel engine or a gearbox-mounted output shaft referred to as power take-off (PTO). Considering price, overload capability, and production simplicity, SyRM can be a viable ePTO solution [2,3]. This paper presents novel SyRM parameterization with a reduced number of parameters (rotor nomenclature according to Figure 1). During the optimization process, infeasible models may occur. Depending on the designer's choice, the infeasible models can be rejected or modified until they become feasible. The idea is that the procedure does not discard randomly generated infeasible geometry but modifies it until feasibility is achieved (hence the name "forced" feasibility).

Peak Performance Requirements
The machine requirements in this paper are derived from Reference [3]. All peak performance requirements are listed in Table 1.

Optimization Method
Most of the requirements for the design of electrical machines are in conflict with each other (reduction of volume or mass, increase of efficiency, etc.). This is evident in the problem of increasing efficiency [4][5][6] through global legislative initiatives [1]. For traction drives, high efficiency in limited packaging space is an absolute imperative [7]. Therefore, a manual design that satisfies all constraints can be an overwhelming task due to a large number of coupled parameters that affect the performance and quality of the machine.
According to Pellegrino [8,9], the computational load increases proportionally with the number of geometric parameters. This is inherently the case for IPM and SyRM machines, leading to a high number of optimization variables and a longer optimization time.
Today, optimization algorithms enjoy great popularity among designers of electrical machines [10][11][12][13][14][15][16]. The personal experience of the designer should not be underestimated, but, due to the non-linearity and complexity of the relationships between the geometry of electrical machines and their performance, it is generally believed that only mathematical optimization can push the boundaries to better designs.
Optimization algorithms can be divided to gradient based methods and stochastic or metaheuristic methods (PyOpt provides several open-source algorithms [17]). Gradient type methods converge fast but have difficulties with global optima. Usually they require feasible starting point which can be a problematic task in complex problems (Quasi Newton method [18]). Stochastic methods are heavily used in electrical machine optimization (Pow-ell's method [18], Nelder-Mead method [19]). The disadvantage is that the convergence can last for days, and global optimum cannot be mathematically proven. On the other hand, from engineer's point of view, these methods can find a satisfying global result. Popular metaheuristic methods are based on natural behavior (Evolutionary algorithm [20], Differential evolution [21], Particle Swarm [22]).
All methods are generally set to solve a single or multi-objective problem. The goal of design optimization is to have a chosen objective function f ( x) reach its minimum or maximum value while keeping other engineering indices within an acceptable range [23].
The use of finite element analysis (FEA) is inevitable in the case of SyRM's because saturation in the rotor bridges and posts significantly affects the final performance. FEA is computationally intensive and optimization can require thousands of calculations trough generations. Significant time savings can be achieved if all calculations are performed using magnetostatic simulations with fixed rotor position. Detailed explanation of different approaches for calculation of IPM machine parameters and performance using only magnetostatic simulations is available in Reference [24]. This paper uses an improved version of DE algorithm proposed by Žarko et al. [25] based on Lampinen's constraint function approach [21,26,27].

Preset Model
The number of slots and poles is chosen to be 36/4 with 4 rotor flux barriers, resulting in a two-layer integer slot winding with distributed overlapping coils. This combination provides a good compromise between the inherent ability to mitigate torque pulsations, susceptibility to noise, and the ability to use multiple parallel paths.
The ideal number of turns per coil (N c ) and parallel paths (a p ) for matching the base speed is automatically calculated based on winding feasibility [28] and ultra-fast scaling laws [29].
The goal of this paper is to prove that the forced feasibility approach (which will be discussed in later sections) yields a shorter optimization time. Considering that the selected machine has a relatively large number of parameters, all non-rotor parameters are taken from the optimized machine design in Reference [3] (Table 2). In addition, the initial design is constrained by peak performance requirements at the base speed (Table 1).
A set of parameters which are subject to optimization are listed in Table 3; colors and numbering in Tables 2 and 3

Automated Geometry Construction
Any optimization requires automatic design generation. As mentioned earlier, a higher number of optimization variables is associated with a longer optimization time. Gamba et al. [30] have shown that three variables per barrier (3 · k, where k is the number of barriers) is the appropriate number of parameters to use for a fast yet accurate description of multi-barrier SyRM. To reduce this even further, the 2 + (2 · k − 1) alternative is proposed in this paper. Several rotor barrier types are commonly used in SyRM design: circular, hyperbolic, fluid (Zhukovsky), segmented, etc. (open-source SyRE project offers more details and instructions on geometry generation [31]).
Within this paper, barrier line profiles ( Figure 2) are derived from conformal mapping theory and the Zhukovsky airflow potential formulation [30,31]. This was originally developed to describe the flow paths of fluids channeled by two infinite plates forming an angle π/p, and a plug centered at the origin of the reference frame (p is the number of pole pairs). In the solid rotor context, the plug represents the non-magnetic shaft with a radius of D sh /2. Equations (1) and (2) express the magnetic field potential lines in parametric form. (1)

Single Barrier Construction
The first step in creating the flux barrier is to select the point E k in polar coordinates. E k consists of two components, the radial component r k = (D b − 2δ)/2 − w bb , which has a fixed value, and the variable angle ϑ k (5). The centerline parameter C k is computed by solving the Equation (6). Virtual center barrier line is then computed by solving Equation (7). It is important to note that the angle vector should be generated in the range ϑ ∈ ϑ k − π 2p , ϑ k ; otherwise, if the angle is close to π 2p , the radial component will be weighted to infinity, leading to a computational error. The point G k (8) is the center barrier coordinate (always lies at angle π 2p ) and a reference point for calculating the inner and outer barrier line.
The barrier is constructed from a virtual centerline, which is modified by adding offsets ∆r in and ∆r out to form inner and outer flux lines ( Figure 3). Offsets in millimeters are calculated from per-unitized input offset parameters according to expressions (4), (3), where d k in and d k out represent distance between G k points according to Figure 4. [mm] [mm] The next step in the construction of the barrier is to compute C k out (9) and C k in (11), which completely define the equations of the inner and outer barriers. Solving the Equation (2) = r k with the arguments C k in , C k out gives the intersection point angles ϑ k in , ϑ k out . The last step is the calculation of barrier lines r out (10) and r in (12) over the given angles.

Construction of All Rotor Flux Barriers
The previous section described the flux barrier construction based on three parameters: variable angle ϑ k , inner (∆r in ) and outer (∆r out ) barrier offsets. If all barriers were constructed according to this principle, the total parameter number would be 3 · k. Considering that center virtual barrier line is not part of the final barrier, there is a way to reduce the total number of parameters. Two angles, ϑ min and ϑ max , indicate the allowed center barrier line range. Depending on the number of barriers k, they define angles ϑ 1..k with equidistant angular offsets ∆ϑ r defining points E 1..k used to construct each virtual center barrier line ( Figure 4). To make the procedure robust, and to avoid clashes with the shaft, final barrier (closest to the shaft) inner offset value is always zero. This procedure reduces the total number of barrier parameters to 2 + (2 · k − 1). The constant 2 refers to the parameters ϑ min and ϑ max , −1 refers to the constant ∆r final in = 0, which can be removed from the total number of parameters.
The choice of equidistant virtual flux line angle offset (∆ϑ r ) might appear to be a design-limiting constraint, but this is not the case. The virtual lines are used only as a reference for calculating the inner and outer barriers. Only in the special case, they represent flux barrier middle line when ∀∆r k out = ∀∆r k in . Optimization results proved that ∀∆r k out = ∀∆r k in (Table 4). Furthermore, ∀∆r k out = 0 and ∀∆r k in = 0 (except in the final barrier, where ∆r final in = 0 by default). This proves that the assumption of equidistant virtual line offset angles has no effect on the final rotor design and can be used in the case of the Zhukovsky flux barrier type.

Optimization Procedure Details
The optimization of the 2D cross-section is set up as a single-objective problem mathematically defined as: Find the vector of parameters (13), subject to D parameter boundary constraints (14) and subject to m inequality constraint functions (15), which will maximize objective function (16).
The applied optimization workflow is shown in Figure 5. The optimization process starts with the problem definition (boundaries, constraints, objectives, etc.) and a preset of constant parameters (slots, poles, active diameter, etc.).
After entering the optimization loop, the following steps are performed iteratively: 1. the optimization algorithm generates the vector x (optimization variables); 2.
variables are converted to model parameters; 3.
model is generated based on the model parameters; 4. model is solved; 5.
performance is extracted from the solution (values for constraints and objective functions are calculated from the solution); 6.
data is passed to the optimization algorithm.
Single-objective optimization has only one objective function, which may have multiple variables and scaling factors. When the optimization converges, the value of the objective function saturates around a certain value. Without stopping criteria, the objective function will improve in infinitesimal steps for a very long time without any significant design improvement. In this paper, the following stopping criteria are used: after each generation, the last four objective function increments f 1..4 are checked; if the difference between f 1 and f 4 is less than 0.2%, the optimization is stopped; otherwise, it continues until the 160th generation.

Model and Solution Feasibility
In optimization problems, the notion of feasibility is related to the acceptance criteria of the solution. A solution is declared feasible if it satisfies certain criteria (15). The feasibility of a solution implies that a solution exists, i.e., the problem or model that provides this solution is solvable. During the optimization process, models may occur that are not solvable. These models have no solution and can be declared infeasible without solving them, which can decrease computational burden in the case of FEA-based solver. The models that are not solvable are declared as infeasible considering the model feasibility criteria.

Geometrical Feasibility
A special subset of model feasibility criteria is geometric feasibility, which characterizes whether the model geometry satisfies certain criteria (e.g., mechanical integrity, physicality, overlaps, etc.). In some cases, trying to solve a geometrically infeasible model with external solvers can lead to a crash of the whole optimization routine, a problem that highlights the need for detection methods to deal with geometrically infeasible models (e.g., barrier 1 (blue) collides with barrier 2 (yellow), the overlap is marked in red, Figure 6c).
To avoid the construction of such an invalid model, a geometric feasibility procedure can be performed within the optimization algorithm for each candidate vector. Regardless of the method of handling and detecting geometric infeasibility, it is always beneficial to reduce the occurrence of infeasible models. The simplest method for reducing the occurrence of geometrically infeasible models is to introduce lower and upper parameter bounds, which can be in the form of linear (Figure 7a), non-linear function bounds (Figure 7b), or complex bounds (Figure 7c).
For simple problems, the introduction of bounds can completely avoid the occurrence of infeasible models, while for complex problems it reduces the probability of the occurrence of the infeasible candidates. For this reason, optimization algorithms must include a method for dealing with geometrically infeasible candidates.
In general, a geometrically infeasible candidate is discarded while the new candidate takes its place. To generate a replacement candidate, Žarko et al. [25,32] randomly initialize the entire parameter set (13) until a geometrically feasible replacement candidate appears. The drawback of this method is a possible rejection of candidates with some good properties. Moreover, this method may lead to slow convergence to the optimal solution if the optimal candidate is on the boundary of the feasible space.
The alternative approach is forced feasibility, where each infeasible design is subjected to parameter modification until feasibility is achieved (i.e., barrier 1 (blue) and barrier 2 (yellow) are modified until the candidate achieves the specified flux carrier width w goal , Figure 6d). This approach requires smart parametrization with minimum feasibility con-straints and can potentially be extremely complex. On the other hand, potential benefits include reduced optimization time (no need to wait for a random feasible design to emerge) and faster convergence to the final result. The substitution of infeasible candidates is implemented in the form of a projection onto the feasible space (Figure 7d).
Projection operator can be mathematically written as: where x orig represents parameter vector of original, geometrically infeasible candidate, x new is its geometrically feasible replacement/alternative, and Q is weighting matrix (18) in which its coefficients control the projection path (1, 2, and 3 in Figure 7d). If all coefficients satisfy q 1..D = 1, projection is orthogonal, and the new geometrically feasible candidate is closest to the original infeasible candidate.

Forced Feasibility Algorithm
For simplicity, the functionality is explained using two barrier SyRM (Figure 6b-d). The first step is geometry generation based on the parameterization approach from Section 3 (Algorithm 1: ln: 1-5). The next step is to compute the variables d k−1,k , w k−1,k , w goal ( Figure 6b) and check whether the infeasibility condition is satisfied (Algorithm 1: ln: 11). The design may be infeasible for one of two reasons: there is a barrier conflict (w k−1,k ≤ 0) or the flux carrier width is too small (w k−1,k ≤ w goal ). If not feasible (Figure 6c), proceed to forcing feasibility (Algorithm 2).
The purpose of the force feasibility function (Algorithm 2: ln: 1) is to provide the current position information of the barrier (w goal , d k−1,k , ∆r k out , ∆r k−1 in , Figure 6b) and provide minimally modified barrier offsets (∆r k out , ∆r k−1 in ), defining a new feasible design (Figure 6d). The minimal deviation is secured via MATLAB function FMINCON function [33]. FMINCON default input parameters are listed in Algorithm 2: ln: 2-7, where x c represents the current barrier offset vector.
During the search, FMINCON iteratively calls COSTFCN, which is responsible for the convergence of the search (calculates the deviation of the generated x and the initial offset vector x c ), and CONSTRAINTFCN which calculates the deviation of the generated flux carrier width (w) from the carrier width (w goal ). When the algorithm converges, FMINCON returns ∆r k out , ∆r k−1 in , and a new feasible geometry is generated (Figure 6d). In summary, when the system detects an infeasible case, such as in the example shown in Figure 6b, where 1st and 2nd barriers overlap, ∆r 2out and ∆r 1in are iteratively modified until flux carrier width w goal is reached (Figure 6d). This prevents the generation of too thin or too wide barrier geometries and improves the optimization procedure. The flux carrier width is randomly generated according to the equation w goal = d k−1,k /R w k in the range R w k = [1.561, 1.88] (Figure 6a). The range is empirically derived from several optimized designs.
Calculate flux carrier width 9: Calculate minimal flux carrier width 10: w goal = d k−1,k /R w k

11:
if w k−1,k ≤ 0 or w k−1,k ≤ w goal then Barriers are not feasible, force feasibility 12: function FORCEFEASIBILITY( w goal , d k−1,k , ∆r k out , ∆r k−1 in ) 13: . . . 14: return ∆r k out , ∆r k−1 in x 0 = 0.5 · ub Initial guess 7: x c = [∆r k out ∆r k−1 in ] Minimization goal Find a solution which satisfies constraints and minimally changes input parameters x c via FMINCON function 8: function FMINCON(COSTFCN, x 0 , B, A eq , B eq , lb, ub,CONSTRAINTFCN, w goal , d k−1,k , x c , opt, k) FMINCON iteratively calls COSTFCN, CONSTRAINTFCN and returns values which satisfy the constraint(s) with minimal deviation from current offset vector x c 9: . . . 10: return ∆r k out , ∆r k−1 in 11: COSTTFCN is responsible for result search convergence 12: function CONSTRAINTFCN(x, w goal , d k−1,k ) CONSTRAINTFCN calculates the deviation of generated flux carrier width from the goal width 16: x(1) → ∆r k out 17: x(2) → ∆r k−1 in 18: return P 22: return ∆r k out , ∆r k−1 in Overall, this system does not discard randomly generated infeasible geometry, but modifies it until feasibility is achieved (hence the name "forced" feasibility). A similar procedure can be implemented for any problem that can be described by smooth analytic functions implemented as user-defined CONSTRAINTFCN.

Handling of Inequality Constraints
Inequality constraints usually arise from various electromagnetic, thermal, mechanical, manufacturing, economic, or normative limits, such as maximum flux density in the stator tooth, maximum magnet temperature, maximum stress in the rotor bridge, minimum magnet dimensions, maximum active material cost, maximum noise, etc. [32].
The traditional approach to constraint handling uses penalty functions to penalize the solutions that violate the constraints. This principle is implemented in terms of weighted sums that modify each objective function. Despite the popularity of penalty functions, they have several drawbacks. The most important is the need for careful fine-tuning of the penalty factors responsible for efficiently approximating the feasible range. Moreover, this method may suffer from problems related to poor choice of weighting factors, which may affect convergence. In this paper, we use an improved constraint function algorithm developed by Žarko et al. [34].
Inequality constraints for this particular case are defined in Table 5. The constraint function g 1 checks rotor structural factor of safety at maximum over-speed (1.2 · n max ). Magnetostatic torque at n b T static ≥600 Nm g 3 Flux density in stator yoke B sy,max ≤1.6 T g 4 Flux density in stator tooth B st,max ≤1.8 T g 5 Power factor cos ϕ min >0.6 g 6 Torque ripple with skewing T ripple,max ≤10 % The procedure related to constraint function g 2 contains several subfunctions designed according to ultra-fast scaling laws [29]. Multiple magnetostatic FEA calculations are performed to find the maximum torque versus current phase advance curve and to determine the optimal maximum torque-per-ampere (MTPA) control angle by polynomial fitting (the input machine has one turn per coil and one parallel path). The number of turns per coil and the number parallel paths of the machine is then matched to the required base speed. The optimization initially assumes a fixed stack length ( Table 2, parameter l s max ). If the calculated torque at l s max is smaller than required, the machine design does not satisfy the constraint. Finally, it is checked whether the stator phase current is smaller than the maximum inverter current; otherwise, the constraint is not satisfied. After this step, all magnetostatic calculations are completed. The results are extracted and evaluated in the following constraint functions. Constraint g 3 checks the maximum stator yoke flux density, while g 4 checks the maximum tooth flux density. The purpose is to penalize the high saturation and the designs with increased iron losses.
Finally, a transient FEA calculation is performed at base speed to determine the power factor, average torque, torque ripple, and terminal voltage. Constraint g 5 checks the power factor. The primary goal in this ePTO optimization case study was maximization of the average torque, while torque ripple was of secondary importance, therefore being handled in the constraint (g 6 ). As a torque-ripple mitigation option, a rotor skew was selected. Skewing angle is one stator slot or 360 • /36 = 10 • mechanical degrees [35]. Without loss of generality, other more affordable, or more practical, torque-ripple mitigation techniques can be applied to this problem [36,37]. In addition, torque ripple can be included as another optimization objective so that the design trade-off is made on the Pareto front of torque versus torque ripple.
If all constraints are satisfied, the final step is to compute the objective function f (19). η is the efficiency, and T trans,skew is the transient torque with applied rotor skew. The constants 700 and 0.96 represent scaling coefficients, which are used for combining torque and efficiency within a single objective function. Furthermore, these values are important for proper optimization convergence and were intentionally chosen to be larger than the peak torque and efficiency at the base speed to constrain the objective function to values ≤ 2 (important for the final comparison of the optimization convergence of forced feasibility and randomly generated geometries).

Optimization Results
Five consecutive optimization runs of the DE algorithm with population size NP = 24 of geometries generated both randomly and with forced feasibility were performed (results in Tables 4 and 6).
To reduce optimization time and compare both approaches, most of the design parameters were taken from a previously optimized design [3]. Twenty-eight parameters were frozen (Table 2), and only 7 parameters were used for optimization (Table 3). This trade-off yields the same average objective function result (difference is +0.05%).
The average number of generations required to achieve convergence in the random generation case is 115.6, and 95.6 in the forced feasibility case, which is a reduction of 17.3%. In addition, in terms of computation time, the average convergence of forced feasibility is 3.34 h shorter (12.3% reduction).
When comparing the average torque results (T trans,skew ), both approaches yield practically the same outcome ( Table 4). The remaining results (η, B sy,max , B st,max , cos ϕ, T ripple ) are identical for forced feasibility and random generation. This is expected for two reasons: most of the parameters were taken from Reference [3], and the optimization algorithm is the same. Finally, the identical results confirm that forced feasibility does not affect negatively the optimization outcome.  Table 4 summarizes the optimized rotor parameters for both optimization approaches. As mentioned earlier, the most important result is the fact that all parameters are greater than zero (neither the inner nor the outer blocking line sticks to the virtual centerline) and ∀∆r k out = ∀∆r k in . All optimized cross sections are shown in Figure 8, which confirms that the position of the centerline of the virtual barriers does not affect the optimization result.

Conclusions
This paper demonstrates a novel automatic design procedure of SyRM rotor with minimum number of geometric parameters, which simplifies the design generation and reduces the optimization time. The presented procedure is implemented into an existing single objective DE optimization algorithm framework which interrupts evaluation of constraint functions when the inequality constraint is violated, thus saving computation time.
The second paper contribution is the introduction of novel forced feasibility concept which improves optimization convergence proved by successive comparative optimization runs with randomly generated rotor barrier geometries. The results show that properly implemented forced feasibility leads to a further reduction in optimization time (12.3% shorter).
The machine design originally presented in Reference [3] has 21 optimization variables. The entire optimization process took 7 days. Without the forced feasibility method, the process would be 12.3% longer (approximately 1 additional day). Considering that the optimization time is proportional to the number of parameters (i.e., a two-layer V-shape PM machine may have 44 parameters), it can be concluded that the total calculation time can be significantly reduced by using forced feasibility approach on different machine topologies.
Author Contributions: This paper is a continuation of the paper "Design and Optimization of Synchronous Reluctance Machine for actuation of Electric Multi-purpose Vehicle Power Take-Off" submitted for ICEM 2020 conference. B.B. created a preset model and automated geometry scripting with a minimal set of parameters. The original multi-objective optimization procedure was adapted to support single-objective optimization, which is easier to evaluate in terms of forced feasibility impact on optimization time reduction. B.B. adapted the original code developed by S.S., simulated the results, and prepared all tables and figures. S.S. supervised the process and formulated the contributions. T.J. worked on the mathematical formulation of the forced feasibility and the practical implementation of the algorithm through the MATLAB fmincon minimization function. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: