- Email: [email protected]
A new filled function applied to global optimization
Available online at www.sciencedirect.com
Computers & Operations Research 31 (2004) 61 – 80
www.elsevier.com/locate/dsw
A new #lled function applied to global optimization Xian Liua;∗ , Wilsun Xub a
Department of Systems Engineering, University of Arkansas at Little Rock, Little Rock, AR 72204-1099, USA Department of Electrical and Computer Engineering, University of Alberta, Edmonton, Alta, Canada T6G 2G7
b
Received 1 March 2002; received in revised form 1 September 2002
Abstract The (lled function method (FFM) is an approach to #nd the global minimizer of multi-modal functions. The numerical applicability of conventional #lled functions is limited as they are de#ned on either exponential or logarithmic terms. This paper proposes a new #lled function that does not have such disadvantages. An algorithm is presented according to the theoretical analysis. A computer program is designed, implemented, and tested. Numerical experiments on typical testing functions show that the new approach is superior to the conventional one. The result of optimization design for an electrical machine is also reported.
Scope and purpose In the context of mathematical programming, global optimization is concerned with the theory and algorithms on minima of multi-modal functions. In general, global optimization approaches can be classi#ed into two categories: probabilistic and deterministic. The former can usually be applied to general multi-modal functions, whereas the latter typically concentrates on some particular classes of functions. The (lled function method is one of a few deterministic approaches which intend to #nd the global minimum for general multi-modal functions. However, the numerical performance of conventional #lled functions is undesirable as they are de#ned on either exponential or logarithmic terms or multiple parameters. This paper proposes a new #lled function that does not have the above disadvantages. The present work consists of theoretical analysis, algorithm design, computer implementation, mathematical validation, and engineering application. ? 2003 Elsevier Ltd. All rights reserved. Keywords: Filled function method; Global optimization; Mathematical programming; Minimization
∗
Corresponding author. Tel.: +1-780-492-3937; fax: +1-780-492-1811. E-mail addresses: [email protected] (X. Liu), [email protected] (W. Xu).
0305-0548/04/$ - see front matter ? 2003 Elsevier Ltd. All rights reserved. PII: S 0 3 0 5 - 0 5 4 8 ( 0 2 ) 0 0 1 5 4 - 5
62
X. Liu, W. Xu / Computers & Operations Research 31 (2004) 61 – 80
1. Introduction Researches in optimization methods applied to multi-modal functions have been active since the 1970s. Early contributions were reported in the two volumes named towards global optimization [1,2]. The progress since then can be found in [3,4]. In general, the reported approaches can be classi#ed into two categories: deterministic and probabilistic. It is generally agreed that the former concerns the covering method [5,6], the (lled function method [7,8], the trajectory method [9,10], and the tunneling method [11], whereas the latter refers to the clustering method [12,13] and other methods [14–16]. The well-known genetic and simulated annealing algorithms are also commonly regarded as the probabilistic approach. It seems that the (lled function method (FFM) has several advantages over others due to its relatively easy actualization with a process that aims at successively #nding smaller local minima. In the research of global optimization, there are two major issues: (1) How to #nd a lower minimizer of the objective function f(X ) from a known local minimizer of f(X ). (2) How to evaluate the convergence and, accordingly, design the stopping criteria. Since the global minimization concerns the whole region rather than the neighborhood of some particular points, there does not exist a general criterion, like the zero gradient criterion for the local minimizer, for the global minimizer. Although a great deal of eEorts has been made to #nd a general criterion to characterize the global minimizer, none is constructive. For this reason, so far none has been accepted to solve practical problems. On the other hand, many global minimization approaches introduced restrictions to the objective functions, such as global convexity, polynomial structure, rational features, etc. Without these restrictions, it is very diFcult, if not impossible, to conduct the convergence analysis. The situation is similar to the stopping criteria. Unlike other deterministic approaches, the FFM is a technique to #nd the global minimizer for general functions and does not assume any particular features of the objective function. In the literature, all researches on the FFM have mainly been concerning issue (1) above. In its nature, the FFM is a deterministic approach in the sense of issue (1). However, just like any other deterministic approaches for the global optimization, if no any restriction to f(X ) is enforced, the stopping criteria cannot avoid involving probabilistic strategies. Therefore, the research on stopping criteria, as well as the convergence analysis, has evolved to a specialistic area. This paper focuses on the #rst issue mentioned above. The rest of this paper is organized as follows. An overview of the FFM is given in Section 2. Then a class of new #lled functions, called the M -function, is de#ned and analyzed in Section 3. Next, in Section 4, an algorithm is presented and results of numerical experiments are reported. The practice for solving an engineering problem, optimization design of electrical machines, is addressed in Section 5. Conclusions are included in Section 6. Finally, the testing functions used for numerical experiments are described in Appendix A and the modeling issues of electrical machines are detailed in Appendix B. 2. Essentials of the FFM The FFM is an approach to #nd the global minimizer of a multi-modal function f(X ) on Rn , under the assumptions that f(X ) is continuously diEerentiable and f(X ) → +∞ as X → +∞.
X. Liu, W. Xu / Computers & Operations Research 31 (2004) 61 – 80
63
Thus, there exists a closed bounded domain that contains all global minimizers of f(X ). In the following, for the presentation purpose, will be called the operating region of the concerned global optimization problem. The FFM also assumes that f(X ) has only a #nite number of minimizers. To introduce the concept of the FFM, let us de#ne the following concepts: Denition 2.1. A basin of f(X ) at an isolated minimizer X1 is a connected domain B1 which contains X1 and in which starting from any point the steepest descent trajectory of f(X ) converges to X1 , but outside which the steepest descent trajectory of f(X ) does not converge to X1 . Denition 2.2. A hill of f(X ) at X1 is the basin of −f(X ) at its minimizer X1 , if X1 is a maximizer of f(X ). Denition 2.3. A local minimizer X2 is said to be higher than X1 if and only if f(X2 ) ¿ f(X1 ), and, for this case, B2 is said to be a higher basin than B1 . In this paper, Bh and Bl denote all higher and lower basins than current basin B1 of f(X ), respectively. Denition 2.4. A function P(X ) is called a #lled function of f(X ) at X1 if (1) X1 is a maximizer of P(X ) and the whole basin B1 becomes a part of a hill of P(X ); (2) P(X ) has no stationary points in any Bh s; and (3) There is a point X in a Bl (if such a basin exists) that minimizes P(X ) on the line through X and X1 . In this paper, we will allow an in#nite maximizer of P(X ). The FFM consists of two phases, local minimization and #lling: Phase 1: In this phase, a local minimizer X1 of f(X ) is found. Any eEective down-hill technique, for instance, the variable metric method, can be employed in phase one. Phase 2: In this phase, a #lled function is constructed and is minimized by the down-hill procedure. Phase 2 ends when such an Xs is found that Xs is in a Bl . Then, the FFM reenters phase 1, with Xs as the starting point, to #nd a new local minimizer X2 of f(X ) (if such one exists), and so on. The above process is repeated until the global minimizer is found. Several #lled functions have been proposed in the literature. Three popular ones are [7,8]: P(X; r; ) = exp(−X − X1 2 =2 )=[r + f(X )]; G(X; r; ) = −{2 ln[r + f(X )] + X − X1 p }; Q(X; a) = −[f(X ) − f(X1 )] exp(aX − X1 p );
(1)
where p = 1 or 2, r and are adjustable parameters, and a is an adjustable positive weight factor. Both P- and G-functions require two adjustable parameters, which need to be appropriately iterated and coordinated each other, hence their algorithmic realization is fairly complicated. For this reason, it is usually agreed that the Q-function is somehow better than the other two, since it involves only one adjustable parameter. However, the Q-function includes an exponential term whose argument
64
X. Liu, W. Xu / Computers & Operations Research 31 (2004) 61 – 80
is the product of the weight factor a and the norm. As a becomes larger and larger, as required to preserve the #lling property, the rapid increasing value of the exponential term will result in failure of computation even if the size of the feasible region is moderate. In practice, this kind of ill-conditioning problem frequently occurs. To make the Q-function work, many additional cares must be incorporated into the algorithm. It is obvious that the exponential term in the Q-function has seriously limited its applicability to the practical global optimization problems, especially those raised from engineering. 3. The M -function The M -function is de#ned as 1 M (X; a) = − aX − X1 2 ; (2) (f − f1 )1=m where f − f1 stands for f(X ) − f(X1 ), m is a natural number greater than 1, and a is a positive real used as the weight factor. For simplicity M (X; a) will be denoted as M (X ) throughout the paper. The #lling properties of M (X ) are exhibited by the following theorems. Theorem 3.1. Given d ∈ Rn and f(X ) ¿ f(X1 ), if dT ∇f(X ) ¿ 0 and dT (X − X1 ) ¿ 0, or dT ∇ f(X ) ¿ 0 and dT (X − X1 ) ¿ 0, then d is a descent direction of M (X ) at point X . Proof. It follows from (2) that dT ∇f(X ) T T d ∇M (X ) = − + 2ad (X − X1 ) : m(f − f1 )1+1=m Therefore, the conditions given guarantee dT ∇M (X ) ¡ 0. Consider the following de#nition: dT ∇f(X ) al (X ) = − 2m(f − f1 )1+1=m dT (X − X1 ) then there is another result: Theorem 3.2. If f(X ) ¿ f(X1 ), dT ∇f(X ) ¡ 0, and dT (X − X1 ) ¿ 0, then the sign of dT ∇M (X ) is that of al (X ) − a. Note that it is possible to have al (X )−a ¿ 0, because al (X ) → +∞ as f(X ) ¿ f(X1 ) and f(X ) → f(X1 ). These results clearly characterize the #lling property of the M (X ). Particularly, in the ascent region of the current basin B1 or a higher basin than B1 , d is always a descent direction of M (X ). On the other hand, in the descent region of a higher basin than B1 , d is still a descent direction of M (X ) provided that the weight factor a is suFciently large. Furthermore, in a lower basin than B1 , d may become an ascent direction of M (X ) and this possibility does exist. Therefore, under the continuously diEerentiable assumption for f(X ), M (X ) must have a stationary point along d.
X. Liu, W. Xu / Computers & Operations Research 31 (2004) 61 – 80
65
3.0 2.5 M(X), a=0.004, m = 2 2.0 1.5 1.0 f(X), M(X)
f(X) 0.5 0 −0.5 −1.0 M(X), a=0.03, m = 2
−1.5 b1
−2.0 2
4
b2
X1 6
8
b3 10
b4 12
14
b5 16
b6 18
20
22
X
Fig. 1. The behavior of M (X ), diEerent a.
The behavior of M (X ) may be well interpreted graphically when X is a one-dimensional variable. In this case, there are only two possible search directions, the positive and negative directions of the X -axis. Consider an arbitrary multi-modal function f(X ) that has a local minimizer X = X1 (Fig. 1). When X ∈ (X1 ; b2 ), M (X ) descends along the positive X -axis. When X ∈ (b1 ; X1 ), M (X ) descends along the negative X -axis. Let ˜i be the unit vector that points to the positive direction of the X -axis. Then the above features can be rephrased as follows: M (X ) always descends along the direction (X − X1 )˜i if X = X1 . The situation is similar in interval (b3 ; b4 ). This property holds for any positive weight factor a, according to the theory of M (X ). However, when X ∈ (b2 ; b3 ), weight factor a must be suFciently large to make M (X ) descending along the positive X -axis. On the other hand, when X ∈ (b4 ; b5 ), M (X ) will have a stationary point and the second phase may stop at that point. The global minimizer b6 can be reached by resuming phase 1 with b5 as the starting point. The eEect of m is shown in Fig. 2. From its de#nition, the M -function appears more applicable to computational assignments, because it retains such advantages of conventional #lled functions that the exponential terms and the second parameter are not needed. Other desirable features include that the formulation is of the sum of two terms, rather than a product. Also, the M -function has only the linear dependency on the norm, although the norm itself may be quadratical, whereas in the P- and Q-functions, the norm is used as the argument of a exponential function. In addition, as described below,
66
X. Liu, W. Xu / Computers & Operations Research 31 (2004) 61 – 80 2.5 M(X), a=0.01, m=2 2
1.5 M(X), a=0.01, m=3
1 f(X) f(X), M(X)
0.5
0 −0.5 −1 −1.5 b1
−2 2
4
X1
b2 6
8
b3 10
b4 12
14
b6
b5 16
18
20
22
X
Fig. 2. The behavior of M (X ), diEerent m.
the lower bound of weight factor a of the M -function is usually smaller than the one of the Q-function. The weight factor a plays a crucial role in a #lled function. Theoretically, the value of a must be suFciently large to preserve a desirable #lling capability. Computationally, the value of a should be small to make the numerical procedures healthy, because there is a product of a and the norm in the formulation. Therefore, a #lled function is a robust one if the value of al (X ) is small given a particular X . In the following, we compare the M -function with the Q-function in terms of the lower bound of a. Consider the de#nition below aq (X ) = −
dT ∇f(X ) : 2(f − f1 )dT (X − X1 )
It is easy to show from (1) that, with the same conditions as given in Theorem 3.2, if a ¿ aq (X ), then d is a descent direction of Q(X ) at point X . Furthermore, we have the following theorem: Theorem 3.3. Under the same conditions as given in Theorem 3.2, if f(X ) ¿ f(X1 ) + m−m , then al (X ) ¡ aq (X ) at X . Proof. f(X ) ¿ f(X1 ) + m−m leads to m(f − f1 )1=m ¿ 1. This in turn implies al (X ) ¡ aq (X ).
X. Liu, W. Xu / Computers & Operations Research 31 (2004) 61 – 80
67
The above theorem is important as it presents the fact that, even with a small weight factor a, the M -function preserves the desired #lling property as long as a proper m is chosen. We give some remarks to end this section. One might argue that M -function has a weakness because of the root term, such that it has no de#nition in the region corresponding to f(X ) ¡ f(X1 ) if m is even. This feature, however, does not aEect its application, because in the numerical iteration process of minimization, we always check the function value of current iteration point #rst. If at some Xs we obtain f(Xs ) ¡ f(X1 ), then Xs is in a lower basin than B1 already. Phase 2 can ends right here and the FFM can reenter phase 2 by starting from Xs . Therefore, the root term in the formulation of M -function is never an obstacle in the arithmetic realization of the FFM. Moreover, an odd m may be chosen if the above problem really needs to be avoided. Yet another concern might be the situation at such an iterating point X that the value of f(X ) − f(X1 ) is positive but very small, resulting in a relatively large ∇M (X ). As a matter of fact, this feature is just what is needed for any #lled functions: in the neighborhood of X2 such that {X2 = X | f(X ) = f(X1 ) and X2 = X1 }, a steep “barrier” is necessary in order to stop phase 2, as revealed by Theorem 3.2. This feature is diEerent in nature than the steep behavior inherited in the Q-function due to the eEect of the sequentially increased exponential term.
4. An algorithm and the numerical experiment In optimization studies, modeling, algorithmization, and numerical realization are three relevant yet independent aspects. It has been well known that the numerical realization of an established formulation for nonlinear programming usually involves more diFculties. The presented FFM is not an exception and even brings an additional layer of complexity. In the modeling aspect of local optimization studies, the function to be minimized is usually implicitly assumed to have the globally convex property: f(X ) → +∞ as X → +∞. Although in theory it is not impossible to develop a globally convex #lled function, there remain several analytical diFculties in the modeling aspect. As a matter of fact, almost all #lled functions reported so far are globally concave, including the M (X ) proposed in the present work. As a result, numerically speaking, many well-developed local minimization techniques may not be eEective in phase 2 of the FFM. In general, it seems that the algorithms not relying on the gradient, Hessian matrix, or conjugated directions are more appropriate. In early studies on the FFM, several local minimization methods, which compare the numerical values of the function only, were investigated and customized to minimize the #lled function. As demonstrated from extensive numerical experiments for various applications, the method adapted in the present work, the linear search approach, has been giving the most expected results. In this section, we present an algorithm for the n-dimensional case. The involved nomenclature in the algorithm is listed as follows: a: fmin : fmin0 : : L:
weight factor in the formulation of the #lled function value of the objective function at the local minimizer initial value of fmin , used for the iteration purpose. It is advisable to select fmin0 as a large positive value operating region of the concerned global optimization problem (de#ned in Section 2) counter recording the number of linear search operations in phase 2 of the FFM
68
X. Liu, W. Xu / Computers & Operations Research 31 (2004) 61 – 80
Table 1 Numerical experiments M -function
Six-hump camel-back Branin G–P Rastrigin Shubert III Sine-square I Sine-square II Sine-square III
Lm : X0 : #:
Q-function
G-function
Kf
Kt
Kf
Kt
Kf
Kt
204 97 354 226 406 784 113 997
346 252 498 572 569 2387 1716 3564
204 97 280 118 NF NF NF NF
346 3640 1059 3589
204 97 261 117 444 1063 113 1562
346 252 404 619 685 2665 1716 3784
speci#ed upper bound of L. An estimation is Lm ≈ 2n dmin D Gmax , where dmin is the smallest search step, D is a metric of , and Gmax = supx∈ ∇f(X ). initial point to start phase 1 of the FFM (local minimization) multiplier to update the weight factor a (# ¿ 1, usually taking 2.0)
The algorithm is detailed as follows: Step 0: Specify X0 , a, and Lm . fmin0 → fmin . Step 1: 0 → M . Enter phase 1 of the FFM. Activate the minimization procedure to minimize the objective function f(X ), starting from X0 . Find a local minimizer X1 . If f(X1 ) 6 fmin , then f(X1 ) → fmin ; otherwise. #a → a. Step 2: Enter phase 2 of the FFM. Step 3: Construct a search direction s using a random number generating procedure. L + 1 → M . Starting from X1 , activate the linear search procedure to minimize the #lled function. Step 4: Continue the linear search. Arrive at point X . Step 5: 1. If X ∈ , then go to Step 6; 2. If f(X ) ¡ f(X1 ), then X → X0 ; go to Step 1. 3. If X is the minimizer of the #lled function, then X → X0 ; go to Step 1. 4. Otherwise, go to Step 4. Step 6: If L ¡ Lm , then go to Step 2; otherwise taking the smallest minimum as the global one. Stop. Several testing functions were proposed in literature to evaluate the numerical performance of a new approach for global optimization. To estimate its eEectiveness, the M -function has been tested, with m = 2, with a set of testing functions detailed in Appendix A. Two conventional #lled functions, the Q- and G-functions, have also been tested. Identical starting points are selected for all #lled functions. The results of numerical experiments are presented in Table 1, where Kf is the total number of evaluations for the objective function and the #lled function when the global minimizer is found; Kt the total number of evaluations for the objective function and the #lled function when the algorithm terminated; and NF the failed to #nd the global minimizer.
X. Liu, W. Xu / Computers & Operations Research 31 (2004) 61 – 80
69
The evaluation of an algorithm or a formulation of the #lled function may involve several layers of the concerned numerical procedures. It is believed, however, that the number of iterations should not be regarded as an appropriate index [17]. This is because, in the down-hill search process, the number of iterations is dependent on what particular methods are used for the lower layers, such as unconstrained minimization or linear search. Typically, there may be several function evaluations per iteration, and it is these evaluations that consume the major part of total execution time in practical problems. For this reason, we used the total number of evaluations of the objective function and the #lled function as an appropriate measure for the performance of diEerent #lled functions. Note that Kt depends not only on Lm , but also on several implementation issues such as the particular linear search algorithm. For example, the number of evaluations of the objective function per linear search is generally diEerent with the Fibonacci, gold ratio, parabolic, or spine line interpolation procedure. In the present work, the parabolic interpolation procedure is used. In testing, identical initial weight factor and multiplier were used for all #lled functions. It is noted that the results of numerical experiments shown in Table 1, in terms of Kf and Kt , imply that the M -function is superior to the Q-function, especially, for the complicated functions like the Shubert III function, or the high-dimensional functions like the Sine-square functions. On the other hand, we also included the G-function as a reference. Strictly speaking, it is inappropriate to compare the G-function with the Q-function or the M -function, because the former has two adjustable parameters whereas the latter have only one. This means that the G-function has an extra degree of freedom. Nevertheless, it might be interesting to see how the G-function behaves, since no testing results for the G- or Q-function were reported in the original paper [7]. To make the G-function have the same degree of freedom as the G- or Q-function does, we used r = f(X1 ) and 2 = 1=a, as implied in [7]. For this particular setting, the M -function is superior to the G-function. We also observed that the value of a had little inNuence to the M - or G-functions except it had been frequently updated when the Shubert III function or the Sine-square functions was tested, because of the large number of local minimums or high dimensionality of these testing functions. This was not the case, however, for the Q-function. Since the Q-function includes an exponential term, in testing a moderate value a often resulted in overNow when the current iterative point X was far away from X1 , while the M - or G-functions still worked well in the same situation. For example, in testing the Q-function failed when a = 4:4, while the M - or G-functions worked well with a = 35:2. Therefore, we conclude that the main disadvantage of the Q-function is its exponential term. Perhaps this is why no any testing was reported in the original paper [7]. In summary, according to the theoretical analysis and numerical testing, it is reasonable to expect that the M -function will be well applicable to the usually more complicated engineering optimization programs. This is actually the practice described in the next section. 5. Optimization design of an electrical machine In the last decade there was a growing interest in design and manufacturing optimization of electrical machines due to increased cost of materials plus the cost of losses on their assumed lifetime [18–22]. A motor’s design can be described by a vector X of n components stating physical characteristics, such as dimensions of core and windings, as well as electromagnetic characteristics, such as densities of current and Nux. The design is subject to a set of m constraints that
70
X. Liu, W. Xu / Computers & Operations Research 31 (2004) 61 – 80
may include speci#cations arising from electrical, mechanical, thermal, manufacturing or standards limits. The task of the design optimization is to make a speci#ed objective function f(X ) reach its optimal value while keeping other technical indices within acceptable ranges. The complexity of electromagnetic characteristics incurred in motor design almost exclusively leads to multi-modal objective functions that have many local minimizers. The well-performed procedures for #nding a local minimizer usually lose their eFciency when applied to #nd the global minimizer of multi-modal functions. Therefore, it is really desirable to employ an approach based on the global optimization methodology to solve the concerned design problems. As examples to demonstrate the applicability of the proposed M -function to practical problems, optimized designs are produced for a single-phase induction motor. The original version of the FFM was developed for solving unconstrained problems under the assumption that operating region is known, where is a closed bounded domain that contains all global minimizers of f(X ). Since the design of electrical machines is a constrained optimization problem, our #rst challenge is to customize the FFM to the constrained environment. This process consists of two issues: con#guration and transformation. The former involves how to establish , whereas the latter deals with the constraints associated with the established . We discuss the con#guration issue #rst. By de#nition, there is a feasible region D in any constrained problem. Any solutions outside D are not useful. Therefore, for a motor design problem, we just look for its global optimal point inside D. This implies that the FFM only needs to be active throughout D. A natural approach seems to let coincide with D. This approach, however, has pitfalls. Recall that the FFM consists of two phases. In phase 1, a local minimizer X1 is found. At the beginning of phase 2, the M -function is constructed and it has an impulse at X1 . Then a down-hill search process is activated in the neighborhood of X1 . Accordingly, the position of X1 will directly aEect the behavior of phase 2. For motor design models that take cost as objective function subject to electromagnetic performance limits or vice versa, the local optimal design is usually located on or near boundaries of feasible region D. If let = D, then the impulse is just on or near one of the boundaries of . Consequently, starting near X1 the down-hill iterations intend to leave very soon and terminate one cycle of phase 2 without any meaningful results. For this reason, An applicable con#guration should have D ⊂ , and this depends on what particular technique is employed in phase 1 to integrate the constraints. The transformation issue cares about handling the constraints in optimization. One of the eEective approaches for solving constrained problems is the indirect method (IDM), which uses the augmented objective function rather than the original one. A typical augmented objective function consists of the weighted sum of all constraint terms and the original objective function. Here, the constraint term is de#ned to be a transformation of the constraint of interest. For example, one of the most popular IDMs, the exterior penalty function method (EPFM), takes the following form [23]: m F(X; r) = f(X ) + r {min[0; gi (X )]}2 ; i=1
where r is the penalty factor with an adjustable positive value. In the case of optimal points located on the boundaries of D, some IDMs may introduce a deviation. For instance, the minimizer found with the (EPFM) are outside D when r is reasonably large. As a result, a new region DF can be formed by enlarging region D with a small increment %.
X. Liu, W. Xu / Computers & Operations Research 31 (2004) 61 – 80
71
With an appropriate % value, all of the minima of F(X; r) are contained in DF . Therefore, we may make = DF . Note that D ⊂ has been satis#ed. Also, a small % is unlikely to lead to too much unnecessary search work. Based on the forgoing analysis, it is possible to use the convergence error &p to determine %. In the EPFM case, one of the most popularly adopted convergence criteria is the closeness measure of the minimizer to the boundaries. For instance, phase 1 can be terminated if gi (X ) ¿ − &p , where i = 1; 2; : : : ; m, and &p ¿ 0. Thus it is reasonable to choose % ¿ &p . The present work has integrated the EPFM into the FFM. For the single-phase induction motor of interest, the involved nomenclature is listed as follows: Lc C Dm Da Zm Za W H F X0 X∗ Cum Cua Al Si cm ca csi cal
length of core (cm), capacity of the operating capacitor (F), diameter of conductor in main winding (mm), diameter of conductor in auxiliary winding (mm), turns of main winding, turns of auxiliary winding, width of rotor end-ring (cm), height of rotor end-ring (cm), value of the objective function (cost of eEective materials), value of original design variables, value of optimal design variables, amount of copper of the main winding, amount of copper of the auxiliary winding, amount of aluminium of the rotor winding, amount of silicon sheet of the core, unit cost of copper of the main winding, unit cost of copper of the auxiliary winding, unit cost of silicon sheet, unit cost of aluminium,
The objective function is chosen as the cost of eEective materials and formulated as follows: f(X ) = Cum cm + Cua ca + Sicsi + Alcal : The details of electromagnetic modeling of concerned motors are included in Appendix B, since this section focuses on how to customize a theoretical optimization method and integrate it into a well-developed electromagnetic model. In the concerned design, eight variables are chosen to adjust the following parameters: length of core, operating capacitor, diameters of conductor in main and auxiliary windings, turns of main and windings, and dimensions of the rotor end-ring. The speci#ed constraints included the output power, eFciency (%), ratio of backward rotating #eld to forward rotating #eld (%), and dimensions of the rotor end-ring. Two sets of optimized results, D1 and D2, are shown in Tables 2 and 3, respectively. The original designs themselves had been produced by a team of experienced engineers using a computer program of motor analysis. Many diEerent combinations of design variables had been carefully selected to obtain a good solution. Therefore, the original designs of this motor were sub-optimal solutions
72
X. Liu, W. Xu / Computers & Operations Research 31 (2004) 61 – 80
Table 2 Design D1
X0 X∗
Lc
C
Dm
Da
Zm
Za
W
H
1:80 2:99
2:00 2:20
0:150 0:135
0:140 0:131
9647 4956
6711 4153
0:353 0:233
0:675 0:554
Lc
C
Dm
Da
Zm
Za
W
H
1:80 3:06
2:00 2:36
0:160 0:131
0:150 0:132
9647 5009
6711 4056
0:300 0:246
0:600 0:548
Table 3 Design D2
X0 X∗
Table 4 Minimizers of design D1 Index of minimizers Object function F
1 1:0000
2 0:9594
3 0:9490
4 0:9416
5 0:9240
6 0:9259
Index of minimizers Object function F
7 0:9218
8 0:9174
9 0:8963
10 0:8733
11 0:8708
12 0:8508
Index of minimizers Object function F
13 0:8668
14 0:8270
15 0:8140
16 0:8039
17 0:7988
Index of minimizers Object function F
1 1:0000
2 0:9885
3 0:9495
4 0:9270
5 0:9022
6 0:9079
Index of minimizers Object function F
7 0:8893
8 0:8761
9 0:8570
10 0:8399
11 0:8272
12 0:8166
Table 5 Minimizers of design D2
already. As an eEort to solve a practical problem, the M -function is applied to improve the original designs. The results are listed in Tables 4 and 5, where the objective function is normalized on the value of the original designs. For design D1, up to 17 local minimizers were found. while for design D2, the total number of identi#ed local minimizers is 12. The merit is obvious: The global optimization program made about 20% reduction in the cost of both designs while satisfying all speci#ed constraints. Thus the M -function indeed exhibited the desirable power. This is a signi#cant result since, especially, the starting points correspond to previously screened designs.
X. Liu, W. Xu / Computers & Operations Research 31 (2004) 61 – 80
73
6. Conclusions The #lled function method (FFM) is an approach to #nd the global minimizer of multi-modal and multidimensional functions. Conventional functions, however, rely on either an exponential term or multiple parameters. In numerical applications, these characteristics are strongly undesirable. In this paper, a new class of #lled functions called the M -function is proposed. The M -function has several desirable features such as (1) no any exponential or logarithmic terms; (2) only one parameter; and (3) formulated as a sum of two terms, rather than a product. The eEectiveness of the M-function applied to practical problems is exhibited by results of numerical experiments on typical testing functions. A computer program is designed based on the new approach and has been successfully applied to the optimization design of electrical machines. Appendix A. Testing functions Six-hump camel-back (n = 2) [9]: fC (X ) = 4x12 − 2:1x14 + x16 =3 + x1 x2 − 4x22 + 4x24
(−3 6 x1 ; x2 6 3):
The global minimizers are (0:08983; −0:7126) and (−0:08983; 0:7126). Branin (n = 2) [24]: fB (X ) = (x2 − 1:275x12 =,2 + 5x1 =, − 6)2 + 10(1 − 0:125=,) cos(x1 ) + 10 (−5 6 x1 6 10; 0 6 x2 6 15): The global minimizers are (−3:142; 12:275), (3:142; 2:275), and (9:425; 2:425). Goldstein–Price (G–P) (n = 2) [2]: fG (x) = [1 + (x1 + x2 + 1)2 (19 − 14x1 + 3x12 − 14x2 + 6x1 x2 + 3x22 ] ×[30 + (2x1 − 3x2 )2 (18 − 32x1 + 12x12 + 48x2 − 36x1 x2 + 27x22 )]
(−3 6 x1 ; x2 6 3):
The global minimizer is (0; −1). Rastrigin (n = 2) [25]: fR (X ) = x12 + x22 − cos(18x1 ) − cos(18x2 )
(−1 6 x1 ; x2 6 1):
This function has about 50 minimizers. The global minimizer is (0, 0). Shubert III (n = 2) [8]: 5 5 i cos[(i + 1)x1 + 1] i cos[(i + 1)x2 + 1] fS (X ) = i=1
i=1 2
+ [(x1 + 1:42513) + (x2 + 0:80032)2
(−10 6 x1 ; x2 6 10):
This function has 760 minimizers. The global minimizer is (−1:42513; −0:80032). Because of the large number of local minimizers and the steep slope around the global minimizer, the Shubert
74
X. Liu, W. Xu / Computers & Operations Research 31 (2004) 61 – 80
300 250 200 150 100 50 0 −50 −100 −150 10 10
5 5
0 x2
0
−5
−5 −10
x1
−10
Fig. 3. The Shubert function III.
function III has widely been recognized as an important testing function. An illustration is given in Fig. 3. Sine-square I (n = 6) [17]: n− 1 2 2 2 2 fsq1 (X ) = 10 sin (,x1 ) + (x n − 1) + (xi − 1) [1 + 10 sin (,xi+1 )] ,=n i=1
(−10 6 xi 6 10): This function has about 60 minimizers. The global minimizer is (1, 1, 1, 1, 1, 1). Sine-square II (n = 6) [17]: fsq2 (X ) =
2
2
10 sin (,y1 ) + (yn − 1) +
n− 1 i=1
yi = 1 + (xi − 1)=4
(−10 6 xi 6 10):
This function has about 30 minimizers. The global minimizer is (1, 1, 1, 1, 1, 1).
2
2
(yi − 1) [1 + 10 sin (,yi+1 )] ,=n;
X. Liu, W. Xu / Computers & Operations Research 31 (2004) 61 – 80
75
Sine-square III (n = 6) [17]: fsq3 (X ) = sin2 (3,x1 ) + (x n − 1)2 [1 + sin2 (2,x n )]
+
n− 1
2
2
(xi − 1) [1 + sin (3,xi+1 )]
10
(−10 6 xi 6 10):
i=1
This function has about 180 minimizers. The global minimizer is (1, 1, 1, 1, 1, 1).
Appendix B. An electromagnetic model of single-phase induction (SPI) motors Most SPI motors use an auxiliary winding with a capacitor to produce the start torque, since the magnetic #eld generated by a single-phase source remains stationary in space and pulses with time. In the following, the main winding and auxiliary winding are denoted as the M and A winding, respectively. One of the representative classes of SPI motors is the capacitor-run type, in which the A winding remains being connected to the source after starting. The present work is devoted to capacitor-run SPI motors. SPI motors in nature are asymmetrical two-phase motors, because both the eEective turns and the current magnitude of the M winding are usually diEerent than those of the A winding, and the phase diEerence between the two currents is not equal 90 electrical degrees. Most electromagnetic models of SPI motors reported in the literature can be considered being variants of the one described in [26], where the equivalent circuits of SPI motors are derived from the symmetrical component method. The circuits corresponding to the forward and backward phase sequences of the M winding and the A winding are shown in Fig. 4(a) – (d), respectively. In the present work, an electromagnetic analysis procedure is developed based on these equivalent circuits. rM1 XM1
XmM
ZM+
rM1 XM1
X’M2 r’M2
X’M2 XmM
ZM-
(a)
rA1 XA1 ZA+
(b)
rA1 XM1
X’A2 XmA
(c)
r’M2 2-S
S
r’A2 S
X’A2 XmA
ZA-
r’A2 2-S
(d)
Fig. 4. The equivalent circuits of SPI motors.
76
X. Liu, W. Xu / Computers & Operations Research 31 (2004) 61 – 80
The principal variables and parameters used in the electromagnetic analysis procedure are listed as follows: cos . power factor Ist ratio of starting current to rated current rA 1 resistance of the A winding rA 2 rotor resistance referred to the A winding rM 1 resistance of the M winding rM rotor resistance referred to the M winding 2 xA 1 leakage reactance of the A winding xA 2 rotor leakage reactance referred to the A winding xM 1 leakage reactance of the M winding xM rotor leakage reactance referred to the M winding 2 xmA magnetizing reactance of the A winding xmM magnetizing reactance of the M winding xa xmM + xM 2 Ka xmM =xa n rotating speed (rpm) pFe1 stator core losses pFe2 rotor core losses pmech friction and windage losses ps stray losses P1 input power P2 output power U voltage of the single phase power supply S slip Sf A slot space #lling factor of the A winding Sf M slot space #lling factor of the M winding Tm ratio of maximum torque to rated torque Tst ratio of torque to rated torque 4 eFciency B.1. Formulas for objective function and constraints The objective function and constraints in the design procedure are expressed in terms of a set of analytical formulas. According to Fig. 4(a), the forward sequence impedance of the M winding is ZM+ = rM1 + jxM1 + 2Zf : The 2Zf is the load impedance: 2Zf =
+ r =S) jxmM (jxM M2 2 ); rM2 =S + j(xmM + xM 2
= 2(Rf + jXf );
X. Liu, W. Xu / Computers & Operations Research 31 (2004) 61 – 80
77
where Rf =
=S 0:5Ka2 rM 2 ; =(Sx )]2 1 + [rM a 2
Xf =
=S)2 + x 0:5xmM (rM mM xa xM2 2 : (rM2 =S)2 + xa2
Similarly, from Fig. 4(b), we have ZM− = rM1 + jxM1 + 2Zb : The 2Zb is the load impedance: 2Zb =
+ r =(2 − S)] jxmM [jxM 2 M2 =(2 − S) + j(x rM2 mM + xM2 )
= 2(Rb + jXb ); where Rb =
=(2 − S) 0:5Ka2 rM 2 ; 1 + {rM2 =[(2 − S)xa ]}2
Xb =
=(2 − S)]2 + x 0:5xmM [rM mM xa xM2 2 : =(2 − S)]2 + x 2 [rM a 2
For the convenience of programming, in the design procedure a set of intermediary parameters is introduced to derive the armature currents, power factor, and the input power: RTM = rM1 + Rf + Rb ; XTM = xM1 + Xf + Xb ; RTA = rA1 + a2 (Rf + Rb ) + rc ; XTA = xA1 + a2 (Xf + Xb ) − xc ; Rfb = Rf − Rb ; Xfb = Xf − Xb ; 2 ); R3 = RTM RTA − XTM XTA − a2 (R2fb − Xfb
X3 = RTM XTA + XTM RTA − 2a2 Rfb Xfb ; R4 = R3 U=(R23 + X32 );
78
X. Liu, W. Xu / Computers & Operations Research 31 (2004) 61 – 80
X4 = R4 X3 =R3 ; R5 = R4 (RTA − aXfb ) + X4 (XTA + aRfb ); X5 = R4 (XTA + aRfb ) − X4 (RTA − aXfb ); R6 = R4 (RTM + aXfb ) + X4 (XTM − aRfb ); X6 = R4 (XTM − aRfb ) − X4 (RTM + aXfb ): The M winding current, the A winding current, and the total input current are: IM = R25 + X52 ; IA = R26 + X62 ; I1 =
(R5 + R6 )2 + (X5 + X6 )2 ;
respectively. The power factor is cos . = (R5 + R6 )=I1 : The input power is P1 = U (R5 + R6 ) + pFe1 : The air-gap power (transformed from stator to rotor) is 2 Pg = (IM + a2 IA2 )Rfb + 2a(Rf + Rb )(R5 X6 − R6 X5 ):
Finally, the output power is determined as follows: P2 = Pg (1 − S) − pmech − pFe2 − ps : B.2. The design procedure Given the speci#cations (ratings and performance constraints) and the core dimensions, a typical design procedure for SPI motors consists of the following steps: 1. 2. 3. 4. 5. 6. 7. 8.
derive the dependent core data; design the main winding; determine the resistance and reactances of the main winding; estimate the magnetomotive forces in multiple core regions; estimate the iron, mechanical, and stray losses; design the auxiliary winding; calculate the full-load performance; estimate the starting performance;
X. Liu, W. Xu / Computers & Operations Research 31 (2004) 61 – 80
79
9. determine the maximum torque; and 10. evaluate the cost of eEective materials. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24]
Dixon L, SzegRo G, editors. Towards global optimization. Amsterdam: North-Holland, 1975. Dixon L, SzegRo G, editors. Towards global optimization, vol. 2. Amsterdam: North-Holland, 1978. Horst R, Tuy H. Global optimization (deterministic approaches), 3rd ed. Berlin: Springer, 1996. S TRorn A, Zilinskas A. Global optimization. Berlin: Springer, 1989. Basso P. Iterative methods for the localization of the global maximum. SIAM Journal on Numerical Analysis 1982;19:781–92. Mladineo R. An algorithm for #nding the global maximum of a multimodal, multivariate function. Mathematical Programming 1986;34:188–200. Ge R, Qin Y. A class of #lled functions for #nding global minimizers of a function of several variables. Journal of Optimization Theory and Applications 1987;54:241–52. Ge R. A #lled function method for #nding a global minimizer of a function of several variables. Mathematical Programming 1990;46:191–204. Branin F. Widely convergent methods for #nding multiple solutions of simultaneous nonlinear equations. IBM Journal of Research Developments 1972;504 –22. Snyman J, Fatti L. A multi-start global minimization algorithm with dynamic search trajectories. Journal of Optimization Theory and Applications 1987;54:121–41. Levy A, Montalvo A. The tunneling algorithm for the global minimization of functions. SIAM Journal on Scienti#c and Statistical Computing 1985;6:15–29. Becker R, Lago G. A global optimization algorithm. In: Proceedings of the Eighth Allerton Conference on Circuits and Systems Theory, Monticello, IL, USA, 1970. p. 3–12. TRorn A. Cluster analysis using seed points and density determined hyperspheres as an aid to global optimization. IEEE Transactions on Systems, Man, and Cybernetics 1977;7:610–6. Archetti F. A probabilistic algorithm for global optimization problems with a dimensionality reduction technique. In: Balakrishnon A, Toma M, editors. Lecture notes in control and information sciences, vol. 23. Berlin: Springer, 1980. p. 36 – 42. Baba N. Global optimization of functions by the random optimization method. International Journal of Control 1979;30:1061–5. Dorea C. Limiting distribution for random optimization methods. SIAM Journal on Control and Optimization 1986;24:76–82. Ge R. The globally convexized #lled functions for global optimization. Applied Mathematics and Computation 1990;35:131–58. Bianchi N, Bolognani S. Design optimization of electric motors by genetic algorithms. IEE Proceedings—Electric Power Applications 1998;145:475–83. Faiz J, Finch J. Aspects of design optimization for switched reluctance motors. IEEE Transactions on Energy Conversion 1993;8:704–13. Liu X, Slemon G. An improved method of optimization for electrical machines. IEEE Transactions on Energy Conversion 1991;6:492–6. Liu X, Xu W. A novel optimization method for electrical machines. In: Proceedings of the 32nd Annual North America Power Symposium, Waterloo, Canada, 2000. p. 8–14 to 8–19. Slemon G, Liu X. Modeling and design optimization of permanent magnet motors. Electric Machines and Power Systems 1992;20:71–92. Courant R. Variational methods for the solution of problems of equilibrium and vibration. Bulletin of the American Mathematical Society 1943;49:1–23. Branin F, Hoo S. A method for #nding multiple extrema of a function of n variables. In: Lootsma FA, editor. Numerical methods of nonlinear optimization. London: Academic Press, 1972. p. 231–7.
80
X. Liu, W. Xu / Computers & Operations Research 31 (2004) 61 – 80
[25] Rastrigin L. Systems of extremal control. Moscow: Nauka, 1974 (in Russian). [26] Veinott C. Theory and design of small induction motors. New York: McGraw-Hill, 1959.
Xian Liu received a Master degree in operations research from the University of Waterloo, Canada, and a Ph.D. degree in computer engineering from the University of British Columbia, Canada. Now he is an assistant professor in the Department of Systems Engineering, University of Arkansas at Little Rock, USA. His primary research interests include operations research methods applied to computer engineering, global optimization, and communication networks. He has published 20 journal papers and several conference papers in above research areas. Wilsun Xu received a Master degree and a Ph.D. degree in Electrical Engineering from the University of Saskatchewan and the University of British Columbia, Canada. At present he is a professor in the Department of Electrical and Computer Engineering at the University of Alberta, Canada. His primary research interests include operations research methods applied to power systems and electrical machines. He has published more than 50 academic papers in above research areas.
Computers & Operations Research 31 (2004) 61 – 80
www.elsevier.com/locate/dsw
A new #lled function applied to global optimization Xian Liua;∗ , Wilsun Xub a
Department of Systems Engineering, University of Arkansas at Little Rock, Little Rock, AR 72204-1099, USA Department of Electrical and Computer Engineering, University of Alberta, Edmonton, Alta, Canada T6G 2G7
b
Received 1 March 2002; received in revised form 1 September 2002
Abstract The (lled function method (FFM) is an approach to #nd the global minimizer of multi-modal functions. The numerical applicability of conventional #lled functions is limited as they are de#ned on either exponential or logarithmic terms. This paper proposes a new #lled function that does not have such disadvantages. An algorithm is presented according to the theoretical analysis. A computer program is designed, implemented, and tested. Numerical experiments on typical testing functions show that the new approach is superior to the conventional one. The result of optimization design for an electrical machine is also reported.
Scope and purpose In the context of mathematical programming, global optimization is concerned with the theory and algorithms on minima of multi-modal functions. In general, global optimization approaches can be classi#ed into two categories: probabilistic and deterministic. The former can usually be applied to general multi-modal functions, whereas the latter typically concentrates on some particular classes of functions. The (lled function method is one of a few deterministic approaches which intend to #nd the global minimum for general multi-modal functions. However, the numerical performance of conventional #lled functions is undesirable as they are de#ned on either exponential or logarithmic terms or multiple parameters. This paper proposes a new #lled function that does not have the above disadvantages. The present work consists of theoretical analysis, algorithm design, computer implementation, mathematical validation, and engineering application. ? 2003 Elsevier Ltd. All rights reserved. Keywords: Filled function method; Global optimization; Mathematical programming; Minimization
∗
Corresponding author. Tel.: +1-780-492-3937; fax: +1-780-492-1811. E-mail addresses: [email protected] (X. Liu), [email protected] (W. Xu).
0305-0548/04/$ - see front matter ? 2003 Elsevier Ltd. All rights reserved. PII: S 0 3 0 5 - 0 5 4 8 ( 0 2 ) 0 0 1 5 4 - 5
62
X. Liu, W. Xu / Computers & Operations Research 31 (2004) 61 – 80
1. Introduction Researches in optimization methods applied to multi-modal functions have been active since the 1970s. Early contributions were reported in the two volumes named towards global optimization [1,2]. The progress since then can be found in [3,4]. In general, the reported approaches can be classi#ed into two categories: deterministic and probabilistic. It is generally agreed that the former concerns the covering method [5,6], the (lled function method [7,8], the trajectory method [9,10], and the tunneling method [11], whereas the latter refers to the clustering method [12,13] and other methods [14–16]. The well-known genetic and simulated annealing algorithms are also commonly regarded as the probabilistic approach. It seems that the (lled function method (FFM) has several advantages over others due to its relatively easy actualization with a process that aims at successively #nding smaller local minima. In the research of global optimization, there are two major issues: (1) How to #nd a lower minimizer of the objective function f(X ) from a known local minimizer of f(X ). (2) How to evaluate the convergence and, accordingly, design the stopping criteria. Since the global minimization concerns the whole region rather than the neighborhood of some particular points, there does not exist a general criterion, like the zero gradient criterion for the local minimizer, for the global minimizer. Although a great deal of eEorts has been made to #nd a general criterion to characterize the global minimizer, none is constructive. For this reason, so far none has been accepted to solve practical problems. On the other hand, many global minimization approaches introduced restrictions to the objective functions, such as global convexity, polynomial structure, rational features, etc. Without these restrictions, it is very diFcult, if not impossible, to conduct the convergence analysis. The situation is similar to the stopping criteria. Unlike other deterministic approaches, the FFM is a technique to #nd the global minimizer for general functions and does not assume any particular features of the objective function. In the literature, all researches on the FFM have mainly been concerning issue (1) above. In its nature, the FFM is a deterministic approach in the sense of issue (1). However, just like any other deterministic approaches for the global optimization, if no any restriction to f(X ) is enforced, the stopping criteria cannot avoid involving probabilistic strategies. Therefore, the research on stopping criteria, as well as the convergence analysis, has evolved to a specialistic area. This paper focuses on the #rst issue mentioned above. The rest of this paper is organized as follows. An overview of the FFM is given in Section 2. Then a class of new #lled functions, called the M -function, is de#ned and analyzed in Section 3. Next, in Section 4, an algorithm is presented and results of numerical experiments are reported. The practice for solving an engineering problem, optimization design of electrical machines, is addressed in Section 5. Conclusions are included in Section 6. Finally, the testing functions used for numerical experiments are described in Appendix A and the modeling issues of electrical machines are detailed in Appendix B. 2. Essentials of the FFM The FFM is an approach to #nd the global minimizer of a multi-modal function f(X ) on Rn , under the assumptions that f(X ) is continuously diEerentiable and f(X ) → +∞ as X → +∞.
X. Liu, W. Xu / Computers & Operations Research 31 (2004) 61 – 80
63
Thus, there exists a closed bounded domain that contains all global minimizers of f(X ). In the following, for the presentation purpose, will be called the operating region of the concerned global optimization problem. The FFM also assumes that f(X ) has only a #nite number of minimizers. To introduce the concept of the FFM, let us de#ne the following concepts: Denition 2.1. A basin of f(X ) at an isolated minimizer X1 is a connected domain B1 which contains X1 and in which starting from any point the steepest descent trajectory of f(X ) converges to X1 , but outside which the steepest descent trajectory of f(X ) does not converge to X1 . Denition 2.2. A hill of f(X ) at X1 is the basin of −f(X ) at its minimizer X1 , if X1 is a maximizer of f(X ). Denition 2.3. A local minimizer X2 is said to be higher than X1 if and only if f(X2 ) ¿ f(X1 ), and, for this case, B2 is said to be a higher basin than B1 . In this paper, Bh and Bl denote all higher and lower basins than current basin B1 of f(X ), respectively. Denition 2.4. A function P(X ) is called a #lled function of f(X ) at X1 if (1) X1 is a maximizer of P(X ) and the whole basin B1 becomes a part of a hill of P(X ); (2) P(X ) has no stationary points in any Bh s; and (3) There is a point X in a Bl (if such a basin exists) that minimizes P(X ) on the line through X and X1 . In this paper, we will allow an in#nite maximizer of P(X ). The FFM consists of two phases, local minimization and #lling: Phase 1: In this phase, a local minimizer X1 of f(X ) is found. Any eEective down-hill technique, for instance, the variable metric method, can be employed in phase one. Phase 2: In this phase, a #lled function is constructed and is minimized by the down-hill procedure. Phase 2 ends when such an Xs is found that Xs is in a Bl . Then, the FFM reenters phase 1, with Xs as the starting point, to #nd a new local minimizer X2 of f(X ) (if such one exists), and so on. The above process is repeated until the global minimizer is found. Several #lled functions have been proposed in the literature. Three popular ones are [7,8]: P(X; r; ) = exp(−X − X1 2 =2 )=[r + f(X )]; G(X; r; ) = −{2 ln[r + f(X )] + X − X1 p }; Q(X; a) = −[f(X ) − f(X1 )] exp(aX − X1 p );
(1)
where p = 1 or 2, r and are adjustable parameters, and a is an adjustable positive weight factor. Both P- and G-functions require two adjustable parameters, which need to be appropriately iterated and coordinated each other, hence their algorithmic realization is fairly complicated. For this reason, it is usually agreed that the Q-function is somehow better than the other two, since it involves only one adjustable parameter. However, the Q-function includes an exponential term whose argument
64
X. Liu, W. Xu / Computers & Operations Research 31 (2004) 61 – 80
is the product of the weight factor a and the norm. As a becomes larger and larger, as required to preserve the #lling property, the rapid increasing value of the exponential term will result in failure of computation even if the size of the feasible region is moderate. In practice, this kind of ill-conditioning problem frequently occurs. To make the Q-function work, many additional cares must be incorporated into the algorithm. It is obvious that the exponential term in the Q-function has seriously limited its applicability to the practical global optimization problems, especially those raised from engineering. 3. The M -function The M -function is de#ned as 1 M (X; a) = − aX − X1 2 ; (2) (f − f1 )1=m where f − f1 stands for f(X ) − f(X1 ), m is a natural number greater than 1, and a is a positive real used as the weight factor. For simplicity M (X; a) will be denoted as M (X ) throughout the paper. The #lling properties of M (X ) are exhibited by the following theorems. Theorem 3.1. Given d ∈ Rn and f(X ) ¿ f(X1 ), if dT ∇f(X ) ¿ 0 and dT (X − X1 ) ¿ 0, or dT ∇ f(X ) ¿ 0 and dT (X − X1 ) ¿ 0, then d is a descent direction of M (X ) at point X . Proof. It follows from (2) that dT ∇f(X ) T T d ∇M (X ) = − + 2ad (X − X1 ) : m(f − f1 )1+1=m Therefore, the conditions given guarantee dT ∇M (X ) ¡ 0. Consider the following de#nition: dT ∇f(X ) al (X ) = − 2m(f − f1 )1+1=m dT (X − X1 ) then there is another result: Theorem 3.2. If f(X ) ¿ f(X1 ), dT ∇f(X ) ¡ 0, and dT (X − X1 ) ¿ 0, then the sign of dT ∇M (X ) is that of al (X ) − a. Note that it is possible to have al (X )−a ¿ 0, because al (X ) → +∞ as f(X ) ¿ f(X1 ) and f(X ) → f(X1 ). These results clearly characterize the #lling property of the M (X ). Particularly, in the ascent region of the current basin B1 or a higher basin than B1 , d is always a descent direction of M (X ). On the other hand, in the descent region of a higher basin than B1 , d is still a descent direction of M (X ) provided that the weight factor a is suFciently large. Furthermore, in a lower basin than B1 , d may become an ascent direction of M (X ) and this possibility does exist. Therefore, under the continuously diEerentiable assumption for f(X ), M (X ) must have a stationary point along d.
X. Liu, W. Xu / Computers & Operations Research 31 (2004) 61 – 80
65
3.0 2.5 M(X), a=0.004, m = 2 2.0 1.5 1.0 f(X), M(X)
f(X) 0.5 0 −0.5 −1.0 M(X), a=0.03, m = 2
−1.5 b1
−2.0 2
4
b2
X1 6
8
b3 10
b4 12
14
b5 16
b6 18
20
22
X
Fig. 1. The behavior of M (X ), diEerent a.
The behavior of M (X ) may be well interpreted graphically when X is a one-dimensional variable. In this case, there are only two possible search directions, the positive and negative directions of the X -axis. Consider an arbitrary multi-modal function f(X ) that has a local minimizer X = X1 (Fig. 1). When X ∈ (X1 ; b2 ), M (X ) descends along the positive X -axis. When X ∈ (b1 ; X1 ), M (X ) descends along the negative X -axis. Let ˜i be the unit vector that points to the positive direction of the X -axis. Then the above features can be rephrased as follows: M (X ) always descends along the direction (X − X1 )˜i if X = X1 . The situation is similar in interval (b3 ; b4 ). This property holds for any positive weight factor a, according to the theory of M (X ). However, when X ∈ (b2 ; b3 ), weight factor a must be suFciently large to make M (X ) descending along the positive X -axis. On the other hand, when X ∈ (b4 ; b5 ), M (X ) will have a stationary point and the second phase may stop at that point. The global minimizer b6 can be reached by resuming phase 1 with b5 as the starting point. The eEect of m is shown in Fig. 2. From its de#nition, the M -function appears more applicable to computational assignments, because it retains such advantages of conventional #lled functions that the exponential terms and the second parameter are not needed. Other desirable features include that the formulation is of the sum of two terms, rather than a product. Also, the M -function has only the linear dependency on the norm, although the norm itself may be quadratical, whereas in the P- and Q-functions, the norm is used as the argument of a exponential function. In addition, as described below,
66
X. Liu, W. Xu / Computers & Operations Research 31 (2004) 61 – 80 2.5 M(X), a=0.01, m=2 2
1.5 M(X), a=0.01, m=3
1 f(X) f(X), M(X)
0.5
0 −0.5 −1 −1.5 b1
−2 2
4
X1
b2 6
8
b3 10
b4 12
14
b6
b5 16
18
20
22
X
Fig. 2. The behavior of M (X ), diEerent m.
the lower bound of weight factor a of the M -function is usually smaller than the one of the Q-function. The weight factor a plays a crucial role in a #lled function. Theoretically, the value of a must be suFciently large to preserve a desirable #lling capability. Computationally, the value of a should be small to make the numerical procedures healthy, because there is a product of a and the norm in the formulation. Therefore, a #lled function is a robust one if the value of al (X ) is small given a particular X . In the following, we compare the M -function with the Q-function in terms of the lower bound of a. Consider the de#nition below aq (X ) = −
dT ∇f(X ) : 2(f − f1 )dT (X − X1 )
It is easy to show from (1) that, with the same conditions as given in Theorem 3.2, if a ¿ aq (X ), then d is a descent direction of Q(X ) at point X . Furthermore, we have the following theorem: Theorem 3.3. Under the same conditions as given in Theorem 3.2, if f(X ) ¿ f(X1 ) + m−m , then al (X ) ¡ aq (X ) at X . Proof. f(X ) ¿ f(X1 ) + m−m leads to m(f − f1 )1=m ¿ 1. This in turn implies al (X ) ¡ aq (X ).
X. Liu, W. Xu / Computers & Operations Research 31 (2004) 61 – 80
67
The above theorem is important as it presents the fact that, even with a small weight factor a, the M -function preserves the desired #lling property as long as a proper m is chosen. We give some remarks to end this section. One might argue that M -function has a weakness because of the root term, such that it has no de#nition in the region corresponding to f(X ) ¡ f(X1 ) if m is even. This feature, however, does not aEect its application, because in the numerical iteration process of minimization, we always check the function value of current iteration point #rst. If at some Xs we obtain f(Xs ) ¡ f(X1 ), then Xs is in a lower basin than B1 already. Phase 2 can ends right here and the FFM can reenter phase 2 by starting from Xs . Therefore, the root term in the formulation of M -function is never an obstacle in the arithmetic realization of the FFM. Moreover, an odd m may be chosen if the above problem really needs to be avoided. Yet another concern might be the situation at such an iterating point X that the value of f(X ) − f(X1 ) is positive but very small, resulting in a relatively large ∇M (X ). As a matter of fact, this feature is just what is needed for any #lled functions: in the neighborhood of X2 such that {X2 = X | f(X ) = f(X1 ) and X2 = X1 }, a steep “barrier” is necessary in order to stop phase 2, as revealed by Theorem 3.2. This feature is diEerent in nature than the steep behavior inherited in the Q-function due to the eEect of the sequentially increased exponential term.
4. An algorithm and the numerical experiment In optimization studies, modeling, algorithmization, and numerical realization are three relevant yet independent aspects. It has been well known that the numerical realization of an established formulation for nonlinear programming usually involves more diFculties. The presented FFM is not an exception and even brings an additional layer of complexity. In the modeling aspect of local optimization studies, the function to be minimized is usually implicitly assumed to have the globally convex property: f(X ) → +∞ as X → +∞. Although in theory it is not impossible to develop a globally convex #lled function, there remain several analytical diFculties in the modeling aspect. As a matter of fact, almost all #lled functions reported so far are globally concave, including the M (X ) proposed in the present work. As a result, numerically speaking, many well-developed local minimization techniques may not be eEective in phase 2 of the FFM. In general, it seems that the algorithms not relying on the gradient, Hessian matrix, or conjugated directions are more appropriate. In early studies on the FFM, several local minimization methods, which compare the numerical values of the function only, were investigated and customized to minimize the #lled function. As demonstrated from extensive numerical experiments for various applications, the method adapted in the present work, the linear search approach, has been giving the most expected results. In this section, we present an algorithm for the n-dimensional case. The involved nomenclature in the algorithm is listed as follows: a: fmin : fmin0 : : L:
weight factor in the formulation of the #lled function value of the objective function at the local minimizer initial value of fmin , used for the iteration purpose. It is advisable to select fmin0 as a large positive value operating region of the concerned global optimization problem (de#ned in Section 2) counter recording the number of linear search operations in phase 2 of the FFM
68
X. Liu, W. Xu / Computers & Operations Research 31 (2004) 61 – 80
Table 1 Numerical experiments M -function
Six-hump camel-back Branin G–P Rastrigin Shubert III Sine-square I Sine-square II Sine-square III
Lm : X0 : #:
Q-function
G-function
Kf
Kt
Kf
Kt
Kf
Kt
204 97 354 226 406 784 113 997
346 252 498 572 569 2387 1716 3564
204 97 280 118 NF NF NF NF
346 3640 1059 3589
204 97 261 117 444 1063 113 1562
346 252 404 619 685 2665 1716 3784
speci#ed upper bound of L. An estimation is Lm ≈ 2n dmin D Gmax , where dmin is the smallest search step, D is a metric of , and Gmax = supx∈ ∇f(X ). initial point to start phase 1 of the FFM (local minimization) multiplier to update the weight factor a (# ¿ 1, usually taking 2.0)
The algorithm is detailed as follows: Step 0: Specify X0 , a, and Lm . fmin0 → fmin . Step 1: 0 → M . Enter phase 1 of the FFM. Activate the minimization procedure to minimize the objective function f(X ), starting from X0 . Find a local minimizer X1 . If f(X1 ) 6 fmin , then f(X1 ) → fmin ; otherwise. #a → a. Step 2: Enter phase 2 of the FFM. Step 3: Construct a search direction s using a random number generating procedure. L + 1 → M . Starting from X1 , activate the linear search procedure to minimize the #lled function. Step 4: Continue the linear search. Arrive at point X . Step 5: 1. If X ∈ , then go to Step 6; 2. If f(X ) ¡ f(X1 ), then X → X0 ; go to Step 1. 3. If X is the minimizer of the #lled function, then X → X0 ; go to Step 1. 4. Otherwise, go to Step 4. Step 6: If L ¡ Lm , then go to Step 2; otherwise taking the smallest minimum as the global one. Stop. Several testing functions were proposed in literature to evaluate the numerical performance of a new approach for global optimization. To estimate its eEectiveness, the M -function has been tested, with m = 2, with a set of testing functions detailed in Appendix A. Two conventional #lled functions, the Q- and G-functions, have also been tested. Identical starting points are selected for all #lled functions. The results of numerical experiments are presented in Table 1, where Kf is the total number of evaluations for the objective function and the #lled function when the global minimizer is found; Kt the total number of evaluations for the objective function and the #lled function when the algorithm terminated; and NF the failed to #nd the global minimizer.
X. Liu, W. Xu / Computers & Operations Research 31 (2004) 61 – 80
69
The evaluation of an algorithm or a formulation of the #lled function may involve several layers of the concerned numerical procedures. It is believed, however, that the number of iterations should not be regarded as an appropriate index [17]. This is because, in the down-hill search process, the number of iterations is dependent on what particular methods are used for the lower layers, such as unconstrained minimization or linear search. Typically, there may be several function evaluations per iteration, and it is these evaluations that consume the major part of total execution time in practical problems. For this reason, we used the total number of evaluations of the objective function and the #lled function as an appropriate measure for the performance of diEerent #lled functions. Note that Kt depends not only on Lm , but also on several implementation issues such as the particular linear search algorithm. For example, the number of evaluations of the objective function per linear search is generally diEerent with the Fibonacci, gold ratio, parabolic, or spine line interpolation procedure. In the present work, the parabolic interpolation procedure is used. In testing, identical initial weight factor and multiplier were used for all #lled functions. It is noted that the results of numerical experiments shown in Table 1, in terms of Kf and Kt , imply that the M -function is superior to the Q-function, especially, for the complicated functions like the Shubert III function, or the high-dimensional functions like the Sine-square functions. On the other hand, we also included the G-function as a reference. Strictly speaking, it is inappropriate to compare the G-function with the Q-function or the M -function, because the former has two adjustable parameters whereas the latter have only one. This means that the G-function has an extra degree of freedom. Nevertheless, it might be interesting to see how the G-function behaves, since no testing results for the G- or Q-function were reported in the original paper [7]. To make the G-function have the same degree of freedom as the G- or Q-function does, we used r = f(X1 ) and 2 = 1=a, as implied in [7]. For this particular setting, the M -function is superior to the G-function. We also observed that the value of a had little inNuence to the M - or G-functions except it had been frequently updated when the Shubert III function or the Sine-square functions was tested, because of the large number of local minimums or high dimensionality of these testing functions. This was not the case, however, for the Q-function. Since the Q-function includes an exponential term, in testing a moderate value a often resulted in overNow when the current iterative point X was far away from X1 , while the M - or G-functions still worked well in the same situation. For example, in testing the Q-function failed when a = 4:4, while the M - or G-functions worked well with a = 35:2. Therefore, we conclude that the main disadvantage of the Q-function is its exponential term. Perhaps this is why no any testing was reported in the original paper [7]. In summary, according to the theoretical analysis and numerical testing, it is reasonable to expect that the M -function will be well applicable to the usually more complicated engineering optimization programs. This is actually the practice described in the next section. 5. Optimization design of an electrical machine In the last decade there was a growing interest in design and manufacturing optimization of electrical machines due to increased cost of materials plus the cost of losses on their assumed lifetime [18–22]. A motor’s design can be described by a vector X of n components stating physical characteristics, such as dimensions of core and windings, as well as electromagnetic characteristics, such as densities of current and Nux. The design is subject to a set of m constraints that
70
X. Liu, W. Xu / Computers & Operations Research 31 (2004) 61 – 80
may include speci#cations arising from electrical, mechanical, thermal, manufacturing or standards limits. The task of the design optimization is to make a speci#ed objective function f(X ) reach its optimal value while keeping other technical indices within acceptable ranges. The complexity of electromagnetic characteristics incurred in motor design almost exclusively leads to multi-modal objective functions that have many local minimizers. The well-performed procedures for #nding a local minimizer usually lose their eFciency when applied to #nd the global minimizer of multi-modal functions. Therefore, it is really desirable to employ an approach based on the global optimization methodology to solve the concerned design problems. As examples to demonstrate the applicability of the proposed M -function to practical problems, optimized designs are produced for a single-phase induction motor. The original version of the FFM was developed for solving unconstrained problems under the assumption that operating region is known, where is a closed bounded domain that contains all global minimizers of f(X ). Since the design of electrical machines is a constrained optimization problem, our #rst challenge is to customize the FFM to the constrained environment. This process consists of two issues: con#guration and transformation. The former involves how to establish , whereas the latter deals with the constraints associated with the established . We discuss the con#guration issue #rst. By de#nition, there is a feasible region D in any constrained problem. Any solutions outside D are not useful. Therefore, for a motor design problem, we just look for its global optimal point inside D. This implies that the FFM only needs to be active throughout D. A natural approach seems to let coincide with D. This approach, however, has pitfalls. Recall that the FFM consists of two phases. In phase 1, a local minimizer X1 is found. At the beginning of phase 2, the M -function is constructed and it has an impulse at X1 . Then a down-hill search process is activated in the neighborhood of X1 . Accordingly, the position of X1 will directly aEect the behavior of phase 2. For motor design models that take cost as objective function subject to electromagnetic performance limits or vice versa, the local optimal design is usually located on or near boundaries of feasible region D. If let = D, then the impulse is just on or near one of the boundaries of . Consequently, starting near X1 the down-hill iterations intend to leave very soon and terminate one cycle of phase 2 without any meaningful results. For this reason, An applicable con#guration should have D ⊂ , and this depends on what particular technique is employed in phase 1 to integrate the constraints. The transformation issue cares about handling the constraints in optimization. One of the eEective approaches for solving constrained problems is the indirect method (IDM), which uses the augmented objective function rather than the original one. A typical augmented objective function consists of the weighted sum of all constraint terms and the original objective function. Here, the constraint term is de#ned to be a transformation of the constraint of interest. For example, one of the most popular IDMs, the exterior penalty function method (EPFM), takes the following form [23]: m F(X; r) = f(X ) + r {min[0; gi (X )]}2 ; i=1
where r is the penalty factor with an adjustable positive value. In the case of optimal points located on the boundaries of D, some IDMs may introduce a deviation. For instance, the minimizer found with the (EPFM) are outside D when r is reasonably large. As a result, a new region DF can be formed by enlarging region D with a small increment %.
X. Liu, W. Xu / Computers & Operations Research 31 (2004) 61 – 80
71
With an appropriate % value, all of the minima of F(X; r) are contained in DF . Therefore, we may make = DF . Note that D ⊂ has been satis#ed. Also, a small % is unlikely to lead to too much unnecessary search work. Based on the forgoing analysis, it is possible to use the convergence error &p to determine %. In the EPFM case, one of the most popularly adopted convergence criteria is the closeness measure of the minimizer to the boundaries. For instance, phase 1 can be terminated if gi (X ) ¿ − &p , where i = 1; 2; : : : ; m, and &p ¿ 0. Thus it is reasonable to choose % ¿ &p . The present work has integrated the EPFM into the FFM. For the single-phase induction motor of interest, the involved nomenclature is listed as follows: Lc C Dm Da Zm Za W H F X0 X∗ Cum Cua Al Si cm ca csi cal
length of core (cm), capacity of the operating capacitor (F), diameter of conductor in main winding (mm), diameter of conductor in auxiliary winding (mm), turns of main winding, turns of auxiliary winding, width of rotor end-ring (cm), height of rotor end-ring (cm), value of the objective function (cost of eEective materials), value of original design variables, value of optimal design variables, amount of copper of the main winding, amount of copper of the auxiliary winding, amount of aluminium of the rotor winding, amount of silicon sheet of the core, unit cost of copper of the main winding, unit cost of copper of the auxiliary winding, unit cost of silicon sheet, unit cost of aluminium,
The objective function is chosen as the cost of eEective materials and formulated as follows: f(X ) = Cum cm + Cua ca + Sicsi + Alcal : The details of electromagnetic modeling of concerned motors are included in Appendix B, since this section focuses on how to customize a theoretical optimization method and integrate it into a well-developed electromagnetic model. In the concerned design, eight variables are chosen to adjust the following parameters: length of core, operating capacitor, diameters of conductor in main and auxiliary windings, turns of main and windings, and dimensions of the rotor end-ring. The speci#ed constraints included the output power, eFciency (%), ratio of backward rotating #eld to forward rotating #eld (%), and dimensions of the rotor end-ring. Two sets of optimized results, D1 and D2, are shown in Tables 2 and 3, respectively. The original designs themselves had been produced by a team of experienced engineers using a computer program of motor analysis. Many diEerent combinations of design variables had been carefully selected to obtain a good solution. Therefore, the original designs of this motor were sub-optimal solutions
72
X. Liu, W. Xu / Computers & Operations Research 31 (2004) 61 – 80
Table 2 Design D1
X0 X∗
Lc
C
Dm
Da
Zm
Za
W
H
1:80 2:99
2:00 2:20
0:150 0:135
0:140 0:131
9647 4956
6711 4153
0:353 0:233
0:675 0:554
Lc
C
Dm
Da
Zm
Za
W
H
1:80 3:06
2:00 2:36
0:160 0:131
0:150 0:132
9647 5009
6711 4056
0:300 0:246
0:600 0:548
Table 3 Design D2
X0 X∗
Table 4 Minimizers of design D1 Index of minimizers Object function F
1 1:0000
2 0:9594
3 0:9490
4 0:9416
5 0:9240
6 0:9259
Index of minimizers Object function F
7 0:9218
8 0:9174
9 0:8963
10 0:8733
11 0:8708
12 0:8508
Index of minimizers Object function F
13 0:8668
14 0:8270
15 0:8140
16 0:8039
17 0:7988
Index of minimizers Object function F
1 1:0000
2 0:9885
3 0:9495
4 0:9270
5 0:9022
6 0:9079
Index of minimizers Object function F
7 0:8893
8 0:8761
9 0:8570
10 0:8399
11 0:8272
12 0:8166
Table 5 Minimizers of design D2
already. As an eEort to solve a practical problem, the M -function is applied to improve the original designs. The results are listed in Tables 4 and 5, where the objective function is normalized on the value of the original designs. For design D1, up to 17 local minimizers were found. while for design D2, the total number of identi#ed local minimizers is 12. The merit is obvious: The global optimization program made about 20% reduction in the cost of both designs while satisfying all speci#ed constraints. Thus the M -function indeed exhibited the desirable power. This is a signi#cant result since, especially, the starting points correspond to previously screened designs.
X. Liu, W. Xu / Computers & Operations Research 31 (2004) 61 – 80
73
6. Conclusions The #lled function method (FFM) is an approach to #nd the global minimizer of multi-modal and multidimensional functions. Conventional functions, however, rely on either an exponential term or multiple parameters. In numerical applications, these characteristics are strongly undesirable. In this paper, a new class of #lled functions called the M -function is proposed. The M -function has several desirable features such as (1) no any exponential or logarithmic terms; (2) only one parameter; and (3) formulated as a sum of two terms, rather than a product. The eEectiveness of the M-function applied to practical problems is exhibited by results of numerical experiments on typical testing functions. A computer program is designed based on the new approach and has been successfully applied to the optimization design of electrical machines. Appendix A. Testing functions Six-hump camel-back (n = 2) [9]: fC (X ) = 4x12 − 2:1x14 + x16 =3 + x1 x2 − 4x22 + 4x24
(−3 6 x1 ; x2 6 3):
The global minimizers are (0:08983; −0:7126) and (−0:08983; 0:7126). Branin (n = 2) [24]: fB (X ) = (x2 − 1:275x12 =,2 + 5x1 =, − 6)2 + 10(1 − 0:125=,) cos(x1 ) + 10 (−5 6 x1 6 10; 0 6 x2 6 15): The global minimizers are (−3:142; 12:275), (3:142; 2:275), and (9:425; 2:425). Goldstein–Price (G–P) (n = 2) [2]: fG (x) = [1 + (x1 + x2 + 1)2 (19 − 14x1 + 3x12 − 14x2 + 6x1 x2 + 3x22 ] ×[30 + (2x1 − 3x2 )2 (18 − 32x1 + 12x12 + 48x2 − 36x1 x2 + 27x22 )]
(−3 6 x1 ; x2 6 3):
The global minimizer is (0; −1). Rastrigin (n = 2) [25]: fR (X ) = x12 + x22 − cos(18x1 ) − cos(18x2 )
(−1 6 x1 ; x2 6 1):
This function has about 50 minimizers. The global minimizer is (0, 0). Shubert III (n = 2) [8]: 5 5 i cos[(i + 1)x1 + 1] i cos[(i + 1)x2 + 1] fS (X ) = i=1
i=1 2
+ [(x1 + 1:42513) + (x2 + 0:80032)2
(−10 6 x1 ; x2 6 10):
This function has 760 minimizers. The global minimizer is (−1:42513; −0:80032). Because of the large number of local minimizers and the steep slope around the global minimizer, the Shubert
74
X. Liu, W. Xu / Computers & Operations Research 31 (2004) 61 – 80
300 250 200 150 100 50 0 −50 −100 −150 10 10
5 5
0 x2
0
−5
−5 −10
x1
−10
Fig. 3. The Shubert function III.
function III has widely been recognized as an important testing function. An illustration is given in Fig. 3. Sine-square I (n = 6) [17]: n− 1 2 2 2 2 fsq1 (X ) = 10 sin (,x1 ) + (x n − 1) + (xi − 1) [1 + 10 sin (,xi+1 )] ,=n i=1
(−10 6 xi 6 10): This function has about 60 minimizers. The global minimizer is (1, 1, 1, 1, 1, 1). Sine-square II (n = 6) [17]: fsq2 (X ) =
2
2
10 sin (,y1 ) + (yn − 1) +
n− 1 i=1
yi = 1 + (xi − 1)=4
(−10 6 xi 6 10):
This function has about 30 minimizers. The global minimizer is (1, 1, 1, 1, 1, 1).
2
2
(yi − 1) [1 + 10 sin (,yi+1 )] ,=n;
X. Liu, W. Xu / Computers & Operations Research 31 (2004) 61 – 80
75
Sine-square III (n = 6) [17]: fsq3 (X ) = sin2 (3,x1 ) + (x n − 1)2 [1 + sin2 (2,x n )]
+
n− 1
2
2
(xi − 1) [1 + sin (3,xi+1 )]
10
(−10 6 xi 6 10):
i=1
This function has about 180 minimizers. The global minimizer is (1, 1, 1, 1, 1, 1).
Appendix B. An electromagnetic model of single-phase induction (SPI) motors Most SPI motors use an auxiliary winding with a capacitor to produce the start torque, since the magnetic #eld generated by a single-phase source remains stationary in space and pulses with time. In the following, the main winding and auxiliary winding are denoted as the M and A winding, respectively. One of the representative classes of SPI motors is the capacitor-run type, in which the A winding remains being connected to the source after starting. The present work is devoted to capacitor-run SPI motors. SPI motors in nature are asymmetrical two-phase motors, because both the eEective turns and the current magnitude of the M winding are usually diEerent than those of the A winding, and the phase diEerence between the two currents is not equal 90 electrical degrees. Most electromagnetic models of SPI motors reported in the literature can be considered being variants of the one described in [26], where the equivalent circuits of SPI motors are derived from the symmetrical component method. The circuits corresponding to the forward and backward phase sequences of the M winding and the A winding are shown in Fig. 4(a) – (d), respectively. In the present work, an electromagnetic analysis procedure is developed based on these equivalent circuits. rM1 XM1
XmM
ZM+
rM1 XM1
X’M2 r’M2
X’M2 XmM
ZM-
(a)
rA1 XA1 ZA+
(b)
rA1 XM1
X’A2 XmA
(c)
r’M2 2-S
S
r’A2 S
X’A2 XmA
ZA-
r’A2 2-S
(d)
Fig. 4. The equivalent circuits of SPI motors.
76
X. Liu, W. Xu / Computers & Operations Research 31 (2004) 61 – 80
The principal variables and parameters used in the electromagnetic analysis procedure are listed as follows: cos . power factor Ist ratio of starting current to rated current rA 1 resistance of the A winding rA 2 rotor resistance referred to the A winding rM 1 resistance of the M winding rM rotor resistance referred to the M winding 2 xA 1 leakage reactance of the A winding xA 2 rotor leakage reactance referred to the A winding xM 1 leakage reactance of the M winding xM rotor leakage reactance referred to the M winding 2 xmA magnetizing reactance of the A winding xmM magnetizing reactance of the M winding xa xmM + xM 2 Ka xmM =xa n rotating speed (rpm) pFe1 stator core losses pFe2 rotor core losses pmech friction and windage losses ps stray losses P1 input power P2 output power U voltage of the single phase power supply S slip Sf A slot space #lling factor of the A winding Sf M slot space #lling factor of the M winding Tm ratio of maximum torque to rated torque Tst ratio of torque to rated torque 4 eFciency B.1. Formulas for objective function and constraints The objective function and constraints in the design procedure are expressed in terms of a set of analytical formulas. According to Fig. 4(a), the forward sequence impedance of the M winding is ZM+ = rM1 + jxM1 + 2Zf : The 2Zf is the load impedance: 2Zf =
+ r =S) jxmM (jxM M2 2 ); rM2 =S + j(xmM + xM 2
= 2(Rf + jXf );
X. Liu, W. Xu / Computers & Operations Research 31 (2004) 61 – 80
77
where Rf =
=S 0:5Ka2 rM 2 ; =(Sx )]2 1 + [rM a 2
Xf =
=S)2 + x 0:5xmM (rM mM xa xM2 2 : (rM2 =S)2 + xa2
Similarly, from Fig. 4(b), we have ZM− = rM1 + jxM1 + 2Zb : The 2Zb is the load impedance: 2Zb =
+ r =(2 − S)] jxmM [jxM 2 M2 =(2 − S) + j(x rM2 mM + xM2 )
= 2(Rb + jXb ); where Rb =
=(2 − S) 0:5Ka2 rM 2 ; 1 + {rM2 =[(2 − S)xa ]}2
Xb =
=(2 − S)]2 + x 0:5xmM [rM mM xa xM2 2 : =(2 − S)]2 + x 2 [rM a 2
For the convenience of programming, in the design procedure a set of intermediary parameters is introduced to derive the armature currents, power factor, and the input power: RTM = rM1 + Rf + Rb ; XTM = xM1 + Xf + Xb ; RTA = rA1 + a2 (Rf + Rb ) + rc ; XTA = xA1 + a2 (Xf + Xb ) − xc ; Rfb = Rf − Rb ; Xfb = Xf − Xb ; 2 ); R3 = RTM RTA − XTM XTA − a2 (R2fb − Xfb
X3 = RTM XTA + XTM RTA − 2a2 Rfb Xfb ; R4 = R3 U=(R23 + X32 );
78
X. Liu, W. Xu / Computers & Operations Research 31 (2004) 61 – 80
X4 = R4 X3 =R3 ; R5 = R4 (RTA − aXfb ) + X4 (XTA + aRfb ); X5 = R4 (XTA + aRfb ) − X4 (RTA − aXfb ); R6 = R4 (RTM + aXfb ) + X4 (XTM − aRfb ); X6 = R4 (XTM − aRfb ) − X4 (RTM + aXfb ): The M winding current, the A winding current, and the total input current are: IM = R25 + X52 ; IA = R26 + X62 ; I1 =
(R5 + R6 )2 + (X5 + X6 )2 ;
respectively. The power factor is cos . = (R5 + R6 )=I1 : The input power is P1 = U (R5 + R6 ) + pFe1 : The air-gap power (transformed from stator to rotor) is 2 Pg = (IM + a2 IA2 )Rfb + 2a(Rf + Rb )(R5 X6 − R6 X5 ):
Finally, the output power is determined as follows: P2 = Pg (1 − S) − pmech − pFe2 − ps : B.2. The design procedure Given the speci#cations (ratings and performance constraints) and the core dimensions, a typical design procedure for SPI motors consists of the following steps: 1. 2. 3. 4. 5. 6. 7. 8.
derive the dependent core data; design the main winding; determine the resistance and reactances of the main winding; estimate the magnetomotive forces in multiple core regions; estimate the iron, mechanical, and stray losses; design the auxiliary winding; calculate the full-load performance; estimate the starting performance;
X. Liu, W. Xu / Computers & Operations Research 31 (2004) 61 – 80
79
9. determine the maximum torque; and 10. evaluate the cost of eEective materials. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24]
Dixon L, SzegRo G, editors. Towards global optimization. Amsterdam: North-Holland, 1975. Dixon L, SzegRo G, editors. Towards global optimization, vol. 2. Amsterdam: North-Holland, 1978. Horst R, Tuy H. Global optimization (deterministic approaches), 3rd ed. Berlin: Springer, 1996. S TRorn A, Zilinskas A. Global optimization. Berlin: Springer, 1989. Basso P. Iterative methods for the localization of the global maximum. SIAM Journal on Numerical Analysis 1982;19:781–92. Mladineo R. An algorithm for #nding the global maximum of a multimodal, multivariate function. Mathematical Programming 1986;34:188–200. Ge R, Qin Y. A class of #lled functions for #nding global minimizers of a function of several variables. Journal of Optimization Theory and Applications 1987;54:241–52. Ge R. A #lled function method for #nding a global minimizer of a function of several variables. Mathematical Programming 1990;46:191–204. Branin F. Widely convergent methods for #nding multiple solutions of simultaneous nonlinear equations. IBM Journal of Research Developments 1972;504 –22. Snyman J, Fatti L. A multi-start global minimization algorithm with dynamic search trajectories. Journal of Optimization Theory and Applications 1987;54:121–41. Levy A, Montalvo A. The tunneling algorithm for the global minimization of functions. SIAM Journal on Scienti#c and Statistical Computing 1985;6:15–29. Becker R, Lago G. A global optimization algorithm. In: Proceedings of the Eighth Allerton Conference on Circuits and Systems Theory, Monticello, IL, USA, 1970. p. 3–12. TRorn A. Cluster analysis using seed points and density determined hyperspheres as an aid to global optimization. IEEE Transactions on Systems, Man, and Cybernetics 1977;7:610–6. Archetti F. A probabilistic algorithm for global optimization problems with a dimensionality reduction technique. In: Balakrishnon A, Toma M, editors. Lecture notes in control and information sciences, vol. 23. Berlin: Springer, 1980. p. 36 – 42. Baba N. Global optimization of functions by the random optimization method. International Journal of Control 1979;30:1061–5. Dorea C. Limiting distribution for random optimization methods. SIAM Journal on Control and Optimization 1986;24:76–82. Ge R. The globally convexized #lled functions for global optimization. Applied Mathematics and Computation 1990;35:131–58. Bianchi N, Bolognani S. Design optimization of electric motors by genetic algorithms. IEE Proceedings—Electric Power Applications 1998;145:475–83. Faiz J, Finch J. Aspects of design optimization for switched reluctance motors. IEEE Transactions on Energy Conversion 1993;8:704–13. Liu X, Slemon G. An improved method of optimization for electrical machines. IEEE Transactions on Energy Conversion 1991;6:492–6. Liu X, Xu W. A novel optimization method for electrical machines. In: Proceedings of the 32nd Annual North America Power Symposium, Waterloo, Canada, 2000. p. 8–14 to 8–19. Slemon G, Liu X. Modeling and design optimization of permanent magnet motors. Electric Machines and Power Systems 1992;20:71–92. Courant R. Variational methods for the solution of problems of equilibrium and vibration. Bulletin of the American Mathematical Society 1943;49:1–23. Branin F, Hoo S. A method for #nding multiple extrema of a function of n variables. In: Lootsma FA, editor. Numerical methods of nonlinear optimization. London: Academic Press, 1972. p. 231–7.
80
X. Liu, W. Xu / Computers & Operations Research 31 (2004) 61 – 80
[25] Rastrigin L. Systems of extremal control. Moscow: Nauka, 1974 (in Russian). [26] Veinott C. Theory and design of small induction motors. New York: McGraw-Hill, 1959.
Xian Liu received a Master degree in operations research from the University of Waterloo, Canada, and a Ph.D. degree in computer engineering from the University of British Columbia, Canada. Now he is an assistant professor in the Department of Systems Engineering, University of Arkansas at Little Rock, USA. His primary research interests include operations research methods applied to computer engineering, global optimization, and communication networks. He has published 20 journal papers and several conference papers in above research areas. Wilsun Xu received a Master degree and a Ph.D. degree in Electrical Engineering from the University of Saskatchewan and the University of British Columbia, Canada. At present he is a professor in the Department of Electrical and Computer Engineering at the University of Alberta, Canada. His primary research interests include operations research methods applied to power systems and electrical machines. He has published more than 50 academic papers in above research areas.