On a smoothed penalty-based algorithm for global optimization

Journal of Global Optimization On a smoothed penalty-based algorithm for global optimization --Manuscript Draft-- Manuscript Number: JOGO-D-16-00039R2 Full Title: On a smoothed penalty-based algorithm for global optimization Article Type: S.I. : EURO-2015 Keywords: Global optimization; Penalty function; Artificial fish swarm; Markov chains. Corresponding Author: Ana Maria A.C. Rocha, Ph.D. Algoritmi Research Centre, University of Minho Braga, PORTUGAL Corresponding Author Secondary Information: Corresponding Author's Institution: Algoritmi Research Centre, University of Minho Corresponding Author's Secondary Institution: First Author: Ana Maria A.C. Rocha, Ph.D. First Author Secondary Information: Order of Authors: Ana Maria A.C. Rocha, Ph.D. M. Fernanda P. Costa, Ph.D. Edite M.G.P. Fernandes, Aggr. Order of Authors Secondary Information: Funding Information: Abstract: This paper presents a coercive smoothed penalty framework for nonsmooth and nonconvex constrained global optimization problems. The properties of the smoothed penalty function are derived. Convergence to an $\epsilon$-global minimizer is proved. At each iteration $k$, the framework requires the $\epsilon^{(k)}$-global minimizer of a subproblem, where $\epsilon^{(k)} \rightarrow \epsilon$. We show that the subproblem may be solved by well-known stochastic metaheuristics, as well as by the artificial fish swarm (AFS) algorithm. In the limit, the AFS algorithm convergence to an $\epsilon^{(k)}$- global minimum of the real-valued smoothed penalty function is guaranteed with probability one, using the limiting behavior of Markov chains. In this context, we show that the transition probability of the Markov chain produced by the AFS algorithm, when generating a population where the best fitness is in the $\epsilon^{(k)}$-neighborhood of the global minimum, is one when this property holds in the current population, and is strictly bounded from zero when the property does not hold. Preliminary numerical experiments show that the presented penalty algorithm based on the coercive smoothed penalty gives very competitive results when compared with other penalty-based methods. Response to Reviewers: Dear Editor All the issues have been addressed. Best regards, The authors Powered by Editorial Manager® and ProduXion Manager® from Aries Systems Corporation Journal of Global Optimization manuscript No. (will be inserted by the editor) On a smoothed penalty-based algorithm for global optimization Ana Maria A.C. Rocha ·M. Fernanda P. Costa · Edite M.G.P. Fernandes Received: date / Accepted: date Abstract This paper presents a coercive smoothed penalty framework for nonsmooth and nonconvex constrained global optimization problems. The properties of the smoothed penalty function are derived. Convergence to an ε-global minimizer is proved. At each iteration k, the framework requires the ε(k)-global minimizer of a subproblem, where ε(k)→ε. We show that the subproblem may be solved by well-known stochastic metaheuristics, as well as by the artificial fish swarm (AFS) algorithm. In the limit, the AFS algorithm convergence to an ε(k)-global minimum of the real-valued smoothed penalty function is guaranteed with probability one, using the limiting behavior of Markov chains. In this context, we show that the transition probability of the Markov chain produced by the AFS algorithm, when generating a population where the best fitness is in the ε(k)-neighborhood of the global minimum, is one when this property holds in the current population, and is strictly bounded from zero when the property does not hold. Preliminary numerical experiments show that the presented penalty algorithm based on the coercive smoothed penalty gives very competitive results when compared with other penalty-based methods. Keywords Global optimization ·Penalty function ·Artificial fish swarm ·Markov chains Ana Maria A.C. Rocha Algoritmi Research Centre, Department of Production and Systems, University of Minho, Campus de Gualtar, 4710-057 Braga, Portugal E-mail: [email protected] M. Fernanda P. Costa Centre of Mathematics, Department of Mathematics and Applications, University of Minho, Campus de Gualtar, 4710-057 Braga, Portugal E-mail: [email protected] Edite M.G.P. Fernandes Algoritmi Research Centre E-mail: [email protected] Manuscript Click here to download Manuscript Smoothed-penaltyAFS_final_2017.tex Click here to view linked References 2 A.M.A.C. Rocha et al. 1 Introduction This paper aims to contribute to the research area concerned with smoothed penalty function methods in constrained global optimization (CGO), by using well-established metaheuristics to compute a sequence of approximations to the solution of the CGO problem. The problem to be addressed is of the form: min x∈Ωf(x) subject to gi(x)≤0,i=1,...,p(1) where f:Rn→Rand gi:Rn→R,i=1,...,pare nonlinear continuous functions, possibly nondifferentiable, and Ω={x∈Rn:−∞<l≤x≤u<∞}. When the problem has equality constraints, h(x) = 0, they may be reformulated into the above form by using a couple of inequality constraints h(x)−υ≤0 and −h(x)−υ≤0, where υis a small positive relaxation parameter. Functions fand gmay be nonconvex and many local minima may exist in the feasible region. When solving global optimization problems, penalty-type methods combined with stochastic techniques have been appearing in the literature [1–7]. They are simple to implement and provide in general high quality solutions. Penalty functions require in general the use of a positive penalty parameter that aims to balance function and constraint violation values. Setting the initial value for the penalty parameter and tuning its values throughout the iterative process are not easy tasks. Furthermore, the performance of the algorithm is greatly affected by their values. In some cases, the optimal solution of the problem is attained only when the penalty parameter approaches infinity. This is the case with the quadratic penalty [8]. Popular penalties when combined with stochastic heuristics are the dynamic penalty, in which the parameter depends continuously on the iteration counter [1,3,5,6], and the adaptive penalty function. These make the performance of the algorithm less sensitive to penalty variations [2,9,7]. Penalties that rely on two parameters have appeared in the literature. A two-parameter hyperbolic penalty is presented and the convergence properties are discussed in [10]. A family of non-coercive smoothed penalty functions and an algorithm for nonlinear constrained optimization are presented in [11]. The convergence properties are also shown. This last work is based on [12] where the therein studied penalty functions can be rewritten like those in [11]. Exact penalty functions are very efficient since every global minimizer of the penalty function is a global solution of the related constrained problem, and conversely, for some finite values of the penalty. In [13], an improved version of the well-known global deterministic DIRECT algorithm, tailored for global optimization with simple bounds, is incorporated into an exact penalty approach. It is proved that every accumulation point of a sequence of iterates produced by the algorithm is a global minimizer of the constrained problem. Furthermore, under a weak assumption, it is proved that the penalty parameter is updated a finite number of times. In this paper, we aim to further explore penalty-based approaches for solving CGO problems. The contribution of the present study is a coercive smoothed penalty framework for nonsmooth and nonconvex CGO problems. The penalty added to the objective function is a smoothed penalty function and depends on two parameters, one is a penalty weight for constraint violation and the other is a smoothing parameter. Further, convergence of a sequence of iterates to an ε-global minimizer is proved. At each iteration k, the penalty framework requires an ε(k)-global minimizer of a bound constrained optimization subproblem, where ε(k)→ε. The subproblems are solved by well-known stochastic metaheuristics. One particular metaheuristic – the artificial fish swarm (AFS) algorithm [14] – is further analyzed Smoothed penalty-based algorithm 3 and its convergence to an ε(k)-global minimum of the real-valued smoothed penalty function is guaranteed with probability one, using the limiting behavior of Markov chains. We also show that the proposed smoothed penalty framework, when a selected set of well-known metaheuristics (including the AFS algorithm) are used to solve the subproblems, is effective in finding global optimal solutions to the CGO problem. Metaheuristics are general-purpose, usually stochastic and population-based, and define a set of solutions using nature-inspired operations. They are important and effective strategies to solve hard optimization problems. This type of problems cannot be solved to optimality by any exact method within a reasonable time limit [15]. Metaheuristics are specially tailored for black-box optimization problems, where there is no explicit form for the involved functions. Objective and constraint functions are evaluated through physical or computational simulation, thus preventing the use of automatic differentiation techniques to obtain the derivatives. Numerical differentiation is also avoided due to high cost of function evaluations. Challenging problems of this class are optimal control problems, simulationbased optimization and parameter estimation over differential equations, where a function model has to be obtained by sampling or by a third party software library. In general, metaheuristics are simple to implement and use, do not require transformation of the original problem, perform quite well and generate good quality solutions in less time than the traditional optimization techniques. However, when derivative information is available, gradientbased algorithms are substantially superior in terms of convergence speed since fast global convergence may be guaranteed. With stochastic metaheuristics only stochastic convergence can be ensured. A state-of-art concerning metaheuristics in optimization is presented in [16]. Interesting and challenging issues related to the use of metaheuristics are addressed in [17]. The paper is organized as follows. In Section 2, the new smoothed penalty function is presented, and the penalty-based algorithm with its convergence properties are derived. In Section 3, we refer to the metaheuristics that have been selected to solve the bound constrained optimization subproblem and discuss the asymptotic convergence properties of the AFS algorithm, Section 4 presents some numerical experiments and we conclude the paper in Section 5. 2 Penalty Method The feasible set of the problem (1) is defined as F={x∈Ω⊂Rn:gi(x)≤0,i=1,...,p}.(2) In this paper, we aim at converging to a global solution of the nonlinear optimization problem (1) using a penalty framework. In [11], a smoothed penalty function that depends on a smooth parameter µ∈(0,1], and is defined as Pµ(t) = µPt µ, where P(t) = log(1+exp(t)) (3) is used, in the sense that a differentiable function is used by smoothing the function t+= max{0,t}(see also [12,18]). Thus, in [11], the penalty term associated with the constraint gi(x)≤0 has the form Pµ(gi(x)) = µlog1+expgi(x) µ (4) 4 A.M.A.C. Rocha et al. and aims to penalize constraints violation of the problem (1). Thus, at each iteration k, where kdenotes the iteration counter of the outer cycle, and for fixed values of the smoothing parameter, µ(k)∈(0,1], and of a penalty weight, τ(k)≥1, the subproblem takes the form: min x∈Ωφ(x;τ(k),µ(k))≡f(x)+τ(k)µ(k) p ∑ i=1 Pgi(x) µ(k).(5) Hence, this penalty method consists in solving a family of bound constrained optimization subproblems of the form (5) where µ(k)→0 and the function Psatisfies certain properties. It is proved in [12] that a path of optimal solutions {xµ(k)}µ(k)>0defined by xµ(k)≡ argminx∈Ωφ(x;τ(k),µ(k))converges towards the optimal set of the original convex constrained optimization problem (likewise (1)) when µ(k)→0. We note that the parameter µ aims to control the precision of the smoothing and τaims to change the inclination of τPµ(t). Thus, the method used for solving the subproblems should ensure that the bound constraints are always satisfied and global optimal solutions are obtained. Some interesting differences between penalty algorithms are located on the framework used to find an approximate solution to the subproblem (5). Our work is inspired by the smoothed penalty functions studied in [11]. However, we consider nonconvex optimization problems and we do not assume that the objective fand the constraints gi,i=1,...,pare differentiable. Thus, we do not assume that the Mangasarian-Fromovitz condition holds. We use the following definition: Definition 1 (ε(k)-global minimizer) Let φk(x)≡φ(x;τ(k),µ(k))be a continuous objective function defined over a bounded space Ω⊂Rn. The point x(k)∈Ωis an ε(k)-global minimizer of the subproblem (5) if φk(x(k))≤miny∈Ωφk(y) + ε(k), where ε(k)>0 is the error bound which reflects the accuracy required for the approximation. Our penalty-based algorithm requires that at each iteration kof the outer cycle, an ε(k)-global solution of the subproblem (5) is obtained. When the optimization problem is nonconvex, a global optimization method is required to approximately solve the subproblem (5), so that the algorithm has some guarantee to converge to a global solution instead of being trapped in a local one. To compute an ε(k)-global minimizer of the subproblem (5), for fixed values of τ(k)and µ(k), we propose the use of well-known stochastic metaheuristics. 2.1 Smoothed penalty function We now present a new smoothed penalty term for the subproblem described in (5): Pµ(gi(x)) ≡µPgi(x) µ=µ1+θgi(x) µ (6) where the function θ(t)is defined, for t∈R, by θ(t) = tanh(t),if t≤0, sinh(t),otherwise, being P(t) = 1+θ(t).(7) Figure 1 displays plots of P(t),P0(t) = θ0(t)as well as of the penalty term P(t)defined in (3), and its first derivative, for comparison. As previously defined, we have Pµ(t) = µPt µwhere P:R→R+. For any µ∈(0,1] and t∈R, the proposed Pµ(t)defined by (6) with θ(t)as in (7) satisfy the following properties. Smoothed penalty-based algorithm 5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 0 0.5 1 1.5 2 2.5 3 3.5 4 t P(t) = (1+θ(t)) defined in (7) P’(t) P(t) =log(1+exp(t)) defined in (3) exp(t)/(1+exp(t)) Fig. 1 Plots of P(t)in (3) and (7), and first derivatives Property 1 Pµ(t)is convex, continuous and differentiable as a function of t, and P0 µ(0)>0. Proof Function θ(t)is convex, P(t)is convex and since µ∈(0,1],Pµ(t)is also convex. Similarly Pµ(t)is differentiable since θ(t)is differentiable. Since P0 µ(t) = µθ0t µ, P0 µ(t) =    1/cosh2t µ,if t≤0, cosht µ,otherwise, (8) and we get P0 µ(0) = 1. Furthermore, P0 µ(0+) = 1. ut Property 2 limt→−∞P(t) = 0, limt→−∞Pµ(t) = 0 and limt→0−Pµ(t) = µ. Proof lim t→−∞P(t) = lim t→−∞(1+tanh(t)) = 0 and lim t→−∞Pµ(t) = lim t→−∞µ1+tanht µ=0. Furthermore, lim t→0−Pµ(t) = lim t→0−µ1+tanht µ=µ. ut Property 3 limt→∞ Pµ(t) t= +∞and limt→∞ Pµ(−t) t=0. Proof Recalling the definition of θin (7), limt→∞ Pµ(t) t=limt→∞ µ t1+sinht µ=limt→∞cosht µ =ν>0 where ν= +∞, defining a coercive penalty. To prove the second part of the property, we use (7) and: limt→∞ Pµ(−t) t=limt→∞ µ t1+tanh−t µ=0. ut 6 A.M.A.C. Rocha et al. Property 4 For µ∈(0,1],P0 µ(t)≤P0(t)when t≤0 and P0 µ(t)≥P0(t)when t>0. Proof For µ∈(0,1], we have t µ≤twhen t≤0, and t µ≥twhen t>0. From (8), we have (recalling that the hyperbolic cosine is strictly decreasing for negative arguments and strictly increasing for positive arguments) for t≤0: P0 µ(t) = 1/cosh2t µ≤1/cosh2(t) = P0(t) and for t>0 P0 µ(t) = cosht µ≥cosh(t) = P0(t). ut Figure 2 displays the plots of Pµ(t),P(t)and P(t)−P(0)as well as their first derivatives, to check Properties 4 and 5. −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −1 0 1 2 3 4 t P(t) Pµ (t) (for µ=0.5) P’(t) P’µ(t) (for µ=0.5) P(t) − P(0) Fig. 2 Plots of Pµ(t),P(t)−P(0),P0 µ(t)and P0(t) Property 5 For t∈Rand µ∈(0,1],Pµ(t) = µPt µ>P(t)−P(0). Proof We note that P(0) = 1 and P(t)−P(0) = θ(t). Consider µ∈(0,1]. When t≤0, t µ≤t; on the other hand, when t>0, t µ≥t. First, we consider the case t≤0, where 1+tanht µ>0≥tanh(t) (recalling that −1<tanh(t)≤0 for t≤0) and thus µ1+tanht µ>tanh(t). Now, to prove that Pµ(t)>P(t)−P(0) = θ(t)is true for t>0, we use P0 µ(t)≥P0(t)(when t>0) from Property 4, and get Zt 0 coshx µdx ≥Zt 0 cosh(x)dx Smoothed penalty-based algorithm 7 which implies µsinht µ≥sinh(t) and therefore µsinht µ+µ>sinh(t). ut Property 6 For t∈R,Pµ(t)converges pointwise to νt+=νmax{0,t}when µ→0+where ν= +∞(using the convention that ∞×0=0). Proof First, we consider the case t≤0, where lim µ→0Pµ(t) = lim µ→0µ1+tanht µ=0=νt+. However, when t>0, using the change of variable z=t µ, and Property 3 we get limµ→0Pµ(t) = limµ→0µ1+sinht µ=limz→∞ t z(1+sinh(z)) =limz→∞ t z+tlimz→∞ sinh(z) z =tlimz→∞cosh(z) =tν |{z} with ν=+∞ =νt+ ut Figure 3 shows the plots of Pµ(t)for values of µ=1,0.5,0.25 and 0.125, as well as the plots of penalty (3) (identified with Plog in the legend) for µ=1 and 0.125. −1.5 −1 −0.5 0 0.5 1 1.5 0 0.5 1 1.5 2 2.5 3 t Pµ(t) µ=1 (corresponds to P(t)) µ=0.5 µ=0.25 µ=0.125 µ=1 (Plog) µ=0.125 (Plog) Fig. 3 Behavior of Pµ(t)for four different values of µ We now define the penalty term in (5) as pµ(x) = p ∑ i=1 Pµ(gi(x)) = µ p ∑ i=1 Pgi(x) µ(9) and if the algorithm never increases the value of the product τµ, we have: 8 A.M.A.C. Rocha et al. Property 7 If τµ is bounded, then limµ→0+τpµ(x) = 0, for any xin the interior of F. Proof For any xin the interior of F,gi(x)<0, i=1,...,pand lim µ→0+τpµ(x) = lim µ→0+τµ p ∑ i=11+θgi(x) µ=lim µ→0+τµ p ∑ i=1 Pgi(x) µ=0 by Property 2 and boundedness of τµ.ut Property 8 If the sequence {τµ} → 0, then limµ→0+τpµ(x) = 0, for any x∈F. Proof For any x∈F: lim µ→0+τpµ(x) = lim µ→0+τµ p ∑ i=11+θgi(x) µ=lim µ→0+τµ p ∑ i=1 Pgi(x) µ=0, since {τµ} → 0 and Pgi(x) µis bounded when gi(x)≤0. ut 2.2 The penalty-based algorithm A formal description of the penalty algorithm for solving the original problem (1) is presented in Algorithm 2.1. We remark the following: Remark 1 When finding the global minimum of a continuous objective function f(x)over a bounded space F⊂Rn, the point ¯x∈Fis an ε-global minimizer if f(¯x)≤miny∈Ff(y)+ε, where εis the error bound which reflects the accuracy required for the solution. We analyze feasibility at iteration kusing the vector g+, componentwise defined by g+ i(x(k)) = maxn0,gi(x(k))o,i=1,...,p.(10) Algorithm 2.1 Penalty algorithm Require: τ(1)≥1, 0 <µ(1)≤1, γτ>1, 0 <γµ<(γτ)−1,γg<1, γε<1, 0 <εmin 1, ηmin 1, ε(1)>εmin, kmax,LB 1: Set k=1, randomly generate a set of solutions in Ωand select x(0) 2: while kg+(x(k−1))k>ηmin or f(x(k−1))>LB +εminand k≤kmax do 3: Find an ε(k)-global minimizer of subproblem (5), x(k), so that φ(x(k);τ(k),µ(k))≤φ(x;τ(k),µ(k))+ ε(k)for all x∈Ω(11) 4: Set µ(k+1)=γµµ(k) 5: if kg+(x(k))k ≤ γgkg+(x(k−1))kthen 6: Set τ(k+1)=τ(k),ε(k+1)=maxnεmin,γεε(k)o 7: else 8: Set τ(k+1)=γττ(k),ε(k+1)=ε(k) 9: end if 10: Set k=k+1 11: end while Smoothed penalty-based algorithm 15 where 1IA(y(t))denotes the indicator function for the set A, which returns 1 if y(t)∈A, and returns 0 otherwise. The selection kernel in (17) is interpreted as follows. If the trial point y(t)∈Eis better than or equal to the current point x(t), i.e. if y(t)∈B(x(t)), and is also in the set A, then it will be a current point for the next population x(t+1)←y(t)(y(t)will transition to the set A, in particular to the set A∩B(x(t))with probability one). However, if y(t)is worse than x(t), i.e. if y(t)∈Bc(x(t)), then the current point will be preserved: x(t+1)←x(t). The stochastic kernel of the algorithm is then: K(x(t),A) = ZE Km(x(t),dy)1IA∩B(x(t))(y(t))+1IA(x(t))ZE Km(x(t),dy)1IBc(x(t))(y(t)) =ZA∩B(x(t)) Km(x(t),dy)+ 1IA(x(t))ZBc(x(t)) Km(x(t),dy) =Km(x(t),A∩B(x(t)))+1IA(x(t))Km(x(t),Bc(x(t))). (18) We now address the case of a population with more than one point (with E=Ωmand m>1). The set of states better than or equal to the state X(t)is re-defined, according to (14), as B(X(t)) = {Y(t)∈E:φk best(Y(t))≤φk best(X(t))}. If Y(t)∈Eis in A∩B(X(t)), the population transitions to A. However, if Y(t)∈Bc(X(t)), meaning that the best point of the population Y(t)is worse than the best point of population X(t), the entire population Y(t)may be rejected and all points in X(t)are selected. This is equivalent to the above structure with one point in the population (see (18)). In the other case, where the selected population X(t+1)will have some trial points from Y(t)as well as points from X(t), the best point in X(t)is in the selected population and consequently X(t+1)∈B(X(t)). In the sequence of expression (18), the kernel has the following structure: K(X(t),A) = Km(X(t),A∩B(X(t))) + 1IA(X(t))ZBc(X(t)) Km(X(t),dY )1IA(X(t+1)). (19) If we restrict the analysis to the set Aε(k)of ε(k)-optimal states, then using (19) we obtain for the general case: i) if Aε(k)⊂B(X(t))then X(t)/∈Aε(k),Aε(k)∩B(X(t)) = Aε(k)and 1IAε(k)(X(t)) = 0; thus K(X(t),Aε(k)) = Km(X(t),Aε(k)); ii) if B(X(t))⊆Aε(k)then X(t)∈Aε(k),Aε(k)∩B(X(t)) = B(X(t))and 1IAε(k)(X(t)) = 1; thus K(X(t),Aε(k)) = Km(X(t),B(X(t)))+ZBc(X(t)) Km(X(t),dY )1IAε(k)(X(t+1)) =Km(X(t),B(X(t)))+Km(X(t),Bc(X(t))) = 1 since X(t+1)∈B(X(t))⊆Aε(k)and 1IAε(k)(X(t+1)) = 1. Therefore the stochastic kernel restricted to the set Aε(k)is given by K(X(t),Aε(k)) = Km(X(t),Aε(k))1IAc ε(k)(X(t))+ 1IAε(k)(X(t)) satisfying the preconditions of Lemma 1 if Km(X(t),Aε(k))≥λ>0, a value strictly bounded from zero, although it may be ε(k)dependent, for all X(t)∈Ac ε(k). 16 A.M.A.C. Rocha et al. Thus, we must guarantee that the set Aε(k)can be reached from everywhere outside Aε(k) with some probability λ>0. The argument is the following. Each point of the population is moved at random and independently. Let Mε(k)={x(t)∈Ω:φk(x(t))−φk,∗≤ε(k)} be the set of ε(k)-optimal solutions, for some ε(k)>0, where φk(x(t))is the penalty function value computed at x(t). If the probability that a point is moved to a point in the set Mε(k)⊆Ω is larger than β≡β(ε(k))>0, then the probability that a point is moved to another point outside Mε(k)is smaller or equal to 1−β. Now, the probability that a population of mpoints moves to another one which is outside Aε(k), i.e., to a population whose best point satisfies φk best(X(t))−φk,∗>ε(k), is smaller or equal to ∏m i=1(1−β)=(1−β)m, since each point moves at random and independently. Therefore, the probability that a population enters the set Aε(k)∈Ais at least λ≡1−(1−β)m>0. Since β>0 implies λ>0, we must have Km(x(t),Mε(k))≥β>0. We note here that the uniform distribution used to create the trial point or, for that matter, any other distribution with density function that is nowhere zero over the set Esatisfy the condition. ut Using the results of Lemma 1 and both Theorems 3 and 4, the convergence property of the AFS algorithm follows. (This result corresponds to Theorem 2 in [31].) Theorem 5 [31] Assume that the objective function of subproblem (5),φk:Ω→R, is bounded below. Then, the stochastic AFS algorithm will converge to an ε(k)-global minimizer of φk, in the sense of Definition 4, regardless of the initial distribution. 4 Numerical Experiments For a preliminary practical validation of the proposed algorithm based on the smoothed penalty function, two sets of benchmark constrained global optimization problems are used. The C programming language is used in this real-coded algorithm and the numerical experiments were performed on a PC with a a 2.7 GHz Core i7-4600U and 8 Gb of memory. 4.1 Using metaheuristics to solve subproblem (5) First, we analyze the performance of Algorithm 2.1 when different metaheuristics are used to solve subproblem (5). For this experiment, the following parameter values are used. Due to their restrictive conditions, the parameters τ(1)and µ(1)have been naturally set both to one. The parameter γgis set to 0.5, the most used value in the literature. Parameters γτand γµmust be related so that Property 8 holds and after fixing γτ=2, γµis set to 0.25. Other parameters are set as follows: ε(1)=1, γε=0.1, εmin =10−5,ηmin =10−6 and kmax =20. For population-based algorithms, we set m=50 and a maximum number of iterations, tmax =100, is allowed; in the SA we set tmax =2500. Each problem was solved 30 times. The metaheuristics GA, jDE, CMA-ES, ACO, SA and AFS are tested using the parameter values as suggested in the above cited papers. We use 11 problems that have only inequality constraints and are a subset of the well-known g-suite [32]. We note that g02, g08 and g12 are maximization problems that were converted into minimization ones. We report the number of the problem, ‘Prob.’, the median (as a measure of the central tendency of the distribution) of the 30 solutions, ‘ fmedian’, and the average number of function evaluations, Smoothed penalty-based algorithm 17 Table 1 Results produced by Algorithm 2.1 using GA, jDE, CMA-ES, ACO, SA and AFS algorithms to solve subproblem (5) Prob. LB GA jDE CMA-ES ACO SA AFS g01 -15.000000 fmedian -14.999992 -14.998716 -14.999944 -14.999991 -14.852036 -14.999994 N f eavg 75969 101001 102021 74023 90210 88770 g02 -0.803619 fmedian -0.757711 -0.800193 -0.717933 -0.754064 -0.360098 -0.540234 N f eavg 98199 96214 102021 96061 87034 171397 g04 -30665.53867 fmedian -30665.52283 -30665.53568 -30662.61083 -30665.53798 -30664.91162 -30665.47540 N f eavg 102933 100994 44955 97528 100061 112993 g06 -6961.813876 fmedian -6912.430838 -6943.599554 -6958.370917 -6961.673093 -6961.312079 -6961.075018 N f eavg 103001 101001 41083 101001 100061 106871 g07 24.306209 fmedian 25.136433 26.269085 24.325982 26.162025 24.588419 24.337138 N f eavg 103001 101001 102021 101001 100061 136988 g08 -0.095825 fmedian -0.095824 -0.095823 -0.095825 -0.095824 -0.095823 -0.095822 N f eavg 11534 11916 21085 10736 15367 10951 g09 680.630057 fmedian 680.666794 680.655094 680.630513 680.632149 680.690420 680.721112 N f eavg 103001 101001 102021 101001 100061 125957 g10 7049.248021 fmedian 7243.461699 7514.086284 7135.130208 8550.501792 7541.187987 7826.758564 N f eavg 103001 101001 102021 101001 100061 129351 g12 -1.000000 fmedian -1.000000 -1.000000 -1.000000 -1.000000 -0.999990 -1.000000 N f eavg 7599 5269 19045 7013 63425 6661 g18 -0.866025 fmedian -0.865728 -0.828130 -0.866001 -0.833600 -0.855751 -0.866016 N f eavg 103001 101001 102021 101001 100061 101619 g24 -5.508013 fmedian -5.508009 -5.508007 -5.508011 -5.508008 -5.508006 -5.508005 N f eavg 66703 64454 70225 52894 83331 64589 ‘N f eavg’, after the 30 runs. For comparative purposes, the best-known solution available in the literature [32], ‘LB’, is shown. From the results summarized in Table 1, we conclude that CMA-ES comes five times in 1st place (results presented “underlined”) and once in 2nd place (presented in “italic” style), both ACO and AFS come twice in 1st place and twice in 2nd place, jDE comes twice in 1st place and once in 2nd place, GA comes four times in 2nd place and SA comes once in 2nd place. To analyze the statistical significance of the results, we perform a Friedman test. The MatlabTM (Matlab is a registered trademark of the MathWorks, Inc.) function friedman is used. This is a non-parametric statistical test for multiple comparisons to determine significant differences in mean for one independent variable with two or more levels (the results achieved by two or more methods) and a dependent variable (the matched groups taken as the problems) [33]. When the assumption of normality of the data is absent, non-parametric tests are the preferred ones. The null hypothesis in this test is that the mean ranks assigned to the results of the methods under testing are the same. We apply the statistical procedure in the joint analysis of the results of the six distributions of fmedian values. The Friedman’s chi-square value is 11.06 and the corresponding p-value = 0.0502 indicates for a significance level of 5% (or even 1%) that we do not have enough evidence to reject the null hypothesis of “no significant differences on mean ranks”, i.e., the observed differences between the six distributions of fmedian values are not statistically significant. From the results in the table, we can also conclude that the proposed penalty algorithm is effective in finding the global optimal solutions to CGO problems. 18 A.M.A.C. Rocha et al. 4.2 Sensitivity analysis of parameters The Algorithm 2.1 using the AFS metaheuristic for solving the subproblem (5), undergoes a sensitivity analysis on the parameters ε(1)>0, 0 <γε<1 and 0 <γµ<1, which is related to γτ>1 by the condition γµ<(γτ)−1. For this set of experiments, we set kmax =20, m= min{5n,50},tmax =50 and solve each problem 30 times. The effect of the pair (ε(1),γε)on the performance of the Algorithm 2.1 is analyzed by using four different pairs of values (0.1, 0.5), (1, 0.1), (10, 0.05) and (100, 0.01), after fixing γτ=2 and γµ=0.25. We use a graphical procedure to visualize the performance differences among the results produced by these four sets of parameters, in relative terms on the previously selected 11 problems, known as performance profiles [34]. These profiles correspond to (cumulative) distribution functions for a performance metric. The plot reports (on the vertical axis) the percentage of problems solved with the competing variants of the algorithm that is within a certain threshold, t, (on the horizontal axis) of the best result. The performance profiles of the four distributions of the metric fmedian −LB (error of the median fvalue to the LB), for the selected problems, are shown in the Figure 4(a). Figure 4(b) corresponds to the metric N f eavg. The higher the percentage the better. A higher value for t=1 means that the algorithm based on that pair of values achieves the best median values (when fmedian is analyzed) and takes the lowest computational efforts (when N f eavg is analyzed) mostly. Therefore, the pair (1, 0.1) is more successful since it has the highest probability of being the best setting for (ε(1),γε)as far as the median of fvalues is concerned. From the plot, the probability that the pair (1, 0.1) is the winner on a given problem is about 0.45 since it produces the best fmedian in five of the 11 problems. (Each one of the pairs (0.1, 0.5) and (10, 0.05) produces the best fmedian in three out of 11 problems and the pair (100, 0.01) does not produce any best result.) On the other hand, the right end of the profile gives the percentage of the problems that are successfully solved by each variant. From the results relative to the metric N f eavg, we may conclude that the par (1, 0.1) is the best in one of the 11 problems, being the pair (100, 0.01) the best in six out of 11 problems. Thus, the pair that is able to reach better fvalues is on average computationally more expensive. 1 1.5 2 2.5 3 3.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t fmedian − LB ε(1)=0.1; γε=0.5 ε(1)=1; γε=0.1 ε(1)=10; γε=0.05 ε(1)=100; γε=0.01 10 20 30 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t (a) fmedian −LB metric 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t Nfeavg ε(1)=0.1; γε=0.5 ε(1)=1; γε=0.1 ε(1)=10; γε=0.05 ε(1)=100; γε=0.01 (b) N f eavg metric Fig. 4 Performance profiles for four different pairs of (ε(1),γε), using γτ=2 and γµ=0.25 We now aim to analyze the effect of the reduction factor γµon the performance of Algorithm 2.1. Besides 0.25, we also test the values 0.1 and 0.05 (with γτfixed at 2) and γµ=0.05 with γτ=10. From the profiles in the Figure 5(a) we may conclude that γµ=0.25 Smoothed penalty-based algorithm 19 with γτ=2 gives the best median solutions in seven of the 11 problems. The variant with γτ=2 and γµ=0.05 is not able to find a feasible solution to one of the problems and this is noted in its profile that does not reach the top at the right end of the plot. From the Figure 5(b) it is possible to conclude that the variant with γµ=0.25 and γτ=2 also is on average computationally less expensive. 1 1.5 2 2.5 3 3.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t fmedian − LB γτ=2; γµ=0.05 γτ=2; γµ=0.1 γτ=2; γµ=0.25 γτ=10; γµ=0.05 200 400 600 800 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t (a) fmedian −LB metric 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t Nfeavg γτ=2; γµ=0.05 γτ=2; γµ=0.1 γτ=2; γµ=0.25 γτ=10; γµ=0.05 (b) N f eavg metric Fig. 5 Performance profiles for different γµand γτ, using ε(1)=1 and γε=0.1 4.3 Comparison with penalty-based algorithms We now compare our results with those obtained by recent penalty-based algorithms [13, 29,35,36]. To solve the subproblem (5), in the context of the proposed Algorithm 2.1, the AFS algorithm is used. A genetic algorithm based augmented Lagrangian (GAAL) method is proposed in [35], a hyperbolic augmented Lagrangian (HAL) algorithm based on the improved AFS metaheuristic is presented in [29], a shifted hyperbolic augmented Lagrangian (sh-HAL) combined with an enhanced two-swarm AFS algorithm is analyzed in [36], and in [13], an exact penalty (ex-PEN) method based on the DIRECT solver is proposed. Table 2 contains our results and those obtained by GAAL, available in [35]. During these experiments we use γτ=2, γµ=0.25, ε(1)=1, γε=0.1, kmax =20, tmax =500 and consider similar conditions to those reported in the cited paper, as follows: m=max{10n,50}, each problem was solved 25 times and the algorithm terminates when the absolute difference between the function values of two consecutive iterations is less than or equal to 10−4. The other parameters of Algorithm 2.1 are set as previously defined. We report ‘ fbest’ (the best solution obtained after the 25 runs), ‘ fmedian’, ‘ fworst’ (the worst solution of the 25 runs) and the percentage of successful runs, ‘SR’, where a run is considered to be successful if the obtained feasible solution, fsol, satisfies |fsol −LB| ≤ 10−4. The solutions produced by GAAL have high quality [35]. This is an expectable behavior since the GAAL algorithm employs, at the end of each set of five iterations, the gradient-based fmincon local search function (from MatlabTM Optimization Toolbox) starting from the best found solution. We have also invoked a local search procedure, at each iteration, to be able to obtain higher quality solutions. However, a derivative-free local solver named Hooke-and-Jeeves (HJ) is used instead [37]. The search HJ is allowed to run for 50 iterations. We conclude that the results produced by our algorithm based on the smoothed penalty function are very competitive. Applying the Friedman test for the comparison of the distributions of fbest values relative to the Algorithm 2.1 and GAAL, the chi-square statistical value of 2.78 and a p-value of 0.0956 20 A.M.A.C. Rocha et al. are obtained. Thus, there is not enough evidence to reject the null hypothesis of “no significant differences on mean ranks”, at a significance level of 5%. When the Friedman test is applied to the distributions of fmedian values, the chi-square statistical value is 7.36 and the p-value = 0.0067. Hence, the null hypothesis is rejected, at a significance level of 5%, and we conclude that the distributions of fmedian values have statistically significant differences. Table 2 Comparison of Algorithm 2.1 with GAAL Algorithm 2.1 GAAL in [35] Prob. fbest fmedian fworst SR fbest fmedian fworst SR g01 -15.000001 -14.999999 -14.999990 100 -15 -15 -12.453125 96 g02 -0.798167 -0.712454 -0.643880 0 -0.803619 -0.803619 -0.744767 92 g04 -30665.538671 -30665.536658 -30665.424037 20 -30665.538672 -30665.538672 -30665.538672 100 g06 -6961.813263 -6961.687392 -6961.641895 0 -6961.813876 -6961.813876 -6961.813876 100 g07 24.306302 24.306700 24.310923 8 24.306209 24.306209 24.306209 100 g08 -0.095825 -0.095823 -0.095802 100 -0.095825 -0.095825 0 32 g09 680.630143 680.630792 680.636586 8 680.630057 680.630057 680.630057 100 g10 7073.392696 7452.558621 8110.527953 0 7049.24802 7049.24802 7049.24802 100 g12 -1.000000 -1.000000 -1.000000 100 -0.999375 -0.999375 -0.999375 100 g18 -0.866015 -0.866004 -0.865955 100 -0.866025 -0.866025 -0.674981 96 g24 -5.508003 -5.508001 -5.507999 100 -5.508013 -5.508013 -5.508013 100 For the comparison of the Algorithm 2.1 with the HAL algorithm in [29], the sh-HAL algorithm in [36] and the ex-PEN method [13], a different set of 20 small constrained global optimization problems is tested (the number of variables ranges from 2 to 6 and the number of constraints ranges from 1 to 12) [19]. Table 3 lists the number of the problem and the best-known solution available in the literature, ‘LB’, as shown in [19]; the best solution obtained by Algorithm 2.1 (during the 30 runs), ‘ fbest’; the median, ‘ fmedian’; and the number of function evaluations required by the best solution, ‘N f ebest’. The table also displays the results produced by HAL, by sh-HAL, the solution found by the ex-PEN, ‘fsol’ and the number of function evaluations, ‘N f e’ available in [13]. The parameters for this experiment are set as follows: m=max{50,10n},tmax =20, kmax =10 and the parameter υfor the equality constraints of the problems was set to 10−15. From the results we may conclude that Algorithm 2.1 based on the proposed smoothed penalty function performs reasonably well. The comparison with the results in [13] is very favorable. From the comparison with the results in [29], we conclude that Algorithm 2.1 is able to reach best and median solutions similar to its competitor using a lower computational effort. The comparison with the results in [36] is slightly favorable to sh-HAL. The statistical test, for multiple comparisons based on the Friedman test, on the three distributions of fmedian values relative to the Algorithm 2.1, HAL and sh-HAL (penalty-based algorithms combined with variants of the AFS metaheuristic), gives a chi-square statistical value of 6.4, with p-value = 0.0408. Hence, the null hypothesis of “no significant differences on mean ranks” is rejected, at a significance level of 5%, and there is evidence that the three distributions of fmedian values have statistically significant differences. (However, at a significance level of 1% – a smaller probability of making a wrong decision – there is no evidence to reject the null hypothesis.) Smoothed penalty-based algorithm 21 Table 3 Comparison of Algorithm 2.1 with HAL, sh-HAL and ex-PEN results of the Algorithm 2.1 HAL in [29] sh-HAL in [36] ex-PEN in [13] Prob. LB fbest N f ebest fmedian fbest N f ebest fmedian fbest N f ebest fmedian fsol N f e 1 0.029313 0.0293 4719 0.0295 0.0342 9608 0.1204 0.0294 9723 0.0450 0.0625 39575 2(a) -400.00 -362.9380 7565 -313.5667 -380.674 15813 -369.111 -400.0000 2434 -400.0000 -134.1127 115107 2(b) -600.00 -530.9222 7472 -306.3665 -385.051 15808 -360.786 -600.0000 4038 -400.0000 -768.4569 120057 2(c) -750.00 -749.9787 7662 -740.9088 -743.416 15612 -693.743 -750.0000 2433 -750.0000 -82.9774 102015 2(d) -400.00 -399.9999 7627 -399.2205 -399.910 15394 -399.492 -400.0000 2711 -400.0000 -385.1704 229773 3(a) -0.38880 -0.3853 9382 -0.3552 -0.3880 18928 -0.3849 -0.3840 19144 -0.3691 -0.3861 48647 3(b) -0.38881 -0.3888 2295 -0.3888 -0.3888 2589 -0.3888 -0.3888 548 -0.3888 -0.3888 3449 4 -6.6666 -6.6667 2537 -6.6667 -6.6667 2242 -6.6667 -6.6667 698 -6.6667 -6.6666 3547 5 201.16 201.1593 3099 201.1593 201.159 2926 201.159 201.1593 2717 201.1593 201.1593 14087 6 376.29 376.2925 2689 376.4515 376.292 5617 376.293 376.2919 1578 376.2919 0.4701 1523 7 -2.8284 -2.8284 2436 -2.8284 -2.8284 3434 -2.8284 -2.8284 886 -2.8284 -2.8058 13187 8 -118.70 -118.7048 2908 -118.7036 -118.705 2884 -118.705 -118.7049 995 -118.7049 -118.7044 7621 9 -13.402 -13.4019 4508 -13.4017 -13.4018 5732 -13.4017 -13.4019 1437 -13.4019 -13.4026 68177 10 0.74178 0.7418 2646 0.7419 0.7418 6342 0.7418 0.7418 1155 0.7418 0.7420 6739 11 -0.50000 -0.5000 2228 -0.5000 -0.5000 3313 -0.5000 -0.5000 1043 -0.5000 -0.5000 3579 12 -16.739 -16.7389 227 -16.7389 -16.7389 98 -16.7389 -16.7389 267 -16.7389 -16.7389 3499 13 189.35 189.3421 4432 189.3445 189.345 9230 189.347 189.3449 9703 189.3449 195.9553 8085 14 -4.5142 -4.5142 4117 -4.5142 -4.5142 6344 -4.5142 -4.5142 1170 -4.5142 -4.3460 19685 15 0.0000 0.0000 4700 0.0000 0.0000 2546 0.0000 0.0000 3187 0.0000 0.0000 1645 16 0.70492 0.7049 366 0.7049 0.7049 1850 0.7049 0.7049 371 0.7049 0.7181 22593 22 A.M.A.C. Rocha et al. Pairwise comparisons may be carried out to determine which mean ranks are significantly different. The Matlab function multcompare is applied. The estimates of the 95% confidence intervals are shown in the Figure 6(a) for each case under testing. Two compared distributions of fmedian are significantly different if their intervals are disjoint and are not significantly different if their intervals overlap. Hence, we conclude that the mean ranks produced by Algorithm 2.1 is significantly different from that of sh-HAL. For the other two pairs of comparison there are no significant differences on the mean ranks. 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 sh−HAL HAL Algorithm 2.1 fmedian The mean column ranks of Algorithm 2.1 and sh−HAL are significantly different (a) fmedian metric 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 sh−HAL HAL Algorithm 2.1 Nfebest The mean column ranks of HAL and sh−HAL are significantly different (b) N f ebest metric Fig. 6 Estimates of the 95% confidence intervals When the Friedman test is applied to the distributions of fbest values, the obtained chisquare statistical value is 2.4762 with p-value = 0.2899. Thus, there is no evidence that the three distributions of fbest have statistically significant differences. Further, the Friedman test applied to the other metric N f ebest, gives a chi-square statistical value of 10 and a p-value = 0.0067, and we conclude that there is evidence that the three sets of N f ebest have statistically significant differences, at a significance level of 5% (or even 1%). Using multcompare, the estimates of the 95% confidence intervals are shown in the Figure 6(b). It is possible to conclude that the statistically significant differences in the mean ranks are only between HAL and sh-HAL. 5 Conclusions We have demonstrated that a penalty-based algorithm, that uses the coercive smoothed penalty (6), has guaranteed convergence to an ε-global minimum of a CGO problem. The newly developed smoothed penalty is continuous, differentiable and relies on the hyperbolic tangent and hyperbolic sine functions. Approximations to the global optimal solutions of the subproblems are obtained by some well-established and well-known stochastic metaheuristics, as well as by the AFS algorithm. Further, using a relation between the limiting behavior of a Markov chain and the stochastic convergence of random sequences of iterates, we have also proved that the transition probability of the Markov chain produced by the AFS algorithm satisfies the property of Lemma 1 (see also [31]). This property is guaranteed by the movement kernel positiveness and the selection kernel elitism of the AFS algorithm. Then, we can conclude that in the limit, the AFS algorithm convergence to an ε(k)-global minimum of the real-valued smoothed penalty function is guaranteed with probability one, where ε(k)→ε. Smoothed penalty-based algorithm 23 Some preliminary numerical experiments and a comparison with other penalty-based frameworks show that our proposed penalty algorithm is very competitive. A more exhaustive comparison with other techniques for CGO remains to be done. Large dimensional problems will be considered in a near future. Acknowledgments The authors would like to thank two anonymous referees for their valuable comments and suggestions to improve the paper. This work has been supported by COMPETE: POCI-01-0145-FEDER-007043 and FCT - Fundac¸˜ ao para a Ciˆ encia e Tecnologia within the projects UID/CEC/00319/2013 and UID/MAT/00013/2013. References 1. Ali, M.M., Golalikhani, M., Zhuang, J.: A computational study on different penalty approaches for solving constrained global optimization problems with the electromagnetism-like method, Optimization, 63(3), 403–419 (2014) 2. Ali, M.M., Zhu, W.X.: A penalty function-based differential evolution algorithm for constrained global optimization. Comput. Optim. Appl. 54(3), 707–739 (2013) 3. Barbosa, H.J.C., Lemonge, A.C.C.: An adaptive penalty method for genetic algorithms in constrained optimization problems. In: H. Iba (Ed.), Frontiers in Evolutionary Robotics pp. 9–34, I-Tech Education Publ., Austria (2008) 4. Coello, C.A.C.: Theoretical and numerical constraint-handling techniques used with evolutionary algorithms: a survey of the state of the art. Comput. Method. Appl. M. 191(11–12), 1245–1287 (2002) 5. Liu, J.-L., Lin, J.-H.: Evolutionary computation of unconstrained and constrained problems using a novel momentum-type particle swarm optimization. Eng. Optimiz. 39(3), 287–305 (2007) 6. Petalas, Y.G., Parsopoulos, K.E., Vrahatis, M.N.: Memetic particle swarm optimization. Ann. Oper. Res. 156(1), 99–127 (2007) 7. Silva E.K., Barbosa H.J.C., Lemonge A.C.C.: An adaptive constraint handling technique for differential evolution with dynamic use of variants in engineering optimization. Optimization and Engineering, 12(1–2), 31–54 (2011) 8. Bertsekas, D.P.: Nonlinear Programming, 2nd edn. Athena Scientific, Belmont (1999) 9. Mezura-Montes, E., Coello, C.A.C.: Constraint-handling in nature-inspired numerical optimization: Past, present and future. Swarm Evolut. Comput. 1(4), 173–194 (2011) 10. Xavier, A.E.: Hyperbolic penalty: a new method for nonlinear programming with inequalities. Intl. Trans. in Op. Res. 8(6), 659–671 (2001) 11. Gonzaga, C.C., Castillo, R.A.: A nonlinear programming algorithm based on non-coercive penalty functions. Math. Program. Ser. A, 96(1), 87–101 (2003) 12. Auslender, A., Cominetti, R., Haddou, M.: Asymptotic analysis for penalty and barrier methods in convex and linear programming. Math. Oper. Res. 22(1), 43–62 (1997) 13. Di Pillo, G., Lucidi, S., Rinaldi, F.: An approach to constrained global optimization based on exact penalty functions. J. Glob. Optim. 54(2), 251–260 (2012) 14. Rocha, A.M.A.C., Fernandes, E.M.G.P., Martins, T.F.M.C.: Novel fish swarm heuristics for bound constrained global optimzation problems. In: B. Murgante et al. (Eds.) Computational Science and Its Applications – ICCSA 2011, LNCS, Vol. 6784, Part III pp. 185–199 (2011) 15. Boussa¨ ıd, I., Lepagnot, J., Siarry, P.: A survey on optimization metaheuristics. Inf. Sci. 237, 82–117 (2013) 16. Gendreau, M., Potvin, J.-Y. (Eds.): Handbook of Metaheuristics, 2nd edn. International Series in Operations Research and Management Science, Springer, (2010) 17. S¨ orensen, K.: Metaheuristics – the metaphor exposed. Intl. Trans. in Op. Res. 22, 3–18 (2015) 18. J. D. Griffin, J.D., and T. G. Kolda, T.G.: Nonlinearly constrained optimization using heuristic penalty methods and asynchronous parallel generating set search. Appl. Math. Res. Express 2010(1), 36–62 (2010) 24 A.M.A.C. Rocha et al. 19. Birgin, E.G., Floudas, C.A., Mart´ ınez, J.M.: Global minimization using an Augmented Lagrangian method with variable lower-level constraints. Math. Program. Ser. A, 125(1), 139–162 (2010) 20. Holland, J.H.: Adaptation in Natural and Artificial Systems, University of Michigan Press, Ann Arbor, MI, USA, (1975) 21. Brest, J., Greiner, S., Boˇ skovi´ c, B., Mernik, M., ˇ Zumer, V.: Self-adapting control parameters in differential evolution: a comparative study on numerical benchmark problems. IEEE Trans. Evol. Comput. 10, 646–657 (2006) 22. Hansen, N., Ostermeier, A.: Completely derandomized self-adaptation in evolution strategies. Evol. Comput. 9(2), 159–195 (2001) 23. Hansen, N.: The CMA evolution strategy: a comparing review. In Lozano, J.A., Larranaga, P., Inza, I., Bengoetxea, E. (Eds.), Towards a New Evolutionary Computation. Advances on Estimation of Distribution Algorithms, Springer, pp. 75–102 (2006) 24. Socha, K., Dorigo, M.: Ant colony optimization for continuous domains. Eur. J. Oper. Res. 185(3), 1155–1173 (2008) 25. Kirkpatrick, S., Gelatt, C., Vecchi, M.: Optimization by simulated annealing. Science 220, 671–680 (1983) 26. Jiang, M., Wang, Y., Pfletschinger, S., Lagunas, M.A., Yuan, D.: Optimal multiuser detection with artificial fish swarm algorithm. In D.-S. Huang, L. Heutte, M. Loog (Eds.) CCIS 2, ICIC 2007, SpringerVerlag, pp. 1084–1093 (2007) 27. Neshat, M., Sepidnam, G., Sargolzaei, M., Toosi, A.N.: Artificial fish swarm algorithm: a survey of the state-of-the-art, hybridization, combinatorial and indicative applications. Artif. Intell. Rev. 42(4), 965– 997 (2014) 28. Rocha, A.M.A.C., Martins, T.F.M.C., Fernandes, E.M.G.P.: An augmented Lagrangian fish swarm based method for global optimization. J. Comput. Appl. Math. 235(16), 4611–4620 (2011) 29. Costa, M.F.P., Rocha, A.M.A.C., Fernandes, E.M.G.P.: An artificial fish swarm algorithm based hyperbolic augmented Lagrangian method. J. Comput. Appl. Math. 259(Part B), 868–876 (2014) 30. Karr, A. F., Probability, Springer Texts in Statistics, Springer-Verlag, New York, (1993) 31. Rudolph, G.: Convergence of evolutionary algorithms in general search spaces. Proceedings of Third IEEE Conference on Evolutionary Computation, IEEE Press, NJ, pp. 50–54 (1996) 32. Liang, J.J., Runarsson, T.P., Mezura-Montes, E., Clerc, M., Suganthan, P.N., Coello Coello, C.A., Deb, C.: In: Problem Definition and Evolution Criteria for the CEC 2006 Special Session on Constrained Real-Parameter Optimization, IEEE Congress on Evolutionary Computation, Vancouver, Canada, 17–21 July (2006) 33. Derrac, J., Garc´ ıa, S., Molina, D., Herrera, F.: A practical tutorial on the use of nonparametric statistical tests as a methodology for comparing evolutionary and swarm intelligence algorithms. Swarm Evolut. Comput. 1(1), 3–18 (2011) 34. Dolan, E.D., Mor´ e, J.J.: Benchmarking optimization software with performance profiles. Math. Program. Ser. A, 91(2), 201–213 (2002) 35. Deb, K., Srivastava, S.: A genetic algorithm based augmented Lagrangian method for constrained optimization. Comput. Optim. Appl. 53(3), 869–902 (2012) 36. Rocha, A.M.A.C., Costa, M.F.P., Fernandes, E.M.G.P.: A shifted hyperbolic augmented Lagrangianbased artificial fish two-swarm algorithm with guaranteed convergence for constrained global optimization. Eng. Optimiz. 48(12), 2114–2140 (2016) 37. Hooke, R., Jeeves, T.A.: Direct search solution of numerical and statistical problems. J. Ass. Comput. Mach. 8(2), 212–229 (1961)