# A Novel BBO Algorithm Based on Local Search and Nonuniform Variation for Iris Classification

Authors: Lisheng Wei; Ning Wang; Huacai Lu
Journal: Complexity
 (2021)

---

## Abstract

In order to improve the iris classification rate, a novel biogeography-based optimization algorithm (NBBO) based on local search and nonuniform variation was proposed in this paper. Firstly, the linear migration model was replaced by a hyperbolic cotangent model which was closer to the natural law. And, the local search strategy was added to traditional BBO algorithm migration operation to enhance the global search ability of the algorithm. Then, the nonuniform variation was introduced to enhance the algorithm in the later iteration. The algorithm could achieve a stronger iris classifier by lifting weaker similarity classifiers during the training stage. On this base, the convergence condition of NBBO was proposed by using the Markov chain strategy. Finally, simulation results were given to demonstrate the effectiveness and efficiency of the proposed iris classification method.


---

## Body

## 1. Introduction

Iris classification is especially suitable for recognition with uniqueness, stability, inviolability, and reliability characteristics and has been one of the hottest biological characteristic recognition research spot recently [1]. It has been widely used technology in national defence, financial industry, and entrance guard system [2, 3]. Alfred Wallace and Charles Darwin proposed the theory of biogeography to study the distribution, migration, and extinction of habitats of biological species in the 19th century. Inspired by biogeography, Simon proposed the BBO algorithm in 2008 [4]. BBO is a new swarm intelligence optimization algorithm and solves the optimization problem by simulating the mathematical model of species migration in biogeography. In general, BBO can find a good initial candidate solution and generate an acceptable optimization solution, which has been widely used in image processing [5], scheduling optimization [6, 7], parameter estimation [8], power flow calculation [9], and load analysis [10]. However, BBO has some shortcomings, such as weak search ability and easy to fall into local optimum in the later period of operation. How to use BBO global search ability to improve the iris classification rate has important theoretical and practical value.In recent years, many researchers have worked on the improvement of the BBO performance. For instance, Wang et al. [11] combined the chaotic mapping strategy with the BBO optimal migration model and proposed a biogeography optimization algorithm based on based on the adaptive population migration mechanism. Feng et al. [12] proposed an improved BBO with random ring topology and Powell’s method, in which the local loop topology was used instead of the global topology. The adaptive Powell method was modified to adapt the evolutionary process to improve the accuracy of the solution. Li et al. [13] designed a hybrid algorithm by combining the artificial bee colony algorithm (ABC) with BBO to obtain ABC’s exploring ability and BBO’s developing ability. Chen et al. [14] combined Cuckoo Search (CS) and bBBO using two search strategies: heterogeneous rhododendron search and biogeography-based discovery. Zhao et al. [15] proposed an optimization method based on two-stage differential biogeography to solve the problem of premature convergence and reduce rotation variance.In order to improve the iris classification rate, the NBBO based on local search and nonuniform variation was proposed in this paper. Firstly, the hyperbolic cotangent migration model was used to replace the original linear migration model to improve the adaptability of the algorithm. Then, the local search strategy was added to the BBO migration operation in order to enhance the algorithm search ability and improve the convergence speed of the algorithm. The nonuniform variation was used in the operation to increase the search ability of the algorithm and prevent the local optimum, which could achieve a stronger iris classifier by lifting weaker similarity classifiers during the training stage.The rest of the paper was organized as follows. Section2 presented the improved biogeography-based optimization algorithm, including migration model, migration operations, and nonuniform variation operation. Section 3 presented NBBO test on different test functions and iris classification samples. Finally, a brief conclusion and future work are given in Section 4.


## 2. Methods

In nature, biological species live in different areas with obvious boundaries, which are called habitats. Habitat suitability index (HSI) is used to describe the survival degree of species suitable for habitat. HSI is affected by humidity, temperature, habitat area, and vegetation diversity. These factors are called suitability variables (SIV). Habitats with higher HSI can accommodate more biological species, whereas habitats with lower HSI can only accommodate fewer biological species. BBO is a population-based algorithm, which regards each solution as a habitat, fitness as HSI, and each component of solution as a SIV. The optimization problem can be solved by simulating the migration and mutation process in biogeography. The two main operations of BBO are migration operation and mutation operation [16, 17].Like other intelligent algorithms, BBO is prone to oscillation in the iterative process, resulting in slow and premature convergence, poor local search ability, and inaccuracy. This section presents the NBBO algorithm in terms of migration model and operation. The flow chart of the improved algorithm is shown in Figure1.Figure 1
NBBO algorithm flowchart.


### 2.1. Migration Model

There are different migration models in biogeography. As shown in formula (1), the original BBO adopts the simplest linear migration model. But, the ecosystem is nonlinear in nature; small changes in a certain part of the system may have complex effects on the whole system. For the actual species, the migration process is more complicated than what is described by the linear model. Among the related improved algorithms, the cosine mobility model proposed by Ma performs best under most of the functions tested [18], and the mobility model close to the natural law is better than the simple linear model. In order to better adapt to the nonlinear migration problem, this paper proposes a hyperbolic cotangent nonlinear migration model to improve the performance of the basic BBO algorithm.(1)λk=I⋅1−kn,μk=E⋅kn,(2)λk=I2⋅−ekπ/n−e−kπ/nekπ/n+e−kπ/n+1,(3)μk=E2⋅ekπ/n−e−kπ/nekπ/n+e−kπ/n+1,where λk and μk are related to the number of species k and n is the maximum number of species. λk represents the immigration rate which refers to the probability of species moving into the habitat, μk represents the emigration rate which refers to the probability of species moving out of the habitat, I represents the maximum immigration rate, and E represents the maximum emigration rate, all of which are functions of the number of habitat populations k. When the number of habitat species is 0, the immigration rate is the highest, that is, λk=I, the migration rate is the smallest, that is, μk=0. When the number of species reaches the maximum number of species n that can be accommodated in the habitat, the immigration rate λk drops to 0, and the migration rate reaches the maximum, that is, μk=E. When the number of species is k0, the immigration rate and the emigration rate of the species are equal, that is, λk=μk, and the dynamic equilibrium state is reached at this time.In the hyperbolic tangential mobility model, the trend of mobility changing with species number is like that of cosine model, but the amplitude is more moderate than that of cosine model.


### 2.2. Local Search Migration Operations

In the original BBO, a generalized search was carried out based on the immigration rate and the emigration rate, a habitat with a higher emigration rate was randomly selected, and the characteristics of the corresponding position (SIV) were moved to the habitat with low HSI. However, in some habitats with high migration rates, the characteristics of the corresponding locations SIV are not necessarily the best for the migrated habitats. If migrated directly, it may result in a decrease in the suitability index of the move into habitat.Therefore, in order to enhance the global search ability of the algorithm and improve the convergence speed of the algorithm, the selected migration characteristics can be modified by using equation (6), which can reduce the poor migration characteristics from the habitat with higher migration rate to the low migration rate. Similarly, if you move out of the same habitat with a high migration rate each time, as the number of evolutions increases, the population diversity decreases, and the convergence speed will slow down:(4)C=MinVal+rand∗MaxVal−MinVal,(5)randTemp=rand1,2,(6)X′k,j=Xk,j+C∗randTemp∗Xk,j,where C represents a space allocated to each individual in the population for storing algorithm parameters. Xk,j represents the current individual.X′k,j represents the new individual generated by the local search. MinVal=−1 and MaxVal=1.


### 2.3. Nonuniform Variation Operation

In the early stage of iteration, the population position is usually far from the optimal solution and needs to be iterated many times to find a feasible search space when the search radius is too small. At the end of iteration, the population position is close to or near the optimal solution. A small range of searches is used as a fine-tuning of the solution vector. It is shown by equation (7), as the traditional mutation operator cannot search the optimal solution efficiently, the nonuniform mutation operator is introduced. The whole population is searched at a large range during the early stage of search. As the number of iterations increases, the search range is gradually narrowed, and the algorithm is effectively balanced to escape from local minimum.(7)ms=mmax⋅1−PsPmax,where Ps denotes the probability for the S number of species; S denotes the species number; and mmax denotes the maximum variation rate given; Pmax=maxPs.The specific operation of nonuniform variation is like uniform variation, but it focuses on searching for small areas near the original individual. Assuming that the variation of thekth component Xk of the t generation (t is the current number of iterations) of Habitat Xit=X1,X2,…,Xk,…,Xn is performed, the variation component is as follows:(8)Xk′=Xk+Δt,Umax−Xk,if rand<0.5,(9)Xk′=Xk−Δt,Xk−Umin,if rand≥0.5,where Umax and Umin are the upper and lower limits of component Xk, respectively. Umax−Xk and Xk−Umin represent random numbers in the range of 0,y and meet the nonuniform distribution:(10)limt⟶∞Δt,y=0.This feature allows the mutation operator to search the entire space evenly in the initial stage of the algorithm (whent is very small), and the latter part of the algorithm focuses on several local range exact searches. Δt,y can be defined as follows:(11)Δt,y=y⋅1−r1−t/T⋅b,where r is a random number that fits evenly within the 0,1 range, T is the maximum number of iterations, and b is the system parameter that determines the degree of dependency on the number of iterations, select here b=2.Then, the specific steps of NBBO algorithm are as follows: (Algorithm1)Algorithm 1:The main procedure of NBBO.
BeginAlgorithmic parameters setting:Population size = 50,Maximum number of iterationsG = 500,Maximum immigration rateI=1,Maximum emigration rateE=1,Mutation probabilityPm=0.05,Number of elites retainedKeep=2./∗ Initialization ∗/Random generation of a set of initial habitats constitutes the initial population.Calculate the fitness value of each habitat in the population.While notT /∗T is the condition for the end of the iteration ∗/Sort the habitat in descending order according to the fitness value.The immigration rateλi and emigration rate μi of each habitat were calculated according to formula (1)./∗ End of initialization∗//∗ Migration ∗/piSelectHi according to the immigration probability λiforj=1:D /∗D is dimension ∗/if rand<λiSelectHj according to the migration probability μiRandomly select a characteristic variableSIV from HjReplace a random characteristic variableSIV in Hi with SIVelsePerform a local search on the characteristic variableSIVaccording to formula (6) can get SIV′.Replace random characteristic variableSIV in Hi with SIV′.end ifend forend for/∗ End of migration ∗//∗ Mutation ∗/fori=1:SizeCalculate the probability of mutationPs according to λi and μiSelection of unmutated habitatHi based on Psif rand0,1<msthenif rand<0.5ReplaceHiSIV with a randomly generated SIV by formula (8)if rand≥0.5ReplaceHiSIV with a randomly generated SIV by (9)end ifend ifend ifend for/∗ End of mutation ∗/Recalculate the habitat fitness valuesend whileend


### 2.4. Convergence Analysis

By using Markov chain strategy, we can analyze the convergence of the proposed NBBO algorithms [19, 20]. Let Xt=x1t,x2t,…,xnt be a population with improved algorithm algebra of t. We divide Xt into n subsets based on the fitness of the population. xiti=1,2,…,n is a subset of them; next, we use Markov chain to analyze the global convergence of the improved BBO algorithm.The evolution process of the NBBO algorithm is mainly composed of selection, migration and mutation operations, which are independent of evolution algebra. Considering the number of individuals of BBO algorithm, the evolution process of the NBBO algorithm satisfies the finite homogeneous Markov chain. Population subsetxiti=1,2,…,n can be equivalent to a state on a finite homogeneous Markov chain, pi means the probability of being in state xit, and then pi>0 and ∑i=1npi=1.The transition probabilitypij=pxi⟶xj of state xit is as follows:(12)pij=0,i<jpij>0,i≥j.Then, the state transition matrixP of the Markov chain is(13)P=p110⋯0p21p22⋯0⋮⋮⋱⋮pn1pn2⋯pnn=p110Q2Q1,where(14)Q1=p22⋯0⋮⋱⋮pn2⋯pnn,Q2=p21,p31,⋯,pn1T.It can be found by formulas (12) and (13) Q2≠0, Q1≠0. Since the sum of the probabilities of each row in the Markov transition matrix is 1, p11=1≠0. It can be concluded that the state transition matrix P of Markov chain satisfies the reducible random matrix. According to [21], the definition of a reductive random matrix is as follows:Let matrixZ of order N be a reducible random matrix which can be obtained by the same row transformation and column transformation, namely,(15)Z=C0RT,where C is the mm≤N order primitive matrix and R and T are nonzero matrices.Then, we have(16)Z∞=limk⟶∞Zk=limk⟶∞Ck0∑i=1k−1TiRCk−iTk=C∞0R∞0.The above matrix is a stable random matrix, andZ∞ satisfies the following conditions:(17)Z∞=ZijN∗N,Zij>0,1≤i≤N,1≤j≤m,Zij=0,1≤i≤N,m<j≤N.Therefore, state transition matrixP is a reductive random matrix, then(18)P∞=limt⟶∞p110Q2Q1t=limt⟶∞p11t0∑i=1t−1Q1iQ2p11t−iQ1t=limt⟶∞p11∞0Q2∞0.In formula (12) p11∞=p11=1 and Q1∞=0, since the sum of the probabilities of each line in the Markov transition matrix is 1; there must be Q2∞=1,1,…,1T, then(19)P∞=10⋯010⋯0⋮⋮⋱⋮10⋯0.Formula (19) shows that when the number of iterations t⟶∞, probability pi1∞=1i=1,2,…,n, That is, no matter how the initial state can finally converge to the global optimal solution with probability 1, it can be concluded that(20)limt⟶∞pYt=f∗=1,where Yt is the optimal fitness value of population Xt, f∗ represents the global optimal fitness value, and pYt=f∗ indicates the probability that the optimal individual of the t generation is the global optimal value. It is proved that the NBBO algorithm has global convergence.


## 2.1. Migration Model

There are different migration models in biogeography. As shown in formula (1), the original BBO adopts the simplest linear migration model. But, the ecosystem is nonlinear in nature; small changes in a certain part of the system may have complex effects on the whole system. For the actual species, the migration process is more complicated than what is described by the linear model. Among the related improved algorithms, the cosine mobility model proposed by Ma performs best under most of the functions tested [18], and the mobility model close to the natural law is better than the simple linear model. In order to better adapt to the nonlinear migration problem, this paper proposes a hyperbolic cotangent nonlinear migration model to improve the performance of the basic BBO algorithm.(1)λk=I⋅1−kn,μk=E⋅kn,(2)λk=I2⋅−ekπ/n−e−kπ/nekπ/n+e−kπ/n+1,(3)μk=E2⋅ekπ/n−e−kπ/nekπ/n+e−kπ/n+1,where λk and μk are related to the number of species k and n is the maximum number of species. λk represents the immigration rate which refers to the probability of species moving into the habitat, μk represents the emigration rate which refers to the probability of species moving out of the habitat, I represents the maximum immigration rate, and E represents the maximum emigration rate, all of which are functions of the number of habitat populations k. When the number of habitat species is 0, the immigration rate is the highest, that is, λk=I, the migration rate is the smallest, that is, μk=0. When the number of species reaches the maximum number of species n that can be accommodated in the habitat, the immigration rate λk drops to 0, and the migration rate reaches the maximum, that is, μk=E. When the number of species is k0, the immigration rate and the emigration rate of the species are equal, that is, λk=μk, and the dynamic equilibrium state is reached at this time.In the hyperbolic tangential mobility model, the trend of mobility changing with species number is like that of cosine model, but the amplitude is more moderate than that of cosine model.


## 2.2. Local Search Migration Operations

In the original BBO, a generalized search was carried out based on the immigration rate and the emigration rate, a habitat with a higher emigration rate was randomly selected, and the characteristics of the corresponding position (SIV) were moved to the habitat with low HSI. However, in some habitats with high migration rates, the characteristics of the corresponding locations SIV are not necessarily the best for the migrated habitats. If migrated directly, it may result in a decrease in the suitability index of the move into habitat.Therefore, in order to enhance the global search ability of the algorithm and improve the convergence speed of the algorithm, the selected migration characteristics can be modified by using equation (6), which can reduce the poor migration characteristics from the habitat with higher migration rate to the low migration rate. Similarly, if you move out of the same habitat with a high migration rate each time, as the number of evolutions increases, the population diversity decreases, and the convergence speed will slow down:(4)C=MinVal+rand∗MaxVal−MinVal,(5)randTemp=rand1,2,(6)X′k,j=Xk,j+C∗randTemp∗Xk,j,where C represents a space allocated to each individual in the population for storing algorithm parameters. Xk,j represents the current individual.X′k,j represents the new individual generated by the local search. MinVal=−1 and MaxVal=1.


## 2.3. Nonuniform Variation Operation

In the early stage of iteration, the population position is usually far from the optimal solution and needs to be iterated many times to find a feasible search space when the search radius is too small. At the end of iteration, the population position is close to or near the optimal solution. A small range of searches is used as a fine-tuning of the solution vector. It is shown by equation (7), as the traditional mutation operator cannot search the optimal solution efficiently, the nonuniform mutation operator is introduced. The whole population is searched at a large range during the early stage of search. As the number of iterations increases, the search range is gradually narrowed, and the algorithm is effectively balanced to escape from local minimum.(7)ms=mmax⋅1−PsPmax,where Ps denotes the probability for the S number of species; S denotes the species number; and mmax denotes the maximum variation rate given; Pmax=maxPs.The specific operation of nonuniform variation is like uniform variation, but it focuses on searching for small areas near the original individual. Assuming that the variation of thekth component Xk of the t generation (t is the current number of iterations) of Habitat Xit=X1,X2,…,Xk,…,Xn is performed, the variation component is as follows:(8)Xk′=Xk+Δt,Umax−Xk,if rand<0.5,(9)Xk′=Xk−Δt,Xk−Umin,if rand≥0.5,where Umax and Umin are the upper and lower limits of component Xk, respectively. Umax−Xk and Xk−Umin represent random numbers in the range of 0,y and meet the nonuniform distribution:(10)limt⟶∞Δt,y=0.This feature allows the mutation operator to search the entire space evenly in the initial stage of the algorithm (whent is very small), and the latter part of the algorithm focuses on several local range exact searches. Δt,y can be defined as follows:(11)Δt,y=y⋅1−r1−t/T⋅b,where r is a random number that fits evenly within the 0,1 range, T is the maximum number of iterations, and b is the system parameter that determines the degree of dependency on the number of iterations, select here b=2.Then, the specific steps of NBBO algorithm are as follows: (Algorithm1)Algorithm 1:The main procedure of NBBO.
BeginAlgorithmic parameters setting:Population size = 50,Maximum number of iterationsG = 500,Maximum immigration rateI=1,Maximum emigration rateE=1,Mutation probabilityPm=0.05,Number of elites retainedKeep=2./∗ Initialization ∗/Random generation of a set of initial habitats constitutes the initial population.Calculate the fitness value of each habitat in the population.While notT /∗T is the condition for the end of the iteration ∗/Sort the habitat in descending order according to the fitness value.The immigration rateλi and emigration rate μi of each habitat were calculated according to formula (1)./∗ End of initialization∗//∗ Migration ∗/piSelectHi according to the immigration probability λiforj=1:D /∗D is dimension ∗/if rand<λiSelectHj according to the migration probability μiRandomly select a characteristic variableSIV from HjReplace a random characteristic variableSIV in Hi with SIVelsePerform a local search on the characteristic variableSIVaccording to formula (6) can get SIV′.Replace random characteristic variableSIV in Hi with SIV′.end ifend forend for/∗ End of migration ∗//∗ Mutation ∗/fori=1:SizeCalculate the probability of mutationPs according to λi and μiSelection of unmutated habitatHi based on Psif rand0,1<msthenif rand<0.5ReplaceHiSIV with a randomly generated SIV by formula (8)if rand≥0.5ReplaceHiSIV with a randomly generated SIV by (9)end ifend ifend ifend for/∗ End of mutation ∗/Recalculate the habitat fitness valuesend whileend


## 2.4. Convergence Analysis

By using Markov chain strategy, we can analyze the convergence of the proposed NBBO algorithms [19, 20]. Let Xt=x1t,x2t,…,xnt be a population with improved algorithm algebra of t. We divide Xt into n subsets based on the fitness of the population. xiti=1,2,…,n is a subset of them; next, we use Markov chain to analyze the global convergence of the improved BBO algorithm.The evolution process of the NBBO algorithm is mainly composed of selection, migration and mutation operations, which are independent of evolution algebra. Considering the number of individuals of BBO algorithm, the evolution process of the NBBO algorithm satisfies the finite homogeneous Markov chain. Population subsetxiti=1,2,…,n can be equivalent to a state on a finite homogeneous Markov chain, pi means the probability of being in state xit, and then pi>0 and ∑i=1npi=1.The transition probabilitypij=pxi⟶xj of state xit is as follows:(12)pij=0,i<jpij>0,i≥j.Then, the state transition matrixP of the Markov chain is(13)P=p110⋯0p21p22⋯0⋮⋮⋱⋮pn1pn2⋯pnn=p110Q2Q1,where(14)Q1=p22⋯0⋮⋱⋮pn2⋯pnn,Q2=p21,p31,⋯,pn1T.It can be found by formulas (12) and (13) Q2≠0, Q1≠0. Since the sum of the probabilities of each row in the Markov transition matrix is 1, p11=1≠0. It can be concluded that the state transition matrix P of Markov chain satisfies the reducible random matrix. According to [21], the definition of a reductive random matrix is as follows:Let matrixZ of order N be a reducible random matrix which can be obtained by the same row transformation and column transformation, namely,(15)Z=C0RT,where C is the mm≤N order primitive matrix and R and T are nonzero matrices.Then, we have(16)Z∞=limk⟶∞Zk=limk⟶∞Ck0∑i=1k−1TiRCk−iTk=C∞0R∞0.The above matrix is a stable random matrix, andZ∞ satisfies the following conditions:(17)Z∞=ZijN∗N,Zij>0,1≤i≤N,1≤j≤m,Zij=0,1≤i≤N,m<j≤N.Therefore, state transition matrixP is a reductive random matrix, then(18)P∞=limt⟶∞p110Q2Q1t=limt⟶∞p11t0∑i=1t−1Q1iQ2p11t−iQ1t=limt⟶∞p11∞0Q2∞0.In formula (12) p11∞=p11=1 and Q1∞=0, since the sum of the probabilities of each line in the Markov transition matrix is 1; there must be Q2∞=1,1,…,1T, then(19)P∞=10⋯010⋯0⋮⋮⋱⋮10⋯0.Formula (19) shows that when the number of iterations t⟶∞, probability pi1∞=1i=1,2,…,n, That is, no matter how the initial state can finally converge to the global optimal solution with probability 1, it can be concluded that(20)limt⟶∞pYt=f∗=1,where Yt is the optimal fitness value of population Xt, f∗ represents the global optimal fitness value, and pYt=f∗ indicates the probability that the optimal individual of the t generation is the global optimal value. It is proved that the NBBO algorithm has global convergence.


## 3. Simulation Experiments

### 3.1. Test Function Simulation

In order to test the optimization performance of the proposed IMBBO algorithm, 13 standard test functions were selected for simulation experiments [22, 23]. Table 1 shows the test function expression, the range of values for each independent variable, and the theoretical optimal value. We select the optimal value, average value, and standard deviation of the test results as the evaluation criteria to test the algorithm performance.Table 1
Test function.
FunctionExpressionRangesTheory optimalf1:Quarticfx1=∑i=1nix4+random0,1(−1.28, 1.28)0f2:Spherefx4=∑i=1nxi2{−100, 100)0f3:Powellfx5=∑i=1d/4x4i−3+10x4i−22+5x4i−1−x4i2+x4i−2−2x4i−14+10x4i−3−x4i4(−4, 5)0f4:Easomfx10=−cosx1cosx2exp−x1−π2−x2−π2(−100, 100)−1f5:Alpinefx6=∑i=1nxisinxi+0.1xi(10, 10)0f6:Griewankfx7=1/4000∑i=1nxi2−∏i=1ncosxi/i+1(−600, 600)0f7:Ackleyfx8=−20exp−0.21/n∑i=1nxi2−exp1/n∑i=1ncos2πxi+20+e(−32, 32)0f8:Kowalikfx13=∑i=111ai−x1bi2+bix2/bi2+bix3+x42(−5, 5)3.075E − 04The experimental function of Table1 is numerically tested using the proposed NBBO algorithm proposed in this paper, and the basic BBO algorithm, grey wolf optimization algorithm [24] (DE), elephant herding optimization [25] (EHO), salp swarm algorithm [26] (DEBBO), and particle swarm optimization biogeography optimization algorithm [27] and (BBOPSO) for comparison of results. The experimental parameter settings of the algorithm are shown in Table 2.Table 2
Algorithm parameter settings.
Setting objectParameterValueCommon partPopulation sizeThe maximum number of iterationsElite reservesSize=50G=500Keep=2BBO/BBOPSOHabitat modification probabilityMaximum immigrationMaximum migrationMutation probabilityPmod=1I=1E=1Pm=0.05DEWeighting factorCrossover probabilityF=0.5CR=0.5EHOThe scale factorsThe number of clansα=0.5, β=0.5nClan=5DEBBOWeighting factorCrossover probabilityDE mutation schemeF=0.5CR=0.9DE/rand/1/binThe proposed NBBOHabitat modification probabilityMaximum immigrationMaximum migrationMutation probabilitySystem parametersPmod=1I=1E=1Pm=0.05b=2To study the performance of the algorithm in different dimensions, the dimensions of test functionf1−f3 and f5−f7 are D=10,30,50, the dimension of test functions f4 is D=2, and the dimension of test functions f8 is D=4. In order to ensure that the comparison between each method is fair, we ensure that the experimental environment of each method is consistent during the experiment: operating system: window 10, memory: 8 GB, programming language: MATLAB 2016, and ensure that the setting of common parameters of the algorithm is consistent, and the parameters of different algorithms are set according to the parameters of specific references. As the algorithm is running randomly, the six algorithms independently run each test function 30 times to avoid unnecessary errors caused by the random operation of the algorithm.Tables3 and 4 show the recorded optimal function value, average value, and standard deviation obtained by the algorithm. In the iterative process of the algorithm, we record the best fitness of each evaluation, find the average optimal value of each evaluation and draw a curve, which reflects the convergence trend of the algorithm. The convergence curves of some test functions are shown in Figures 2 and 3.Table 3
Unimodal function test results.
FunctionDimensionEvaluating indicatorReference [4]BBOReference [24]DEReference [25]EHOReference [26]DEBBOReference [27]BBOPSOThe proposed methodf1xQuarticD = 10BestAve|Std3.40E − 031.37 E − 029.80 E − 031.70E − 034.10 E − 031.70 E − 032.10E − 036.60 E − 033.70 E − 031.81E − 049.45 E − 046.91 E − 047.92E − 051.20 E − 038.62 E − 045.78 E − 054.51 E − 043.25 E − 04D = 30BestAveStd7.58E − 021.47 E − 014.41 E − 023.50E − 025.46 E − 021.03 E − 028.44E − 021.33 E − 014.13 E − 021.05E − 021.46 E − 024.20 E − 036.20E − 032.23 E − 029.10 E − 031.10 E − 032.70 E − 032.00 E − 03D = 50BestAveStd2.51E − 015.67 E − 012.08 E − 011.21E − 012.29 E − 015.16 E − 021.98E − 015.37 E − 011.38 E − 013.47E − 025.61 E − 021.56 E − 024.63E − 027.82 E − 022.22 E − 023.00 E − 038.50 E − 034.20 E − 03f2xSphereD = 10BestAveStd1.94E + 008.38 E + 006.12 E + 001.18E − 203.52 E − 202.77 E − 202.35E − 043.74 E − 046.88 E − 056.72E − 321.69 E − 292.06 E − 292.02E − 281.37 E − 271.49 E − 273.66 E − 455.79 E − 439.78 E − 43D = 30BestAveStd7.13E + 011.43 E + 024.64 E + 011.20E − 032.60 E − 032.20 E − 033.00E − 033.60 E − 030.31 E − 041.40E − 105.83 E − 104.90 E − 106.73E − 071.24 E − 053.35 E − 052.09 E − 241.43 E − 223.91 E − 22D = 50BestAveStd5.93E + 027.27 E + 021.34 E + 021.58E + 012.30 E + 014.63 E + 007.90E − 039.00 E − 031.20 E − 032.38E − 051.62 E − 041.54 E − 044.11E − 021.64 E − 019.32 E − 024.66 E − 186.83 E − 178.48 E − 17f3xPowellD = 10BestAveStd2.47E − 023.01 E − 012.54 E − 012.04E − 051.52 E − 041.04 E − 044.64E − 067.51 E − 062.22 E − 063.89E − 042.40 E − 031.80 E − 032.63E − 041.45 E − 031.24 E − 038.88 E − 135.16 E − 077.21 E − 07D = 30BestAveStd2.47E + 005.77 E + 003.75 E + 006.19E + 011.36 E + 014.70 E + 012.36E − 044.10 E − 041.50 E − 041.21E − 013.92 E − 011.61 E − 015.35E − 022.90 E − 011.25 E − 011.13 E − 101.26 E − 062.10 E − 06D = 50BestAveStd1.73E + 015.16 E + 012.44 E + 011.04E + 031.87 E + 034.83 E + 022.20E − 033.30 E − 036.48 E − 042.02E + 004.55 E + 002.61 E + 001.63E + 002.55 E + 009.34 E − 014.07 E − 094.69 E − 066.33 E − 09f4xEasomD = 2BestAveStd−9.13E − 01–5.00 E − 014.45 E − 01−1.00 E + 00−1.00 E + 000.00 E + 00−9.95E − 01–9.60 E − 013.55EE − 02−1.00 E + 00−1.00 E + 000.00 E + 00−1.00 E + 00−1.00 E + 000.00 E + 00−1.00 E + 00−1.00 E + 000.00 E + 00Table 4
Multimodal function test results.
FunctionDimensionEvaluating indicatorReference [4] BBOReference [24] DEReference [25] EHOReference [26] DEBBOReference [27] BBOPSOThe proposed methodf5xAlpineD = 10BestAveStd2.59E − 025.18 E − 022.26 E − 021.41E − 192.97 E − 165.21 E − 163.59E − 064.66 E − 068.50 E − 074.91E − 068.09 E − 055.62 E − 059.43E − 091.58 E − 051.77 E − 051.26 E − 268.04 E − 251.73 E − 24D = 30BestAveStd3.24E − 017.83 E − 012.94 E − 017.41E + 001.11 E + 011.97 E + 005.22E − 067.12 E − 051.75 E − 054.50E − 036.21 E − 025.85 E − 024.43E − 039.06 E − 033.07 E − 039.64 E − 161.05 E − 133.65 E − 13D = 50BestAveStd2.32E + 003.26 E + 005.16 E − 013.12E + 013.62 E + 012.45 E + 001.68E − 042.29 E − 045.69 E − 053.33E − 011.15 E + 005.32 E − 014.26E − 021.58 E − 011.41 E − 013.55 E − 131.03 E − 102.92 E − 10f6xGriewankD = 10BestAveStd1.00E + 001.07 E + 004.12 E − 021.00E + 001.00 E + 000.00 E + 001.10E + 001.14 E + 002.33 E − 020.00 E + 002.48 E − 021.57 E − 020.00 E + 002.95 E − 022.54 E − 020.00 E + 000.00 E + 000.00 E + 00D = 30BestAveStd1.88E + 002.40 E + 003.36 E − 011.00E + 001.00 E + 007.20 E − 042.90E + 003.61 E + 005.98 E − 016.03E − 101.50 E − 033.20 E − 033.66E − 071.24 E − 022.40 E − 020.00 E + 000.00 E + 000.00 E + 00D = 50BestAveStd4.86E + 008.00 E + 002.26 E + 001.12E + 001.18 E + 003.33 E − 027.72E + 009.82 E + 001.01 E + 004.22E − 051.90 E − 035.00 E − 039.11E − 021.78 E − 016.79 E − 0202.22 E − 177.02 E − 17f7xAckleyD = 10BestAveStd8.70E − 012.05 E + 006.17 E − 013.43E − 118.75 E − 116.33 E − 117.20E − 038.10 E − 036.43 E − 042.66 E − 152.66 E − 150.00 E + 006.22E − 158.85 E − 132.54 E − 122.66 E − 152.66 E − 150.00 E + 00D = 30BestAveStd3.58E + 004.30 E + 004.09 E − 011.28E − 021.63 E − 022.30 E − 031.44E − 021.56 E − 027.98 E − 045.12E − 066.50 E − 062.13 E − 061.08E − 042.66 E − 041.25 E − 041.23 E − 131.23 E − 131.07 E − 12D = 50BestAveStd5.43E + 006.18 E + 005.51 E − 012.14E + 002.57 E + 002.98 E − 011.71E − 021.88 E − 029.30 E − 028.84E − 041.70 E − 031.20 E − 032.74E − 027.18 E − 026.75 E − 021.75 E − 109.85 E − 109.57 E − 10f8xKowalikD = 4BestAveStd6.13E − 042.70 E − 034.30 E − 033.17E − 046.10 E − 044.01 E − 043.24E − 043.52 E − 042.27 E − 053.16E − 049.99 E − 041.20 E − 034.91E − 047.38 E − 043.14 E − 043.08 E − 043.23 E − 042.26 E − 05Figure 2
Unimodal function simulation results diagram. (a)f1:Quarticd=10, (b) f1:Quarticd=30, (c) f1:Quarticd=50, (d) f2:Sphered=10, (e) f2:Sphered=30, (f) f2:Sphered=50, (g) f3:Powelld=10, (h) f3:Powelld=30, (i) f3:Powelld=50, (j) f4:Easomd=2.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)Figure 3
Multimodal function simulation results diagram, (a)f5:Alpined=10, (b) f5:Alpined=30, (c) f5:Alpined=50, (d) f6:Griewankd=10, (e) f6:Griewankd=30, (f) f6:Griewankd=50, (g) f7:Ackleyd=10, (h) f7:Ackleyd=30, (i) f7:Ackleyd=50, (j) f8:Kowalikd=4.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)In order to compare the optimization performance of each algorithm more easily and clearly, the algorithm with the best optimization performance of the optimal value, average value and standard deviation of each test function is boldly marked in Tables3 and 4. From Tables 3 and 4 and the simulation diagram, we find that the convergence speed of the analytical algorithm can be qualitatively and quantitatively analyzed. For the simple multidimensional unimodal function in Figure 2, the NBBO algorithm has faster convergence speed and solution precision in the later iteration although its advantage is not obvious in the initial stage of iteration. For the more complex multidimensional multimodal function in Figure 3, the other five algorithms have the phenomenon of search stagnation or slow search, i.e., they fall into local optimum and are difficult to jump out. Then, by comparing the optimal value of each test function and the average value of the optimization results, we can see that IMBBO algorithm has achieved a smaller value. Although the optimization results of DE, DEBBO, and BBOPSO algorithms also reach the global optimal value in the optimization of some test functions with different dimensions, it can be seen from the optimization results of single peak function and multi peak function that NBBO algorithm has a great advantage in the optimization ability. Figures 2 and 3 show that the proposed NBBO algorithm can not only improve the convergence of the algorithm but also improve the evolution ability and the ability to jump out of the local optimal solution. It is superior to the other four comparison algorithms in the stability of the optimization result. In summary, when solving the optimization problem of single-peak or multi-peak functions, the proposed NBBO can quickly and effectively find the optimal solution under a short number of iterations, and its optimal solution has higher precision and stability.Considering the fixed dimension test function has few dimensions to be adjusted, the accuracy and stability of each algorithm are relatively high. For the multidimensional test function, when its dimensionD is increased from 10 to 30 or 50, the accuracy and stability of the algorithm have decreased. This is because the objective function is more complicated and the algorithm needs to make more adjustments. However, no matter how the dimension changes, the proposed NBBO algorithm still maintains a stable and high-precision optimization performance for high-dimensional complex functions. In reference [28], derrick proposed that for the evaluation of the performance of the improved evolutionary algorithm, only based on the average value and the standard deviation value to compare the advantages and disadvantages of the algorithm is not rigorous enough, so statistical test should be carried out. The statistical test shows that the improved algorithm has significant improvement advantages over other existing algorithms. Therefore, in order to verify the significant difference between the IMBBO algorithm and other algorithms, a nonparametric estimation method, Wilcoxon rank sum test, is used for statistical analysis. For each test function, the 30 times solution results of IMBBO algorithm are tested with those of BBO, DE, EHO, DEBBO, and BBOPSO algorithms respectively, and H0 is assumed: the solution results of the two algorithms are generally the same; H1: the solution results of the two algorithms are generally different. Set the threshold value of significance level of zero hypothesis test α = 0.05. When the test result P<0.05, reject the zero hypothesis; the experimental results of the two algorithms have significant difference; when the test result P>0.05, accept the zero hypothesis; the experimental results of the two algorithms have no significant difference. The specific inspection results are shown in Tables 5 and 6.Table 5
Wilcoxon rank sum test results of unimodal function.
FunctionDimensionNBBO vs BBONBBO vs DENBBO vs EHONBBO vs DEBBONBBO vs BBOPSOPWPWPWPWPWf1xQuarticD = 102.87e − 11+2.87e − 11+2.87e − 11+1.40e − 03+1.68e − 06+D = 302.82e − 11+2.82e − 11+2.82e − 11+2.82e − 11+6.93e − 11+D = 502.87e − 11+2.87e − 11+2.87e − 11+2.82e − 11+2.87e − 11+f2xSphereD = 102.87e − 11+2.87e − 11+2.87e − 11+2.87e − 11+2.78e − 11+D = 302.87e − 11+2.87e − 11+2.82e − 11+2.87e − 11+2.87e − 11+D = 502.87e − 11+2.87e − 11+2.82e − 11+2.87e − 11+2.87e − 11+f3xPowellD = 102.87e − 11+2.87e − 11+2.87e − 11+2.87e − 11+2.87e − 11+D = 302.87e − 11+2.87e − 11+2.93e − 11+2.93e − 11+2.93e − 11+D = 502.87e − 11+2.87e − 11+2.78e − 11+2.87e − 11+2.87e − 11+f4xEasomD = 21.17e − 12+N/AN/A1.17e − 12+N/AN/AN/AN/ATable 6
Wilcoxon rank sum test results of multimodal function.
FunctionDimensionNBBO Vs BBONBBO vs DENBBO vs EHONBBO vs DEBBONBBO vs BBOPSOPWPWPWPWPWf5xAlpineD = 102.87e − 11+2.87e − 11+2.87e − 11+2.87e − 11+2.87e − 11+D = 302.87e − 11+2.87e −  11+2.87e − 11+2.87e − 11+2.87e − 11+D = 502.87e − 11+2.87e − 11+2.87e − 11+2.87e − 11+2.87e − 11+f6xGriewankD = 101.17e − 12+1.69e − 14+1.17e − 12+5.62e − 11+5.62e − 11+D = 301.17e − 12+1.69e − 14+1.17e − 12+1.17e − 12+1.17e − 12+D = 503.06e − 12+3.06e − 12+3.06e − 12+3.06e − 12+3.06e − 12+f7xAckleyD = 101.17e − 12+1.17e − 121.03e − 13+N/AN/A1.15e − 12+D = 302.87e − 11+2.87e − 11+2.87e − 11+2.87e − 11+2.87e − 11+D = 502.87e − 11+2.87e − 11+2.87e − 11+2.87e − 11+2.87e − 11+f8xKowalikD = 22.87e − 11+4.27e − 07+1.68e − 06+2.31e − 08+2.87e − 11+Tables5 and 6 show the test results of test function f1−f8. Among them, for the probability p value less than the significance level threshold value of 0.05, bold display; “W” is the significance test judgment result, combined with Tables 3 and 4, use “+” to indicate that NBBO is more accurate and significant than other algorithms, use “−” to indicate that NBBO is less accurate and significant than other algorithms, and use “ = ” to indicate that NBBO algorithm. There is no significant difference in the accuracy of other algorithms; “N/A” means that significance judgment cannot be made. From the test results, most of the test p values are less than 0.05 and “W” is “+,” so zero hypothesis is rejected. Therefore, there are significant differences between the calculation results of IMBBO and BBO, DE, EHO, DEBBO, and BBOPSO, and NBBO is significantly better.To sum up, the biogeography-based optimization algorithm based on local search and nonuniform variation has good results for most of the benchmark functions, which effectively solves the slow convergence speed of the original BBO algorithm and easy to fall into the local optimization problem in the later stage of the search and greatly improves the convergence speed and search accuracy of the BBO algorithm without adding too many programming steps.


### 3.2. Iris Classification Based on NBBO Multilayer Perceptron

In order to test the training performance of NBBO multilayer perceptron for iris classification, the mean square error (MSE) of the actual output value of multilayer perceptron and the expected output value is used as an indicator to measure the performance of multilayer perceptron [29]. MSE of all training samples is calculated as follows:(21)MSE=1s∑k=1s∑i=1moik−dik2,where s is the number of training samples, m is the number of outputs, dik is the expected output of the i input unit under k training samples, and oik is the actual output of the i input unit under k training samples. The MSE value is smaller is equal to the performance of the multilayer perceptron is better, and vice versa. The HIS of the i habitat is calculated by(22)HSIHabitati=EHabitati.Figure4 shows the schematic diagram of a multilayer perceptron based on NBBO algorithm. NBBO algorithm receives the average MSE of all training samples and all expected output samples as the objective function for iterative training and adjusts the link weights and deviations through continuous iterative evolution, so as to provide the optimal connection weight and deviation for iris classification samples after iterative training.Figure 4
NBBO-based Iris classification diagram.The iris dataset has 150 training samples divided into 3 categories: Setosa, Versicolor, and Virginica, with a total of 4 basic features: sepal length, sepal width, petal length, and petal width. We used the structure of multilayer perceptron as 4-9-3 to classify the problem. Table7 shows the experimental results of the data set algorithm. Figure 5 shows the average convergence trend and classification accuracy.Table 7
Simulation results of iris classification.
AlgorithmOptimal valueAverage valueStandard deviationReference [4]3.15E − 025.34E − 021.52E − 02Reference [24]6.12E − 021.31E − 013.48E − 02Reference [25]2.37E − 013.92E − 017.36E − 02Reference [26]2.73E − 024.72E − 021.29E − 02Reference [27]1.51E − 012.56E − 015.34E − 02The proposed method1.34 E − 021.93 E − 023.80 E − 03Figure 5
(a) Convergence curve and (b) classification accuracy.
(a)(b)Compared with other five iris classification algorithms in Table7, we can see that the proposed NBBO has the better optimal value, average value, and standard deviation than the other five algorithms. The optimal and average values of the MSE calculated by NBBO are better than other algorithms, which indicates that the optimization performance of the NBBO is the best among the seven algorithms. At the same time, the minimum standard deviation calculated by NBBO also indicates that the NBBO algorithm inherits the strong robustness of the BBO algorithm and is more robust than other algorithms.Figure5(a) shows the classification convergence curves of the six algorithms for the iris classification problem, among which the convergence effect of the DE algorithm is the worst of the six algorithms, and the EHO and BBOPSO algorithms fall into the local optimum in the late iteration. BBO algorithm and DEBBO algorithm have similar convergence precision but the convergence speed is not ideal. The proposed NBBO algorithm is optimal in convergence precision and convergence speed. Figure 5(b) shows that DE, EHO, and BBOPSO show significant fluctuations in the classification accuracy of 10 independent operations. The NBBO algorithm has the best stability and average, and its classification accuracy is the highest among the three stable algorithms BBO, DEBBO, and NBBO.


## 3.1. Test Function Simulation

In order to test the optimization performance of the proposed IMBBO algorithm, 13 standard test functions were selected for simulation experiments [22, 23]. Table 1 shows the test function expression, the range of values for each independent variable, and the theoretical optimal value. We select the optimal value, average value, and standard deviation of the test results as the evaluation criteria to test the algorithm performance.Table 1
Test function.
FunctionExpressionRangesTheory optimalf1:Quarticfx1=∑i=1nix4+random0,1(−1.28, 1.28)0f2:Spherefx4=∑i=1nxi2{−100, 100)0f3:Powellfx5=∑i=1d/4x4i−3+10x4i−22+5x4i−1−x4i2+x4i−2−2x4i−14+10x4i−3−x4i4(−4, 5)0f4:Easomfx10=−cosx1cosx2exp−x1−π2−x2−π2(−100, 100)−1f5:Alpinefx6=∑i=1nxisinxi+0.1xi(10, 10)0f6:Griewankfx7=1/4000∑i=1nxi2−∏i=1ncosxi/i+1(−600, 600)0f7:Ackleyfx8=−20exp−0.21/n∑i=1nxi2−exp1/n∑i=1ncos2πxi+20+e(−32, 32)0f8:Kowalikfx13=∑i=111ai−x1bi2+bix2/bi2+bix3+x42(−5, 5)3.075E − 04The experimental function of Table1 is numerically tested using the proposed NBBO algorithm proposed in this paper, and the basic BBO algorithm, grey wolf optimization algorithm [24] (DE), elephant herding optimization [25] (EHO), salp swarm algorithm [26] (DEBBO), and particle swarm optimization biogeography optimization algorithm [27] and (BBOPSO) for comparison of results. The experimental parameter settings of the algorithm are shown in Table 2.Table 2
Algorithm parameter settings.
Setting objectParameterValueCommon partPopulation sizeThe maximum number of iterationsElite reservesSize=50G=500Keep=2BBO/BBOPSOHabitat modification probabilityMaximum immigrationMaximum migrationMutation probabilityPmod=1I=1E=1Pm=0.05DEWeighting factorCrossover probabilityF=0.5CR=0.5EHOThe scale factorsThe number of clansα=0.5, β=0.5nClan=5DEBBOWeighting factorCrossover probabilityDE mutation schemeF=0.5CR=0.9DE/rand/1/binThe proposed NBBOHabitat modification probabilityMaximum immigrationMaximum migrationMutation probabilitySystem parametersPmod=1I=1E=1Pm=0.05b=2To study the performance of the algorithm in different dimensions, the dimensions of test functionf1−f3 and f5−f7 are D=10,30,50, the dimension of test functions f4 is D=2, and the dimension of test functions f8 is D=4. In order to ensure that the comparison between each method is fair, we ensure that the experimental environment of each method is consistent during the experiment: operating system: window 10, memory: 8 GB, programming language: MATLAB 2016, and ensure that the setting of common parameters of the algorithm is consistent, and the parameters of different algorithms are set according to the parameters of specific references. As the algorithm is running randomly, the six algorithms independently run each test function 30 times to avoid unnecessary errors caused by the random operation of the algorithm.Tables3 and 4 show the recorded optimal function value, average value, and standard deviation obtained by the algorithm. In the iterative process of the algorithm, we record the best fitness of each evaluation, find the average optimal value of each evaluation and draw a curve, which reflects the convergence trend of the algorithm. The convergence curves of some test functions are shown in Figures 2 and 3.Table 3
Unimodal function test results.
FunctionDimensionEvaluating indicatorReference [4]BBOReference [24]DEReference [25]EHOReference [26]DEBBOReference [27]BBOPSOThe proposed methodf1xQuarticD = 10BestAve|Std3.40E − 031.37 E − 029.80 E − 031.70E − 034.10 E − 031.70 E − 032.10E − 036.60 E − 033.70 E − 031.81E − 049.45 E − 046.91 E − 047.92E − 051.20 E − 038.62 E − 045.78 E − 054.51 E − 043.25 E − 04D = 30BestAveStd7.58E − 021.47 E − 014.41 E − 023.50E − 025.46 E − 021.03 E − 028.44E − 021.33 E − 014.13 E − 021.05E − 021.46 E − 024.20 E − 036.20E − 032.23 E − 029.10 E − 031.10 E − 032.70 E − 032.00 E − 03D = 50BestAveStd2.51E − 015.67 E − 012.08 E − 011.21E − 012.29 E − 015.16 E − 021.98E − 015.37 E − 011.38 E − 013.47E − 025.61 E − 021.56 E − 024.63E − 027.82 E − 022.22 E − 023.00 E − 038.50 E − 034.20 E − 03f2xSphereD = 10BestAveStd1.94E + 008.38 E + 006.12 E + 001.18E − 203.52 E − 202.77 E − 202.35E − 043.74 E − 046.88 E − 056.72E − 321.69 E − 292.06 E − 292.02E − 281.37 E − 271.49 E − 273.66 E − 455.79 E − 439.78 E − 43D = 30BestAveStd7.13E + 011.43 E + 024.64 E + 011.20E − 032.60 E − 032.20 E − 033.00E − 033.60 E − 030.31 E − 041.40E − 105.83 E − 104.90 E − 106.73E − 071.24 E − 053.35 E − 052.09 E − 241.43 E − 223.91 E − 22D = 50BestAveStd5.93E + 027.27 E + 021.34 E + 021.58E + 012.30 E + 014.63 E + 007.90E − 039.00 E − 031.20 E − 032.38E − 051.62 E − 041.54 E − 044.11E − 021.64 E − 019.32 E − 024.66 E − 186.83 E − 178.48 E − 17f3xPowellD = 10BestAveStd2.47E − 023.01 E − 012.54 E − 012.04E − 051.52 E − 041.04 E − 044.64E − 067.51 E − 062.22 E − 063.89E − 042.40 E − 031.80 E − 032.63E − 041.45 E − 031.24 E − 038.88 E − 135.16 E − 077.21 E − 07D = 30BestAveStd2.47E + 005.77 E + 003.75 E + 006.19E + 011.36 E + 014.70 E + 012.36E − 044.10 E − 041.50 E − 041.21E − 013.92 E − 011.61 E − 015.35E − 022.90 E − 011.25 E − 011.13 E − 101.26 E − 062.10 E − 06D = 50BestAveStd1.73E + 015.16 E + 012.44 E + 011.04E + 031.87 E + 034.83 E + 022.20E − 033.30 E − 036.48 E − 042.02E + 004.55 E + 002.61 E + 001.63E + 002.55 E + 009.34 E − 014.07 E − 094.69 E − 066.33 E − 09f4xEasomD = 2BestAveStd−9.13E − 01–5.00 E − 014.45 E − 01−1.00 E + 00−1.00 E + 000.00 E + 00−9.95E − 01–9.60 E − 013.55EE − 02−1.00 E + 00−1.00 E + 000.00 E + 00−1.00 E + 00−1.00 E + 000.00 E + 00−1.00 E + 00−1.00 E + 000.00 E + 00Table 4
Multimodal function test results.
FunctionDimensionEvaluating indicatorReference [4] BBOReference [24] DEReference [25] EHOReference [26] DEBBOReference [27] BBOPSOThe proposed methodf5xAlpineD = 10BestAveStd2.59E − 025.18 E − 022.26 E − 021.41E − 192.97 E − 165.21 E − 163.59E − 064.66 E − 068.50 E − 074.91E − 068.09 E − 055.62 E − 059.43E − 091.58 E − 051.77 E − 051.26 E − 268.04 E − 251.73 E − 24D = 30BestAveStd3.24E − 017.83 E − 012.94 E − 017.41E + 001.11 E + 011.97 E + 005.22E − 067.12 E − 051.75 E − 054.50E − 036.21 E − 025.85 E − 024.43E − 039.06 E − 033.07 E − 039.64 E − 161.05 E − 133.65 E − 13D = 50BestAveStd2.32E + 003.26 E + 005.16 E − 013.12E + 013.62 E + 012.45 E + 001.68E − 042.29 E − 045.69 E − 053.33E − 011.15 E + 005.32 E − 014.26E − 021.58 E − 011.41 E − 013.55 E − 131.03 E − 102.92 E − 10f6xGriewankD = 10BestAveStd1.00E + 001.07 E + 004.12 E − 021.00E + 001.00 E + 000.00 E + 001.10E + 001.14 E + 002.33 E − 020.00 E + 002.48 E − 021.57 E − 020.00 E + 002.95 E − 022.54 E − 020.00 E + 000.00 E + 000.00 E + 00D = 30BestAveStd1.88E + 002.40 E + 003.36 E − 011.00E + 001.00 E + 007.20 E − 042.90E + 003.61 E + 005.98 E − 016.03E − 101.50 E − 033.20 E − 033.66E − 071.24 E − 022.40 E − 020.00 E + 000.00 E + 000.00 E + 00D = 50BestAveStd4.86E + 008.00 E + 002.26 E + 001.12E + 001.18 E + 003.33 E − 027.72E + 009.82 E + 001.01 E + 004.22E − 051.90 E − 035.00 E − 039.11E − 021.78 E − 016.79 E − 0202.22 E − 177.02 E − 17f7xAckleyD = 10BestAveStd8.70E − 012.05 E + 006.17 E − 013.43E − 118.75 E − 116.33 E − 117.20E − 038.10 E − 036.43 E − 042.66 E − 152.66 E − 150.00 E + 006.22E − 158.85 E − 132.54 E − 122.66 E − 152.66 E − 150.00 E + 00D = 30BestAveStd3.58E + 004.30 E + 004.09 E − 011.28E − 021.63 E − 022.30 E − 031.44E − 021.56 E − 027.98 E − 045.12E − 066.50 E − 062.13 E − 061.08E − 042.66 E − 041.25 E − 041.23 E − 131.23 E − 131.07 E − 12D = 50BestAveStd5.43E + 006.18 E + 005.51 E − 012.14E + 002.57 E + 002.98 E − 011.71E − 021.88 E − 029.30 E − 028.84E − 041.70 E − 031.20 E − 032.74E − 027.18 E − 026.75 E − 021.75 E − 109.85 E − 109.57 E − 10f8xKowalikD = 4BestAveStd6.13E − 042.70 E − 034.30 E − 033.17E − 046.10 E − 044.01 E − 043.24E − 043.52 E − 042.27 E − 053.16E − 049.99 E − 041.20 E − 034.91E − 047.38 E − 043.14 E − 043.08 E − 043.23 E − 042.26 E − 05Figure 2
Unimodal function simulation results diagram. (a)f1:Quarticd=10, (b) f1:Quarticd=30, (c) f1:Quarticd=50, (d) f2:Sphered=10, (e) f2:Sphered=30, (f) f2:Sphered=50, (g) f3:Powelld=10, (h) f3:Powelld=30, (i) f3:Powelld=50, (j) f4:Easomd=2.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)Figure 3
Multimodal function simulation results diagram, (a)f5:Alpined=10, (b) f5:Alpined=30, (c) f5:Alpined=50, (d) f6:Griewankd=10, (e) f6:Griewankd=30, (f) f6:Griewankd=50, (g) f7:Ackleyd=10, (h) f7:Ackleyd=30, (i) f7:Ackleyd=50, (j) f8:Kowalikd=4.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)In order to compare the optimization performance of each algorithm more easily and clearly, the algorithm with the best optimization performance of the optimal value, average value and standard deviation of each test function is boldly marked in Tables3 and 4. From Tables 3 and 4 and the simulation diagram, we find that the convergence speed of the analytical algorithm can be qualitatively and quantitatively analyzed. For the simple multidimensional unimodal function in Figure 2, the NBBO algorithm has faster convergence speed and solution precision in the later iteration although its advantage is not obvious in the initial stage of iteration. For the more complex multidimensional multimodal function in Figure 3, the other five algorithms have the phenomenon of search stagnation or slow search, i.e., they fall into local optimum and are difficult to jump out. Then, by comparing the optimal value of each test function and the average value of the optimization results, we can see that IMBBO algorithm has achieved a smaller value. Although the optimization results of DE, DEBBO, and BBOPSO algorithms also reach the global optimal value in the optimization of some test functions with different dimensions, it can be seen from the optimization results of single peak function and multi peak function that NBBO algorithm has a great advantage in the optimization ability. Figures 2 and 3 show that the proposed NBBO algorithm can not only improve the convergence of the algorithm but also improve the evolution ability and the ability to jump out of the local optimal solution. It is superior to the other four comparison algorithms in the stability of the optimization result. In summary, when solving the optimization problem of single-peak or multi-peak functions, the proposed NBBO can quickly and effectively find the optimal solution under a short number of iterations, and its optimal solution has higher precision and stability.Considering the fixed dimension test function has few dimensions to be adjusted, the accuracy and stability of each algorithm are relatively high. For the multidimensional test function, when its dimensionD is increased from 10 to 30 or 50, the accuracy and stability of the algorithm have decreased. This is because the objective function is more complicated and the algorithm needs to make more adjustments. However, no matter how the dimension changes, the proposed NBBO algorithm still maintains a stable and high-precision optimization performance for high-dimensional complex functions. In reference [28], derrick proposed that for the evaluation of the performance of the improved evolutionary algorithm, only based on the average value and the standard deviation value to compare the advantages and disadvantages of the algorithm is not rigorous enough, so statistical test should be carried out. The statistical test shows that the improved algorithm has significant improvement advantages over other existing algorithms. Therefore, in order to verify the significant difference between the IMBBO algorithm and other algorithms, a nonparametric estimation method, Wilcoxon rank sum test, is used for statistical analysis. For each test function, the 30 times solution results of IMBBO algorithm are tested with those of BBO, DE, EHO, DEBBO, and BBOPSO algorithms respectively, and H0 is assumed: the solution results of the two algorithms are generally the same; H1: the solution results of the two algorithms are generally different. Set the threshold value of significance level of zero hypothesis test α = 0.05. When the test result P<0.05, reject the zero hypothesis; the experimental results of the two algorithms have significant difference; when the test result P>0.05, accept the zero hypothesis; the experimental results of the two algorithms have no significant difference. The specific inspection results are shown in Tables 5 and 6.Table 5
Wilcoxon rank sum test results of unimodal function.
FunctionDimensionNBBO vs BBONBBO vs DENBBO vs EHONBBO vs DEBBONBBO vs BBOPSOPWPWPWPWPWf1xQuarticD = 102.87e − 11+2.87e − 11+2.87e − 11+1.40e − 03+1.68e − 06+D = 302.82e − 11+2.82e − 11+2.82e − 11+2.82e − 11+6.93e − 11+D = 502.87e − 11+2.87e − 11+2.87e − 11+2.82e − 11+2.87e − 11+f2xSphereD = 102.87e − 11+2.87e − 11+2.87e − 11+2.87e − 11+2.78e − 11+D = 302.87e − 11+2.87e − 11+2.82e − 11+2.87e − 11+2.87e − 11+D = 502.87e − 11+2.87e − 11+2.82e − 11+2.87e − 11+2.87e − 11+f3xPowellD = 102.87e − 11+2.87e − 11+2.87e − 11+2.87e − 11+2.87e − 11+D = 302.87e − 11+2.87e − 11+2.93e − 11+2.93e − 11+2.93e − 11+D = 502.87e − 11+2.87e − 11+2.78e − 11+2.87e − 11+2.87e − 11+f4xEasomD = 21.17e − 12+N/AN/A1.17e − 12+N/AN/AN/AN/ATable 6
Wilcoxon rank sum test results of multimodal function.
FunctionDimensionNBBO Vs BBONBBO vs DENBBO vs EHONBBO vs DEBBONBBO vs BBOPSOPWPWPWPWPWf5xAlpineD = 102.87e − 11+2.87e − 11+2.87e − 11+2.87e − 11+2.87e − 11+D = 302.87e − 11+2.87e −  11+2.87e − 11+2.87e − 11+2.87e − 11+D = 502.87e − 11+2.87e − 11+2.87e − 11+2.87e − 11+2.87e − 11+f6xGriewankD = 101.17e − 12+1.69e − 14+1.17e − 12+5.62e − 11+5.62e − 11+D = 301.17e − 12+1.69e − 14+1.17e − 12+1.17e − 12+1.17e − 12+D = 503.06e − 12+3.06e − 12+3.06e − 12+3.06e − 12+3.06e − 12+f7xAckleyD = 101.17e − 12+1.17e − 121.03e − 13+N/AN/A1.15e − 12+D = 302.87e − 11+2.87e − 11+2.87e − 11+2.87e − 11+2.87e − 11+D = 502.87e − 11+2.87e − 11+2.87e − 11+2.87e − 11+2.87e − 11+f8xKowalikD = 22.87e − 11+4.27e − 07+1.68e − 06+2.31e − 08+2.87e − 11+Tables5 and 6 show the test results of test function f1−f8. Among them, for the probability p value less than the significance level threshold value of 0.05, bold display; “W” is the significance test judgment result, combined with Tables 3 and 4, use “+” to indicate that NBBO is more accurate and significant than other algorithms, use “−” to indicate that NBBO is less accurate and significant than other algorithms, and use “ = ” to indicate that NBBO algorithm. There is no significant difference in the accuracy of other algorithms; “N/A” means that significance judgment cannot be made. From the test results, most of the test p values are less than 0.05 and “W” is “+,” so zero hypothesis is rejected. Therefore, there are significant differences between the calculation results of IMBBO and BBO, DE, EHO, DEBBO, and BBOPSO, and NBBO is significantly better.To sum up, the biogeography-based optimization algorithm based on local search and nonuniform variation has good results for most of the benchmark functions, which effectively solves the slow convergence speed of the original BBO algorithm and easy to fall into the local optimization problem in the later stage of the search and greatly improves the convergence speed and search accuracy of the BBO algorithm without adding too many programming steps.


## 3.2. Iris Classification Based on NBBO Multilayer Perceptron

In order to test the training performance of NBBO multilayer perceptron for iris classification, the mean square error (MSE) of the actual output value of multilayer perceptron and the expected output value is used as an indicator to measure the performance of multilayer perceptron [29]. MSE of all training samples is calculated as follows:(21)MSE=1s∑k=1s∑i=1moik−dik2,where s is the number of training samples, m is the number of outputs, dik is the expected output of the i input unit under k training samples, and oik is the actual output of the i input unit under k training samples. The MSE value is smaller is equal to the performance of the multilayer perceptron is better, and vice versa. The HIS of the i habitat is calculated by(22)HSIHabitati=EHabitati.Figure4 shows the schematic diagram of a multilayer perceptron based on NBBO algorithm. NBBO algorithm receives the average MSE of all training samples and all expected output samples as the objective function for iterative training and adjusts the link weights and deviations through continuous iterative evolution, so as to provide the optimal connection weight and deviation for iris classification samples after iterative training.Figure 4
NBBO-based Iris classification diagram.The iris dataset has 150 training samples divided into 3 categories: Setosa, Versicolor, and Virginica, with a total of 4 basic features: sepal length, sepal width, petal length, and petal width. We used the structure of multilayer perceptron as 4-9-3 to classify the problem. Table7 shows the experimental results of the data set algorithm. Figure 5 shows the average convergence trend and classification accuracy.Table 7
Simulation results of iris classification.
AlgorithmOptimal valueAverage valueStandard deviationReference [4]3.15E − 025.34E − 021.52E − 02Reference [24]6.12E − 021.31E − 013.48E − 02Reference [25]2.37E − 013.92E − 017.36E − 02Reference [26]2.73E − 024.72E − 021.29E − 02Reference [27]1.51E − 012.56E − 015.34E − 02The proposed method1.34 E − 021.93 E − 023.80 E − 03Figure 5
(a) Convergence curve and (b) classification accuracy.
(a)(b)Compared with other five iris classification algorithms in Table7, we can see that the proposed NBBO has the better optimal value, average value, and standard deviation than the other five algorithms. The optimal and average values of the MSE calculated by NBBO are better than other algorithms, which indicates that the optimization performance of the NBBO is the best among the seven algorithms. At the same time, the minimum standard deviation calculated by NBBO also indicates that the NBBO algorithm inherits the strong robustness of the BBO algorithm and is more robust than other algorithms.Figure5(a) shows the classification convergence curves of the six algorithms for the iris classification problem, among which the convergence effect of the DE algorithm is the worst of the six algorithms, and the EHO and BBOPSO algorithms fall into the local optimum in the late iteration. BBO algorithm and DEBBO algorithm have similar convergence precision but the convergence speed is not ideal. The proposed NBBO algorithm is optimal in convergence precision and convergence speed. Figure 5(b) shows that DE, EHO, and BBOPSO show significant fluctuations in the classification accuracy of 10 independent operations. The NBBO algorithm has the best stability and average, and its classification accuracy is the highest among the three stable algorithms BBO, DEBBO, and NBBO.


## 4. Conclusions and Future Work

This paper proposed a novel biogeography-based optimization algorithm based on local migration strategy and nonuniform mutation for iris classification. Firstly, the hyperbolic cotangent transfer model was used to replace the original transfer model to enhance the adaptability of the algorithm. Secondly, the local transfer strategy was used to improve the global search ability of the algorithm. Finally, the nonuniform mutation operation was used to enhance the development ability of the algorithm in the later stage of iteration for the optimal solution in a small range, to improve the ability of the algorithm to exit the local optimum. The simulation results showed that the proposed NBBO has higher stability, better convergence accuracy, and faster convergence speed than other method for iris classification. How to use NBBO algorithm to solve other multiobjective optimization problem compare with multiobjective optimization method [30–33] and realize better applications for engineering practice will be the focus of our next work.



---

*Source: 6694695-2021-04-15.xml*