A revised particle swarm optimization based discrete Lagrange multipliers method for nonlinear programming problems

Computers & Operations Research 38 (2011) 1164–1174

Contents lists available at ScienceDirect

Computers & Operations Research journal homepage: www.elsevier.com/locate/caor

A revised particle swarm optimization based discrete Lagrange multipliers method for nonlinear programming problems Ali Mohammad Nezhad , Hashem Mahlooji Department of Industrial Engineering, Sharif University of Technology, Tehran, Iran

a r t i c l e i n f o

a b s t r a c t

Available online 18 November 2010

In this paper, a new algorithm for solving constrained nonlinear programming problems is presented. The basis of our proposed algorithm is none other than the necessary and sufficient conditions that one deals within a discrete constrained local optimum in the context of the discrete Lagrange multipliers theory. We adopt a revised particle swarm optimization algorithm and extend it toward solving nonlinear programming problems with continuous decision variables. To measure the merits of our algorithm, we provide numerical experiments for several renowned benchmark problems and compare the outcome against the best results reported in the literature. The empirical assessments demonstrate that our algorithm is efficient and robust. & 2010 Elsevier Ltd. All rights reserved.

Keywords: Particle swarm optimization Discrete Lagrange multipliers Nonlinear programming Priority based feasibility strategy

1. Introduction One of the branches in operations research that due to its wide applications in the fields of management science, engineering optimization, and economics has attracted substantial attention over the years is nonlinear programming (NLP). In the realm of NLP, a general constrained nonlinear programming problem is formulated as follows: Minimize

f ðxÞ

s:t hi ðxÞ ¼ 0

i ¼ 1, . . . ,m1

gj ðxÞ r 0

j ¼ 1,:::,m2

ð1Þ

x A X DRn where x ¼ ½x1 ,x2 , . . . ,xn T is a vector of n decision variables, f(x) represents the objective function, hi(x) for i¼ 1,y,m1 show the equality constraints, gj ðxÞ for j¼1,y,m2 stand for the inequality constraints and X is a closed continuous set. Either of the functions f ðxÞ, hi ðxÞ or gj ðxÞ can be linear or nonlinear. One can cite numerous methods of solving nonlinear programming problems from the rich literature in this area. To name a few, feasible directional methods [1], gradient projection method [2], reduced gradient method [3], and penalty function methods like augmented Lagrangian penalty function [4] can be mentioned. In addition, within the past decade, a number of novel methods have also been proposed like discrete Lagrange multipliers method for NLP [5,6].  Corresponding author. Tel.: +98 21 6616 5717; fax: + 98 21 6602 2702.

E-mail addresses: [email protected], [email protected] (A. Mohammad Nezhad), [email protected] (H. Mahlooji). 0305-0548/$ - see front matter & 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.cor.2010.11.007

As a category of stochastic optimization methods, evolutionary computation techniques have attracted substantial interest on the part of scholars and practitioners in developing more efficient and robust computational procedures for solving complex optimization problems. Genetic algorithms are among the most well-known evolutionary algorithms. As a probabilistic search method, a genetic algorithm is designed to behave in accordance to the mechanism of natural evolution and selection. A genetic algorithm possesses the capability of simultaneously evaluating numerous points in the search space, so it stands a better chance of obtaining the global solution for the type of problems discussed here. Genetic algorithms enjoy the added advantage of ease of implementation, since they only use a scalar measure and not the derivative information in the search process. Another evolutionary algorithm known as particle swarm optimization (PSO) was introduced by Eberhart and Kennedy in 1995 [7], which is based on individual entities called particles. Imitating the social behavior of animals has been the key idea in developing the PSO algorithm. The PSO solution process requires particles to pass (fly) through the search space as their attributes such as positions and velocities are being updated on a constant basis with regard to the best performance of all particles as well as the best performance of each particle and its neighbors. PSO has some interesting properties to offer. For one thing, PSO has memory in the sense that particles retain good solutions (this feature is similar to the one involved in the elitist genetic algorithm). In addition, some important topologies of PSO enjoy the property of sharing information between all the particles (for instance, the global star topology). General quadratic penalty approaches are the commonly used constraint handling methods that the evolutionary algorithms employ. Static penalty methods consider substantial penalties and convert

A. Mohammad Nezhad, H. Mahlooji / Computers & Operations Research 38 (2011) 1164–1174

(1) into an unconstrained optimization problem [4] while dynamic penalty methods, on the other hand, try to increase the penalties in a gradual fashion in order to guard against the problems associated with the static penalty methods. These methods convert (1) into a sequence of unconstrained sub-problems, which tend to converge when all these sub-problems are solved optimally [4]. This means that, even if a single unconstrained sub-problem is not solved optimally, there is no assurance of reaching a constrained global minimum solution. Parsopoulos and Vrahatis [8] came up with a penalty framework for solving constrained optimization problems which employs a particle swarm paradigm (PSP) in the solution of resulting unconstrained sub-problems. Having introduced the notion of co-evolution, Coello [9] developed a penalty based algorithm by using a co-evolution model to adapt the penalty factors into the objective function in a genetic algorithm framework. Akin to the work of Coello [9], He and Wang [10] developed a co-evolutionary particle swarm optimization in which penalty factors and decision variables are optimized in a definite sequence by PSO. In spite of their popularity, the quadratic penalty methods may encounter the difficulty of ill-conditioning, caused by the sensitive effect of penalty factors on the convergence behavior that can lead to a premature termination or slow rate of progress [4]. To remedy this drawback, Eberhard and Sedlaczek [11] developed a particle swarm featured by the augmented Lagrange penalty function which enjoys the benefit of smoothness. Their algorithm implements a dynamic penalty adaptation approach leading to a better quality of solutions with respect to the general quadratic penalty methods. Recently, some interesting constraint handling methods have been proposed for the particle swarm optimization and genetic algorithms. For example, He and Wang [12] presented a hybrid particle swarm optimization in the context of a feasibility based rule that takes advantage of not using the penalty function to handle the constraints. To refuse the dominant effect of local optimum solutions, a simulated annealing algorithm has been incorporated in this framework to mutate the best solution with the aim of a more effective exploration of the solution space. Dong et al. [13] proposed a particle swarm optimization for nonlinear programming problems for which the constraints are handled with a priority based ranking method. Defining the concept of constraint fitness, they dealt with a constraint handling method to quantify the degree of violation for each particle. Chootinan and Chen [14] recommended using the gradient repair method to guide infeasible solutions toward the feasible region. They applied a simple GA configured by the gradient repair method to drive the constraints into satisfaction. With the above points in mind, Zahara and Kao [15] employed the Nelder–Mead simplex method along with PSO towards solving constrained optimization problems. They made use of gradient repair method of Ref. [14] as well as the priority based ranking method of Ref. [13] to manage the constraints. Aside from PSO and GA, one can find some research efforts about the simulated annealing in the context of stochastic optimization methods. Wah and Wang [5,6] offered an improved simulated annealing including the discrete Lagrange multipliers method for discrete constrained optimization problems. Pedamallu and Ozdamar [16] proposed a simulated annealing method supported by a local search algorithm called feasible sequential quadratic programming to solve general constrained optimization problems. The authors examined both versions of penalty and non-penalty based approaches to handle the constraints. The discrete Lagrange multipliers theory [17,18] is a collection of interconnected mathematical facts associated with the discrete space similar to the continuous Lagrange multipliers theory for the continuous space. More importantly, the theorems and lemmas defined in this theory cover nonconvex and nondifferentiable optimization problems. In fact, a strong mathematical background

1165

supports the discrete Lagrange multipliers method which not only does it but does not rely on any restrictive assumption, it is more widely applicable to tackle real life problems. This is rarely true for the other procedures of constraint handling such as the quadratic and augmented Lagrange penalty functions [4]. In this paper, a revised particle swarm optimization is proposed to solve constrained nonlinear optimization problems which employ the discrete Lagrange multipliers method (RPSO-DLM) to satisfy the constraints. We utilize PSO in this algorithm because of the following facts: (1) PSO is easy to implement in the context of programming and parameter setting. (2) PSO originally has been developed to address nonlinear programming problems [7] and more recently, some of its variants have been successfully applied to constrained optimization problems [12,15,27]. (3) The population based algorithms like PSO allow us to provide each entity (particle) with a separate vector of Lagrange multipliers rather than using one vector of Lagrange multipliers for all the entities (particles). The rest of this paper is organized as follows. In Section 2, the continuous and discrete Lagrange multipliers theories are concisely reviewed. In Section 3, we describe particle swarm optimization method and two of its approaches. In Section 4, we present our contribution for solving constrained nonlinear optimization problems called PSO-DLM. In Section 5, we propose a type of priority strategy to promote PSO which can be helpful to eliminate some defects of our algorithm. Section 6 contains a discussion with the aim of shedding light on how to fine tune the essential parameters and finally, in Section 7, we provide several numerical examples to establish the merits of our proposed algorithm.

2. Lagrange multipliers method In this section, we present a discussion on the continuous [19] and discrete Lagrange multipliers theory [17,18] as much as needed for developing our proposed method. Thus, here, we mainly concentrate on definitions, lemmas and theorems. As the framework of the discrete Lagrange multipliers theory is akin to its continuous counterpart, it is necessary to remind some important parts of the continuous Lagrange multipliers theory in preparation of a rather detailed discussion on this subject. 2.1. Continuous Lagrange multipliers Consider the mathematical program (1) and assume that the inequality constraints are somehow converted to equality constraints by a one to one correspondence (this will be shown in the concluding paragraph of Section 2). Hence, the formulation (1) can be expressed as the equality-constrained minimization problem (2) Minimize

f ðxÞ

s:t hi ðxÞ ¼ 0

i ¼ 1, . . . ,m

ð2Þ

xA X where f ðxÞ and hi ðxÞ for i¼ 1,y,m stand for the objective function and equality constraints, respectively, x ¼ ½x1 ,x2 , . . . ,xn T is the vector of decision variables and m¼ m1 +m2. Let X DRn unless explicitly stated otherwise. The continuous Lagrange multipliers method was originally developed to address continuous constrained optimization problems of this type. This method, attempts to get a local minimum by transforming an equality constrained

1166

A. Mohammad Nezhad, H. Mahlooji / Computers & Operations Research 38 (2011) 1164–1174

problem into an unconstrained problem for which the objective function called the Lagrangian function, is expressed as follows: Lðx, kÞ ¼ f ðxÞ þ kT hðxÞ where k ¼ ½l1 , l2 , . . . , lm T shows the vector of Lagrange multipliers corresponding to the equality constraints and hðxÞ ¼ ½h1 ðxÞ, . . . , hm ðxÞT . Then, the associated problem is minimized with respect to x in order to acquire an extremum point of the Lagrangian function over x A X. To formally proclaim the basics of the continuous Lagrange multipliers theory, it is necessary to borrow the following theorems from Ref. [19]: Theorem 1. Suppose that f and hi for i¼1,y,m are real valued and continuously differentiable on X. Let x be a local minimum of (2) which (obviously) satisfies hðx Þ ¼ 0 and suppose that the Jacobin matrix of hðx Þ is full rank (x is designated as a regular local minimum). Accordingly, it turns out that the gradient vector of f at x can be represented as a linear combination of the gradient vectors of the equality constraints at x

rx f ðx Þ þ rx hðx ÞT k ¼ 0 where rx f ðx Þ denotes the gradient of f at x , rx hðx Þ stands for the Jacobin matrix of h at x* and k* is the corresponding vector of Lagrange multipliers.

where at shows the step size value and t denotes the index of iteration. Factually, Eq. (3) implements a movement toward minimizing the Lagrangian function with respect to x and Eq. (4) modifies the Lagrange multipliers with the aim of deriving the violated constraints to be satisfied until reaching a feasible solution after which this modifying step comes to a halt. Hence, the former performs minimization in the decision variables space while the latter embarks on maximization in the Lagrange multipliers space. It must be noted that once a feasible solution is found, as the second term of L(x, k) is dropped (since h(x)¼0), Eq. (3) resumes minimizing the real objective function of (2). Following this procedure for a finite number of iterations, yields a local minimum, if it satisfies the second order sufficient conditions of Ref. [19]. Based on our presentation of the first order method, one may expect to certainly get a continuous saddle point which is defined in Ref. [4]. By way of an example, one can show that this perception is not necessarily true (see [18]). On the other hand, it is proved that in the continuous space, the set of constrained local minima encompasses the set of solutions that satisfy the first order necessary (Theorem 2) and second order sufficient [19] conditions of local minimum solutions. That is to say, this undesirable property reflects the fact that the continuous Lagrange multipliers method may be unable to access some local minima and in this way, it easily gets trapped in an inferior local minimum. 2.2. Discrete Lagrange multipliers

Proof. One can refer to Ref. [19] for a formal proof. This theorem states that the gradient of the Lagrangian function with respect to x is expected to vanish at the regular local minimum ðx , k Þ. All this means that rx Lðx , k Þ ¼ 0 should be satisfied. We now present Theorem 2 which specifies the first order necessary conditions for a solution to be a local minimum [19]. Theorem 2. Suppose that the assumptions of Theorem 1 are satisfied. Then, a vector of Lagrange multipliers k can be found such that rx Lðx , k Þ ¼ 0. Proof. See Ref. [19] for a formal proof. As mentioned at the beginning of this section, the continuous Lagrange multipliers method would end up with an extremum solution for the converted unconstrained problem whereas this solution may be infeasible for (2) if the Lagrange multipliers are not appropriately set. In other words, aside from minimizing the Lagrangian function, the equality constraints must be consistent with the generated solution which equivalently confirms the assertion of Theorem 2. In the sequel, it can be deduced that for the regular local minimum x*, the succeeding system of equations should be satisfied ( rx Lðx , k Þ ¼ 0

rk Lðx , k Þ ¼ 0 in which the first set of equations is the same as the one in Theorem 2 and the second one indicates that h(x*) ¼0 should hold. Furthermore, as the name implies, the necessary conditions do not strictly guarantee that x* to really be a local minimum and for this purpose, the second order sufficient conditions of Ref. [19] (not relevant to the main theme of this paper) should be checked. In dealing with this system of equations, numerous solving methods have been developed such as the first order method [18], Newton and quasiNewton methods [4], and sequential quadratic programming method [4]. Due to its simplicity, we elaborate on the first order method which consists of two recursive equations given by xt þ 1 ¼ xt at rx Lðxt , kt Þ

ð3Þ

kt þ 1 ¼ kt þ at hðxt Þ

ð4Þ

The discrete Lagrange multipliers (DLM) theory establishes a strong mathematical foundation to address discrete constrained optimization [17]. This theory extends the concept of a continuous saddle point to the discrete space and through defining a Lagrange multipliers formulation, presents the first order necessary and sufficient conditions for a saddle point to be a local minimum. As a major goal, the DLM attempts to elaborate on a set of generic theorems and lemmas dealing with discrete constrained optimization which are not responsive to the problem structure. In addition, the DLM proves that in the discrete space, the set of constrained local minima and the set of saddle points are equivalent. To follow the discussion on the DLM, we now consider the mathematical program (2) for which X DD. Here D denotes an arbitrary set of discrete values. To present the DLM, we first define the generalized Lagrangian function in the discrete space as Ld ðx, kÞ ¼ f ðxÞ þ kT HðhðxÞÞ

ð5Þ

in which, H is an arbitrary continuous transformation function, k shows the vector of Lagrange multipliers corresponding to the equality constraints and hðxÞ ¼ ½h1 ðxÞ, . . . ,hm ðxÞT . In the same way as the continuous Lagrange multipliers theory, the DLM comes up with a Lagrangian formulation to relax the functional constraints of (2) but within the discrete space. However, unlike the continuous Lagrange multipliers, it is proved that the DLM and the concept of saddle point are intimately interconnected so that the necessary and sufficient conditions for a saddle point can be used as the necessary and sufficient conditions for a local minimum solution interchangeably. Now, a formal definition of discrete saddle point is in order. Definition 1. A discrete saddle point is a solution like ðx , k Þ for the generalized Lagrangian function such that Ld ðx , kÞ r Ld ðx , k Þ r Ld ðx, k Þ

x A Ne ðx Þ,

k A Rm

ð6Þ

where Ne ðx Þ denotes an e- neighborhood set around x and e is a small positive value. 



In tackling a discrete optimization problem, it is quite obvious that the gradient vector can no longer be defined. Instead, we adopt the concept of direction of maximum potential drop as a tool to,

A. Mohammad Nezhad, H. Mahlooji / Computers & Operations Research 38 (2011) 1164–1174

loosely speaking, equivalently specify the improving direction. The direction of maximum potential drop (DMPD) is a direction vector defined on Ld ðx, kÞ at point x toward a point likeysuch that Ld has a minimum value at y A Ne ðxÞ [ fxg. In other words, DMPD is defined as

Dx Ld ðx, kÞ ¼ ½y1 x1 , . . . ,yn xn T where Ld ðy, k ¼ ÞMin Ld ðxu , kÞ and xu A Ne ðxÞ [ fxg. In light of this definition, Lemma 1 looks into the first order necessary and sufficient conditions for a discrete saddle point (for a proof of this lemma, see Ref. [17]). Lemma 1. The solution ðx, kÞ in the discrete space is a saddle point if the following conditions are satisfied:

Dx Ld ðx, kÞ ¼ Dx ðf ðxÞ þ kT HðhðxÞÞÞ ¼ 0 hðxÞ ¼ 0

multipliers vector k0 exceed the ones in the optimal Lagrange multipliers vector. Corollary 1. Let ðx , k Þ be a discrete saddle point for the generalized Lagrangian function. The solution ðx , kuÞ can also be a discrete saddle point even if ku Z k holds. In other words, satisfaction of the left hand side of inequality (6) is sufficient for ðx , kuÞ to be a discrete saddle point. Proof. Refer to Ref. [17]. Now that we have reviewed the theoretical basis of the DLM, we elaborate on the way the inequality constraints can be transformed to equality constraints in formulation (2). To achieve this purpose, the maximum function is applied to convert the inequality constraints as gj ðxÞ r0 3max f0, gj ðxÞg ¼ 0

where D indicates the direction of maximum potential drop as defined above. In fact, the conditions that Lemma 1 sets out are similar to the ones in Theorem 2 with the exception of using DMPD instead of the gradient vector. Therefore, taking the idea of the first order method in continuous space, a sequential search procedure can be developed in terms of the conditions contained in Lemma 1 to locate a saddle point with moderately low cost of the generalized Lagrangian function given by xt þ 1 ¼ xt  Dx Ld ðxt , kt Þ

ð7Þ

kt þ 1 ¼ kt þ cHðhðxt ÞÞ

ð8Þ

where  is the addition operator in the discrete space and c is a positive step size value for increasing the Lagrange multipliers values corresponding to the violated constraints. Here, t stands for the iteration number. Eq. (8) looks into maximizing the generalized Lagrangian function with respect to the Lagrange multipliers to fulfill the violated constraints whereas Eq. (7) minimizes the generalized Lagrangian function to generate low valued objective function solutions. Roughly speaking, this procedure aims at finding a discrete saddle point. Note that the corresponding Lagrange multiplier of a constraint will no longer be changed by the time it becomes satisfied. The next lemma demonstrates the sufficient conditions whereby the set of saddle points exactly matches the set of constrained local minima in the discrete space (for a proof of this lemma, see Ref. [17]). Lemma 2. The set of all discrete saddle points is equivalent to the set of all constrained local minimum solutions provided that the non-negative or non-positive continuous transformation function H(as defined in Eq. (5)) fulfills the property given by HðxÞ ¼ 03x ¼ 0 Lemma 2 ascertains that to access a local minima, it suffices to restrict our attention to the set of discrete saddle points which are generated by the above first order method. With this in mind, the constrained global minimum solution can be obtained provided that a global optimization strategy is employed which is capable of seeking the most eligible saddle point. Since the absolute and square functions are the two common instances for this transformation function, in this work, we choose to employ the absolute function. At times, a difficulty may arise in the sense that when the optimum values of Lagrange multipliers are unknown, we may be unable to exactly specify them in agreement with the definition of a discrete saddle point. To overcome such difficulty, the subsequent corollary concludes that the solution ðx , kuÞ will be a discrete saddle point even when the elements of its corresponding Lagrange

1167

j ¼ 1, . . . ,m2

Therefore, the generalized Lagrangian function for (1) in which X D D, can be expressed as follows: Ld ðx, k, lÞ ¼ f ðxÞ þ

m1 X

li Hðhi ðxÞÞ þ

i¼1

m2 X

mj Hðmaxf0,gj ðxÞgÞ

j¼1

in which, l is the Lagrange multipliers vector corresponding to the original inequality constraints. Obviously, the max function is nondifferentiable, but the preceding theorems and lemmas can still be employed since none of them requires Ld ðx, k, lÞ to be differentiated.

3. Particle swarm optimization 3.1. Canonical approach Particle swarm optimization (PSO) algorithm proposed by Eberhart and Kennedy [7] is one of the most novel meta-heuristic algorithms that was essentially developed for solving nonlinear programming problems. Principally, the main notion of PSO has been derived from the social behavior of animals in the nature. PSO acts similar to the genetic algorithms in the sense of its population based approach and imitating the evolutionary process in the nature but unlike the genetic algorithms, no entity is discarded from the population. In PSO, the social interaction of a group in the nature is imitated in the sense that the group members are inclined to move toward the best member in the group. Hence, the behavior of each member is formed by the personal and social information. The general procedure of PSO can be described step by step as follows: Step 1: Initialize a swarm of members called particles in the solution space. Assign a random position and velocity to each particle. Set the current positions as the particles’ best positions (pi ð0Þ) and the position of the best evaluated particle as the global best position (pg ð0Þ). Here, i nominates the particle number. Set t’0 and go to Step 2. Step 2: Modify the velocities and positions according to Eqs. (9) and (10), respectively vi ðt þ 1Þ ¼ ovi ðtÞ þc1 r1 ðpi ðtÞxi ðtÞÞ þ c2 r2 ðpg ðtÞxi ðtÞÞ xi ðt þ 1Þ ¼ xi ðtÞ þvi ðt þ 1Þ

ð9Þ ð10Þ

where o is representative of the inertia weight, c1and c2 are the weighting coefficients and rj(j ¼1,2) are two uniform random numbers inside (0,1). The velocity of each particle is limited to a certain bound vmax that is usually set to a fraction of the solution

1168

A. Mohammad Nezhad, H. Mahlooji / Computers & Operations Research 38 (2011) 1164–1174

space. The inertia weight is determined based on this equation ðo o Þ o ¼ omax  max min t tmax in which, t and tmax show the iteration number and the maximum number of iterations respectively. Shi and Eberhart [20] state that PSO parameters will be fine tuned by using the following values: c1 ¼ c2 ¼ 2

omax ¼ 0:9 omin ¼ 0:4

Afterwards, go to Step 3. Step 3: Evaluate the current positions of the particles. Compare the current fitness with the best fitness for each particle. If a better position has been found, update the best position. Also update the global best position when needed and go to Step 4. Step 4: Check the termination criterion. If the condition is met, terminate the algorithm and deliver the best obtained solution. Otherwise, set t’t þ 1 and return to Step 2. 3.2. Constriction coefficient approach As the canonical approach of the PSO can be regarded as a kind of difference equation, it is possible to analyze the search procedure by Eigen values of this difference equation. Moreover, it is difficult to set an appropriate value for vmax due to its main effect on the convergence rate. Hence, to omit this obstacle and exploit the special characteristics of the search procedure, the constriction coefficient approach of PSO developed by Clerc and Kennedy [21,22] is applied. In this approach, Eq. (9) is replaced by Eq. (11) to update the velocities vi ðt þ1Þ ¼ wðvi ðtÞ þc1 r1 ðpi ðtÞxi ðtÞÞ þc2 r2 ðpg ðtÞxi ðtÞÞÞ

ð11Þ

where

w¼

2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi j2j j2 4jj

ð12Þ

in which j ¼ c1 +c2 and c1 ¼c2 ¼2.05. Clerc [21] noted that upon using the constriction coefficient approach, the algorithm does not come to a premature convergence before a complete search of the solution space. In this approach, the position of the swarm is modified using Eqs. (10), (11) and (12).

4. PSO based discrete Lagrange multipliers As proclaimed in Section 2, in contrast to the discrete space, a solution that satisfies the necessary and sufficient conditions in the continuous space may not be a suitable constrained local optimum solution [18]. In other words, the discrete Lagrange multipliers theory is stronger than the continuous case in terms of the optimality conditions. Hence, if there exists a way to transform a continuous optimization problem into a discrete one and apply the discrete Lagrange multipliers method, it may be possible to obtain the global optimum solution in an efficient manner. Wah and Wang [5] assert that as digital computers transform continuous data into discrete data, due to using an approximation up to a certain precision, no significant difference exists when the discrete Lagrange multipliers method is employed to solve a continuous optimization problem. Generally speaking, our aim is to implement the continuous particle swarm optimization algorithm based on discrete Lagrange multipliers method toward solving a continuous optimization problem. Toward this end, the first set of equations of first order method (Eq. (7)) is replaced by PSO so as to explore the space X and the vector of Lagrange multipliers is adjusted by means Eq. (8). In an attempt to gain a saddle point, as defined in Eq. (6), our algorithm consists of five steps for a minimization problem as discussed below. The first part of the algorithm (Step 3) carries out minimization of Ld ðx, kÞ with

respect to x for a fixed k and the second part of the algorithm (Step 4) maximizes Ld ðx, kÞ with respect to k for a fixed x. It seems that a reasonable combination of minimizing and maximizing the generalized Lagrangian function (Steps 3 and 4) requires a higher chance of resorting to Step 3 at each iteration as Wah and Wang [5] remark. Observing such a tendency, hints at using a non-uniform distribution (e.g. a right skewed beta distribution with parameters a ¼2 and b ¼3) in accordance with imitating this biased behavior. That is to say, a random deviate, generated from such a nonuniform distribution determines which of these two steps should be taken. In our proposed algorithm, we choose to work with beta (a ¼2, b ¼3). For illustration purposes, Fig. 1 shows a right skewed beta distribution in which the vertical dashed line stands for the threshold bound which is decided by the user and the area on the left side of the dashed line indicates the probability of selecting Step 3. To initialize the algorithm, we randomly generate a swarm of particles and assign an initial velocity to each particle. Also, for each particle, a Lagrange multipliers vector is generated within the predetermined bounds. By this representation, each particle can be shown as an ordered-pair like (x,k). Then, the particles are appraised through evaluating the generalized Lagrangian function (Step 1). Next, a random deviate from the beta distribution is generated (Step 2). If this value is less than the predetermined threshold, we go to Step 3, otherwise we go to Step 4. In Step 3, PSO is employed and one iteration is done to minimize the generalized Lagrangian function with respect to x. In Step 4, if any of the constraints corresponding to a particle is violated, its Lagrange multiplier is amplified. We note that if a constraint is satisfied at this step, we do not need to change its Lagrange multiplier. At the end of any of these two steps (3 or 4), the stopping condition is checked (Step 5). Now, we formally present our algorithm by way of a detailed stepwise fashion as follows: Step 1: Initialize the algorithm. Step 1.1: Initialize the main parameters. Let S be the swarm size, r the threshold value for the beta distribution, c the step size value for penalizing the violated constraints, n the number of variables, m the number of constraints, t the iteration number, FEa counter for the number of function evaluations, FE_Max the maximum number of function evaluations, pbest the best feasible solution, and finally, Fbest the generalized Lagrangian function value of pbest. Set Fbest ¼N, S0 ¼{1,y,S} and M¼{1,...,m}. Substitute the absolute function for the transformation function H. Step 1.2: Distribute the initial particles’ positions (xi ð0Þ) and velocities (vi ð0Þ). Moreover, corresponding to each particle, generate an m-dimensional vector of Lagrange multipliers (ki ð0Þ). Here i refers to the particle number. Set t’0. Step 1.3: Evaluate the particles by the generalized Lagrangian function. Determine the particle’s best solution (pi ðtÞ), the global best solution (pg ðtÞ), and their corresponding fitness values. Define

3 Threshold

2.5 2 1.5 1

Step 3

0.5

Step 4

0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

Fig. 1. Plot of a right skewed beta distribution.

0.9

1

A. Mohammad Nezhad, H. Mahlooji / Computers & Operations Research 38 (2011) 1164–1174

Fgt as the global best objective value or the fitness value for the global best solution and Fit as the particle’s best objective value or the fitness value for the particle’s best solution for i¼1,y,S. Step 1.4: Calculate the degree of violation for each particle which reflects the maximum deviation of equality constraints for the particle’s solution given by

1169

where | | is indicative of the absolute function. Now, go to Step 2. Step 2: Generate a random deviate R from the beta distribution. If Rr r, then go to Step 3, otherwise go to Step 4. Step 3: Modify the particles’ positions using Eqs. (10) and (11). Then, evaluate each particle and determine its degree of violation by Eq. (13). Accept the current particle’s solution as the particle’s best solution if its current position is preferred to its position in the previous iteration in the context of fitness and feasibility. Also, update the global best solution whenever needed. Mathematically speaking, for any particle i¼1,y,S, if Ld ðxi ðt þ 1Þ, ki ðtÞÞ r Fit and violit + 1 r violit hold, then set pi ðt þ 1Þ ¼ xi ðt þ 1Þ as the particle’s best solution and Fit þ 1 ¼ Ld ðxi ðt þ 1Þ, ki ðtÞÞ as its generalized Lagrangian function value. Otherwise, pi ðt þ1Þ ¼ pi ðtÞ andFit + 1 ¼Fit. Moreover, modify the global best solution and its objective as follows: let N ¼{i|violit + 1 rviolit, i¼1,y,S}. If minfFit þ 1 g r Fgt holds, then set

stop even before an exhaustive search is complete. In fact, a suboptimal feasible solution pg(t) can make the algorithm prematurely come to a halt, due to concentration of other particles in the suboptimal region. As an alternative, we develop a type of priority strategy governing the particles for the sake of a concurrent diversification and intensification once a feasible solution is found. The priority based feasibility strategy (PBFS) attracts the more feasible particles (within a predetermined bound) instead of the ones that attain better generalized Lagrangian function values, until a feasible solution is found. Also, in this strategy, infeasible particles are not denied, until the whole swarm of particles is contained in the feasible domain. Hence, in contrast to a typical penalty approach, this strategy does not discard any piece of information on infeasible particles. In the priority based feasibility strategy, the swarm is partitioned into two independent categories of feasible and infeasible particles. In the first category, the best particle is the one that attains the highest fitness value (lowest amount of the generalized Lagrangian function value) while in the second, the best is the one that attains the lowest degree of violation in the swarm. In order to resolve the noted problem in the PSO-DLM algorithm, Steps 3 and 4 of PSO-DLM are revised as follows: Let Fg0 ¼ minfFi0 g if viol0 r 10 e6, otherwise, Fg0 ¼N where

pg ðt þ 1Þ ¼ pI ðt þ1Þ and Fgt + 1 ¼FIt + 1 where I ¼ argminfFit þ 1 g stands

I ¼ argminfFi0 g.

for the index of the global best solution. Else, pg ðt þ 1Þ ¼ pg ðtÞ and Fgt + 1 ¼Fgt. Set ki ðt þ1Þ ¼ ki ðtÞ for i¼1,y,S and go to Step 5. Step 4: Adjust the vectors of Lagrange multipliers according to Eq. (8). Due to penalizing the infeasible particles, the Lagrange multipliers corresponding to the unsatisfied constraints are as much increased as to force them to be satisfied whereas the ones corresponding to the satisfied constraints will remain unchanged in the next iterations. Such adjustment with regard to replacing the Htransformation function by the absolute function is done as follows:

Step 3: The evaluations, substitutions and comparisons are done as before with the exception of updating the global best solution such that if minfFit þ 1 g rFgt , then set pg ðt þ 1Þ ¼ pI ðt þ 1Þ where

violti ¼ maxfjhj ðxi ðtÞÞjg

8i ¼ 1, . . . ,S

jAM

ð13Þ



iAN

iAN

ki ðt þ1Þ ¼ ki ðtÞ þ cwt

wtj ¼ maxfjhj ðxi ðtÞÞjg i A Su

j ¼ 1, . . . ,m

where wt ¼ ½wt1 ,wt2 , . . . ,wtm T . Due to making a change in the Lagrange multipliers, the current attributes such as Fgt + 1 and Fit + 1 should be adapted to the new generalized Lagrangian function for the next iterations by Fit þ 1 ¼ Ld ðpi ðtÞ, ki ðt þ 1ÞÞ

i ¼ 1, . . . ,S

Fgt þ 1 ¼ FIt þ 1 where I is the index of the global best solution defined earlier. Set xi ðt þ 1Þ ¼ xi ðtÞ, pi ðt þ1Þ ¼ pi ðtÞ and go to Step 5. Step 5: Check the stopping rule and deliver pbest as the final solution if FE 4FE_Max. Otherwise, if violIt + 1 r10e 6 and Fgt + 1 r Fbest, then Fbest ¼Fgt + 1 and pbest ¼ pg ðt þ 1Þ. Set t’t þ1 and return to Step 2.

5. Priority based feasibility strategy A number of research efforts have been reported on the issue of the swarm convergence to the global optimum solution [21,23]. In this relation, Van den Bergh [23] noticed a critical problem in the original PSO caused by the particles’ positions coinciding with the global best position. In this case, if the particles’ velocities are close to zero, the PSO algorithm may come to a premature convergence, a condition which is more crucial in our algorithm. To be more specific, if the number of equality constraints is increased, the size of feasible region tends to become too small to find a feasible solution. Thus, the structure of this algorithm may bring about a

i A Su

I

i A Su

iAN

I ¼ argminfFit þ 1 g. Otherwise, set pg ðt þ 1Þ ¼ pg ðtÞ. iAN

Step 4: In this case only the renewing rule of Fgt + 1 is revised: If violIt + 1 r 10 e  6, then Fgt + 1 ¼FIt + 1, otherwise Fgt + 1 ¼N. These modifications make the global best objective not to be updated in Step 3 and this issue is postponed until performing Step 4. In addition, the global best objective is given a substantially large value if it corresponds to an infeasible particle. Hence, during the initial iterations of the algorithm, almost every infeasible particle, succeeded in cutting down the degree of violation, can take the place of the global best solution and in this way temporarily attract the attention of other particles. Upon detecting a feasible particle in the position of the global best solution, the global best objective takes a finite value and then a greater number of particles are filtered because of being compared with a feasible particle. However, infeasible particles are not completely refused from consideration since there may exist an infeasible particle which has a lower value of the generalized Lagrangian function with respect to a feasible particle. If such a particle is found, two cases may occur: (1) The particle is feasible. In this case, as a better feasible particle has been found, the global best objective or Fgt + 1 take a tighter value and the domain of comparison becomes more limited. (2) The particle is infeasible. In this case Fgt + 1 ¼N provides another opportunity for infeasible particles to be considered as the global best solution. In fact, these two cases imply amplifying the PSO-DLM in spirit of exploitation and exploration to avoid getting trapped in local optimum solutions. Incorporating the feature of smartness in the particle swarm optimization by using this procedure can establish a full cooperation among all particles in order not to accept a suboptimal solution (to be tested in Section 7). With this modification, our algorithm is called RPSO-DLM (revised particle swarm optimization) hereafter.

1170

A. Mohammad Nezhad, H. Mahlooji / Computers & Operations Research 38 (2011) 1164–1174

6. Sensitivity analysis and parameter tuning

Average objective value

40

ρ = 0.6 50

ρ = 0.9 60

70

80

90

100

-30365 -30390 -30415 -30440 -30465 -30490 -30515 -30540 -30565 -30590 -30615 -30640 -30665

ρ = 0.6 0.95

1

1.05

1.1

1.15

1.2

1.25

1.3

Average objective value

-30645 -30646 -30647 -30648 -30649 -30650 -30651 -30652 -30653 -30654 C Fig. 3. Average objective values versus different levels of the step size parameter for G04.

c = 0.9 c = 1.1 c = 1.3 500 -30520

1000

1500

2000

2500

3000

3500

4000

4500

5000

Average objective value

-30540 -30560 -30580 -30600 -30620 -30640 -30660 -30680 Iteration Fig. 4. Average improvement of the objective function versus the number of iterations for G04 when S¼ 40.

c = 0.9 c = 1.1 c = 1.3 -30520

500

1000

1500

2000

2500

3000

3500

4000

4500

5000

-30540 Average objective value

In this section, we plan to perform a series of sensitivity analysis experiments to observe which of the parameters can essentially affect the total efficiency of our proposed algorithm. By the end of this section, we aim at screening the significant parameters which can seriously influence the conclusions based upon the quality and the corresponding rate of convergence. The relevant outcomes direct us to fine tune the critical parameters. In the course of running the algorithm, the following parameters of step size (c), swarm size (S) and the threshold value of the beta distribution (r) are treated as the critical parameters. To do so, for each parameter, some levels have been considered as follows: three levels of 0.9, 1.1 and 1.3 for the step size, four levels of 40, 60, 80 and 100 for the swarm size and two levels of 0.6 and 0.9 for the threshold value. To trace the effect of changing the levels of the critical parameters, we chose two test problems from Ref. [24]. These are G03 and G04. G03 is a maximization problem while G04 is a minimization problem. G04 is regarded as a relatively simple problem but G03 is thought of as a rather hard problem. Our algorithm was coded in MATLAB 7 on a CORE 2duo @2.4 GHZ notebook and for each combination of the parameter levels, 10 replications were made based on generating random initial swarms. Fig. 2 illustrates the average objective values of the final solutions for G04 obtained by the RPSO-DLM versus different values of the swarm size, for two threshold values of 0.6 and 0.9. As can be seen, the plots of the average objective values at those two levels are quite different. The 0.9 plot fluctuates rather widely while the 0.6 plot behaves monotonically. Besides, the 0.6 plot is indicative of a superior performance when compared to the 0.9 plot. For r ¼ 0.6, while the impact of the swarm size is almost negligible for S A ½40,80, for S480 the results tend to worsen. More importantly, for S 480, we should expect a higher computation time due to the higher number of particles. Therefore, by way of experimentation, we like to conclude that for simple (which will be defined later) problems like G04, r ¼0.6 leads to a better solution. Fig. 3 reveals how much the average objective values can be affected by the variations in the step size parameter for the case of r ¼0.6. This figure displays a decreasing concave graph indicating that more suitable solutions can be obtained for the higher values of c. However, in spite of being effective on the quality, the step size may weaken the algorithm capability in the context of converging to the global optima due to a larger step size that in turn, requires a larger number of iterations. To curb against such undesirable outcomes, we need to consider the effect of the step size as well as the swarm size on the convergence rate. To achieve this purpose,

0.9 -30644

-30560 -30580 -30600 -30620 -30640 -30660 -30680 Iteration

Fig. 5. Average improvement of the objective function versus the number of iterations for G04 when S¼ 60.

S Fig. 2. Average objective values versus different levels of the swarm size parameter for G04.

Figs. 4, 5 and 6 compare three levels of the step size in terms of the gradual rate of improvement in the objective function for three values of the swarm size (40, 60 and 80) when r is fixed at 0.6. As the figures imply, the swarm size has a more pronounced effect on the convergence rate than the effect exerted by the step size. As a consequence, two values of 60 and 80 for the swarm size denote more promising results when compared to the value of 40. In addition, the step size parameter has a trivial effect on the

A. Mohammad Nezhad, H. Mahlooji / Computers & Operations Research 38 (2011) 1164–1174

r ¼ 0:6, 60 rS r80 and 1:1 rc r 1:3 So far, we have considered a simple test problem with a relatively vast feasible region for which a large number of local optima may exist. In hope of finding improved feasible solutions for such problems, a good idea is to rapidly find a feasible solution before the particles’ velocities approach zero because a futile search around the infeasible region is prone to premature convergence. Setting r ¼0.6 achieves this purpose in that a greater effort ought to be spent to penalize the violated constraints (which in turn directs the algorithm toward the feasible region) in contrast tor ¼ 0.9. However, for a typical hard problem with highly nonlinear constraints, the feasible region may be so small that a feasible solution can hardly be found. Since, in this case, the previous interpretations may not be valid, it is imperative to repeat the same analysis for a typical hard problem as well. To fulfill this task, we pick the test problem G03 for evaluation. We employ the same test conditions including the parameter settings exactly as the ones used for G04. Similar to Fig. 2, Fig. 7 depicts the average objective values of RPSO-DLM versus different values of the swarm size, for G03. Contrary to the findings in dealing with the test problem G04, here, the better results correspond to the higher level of the threshold

500 -30520

c = 0.9 c = 1.1 c = 1.3 1000 1500

2000

2500

3000

3500

4000

4500

value, r ¼0.9. Furthermore, the plot with r ¼0.9 in Fig. 7 displays much less fluctuations with the changes in S than the plot with r ¼0.6. For the values of S A ½60,80, the average objective values of the plot with r ¼0.9 seem to be stabilized at or near 1. Fig. 8, in a similar manner, displays the effect of various levels of the step size on the quality of the average objective values. As can be seen, unlike the case of G04, the average objective values for G03 monotonically decrease as c is increased. In the same way as the first set of experiments, in Figs. 9 and 10, we have plotted the gradual rate of improvement versus the number of iterations, considering two levels of 60 and 80 for the

Average objective value

convergence rate so that 0.9, 1.1 and 1.3 lead to the same performance, particularly for S A ½60,80. Hence, based on the above discussion related to Figs. 2 to 6, it is recommended that for simple problems of this type, the parameters be set as

0.9

Average objective value

-30580 -30600 -30620 -30640

0.95

1

1.05

-30660 -30680

1.1 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 1000

ρ = 0.9 1

Average objective value

Average objective value

0.97 0.91 0.88 0.85 0.82 0.79 0.76 0.73 0.7 40

50

60

70 S

80

90

100

Fig. 7. Average objective values versus different levels of the swarm size parameter for G03.

1.2

1.25

1.3

2000

3000

4000

5000 6000 Iteration

7000

8000

9000

10000

Fig. 9. Average improvement of the objective function versus the number of iterations for G03 when S ¼60.

ρ = 0.6

0.94

1.15

c = 0.9 c = 1.1 c = 1.3

Iteration Fig. 6. Average improvement of the objective function versus the number of iterations for G04 when S ¼80.

1.1 C

Fig. 8. Average objective values versus different levels of the step size parameter for G03.

-30540 Average objective value

ρ = 0.9

0.999 0.998 0.997 0.996 0.995 0.994 0.993 0.992 0.991 0.99 0.989 0.988 0.987 0.986 0.985

5000

-30560

1171

1.1 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 1000

c = 0.9 c = 1.1 c = 1.3

2000

3000

4000

5000 6000 Iteration

7000

8000

9000

10000

Fig. 10. Average improvement of the objective function versus the number of iterations for G03 when S ¼80.

1172

A. Mohammad Nezhad, H. Mahlooji / Computers & Operations Research 38 (2011) 1164–1174

swarm size (because of being superior than others) and three levels of 0.9, 1.1 and 1.3 for the step size. Apparently, S ¼60 outperforms S¼ 80 in the sense that it converges more rapidly to the global optima and is more robust to variations of c. However, as the plots demonstrate, for both cases, c ¼0.9 and c¼1.1 show equal and yet better performance in contrast to c ¼1.3. Thus, according to Figs. 7 and 8, along with the conclusions made just above, it seems reasonable to set 60 rS r80 as well as 0.9 rcr1.1 to avoid a slow algorithm. For the so called hard problems like G03, which consist of highly nonlinear constraints, the feasible region is normally too small and hence the number of local optimal solutions is very low. Besides, a non-smooth complicated objective function may exert a negative effect on the convergence rate of the algorithm. All this means that our algorithm may be trapped in a local optimum solution because of over emphasizing on the obtaining a feasible solution. To overcome this dilemma, we suggest an alternative approach which also examines the qualified infeasible solutions to perform a more detailed search operation. As the number of iterations increases, the criterion for defining the solution space becomes more and more difficult for the qualified infeasible solutions to be satisfied. This approach is equivalent to choosing larger values for the threshold parameter (r).

7. Numerical computations In this section, we make an attempt to justify the constraint handling capability of the RPSO-DLM algorithm by posing this question: How promising the proposed constraint handling method can be when compared to the existing constraint handling methods in standard particle swarm optimization (SPSO 2007 [25]) and to what degree may PBFS be effective? Aside from that, we intend to asses the merits of our algorithm against some of the existing efficient optimization algorithms in the literature. To this end, we have selected ten test problems as given by Michalewicz [24]. These problems are listed in Table 1. To obtain the numerical results in this section, our algorithm was coded in MATLAB 7 and the programs were run on a CORETM 2 duo @2.4 GHZ notebook. For all runs, beta (a ¼2, b ¼3) has been chosen as the non-uniform distribution and in agreement with Lemma 2, the absolute function has been used as the H function. Depending on the problem complexity, the swarm size has been set to a value in the range 60rSr80 and 1.1 has been chosen for the step size parameter. Prior to fine tuning the threshold parameter, we should approximate the degree of the problem complexity. Toward this end, a large number of random solutions are generated within the associated bounds of the decision variables. Let P represent the percent of the solutions satisfying the entire set of functional constraints. In the case at least one highly nonlinear constraint (i.e., nonlinear equality constraint, high degree geometric inequality constraint, etc.) exists, Table 1 The description of Michalewicz test problems. Problems

n

Objective function

P (%)

LI

NE

NI

G01 G02 G03 G04 G05 G06 G07 G08 G09 G10

13 20 10 5 4 2 10 2 7 8

Quadratic Nonlinear Polynomial Quadratic Cubic Cubic Quadratic Nonlinear Polynomial Linear

0.0111 99.8474 0.0000 52.1230 0.0000 0.0066 0.0003 0.8560 0.5121 0.0010

9 0 0 0 2 0 3 0 0 3

0 0 1 0 3 0 0 0 0 0

0 2 0 6 0 2 5 2 4 3

n: number of variables; P: ratio of the feasible region; LI: number of linear inequality constraints; NE: number of nonlinear equality constraints; NI: number of nonlinear inequality constraints.

r is set to 0.9. If Po0.1 and the problem includes a non-polynomial objective function, then r ¼0.9. Otherwise, r ¼0.6 should be used. To answer the first question presented in the opening paragraph in this section, we have provided a test bed to compare RPSO-DLM with the standard particle swarm optimization [25] configured by three different constraint handling methods namely the quadratic penalty function (QP) [4], the augmented Lagrange penalty function (ALAG) [4] and the discrete Lagrange multipliers (DLM). Among these three methods, the last one is very similar to the approach adopted in this paper (RPSO-DLM) except for the fact that it employs SPSO and not our revised PSO. To do this, two well-suited engineering design problems [15] plus three problems code named G04, G07 and G10 have been chosen for the numerical evaluation. The engineering design problems involve the tension/compression design problem containing three decision variables, four nonlinear inequality constraints and polynomial objective function as well as the welded beam design problem with four decision variables, seven nonlinear inequality constraints and polynomial objective function. In spite of the default parameter settings of SPSO, in order to improve the efficiency of SPSO-ALAG, SPSO-QP and SPSO-DLM, it was decided to assign a sufficiently large value to the swarm size parameter in hope of generating the best possible final solution. Aside from the swarm size, SPSO-DLM has been prepared in the same way as RPSO-DLM except for the case G04 for which we utilized the uniform distribution in lieu of the beta distribution. The same stopping rule has been applied to terminate the SPSO at each iteration of SPSO-ALAG and SPSO-QP and each of the last two algorithms have been separately fine tuned up to the highest level. Also, a maximum number of 200,000 function evaluations, along with generating the known optimal solution (determined by the LINGO 8 package) have been considered as the two stopping rules in all cases and for all algorithms whichever comes first. Table 2 presents the outcome of the numerical experiments. Each entry in this table is the (rounded) average of 20 runs. The second column displays the optimal objective value for each of the five test problems. The rest of this table presents the best, average, standard deviation as well as the average number of function evaluations (#FE) for each test problem for all four competing algorithms. Except for the tension/compression problem, RPSO-DLM strictly outperforms the others in terms of the best and the average solutions in all test problems. In four out of five cases (tension/compression, welded beam, G04 and G07), RPSO-DLM has reached the optimal solution in all 20 runs. Just one of the 20 runs of RPSO-DLM for the remaining problem, G10, fails to achieve the optimal solution. Nevertheless, RPSO-DLM relies on much fewer number of function evaluations to converge while SPSO-QP and SPSO-ALAG often fail to converge in spite of crossing the permitted number of function evaluations. SPSO-DLM has succeeded in getting the global optima in three cases but at the cost of more number of function evaluations. Besides, for G10, SPSO-ALAG and SPSO-LDM not only have failed to recover the optimal solution but also are not capable of generating even a feasible solution. Hence, from this point of view, RPSO-DLM has had the best performance too. The main advantage of our algorithm can be attributed to the powerful mathematical background of the discrete Lagrange multipliers due to its independence from the problem structure, though it has been undergone a change to be adapted to the continuous space, as pointed out at the beginning of Section 4. This is in contrast to the ALAG and QP methods for which to obtain the global optima requires providing some preliminary assumptions [4]. Although SPSO-DLM appears to work fine, due to using a Lagrangian function framework similar to penalty function approaches, it may be unsuccessful in finding even a feasible solution as the problem becomes intractable. However, the initial results indicate that our proposed strategy, PBFS, can be helpful in getting around this problem.

91605 95228 197020 200080 200100 0 0.00 0.88 1.49 – 0.0126652 1.724858  30664.979 26.0069 –

(1) Constrained simulated annealing (CSA) [5] (2) Stochastic ranking (SR) [26] (3) Improved particle swarm optimization (IPSO) [27] The first algorithm, CSA, relies on a simulated annealing approach configured by the discrete Lagrange multipliers method. The second algorithm, SR, concerns a non-penalty function constraint handling method incorporated in an evolutionary strategy algorithm which attempts to balance between the constraint satisfaction and objective function evaluation. The last algorithm, IPSO, presents a feasibility based rule embedded in an improved particle swarm optimization which employs a number of new operators to rectify the shortages of the standard PSO. To establish a fair framework for judgment, we chose to limit the number of function evaluations to 350,000 and the number of repeat runs to 30 for each test problem. Again, our algorithm will be terminated upon obtaining the global optima. The reason we have selected these limits has to do with the fact that SR and IPSO’s performance have been reported in the literature [26,27] with 350,000 and 340,000 function evaluations respectively, as well as 30 repeat runs. Unlike the cases of SR and IPSO, the literature just contains this piece of information on CSA that the number of repeat runs has been set to 100 in a study performed in 1999 [5]. Table 3 aims in comparing the RPSO-DLM algorithm against CSA, SR and IPSO based upon such indices the best, average and the standard deviation except for CSA for which the standard deviation has not been reported. Also, the fourth column of RPSO-DLM represents the number of function evaluations (#FE). In each row of this table, the best values have been printed in bold numbers. Based on the results contained in Table 3, in 8 out of 10 cases, RPSO-DLM attains the global optima where for 6 of these cases, all the repeat runs finish with finding the global optima. However, CSA, SR and IPSO generate the optimum solution for 10, 6 and 5 test problems respectively. If we do not consider G01 and G05, none of the existing methods is capable of outperforming RPSO-DLM. In fact, RPSO-DLM comes up with inferior results for these two problems despite passing the mark of 350,000 function evaluations. Literature reports on the optimal solution for G05 are sketchy in the sense that two different values are reported for the optimal solution [28] one of which is none other than 5126.498. Interestingly enough, LINGO along with GAMS produce exactly this value as the optimal solution for G05. With regard to the average outcomes, it can be deduced that our proposed algorithm performs better than SR and IPSO for G04, G06, G07, G08, G09 and G10 on the basis of high quality and more accurate solutions. It is worthy to note that for these six test problems, RPSO-DLM attains the global optima in all repeat runs, before meeting the terminating condition. Generally speaking, CSA and RPSO-DLM are the two promising optimization algorithms in this numerical evaluation though care should be taken about this assertion because of the imperfect information about CSA.

61189 234607 93687 224538 – 0 0.01 14.51 7.89 – 0.0126652 1.731805  30662.007 37.2819 – 0.0126652 1.724892  30665.53 28.417 –* 258690 205322 211847 308293 274296 0.00 0.01 16.44 2.46 2006.50 0.0126809 1.730243  30649.672 28.6370 9181.46 0.0126653 1.724852  30664.225 25.483 7088.23 17415 48890 47408 138830 153160 0 0 0 0 0.16 0.0126652 1.724852  30665.53 24.3062 7049.37 0.0126652 1.724852  30665.53 24.3062 7049.33 0.0126652 1.724852  30665.53 24.3062 7049.33

8. Conclusions

* No feasible solution.

Tension compression Welded beam G04(min) G07(min) G10(min)

1173

Next, we identified three of the best existing algorithms for more numerical experiments. These algorithms are:

0.0126652 1.724852  30665.53 24.48925 –

Std Average Best Best Average Best Best

Average

Std

#FE

SPSO-QP Optimal (LINGO) RPSO-DLM Problem

Table 2 Numerical results of SPSO-QP, SPSO-ALAG, SPSO-DLM, and RPSO-DLM on some selected test problems.

Std

#FE

SPSO-ALAG

Average

Std

#FE

SPSO-DLM

#FE

A. Mohammad Nezhad, H. Mahlooji / Computers & Operations Research 38 (2011) 1164–1174

In this paper, a particle swarm approach was introduced to tackle constrained nonlinear programming problems for which the constraints are managed using the discrete Lagrange multipliers method. In contrast to the mathematical theory of the popular exact methods which are essentially dependent on some sensitive assumptions, the discrete Lagrange multipliers theory can be expanded on nonconvex and nondifferentiable optimization problems in the realm of discrete space. Having been extended to the continuous space, the discrete Lagrange multipliers method was

Std

0 0.0159 0.0000 0 0.0005 0.0000 0.0819 0 0.0012 18.4433  15 0.66353 1  30665.53 5126.5  6961.81 24.9613 0.095825 680.668 7377.286

A. Mohammad Nezhad, H. Mahlooji / Computers & Operations Research 38 (2011) 1164–1174

Average

1174

embedded into a constriction coefficient particle swarm optimizer with the aim of finding the elegant saddle points. To alleviate this weakness of PSO in dealing with highly nonlinear constraints, we revised PSO-DLM by using the feasibility strategy of FBPS. Numerical experiments reiterate that our RPSO-DLM algorithm is efficient and robust for the sake of yielding promising results in spite of depending on half a dozen critical parameters.

 15 0.803603 1  30665.53 5126.498  6961.81 24.3345 0.095825 680.6307 7251.2066 0 0.02 0.0002 0.0000 3.5 160 0.066 0 0.0340 530  15 0.7820 1  30665.54 5128.881  6875.94 24.374 0.095825 680.656 7559.192  15 0.803515 1  30665.53 5126.497  6961.81 24.307 0.095825 680.6301 7054.316  15 0.802829 1  30665.53 4221.956  6961.81 24.3062 0.095825 680.6301 7049.33  15 0.803619 1  30665.53 4221.956  6961.81 24.3062 0.095825 680.6301 7049.33 350100 350200 328532 45158 186347 60080 138607 7398 29544 156392 0.5456 0.06197 0.1076 0 2.9391 0 0 0 0 0  14.7765 0.7401799 0.9488  30665.53 5127.347  6961.81 24.3062 0.095825 680.6301 7049.33

Std

 15 unknown 1  30665.53 unknown  6961.81 24.3062 unknown 680.6301 7049.33

Average

G01(min) G02(max) G03(max) G04(min) G05(min) G06(min) G07(min) G08(max) G09(min) G10(min)

 14.987 0.803619 1  30665.53 5126.498  6961.81 24.3062 0.095825 680.6301 7049.33

Best Best Best Best

#FE

CSA RPSO-DLM Global optimum (LINGO) Problem ID

Table 3 Numerical results of CSA, SR, IPSO, and RPSO-DLM for the Michalewicz classic test problems.

Average

SR

Average

Std

IPSO

References [1] Zoutendijk G. Methods of feasible directions: a study in linear and non-linear programming. Amsterdam: Elsevier; 1960. [2] Rosen JB. The gradient projection method for nonlinear programming. Journal of the Society for Industrial and Applied Mathematics 1960;8(1):181–217. [3] Wolf P. Methods of nonlinear programming. Recent advances in mathematical programming. New York: McGrew-Hill; 1963. [4] Bazaraa MS, Sherali HD, Shetty LM. Nonlinear Programming: theory and algorithms. New York: John Wiley and Sons; 2006. [5] Wah BW, Wang T. Constrained simulated annealing with applications in nonlinear continuous constrained global optimization. In: Proceedings of the 11th IEEE international conference on tools with artificial intelligence, 1999. p. 381–8. [6] Wah BW, Chen Y, Wang T. Simulated annealing with asymptotic convergence for nonlinear constrained optimization. Journal of Global Optimization 2007;39(1):1–37. [7] Eberhart RC, Kennedy J. A new optimizer using particle swarm theory. In: Proceedings of the sixth international symposium on micro machine and human science, 1995. p. 39-43. [8] Parsopoulos KE, Vrahatis MN. Particle swarm optimization method for constrained optimization problems. In: Proceedings of the second Euro-international symposium on computational intelligence, Slovakia, 2002. p. 214-20. [9] Coello CA. Use of a self-adaptive penalty approach for engineering optimization problems. Computers in Industry 2000;41(2):113–27. [10] He Q, Wang L. An effective co-evolutionary particle swarm optimization for constrained engineering design problems. Engineering Applications of Artificial Intelligence 2007;20(1):89–99. [11] Eberhard P, Sedlaczek K. Using augmented Lagrangian particle swarm optimization for constrained problems in engineering. Advanced Design of Mechanical Systems: From Analysis to Optimization 2009:253–71. [12] He Q, Wang L. A hybrid particle swarm optimization with a feasibility-based rule for constrained optimization. Applied Mathematics and Computation 2007;186(2):1407–22. [13] Dong Y, Tang J, Xu B, Wang D. An application of swarm optimization to nonlinear programming. Computers & Mathematics with Applications 2005;49(11–12):1655–68. [14] Chootinan P, Chen A. Constraint handling in genetic algorithms using a gradientbased repair method. Computers & Operations Research 2006;33(8):2263–81. [15] Zahara E, Kao YT. Hybrid Nelder–Mead simplex search and particle swarm optimization for constrained engineering design problems. Expert Systems with Applications 2009;36(2):3880–6. [16] Pedamallu CS, Ozdamar L. Investigating a hybrid simulated annealing and local search algorithm for constrained optimization. European Journal of Operational Research 2008;185(3):1230–45. [17] Wah BW, Wu Z. The theory of discrete Lagrange multipliers for nonlinear discrete optimization. In: Principles and practice of constraint programming– CP’99, 1999. p. 28–42. [18] Zhe W. Discrete Lagrangian methods for solving nonlinear discrete constrained optimization problems. M.Sc. thesis, University of Illinois, USA; 1998. [19] Avriel M. Nonlinear programming: analysis and methods. New Jersey: Prentice-Hall; 1976. [20] Shi Y, Eberhart RC. Modified particle swarm optimizer. In: The 1998 IEEE international conference on evolutionary computation, ICEC’98, 1998. p. 69–73. [21] Clerc M. The swarm and the queen: towards a deterministic and adaptive particle swarm optimization. In: Proceedings of IEEE congress on evolutionary computation, Washington, D.C., USA, 1999. p. 1951–7. [22] Clerc M, Kennedy J. The particle swarm-explosion, stability, and convergence in a multidimensional complex space. IEEE Transactions on Evolutionary Computation 2002;6(1):58–73. [23] Van den Bergh F. An analysis of particle swarm optimizers. Ph.D. thesis, University of Pretoria, South Africa; 2001. [24] Koziel S, Michalewicz Z. Evolutionary algorithms, homomorphous mappings, and constrained parameter optimization. Evolutionary Computation 1999;7(1):19–44. [25] Standard particle swarm optimization: /http://www.particleswarm.info/Pro grams.html/S. [26] Runarsson TP, Yao X. Stochastic ranking for constrained evolutionary optimization. IEEE Transactions on Evolutionary Computation 2000;4(3):284–94. [27] Sun C, Zeng J, Pan J. An improved particle swarm optimization with feasibilitybased rules for constrained optimization problems. Next-Generation Applied Intelligence 2009:202–11. [28] Zhu W, Ali MM. Solving nonlinearly constrained global optimization problem via an auxiliary function method. Journal of Computational and Applied Mathematics 2009;230(2):491–503.