Robust Design Optimization of Electric Machines with Isogeometric Analysis

: In electric machine design, efficient methods for the optimization of the geometry and associated parameters are essential. Nowadays, it is necessary to address the uncertainty caused by manufacturing or material tolerances. This work presents a robust optimization strategy to address uncertainty in the design of a three-phase, six-pole permanent magnet synchronous motor (PMSM). The geometry is constructed in a two-dimensional framework within MATLAB ® , employing isogeometric analysis (IGA) to enable flexible shape optimization. The main contributions of this research are twofold. First, we integrate shape optimization with parameter optimization to enhance the performance of PMSM designs. Second, we use robust optimization, which creates a min–max problem, to ensure that the motor maintains its performance when facing uncertainties. To solve this bilevel problem, we work with the maximal value functions of the lower-level maximization problems and apply a version of Danskin’s theorem for the computation of generalized derivatives. Additionally, the adjoint method is employed to efficiently solve the lower-level problems with gradient-based optimization. The paper concludes by presenting numerical results showcasing the efficacy of the proposed robust optimization framework. The results indicate that the optimized PMSM designs not only perform competitively compared to their non-robust counterparts but also show resilience to operational and manufacturing uncertainties, making them attractive for industrial applications.


Introduction
The growing importance of electric machines in transportation and other applications necessitates highly optimized designs.This trend is clearly seen in the shift towards electric vehicles in the automotive industry.This has become especially significant in efforts to mitigate climate change and promote environmental sustainability.Over the past decade, the development of electric machines has shifted from physical to virtual prototypes, driven by advancements in numerical simulation.This change has resulted in cost savings and reduced material usage, improving sustainability.To further exploit these benefits, the current trend is to move from a purely numerical simulation-based procedure to one that improves the performance through optimization methods that are able to handle many design variables.
Such optimization processes in the design of electric machines aim to balance competing objectives such as high torque capacity, operational smoothness and cost-effectiveness [1,2].However, one major challenge is uncertainties during manufacturing, which can affect the performance and efficiency in a negative way or even lead to the failure [3] of the electric machine.In order to ensure a reliable design, it is advisable to account for these uncertainties in the design optimization process using robust optimization methods that account for these variations.
The field of electric machine optimization encompasses a wide range of specialized studies.In the optimization of electric machines, there exist two main approaches: gradientbased algorithms and gradient-free methods, like particle swarm optimization (PSO).The latter are non-intrusive but require many simulation runs and are thus time-consuming [4].In this paper, we choose a gradient-based optimization approach and seek to obtain robust optima with a moderate number of optimization runs.In the field of gradient-based optimization, much research has been devoted to advancing free-form optimization using the finite element method (FEM) and variants thereof (such as isogeometric analysis (IGA)), which has led to a reduction in torque ripple [1,5].Alongside these, there has been a focus on parameter optimization with different performance metrics, such as the magnet mass [6,7].However, there is a noticeable gap in research that concurrently addresses both free-shape and parameter optimization and their combined effects on electric machine performance.Only a few studies have appeared at this intersection, using FEM for their investigations [8].In contrast to these methods, our study [9] introduces a gradient-based optimization framework within the IGA context, combining both parameter and shape optimization.This approach, supported by numerical evidence, outperforms the sequential optimization, which first addresses the parameters and then the shape.Additionally, the discretization technique employed in [9] provides several benefits, particularly for free-form optimization [10].
The objective of this paper is to build on these results and extend our previously established gradient-based optimization methodology, which combines parameter and free-shape design variables in IGA, by introducing a robustoptimization formulation.We deal with uncertainties using the worst-case optimization method originally proposed by Ben-Tal [11].In the literature, the worst-case approach for robustification against uncertainties has been already successfully applied for a three-phase, six-pole PMSM, where both uncertainties and optimization variables are parameters [3].Notably, these studies employed the standard FEM for their discretization and only addressed parameter optimization; they did not consider free-form optimization.In contrast, the present study generalizes this robust optimization approach by allowing a wider set of design variables on an electric motor that is simulated with state of the art techniques to handle rotation, i.e., mortaring [12].
A novel aspect of this research is the inclusion of uncertainties in both the motor parameters and its shape, and the consideration of both as variables in the optimization process.To improve the efficiency, we utilize IGA instead of the conventional FEM in discretizing the PDE.The robust approach results in a min-max optimization problem, where objective function and critical constraints are replaced by the maximum over an uncertainty set.To solve this, we employ a generalization of Danskin's theorem [13], which yields an efficient tool to compute generalized derivatives of the resulting robustified objective function and constraints under mild assumptions.The approach's advantage is that it only requires first derivatives.However, it is crucial to note that the outer minimization problem can become non-smooth.This property necessitates the use of a non-smooth optimization solver to ensure convergence to a non-smooth optimum; see also the discussion in Bertsimas et al. [14].
The rest of this paper is structured as follows.Section 2 introduces the mathematical and physical foundations with a focus on the relevant aspects for the optimization.In Section 3, the robust optimization method and its application to the optimal design of a PMSM are described.Section 4 presents numerical results for the robust shape and parameter optimization of a PMSM.Section 5 summarizes the main results and explores possible directions for future research.

Geometry Representation
This section introduces the physical model and the computational methods that we use for our optimization techniques.Figure 1 shows the motor model that will be used in this paper.The colors of the motor represent the materials used, i.e., gray for iron, blue for air, red for copper and green for the magnet, where the magnetization direction is shown by the arrows.The copper represents stranded conductor windings, where the phase and current direction are indicated.Note that the motor is modeled with two small air gaps at the rotor iron bridges to represent saturation effects and allowing linear material properties to be used.In 2005, Hughes et al. introduced IGA, which we use to discretize the physical domain and the PDE [15].Unlike the standard FEM, which uses polynomials as basis functions, IGA uses B-splines and non-uniform rational B-splines (NURBS) [16].The flexibility provided by spline functions is essential in accurately representing conic sections while maintaining computational efficiency.This feature is particularly useful in simulating rotating machines.The formulation of the univariate B-spline basis functions N p i of degree p is expressed by a knot vector:

M H α
where the knots ξ i define the support of the basis functions.These functions are computed through a recursive process, utilizing Cox-de Boor's recursion formula [17], which is given by These basis functions have useful mathematical properties, such as the partition of unity, compact support of the basis function and pointwise positivity, which make them suitable for approximating geometries.A comprehensive discussion of these properties can be found in the work of Hughes et al. [15].The two-dimensional domain in our work is constructed using the tensor-product method, which allows for the creation of multivariate B-splines or NURBS.This method involves the multiplication of individual basis functions from each parametric coordinate, thus forming a set of multivariate basis functions.This approach provides an efficient framework and can be also extended for the construction of three-dimensional domains [15].To precisely map the parametric domain Ω ∈ (0, 1) d to the physical domain Ω ⊆ R r , where r ≥ d, we use geometric basis functions Ĝk , which can be either B-splines or NURBS.Each of these basis functions Ĝk is multiplied with the corresponding control point C k .This leads to the mapping F : Ω → Ω defined by Here, x is located within the parametric domain Ω and F(x) is the corresponding point in the physical domain Ω.The mapping defined in (2) must satisfy some regularity conditions, e.g., it must be a diffeomorphism, which is required for the application of the transformation formula for integrals [18].A visual representation of the relationship between the control points and the model's geometry is shown in Figure 2. Shifting the position of a control point C k results in a corresponding change in the geometry according to the mapping (2).Therefore, the transformation (2) can be used directly for shape optimization, which is a convenient feature of IGA.For geometries that are complex in their topology or material distribution, such as the rotor and stator in Figure 1, multipatch decomposition is necessary.This method divides the overall domain into several contiguous patches and couples the control points of neighboring patches to ensure continuity; see [20].The grey lines in Figure 1 indicate the boundaries of the images of the patches.

Governing Equations
In this subsection, we describe the physical system governing PMSMs [12].It is based on the magnetostatic approximation of Maxwell's equations.This simplification is justified and often used when electric machines operate at a low frequency, since, in this case, the displacement currents and eddy currents are negligible compared to the source current.This assumption simplifies our model to an elliptic PDE, which can be expressed in terms of the magnetic vector potential A as follows: In ( 3), the term B = ∇ × A denotes the magnetic flux density, while ν represents the material's linear reluctivity.The term J source is the density of the current source, and B r is the remanent magnetic flux density from the permanent magnets, which follows from the constitutive material law [21].For a more detailed understanding of these physical principles, see [22,23].According to Helmholtz's theorem, A is not uniquely determined by the PDE.To simplify the model, ensure its uniqueness and reduce the computational cost, it is common to reduce (3) to a two-dimensional problem by assuming a substantial axial length of the motor.This step divides the computational area into the 2D rotor and 2D stator sections, Ω R and Ω S , respectively, and leads to the elliptic PDE in 2D, where only the A z -component of the magnetic vector potential is considered.This component is separated for the rotor and stator, denoted as A z,R and A z,S .In electric machine modeling, boundary conditions play a critical role.For the Dirichlet boundary, denoted as Γ d , which is located at the outer boundaries of the stator and rotor-see Figure 1-the condition A z = 0 is imposed.On the boundary Γ ap , antiperiodic conditions are implemented due to the simulation and optimization being applied to only one pole of the entire motor.
For the coupling interface Γ ag between the rotor and stator, the following conditions must be met: where H = νB denotes the magnetic field strength, with B as the magnetic flux density, and both A z and H θ are evaluated in their respective rotor or stator local coordinate systems in dependence on the rotation angle α.This ensures the continuity of the magnetic vector potential A z and the azimuthal magnetic field strength H θ across Γ ag [24].A visualization of the the boundary conditions can be found in Figure 1.From the 2D reduction, the reduced remanence B ⊥ r is expressed as where B r is the remanence magnitude and ϕ is the magnet direction angle, defining the orientation of the magnet's remanence in the two-dimensional plane.The source current density in the k-th winding is defined as and is zero outside the k-th winding.The synchronous operation of the motor implies that the time dependency in ( 7) is represented by the rotor angle α, adjusted by the number of pole pairs P. In (7), n wind denotes the number of windings in each coil.I app represents the application current, which is the primary source of electromagnetic field generation in the system, and A coil refers to the cross-sectional area of the coil, which affects both the material usage and electromagnetic performance.The total current density J z,source is the sum over all windings, and is assumed to be evenly distributed within the coils.

Discretization
In this subsection, we focus on the IGA discretization of the governing Equation (4).To apply a Ritz-Galerkin method, we introduce ansatz and test spaces, both spanned by a basis N(x), which is defined in (1): Deriving the weak formulation of ( 4) and inserting the ansatz and test function expansions from (8), we obtain the system of equations Equation ( 9) is a coupled system of discretized elliptic PDEs.The stiffness matrices K R , K S and right-hand sides b R , b S for the rotor and stator are obtained in a standard way from the weak formulation of the elliptic PDEs (4).The coupling matrices G R and G S result from the application of harmonic (Fourier) test functions in a mortaring approach to discretize the coupling conditions (5), as detailed in [12].For distinct rotation angles α ∈ {α 0 , α 1 , . . ., α N−1 }, the stator coupling matrix G S is multiplied by a rotation matrix R α .This matrix contains block diagonal entries of sine and cosine functions, thereby eliminating the need for the reassembly of the coupling matrix for each α.This improves the computational efficiency when evaluating multiple rotation angles [24].For more details about the discretization procedure, see [9].In (9), the right-hand side and the rotation matrix depend on the rotation angle α and consequently also the solution (u R,α , u R,α , λ α ).
Taking advantage of the symmetry of the three-phase current, we consider only the angles α ∈ {0, 1, . . ., 19} degrees.In order to evaluate the performance of the system, we solve (9) and compute the torque T α (u) based on energy conservation principles [24] by Here, R ′ α is the derivative of the rotation matrix R α with respect to the rotation angle α.Since the computation is done in 2D, we include the axial length L directly in the calculation of the torque T α (u).To determine the average torque performance, we need to solve the linear system (9) for each angle α, since satisfying a predefined average torque will be a constraint in the optimization.

Optimization
In the literature on the optimization of electric machines, various performance aspects are important.We summarize the most important objective functions that are considered in [19].A major focus, especially for industrial applications, is the cost-effective use of materials.For machines using permanent magnets, reducing the size of the magnet is a critical goal, because these magnets are made of rare earth materials, which are costly and have large environmental footprints [25].A critical constraint is to keep the average torque at a certain level because the cost and environmental impact increase as more permanent magnet material is used.Another important goal is to reduce torque ripple, which helps to lower vibration and noise.These objectives have, for example, been studied in [2,5,6].However, in the optimization of electric machines, the focus has been on either parameter optimization or free-form optimization.There is less work on combining these two approaches, e.g., [8,9].However, combining the two promises better results because it offers more design freedom.
For practical design optimization problems, often, several objective functions are of interest, which leads to a multiobjective optimization problem.Here, we consider two competing objective functions and apply the weighted sum approach; see [26].The first objective function is the area of the permanent magnets A mag , contributing mainly to parameter optimization, and the second is the standard deviation T of the torque, associated with free-shape optimization.The area A mag of the permanent magnet is computed as the product of the width and the height of the region highlighted in green in Figure 1.For the torque evaluation, the standard deviation T and the mean torque T mean are calculated based on the torque (10) for a discrete set of rotation angles α by In our optimization, the mean torque T mean is constrained to be at least a predetermined minimum level T desired .This constraint affects the required size of the permanent magnet.Altogether, the combined parameter and shape optimization problem is defined as T mean ≥ T desired , (Torque Constraint) Here, m 1 and m 2 are weighting factors introduced to balance the components of the objective function.The state equation e(u, p) = 0 represents the governing equations of the system (9).p are the optimization variables, including both shape and further design parameters.In our case, p will consist of control points located on the air gap, the geometric parameters of the permanent magnet and the electric phase angle.Further explanation of the optimization variables is given in Section 4. Numerical results for another electric machine and without the consideration of uncertainty are presented in [9,27] and demonstrate that the combined parameter and shape optimization is able to reduce the torque ripple significantly while also reducing the size of the magnets.

Robust Optimization
As described in Section 1, the goal of this paper is to make the optimization problem (12) robust to uncertainties.Uncertainties are inevitable in real-world scenarios and can affect the model parameters q as well as the design variables p.In our case, we consider uncertainty in the application current I app and in the position of the control points that define the shape of the machine's air gap.Therefore, we introduce the vector w = (δp, q) ∈ R n p × R n q to represent the uncertainty in the system and apply a robust optimization approach [11].A comprehensive overview of robust nonlinear optimization can be found in [28].
We define a closed, convex and bounded uncertainty set that contains all relevant realizations of the uncertain quantities.This set accounts for uncertainties in both parameters and optimization variables.During optimization without uncertainties, the model parameters are fixed.However, the fixed values impact the optimized design.Individual sets represent uncertainties, and the combined uncertainty set is the Cartesian product of these individual sets.

•
Uncertainty set for parameters: Uncertainty set for design variables: with R q ∈ R n q ,m q , m q ≤ n q , rank(R q ) = m q .
The combined uncertainty set results from the Cartesian product: The uncertainty sets U p,2 and U q,2 are ellipsoids.R p v p and R q v q represent the allowable deviation and can be chosen by prior knowledge or based on statistical considerations.Moreover, there exist approaches based on machine learning algorithms to adaptively learn the optimal uncertainty set for robust optimization.The advantages and disadvantages of using a learned uncertainty set versus a predefined one are discussed in [29].Therefore, our goal in making the optimization problem robust to uncertainty is to ensure that the solution remains feasible and competitive with the optimized solution from (12), regardless of the uncertainties that can occur within the defined uncertainty set W w,2 .
For this, we derive the robust optimization formulation.We consider the following class of optimization problems that include problem (12).We denote by u the state, by p the optimization variables and by q the parameters.The nominal problem is given by min u,p F 0 (u, p, q) s.t.e(u, p, q) = 0, (state equation) Under the assumption that the state equation can be uniquely solved for the state u, we can eliminate u and obtain with f i (p, q) := F i (u(p, q), p, q), i ∈ I ∪ {0}, where e(u(p, q), p, q) = 0 the reduced problem See also [3].Since p, q are uncertain, we consider the following robust formulation [3,14,30,31]: where U p,2 × U q,2 is defined in (13).Compared to [30], where also the robust optimization of electric machines has been considered, the uncertainty here is not only in the model parameters but also in the design variables, including the control points.
We propose next an algorithm to solve the robust counterpart (16).To this end, it is convenient to introduce the value functions of the robustified objective function and constraints given by Then, one can rewrite the robust counterpart (16) as Since all arising functions are potentially non-convex, a global solution of ( 18) is usually intractable.
To obtain tractability, (δp, q) → f i (p + δp, q) in ( 17) is usually approximated by a linear or quadratic approximation; see, for example, [3,31,32].However, it can be shown that a linear approximation as applied in [31] cannot be used if the uncertainty extends to the optimization variables.This problem could be fixed by using a quadratic approximation, but this requires derivatives up to the third order.In the present case of PDE constraints, this is expensive to implement.
Therefore, we work directly with the problem (18), apply a non-smooth optimization algorithm and use the fact that a generalization of Danskin's theorem yields Clarke's generalized gradient of the worst-case functions φ i under quite general assumptions requiring only first derivatives [13,33].A related approach is used in [14,34], where the descent directions of the worst-case functions are used.
Theorem 1.Let U p,2 × U q,2 ⊂ R n p × R n q be compact and f i : R n p × U q,2 → R be continuous and continously differentiable with respect to p.Then, the following holds for φ i given by (17).

1.
φ i is locally Lipschitz continuous and directionally differentiable.

2.
If W * w,2 (p) denotes the set of maximizers (δp * , q * ) in ( 17) and Clarke's generalized gradient of φ i is given by

3.
If W * w,2 (p) contains a single element, then φ i (p) is differentiable at p, and (19) yields the classical gradient.
Hence, the advantage of solving (18) by using a non-smooth optimization method using generalized derivatives obtained by (19) has the advantage that only first derivatives are required.However, the computation of a maximizer (δp * , q * ) ∈ W * w,2 (p) for is required.
This motivates the following strategy.We apply a non-smooth optimization method to (18).If the algorithm requests in iteration j value φ i (p j ) and generalized derivative g j i ∈ ∂φ i (p j ) of φ i , then an approximate maximizer (δp j i , q j i ) of ( 17) is computed by applying a few iterations of an optimization method to (17) with the previous approximate maximizer (δp The proposed algorithm is summarized in Algorithm 1.

Algorithm 1 Robust Optimization Algorithm Using Danskin's Theorem
Require: Initial point p 0 , initial worst cases (δp 0 i , q 0 i ) ∈ U p × U q , i ∈ I ∪ {0}.1: Apply a non-smooth solver to (18).Initiate iteration counter j = 0. 2: while Outer stopping criterion not satisfied by p j do 3: Apply an iteration of the non-smooth solver for (18) Apply inner solver to (17) with starting point (δp until a stopping criterion is met.Result (δp j i , q j i )

7:
Return approximate values 8: Generate new outer iterate p j+1 11: Update outer counter j = j + 1 12: end while A critical step involves computing the derivatives with respect to the control points to be able to solve (17) efficiently.We provide an adjoint-based approach in Appendix A.

Numerical Results
This section presents a numerical validation of the robust optimization approach.The IGA discretization (9) of the state equation and the corresponding adjoint-based derivative computation is based the open-source package GeoPDEs [35].Harmonic mortaring between the rotor and stator domains is applied.The approximate solution of the inner maximization problems ( 17) is performed with fmincon and the solution of the non-smooth robust formulation (18) by using GRANSO (GRadient-based Algorithm for Non-Smooth Optimization) [36].
Our goal is to optimize the geometry of the rotor of an electric machine, as shown in Figure 1.In (7), the current density is I app = 9 and the winding number n wind = 12.The optimization variables p are the electric phase angle and the geometric parameters, consisting of the width (MW), height (MH) and distance from the air gap (MAG) of the permanent magnet, as shown in Figure 1, as well as five control points located at the air gap.Uncertainties are considered in the application current I app , as well as in the radial offsets of the five optimized control points.The uncertainty set for the parameter I app is onedimensional and given by 5.5% of the nominal value 9.The uncertainty set for the control points is a five-dimensional sphere with radius 0.25 mm, which is 0.56% of the rotor radius.A summary of the starting point for the optimization and the bounds on the optimization variables can be found in Table 1.The parameters required to describe the geometry of the machine can be found in ( [12], Figure A.1).The weights in the objective function are chosen as m 1 = 3000 and m 2 = 130 to balance the competing optimization goals.The torque is evaluated at discrete angles of rotation α ∈ {0, 1, . . ., 19}.The mean torque has the lower bound T desired = 1.55 Nm.The rotor is constructed from 12 distinct patches, and the stator shape is constructed from 79 patches.This results in a total of 160 control points that define the rotor's shape, and another 669 control points that define the stator's shape.We have in total 792 degrees of freedom for the rotor and 2783 for the stator.The magnetic flux density for the initial design is depcted in Figure 3. Figure 4a shows the initial design, Figure 4b the nominal optimized design for the optimization problem (12) and Figure 4c the robust optimized design that accounts for uncertainties in parameters and shape.Starting from an initial magnet geometry with width MW = 19 mm and height MH = 7 mm, the nominal optimization achieves a reduction in the magnet's material usage by approximately 26%, while the desired average torque is maintained.In addition, the torque ripple is decreased by 96%, as seen also in Figure 5.The robust design, which optimizes the worst-case behavior and satisfies the constraints for all realizations of the uncertain quantities in the uncertainty set, still shows an improvement over the initial configuration.Here, a reduction in the magnet material by about 8% is achieved, as well as a torque ripple reduction of approximately 92% (see Figure 5).This is also represented in the values of the objective function, where the nominal optimization reaches a value of 0.766 from initially 13.53, whereas the robust optimization approach reaches a value of 1.464.In addition, the robust optimization ensures the feasibility of the critical constraint, which is the lower bound on the average torque, for all realizations in the uncertainty set.On the other hand, if we compute the worst case of the average torque for the nominal optimized design over the uncertainty set, then it violates the lower bound, in contrast to the robust optimal design, as shown in Table 2.This is because the nominal optimization problem (12) does not take uncertainties into account and chooses the area of the magnet to be too small, so that it does not achieve a sufficient average torque when uncertainties occur.However, the robust problem (18) resolves the uncertainties and accordingly selects a larger area of the magnet, achieving a sufficient average torque that is feasible even in the worst case.In summary, robust optimization can handle uncertainties and still achieve comparable performance.
For the nominal optimization problem, the computation is completed in about 21 min and requires 102 function evaluations.A comparison of the nominal optimization method with an evolutionary algorithm can be found in [9].In [27], it is shown that the nominal optimization method can simultaneously reduce several objectives relevant to electric machines, such as the cost, torque ripple and losses, even with a large number of optimization variables.
The robust optimization requires 79 min for a similar number of iterations.This is primarily due to the iterative determination of worst-case scenarios, which results in an increased number of 405 function evaluations.However, as the optimization progresses, the number of iterations required to compute the worst case decreases.After only a few iterations, the number of iterations is quickly reduced, as observed in Figure 6.This is mainly for the following reason.During the outer loop, the optimization variables of the robust optimization problem (18) converge and thus also the solutions of the worst-case problems (17) converge.Since we start the solution of ( 17) in each iteration with the solution obtained in the previous outer iteration, we have already a good starting point and require only a few iterations to solve the worst-case problems (17).The fact that the robust problem is solved at about four times the cost of the nominal problem underlines the efficiency of the proposed approach.The details are summarized in Tables 2 and 3.   Since the current implementation is done in MATLAB ® , the code is not yet optimized for performance.Calculating the derivatives and solving the problem for different rotation angles is, for example, not yet parallelized.One limitation to note is that the calculations do not account for saturation (nonlinear material behavior), which is essential for accurate design.Therefore, we provide a comparison for the magnetic flux density for the optimized design and the robust optimized design in Figure 7.The magnetic flux density for the initial design is shown in Figure 3.We focus next on the achieved reduction in the torque ripple.By applying free-shape optimization at the air gap by moving the control points, as shown in Figure 4b,c, we are able to reduce the torque ripple significantly.Figure 5 shows the torque profiles during nominal and robust optimization.In Figure 5a, a significant reduction in torque ripple is observed after only 20 iterations, similar to the results in [9].In Figure 5b, we observe at iteration zero that the mean torque from the nominal design is initially insufficient for the robust design due to uncertainties in the optimal design.To address this, the mass of the magnet is increased to satisfy the torque requirements.This early adjustment phase is typical for robust optimization methods and highlights that solutions that are feasible under nominal conditions may not be as feasible under uncertainty.For the first iteration of the robust optimization, assuming the worst case for every parameter even yields torque ripples comparable to the initial design.After several iterations, the optimizer manages to reduce the torque again quite significantly, but not as much as the nominal optimization.Overall, both the nominal and robust optimization methods show a reduction in torque ripple within the first few iterations, which is also the main contributor to the rapid decrease in the objective function value.

Conclusions and Discussion
In this paper, we have proposed an algorithm for the robust optimization of complex PDE-constrained problems.We have applied the method to the combined shape and parameter optimization, with uncertainties incorporated in both.The underlying PDEs are discretized using IGA.This yields a precise representation of the geometry and provides a convenient framework for shape optimization.While the focus has been on the context of electric machine design, the principles and techniques employed are not limited to this application and can be applied to a wide range of robust optimal design problems.The additional costs for robust optimization are small compared to nominal optimization.In future research, we plan to extend the optimization framework to three-dimensional problems, including nonlinearity and hysteresis.Moreover, we plan to apply our methods to robust topology optimization.solution, reducing the computational cost.In order to compute the gradient in (A1), one needs to determine the derivatives of the stiffness matrix dK ij dC kd and the right-hand side db i dC kd with respect to the control points C kd .Theorem A1.Let the mapping function F : Ω → Ω be defined according to (2).Assume that the stiffness matrix K is differentiable with respect to C kd , which denotes coordinate d (d = x or d = y) of the k-th control point.Further, assume that the material coefficient ν is a constant.Then, the derivative of K with respect to C kd can be decomposed into three components as follows: ji,kd + K Here, Dkd and J F are defined as which is derived from (2).
Proof of Theorem A1.The stiffness matrix K ij for an elliptic problem (9) is defined as We transform this equation in the parametric space Ω.Using the mapping F as described in (2), the stiffness matrix is transformed into where the Jacobian matrix J F is given by (A7).To obtain (A8), we use the transformation rule ∇N i = J −⊤ F ∇ Ni for H 1 functions, because the mapping F fulfills the required regularity assumptions; see, for example, [41].To calculate the derivative of (A8) with respect to the control points C kd , we have to differentiate the terms that depend on the control points: Consequently, we apply the product rule to these terms, leading to (A4) and (A5).

Detailing K
(1) ij,kd : The derivative of ∇N i with respect to C kd is given by We can further simplify this to This shows (A5).
The above approach can be used to calculate the derivative of the stiffness matrix with respect to the control points for an elliptic problem, discretized using IGA.To calculate the derivative of the right-hand side db i dC kd with respect to the control points C kd that also appear in (A1), we can use the same approach as in Theorem 2.

Figure 2 .
Figure 2. Geometry model of the PMSM with highlighted control points as described in (2) and patch boundaries.The model was originally based on [19]. .

Figure 4 .
Figure 4.A visual representation of the optimization results: the initial design (a), the optimized design under nominal conditions (b) and the robust optimized design (c).

Figure 5 .
Figure 5. Optimization history of the torque ripple for nominal optimization (a) and robust optimization (b).The nominal optimization and robust optimization show a significant reduction in the torque ripple at the end of the optimization.

Figure 6 .
Figure 6.Number of iterations to determine the worst-case scenario during the robust optimization process.The number of needed iterations decreases with increasing optimization iterations.

Figure 7 .
Figure 7. Resulting magnetic flux density of the the nominal design (a) and the robust design (b).

: 1 F
Since the second term on the right-hand side vanishes, we deduce We employ the Jacobi formula from matrix calculus, which allows us to express the derivative of |J F | with respect to C kd as∂ ∂C kd (|J F |) = |J F | Tr J −∂J F ∂C kd .

Table 2 .
Comparing the iteration numbers for the nominal and robust optimization.

Table 3 .
Comparing the optimization objectives.