# Heuristic-Based Firefly Algorithm for Bound Constrained Nonlinear Binary Optimization

Authors: M. Fernanda P. Costa; Ana Maria A. C. Rocha; Rogério B. Francisco; Edite M. G. P. Fernandes
Journal: Advances in Operations Research
 (2014)

---

## Abstract

Firefly algorithm (FA) is a metaheuristic for global optimization. In this paper, we address the practical testing of a heuristic-based FA (HBFA) for computing optima of discrete nonlinear optimization problems, where the discrete variables are of binary type. An important issue in FA is the formulation of attractiveness of each firefly which in turn affects its movement in the search space. Dynamic updating schemes are proposed for two parameters, one from the attractiveness term and the other from the randomization term. Three simple heuristics capable of transforming real continuous variables into binary ones are analyzed. A new sigmoid “erf” function is proposed. In the context of FA, three different implementations to incorporate the heuristics for binary variables into the algorithm are proposed. Based on a set of benchmark problems, a comparison is carried out with other binary dealing metaheuristics. The results demonstrate that the proposed HBFA is efficient and outperforms binary versions of differential evolution (DE) and particle swarm optimization (PSO). The HBFA also compares very favorably with angle modulated version of DE and PSO. It is shown that the variant of HBFA based on the sigmoid “erf” function with “movements in continuous space” is the best, in terms of both computational requirements and accuracy.


---

## Body

## 1. Introduction

This paper aims to analyze the merit, in terms of performance, of a heuristic-based firefly algorithm (HBFA) for computing the optimal and binary solution of bound constrained nonlinear optimization problems. The problem to be addressed has the form(1)min⁡ffffifff(x)subjecttox∈Ω⊂Rn(acompactconvexset)ffffffffffixl∈{0,1}forl=1,…,n,
where f is a continuous function. Due to the compactness of Ω, we also have Lbl≤xl≤Ubl, l=1,…,n, where Lb and Ub are the vectors of the lower and upper bounds, respectively. We do not assume that f is differentiable and convex. Instead of searching for any local (nonglobal) solution we want the globally best binary point. Direct search methods might be suitable since we do not assume differentiability. However, they are only local optimization procedures and therefore there is no guarantee that a global solution is reached. For global optimization, stochastic methods are generally used and aim to explore the search space and converge to a global solution. Metaheuristics are higher-level procedures or heuristics that are designed to search for good solutions, known as near-optimal solutions, with less computational effort and time than more classical algorithms. They are usually nondeterministic and their behaviors do not depend on problem’s properties. Population-based metaheuristics have been used to solve a variety of optimization problems, from continuous to the combinatorial ones.Metaheuristics are common for solving discrete binary optimization problems [1–10]. Many approaches have been developed aiming to solve nonlinear programming problems with mixed-discrete variables by transforming the discrete problem into a continuous one [11]. The most used and simple approach solves the continuous relaxed problem and then discretizes the obtained solution by using a rounding scheme. This type of approach works well on simple and small dimension academic and benchmark problems but may be somehow limited on some real-world applications.Recently, a metaheuristic optimization algorithm, termed firefly algorithm (FA), that mimics the social behavior of fireflies based on the flashing and attraction characteristics of fireflies, has been developed [12, 13]. This is a swarm intelligence optimization algorithm that is capable of competing with the most well-known algorithms, like ant colony optimization, particle swarm optimization, artificial bee colony, artificial fish swarm, and cuckoo-search.FA performance is controlled by three parameters: the randomization parameterα, the attractiveness β, and the absorption coefficient γ. Authors have argued that its efficiency is due to its capability of subdividing the population into subgroups (since local attraction is stronger than long-distance attraction) and its ability to adapt the search to problem landscape by controlling the parameter γ [14, 15]. Several variants of the firefly algorithm do already exist in the literature. Based on the settings of their parameters, a classification scheme has appeared. Gaussian FA [16], hybrid FA with harmony search [17], hybrid genetic algorithm with FA [18], self-adaptive step FA [15], and modified FA in [19] are just a few examples. Further improvements have been made aiming to accelerate convergence (see, e.g., [20–22]). A practical convergence analysis of FA with different parameter sets is presented in [23]. FA has become popular and widely used in recent years in many applications, like economic dispatch problems [24] and mixed variable optimization problems [25]. The extension of FA to multiobjective continuous optimization has already been investigated [26]. A recent review of firefly algorithms is available in [14].Based on the effectiveness of FA in continuous optimization, it is predicted that it will perform well when solving discrete optimization problems. Discrete versions of the FA are available for solving discrete NP hard optimization problems [27, 28].The main purpose of this study is to incorporate some heuristics aiming to deal with binary variables in the firefly algorithm for solving nonlinear optimization problems with binary variables. The binary dealing methods that were implemented are adaptations of well-known heuristics for defining 0 and 1 bit strings from real variables. Furthermore, a new sigmoid function aiming to constrain a real valued variable to the range[0,1] is also proposed. Three different implementations to incorporate the heuristics for binary variables and adapt FA to binary optimization are proposed. We apply the proposed heuristic strategies to solve a set of benchmark nonlinear problems and show that the newly developed HBFA is effective in binary nonlinear programming.The remaining part of the paper is organized as follows. Section2 reviews the standard FA and presents new dynamic updates for some FA parameters, and Section 3 describes different heuristic strategies and reports on their implementations to adapt FA to binary optimization. All the heuristic approaches are validated with tests on a set of well-known bound constrained problems. These results and a comparison with other methods in the literature are shown in Section 4. Finally, the conclusions and ideas for future work are listed in Section 5.


## 2. Firefly Algorithm

Firefly algorithm is a bioinspired metaheuristic algorithm that is able to compute a solution to an optimization problem. It is inspired by the flashing behavior of fireflies at night. According to [12, 13, 19], the three main rules used to construct the standard algorithm are the following:(i)
all fireflies are unisex, meaning that any firefly can be attracted to any other brighter one;(ii)
the attractiveness of a firefly is determined by its brightness which is associated with the encoded objective function;(iii)
attractiveness is directly proportional to brightness but decreases with distance.Throughout this paper, we let∥·∥ represent the Euclidean norm of a vector. We use the vector x=(x1,x2,…,xn)T to represent the position of a firefly in the search space. The position of the firefly j will be represented by xj∈Rn. We assume that the size of the population of fireflies is m. In the context of problem (1), firefly j is brighter than firefly i if f(xj)<f(xi).


### 2.1. Standard FA

First, in the standard FA, the positions of the fireflies are randomly generated in the search spaceΩ as follows:
(2)xli=Lbl+rand(Ubl-Lbl),forl=1,…,n,
where rand is a uniformly distributed random number in [0,1], hereafter represented by rand~U[0,1]. The movement of a firefly i is attracted to another brighter firefly j and is given by
(3)xi=xi+β(xj-xi)+α(rand-0.5)S,
where α∈[0,1] is the randomization parameter, rand~U[0,1], S∈Rn is a problem dependent vector of scaling parameters, and
(4)β=β0exp⁡(-γ∥xi-xj∥p)forp≥1
gives the attractiveness of a firefly which varies with the light intensity/brightness seen by adjacent fireflies and the distance between themselves and β0 is the attraction parameter when the distance is zero [12, 13, 22, 29]. Besides the presented “exp⁡” function, any monotonically decreasing function could be used. The parameter γ which characterizes the variation of the attractiveness is the light absorption coefficient and is crucial to determine the speed of convergence of the algorithm. In theory, γ could take any value in the set [0,∞). When γ→0, the attractiveness is constant β=β0, meaning that a flashing firefly can be seen anywhere in the search space. This is an ideal case for a problem with a single (usually global) optimum since it can easily be reached. On the other hand, when γ→∞, the attractiveness is almost zero in the sight of other fireflies and each firefly moves in a random way. In particular, when β0=0, the algorithm behaves like a random search method [13, 22]. The randomization term can be extended to the normal distribution N(0,1) or to any other distribution [15]. Algorithm 1 presents the main steps of the standard FA (on continuous space).Algorithm 1:Standard FA.
Data:  kmax⁡, f*, η
Setk=0;
Randomly generatexi∈Ω, i=1,…,m;
Evaluatef(xi),i=1,…,m, rank fireflies (from lowest to largest f);
while  k≤kmax⁡  and|f(x1)-f*|>η  do
forall
x
i  such thati=2,…,m  do
forall
x
j  such thatj=1,…,i-1  do
Compute randomization term;
Compute attractiveness  β;
Move firefly i towards j using (3);
Evaluate f(xi),i=1,…,m, rank fireflies (from lowest to largest f);
Set  k=k+1;


### 2.2. Dynamic Updates ofα and γ

The relative value of the parametersα and γ affects the performance of FA. The parameter α controls the randomness or, to some extent, the diversity of solutions. Parameter γ aims to scale the attraction power of the algorithm. Small values of γ with large values of α can increase the number of iterations required to converge to an optimal solution. Experience has shown that α must take large values at the beginning of the iterative process to enforce the algorithm to increase the diversity of solutions. However, small α values combined with small values of γ in the final iterations increase the fine-tuning of solutions since the effort is focused on exploitation. Thus, it is possible to improve the quality of the solution by reducing the randomness. Convergence can be improved by varying the randomization parameter α so that it decreases gradually as the optimum solution is approaching [22, 24, 26, 29]. In order to improve convergence speed and solution accuracy, dynamic updates of the parameters α and γ of FA, which depend on the iteration counter k of the algorithm, are implemented as follows.Similarly to the factor which controls the amplification of differential variations, in differential evolution (DE) metaheuristic [5], the inertial weight, in particle swarm optimization (PSO) [29, 30], and the pitch adjusting rate, in the harmony search (HS) algorithm [31], we allow the value of α to decrease linearly with k, from an upper level αmax⁡ to a lower level αmin⁡:
(5)α(k)=αmax⁡-kαmax⁡-αmin⁡kmax⁡,
where kmax⁡ is the maximum number of allowed iterations. To increase the attractiveness with k, the parameter γ is dynamically updated by
(6)γ(k)=γmax⁡exp⁡(kkmax⁡ln⁡(γmin⁡γmax⁡)),
where γmin⁡ and γmax⁡ are the minimum variation and maximum variation of attractiveness, respectively.


### 2.3. Lévy Dependent Randomization Term

We remark that our implementation of the randomization term in the proposed dynamic FA considers the Lévy distribution. Based on the attractivenessβ, in (4), the equation for the movement of firefly i towards a brighter firefly j can be written as follows:
(7)xi=xi+yiwithyi=β(xj-xi)+αL(x1)σxi,
where L(x1) is a random number from the Lévy distribution centered at x1, the position of the brightest firefly, with an unitary standard deviation. The vector σxi represents the variation around x1 (and based on real position x)
(8)σxi=(|x1i-x11|,…,|xni-xn1|)T.


## 2.1. Standard FA

First, in the standard FA, the positions of the fireflies are randomly generated in the search spaceΩ as follows:
(2)xli=Lbl+rand(Ubl-Lbl),forl=1,…,n,
where rand is a uniformly distributed random number in [0,1], hereafter represented by rand~U[0,1]. The movement of a firefly i is attracted to another brighter firefly j and is given by
(3)xi=xi+β(xj-xi)+α(rand-0.5)S,
where α∈[0,1] is the randomization parameter, rand~U[0,1], S∈Rn is a problem dependent vector of scaling parameters, and
(4)β=β0exp⁡(-γ∥xi-xj∥p)forp≥1
gives the attractiveness of a firefly which varies with the light intensity/brightness seen by adjacent fireflies and the distance between themselves and β0 is the attraction parameter when the distance is zero [12, 13, 22, 29]. Besides the presented “exp⁡” function, any monotonically decreasing function could be used. The parameter γ which characterizes the variation of the attractiveness is the light absorption coefficient and is crucial to determine the speed of convergence of the algorithm. In theory, γ could take any value in the set [0,∞). When γ→0, the attractiveness is constant β=β0, meaning that a flashing firefly can be seen anywhere in the search space. This is an ideal case for a problem with a single (usually global) optimum since it can easily be reached. On the other hand, when γ→∞, the attractiveness is almost zero in the sight of other fireflies and each firefly moves in a random way. In particular, when β0=0, the algorithm behaves like a random search method [13, 22]. The randomization term can be extended to the normal distribution N(0,1) or to any other distribution [15]. Algorithm 1 presents the main steps of the standard FA (on continuous space).Algorithm 1:Standard FA.
Data:  kmax⁡, f*, η
Setk=0;
Randomly generatexi∈Ω, i=1,…,m;
Evaluatef(xi),i=1,…,m, rank fireflies (from lowest to largest f);
while  k≤kmax⁡  and|f(x1)-f*|>η  do
forall
x
i  such thati=2,…,m  do
forall
x
j  such thatj=1,…,i-1  do
Compute randomization term;
Compute attractiveness  β;
Move firefly i towards j using (3);
Evaluate f(xi),i=1,…,m, rank fireflies (from lowest to largest f);
Set  k=k+1;


## 2.2. Dynamic Updates ofα and γ

The relative value of the parametersα and γ affects the performance of FA. The parameter α controls the randomness or, to some extent, the diversity of solutions. Parameter γ aims to scale the attraction power of the algorithm. Small values of γ with large values of α can increase the number of iterations required to converge to an optimal solution. Experience has shown that α must take large values at the beginning of the iterative process to enforce the algorithm to increase the diversity of solutions. However, small α values combined with small values of γ in the final iterations increase the fine-tuning of solutions since the effort is focused on exploitation. Thus, it is possible to improve the quality of the solution by reducing the randomness. Convergence can be improved by varying the randomization parameter α so that it decreases gradually as the optimum solution is approaching [22, 24, 26, 29]. In order to improve convergence speed and solution accuracy, dynamic updates of the parameters α and γ of FA, which depend on the iteration counter k of the algorithm, are implemented as follows.Similarly to the factor which controls the amplification of differential variations, in differential evolution (DE) metaheuristic [5], the inertial weight, in particle swarm optimization (PSO) [29, 30], and the pitch adjusting rate, in the harmony search (HS) algorithm [31], we allow the value of α to decrease linearly with k, from an upper level αmax⁡ to a lower level αmin⁡:
(5)α(k)=αmax⁡-kαmax⁡-αmin⁡kmax⁡,
where kmax⁡ is the maximum number of allowed iterations. To increase the attractiveness with k, the parameter γ is dynamically updated by
(6)γ(k)=γmax⁡exp⁡(kkmax⁡ln⁡(γmin⁡γmax⁡)),
where γmin⁡ and γmax⁡ are the minimum variation and maximum variation of attractiveness, respectively.


## 2.3. Lévy Dependent Randomization Term

We remark that our implementation of the randomization term in the proposed dynamic FA considers the Lévy distribution. Based on the attractivenessβ, in (4), the equation for the movement of firefly i towards a brighter firefly j can be written as follows:
(7)xi=xi+yiwithyi=β(xj-xi)+αL(x1)σxi,
where L(x1) is a random number from the Lévy distribution centered at x1, the position of the brightest firefly, with an unitary standard deviation. The vector σxi represents the variation around x1 (and based on real position x)
(8)σxi=(|x1i-x11|,…,|xni-xn1|)T.


## 3. Dealing with Binary Variables

The standard FA is a real-coded algorithm and some modifications are needed to enable it to deal with discrete optimization problems. This section describes the implementation of some heuristics with FA for binary nonlinear optimization problems. In the context of the proposed HBFA, three different heuristics to transform a continuous real variable into a binary one are presented. Furthermore, to extend FA to binary optimization, different implementations to incorporate the heuristic strategies into FA are described. We will use the term “discretization” to define the process that transforms a continuous real variable, represented, for example, byx, into a binary one, represented by b.


### 3.1. Sigmoid Logistic Function

This discretization methodology is the most common in the literature when population-based stochastic algorithms are considered in binary optimization, namely, PSO [6, 8, 9], DE [3], HS [1, 32], artificial fish swarm [33], and artificial bee colony [4, 7, 10].Whenxi moves towards xj, the likelihood is that the discrete components of xi change from binary numbers to real ones. To transform a real number into a binary one, the following sigmoid logistic function constrains the real value to the interval [0,1]:
(9)sig(xli)=11+exp⁡(-xli),
where xli, in the context of FA, is the component l of the position vector xi (of firefly i) after movement—recall (7) and (4). Equation (9) interprets the floating-point components of a solution as a set of probabilities. These are then used to assign appropriate binary values by using
(10)bli={1,ifrand≤sig(xli)0,otherwise,
where sig(xli) gives the probability that the component itself is 0 or 1 [28] and rand~U[0,1]. We note that during the iterative process the firefly positions, x, were not allowed to move outside the search space Ω.


### 3.2. Proposed Sigmoiderf⁡ Function

The error function is a special function with a shape that appears mainly in probability and statistics contexts. Denoted by “erf,” the mathematical function defined by the integral,(11)erf⁡(x)=2π∫0x‍exp⁡(-t2)dt,
satisfies the following properties
(12)erf⁡(0)=0,erf⁡(-∞)=-1,erf⁡(+∞)=1,erf⁡(-x)=-erf⁡(x)
and it has a close relation with the normal distribution probabilities. When a series of measurements are described by a normal distribution with mean 0 and standard deviation σ, the erf function evaluated at (x/σ2), for a positive x, gives the probability that the error of a single measurement lies in the interval [-x,x]. The derivative of the erf⁡ function follows immediately from its definition:
(13)ddterf⁡(t)=2πexp⁡(-t2),fort∈R.
The good properties of the erf⁡ function are thus used to define a new sigmoid function, the sigmoid erf⁡ function:
(14)sig(xli)=0.5(1+erf⁡(xli)),
which is a bounded differentiable real function defined for all x∈R and has a positive derivative at each point. A comparison of both functions (9) and (14) is depicted in Figure 1. Note that the slope at the origin of the sigmoid function in (14) is around 0.5641895, while that of function (9) is 0.25, thus yielding a faster growing from 0 to 1.Figure 1
Sigmoid functions.


### 3.3. Rounding to Integer Part

The simplest discretization procedure of a continuous component of a point into 0/1 bit uses the rounding to the integer part function, known as floor function, and is described in [34]. Each continuous value xli∈R is transformed into a binary one, 0 bit or 1 bit, bli, for l=1,…,n in the following way:
(15)bli=⌊|xlimod⁡2|⌋,
where ⌊z⌋ represents the floor function of z and gives the largest integer not greater than z. The floating-point value xli is first divided by 2 and then the absolute value of the remainder is floored. The obtained integer number is the bit value of the component.


### 3.4. Heuristics Implementation

In this study, three methods capable of computing global solutions to binary optimization problems using FA are proposed.


#### 3.4.1. Movement on Continuous Space

In this implementation of the previously described heuristics, denoted by “movement on continuous space” (mCS), the movement of each firefly is made on the continuous space and its attractiveness term is updated considering the real position vector. The real position of fireflyi is discretized only after all movements towards brighter fireflies have been carried out. We note that the fitness evaluation of each firefly, for firefly ranking, is always based on the binary position. Algorithm 2 presents the main steps of HBFA with mCS.Algorithm 2:HBFA with mCS.
Data:  kmax⁡, f*, η
Setk=0;
Randomly generatexi∈Ω, i=1,…,m;
Discretize position of fireflyi: xi→bi, i=1,…,m;
Computef(bi),i=1,…,m, rank fireflies (from lowest to largest  f);
while
k
≤
k
max
⁡  and|f(b1)-f*|>η  do
forall  xi  such thati=2,…,m  do
forall    xj  such thatj=1,…,i-1  do
Compute randomization term;
Compute attractiveness β;
Move position xi of firefly i towards xj using (7);
Discretize positions: xi→bi, i=1,…,m;
Compute f(bi),i=1,…,m, rank fireflies (from lowest to largest f);
Set k=k+1;


#### 3.4.2. Movement on Binary Space

This implementation, denoted by “movement on binary space” (mBS), moves the binary position of each firefly towards the binary positions of brighter fireflies; that is, each movement is made on the binary space although the corresponding position may fail to be 0 or 1 bit string and must be discretized before the updating of attractiveness. Here, fitness is also based on the binary positions. Algorithm3 presents the main steps of HBFA with mBS.Algorithm 3:HBFA with mBS.
Data: kmax⁡, f*, η
Setk=0;
Randomly generatexi∈Ω, i=1,…,m;
Discretize position of fireflyi: xi→bi, i=1,…,m;
Computef(bi),i=1,…,m, rank fireflies (from lowest to largest f);
while  k≤kmax⁡  and|f(b1)-f*|>η  do
forall  bi  such thati=2,…,m  do
forall  bj  such thatj=1,…,i-1  do
Compute randomization term;
Compute attractiveness β based on distance ∥bi-bj∥p;
Move binary position bi of firefly i towards  bj
using xi=bi+β(bj-bi)+αL(b1)σbi;
Discretize position of firefly i: xi→bi;
Compute f(bi),i=1,…,m, rank fireflies (from lowest to largest  f);
Set k=k+1;


#### 3.4.3. Probability for Binary Component

For this implementation, named “probability for binary component” (pBC), we borrow the concept from the binary PSO [6, 9, 35] where each component of the velocity vector is directly used to compute the probability that the corresponding component of the particle position, xli, is 0 or 1. Similarly, in the FA algorithm, we do not interpret the vector yi in (7) as a step size, but rather as a mean to compute the probability that each component of the position vector of firefly i is 0 or 1. Thus, we define
(16)bli={1,ifrand≤sig(yli)0,otherwise,
where sig() represents a sigmoid function. Algorithm 4 is the pseudocode of HBFA with pBC.Algorithm 4:HBFA with pBC.
Data:  kmax⁡, f*, η
Setk=0;
Randomly generatexi∈Ω, i=1,…,m;
Discretize position of fireflyi: xi→bi, i=1,…,m;
Computef(bi),i=1,…,m, rank fireflies (from lowest to largest f);
while  k≤kmax⁡  and|f(b1)-f*|>η  do
forall
b
i  such thati=2,…,m  do
forall  bj  such thatj=1,…,i-1  do
Compute randomization term;
Compute attractiveness β based on distance ∥bi-bj∥p;
Compute yi using binary positions (see (7));
Discretize yi and define bi using (16);
Compute f(bi),i=1,…,m, rank fireflies (from lowest to largest  f);
Setk=k+1;


## 3.1. Sigmoid Logistic Function

This discretization methodology is the most common in the literature when population-based stochastic algorithms are considered in binary optimization, namely, PSO [6, 8, 9], DE [3], HS [1, 32], artificial fish swarm [33], and artificial bee colony [4, 7, 10].Whenxi moves towards xj, the likelihood is that the discrete components of xi change from binary numbers to real ones. To transform a real number into a binary one, the following sigmoid logistic function constrains the real value to the interval [0,1]:
(9)sig(xli)=11+exp⁡(-xli),
where xli, in the context of FA, is the component l of the position vector xi (of firefly i) after movement—recall (7) and (4). Equation (9) interprets the floating-point components of a solution as a set of probabilities. These are then used to assign appropriate binary values by using
(10)bli={1,ifrand≤sig(xli)0,otherwise,
where sig(xli) gives the probability that the component itself is 0 or 1 [28] and rand~U[0,1]. We note that during the iterative process the firefly positions, x, were not allowed to move outside the search space Ω.


## 3.2. Proposed Sigmoiderf⁡ Function

The error function is a special function with a shape that appears mainly in probability and statistics contexts. Denoted by “erf,” the mathematical function defined by the integral,(11)erf⁡(x)=2π∫0x‍exp⁡(-t2)dt,
satisfies the following properties
(12)erf⁡(0)=0,erf⁡(-∞)=-1,erf⁡(+∞)=1,erf⁡(-x)=-erf⁡(x)
and it has a close relation with the normal distribution probabilities. When a series of measurements are described by a normal distribution with mean 0 and standard deviation σ, the erf function evaluated at (x/σ2), for a positive x, gives the probability that the error of a single measurement lies in the interval [-x,x]. The derivative of the erf⁡ function follows immediately from its definition:
(13)ddterf⁡(t)=2πexp⁡(-t2),fort∈R.
The good properties of the erf⁡ function are thus used to define a new sigmoid function, the sigmoid erf⁡ function:
(14)sig(xli)=0.5(1+erf⁡(xli)),
which is a bounded differentiable real function defined for all x∈R and has a positive derivative at each point. A comparison of both functions (9) and (14) is depicted in Figure 1. Note that the slope at the origin of the sigmoid function in (14) is around 0.5641895, while that of function (9) is 0.25, thus yielding a faster growing from 0 to 1.Figure 1
Sigmoid functions.


## 3.3. Rounding to Integer Part

The simplest discretization procedure of a continuous component of a point into 0/1 bit uses the rounding to the integer part function, known as floor function, and is described in [34]. Each continuous value xli∈R is transformed into a binary one, 0 bit or 1 bit, bli, for l=1,…,n in the following way:
(15)bli=⌊|xlimod⁡2|⌋,
where ⌊z⌋ represents the floor function of z and gives the largest integer not greater than z. The floating-point value xli is first divided by 2 and then the absolute value of the remainder is floored. The obtained integer number is the bit value of the component.


## 3.4. Heuristics Implementation

In this study, three methods capable of computing global solutions to binary optimization problems using FA are proposed.


### 3.4.1. Movement on Continuous Space

In this implementation of the previously described heuristics, denoted by “movement on continuous space” (mCS), the movement of each firefly is made on the continuous space and its attractiveness term is updated considering the real position vector. The real position of fireflyi is discretized only after all movements towards brighter fireflies have been carried out. We note that the fitness evaluation of each firefly, for firefly ranking, is always based on the binary position. Algorithm 2 presents the main steps of HBFA with mCS.Algorithm 2:HBFA with mCS.
Data:  kmax⁡, f*, η
Setk=0;
Randomly generatexi∈Ω, i=1,…,m;
Discretize position of fireflyi: xi→bi, i=1,…,m;
Computef(bi),i=1,…,m, rank fireflies (from lowest to largest  f);
while
k
≤
k
max
⁡  and|f(b1)-f*|>η  do
forall  xi  such thati=2,…,m  do
forall    xj  such thatj=1,…,i-1  do
Compute randomization term;
Compute attractiveness β;
Move position xi of firefly i towards xj using (7);
Discretize positions: xi→bi, i=1,…,m;
Compute f(bi),i=1,…,m, rank fireflies (from lowest to largest f);
Set k=k+1;


### 3.4.2. Movement on Binary Space

This implementation, denoted by “movement on binary space” (mBS), moves the binary position of each firefly towards the binary positions of brighter fireflies; that is, each movement is made on the binary space although the corresponding position may fail to be 0 or 1 bit string and must be discretized before the updating of attractiveness. Here, fitness is also based on the binary positions. Algorithm3 presents the main steps of HBFA with mBS.Algorithm 3:HBFA with mBS.
Data: kmax⁡, f*, η
Setk=0;
Randomly generatexi∈Ω, i=1,…,m;
Discretize position of fireflyi: xi→bi, i=1,…,m;
Computef(bi),i=1,…,m, rank fireflies (from lowest to largest f);
while  k≤kmax⁡  and|f(b1)-f*|>η  do
forall  bi  such thati=2,…,m  do
forall  bj  such thatj=1,…,i-1  do
Compute randomization term;
Compute attractiveness β based on distance ∥bi-bj∥p;
Move binary position bi of firefly i towards  bj
using xi=bi+β(bj-bi)+αL(b1)σbi;
Discretize position of firefly i: xi→bi;
Compute f(bi),i=1,…,m, rank fireflies (from lowest to largest  f);
Set k=k+1;


### 3.4.3. Probability for Binary Component

For this implementation, named “probability for binary component” (pBC), we borrow the concept from the binary PSO [6, 9, 35] where each component of the velocity vector is directly used to compute the probability that the corresponding component of the particle position, xli, is 0 or 1. Similarly, in the FA algorithm, we do not interpret the vector yi in (7) as a step size, but rather as a mean to compute the probability that each component of the position vector of firefly i is 0 or 1. Thus, we define
(16)bli={1,ifrand≤sig(yli)0,otherwise,
where sig() represents a sigmoid function. Algorithm 4 is the pseudocode of HBFA with pBC.Algorithm 4:HBFA with pBC.
Data:  kmax⁡, f*, η
Setk=0;
Randomly generatexi∈Ω, i=1,…,m;
Discretize position of fireflyi: xi→bi, i=1,…,m;
Computef(bi),i=1,…,m, rank fireflies (from lowest to largest f);
while  k≤kmax⁡  and|f(b1)-f*|>η  do
forall
b
i  such thati=2,…,m  do
forall  bj  such thatj=1,…,i-1  do
Compute randomization term;
Compute attractiveness β based on distance ∥bi-bj∥p;
Compute yi using binary positions (see (7));
Discretize yi and define bi using (16);
Compute f(bi),i=1,…,m, rank fireflies (from lowest to largest  f);
Setk=k+1;


## 3.4.1. Movement on Continuous Space

In this implementation of the previously described heuristics, denoted by “movement on continuous space” (mCS), the movement of each firefly is made on the continuous space and its attractiveness term is updated considering the real position vector. The real position of fireflyi is discretized only after all movements towards brighter fireflies have been carried out. We note that the fitness evaluation of each firefly, for firefly ranking, is always based on the binary position. Algorithm 2 presents the main steps of HBFA with mCS.Algorithm 2:HBFA with mCS.
Data:  kmax⁡, f*, η
Setk=0;
Randomly generatexi∈Ω, i=1,…,m;
Discretize position of fireflyi: xi→bi, i=1,…,m;
Computef(bi),i=1,…,m, rank fireflies (from lowest to largest  f);
while
k
≤
k
max
⁡  and|f(b1)-f*|>η  do
forall  xi  such thati=2,…,m  do
forall    xj  such thatj=1,…,i-1  do
Compute randomization term;
Compute attractiveness β;
Move position xi of firefly i towards xj using (7);
Discretize positions: xi→bi, i=1,…,m;
Compute f(bi),i=1,…,m, rank fireflies (from lowest to largest f);
Set k=k+1;


## 3.4.2. Movement on Binary Space

This implementation, denoted by “movement on binary space” (mBS), moves the binary position of each firefly towards the binary positions of brighter fireflies; that is, each movement is made on the binary space although the corresponding position may fail to be 0 or 1 bit string and must be discretized before the updating of attractiveness. Here, fitness is also based on the binary positions. Algorithm3 presents the main steps of HBFA with mBS.Algorithm 3:HBFA with mBS.
Data: kmax⁡, f*, η
Setk=0;
Randomly generatexi∈Ω, i=1,…,m;
Discretize position of fireflyi: xi→bi, i=1,…,m;
Computef(bi),i=1,…,m, rank fireflies (from lowest to largest f);
while  k≤kmax⁡  and|f(b1)-f*|>η  do
forall  bi  such thati=2,…,m  do
forall  bj  such thatj=1,…,i-1  do
Compute randomization term;
Compute attractiveness β based on distance ∥bi-bj∥p;
Move binary position bi of firefly i towards  bj
using xi=bi+β(bj-bi)+αL(b1)σbi;
Discretize position of firefly i: xi→bi;
Compute f(bi),i=1,…,m, rank fireflies (from lowest to largest  f);
Set k=k+1;


## 3.4.3. Probability for Binary Component

For this implementation, named “probability for binary component” (pBC), we borrow the concept from the binary PSO [6, 9, 35] where each component of the velocity vector is directly used to compute the probability that the corresponding component of the particle position, xli, is 0 or 1. Similarly, in the FA algorithm, we do not interpret the vector yi in (7) as a step size, but rather as a mean to compute the probability that each component of the position vector of firefly i is 0 or 1. Thus, we define
(16)bli={1,ifrand≤sig(yli)0,otherwise,
where sig() represents a sigmoid function. Algorithm 4 is the pseudocode of HBFA with pBC.Algorithm 4:HBFA with pBC.
Data:  kmax⁡, f*, η
Setk=0;
Randomly generatexi∈Ω, i=1,…,m;
Discretize position of fireflyi: xi→bi, i=1,…,m;
Computef(bi),i=1,…,m, rank fireflies (from lowest to largest f);
while  k≤kmax⁡  and|f(b1)-f*|>η  do
forall
b
i  such thati=2,…,m  do
forall  bj  such thatj=1,…,i-1  do
Compute randomization term;
Compute attractiveness β based on distance ∥bi-bj∥p;
Compute yi using binary positions (see (7));
Discretize yi and define bi using (16);
Compute f(bi),i=1,…,m, rank fireflies (from lowest to largest  f);
Setk=k+1;


## 4. Numerical Experiments

In this section, we present the computational results that were obtained with HBFA—Algorithms2, 3, and 4, using (9), (14), or (15)—aiming to investigate its performance when solving a set of binary nonlinear optimization problems. Two small 0-1 knapsack problems are also used to test the algorithms’ behavior on linear problems with 0/1 variables.The numerical experiments were carried out on a PC Intel Core 2 Duo Processor E7500 with 2.9 GHz and 4 Gb of memory. The algorithms were coded in Matlab Version 8.0.0.783 (R2012b).


### 4.1. Experimental Setting

Each experiment was conducted 30 times. The size of the population is made to depend on the problem’s dimension and is set tom=min⁡{40,2n}. Some experiments have been carried out to tune certain parameters of the algorithms. In the proposed FA with dynamic α and γ, they are set as follows: β0=1, p=1, αmax⁡=0.5, αmin⁡=0.01, γmax⁡=10, and γmin⁡=0.1. In Algorithms 2 (mCS), 3 (mBS), and 4 (pBC), iterations were limited to kmax⁡=500 and the tolerance for finding a good quality solution is η=10-6. Results reported are averaged (over the 30 runs) of best function values, number of function evaluations, and number of iterations.


### 4.2. Experimental Results

First, we use a set of ten benchmark nonlinear functions with different dimensions and characteristics. For example, five functions are unimodal and the remaining multimodal [3, 9, 10, 36]. They are displayed in Table 1. Although they are widely used in continuous optimization, we now aim to converge to a 0/1 bit string solution.Table 1
Problems set.
Ackley
f
(
x
)
=
-
20
exp
⁡
(
-
0.2
1
n
∑
i
=
1
n
‍
x
i
2
)
-
exp
⁡
(
1
n
∑
i
=
1
n
‍
cos
⁡
(
2
π
x
i
)
)
+
20
+
e
n = 30, Ω=[-30,30]30, f*=0  at  x*=(0,…,0)
Foxholes
f
(
x
)
=
1
l
0.002
+
∑
j
=
1
25
‍
(
1
/
(
j
+
(
x
1
-
a
1
j
)
6
+
(
x
2
-
a
2
j
)
6
)
)
[
a
i
j
]
=
[
-
32
-
16
0
16
32
-
32
-
16
⋯
16
32
-
32
-
32
-
32
-
32
-
32
-
16
-
16
⋯
32
32
]
n = 2, Ω=[-65.536,65.536]2, f*≈13 at x*=(0,0)
Griewank
f
(
x
)
=
1
+
1
4000
∑
i
=
1
n
‍
x
i
2
-
∏
i
=
1
n
cos
⁡
(
x
i
i
)
n = 30, Ω=[-300,300]30, f*=0  at  x*=(0,…,0)
Quartic
f
(
x
)
=
∑
i
=
1
n
‍
i
x
i
4
+
U
[
0,1
]
n = 30, Ω=[-1.28,1.28]30, f*=0+noise at x*=(0,…,0)
Rastrigin
f
(
x
)
=
10
n
+
∑
i
=
1
n
‍
(
x
i
2
-
10
cos
⁡
(
2
π
x
i
)
)
n = 30, Ω=[-5.12,5.12]30, f*=0  at  x*=(0,…,0)
Rosenbrock2
f
(
x
)
=
100
(
x
1
2
-
x
2
)
2
+
(
1
-
x
1
)
2
n = 2, Ω=[-2.048,2.048]2, f*=0  at  x*=(1,1)
Rosenbrock
f
(
x
)
=
∑
i
=
1
n
-
1
‍
100
(
x
i
2
-
x
i
+
1
)
2
+
(
1
-
x
i
)
2
n = 30, Ω=[-2.048,2.048]30, f*=0  at  x*=(1,…,1)
Schaffer
f
(
x
)
=
0.5
+
(
sin
⁡
⁡
(
x
1
2
+
x
2
2
)
)
2
-
0.5
(
1
+
0.001
(
x
1
2
+
x
2
2
)
)
2
n = 2, Ω=[-100,100]2, f*=0 at x*=(0,0)
Spherical
f
(
x
)
=
∑
i
=
1
n
‍
x
i
2
n = 3, Ω=[-5.12,5.12]3, f*=0 at x*=(0,0,0)
Step
f
(
x
)
=
6
n
+
∑
i
=
1
n
‍
⌊
x
i
⌋
n = 5, Ω=[-5.12,5.12]5, f*=30 at  x*=(0,0,0,0,0)First, we aim to compare with the results reported in [3, 9, 10]. Due to poor results, the authors in [10] do not recommend the use of ABC to solve binary-valued problems. The other metaheuristics therein implemented are the following:(i)
angle modulated PSO (AMPSO) and angle modulated DE (AMDE) that incorporate a trigonometric function as a bit string generator into the classic PSO and DE algorithms, respectively;(ii)
binary DE and PSO based on the sigmoid logistic function and (10), denoted by binDE and binPSO, respectively.We noticed that the problems Foxholes, Griewank, Rosenbrock, Schaffer, and Step are not correctly described in [3, 9, 10]. Table 2 shows both the averaged best function values (obtained during the 30 runs), favg, with the St.D. in parentheses, and the averaged number of function evaluations, nfe, obtained with the sigmoid logistic function (see in (9)) and (10), while using the three implementations: mCS, mBS, and pBC. Results obtained for these ten functions indicate that our proposal HBFA produces high quality solutions and outperforms the binary versions binPSO and binDE, as well as AMPSO and AMDE. We also note that mCS has the best “nfe” values on 6 problems, mBS is better on 3 problems (one is a tie with mCS), and pBC on 2 problems. Thus, the performance of mCS is the best when compared with those of mBS and pBC. The latter is the least efficient of all, in particular for the large dimensional problems.Table 2
Comparison with AMPSO, binPSO, binDE, and AMDE based onfavg and St.D. (shown in parentheses).
Prob.
mCS based on (9)
mBS based on (9)
pBC based on (9)
AMPSO
binPSO
binDE
AMDE
f
avg
nfe
f
avg
nfe
f
avg
nfe
f
avg
f
avg
f
avg
f
avg
Ackley
8.88
e
-
16
80
8.88
e
-
16
1156
8.88
e
-
16
2168
1.97
e
01
2.01
e
01
1.73
e
01
1.64
e
01
(
0.00
e
00
)
(
0.00
e
00
)
(
0.00
e
00
)
(
0.57
e
-
01
)
(
0.49
e
-
01
)
(
0.31
e
01
)
(
0.76
e
00
)
Foxholes
1.27
e
01
6.1
1.27
e
01
7.2
1.27
e
01
6.3
0.53
e
-
10
0.53
e
-
10
1.29
e
01
5.0
e
02
(
9.03
e
-
15
)
(
9.03
e
-
15
)
(
9.03
e
-
15
)
(
0.00
e
00
)
(
0.97
e
-
14
)
(
0.86
e
00
)
(
0.0
e
00
)
Griewank
0.00
e
00
80
0.00
e
00
1332
0.00
e
00
2300
1.06
e
02
6.79
e
01
2.63
e
02
2.06
e
02
(
0.00
e
00
)
(
0.00
e
00
)
(
0.00
e
00
)
(
0.44
e
01
)
(
0.98
e
01
)
(
0.10
e
02
)
(
0.39
e
01
)
Quartic
4.55
e
-
01
2012
5.37
e
-
01
1771
4.88
e
-
01
2951
4.15
e
01
2.09
e
01
1.49
e
00
3.55
e
00
(
3.01
e
-
01
)
(
3.03
e
-
01
)
(
2.71
e
-
01
)
(
0.19
e
01
)
(
0.19
e
01
)
(
0.66
e
00
)
(
0.81
e
00
)
Rastrigin
0.00
e
00
80
0.00
e
00
1282.7
0.00
e
00
2406.7
2.25
e
02
3.08
e
02
2.14
e
02
9.05
e
01
(
0.00
e
00
)
(
0.00
e
00
)
(
0.00
e
00
)
(
0.35
e
02
)
(
0.28
e
01
)
(
0.45
e
02
)
(
0.31
e
02
)
Rosenbrock2
0.00
e
00
6.5
0.00
e
00
5.7
0.00
e
00
6.1
0.49
e
-
04
0.14
e
-
03
0.20
e
-
05
0.55
e
-
04
(
0.00
e
00
)
(
0.00
e
00
)
(
0.00
e
00
)
(
0.11
e
-
03
)
(
0.88
e
-
04
)
(
0.15
e
-
04
)
(
0.17
e
-
04
)
Rosenbrock
0.00
e
00
180
0.00
e
00
2190.7
0.00
e
00
2088
2.20
e
03
2.24
e
03
1.81
e
03
9.14
e
01
(
0.00
e
00
)
(
0.00
e
00
)
(
0.00
e
00
)
(
0.86
e
02
)
(
0.77
e
02
)
(
0.16
e
02
)
(
0.42
e
02
)
Schaffer
0.00
e
00
6.1
0.00
e
00
8
0.00
e
00
4.9
0.24
e
-
01
0.73
e
-
01
-
0.995
e
00
-
1.0
e
00
(
0.00
e
00
)
(
0.00
e
00
)
(
0.00
e
00
)
(
0.42
e
-
02
)
(
0.11
e
-
01
)
(
0.27
e
-
06
)
(
0.0
e
00
)
Spherical
0.00
e
00
12
0.00
e
00
12
0.00
e
00
11.7
0.30
e
-
03
0.30
e
-
03
0.15
e
-
03
0.20
e
-
04
(
0.00
e
00
)
(
0.00
e
00
)
(
0.00
e
00
)
(
0.00
e
00
)
(
0.0
e
00
)
(
0.41
e
-
04
)
(
0.15
e
-
04
)
Step
3.00
e
01
42.7
3.00
e
01
42.7
3.00
e
01
48
0.00
e
00
0.17
e
-
04
0.15
e
-
01
0.25
e
00
(
0.00
e
00
)
(
0.00
e
00
)
(
0.00
e
00
)
(
0.00
e
00
)
(
0.15
e
-
04
)
(
0.0
e
00
)
(
0.60
e
-
01
)To analyze the statistical significance of the results we perform a Friedman test. This is a nonparametric statistical test to determine significant differences in mean for one independent variable with two or more levels—also denoted as treatments—and a dependent variable (or matched groups taken as the problems). The null hypothesis in this test is that the mean ranks assigned to the treatments under testing are the same. Since all three implementations are able to reach the solutions within theη error tolerance on 9 out of 10 problems, the statistical analysis is based on the performance criterion “nfe.” In this hypothesis testing, we have three treatments and ten groups. Friedman’s chi-square has a value of 2.737 (with a P value of 0.255). For 2 degrees of freedom reference χ2 distribution, the critical value for a significance level of 5% is 5.99. Hence, since 2.737 ≤ 5.99, the null hypothesis is not rejected and we conclude that there is no evidence that the three mean ranks values have statistically significant differences.To further compare the sigmoid functions with the rounding to integer strategy, we include in Table3 the results obtained by the “erf⁡” function in (14), together with (10), and the floor function in (15). Only the implementations mCS and mBS are tested. The table also shows the averaged number of iterations, nit. The results illustrate that implementation mCS (Algorithm 2) works very well with strategies based on (14), together with (10), and (15). The success rate for all the problems is 100%, meaning that the algorithms stop because the f value at the position of the best/brightest firefly is within a tolerance η of the optimal solution f*, in all runs. Further, mBS (Algorithm 3) works better when the discretization of the variables is carried out by (15). Overall, mCS based on (14) produces the best results on 6 problems, mCS based on (15) gives the best results on 7 problems (including 4 ties with the former case), mBS based on (14) wins only on one problem, and mBS based on (15) wins on 3 problems (all are ties with mCS based on (15)).Table 3
Comparison of mCS versus mBS and (14) versus (15), based on favg, nfe, and nit.
Prob.
mCS based on (14)
mBS based on (14)
mCS based on (15)
mBS based on (15)
f
avg
nfe
nit
f
avg
nfe
nit
f
avg
nfe
nit
f
avg
nfe
nit
Ackley
8.88
e
-
16
80
1
8.88
e
-
16
678.7
16
8.88
e
-
16
80
1
8.88
e
-
16
82.7
1.1
Foxholes
1.27
e
01
5.7
0.4
1.27
e
01
8
1
1.27
e
01
6.1
0.5
1.29
e
01
404.8
100.2
Griewank
0.00
e
00
80
1
0.00
e
00
717.3
16.9
0.00
e
00
80
1
0.00
e
00
82.7
1.1
Quartic
2.38
e
-
01
81.3
1.03
4.66
e
-
01
794.7
18.9
9.48
e
-
02
80
1
3.93
e
-
01
80
1
Rastrigin
0.00
e
00
80
1
0.00
e
00
702.7
16.6
0.00
e
00
80
1
0.00
e
00
80
1
Rosenbrock2
0.00
e
00
6.4
0.6
0.00
e
00
5.3
0.3
0.00
e
00
6.3
0.6
2.33
e
-
01
470.8
116.7
Rosenbrock
0.00
e
00
80
1
0.00
e
00
105.3
1.6
0.00
e
00
80
1
2.90
e
01
20040
500
Schaffer
0.00
e
00
6.8
0.7
0.00
e
00
12.9
2.2
0.00
e
00
5.1
0.3
7.08
e
-
02
205.1
50.3
Spherical
0.00
e
00
9.9
0.2
0.00
e
00
22.7
1.8
0.00
e
00
10.1
0.3
0.00
e
00
11.5
0.4
Step
3.00
e
01
45.9
0.4
3.00
e
01
62.9
1
3.00
e
01
40.5
0.3
3.00
e
01
40.5
0.3Further, when performing the Friedman test on the four distributions of “nfe” values, the chi-square statistical value is 13.747 (and the P value is 0.0033). From the χ2 distribution table, the critical value for a significance level of 5% and 3 degrees of freedom is 7.81. Since 13.747 > 7.81, the null hypothesis is rejected and we conclude that the observed differences of the four distributions are statistically significant.We now introduce in the statistical analysis the results reported in Tables2 and 3 concerned with both implementations mCS and mBS. Six distributions of “nfe” values are now in comparison. Friedman’s chi-square value is 18.175 (P value = 0.0027). The critical value of the chi-square distribution for a significance level of 5% and 5 degrees of freedom is 11.07. Thus, the null hypothesis of “no significant differences on mean ranks” is rejected and there is evidence that the six distributions of “nfe” values have statistically significant differences. Multiple comparisons (two at a time) may be carried out to determine which mean ranks are significantly different. The estimates of the 95% confidence intervals are shown in the graph of Figure 2 for each case under testing. Two compared distributions of “nfe” are significantly different if their intervals are disjoint and are not significantly different if their intervals overlap. Hence, from the six cases, we conclude that the mean ranks produced by mCS based on (15) are significantly different from those of mBS based on (9) and mBS based on (14). For the remaining pairs of comparison there are no significant differences on the mean ranks.Figure 2
Confidence intervals for mean ranks ofnfe.For comparative purposes we include in Table4 the results obtained by using the proposed Lévy (L) distribution in the randomization term, as shown in (7), and those produced by the Uniform (U) distribution, using rand~U[0,1] as shown in (3). The reported tests use implementation mCS (described in Algorithm 2) with the two heuristics for binary variables: (i) the “erf⁡” function in (14), together with (10), and (ii) the floor function in (15). It is shown that the performance of HBFA with Uniform distribution is very sensitive to the dimension of the problem, since the efficiency is good when n is small but gets worse when n is large. Thus, we have shown that the Lévy distribution is a very good bid.Table 4
Comparison between Lévy and Uniform distributions in the randomization term, based onfavg, nfe, and nit (with St.D. in parentheses).
Prob.
mCS based on (14) + L
mCS based on (14) + U
mCS based on (15) + L
mCS based on (15) + U
f
avg
nfe
nit
f
avg
nfe
nit
f
avg
nfe
nit
f
avg
nfe
nit
Ackley
8.88
e
-
16
80
1
8.88
e
-
16
1096
26.4
8.88
e
-
16
80
1
1.41e00
20040
500
(
0.00
e
00
)
(
0.00
e
00
)
(
0.00
e
00
)
(
1.26
e
-
02
)
Foxholes
1.27
e
01
5.7
0.4
1.27
e
01
8
1
1.27
e
01
6.1
0.5
1.27
e
01
6.8
0.7
(
9.03
e
-
15
)
(
9.03
e
-
15
)
(
9.03
e
-
15
)
(
9.03
e
-
15
)
Griewank
0.00
e
00
80
1
0.00
e
00
2436
59.9
0.00
e
00
80
1
1.27
e
-
01
20040
500
(
0.00
e
00
)
(
0.00
e
00
)
(
0.00
e
00
)
(
2.76
e
-
02
)
Quartic
2.38
e
-
01
81.3
1.03
6.90
e
00
14692
366.2
9.48
e
-
02
80
1
5.10
e
-
01
10581
263.5
(
1.70
e
-
01
)
(
1.16
e
01
)
(
1.03
e
-
01
)
(
2.90
e
-
01
)
Rastrigin
0.00
e
00
80
1
0.00
e
00
1157
27.9
0.00
e
00
80
1
1.00
e
-
01
18104
451.6
(
0.00
e
00
)
(
0.00
e
00
)
(
0.00
e
00
)
(
4.03
e
-
01
)
Rosenbrock2
0.00
e
00
6.4
0.6
0.00
e
00
36.8
8.2
0.00
e
00
6.3
0.6
0.00
e
00
5.9
0.5
(
0.00
e
00
)
(
0.00
e
00
)
(
0.00
e
00
)
(
0.00
e
00
)
Rosenbrock
0.00
e
00
80
1
8.22
e
01
14136
352.4
0.00
e
00
80
1
7.42
e
01
17698
441.2
(
0.00
e
00
)
(
1.02
e
02
)
(
0.00
e
00
)
(
8.43
e
01
)
Schaffer
0.00
e
00
6.8
0.7
0.00
e
00
18.5
3.6
0.00
e
00
5.1
0.3
0.00
e
00
5.6
0.4
(
0.00
e
00
)
(
0.00
e
00
)
(
0.00
e
00
)
(
0.00
e
00
)
Spherical
0.00
e
00
9.9
0.2
0.00
e
00
13.3
0.7
0.00
e
00
10.1
0.3
0.00
e
00
14.7
0.8
(
0.00
e
00
)
(
0.00
e
00
)
(
0.00
e
00
)
(
0.00
e
00
)
Step
3.00
e
01
45.9
0.4
3.00
e
01
46.9
0.5
3.00
e
01
40.5
0.3
3.00
e
01
52.3
0.6
(
0.00
e
00
)
(
0.00
e
00
)
(
0.00
e
00
)
(
0.00
e
00
)We add to some problems withn=30 from Table 1—Ackley, Griewank, Rastrigin, Rosenbrock, and Spherical—three other functions Schwefel 2.22, Schwefel 2.26, and Sum of Different Power to compare our results with those reported in [1]. Schwefel 2.22 is unimodal and for Ω=[-10,10]30, the binary solution is (0,…,0) with f*=0; Schwefel 2.26 is multimodal and in Ω=[-500,500]30, the binary solution is (1,1,…,1) with f*=-25.244129544; Sum of Different Power is unimodal and in Ω=[-1,1]30, the minimum is 0 at (0,…,0). For the results of Table 5, we use HBFA based on mCS, with both “erf⁡” function in (14), together with (10), and the floor function (15). The table reports on the average function values, average number of function evaluations, and success rate (SR). Here, 50 independent runs were carried out to compare with the results shown in [1]. The maximum number of function evaluations therein used was 90000. It is shown that our HBFA outperforms the proposed adaptive binary harmony search (ABHS).Table 5
Comparison of HBFA (with mCS) with ABHS in [1] based on favg, nfe, and SR (%).
Prob.
mCS based on (14)
mCS based on (15)
ABHS in [1]
f
avg
nfe
SR (%)
f
avg
nfe
SR (%)
f
avg
nfe
SR (%)
Ackley
8.88
e
-
16
80
100
8.88
e
-
16
80
100
1.56
e
-
01
62350
90
Griewank
0.00
e
00
80
100
0.00
e
00
80
100
3.30
e
-
02
79758
38
Rastrigin
0.00
e
00
80
100
0.00
e
00
80
100
1.32
e
01
90000
0
Rosenbrock
0.00
e
00
80
100
0.00
e
00
80
100
6.80
e
02
90000
0
Schwefel 2.22
0.00
e
00
80
100
0.00
e
00
80
100
0.00
e
00
59870
100
Schwefel 2.26
-
2.52
e
01
80
100
-
2.47
e
01
10867
87
-
1.195
e
04
90000
0
Spherical
0.00
e
00
80
100
0.00
e
00
80
100
0.00
e
00
62234
100
Sum of Different Power
0.00
e
00
91
100
0.00
e
00
168
100
0.00
e
00
80371
100


### 4.3. Effect of Problem’s Dimension on HBFA Performance

We now consider six problems with varied dimensions from the previous set to analyze the effect of problem’s dimension on the HBFA performance. We test three dimensions:n=50, n=100, and n=200. The algorithm’s parameters are set as previously defined. We remark that the size of the population for all the tested problems and dimensions is 40 points.Table6 contains the results for comparison based on averaged values of f, number of function evaluations, and number of iterations. The “St.D.” of the f values are also displayed. Since the implementation mCS, shown in Algorithm 2, performs better and shows more consistent results than the other two, we tested only mCS based on (14) and mCS based on (15).Table 6
Results for varied dimensions (n=50,100,200), considering m=40.
f
*
n
mCS based on (14)
mCS based on (15)
f
avg
St.D.
nfe
nit
f
avg
St.D.
nfe
nit
Ackley
0.00
e
00
50
8.88
e
-
16
0.00
e
00
80
1
8.88
e
-
16
0.00
e
00
80
1
0.00
e
00
100
8.88
e
-
16
0.00
e
00
80
1
8.88
e
-
16
0.00
e
00
80
1
0.00
e
00
200
8.88
e
-
16
0.00
e
00
80
1
8.88
e
-
16
0.00
e
00
80
1
Griewank
0.00
e
00
50
0.00
e
00
0.00
e
00
80
1
0.00
e
00
0.00
e
00
80
1
0.00
e
00
100
0.00
e
00
0.00
e
00
80
1
0.00
e
00
0.00
e
00
80
1
0.00
e
00
200
0.00
e
00
0.00
e
00
80
1
0.00
e
00
0.00
e
00
80
1
Quartic
0.00
e
00
+ noise
50
2.24
e
-
01
1.95
e
-
01
82.7
1.1
1.43
e
-
01
1.59
e
-
01
80
1
0.00
e
00
+ noise
100
4.32
e
-
01
2.73
e
-
01
146.7
2.7
1.73
e
-
01
1.10
e
-
01
80
1
0.00
e
00
+ noise
200
5.23
e
-
01
2.94
e
-
01
1738.7
42.5
1.96
e
-
01
2.13
e
-
01
81.3
1.03
Rosenbrock
0.00
e
00
50
0.00
e
00
0.00
e
00
80
1
0.00
e
00
0.00
e
00
80
1
0.00
e
00
100
0.00
e
00
0.00
e
00
80
1
0.00
e
00
0.00
e
00
80
1
0.00
e
00
200
0.00
e
00
0.00
e
00
80
1
0.00
e
00
0.00
e
00
80
1
Spherical
0.00
e
00
50
0.00
e
00
0.00
e
00
80
1
0.00
e
00
0.00
e
00
80
1
0.00
e
00
100
0.00
e
00
0.00
e
00
80
1
0.00
e
00
0.00
e
00
80
1
0.00
e
00
200
0.00
e
00
0.00
e
00
80
1
0.00
e
00
0.00
e
00
80
1
Step
3.00
e
02
50
3.00
e
02
0.00
e
00
80
1
3.00
e
02
0.00
e
00
80
1
6.00
e
02
100
6.00
e
02
0.00
e
00
80
1
6.00
e
02
0.00
e
00
80
1
1.20
e
03
200
1.20
e
03
0.00
e
00
80
1
1.20
e
03
0.00
e
00
80
1Besides testing significant differences on the mean ranks produced by the two treatments, mCS based on (14) and mCS based on (15), we also want to determine if the differences on mean ranks produced by problem’s dimension—50, 100, and 200—are statistically significant at a significance level of 5%. Hence, we aim to analyze the effects of two factors “A” and “B.” “A” is the HBFA implementation (with two levels) and “B” is the problem’s dimension (with three levels). For this purpose, the results obtained for the six problems for each combination of the levels of “A” and “B” are considered as replications. When performing the Friedman test for factor “A,” the chi-square statistical value is 1.225 (P value = 0.2685) with 1 degree of freedom. The critical value for a significance level of 5% and 1 degree of freedom in the χ2 distribution table is 3.84, and there is no evidence of statistically significant differences. From the Friedman test for factor “B,” we also conclude that there is no evidence of statistically significant differences, since the chi-square statistical value is 0.746 (P value = 0.6886) with 2 degrees of freedom. (The critical value of the χ2 distribution table for a significance level of 5% and 2 degrees of freedom is 5.99.) Hence, we conclude that the dimension of the problem does not affect the algorithm’s performance. Only with problem Quartic, the efficiency of mCS based on (14) gets worse as dimension increases. Overall, both tested strategies are rather effective when binary solutions are required on small as well as on large nonlinear optimization problems.


### 4.4. Solving 0-1 Knapsack Problems

Finally, we aim to analyze the behavior of our best tested strategies when solving well-known binary and linear optimization problems. For this preliminary experiment, we selected two small knapsack problems. The 0-1 knapsack problem (KP) can be described as follows. Letn be the number of items, from which we have to select some of them to be carried in a knapsack. Let wl and vl be the weight and the value of item l, respectively, and let W be the knapsack’s capacity. It is usually assumed that all weights and values are nonnegative. The objective is to maximize the total value of the knapsack under the constraint of the knapsack’s capacity:
(17)max⁡xV(x)≡∑l=1n‍vlxls.t.F∑l=1n‍wlxl≤W,xl∈{0,1},l=1,…,n.
If item l is selected, xl=1; otherwise, xl=0. Using a penalty function, this problem can be transformed into
(18)min⁡x-∑l=1n‍vlxl+μmax⁡{0,∑l=1n‍wlxl-W},
where μ is the penalty parameter which was set to be 100 in this experiment.Case 1 (an instance of a 0-1 KP with 4 items).
Knapsack’s capacity isW=6 and the vectors of values and weights are v=(40,15,20,10) and w=(4,2,3,1). Based on the above-mentioned parameters, the HBFA with mCS based on (14) was run 30 times and the averaged results were the following. With a success rate of 100%, items 1 and 2 are included in the knapsack and items 3 and 4 are excluded, with a maximum value of 55 (St.D. = 0.0e00). On average, the runs required 0.8 iterations and 29.3 function evaluations. With a success rate of 23%, the heuristic based on the floor function, thus mCS based on (15), reached favg=49 (St.D. = 4.0e00) after an average of 6161.1 function evaluations and an average of 384.1 iterations.Case 2 (an instance of a 0-1 KP with 8 items).
The maximum capacity of the knapsack is set to 8 and the vectors of values and weights arev=(83,14,54,79,72,52,48,62) and w=(3,2,3,2,1,2,2,3). The results are averaged over the 30 runs. After 8.7 iterations and 386.7 function evaluations, the maximum value produced by the strategy mCS based on (14) is 286 (St.D. = 0.0e00), with a success rate of 100%. Items 1, 4, 5, and 6 are included in the knapsack and the others are excluded. The heuristic mCS based on (15) did not reach the optimal solution. All runs required 500 iterations and 20040 function evaluations and the average function values were favg=227 with St.D. = 3.14e01.


## 4.1. Experimental Setting

Each experiment was conducted 30 times. The size of the population is made to depend on the problem’s dimension and is set tom=min⁡{40,2n}. Some experiments have been carried out to tune certain parameters of the algorithms. In the proposed FA with dynamic α and γ, they are set as follows: β0=1, p=1, αmax⁡=0.5, αmin⁡=0.01, γmax⁡=10, and γmin⁡=0.1. In Algorithms 2 (mCS), 3 (mBS), and 4 (pBC), iterations were limited to kmax⁡=500 and the tolerance for finding a good quality solution is η=10-6. Results reported are averaged (over the 30 runs) of best function values, number of function evaluations, and number of iterations.


## 4.2. Experimental Results

First, we use a set of ten benchmark nonlinear functions with different dimensions and characteristics. For example, five functions are unimodal and the remaining multimodal [3, 9, 10, 36]. They are displayed in Table 1. Although they are widely used in continuous optimization, we now aim to converge to a 0/1 bit string solution.Table 1
Problems set.
Ackley
f
(
x
)
=
-
20
exp
⁡
(
-
0.2
1
n
∑
i
=
1
n
‍
x
i
2
)
-
exp
⁡
(
1
n
∑
i
=
1
n
‍
cos
⁡
(
2
π
x
i
)
)
+
20
+
e
n = 30, Ω=[-30,30]30, f*=0  at  x*=(0,…,0)
Foxholes
f
(
x
)
=
1
l
0.002
+
∑
j
=
1
25
‍
(
1
/
(
j
+
(
x
1
-
a
1
j
)
6
+
(
x
2
-
a
2
j
)
6
)
)
[
a
i
j
]
=
[
-
32
-
16
0
16
32
-
32
-
16
⋯
16
32
-
32
-
32
-
32
-
32
-
32
-
16
-
16
⋯
32
32
]
n = 2, Ω=[-65.536,65.536]2, f*≈13 at x*=(0,0)
Griewank
f
(
x
)
=
1
+
1
4000
∑
i
=
1
n
‍
x
i
2
-
∏
i
=
1
n
cos
⁡
(
x
i
i
)
n = 30, Ω=[-300,300]30, f*=0  at  x*=(0,…,0)
Quartic
f
(
x
)
=
∑
i
=
1
n
‍
i
x
i
4
+
U
[
0,1
]
n = 30, Ω=[-1.28,1.28]30, f*=0+noise at x*=(0,…,0)
Rastrigin
f
(
x
)
=
10
n
+
∑
i
=
1
n
‍
(
x
i
2
-
10
cos
⁡
(
2
π
x
i
)
)
n = 30, Ω=[-5.12,5.12]30, f*=0  at  x*=(0,…,0)
Rosenbrock2
f
(
x
)
=
100
(
x
1
2
-
x
2
)
2
+
(
1
-
x
1
)
2
n = 2, Ω=[-2.048,2.048]2, f*=0  at  x*=(1,1)
Rosenbrock
f
(
x
)
=
∑
i
=
1
n
-
1
‍
100
(
x
i
2
-
x
i
+
1
)
2
+
(
1
-
x
i
)
2
n = 30, Ω=[-2.048,2.048]30, f*=0  at  x*=(1,…,1)
Schaffer
f
(
x
)
=
0.5
+
(
sin
⁡
⁡
(
x
1
2
+
x
2
2
)
)
2
-
0.5
(
1
+
0.001
(
x
1
2
+
x
2
2
)
)
2
n = 2, Ω=[-100,100]2, f*=0 at x*=(0,0)
Spherical
f
(
x
)
=
∑
i
=
1
n
‍
x
i
2
n = 3, Ω=[-5.12,5.12]3, f*=0 at x*=(0,0,0)
Step
f
(
x
)
=
6
n
+
∑
i
=
1
n
‍
⌊
x
i
⌋
n = 5, Ω=[-5.12,5.12]5, f*=30 at  x*=(0,0,0,0,0)First, we aim to compare with the results reported in [3, 9, 10]. Due to poor results, the authors in [10] do not recommend the use of ABC to solve binary-valued problems. The other metaheuristics therein implemented are the following:(i)
angle modulated PSO (AMPSO) and angle modulated DE (AMDE) that incorporate a trigonometric function as a bit string generator into the classic PSO and DE algorithms, respectively;(ii)
binary DE and PSO based on the sigmoid logistic function and (10), denoted by binDE and binPSO, respectively.We noticed that the problems Foxholes, Griewank, Rosenbrock, Schaffer, and Step are not correctly described in [3, 9, 10]. Table 2 shows both the averaged best function values (obtained during the 30 runs), favg, with the St.D. in parentheses, and the averaged number of function evaluations, nfe, obtained with the sigmoid logistic function (see in (9)) and (10), while using the three implementations: mCS, mBS, and pBC. Results obtained for these ten functions indicate that our proposal HBFA produces high quality solutions and outperforms the binary versions binPSO and binDE, as well as AMPSO and AMDE. We also note that mCS has the best “nfe” values on 6 problems, mBS is better on 3 problems (one is a tie with mCS), and pBC on 2 problems. Thus, the performance of mCS is the best when compared with those of mBS and pBC. The latter is the least efficient of all, in particular for the large dimensional problems.Table 2
Comparison with AMPSO, binPSO, binDE, and AMDE based onfavg and St.D. (shown in parentheses).
Prob.
mCS based on (9)
mBS based on (9)
pBC based on (9)
AMPSO
binPSO
binDE
AMDE
f
avg
nfe
f
avg
nfe
f
avg
nfe
f
avg
f
avg
f
avg
f
avg
Ackley
8.88
e
-
16
80
8.88
e
-
16
1156
8.88
e
-
16
2168
1.97
e
01
2.01
e
01
1.73
e
01
1.64
e
01
(
0.00
e
00
)
(
0.00
e
00
)
(
0.00
e
00
)
(
0.57
e
-
01
)
(
0.49
e
-
01
)
(
0.31
e
01
)
(
0.76
e
00
)
Foxholes
1.27
e
01
6.1
1.27
e
01
7.2
1.27
e
01
6.3
0.53
e
-
10
0.53
e
-
10
1.29
e
01
5.0
e
02
(
9.03
e
-
15
)
(
9.03
e
-
15
)
(
9.03
e
-
15
)
(
0.00
e
00
)
(
0.97
e
-
14
)
(
0.86
e
00
)
(
0.0
e
00
)
Griewank
0.00
e
00
80
0.00
e
00
1332
0.00
e
00
2300
1.06
e
02
6.79
e
01
2.63
e
02
2.06
e
02
(
0.00
e
00
)
(
0.00
e
00
)
(
0.00
e
00
)
(
0.44
e
01
)
(
0.98
e
01
)
(
0.10
e
02
)
(
0.39
e
01
)
Quartic
4.55
e
-
01
2012
5.37
e
-
01
1771
4.88
e
-
01
2951
4.15
e
01
2.09
e
01
1.49
e
00
3.55
e
00
(
3.01
e
-
01
)
(
3.03
e
-
01
)
(
2.71
e
-
01
)
(
0.19
e
01
)
(
0.19
e
01
)
(
0.66
e
00
)
(
0.81
e
00
)
Rastrigin
0.00
e
00
80
0.00
e
00
1282.7
0.00
e
00
2406.7
2.25
e
02
3.08
e
02
2.14
e
02
9.05
e
01
(
0.00
e
00
)
(
0.00
e
00
)
(
0.00
e
00
)
(
0.35
e
02
)
(
0.28
e
01
)
(
0.45
e
02
)
(
0.31
e
02
)
Rosenbrock2
0.00
e
00
6.5
0.00
e
00
5.7
0.00
e
00
6.1
0.49
e
-
04
0.14
e
-
03
0.20
e
-
05
0.55
e
-
04
(
0.00
e
00
)
(
0.00
e
00
)
(
0.00
e
00
)
(
0.11
e
-
03
)
(
0.88
e
-
04
)
(
0.15
e
-
04
)
(
0.17
e
-
04
)
Rosenbrock
0.00
e
00
180
0.00
e
00
2190.7
0.00
e
00
2088
2.20
e
03
2.24
e
03
1.81
e
03
9.14
e
01
(
0.00
e
00
)
(
0.00
e
00
)
(
0.00
e
00
)
(
0.86
e
02
)
(
0.77
e
02
)
(
0.16
e
02
)
(
0.42
e
02
)
Schaffer
0.00
e
00
6.1
0.00
e
00
8
0.00
e
00
4.9
0.24
e
-
01
0.73
e
-
01
-
0.995
e
00
-
1.0
e
00
(
0.00
e
00
)
(
0.00
e
00
)
(
0.00
e
00
)
(
0.42
e
-
02
)
(
0.11
e
-
01
)
(
0.27
e
-
06
)
(
0.0
e
00
)
Spherical
0.00
e
00
12
0.00
e
00
12
0.00
e
00
11.7
0.30
e
-
03
0.30
e
-
03
0.15
e
-
03
0.20
e
-
04
(
0.00
e
00
)
(
0.00
e
00
)
(
0.00
e
00
)
(
0.00
e
00
)
(
0.0
e
00
)
(
0.41
e
-
04
)
(
0.15
e
-
04
)
Step
3.00
e
01
42.7
3.00
e
01
42.7
3.00
e
01
48
0.00
e
00
0.17
e
-
04
0.15
e
-
01
0.25
e
00
(
0.00
e
00
)
(
0.00
e
00
)
(
0.00
e
00
)
(
0.00
e
00
)
(
0.15
e
-
04
)
(
0.0
e
00
)
(
0.60
e
-
01
)To analyze the statistical significance of the results we perform a Friedman test. This is a nonparametric statistical test to determine significant differences in mean for one independent variable with two or more levels—also denoted as treatments—and a dependent variable (or matched groups taken as the problems). The null hypothesis in this test is that the mean ranks assigned to the treatments under testing are the same. Since all three implementations are able to reach the solutions within theη error tolerance on 9 out of 10 problems, the statistical analysis is based on the performance criterion “nfe.” In this hypothesis testing, we have three treatments and ten groups. Friedman’s chi-square has a value of 2.737 (with a P value of 0.255). For 2 degrees of freedom reference χ2 distribution, the critical value for a significance level of 5% is 5.99. Hence, since 2.737 ≤ 5.99, the null hypothesis is not rejected and we conclude that there is no evidence that the three mean ranks values have statistically significant differences.To further compare the sigmoid functions with the rounding to integer strategy, we include in Table3 the results obtained by the “erf⁡” function in (14), together with (10), and the floor function in (15). Only the implementations mCS and mBS are tested. The table also shows the averaged number of iterations, nit. The results illustrate that implementation mCS (Algorithm 2) works very well with strategies based on (14), together with (10), and (15). The success rate for all the problems is 100%, meaning that the algorithms stop because the f value at the position of the best/brightest firefly is within a tolerance η of the optimal solution f*, in all runs. Further, mBS (Algorithm 3) works better when the discretization of the variables is carried out by (15). Overall, mCS based on (14) produces the best results on 6 problems, mCS based on (15) gives the best results on 7 problems (including 4 ties with the former case), mBS based on (14) wins only on one problem, and mBS based on (15) wins on 3 problems (all are ties with mCS based on (15)).Table 3
Comparison of mCS versus mBS and (14) versus (15), based on favg, nfe, and nit.
Prob.
mCS based on (14)
mBS based on (14)
mCS based on (15)
mBS based on (15)
f
avg
nfe
nit
f
avg
nfe
nit
f
avg
nfe
nit
f
avg
nfe
nit
Ackley
8.88
e
-
16
80
1
8.88
e
-
16
678.7
16
8.88
e
-
16
80
1
8.88
e
-
16
82.7
1.1
Foxholes
1.27
e
01
5.7
0.4
1.27
e
01
8
1
1.27
e
01
6.1
0.5
1.29
e
01
404.8
100.2
Griewank
0.00
e
00
80
1
0.00
e
00
717.3
16.9
0.00
e
00
80
1
0.00
e
00
82.7
1.1
Quartic
2.38
e
-
01
81.3
1.03
4.66
e
-
01
794.7
18.9
9.48
e
-
02
80
1
3.93
e
-
01
80
1
Rastrigin
0.00
e
00
80
1
0.00
e
00
702.7
16.6
0.00
e
00
80
1
0.00
e
00
80
1
Rosenbrock2
0.00
e
00
6.4
0.6
0.00
e
00
5.3
0.3
0.00
e
00
6.3
0.6
2.33
e
-
01
470.8
116.7
Rosenbrock
0.00
e
00
80
1
0.00
e
00
105.3
1.6
0.00
e
00
80
1
2.90
e
01
20040
500
Schaffer
0.00
e
00
6.8
0.7
0.00
e
00
12.9
2.2
0.00
e
00
5.1
0.3
7.08
e
-
02
205.1
50.3
Spherical
0.00
e
00
9.9
0.2
0.00
e
00
22.7
1.8
0.00
e
00
10.1
0.3
0.00
e
00
11.5
0.4
Step
3.00
e
01
45.9
0.4
3.00
e
01
62.9
1
3.00
e
01
40.5
0.3
3.00
e
01
40.5
0.3Further, when performing the Friedman test on the four distributions of “nfe” values, the chi-square statistical value is 13.747 (and the P value is 0.0033). From the χ2 distribution table, the critical value for a significance level of 5% and 3 degrees of freedom is 7.81. Since 13.747 > 7.81, the null hypothesis is rejected and we conclude that the observed differences of the four distributions are statistically significant.We now introduce in the statistical analysis the results reported in Tables2 and 3 concerned with both implementations mCS and mBS. Six distributions of “nfe” values are now in comparison. Friedman’s chi-square value is 18.175 (P value = 0.0027). The critical value of the chi-square distribution for a significance level of 5% and 5 degrees of freedom is 11.07. Thus, the null hypothesis of “no significant differences on mean ranks” is rejected and there is evidence that the six distributions of “nfe” values have statistically significant differences. Multiple comparisons (two at a time) may be carried out to determine which mean ranks are significantly different. The estimates of the 95% confidence intervals are shown in the graph of Figure 2 for each case under testing. Two compared distributions of “nfe” are significantly different if their intervals are disjoint and are not significantly different if their intervals overlap. Hence, from the six cases, we conclude that the mean ranks produced by mCS based on (15) are significantly different from those of mBS based on (9) and mBS based on (14). For the remaining pairs of comparison there are no significant differences on the mean ranks.Figure 2
Confidence intervals for mean ranks ofnfe.For comparative purposes we include in Table4 the results obtained by using the proposed Lévy (L) distribution in the randomization term, as shown in (7), and those produced by the Uniform (U) distribution, using rand~U[0,1] as shown in (3). The reported tests use implementation mCS (described in Algorithm 2) with the two heuristics for binary variables: (i) the “erf⁡” function in (14), together with (10), and (ii) the floor function in (15). It is shown that the performance of HBFA with Uniform distribution is very sensitive to the dimension of the problem, since the efficiency is good when n is small but gets worse when n is large. Thus, we have shown that the Lévy distribution is a very good bid.Table 4
Comparison between Lévy and Uniform distributions in the randomization term, based onfavg, nfe, and nit (with St.D. in parentheses).
Prob.
mCS based on (14) + L
mCS based on (14) + U
mCS based on (15) + L
mCS based on (15) + U
f
avg
nfe
nit
f
avg
nfe
nit
f
avg
nfe
nit
f
avg
nfe
nit
Ackley
8.88
e
-
16
80
1
8.88
e
-
16
1096
26.4
8.88
e
-
16
80
1
1.41e00
20040
500
(
0.00
e
00
)
(
0.00
e
00
)
(
0.00
e
00
)
(
1.26
e
-
02
)
Foxholes
1.27
e
01
5.7
0.4
1.27
e
01
8
1
1.27
e
01
6.1
0.5
1.27
e
01
6.8
0.7
(
9.03
e
-
15
)
(
9.03
e
-
15
)
(
9.03
e
-
15
)
(
9.03
e
-
15
)
Griewank
0.00
e
00
80
1
0.00
e
00
2436
59.9
0.00
e
00
80
1
1.27
e
-
01
20040
500
(
0.00
e
00
)
(
0.00
e
00
)
(
0.00
e
00
)
(
2.76
e
-
02
)
Quartic
2.38
e
-
01
81.3
1.03
6.90
e
00
14692
366.2
9.48
e
-
02
80
1
5.10
e
-
01
10581
263.5
(
1.70
e
-
01
)
(
1.16
e
01
)
(
1.03
e
-
01
)
(
2.90
e
-
01
)
Rastrigin
0.00
e
00
80
1
0.00
e
00
1157
27.9
0.00
e
00
80
1
1.00
e
-
01
18104
451.6
(
0.00
e
00
)
(
0.00
e
00
)
(
0.00
e
00
)
(
4.03
e
-
01
)
Rosenbrock2
0.00
e
00
6.4
0.6
0.00
e
00
36.8
8.2
0.00
e
00
6.3
0.6
0.00
e
00
5.9
0.5
(
0.00
e
00
)
(
0.00
e
00
)
(
0.00
e
00
)
(
0.00
e
00
)
Rosenbrock
0.00
e
00
80
1
8.22
e
01
14136
352.4
0.00
e
00
80
1
7.42
e
01
17698
441.2
(
0.00
e
00
)
(
1.02
e
02
)
(
0.00
e
00
)
(
8.43
e
01
)
Schaffer
0.00
e
00
6.8
0.7
0.00
e
00
18.5
3.6
0.00
e
00
5.1
0.3
0.00
e
00
5.6
0.4
(
0.00
e
00
)
(
0.00
e
00
)
(
0.00
e
00
)
(
0.00
e
00
)
Spherical
0.00
e
00
9.9
0.2
0.00
e
00
13.3
0.7
0.00
e
00
10.1
0.3
0.00
e
00
14.7
0.8
(
0.00
e
00
)
(
0.00
e
00
)
(
0.00
e
00
)
(
0.00
e
00
)
Step
3.00
e
01
45.9
0.4
3.00
e
01
46.9
0.5
3.00
e
01
40.5
0.3
3.00
e
01
52.3
0.6
(
0.00
e
00
)
(
0.00
e
00
)
(
0.00
e
00
)
(
0.00
e
00
)We add to some problems withn=30 from Table 1—Ackley, Griewank, Rastrigin, Rosenbrock, and Spherical—three other functions Schwefel 2.22, Schwefel 2.26, and Sum of Different Power to compare our results with those reported in [1]. Schwefel 2.22 is unimodal and for Ω=[-10,10]30, the binary solution is (0,…,0) with f*=0; Schwefel 2.26 is multimodal and in Ω=[-500,500]30, the binary solution is (1,1,…,1) with f*=-25.244129544; Sum of Different Power is unimodal and in Ω=[-1,1]30, the minimum is 0 at (0,…,0). For the results of Table 5, we use HBFA based on mCS, with both “erf⁡” function in (14), together with (10), and the floor function (15). The table reports on the average function values, average number of function evaluations, and success rate (SR). Here, 50 independent runs were carried out to compare with the results shown in [1]. The maximum number of function evaluations therein used was 90000. It is shown that our HBFA outperforms the proposed adaptive binary harmony search (ABHS).Table 5
Comparison of HBFA (with mCS) with ABHS in [1] based on favg, nfe, and SR (%).
Prob.
mCS based on (14)
mCS based on (15)
ABHS in [1]
f
avg
nfe
SR (%)
f
avg
nfe
SR (%)
f
avg
nfe
SR (%)
Ackley
8.88
e
-
16
80
100
8.88
e
-
16
80
100
1.56
e
-
01
62350
90
Griewank
0.00
e
00
80
100
0.00
e
00
80
100
3.30
e
-
02
79758
38
Rastrigin
0.00
e
00
80
100
0.00
e
00
80
100
1.32
e
01
90000
0
Rosenbrock
0.00
e
00
80
100
0.00
e
00
80
100
6.80
e
02
90000
0
Schwefel 2.22
0.00
e
00
80
100
0.00
e
00
80
100
0.00
e
00
59870
100
Schwefel 2.26
-
2.52
e
01
80
100
-
2.47
e
01
10867
87
-
1.195
e
04
90000
0
Spherical
0.00
e
00
80
100
0.00
e
00
80
100
0.00
e
00
62234
100
Sum of Different Power
0.00
e
00
91
100
0.00
e
00
168
100
0.00
e
00
80371
100


## 4.3. Effect of Problem’s Dimension on HBFA Performance

We now consider six problems with varied dimensions from the previous set to analyze the effect of problem’s dimension on the HBFA performance. We test three dimensions:n=50, n=100, and n=200. The algorithm’s parameters are set as previously defined. We remark that the size of the population for all the tested problems and dimensions is 40 points.Table6 contains the results for comparison based on averaged values of f, number of function evaluations, and number of iterations. The “St.D.” of the f values are also displayed. Since the implementation mCS, shown in Algorithm 2, performs better and shows more consistent results than the other two, we tested only mCS based on (14) and mCS based on (15).Table 6
Results for varied dimensions (n=50,100,200), considering m=40.
f
*
n
mCS based on (14)
mCS based on (15)
f
avg
St.D.
nfe
nit
f
avg
St.D.
nfe
nit
Ackley
0.00
e
00
50
8.88
e
-
16
0.00
e
00
80
1
8.88
e
-
16
0.00
e
00
80
1
0.00
e
00
100
8.88
e
-
16
0.00
e
00
80
1
8.88
e
-
16
0.00
e
00
80
1
0.00
e
00
200
8.88
e
-
16
0.00
e
00
80
1
8.88
e
-
16
0.00
e
00
80
1
Griewank
0.00
e
00
50
0.00
e
00
0.00
e
00
80
1
0.00
e
00
0.00
e
00
80
1
0.00
e
00
100
0.00
e
00
0.00
e
00
80
1
0.00
e
00
0.00
e
00
80
1
0.00
e
00
200
0.00
e
00
0.00
e
00
80
1
0.00
e
00
0.00
e
00
80
1
Quartic
0.00
e
00
+ noise
50
2.24
e
-
01
1.95
e
-
01
82.7
1.1
1.43
e
-
01
1.59
e
-
01
80
1
0.00
e
00
+ noise
100
4.32
e
-
01
2.73
e
-
01
146.7
2.7
1.73
e
-
01
1.10
e
-
01
80
1
0.00
e
00
+ noise
200
5.23
e
-
01
2.94
e
-
01
1738.7
42.5
1.96
e
-
01
2.13
e
-
01
81.3
1.03
Rosenbrock
0.00
e
00
50
0.00
e
00
0.00
e
00
80
1
0.00
e
00
0.00
e
00
80
1
0.00
e
00
100
0.00
e
00
0.00
e
00
80
1
0.00
e
00
0.00
e
00
80
1
0.00
e
00
200
0.00
e
00
0.00
e
00
80
1
0.00
e
00
0.00
e
00
80
1
Spherical
0.00
e
00
50
0.00
e
00
0.00
e
00
80
1
0.00
e
00
0.00
e
00
80
1
0.00
e
00
100
0.00
e
00
0.00
e
00
80
1
0.00
e
00
0.00
e
00
80
1
0.00
e
00
200
0.00
e
00
0.00
e
00
80
1
0.00
e
00
0.00
e
00
80
1
Step
3.00
e
02
50
3.00
e
02
0.00
e
00
80
1
3.00
e
02
0.00
e
00
80
1
6.00
e
02
100
6.00
e
02
0.00
e
00
80
1
6.00
e
02
0.00
e
00
80
1
1.20
e
03
200
1.20
e
03
0.00
e
00
80
1
1.20
e
03
0.00
e
00
80
1Besides testing significant differences on the mean ranks produced by the two treatments, mCS based on (14) and mCS based on (15), we also want to determine if the differences on mean ranks produced by problem’s dimension—50, 100, and 200—are statistically significant at a significance level of 5%. Hence, we aim to analyze the effects of two factors “A” and “B.” “A” is the HBFA implementation (with two levels) and “B” is the problem’s dimension (with three levels). For this purpose, the results obtained for the six problems for each combination of the levels of “A” and “B” are considered as replications. When performing the Friedman test for factor “A,” the chi-square statistical value is 1.225 (P value = 0.2685) with 1 degree of freedom. The critical value for a significance level of 5% and 1 degree of freedom in the χ2 distribution table is 3.84, and there is no evidence of statistically significant differences. From the Friedman test for factor “B,” we also conclude that there is no evidence of statistically significant differences, since the chi-square statistical value is 0.746 (P value = 0.6886) with 2 degrees of freedom. (The critical value of the χ2 distribution table for a significance level of 5% and 2 degrees of freedom is 5.99.) Hence, we conclude that the dimension of the problem does not affect the algorithm’s performance. Only with problem Quartic, the efficiency of mCS based on (14) gets worse as dimension increases. Overall, both tested strategies are rather effective when binary solutions are required on small as well as on large nonlinear optimization problems.


## 4.4. Solving 0-1 Knapsack Problems

Finally, we aim to analyze the behavior of our best tested strategies when solving well-known binary and linear optimization problems. For this preliminary experiment, we selected two small knapsack problems. The 0-1 knapsack problem (KP) can be described as follows. Letn be the number of items, from which we have to select some of them to be carried in a knapsack. Let wl and vl be the weight and the value of item l, respectively, and let W be the knapsack’s capacity. It is usually assumed that all weights and values are nonnegative. The objective is to maximize the total value of the knapsack under the constraint of the knapsack’s capacity:
(17)max⁡xV(x)≡∑l=1n‍vlxls.t.F∑l=1n‍wlxl≤W,xl∈{0,1},l=1,…,n.
If item l is selected, xl=1; otherwise, xl=0. Using a penalty function, this problem can be transformed into
(18)min⁡x-∑l=1n‍vlxl+μmax⁡{0,∑l=1n‍wlxl-W},
where μ is the penalty parameter which was set to be 100 in this experiment.Case 1 (an instance of a 0-1 KP with 4 items).
Knapsack’s capacity isW=6 and the vectors of values and weights are v=(40,15,20,10) and w=(4,2,3,1). Based on the above-mentioned parameters, the HBFA with mCS based on (14) was run 30 times and the averaged results were the following. With a success rate of 100%, items 1 and 2 are included in the knapsack and items 3 and 4 are excluded, with a maximum value of 55 (St.D. = 0.0e00). On average, the runs required 0.8 iterations and 29.3 function evaluations. With a success rate of 23%, the heuristic based on the floor function, thus mCS based on (15), reached favg=49 (St.D. = 4.0e00) after an average of 6161.1 function evaluations and an average of 384.1 iterations.Case 2 (an instance of a 0-1 KP with 8 items).
The maximum capacity of the knapsack is set to 8 and the vectors of values and weights arev=(83,14,54,79,72,52,48,62) and w=(3,2,3,2,1,2,2,3). The results are averaged over the 30 runs. After 8.7 iterations and 386.7 function evaluations, the maximum value produced by the strategy mCS based on (14) is 286 (St.D. = 0.0e00), with a success rate of 100%. Items 1, 4, 5, and 6 are included in the knapsack and the others are excluded. The heuristic mCS based on (15) did not reach the optimal solution. All runs required 500 iterations and 20040 function evaluations and the average function values were favg=227 with St.D. = 3.14e01.


## 5. Conclusions and Future Work

In this work we have implemented several heuristics to compute a global optimal binary solution of bound constrained nonlinear optimization problems, which have been incorporated into FA, yielding the herein called HBFA. The problems addressed in this study have bounded continuous search space. Our FA proposal uses dynamic updating schemes for two parameters,γ from the attractiveness term and α from the randomization term, and considers the Lévy distribution to create randomness in firefly movement. The performance behavior of the proposed heuristics has been investigated. Three simple heuristics capable of transforming real continuous variables into binary ones are implemented. A new sigmoid “erf⁡” function is proposed. In the context of the firefly algorithm, three different implementations aiming to incorporate the heuristics for binary variables into FA are proposed (mCS, mBS, and pBC). Based on a set of benchmark problems, a comparison is carried out with other binary dealing metaheuristics, namely, AMPSO, binPSO, binDE, and AMDE. The experimental results show that the implementation denoted by mCS when combined with either the new sigmoid “erf⁡” function or the rounding scheme based on the floor function is quite efficient and superior to the other methods in comparison. The statistical analysis carried out on the results shows evidence of statistically significant differences on efficiency, measured by the number of function evaluations, between the implementation mCS based on the floor function approach and the mBS based on both tested sigmoid functions schemes. We have also investigated the effect of problem’s dimension on the performance of our algorithm. Using the Friedman statistical test we conclude that the differences on efficiency are not statistically significant. Another simple experiment has shown that the implementation mCS with the sigmoid “erf⁡” function is effective in solving two small 0-1 KP. The performance of this simple heuristic strategy will be further analyzed to solve large and multidimensional 0-1 KP. Future developments concerning the HBFA will consider its extension to deal with integer variables in nonlinear optimization problems. Different heuristics to transform continuous real variables into integer variables will be investigated. Challenging mixed-integer nonconvex nonlinear problems will be solved.



---

*Source: 215182-2014-10-08.xml*