|Home | Articles | Forum | Glossary | Books|
Optimization methods try to find the maximum or minimum of a function, where there may exist restrictions or constraints on the independent variables. Finding the maximum or minimum of a function which is also the global maximum or minimum of the function has considerable difficulty and is the complex part of any optimization method. In engineering, it is generally considered practical to search only for local solutions.
In the optimization of electrical machines the objective function and constraints can be computed using either the classical (circuital) approach or a numerical field computation approach, such as the FEM. The FEM is more accurate than the classical approach but requires substantially more sophisticated software and computational time. In numerical field computation problems the standard prerequisites of local optimization (convexity, differentiability, accuracy of the objective function) are usually not guaranteed. Deterministic optimization tools for solving local optimization, such as steepest-descent, conjugate gradient and quasi-Newton methods, are not ideally suited to numeric electromagnetic problems. This is due to the difficulties with numerical calculation of derivatives because of discretization errors and numerical inaccuracy. Recently, noniterative optimization schemes have been proposed using artificial neural networks. Stochastic methods of optimization have, however, become more popular over the last few years for optimization of electrical machines due to their high probability of finding global minima and their simplicity. Stochastic methods such as simulated annealing, genetic algorithm (GA) and evolution strategies have been successfully used for different aspects of electrical machine design.
The population-based incremental learning (PBIL) method is a stochastic nonlinear programming method which shows many advantages over the present stochastic methods. Stochastic optimization has the disadvantage that it is not very efficient. This problem is compounded by the large computation time of the numeric field computations. A solution to this problem is the use of response surface methodology, which has been successfully used for optimization of e.g., PM d.c. commutator motors.
1. Mathematical formulation of optimization problem
The optimization of an electrical machine can be formulated as a general con strained optimization problem with more than one objective, i.e., minimization of costs, minimization of the amount of PM material, maximization of the efficiency and output power, etc. Finding the extremum (Extr) of vector optimization problems is defined as:
Extr F(x)= Extr [f1(x),f2(x),...,fk(x)] (eqn.1)
... where F : _n ?_k gi,hj : _n ?_ x ?_n x =(x1,x2,...,xn) subject to specified equality and inequality constraints:
gi(x) = 0,i =1, 2,...,m (eqn.2) hj (x)=0,j =1, 2,...,p (eqn.3) and specified limits for the independent variables:
xmin = x = xmax (eqn.4)
In eqns (1) to (3) F(x) is the vector objective function with the objectives fi(x) to be minimized, x is the vector of design variables used in the optimization, gi are the nonlinear inequality constraints, hj are the equality constraints, and xmin and xmax are vectors of lower and upper bounds for the design variables.
In vector optimization problems there is a conflict between the individual objective functions fi(x) since there exists no solution vector x for which all objectives gain their individual minimum. Vector optimization problems can be transformed from multi-objective optimization into single objective optimization using the method of objective weighting. Although objective weighting always leads to a non-inferiority (Pareto-optimal) feasible solution, the estimation of the weighting factors and the optimization starting point are subjective choices and their influence can rarely be estimated in advance.
A more practical method of optimization is to minimize only one objective function while restricting the others with appropriate constraints. Most constraints will be upper or lower bound inequality constraints, which means constrained optimization procedures are required. The optimization is thus done for a feasible region, in which all the constraints are satisfied for the design variables.
2. Nonlinear programming methods
In nonlinear programming both the objective and constraint functions may be nonlinear. There is no general agreement to the best optimization method or approach. The extremely vast subject of nonlinear programming has been divided into direct search methods, stochastic methods and gradient methods. A short review of these three classes is described below.
Most numerical field problems have constraints of one form or another. A summary of the main constrained optimization techniques is also given below.
2.1 Direct search methods
Direct search methods are minimization techniques that do not require the explicit evaluation of any partial derivatives of the function, but instead rely solely on values of the objective function, plus information gained from earlier iterations. Direct search methods can loosely be divided into three classes: tabulation, sequential and linear methods.
Tabulation methods assume that the minimum lies within a known region.
The methods of locating the optimum are: (i) evaluation of the function at grid points covering the region given by the inequalities, (ii) random searching assuming that the minimum would be found within a sufficiently large number of evaluations, or (iii) a generalized Fibonacci search finding the solution of the multivariate minimization problem by using a sequence of nested univariate searches.
Sequential methods investigate the objective function by evaluating the function at the vertices of some geometric configuration in the space of the independent variables. This method originated from the evolutionary operation (EVOP). EVOP is based on factorial designs. The objective function is evaluated at the vertices of a hypercube in the space of the independent variables.
The vertex with the minimum function value becomes the center point of the next iteration and a new design is constructed about this point. This is a mutation type search mechanism to direct the search towards the optimal. Fractional factorial experimentation assumes systematic and symmetric vertices to reduce the number of objective function evaluations.
The simplex method evaluates the objective function of n independent variables at n + 1 mutually equidistant points, forming the vertices of a regular simplex. The vertex with the highest value is reflected in the centroid of the remaining n vertices, forming a new simplex. The size of the simplex is reduced if one vertex remains unchanged for more than M consecutive iterations, thus narrowing the search to the minimum.
Linear methods use a set of direction vectors to direct the exploration.
There is a large number of linear methods available such as:
(i) The alternating variable search method which considers each independent variable in turn and alters it until a minimum of the function is located, while the remaining (n - 1) variables remain fixed.
(ii)The method of Hooke and Jeeves which uses exploratory moves and pattern moves to direct the search towards the minimum by attempting to align the search direction with the principal axis of the objective function.
(iii)Rosenbrock's method which uses n mutually ortho-normal direction vectors.
Perturbations along each search direction are done in turn and if the result is no greater than the current best value, this trial point replaces the current point. This is repeated until the minimum is obtained.
(iv)Davies, Swann and Campey method uses n mutually ortho-normal direction vectors and a linear univariate search algorithm on each direction in turn.
After each stage is complete, the direction vectors are redefined.
(v) Quadratic convergent methods minimize functions which are quadratic in the independent variables making use of conjugate directions.
(vi)Powell's method is based on mutually conjugate directions and ensures that a direction is replaced only if by doing so a new set of direction vectors, at least as efficient as the current set, is likely to be obtained.
2.2 Stochastic methods
Simulated annealing (SA). This method generates a sequence of states based on an analogy from thermodynamics where a system is slowly cooled in order to achieve its lowest energy state. This is done using nature's own minimization algorithm based on the Boltzmann probability distribution. Thus, the design configuration is changed from one to two with objective functions f1 and f2 with a probability of P = exp[-(f2 - f1)/kT]. If f2 >f1, the state with probability P is accepted. If f2 <f1, the probability is greater than one and the new state is accepted. At a given temperature, the configurations are arbitrarily changed using a random number generator and the designs are also changed with a dictated probability greater than one. The temperature is thus lowered for the next round of searches. This makes uphill excursions less likely and limits the search space. SA ultimately converges to the global optimum.
Multiple-restart stochastic hill-climbing (MRSH). This method initially generates a random list of solution vectors of the independent variables, using binary vectors. The solution vector corresponding to the minimum result of the objective function is used in an iterative loop. A bit of the solution vector is toggled and evaluated. The minimum for a sufficiently large number of iterations is assumed to be the minimum of the objective function.
Genetic algorithm (GA). This search method is rooted in the mechanisms of evolution and natural genetics. A GA combines the principles of survival of the fittest with a randomized information exchange. GAs generate a sequence of populations by using a selection mechanism, and use crossover as the search mechanism to guide the search towards the optimal solution.
Population-based incremental learning (PBIL). This is a combination of evolutionary optimization and hill-climbing. PBIL is an abstraction of the GA that maintains the statistics of the GA, but abstracts away the crossover operation and redefines the role of the population.
2.3 Gradient methods
Gradient methods select the direction si, of the n dimensional direction vector, using values of the partial derivatives of the objective function F with respect to the independent variables, as well as values of F itself, together with in formation gained from earlier iterations. The solution is thus improved, that is:
F(xi+1) = F(xi) xi+1 = xi + hisi (eqn. 5)
…where hi is the step increment and si is the search direction. Types of gradient optimization methods are:
Methods of steepest descent use the normalized gradient vector at the current point to obtain a new point using a specified step length.
Newton's methods use a second order truncated Taylor series expansion of the objective function F(x). The method requires zero, first and second derivatives of the function at any point.
Quasi-Newton methods use an approximation of the second derivative of the function which is updated after each iteration.
2.4 Constrained optimization techniques
Constrained optimization problems are generally transformed to unconstrained ones and are then optimized using one of the nonlinear programming methods described above. Some of the techniques used on constrained problems are:
Feasible direction method attempts to maintain feasibility by searching from one feasible point to another along feasible arcs. This method assumes a feasible point can be found when the procedure starts.
Penalty function transforms the optimization problem to include the constraints which enable F to be maintained while controlling constraint violations by penalizing them. Exact penalty function is similar to the classical penalty function except that the absolute value of the constraints are used.
Sequential unconstrained minimization technique is also similar to the classical penalty function except the penalty coefficient is increased after each step of the algorithm.
Augmented Lagrangian function or multiplier penalty function uses a Lagrangian function to which the penalty term is added. The problem is trans formed to an augmented Lagrangian function which is minimized.
3. Population-based incremental learning
In terms of Darwinian models of natural selection and evolution, life is a struggle in which only the fittest survive to reproduce. GAs based on natural selection and genetic recombination were first proposed by Holland. GAs generate a sequence of populations by using a selection mechanism, crossover and mutation as search mechanisms. In nature, competition for resources such as food means that the fittest individuals of a species dominate over the weaker ones. This natural phenomenon is called "the survival of the fittest." The fittest individuals thus get a chance to reproduce ensuring, implicitly, the survival of the fittest genes. The reproduction process combines the genetic material (chromosome) from the parents into a new gene. This exchange of part of the genetic material among chromosomes is called crossover.
GAs encode the solution as a population of binary strings. Each solution is associated with a fitness value determined from the objective function. The main operation in GA is crossover, although mutation plays the role of regenerating lost genetic material by causing sporadic and random alteration of the bits of the binary strings. GAs are generally characterized by their population size, crossover type, crossover rate and elitist selection. These control parameters affect how well the algorithm performs. The optimum set of parameters is dependent on the application being optimized.
PBIL is an abstraction of the GA that explicitly maintains the statistics contained in a GA's population, but abstracts away the crossover operation. PBIL is in fact a combination of evolutionary optimization and hill-climbing. The PBIL algorithm uses a real valued probability vector which, when sampled, reveals a high evaluation solution vector with high probability.
The PBIL algorithm creates a probability vector from which samples are drawn to produce the next generation's population. As in the GA, the solution is encoded into a binary vector of fixed length. Initially the values of the probability vector are set to 0.5. A number of solution vectors, analogous to the population in GAs, are generated based upon the probabilities of the probability vector. The probability vector is pushed towards the generated solution vector with the highest evaluation (fitness value). This probability vector can thus be considered a prototype for high evaluation vectors for the function space being explored. Each bit of the probability vector is updated using:
Pi =[Pi × (1.0 - dl)]+(dl × si) (eqn. 6)
...where Pi is the probability of generating a one in the bit position i, si is the ith position in the solution vector for which the probability vector is being changed and dl is the learning rate. The learning rate is the amount the probability vector is changed after each cycle. A new set of solution vectors is produced after each update of the probability vector. The entries in the probability vector start to drift towards either 0.0 or 1.0 as the search progresses to represent a high evaluation solution vector.
Mutation is used in PBIL for the same reasons as in the GA, to inhibit premature convergence. Mutations perturb the probability vector with a small probability in a random direction. PBILs are generally characterized by their number of samples, learning rate, number of vectors to update from and mutation rate. Fig. 1 shows a flow chart representation of PBIL.
PBIL has been shown to work as well, or better, than the GA. The main advantage of the PBIL over the GA is that since PBIL is characterized by fewer parameters and their values are less problem-related, as little problem specific knowledge as possible is needed.
4. Response surface methodology
When using a numeric field computation based on a FEM analysis, the objective function being optimized will require a lengthy computation time. Owing to this large computation time a faster approach than PBIL is required. Response surface methodology (RSM) is a collection of mathematical and statistical techniques useful for analyzing problems of several independent variables and model exploitation.
By careful design of the FEM experiments, RSM seeks to relate the output variable to the input variables that affect it. The computer experimental result y as a function of the input independent variables x1,x2,...,xn is y = f(x1,x2,...,xn)+ d(x1,x2,...,xn)+ d_ (eqn.7)
... where d(x) is the bias error and d_ is a random error component. If the expected result is denoted by E(y)= S, then the surface represented by S = f(x1,x2,...,xn) is called the response surface.
The form of the relationship between the result and the independent variables is unknown. A polynomial expression is used as a suitable approximation for the true functional relationship between y and the independent variables.
A polynomial expression of degree d can be thought of as a Taylor's series expansion of the true underlying theoretical function f(x) truncated after terms of dth order. The assumption that low order polynomial models can be used to approximate the relationship between the result and the independent variables within a restricted region of the operability space is essential to RSM.
Second order polynomial models are used to model the response surface.
f(x)= bo + bijxixj (eqn. 8)
...where i<j and the coefficients bo, bi, bii and bij are found using the method of least squares. These unknown coefficients can be estimated most effectively if proper computer experimental designs are used to collect the data. Designs for fitting response surfaces are called response surface designs.
4.1 Response surface designs
An experimental design for fitting a second order model must have at least three levels of each factor, so that the model parameters can be estimated.
Rotatable designs are the preferred class of second order response surface designs. The variance of the predicted response at a point x, in a rotatable design, is a function only of the distance of the point from the design center and not a function of direction.
The collecting of the sample results is essential since a sufficiently accurate approximation has to be found with a minimum of experiments. If not all the factorial combinations are employed, the design is called an incomplete factorial design. The Box-Behnken three level design has been chosen for investigating the response surface. This is an incomplete factorial design which is a reasonable compromise between accuracy of the function and the required number of computations. The design generates second order rotatable designs or near-rotatable designs, which also possess a high degree of orthogonality.
4.2 Estimation of errors in response surface fitting
The errors in design of experimental methods are generally divided into two types: (i) systematic or bias errors d(x) which are the difference between the expected value of the response E(y)= S and the approximate objective function f(x), and (ii) the random or experimental error _ in sampling.
In numeric computer experiments, replicated experiments result in the same result, so random errors cannot be defined. Only the bias error from systematic departure of the fixed polynomial from the real response, due to an insufficient polynomial order, can be calculated. An estimate of the variance of the bias error is (eqn.9) where se is the standard error of estimate, m is the number of observations, p is the number of coefficients in the polynomial, yi is the observed response and ˆ yi is the predicted response. The normalized error is:
d = se yo (eqn.10)
where yo =(y1 +y2 +...+yn)/m. The accuracy of the response surface E(y) is varied by changing the size of the investigated region, since only second order polynomials are used.
5. Modern approach to optimization of PM motors
The aim of the optimization procedure is to minimize the cost of the active material of the motor, while ensuring a rated power and high efficiency. It is proposed to optimize a PM motor using the PBIL with RSM. The PBIL does not optimize the performance characteristics directly, but uses polynomial fits of the performance characteristics of the motor, created using RSM.
The motor characteristics are calculated using FEM ( Section 3) and classical machine theory for a number of combinations of input parameters. Fig. 2 shows the logic sequence followed by the control program. The output characteristics and input parameters are used to fit polynomial equations of second order to every output characteristic. These polynomials are used as objective functions and constraints in the PBIL optimization.
The following simplifications are made for the FEM: (a) the end leakage reactance is calculated using classical theory (Appendix A, eqn (A.10)) since it would be very difficult to calculate using a two-dimensional FEM; (b) the induced EMF and inductive reactance are assumed constant throughout the load range, and equal to the value obtained at rated current; (c) the model is independent of rotor position.
5.1 PM d.c. commutator motors
Consider a PM d.c. commutator motor with segmental magnets. The optimization problem can be stated as minimize the cost:
C(x)= VMcPM + Vwcw + Vccc (eqn.11)
…where VM, Vw, Vc are the volume of the PM, copper conductors and steel core, respectively, and cPM, cw, cc are the price per m3 for the PM, copper conductors and electrotechnical steel, respectively. The optimization is subject to the constraints:
5.2 PM synchronous motors
A surface PM rotor and a buried PM rotor synchronous motor have been considered. Only the rotor is being optimized since the stator of an existing induction motor has been used. The objective function attempts to minimize the volume of PM material used. The optimization problem can be expressed as: minimize VM(x) subject to the constraints:
Pelm(x) = Pelm(d) Ja(x) = Jath
?(x) = ?d g = gmin (eqn.14) hM,wM = hMmin Dmax <D1in
where VM is the magnet volume, Pelm(d) is the desired electromagnetic power, Jath is the current density thermal limit at Pelm(d), ?d is the desired electrical efficiency at Pelm(d), g and hM are the mechanical minimum sizes for the air gap and PM respectively, and Dmax is the maximum diameter of the outer edge of the PMs.
The independent variables used in the surface PM rotor design are the air gap g, PM thickness hM and the overlap angle ß, and in the case of the buried PM rotor design the air gap g, PM thickness hM and PM width wM (Fig 3). The motor performance obtained from the FEM, using flux linkage and magnetizing reactance.
The objective function is not easily expressed in terms of the independent variables. RSM is thus not used for modeling the objective function directly, but rather for modeling the performance characteristics of the PM synchronous motor, used in the constraints, in terms of the independent variables.
The characteristics for a number of combinations of the independent variables (factors) are calculated using the FEM. The factors range over the whole problem space and their combinations are determined by the Box-Behnken three level design method, with 3 factors and 3 levels the design needs 15 runs. Response surfaces are created for the output characteristics Ef , Xsd, Xsq, Xad and Xaq. A second order polynomial is fitted to each of these five characteristics using the least squares method. These five polynomial equations are then used to calculate the electromagnetic power, efficiency and stator current in the optimization procedure. The PBIL optimization is used for the minimization of the PM volume with the constraints stated in eqn (eqn.14).
The volume of a surface magnet per pole is VM/(2p)=(0.25ßp/3600) ×(D2 2out - D2 Min)lM, where lM is the magnet axial length, D2out and DMin are magnet diameters according to Fig. 3a. The volume of a buried magnet per pole is VM/(2p)=2hMwMlM (Fig. 3b).
Numerical example 1
Optimization of a surface PM rotor. The stator of a commercial 380 V, 50-Hz, four-pole 1.5-kW induction motor has been used to design a synchronous motor with surface type PM rotor. Table 1 shows the results for different constraints in electromagnetic output power, efficiency and power factor. The results show that at low power ratings (1.5 kW) the design needs a substantial increase in PM material if the power factor is constrained to a minimum of 0.9 at rated power. At the high power rating of 2.5 kW the power factor well exceeds the design minimum.
The most feasible power rating for this motor should ensure that the stator winding current density remains well below the specification limits, but also maximizes the power output. A power rating of 2.2 kW is thus considered appropriate due to the excellent power factor at rated power and the maximum use of the stator windings. The stator used in this motor is from a 1.5 kW induction motor. The increase in power rating to 2.2 kW, when operating as a synchronous motor, is due to the higher stator winding current density possible in PM synchronous machines, and the improved efficiency and power factor. The stator winding current density has increased from 9.93 A/mm^2, for the induction motor, to 10.3 A/mm^2 for the PM synchronous motor (forced ventilation).
Higher time and space harmonics have been neglected and the effects of cogging torque have also been neglected throughout this optimization procedure. The cogging torque can significantly be reduced using an appropriate PM overlap angle. Skewing of the stator teeth is not possible as an existing stator is being used. Using a PM ß =(k +0.14)t1 where k is an integer and t1 = 3600/s1 is the stator slot pitch (angle), the minimum cogging torque can be achieved. The overlap angle is thus decreased from ß =73.80 to 71.40. The air gap and magnet thickness are again optimized for this fixed overlap angle. Table 1 shows the final optimized rotor details.
The performance of this optimized motor is compared with an initial surface PM design. Fig. 4 compares the electromagnetic power and Fig. 5 compares the efficiency for the two motors. The optimized surface PM synchronous motor has superior efficiency at the desired rated power with a reduction in the PM volume from 15,000 mm3 per pole for the initial design to 6838 mm3 per pole for the optimized motor.
Numerical example 1
Numerical example 2
Optimization of a buried PM rotor. The stator of a commercial 380-V, 50 Hz, four-pole 1.5-kW induction motor has been used to design a synchronous motor with buried-type PM rotor. Table 2 shows the results for different constraints in electromagnetic output power, efficiency and power factor.
The most feasible power rating is again chosen to be 2.2 kW. The current density in the stator winding is Ja ˜ 10.4 A/mm2 (forced ventilation).
Numerical example 2
The cogging torque can be minimized in the final design by creating asymmetry of the rotor magnetic circuit. This has a marginal effect on the optimization point since it does not change any of the optimization parameters.
The performance of this optimized motor is compared with an initial buried PM motor design. The optimized buried PM motor uses 6627 mm3 of PM material, while the initial design uses 16, 200 mm3. It has also superior efficiency at rated output power.
The optimization of the surface and buried PM synchronous motors both showed improved performance over their initial designs. The volume of PM material was also reduced in both designs. The buried PM motor design is considered the superior design since it uses the minimum of PM material and has a high efficiency over a wide power range (Fig. 7).
The RSM using the PBIL is thus seen as an appropriate method for optimization of PM synchronous motors with the aid of the FEM. This optimization technique can easily be extended for the optimum design of the whole synchronous motor.