This allows us to start the optimization with an ini- ... the optimum. 3) The optimization may start from an ini- tially feasible ..... 141 K. Preis, C. Magele, and 0. Biro ...
1948
IEEE TRANSACTIONS ON MAGNETICS, VOL. 28, NO. 4, JULY 1992
Geometric Parametrization and Constrained Optimization Techniques in the Design of Salient Pole Synchronous Machines Konrad Weeber, Student Member, IEEE, and S . Ratnajeevan H. Hoole, Senior Member, IEEE
Abstract-The shape optimization of magnetic devices is efficiently performed with field calculation and sensitivity analysis based on the finite element method. Several sequential unconstrained optimization techniques are discussed and evaluated with respect to their application in engineering design. The optimization of the geometry of a salient pole generator so as to achieve a desired field configuration in the airgap is used as an illustrative numerical example to demonstrate the geometric parametrization technique, emphasize the importance of constraints in engineering design, and highlight the advantageous features of the augmented Lagrangian multiplier method for nonlinear constrained optimization. For the required geometric parametrizationa recent novel use of structural mapping is extended to incorporate constrained optimization. For the benefit of the electrical engineering readership, the associated equations of structural mapping are presented.
INTRODUCTION
T
HE ANALYSIS of electromagnetic systems is now a well developed and demonstrably accurate art. Using numerical techniques such as the finite element or the boundary element method allows the prediction of the magnetic field patterns and thus the performance of devices and, in particular, electric machines. In practical engineering, however, such an analysis is rarely required exclusively but rather in the context of synthesizing the best engineering solution to a given physical problem. For example, the question in engineering is not ‘‘what torque characteristics will this motor produce?”, but rather “what motor will give me these torque characteristics?” Thus, the task is to optimize the device by determining its dimensions, material properties and field sources to meet specified performance criteria. In recent years there has begun to emerge a methodology for formalizing this optimization process based on the finite element method for magnetic field calculations [ 11-[ 121. The field analysis is part of an iterative cycle that attempts to find the optimum design that fulfills defined performance criteria. Manuscript received November 1I , 1991; revised March 13, 1992. This work is supported by the National Science Foundation’s Division of Electrical, Communications and Systems Engineering and RUI Program under Contract ECS-8809758, the Harvey Mudd College-Southem California Edison Center of Excellence in Electrical System, as well as the Kurt Godel Stipend of the Ministry of Science and Research, Vienna, Austria. The authors are with the Department of Engineering, Harvey Mudd College, Claremont, CA 91711. IEEE Log Number 9200741.
The task of optimally designing electric machines involves determining the dimensions of the magnetic circuit, of windings, and of permanent magnets to obtain a desired performance criterion. In the different optimization techniques presented in the past [13]-[17] the machine analysis is mainly based on algebraic sizing equations and reluctivity methods. In these cases the machine analysis is computationally inexpensive but based on crude approximations for the field distributions. The objective is to minimize the material and operating costs or to maximize the efficiency with the desired performance characteristics imposed as separate constraint functions. The resulting nonlinear constrained optimization problems are solved by different mathematical programming methods available in the literature, such as SUMT techniques [ 131, [ 161 and boundary search methods [ 141, [ 151. Optimization of electric machines based on computationally more expensive but more accurate numerical field calculations has been investigated just recently. In [ 181 the field distribution is determined by the finite difference method and the sensitivity analysis in the optimization procedure is performed by finite difference steps. Field calculations by finite elements and sensitivity analysis by finite difference steps are reported in [19]. The success of shape optimization procedures in general depends strongly on the method used for the parametrization of geometric contours and on the corresponding finite element discretization. The iterative procedures are often hampered by discontinuities in the object and constraint functions as well as by nonsmooth geometric shapes that are obtained as the result of the optimization [ 11, [ 101. However, the application of powerful gradient methods for optimization requires the continuity and differentiability of object and constraint functions. We will present a novel extension of a structural mapping technique for the geometric optimization of magnetic devices. With this technique the displacements of surface grids are used as parameters and are mapped onto the finite element mesh. Different geometric parameter values then result in a finite element mesh of unmodified topology and evenly distorted elements. The implementation of this new technique results in continuous and differentiable object functions and smooth geometric contours of the shapes to be optimized. Thus, the gradient methods employed converge reliably and the geometric shapes obtained are sat-
0018-9464/92$03.00 0 1992 IEEE
1949
WEEBER A N D HOOLE: GEOMETRIC PARAMETRIZATION A N D CONSTRAINED OPTIMIZATION TECHNIQUES
isfactory and, most importantly, manufacturable. In the body of this paper several transformation techniques are discussed for the solution of the nonlinear constrained optimization problem. Finally, the augmented Lagrangian multiplier method is used to solve a numerical example. In these calculations we emphasize the requirement of design constraints to achieve unique optimization results. BY FINITEELEMENTS DESIGNSENSITIVITY Due to the characteristic geometry of electric machines the magnetic fields mainly vary in the cross-sectional plane and are constant in the axial direction, except in the end regions and radial cooling slits. Two-dimensional finite element analysis is now a well developed and demonstrably accurate art for the field calculation in the design cycle of electric machines. Classically, in finite element computation, the magnetic field problem is formulated in terms of the magnetic vector potential A as the state variable which is computed from the matrix equation
{ A ) = {Q(x,Y , J > > , [P(x,Y , CL)]
(1) where the matrix [PI and the column vector { Q } are assembled from permeabilities p , current densities J , and the geometry of the solution region [20]. The resulting knowledge of { A } is used subsequently to determine quantities such as the magnetic flux densities, forces and torques. The solution of (1) represents a direct analysis, where the field calculation is performed for given configurations of materials, geometries and excitations. For the purpose of device optimization we select key parameters, such as geometric dimensions, material properties and current distributions, that define the design. These n design parameters are represented in vector notation as { p } = { pi * * * p i * * * P , , } ~ .Both the coefficient matrix [PI and the field excitation vector { Q } are thus functions of the design variables p i , which implies that the solution for the potential { A } also becomes design dependent. The desired performance of the machine is formulated in terms of an object or quality function F and several constraint functions g, which impose conditions on the minimum of the object function. In the most general case F and g are design dependent both explicitly and implicitly through the design dependence of the potential distribution: F = F({ P I , A ( { P}>)g= 8({ P I ? A ( { P I ) ) . (2) As a major segment of the optimization algorithm, quantitative information on how the performance of the device is affected by changes in the variables p i , is obtained through design sensitivity analysis [3], [6], [21]. For that purpose the total design dependence of object and constraint functions is determined by the total derivatives of F and g with respect to each parameter p i as
i = l * - - n (3a)
i= 1
. . . n.
(3b)
The partial derivatives of F and g represent the explicit dependence on parameters and field distributions and follow directly from the algebraic expressions for F and g. To obtain the design dependence of the field distributions, represented by the partial derivative of the state variable A , both sides of the governing (1) are differentiated with respect to each optimization parameter p i to obtain the following identities
i = l
. . . 12.
(4)
Since the matrix [PI and the vector { Q } are expressed in terms of the physical geometry of the model, their derivatives with respect to any physical parameter are expressed by direct differentiation of the element matrices without need for an additional field computation, as presented in [3], [5], [7], [8], [21], [22] for linear and in [5] for nonlinear magnetic materials. The known potential solution of { A } from (1) is then used in the assembly of the right hand side of (2). It is to be noted that both (1) for the state variable { A } and (2) for its partial derivative are characterized by the same coefficient matrix [3], [9]. Therefore, we preferably use Choleski-decomposition in the solution of both equations: the computational effort for the calculation of a single additional gradient vector { a A / a p i } is then reduced to forward and backward substitutions, using the already decomposed cOel€icient matrix from the solution for { A } [9]. With the. knowledge of { aA /api } the sensitivity of the two-dimensional planar components and the magnitude IB( of the magnetic flux density follow directly as
a’B’ - a ({Ae}T[P,]{A,})”’ aPi api
(5c) In these equations subscripts e denote element quantities and { b } and { c } are the differentiation vectors of first order triangular elements [20]. The derivatives of the magnetic potential and the flux density give a measure for how much the field solution changes due to variations in the
1950
IEEE TRANSACTIONS ON MAGNETICS, VOL. 28, NO. 4, JULY 1992
ith parameterp,. It is this sensitivity information which is available for both linear and nonlinear materials, that provides the essential guideline for the direction in which the parameter values are modified during the optimization procedure. It may be noted that even a design sensitivity analysis that is not incorporated in an optimization algorithm has its own right. It provides designers with information on the trends of the device performance as they change single parameter values and it may thus be used directly in any interactive computer aided design procedure.
where the gradient vector of the object function is given by 7
VFr =
dF({ P},, A ( { PIr)) dF({ P>r, A ( { PI,)) dP2
(9)
An improved quadratic convergence rate is obtained in the Newton method by utilizing second order gradient information in the form of the Hessian matrix. However, the OF GENERAL CONSTRAINED computational effort for obtaining the required second orNUMERICAL OPTIMIZATION der gradients with respect to all parameters may become PROBLEMS prohibitively large, so that the Hessian is frequently The determination of the optimal shape, material propmerely approximated by collecting first order informaerties and current distributions as the device descriptive tion. Those methods are referred to as quasi-Newton or parameters can be reduced to the minimization of a qualvariable metric methods, and the most frequently used ity or object function subject to constraints for fields, flux ones are the Broydon-Fletcher-Goldfarb-Shanno (BFGS) linkages, and parameter values. This is stated in the mathand the Davidon-Fletcher-Powell (DFP) methods [24]. ematical formulation of a general constrained optimiza- Once the search direction is known, the rth line search tion problem: determines the new parameter vector by minimization of minimize the object function the function F along the present search direction. This F({ P I ? A({ PI)) ( 6 4 distance of travel in the direction of {s}, is given by the scalar factor CY subject to inequality constraints (10) {PI, = { P > r - l + CY{slr. s,, A{(P})) = 0 , (6c) tion methods [23]. In these methods the violated conand side constraints straint functions are weighted by a penalty factor w q and added as a penalty function R to the original object func= 1 * * * n. i PL,, 5 PI 5 PU,,, (64 tion F to give the unconstrained modified or pseudo-obThe object and the constraint functions may be nonlinear ject function +: and may depend implicitly on the design variables through (11) +({P>, W J = F({PI) + R({PI7 wq). the state variable A of the magnetic field distribution. The side constraints give the lower and upper bounds of the The constrained optimization of (6)(a-d) is thus reduced parameters and thus define the region of search in the to a series of unconstrained problems for the minimization n-dimensional parameter space. For the remainder of this of which leads to the term sequential unconstrained paper the side constraints are included in the m inequality minimization techniques (SUMT). The pseudo-object constraints by defining the following normalized confunction tends to be numerically ill-conditioned for large straint function penalty factors. Therefore, in the initial stages of the optimization, a small penalty factor w q is used and increased in subsequent iterations q , imposing a more rigorous pen(7) alty for constraint violations. As a fundamental part in constrained optimization let In the external penalty function method (EPFM), R is us first consider unconstrained optimization methods [23] assembled by penalizing the object function only when which are separated into two basic tasks. The first is the constraints are violated ( g j ( { p}) > 0 or hk({ p}) # 0), determination of a search direction {s} in the parameter m space, which will improve the object function. In the conjugate gradient method the new search direction at iteration r consists of a term equal to the steepest descent direction updated by and additional term based on the + wq k = 1 [hk({P))l2. (12) information of the previous line search: The constraints are squared to ensure a continuous slope of the penalty and pseudo-object functions at the constraint boundary. For low values of the penalty factor the pseudo-object function is well behaved, but the constrained optimum is reached only in the limit of w q -+ 00
-
+,
1951
WEEBER AND HOOLE: GEOMETRIC PARAMETRIZATION AND CONSTRAINED OPTIMIZATION TECHNIQUES
and is approached from the infeasible side. From the engineering perspective this drawback is significant, because the design will not be feasible if the optimization is terminated prematurely. However, the advantages of the EPFM are that 1) the penalty function is continuous at the constraint boundary, 2 ) the original object function is not modified by penalty terms inside the feasible region, and 3) The penalty function is defined outside the feasible region. This allows us to start the optimization with an initially infeasible design. The interior penalty function method (IPFM) for inequality constraints leads to a sequence of improving designs, where the constrained optimum is approached from inside the feasible region. The penalty function is generally of the form
f
p2 gl=O
g2=0
Pl
Fig. 1. Geometric interpretation of the Lagrange multipliers at the constrained optimum.
with
m
leading to a singularity at the boundary of the feasible region when g j ( { p } ) 0 . Again, caution has to be exercised in the choice of large penalty factors, which result in steeper slopes at the constraint boundaries and in smaller penalty terms. Additionally, we have to ensure that the search procedure starts from within and never leaves the feasible range. This IPFM is frequently applied for imposing parameter side constraints, where one can easily pick feasible starting values. The advantages of both EPFM and IPFM are combined in the augmented Lagrange multiplier method (ALMM), which is the most commonly used method for transforming the constrained problem into an unconstrained one for sequential minimization [ 181, [ 2 3 ] . By including the conditions for optimality, the dependency of the algorithm on the choice of the penalty factors and their updating is reduced. The necessary Kuhn-Tucker conditions for a constrained optimum at { p* } are equivalent to the minimization of the Lagrangian function L( { p } , { X}) -+
m
VU{P*I, {XI) = VQ{P*I) + J,c XjVgj({p*)) =1 1
+ kc XkVhk({p*}) = 0. = 1
(14)
The Lagrange multipliers X are larger than zero for all active inequality constraints (g,({ p * } ) > 0), and unequal to zero for all active equality constraints (hk ({ p* }) # 0). A geometric interpretation of (14) is that the gradient vectors of the object function and of the active constraints have to point in opposite directions at a constrained optimum, as indicated in Fig. 1 for two active constraints. It has been shown in [26] that the optimization of the constrained problem is identical to the minimization of the augmented Lagrangian m
1
This augmented Lagrangian is basically an exterior penalty function formulation of (12) augmented by additional Lagrange-multiplier terms. The purpose of the sequential minimizations is now to determine the values of the unknown Lagrange multipliers as well as the penalty factor. Although this increases the number of unknowns in the SUMT iterations, some significant computational advantages are gained. For the case of all Xi, X k = 0, (15) reduces to the classical EPFM, which reaches the constrained optimum from the infeasible region, and only for wq 03. However, if the exact values of the Lagrange multipliers that minimize L A ( {p } , {A}, w q ) are known at the beginning, only a single unconstrained minimization is required, which yields an exact constraint satisfaction for a positive, but finite value of wq. As these exact Lagrange multipliers are not known a priori, they are updated at each SUMT iteration q according to -+
AT+'
= A4
+ 2w9 q rJa
( 16a)
A$+' = + 2wqhk({ p } ) . ( 16b) In addition to updating the Lagrange multipliers, the penalty factor wq is increased only modestly at each SUMT iteration, as it does not have to approach large values at the design optimum. The ALMM offers the following four advantages: 1) It is relatively insensitive to the value of wq. 2 ) Precise constraint satisfaction can be achieved at the optimum. 3 ) The optimization may start from an initially feasible or infeasible design, and finally, 4) Nonzero values of the Lagrange multipliers automatically identify the active constraints at the optimum. In addition to the transformation methods discussed in this section, primal methods are used in constrained nonlinear optimization [ 1 2 ] , [19], [ 2 7 ] . In the course of these methods no penalty functions are imposed on constraint violations, but the constraints are dealt with directly, most often by linear or quadratic approximation techniques. GEOMETRIC PARAMETRIZATION: A STRUCTURAL MAPPINGTECHNIQUE The equation outlined in the previous sections represent the framework within which the optimization procedures attempt to approach the optimal design. These methods
1952
rely heavily on the gradients of object and constraint functions, which are obtained by the design sensitivity analysis based on the finite element formulation. The continuity and differentiability of F and g are basic requirements for the application of these gradient methods and the accuracy of the gradient calculations is of crucial importance to the stability and robustness of the optimization algorithm. In recent papers [lo], [28], two key difficulties in geometric parametrization have been recognized. First, as the geometry changes, a new finite element mesh is necessitated for each new value of a geometric parameter. In using free meshing based on Delauney triangulation [20], [29], [30], topologically different meshes arise, each one representing a possible domain discretization with an inherent error. This change in mesh topology and its associated discretization error occurs in discrete steps, leading to discrete changes in the potential and field solutions which automatically implies discontinuities in the object function. Through the finite element discretization, the actually continuous object function becomes only piecewise continuous. This jeopardizes the smooth convergence and restricts the reliable application of gradient methods. Although tunneling [3 11 and stochastic methods [4] are effective in overcoming this problem, it is far simpler to remove these artificial discontinuities through proper meshing as suggested in [ 101, [28]. It is interesting to note that the number of object function discontinuities which are large enough to terminate numerical search methods increases with the fineness of the mesh. Therefore, for finer discretizations and more accurate finite elements meshes it is more likely that line searches terminate prematurely. Furthermore a sensitivity analysis based on finite difference steps becomes inaccurate and often meaningless, because the gradient values can vary significantly due to the changes in mesh topology. On the contrary, when the gradient dF/dpi is calculated by the design sensitivity analysis of (3)-(5) by direct differentiation of the coefficient matrices, it relies on finite element information of a single given mesh topology. It is thus accurate and not affected by topological changes. As a second requirement, geometric parametrization has to impose constraints on the smoothness of the geometric contour. When selecting nodal coordinates at the shape surface as parameters, the boundary is modeled by piecewise linear segments which do not satisfy continuity conditions. In this case the geometric shapes obtained through optimization are characterized by jagged contours that are not manufacturable, although the object function is truly minimized [l], [4], [lo], [28], [32]. Using coefficients of polynomials as parameters in the boundary description results in smooth shape outlines that however tend to undulating contours for higher order polynomials. Therefore, preferably, parametric curves such as cubic spline functions with two continuous derivatives [32] and Bezier curves [33], [34] are applied in shape representation.
IEEE TRANSACTIONS ON MAGNETICS, VOL. 28, NO. 4, JULY 1992
Therefore the two requirements for a robust and reliable geometric parameterization can be defined as: 1) A preserved mesh topology during optimization, and simultaneously 2) The use of C’-continuous curves in the shape description. As we want to optimize the shape of a device, it is reasonable to modify an initial guess by deflecting it according to structural laws of elasticity [35]. Fig. 2 shows how the geometry of a rectangular plate is deformed by three different applied forces. The structural consistency of the plate guarantees smooth shape contours, which depend on the value of the applied force and the structural material characteristics of stiffness and Poisson’s ratio. The main advantage of using this structural deformation as a parametrization method is that the mesh topology does not change: As geometric parameter values are modified, the new geometric shapes are represented by distorted finite element meshes. The distortion of the finite element mesh is proportional to the changes of the parameter values and thus increases continuously. This implies continuous changes of discretization error, field solution and object function. In the two dimensional analysis of plane stress, the constant strain element is one of the earliest and best known elements [36]-[38]. It is used for the discretization of the continuum of the structural design model, because of its nodal compatibility with the first order triangular element now commonly used for the electromagnetic field analysis. Thus, the finite element models for structural and electromagnetic computations consist of identical triangular elements in the common region. In the structural finite element equation, [KI{u) =
{f),
(17)
the stiffness matrix [K], depending on the material properties of elasticity modulus E and Poisson’s ratio U , relates the displacements in the x - and y-directions as nodal degrees of freedom { U } to the loads { f } , which are the nodal forces acting in the two coordinate directions. The deflections of a structure tend to be centered around the point where the loads are applied, and we therefore have to use methods to distribute the local actions of applied point forces and specified nodal displacements more evenly throughout the solution region. For this purpose the two-dimensional structural finite element model is augmented by additional one-dimensional elements. They are placed mainly at the surface contour of the geometry to be optimized. Selecting a higher stiffness for these elements than for the planar elements renders surface contours whose shapes depend mostly on the deformation of the one-dimensional elements. The softer planar elements thus react like a “soft sponge” and are then merely used to map the applied surface deformation evenly onto the solution region. The one-dimensional element used for this purpose is the beam element of Fig. 3 with three degrees of freedom
1953
WEEBER A N D HOOLE: GEOMETRIC PARAMETRIZATION AND CONSTRAINED OPTIMIZATION TECHNIQUES
Fig. 2. Example of a structurally distorted finite element mesh with three different parameter values.
'A
L
40 Fig. 3 . One-dimensional beam element for axial and flexural behavior with 3 degrees of freedom per node.
per node: the axial displacement U , the transverse displacement U , as well as the rotation 8 in the z-sense. The rotation itself is identical to the slope of the bending curve 0 = dv/d.x [35],[36],[38].Thus, the element vector of the nodal degrees of freedom for both axial and bending behavior is written as (U} = (uI
v1
el
u2 u2
e2}T, with
The corresponding column vector for the nodal actions contains point forces in the axial and transverse directions as well as the bending moments in the z-sense,
{f)= {fxl
4
(19) The exact element stiffness matrix within beam theory [38] represents both the axial and bending behavior fyl
fx2
W 2 Y .
fy2
-
bhL2
0 E [KI, = 3
the bending behavior. The axial behavior is characterized by the product of elasticity modulus E and beam cross section bh, whereas the bending behavior is characterized by the flexural rigidity EI. For the assembly of the equations of the structural design model, this decoupling creates the option of selecting material properties and beam cross sections in such a combination as to modify and influence the axial and bending behavior independently. We can thus choose the beam characteristics such that the flexural rigidity is small enough to give locally yet smoothly deflected beams. Furthermore the axial stiffness can be set large enough to distribute the axial deflections evenly throughout the finite element model. The nodal displacements and forces in a local beam coordinate system are transformed by rotation-of-axis transformations into the global coordinate system, in which the nodal displacements of triangular elements are defined (Fig. 4 ) . With these transformations the global equilibrium equations in the nodal degrees of freedom resulting from element assembly pertain to the global coordinate system [36].The nodal connection of beam elements requires some further discussion. If the rotation 8 of beam elements is set equal at common nodes, the &distribution is forced to be continuous along beam elements which leads to a C1-continuous beam deflection curve and geometric contour. This inherent smoothness of the shape is, however, not desired when sharp edges and comers have to be modeled. In this case the &components of neighboring beam elements are not set equal at the common node, resulting in a discontinuous beam deflection curve. These hinged beams are connected by multipoint constraints applied to common nodal wand u-components (Fig. 5). The following two examples briefly demonstrate the characteristics of the combined beam and triangle elements for the purpose of obtaining uniformly deflected finite element meshes. A rectangular plate is augmented by beam elements at the upper boundary, where geomet0
-bhL2
0
0
121
61L
0
-121
61L
0
-6IL
21L2
61L
41L2
-bhL2
0
0
0
-121
0
- 0
-
0
61L
In this expression L denotes the length, b the width, h the height and I the moment of inertia of the cross section of a beam element. This stiffness matrix relates the axial and transverse displacements U , U , and the rotations 8, to applied axial and transverse forces and moments in the z-sense, respectively. As the influence of shearing deflections are neglected, the axial behavior is decoupled from
-61L 21L2
bhL2 0
0
0
0
121 -61L -61L
4ILz
ric deflections are imposed by specified displacements at a single node. The plate is fixed at the bottom, and the horizontal displacements are constrained at the right boundary. Fig. 6 gives the deflection curve for three different values of the elasticity modulus E of the beams, which determines their flexural rigidity and bending be-
IEEE TRANSACTIONS ON MAGNETICS, VOL. 28, NO. 4, JULY 1992
1954
X
beans
b
CONTINUOUS BEAMS
=
IO 0. b
= 1
0.
h =
2 0
(a)
CONTlNUoUS BEAMS
I
82
E
€1 = 6 666 ER = 20 0
Fig. 4. Displacements of triangle in global coordinate system (unprimed) and of beam element in local coordinate system (primed).
I
8 3 04
05
I U3 = U4
HINGED BEAMS v3 = v4 e3
beams
#
64
Fig. 5 . Joint of hinged beams.
E-beans
= 10 0
-
E-beans = 30 0
E-beanr
=
mo 0
(a) 6) (C) Fig. 6 . Deformed mesh of rectangular plate augmented by beam elements at upper boundary.
havior. The curvature of the deflection depends on the flexural rigidity in the same way as coefficients of a spline function govern its distribution. To obtain the deflected shapes of Fig. 7, further beam elements are added at the left boundary, and the joint of vertical and horizontal beams is modeled by an additional hinge node. Fig. 7(a) shows a model with axially soft beams: the axial deflection of the vertical beams is concentrated around the node where the vertical displacement is applied. Due to the axially stiffer beams in Fig. 7(b), the specified displacement is distributed uniformly over the whole length of the vertical plate boundary, resulting in a more evenly deformed mesh throughout the model. It can be seen that the use of hinged beams of high axial stiffness and relatively soft flexural behavior yields a uniformly distorted finite element mesh in the whole region of the structural design model. It is important to emphasize that the structural design model is not used to obtain a numerical solution of an existing structural problem. We rather want to establish a reliable mapping technique that yields a uniformly deflected finite element mesh mainly for complex geometries. For this purpose we are free to assemble the struc-
E i 100 0. b E l = 6 666 ER = 2w 0
i
3 162.
h = 0 632
(b) Fig. 7 . Even distribution of axial beam displacements. (a) Axially soft beams: deflections concentrate in vicinity of applied displacement. (b) Axially stiff beams: deflections uniformly distributed over full beam length.
tural design model by augmenting sections of the electromagnetic finite element model appropriately with beam elements and hinge nodes. In the following we will outline how this structural design model fits into the electromagnetic shape optimization procedure. The structural design model may be deflected either by applying forces { f } as loads, or by specifying nodal displacement components in the vector of state variables {U}. The analogy in the analysis of electromagnetic fields is to “drive” magnetic fields either by specified current and charge distributions, or by defined nonzero potentials. Point forces applied to the surface of the structure have been selected as parameters in a previous publication [39]. However, in using beam elements with different axial and flexural behavior, axial and transversal forces result in displacements of different size. The design model thus becomes more sensitive to transverse forces. This unwanted scaling effect is avoided by selecting specified surface displacements as the design parameters. Then the same increment in a parameter value (that may describe axial or transverse deflections) results in the same amount of change in geometry. The displacement vector may then be partitioned into the unknown U, and the specified U, components, to yield the finite element equation of the structural model without force components as loads in partitioned form
WEEBER AND HOOLE: GEOMETRIC PARAMETRIZATION A N D CONSTRAINED OPTIMIZATION TECHNIQUES
The unknown displacements are then determined by solution of the first row
This leads to the term “structural mapping”: the specified surface displacements U, as parameters are mapped onto the remaining, unknown displacements U, by the mapping matrices K,, and K,,. This mapping relationship is governed by the physical properties of the structural design model which may be chosen in such a combination as to obtain the desired geometric shapes. The solution of (22) is performed for two different purposes in the process of shape optimization. First, the initial geometry of the starting guess is subjected to the specified surface deflections { u , } ~rendering , a set of basis ( i1 * n) deflection shapes { ~ ~ } ~ =
1955
nated. The smoothness of the geometric contours is guaranteed by the use of sufficiently stiff beam elements in the design model at the surface contours. The second application of the mapping (22) arises in the calculation of the magnetic potential gradient { a A / d p i } in sensitivity analysis, which is generally required at the beginning of the rth line search during unconstrained minimization. For the assemblage of the right hand side of (4), the derivatives of the element matrices [PI, and vectors { Q } , with respect to the geometric parameters are required. These derivatives are assembled from the partial derivatives of the element matrices with respect to the coordinates of the triangle nodes and the derivatives of these nodal coordinates with respect to the parameters as
-
The parameters p i during the optimization are then the factors of the surface deflections { u , } ~applied to the structure. Each vector may consist of a single displacement component, whose action is distributed over the complete surface by beam elements. The resulting smooth contour is determined by the flexural characteristics of the beam elements. On the other hand, each vector { u , } may ~ as well consist of a group of specified surface displacements. In this case, it may define the deformation of whole geometric boundaries either by a given constant, linear, or polynomial distribution. The surface displacements specified in each { u , } ~may even resemble the geometric contours of available shape profiles, so that the optimization actually determines the coefficients of an optimal superposition of those different profiles. Therefore quite complex changes of the geometry are essentially represented by a single parameter, which may be effectively applied to reduce the dimensionality of the optimization problem and the involved computational effort. This way the structural mapping technique automatically incorporates the reduced basis concept [40]. The set of basis vectors { ub } i and the nodal coordinates of the starting geometry { x } , are used throughout the optimization to obtain the parameter dependent geometry b < P ) >as
This procedure thus replaces the need for free meshing: For new parameter values the finite element mesh is obtained as a linear superposition of the basis deflection shapes. The resulting mesh never changes its topology, and therefore changes in the discretization error and resulting discontinuities in the object function are elimi-
The detailed expressions for the sensitivities of the element matrices for geometric, material and excitation parameters are given in references [7], [8]. It is the calculation of the gradient of the nodal locations with respect to the geometric parameters dxn/api‘and ayn/api that again involves the solution of the structural mapping (22). The same set of surface deflections as used in (23) is now applied to the present shape { x } of the structure. This is equivalent to applying a unit step of the parameter p i (recall: p i is the coefficient of {U, } i and { ub } i ) to obtain , in turn are the corresponding displacements { u , } ~which the desired derivatives of the nodal coordinate with respect to the parameter p i ,
for a unit change in p i . Thus, the gradients of the nodal coordinates are obtained by a structural solution for the geometry of the present parameter vector. The structural mapping method replaces a geometric parametrization involving free meshing as follows: At the beginning of the optimization procedure n different displacement vectors are defined which characterize the parameterized surface deformations. An n-dimensional basis of deflection vectors is calculated by the structural solution of the initial geometry loaded with the specified displacements (23). Whenever new values of geometric parameters require a new finite element discretization, this is obtained by using the present parameter values as coefficients in the linear superposition of the basis vectors
1956
IEEE TRANSACTIONS ON MAGNETICS, VOL. 28, NO. 4. JULY 1992
T
I
'
(24). The gradient of the object function requires the derivatives of the nodal coordinates with respect to the parameters, which in turn are obtained from a structural solution of unit parameter displacements applied to the present shape (26). Therefore, at the cost of an additional structural solution per line search with a model containing just the regions of variable geometry, the prerequisites of continuity and differentiability of the object function are achieved. Furthermore, smooth shape outlines are obtained even with only a few optimization parameters, as will be shown in the following applications.
ROTOR: pr = 2OO0.0
A NUMERICAL EXAMPLE: THE SHAPEOPTIMIZATION OF A SALIENT POLE
I 0.1
The following very familiar numerical example of determining the pole face contour of a salient pole synchronous machine is considered to demonstrate the structural mapping technique and its contribution to constrained optimization. It is not viewed as a definite design optimization of a given synchronous machine. For this reason we have simplified the stated problem by neglecting stator slotting, three-dimensional field components and nonlinear saturation effects. The current density in the excitation coil and the geometric parameters that define the shape of the pole piece have to be predicted in order to achieve a sinusoidal distribution of the airgap flux density with a peak value of 1.O T (the time-dependency of the flux density wave and thus the induced voltage wave may be easily constructed from the spatial flux density distribution with the knowledge of the rotational speed of the machine). Fig. 8 gives the parameterized geometry of half a pole pitch that is used as solution region in the no-load field analysis. The stator is idealized as a solid steel region without slots, and both stator and rotor are made of linear steel with a relative permeability of 2000. To determine the deflection shape functions by structural analysis, a subregion is used which contains only the areas of variable geometry and extends from the upper edge of the coil to the airgap interface of the stator, as indicated in Fig. 8. This structural model represents thus only a fraction of the electromagnetic model and hence the computational effort to determine the shape functions and the partial derivatives of (25) is kept at a minimum. The parameters pl-p3 define uniform deflections of the left, lower and upper pole piece contour, whereas the parameters p4-p6 define deflections of a single grid point each, resulting in contours that are governed by the flexural behavior of the beam elements at the pole piece surfaces. The material properties of the structural triangular elements ( E = 0.5, ZI = 0.3) and the beam elements ( E = 5.0, width = 10.0, height = 0.05) guarantee that the characteristics of the mesh distortion depend mainly on the flexurally soft beams. The objective of achieving the specified sine-distribution of the airgap flux density just below the stator (shaded area in Fig. 8) is formulated with
,
0.09
0.01
0.1
Fig. 8. Parametrized geometry of salient pole. Dimensions in [m].
the following least square object function which describes the deviation of the numerically computed y-component of the flux density Bj,,sfrom the desired values Bopt,$at defined sampling points s in the airgap F(Pl9
*
'
9
Pn) = F ( { P } ) =
c (By,s
-
Bopt,d2. (27)
Each component of the gradient vector V F is obtained as a summation of terms of the following form
Here, the first term of the partial derivatives vanishes due to the nature of the object function, which does not explicitly depend on p i . The second term involves the sensitivity of the y-component of the flux density with respect to the design parameters which is given by (5b). To prevent mutual penetration of different geometric components the parameter side constraints are implemented in the form of (1 1). Additional constraints for the geometric parameters of the pole face contour impose a minimum airgap of 6, = 0.02 m. A crucial constraint on the operating characteristics imposes bounds on the leakage flux 9, of the field winding, which has to be smaller than a specified fraction of the main flux CP, entering the stator. These flux quantities are calculated with the magnetic vector potential at the points E , F, G as indicated in Fig. 8. The resulting constraint function is given by a linear combination of magnetic potential terms and exhibits an implicit and nonlinear dependence on the design parameters { p } ay
d{PH
