# GIS-Based Niche Hybrid Bat Algorithm for Solving Optimal Spatial Search

Authors: Guoming Du; Yangbo Chen; Wei Sun
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

Complex nonlinear optimization problems are involved in optimal spatial search, such as location allocation problems that occur in multidimensional geographic space. Such search problems are generally difficult to solve by using traditional methods. The bat algorithm (BA) is an effective method for solving optimization problems. However, the solution of the standard BA is easily trapped at one of its local optimum values. The main cause of premature convergence is the loss of diversity in the population. The niche technique is an effective method to maintain the population diversity, to enhance the exploration of the new search domains, and to avoid premature convergence. In this paper, a geographic information system- (GIS-) based niche hybrid bat algorithm (NHBA) is proposed for solving the optimal spatial search. The NHBA is able to avoid the premature convergence and obtain the global optimal values. The GIS technique provides robust support for processing a substantial amount of geographical data. A case in Fangcun District, Guangzhou City, China, is used to test the NHBA. The comparative experiments illustrate that the BA, GA, FA, PSO, and NHBA algorithms outperform the brute-force algorithm in terms of computational efficiency, and the optimal solutions are more easily obtained with NHBA than with BA, GA, FA, and PSO. Moreover, the precision of NHBA is higher and the convergence of NHBA is faster than those of the other algorithms under the same conditions.


---

## Body

## 1. Introduction

The general goal of the optimal spatial search is to find a set of contiguous places that meet specific optimization objectives, such as minimizing the total cost and maximizing the proximity to certain facilities [1]. Optimal spatial search problems often exhibit considerable complexity, especially when they are involved in large, multiple objective location allocation problems, many of which contain combinatorial optimization under complicated, highly nonlinear constraints, and thus, they are often NP-hard [2, 3]. The traditional solutions include the fixed charge location model, the covering model in noncompetitive location theory, and the medianoid and centroid models from competitive location theory. Because of the exponential increase in the search space with problem size, the combinatorial optimal problems where the reliability goal is achieved by discrete choices made from available parts encounter computational difficulties. Therefore, the methods are difficult to run in the acceptable computational time. Consequently, many research efforts have been devoted to improving or finding efficient algorithms. Metaheuristic algorithms, derived from the behavior of biological or physical systems found in nature, have become powerful methods to solve the tough spatial optimization problems [4]. Inspired by this idea, the genetic algorithm (GA) [5], ant colony optimization (ACO) [6], particle swarm optimization (PSO) [7], artificial fish swarm algorithm (AFSA) [8], cuckoo search (CS) [9], monkey algorithm (MA) [10], and firefly algorithm (FA) [11] have been proposed and applied widely. Every algorithm has its own advantages and disadvantages and is also suited for solving different problems. The bat algorithm (BA) is also a metaheuristic algorithm inspired by a property known as echolocation, which is a type of sonar that guides bats in their flying and hunting behavior [12]. Not only can the bats move, but also they can distinguish different types of insects even in complete darkness, thanks to such capabilities [12, 13]. The BA can solve optimization problems completely stronger than the PSO, GA, and harmony search algorithms [14]. The main reason is that the BA uses an appropriate combination of significant advantages of the mentioned algorithms.However, the solution of the BA is sometimes trapped at one of its local optima because of its premature convergence and weak exploitation capabilities. Therefore, it cannot always obtain the ideal results when BA is used alone [15]. The main cause of premature convergence is the loss of diversity in the population. The niche techniques are regarded as an effective method to maintain the population diversity to enhance the exploration of the new search domains [16]. The niche techniques have been proven particularly useful in problems that require multiple solutions to be found by a search algorithm, such as multimodal and multiobjective optimization problems [17]. As a consequence, the aim of this study is to propose the niche hybrid bat algorithm (NHBA) to enhance the performance of BA on nonlinear optimization problems, avoid becoming trapped in a locally optimal value with the subsequent premature convergence problem of BA, and solve the complex problems of optimal spatial search more efficiently, precisely, and reliably.The main contributions of this study are as follows:(1)
It provides a robust solution by hybridizing the niche techniques with BA in combination with GIS techniques(2)
It contributes to maintaining the diversity of the population properly, improving the global search ability, and overcoming the deficiency in the premature convergence on only one solution(3)
It provides an effective method for solving optimal spatial search and multiobjective nonlinear optimization problemsThe rest of the paper is structured as follows. Section2 gives the literature review of the related works. Section 3 introduces the basic concepts, theory, and procedure of BA. Section 4 proposes the GIS-based niche hybrid bat algorithm (NHBA). Section 5 analyzes and discusses the experimental results by implementing comparative studies. Section 6 concludes the paper.


## 2. Literature Review

Optimal spatial search is involved in combinatorial optimization problems such as location allocation. The genetic algorithm (GA), which searches over regions, was first proposed to identify a final, feasible, optimal, or near-optimal solution to a relaxed version of the redundancy allocation problem [18]. Because GA involves extensive coding and decoding, it will take an excessive amount of time to compute. Simulated annealing was used to solve high-dimensional nonlinear optimization problems for multisite land use allocation (MLUA) problems [2]. Simulated annealing can almost guarantee finding the optimal solution when the cooling process is slow enough, and the simulation is running long enough. However, the fine adjustment of parameters affects the convergence rate of the optimization process. The multiobjective harmony search algorithm was used to optimize the multiserver location allocation problem in congested systems [19]. A modified fruit fly optimization algorithm (MFOA) was proposed to solve a location allocation-inventory problem in a two-echelon supply chain network [20]. Ant colony optimization (ACO) with stochastic demand was used to solve the stochastic uncapacitated location allocation problem with an unknown number of facilities [21]. A metaheuristic algorithm based on the harmony search algorithm was presented for the location, allocation, and routing of temporary health centers in rural areas in a crisis situation [22]. An improved artificial bee colony was proposed for the facility location allocation problem of the end-of-life vehicles recovery network [23]. A self-learning particle swarm optimization (SLPSO) algorithm was developed to solve the multiechelon capacitated location-allocation-inventory problem [24]. It is difficult to avoid becoming stuck at the local minimum value when using the single algorithm. Furthermore, these algorithms lack the support of GIS techniques.The BA was proposed by Yang in 2010. There are many variants of BA in the related literature. A multiobjective bat-inspired algorithm was proposed for solving nonlinear global optimization problems [14]. Simulation results suggested that the proposed algorithm worked efficiently. Yang’s original work was expanded to solve complicated constrained nonlinear optimization problems [25]. A modified optimization methodology called opposition-based BA (OBA) was proposed for the infinite impulse response system identification problem [26]. The simulation results revealed OBA to be a more competent candidate than other evolutionary algorithms (EA), such as real-coded GA, differential evolution (DE), and PSO, in terms of accuracy and convergence speed. The discrete real BA (RBA) was proposed for route search optimization of a graph-based road network [27]. Compared with ACO and the intelligent water drops (IWD) algorithm, RBA achieved a better convergence rate. A novel BA was developed in order to implement the loading pattern optimization of a nuclear reactor core [28]. The test results showed that BA was a very promising algorithm for LPO problems and has the potential to be used in other nuclear engineering optimization problems. A new modified BA was proposed to solve the optimal management of the multiobjective reconfiguration problem [29]. This algorithm was tested on the 32-bus IEEE radial distribution system and demonstrated its feasibility and effectiveness. The bat-inspired algorithm was combined with differential evolution (BA-DE) into a new hybrid algorithm to solve constrained optimization problems [15]. Comparisons showed that the BA-DE outperformed the most advanced methods in terms of the final solution’s quality. A new multiobjective optimization based on the modified BA and Pareto front was developed for solving passive power filter design problems [30]. Local and global search characteristics of BA (EBA) through three different methods were enhanced for solving optimization problems [31]. It was proven that EBA was more effective than other algorithms such as the DE, EA, FA, artificial immune system, and charged system search. BA was used to solve the reliability redundancy allocation problem (RAP) and showed that the BA was competitive with the best known heuristics for redundancy allocation [4]. BA was used to design and optimize the concurrent tolerance in mechanical assemblies [32]. It was found that the BA produced better results than other methods (such as EA and GA) in the initial generations of concurrent tolerance problems. The multidirectional bat algorithm (MDBAT), which hybridized the bat algorithm with the multidirectional search algorithm, was proposed for solving unconstrained global optimization problems [33]. It could accelerate the convergence to the region of the optimal response. The multiobjective bat algorithm was proposed for association rule mining [34]. A novel multiobjective BA was proposed for community detection on dynamic social networks [35]. The algorithm used the mean shift algorithm to generate the new solutions and avoid the random process by defining a new mutation operator. A combination of bat and scaled conjugate gradient algorithms is proposed to improve neural network learning capability [36]. The bat algorithm was first practically implemented in swarm robotics [37]. An autonomous binary version of the bat algorithm was used to solve the patient bed assignment problem [38]. The swarm bat algorithm with improved search was introduced in [39]. The adaptive multiswarm bat algorithm (AMBA) was superior to the others over 20 benchmark functions in a comparison with six algorithms [40]. A hybrid bat algorithm with a genetic crossover operation and smart inertia weight (SGA-BA) was used to solve multilevel image segmentation [41]. A novel bat algorithm with double mutation operators was applied to the low-velocity impact localization problem [42].


## 3. Standard Bat Algorithm

The standard BA is a bioinspired algorithm based on the echolocation behavior of bats [12]. Bats emit a very loud sound pulse and listen for the echo that bounces back from the surroundings. Bats fly randomly using frequency, velocity, and position to search for prey [43, 44]. In the BA, the frequency, velocity, and position of each bat in the population are updated for further movements. This algorithm is formulated to imitate the behavior of bats finding their prey such that it serves to solve both single objective and multiobjective optimization problems in the continuous solution domain [14]. The implementation of the standard BA is as follows [12–14]:(1)
All bats use echolocation to sense distance, and they also “know” the difference between prey and background barriers in an uncanny way.(2)
With a varied frequencyfi, varied wavelength λ, and loudness A0, bats fly randomly with velocity vi at position xi to search for prey. They can automatically adjust the wavelength (or frequency) according to the distance from the prey.(3)
Although the loudness can vary in different ways, it is assumed that the loudness varies from largeA0 to a minimum value Amin.In the BA, each bat is defined by its positionxit, velocity vit, frequency fi, loudness Ait, and the emission pulse rate rit at time t in a D-dimensional search space. The new solution at time t is given by(1)xit=xit−1+vit.The velocity is updated by(2)vit=vit−1+xit−1−xgbesttfi,where xgbestt is the current global best location found by all the bats in the past generations. The pulse frequency fi is dominated by(3)fi=fmin+fmax−fminβ,where S is a random number drawn from a uniform distribution. The two parameters fmin and fmax are the minimum and maximum frequency values, respectively.Whenever a solution is selected as the current optimal solution (xgbest) in the local search areas, the new solution (xi) is generated by a random walk as follows:(4)xi=xgbest+rand∗At,where rand is a random number uniformly distributed within [−1, 1] and controls the direction and power of the random walk. At is the average loudness of all the bats.With theith bat approaching the prey, its loudness Ai will decrease and its rate ri of pulse emission will increase. Such a change can be achieved by the following formulas:(5)Ait+1=αAit,(6)rit+1=ri01−exp−γt,where γ and γ are constants. For any 0<α<1 and γ>0, we have(7)Ait+1⟶0,rit⟶ri0,ast⟶∞.The pseudocode of the original BA is given in Algorithm1Algorithm 1: The pseudocode of the original BA.
Objective functionF(x), x = (x1,x2,...,xd)T
Initialize the bat populationxi (i = 1,2, …, n) and vi
Define the pulse frequencyfi at xi
Initialize the pulse ratesri and the loudness Ai
While (t < Imax) (where Imax is the max number of iterations)
Generate new solutions by adjusting the frequency, and updating velocities and locations using equations (1)–(3)
If (rand1 > ri)
Select a solution among the best solutions
Generate a local solution around the selected best solution using equation (4)
End if
Generate a new solution by flying randomly
If (rand2 < Ai and f(xi) < f(xgbest))
Accept the new solutions
Increaseri and reduce Ai using equations (5) and (6)
End if
Rank the bats and find the current bestxgbest
End while


## 4. GIS-Based Niche Hybrid Bat Algorithm (NHBA)

### 4.1. Brief Description of NHBA

The standard BA has a premature convergence rate and is easily trapped in a local optimum. To overcome these problems, NHBA is proposed by means of the niche technique hybridizing BA. The idea of NHBA is as follows.First, the Euclidean distancedij is calculated between the arbitrary two bats Bati and Batj among the population. If the Euclidean distance dij, denoted as || Bati − Batj ||, is less than the niche distance L which is a threshold value and is given previously based on experience, it signifies that the bats Bi and Bj have higher similarity. The fitness values of two individuals are compared. Then, the lower fitness value is given a penalty factor in order to decrease its fitness value. Regarding the two individuals within the niche distance L, the individuals with lower fitness values will become worse and be eliminated by a large probability in the next generations. It means that there exist only elite individuals within the niche distance L. Therefore, the population diversity is maintained. Meanwhile, the two different individuals maintain a certain distance, and all the individuals are scattered over the whole constrained space. Finally, the niche hybrid bat algorithm is implemented in this way.


### 4.2. GIS-Based NHBA for Optimal Spatial Search

The GIS-based NHBA has the advantage of being applied to optimal spatial search with combined GIS techniques during a sequence of iterations. Optimal spatial search is involved in the resource or location allocation of facilities, which include hospitals, supermarkets, schools, cinemas, and fire stations. Hence, optimal spatial search plays an important role in geographic information science. This paper uses the optimal location allocation of facilities as an example. This type of problem is involved in complex nonlinear optimization scenarios, especially in multiobjective location allocation. If we use the brute-force algorithm, i.e., list all the combinatorial results to obtain the optimal value, the computational volume is so enormous that it is impossible to finish the computing work in an acceptable time. The GIS-based NHBA is an effective algorithm that can easily optimize the complex spatial search. The performance procedure is described as follows.


#### 4.2.1. Define the Objective Function

The objective function is strongly associated with the population density, traffic accessibility, and competitive factors of facilities.(1) Population Density. Population density is an important factor affecting the spatial distribution of facilities. When a facility is located in a position with a large population density, it will undoubtedly increase the passenger volume. The influence of the facilities on the population rapidly attenuates with the increase in distance. Suppose the central position of the facility is P(x, y), the formula of which has an effect on a certain point Pj(xj′, yj′) in the 2-dimensional geographical space is expressed by the following equation:(8)Dj=e−kx−xj′2+y−yj′2,where k is the decay coefficient. The greater the value k, the faster the decay is, and vice versa. When the facility scale is large, its influence is also large, and when the value k is relatively low, it means there is a wider scope of influence. In addition, traffic accessibility has a greater impact on the value k. When the transportation is developed, the value k is relatively lower and the influence scope of the facilities is broader. In this paper, we suppose that every facility scale is fixed, which means the value k is constant. The facility with a good position should try to cover the maximum population density, i.e., the formula: max Dj.(2) Traffic Accessibility. Traffic accessibility is another important factor affecting the spatial distribution of facilities. Convenient traffic can attract a larger passenger volume. As a consequence, the population number will be reduced accordingly with the increasing space distance. Suppose the distance d(x, y) between an arbitrary point Pj(xj′, yj′) and the central position of a certain facility P(x, y) is described by the Euclidean distance and expressed by the following equation:(9)dx,y=x−xj′2+y−yj′2.Traffic accessibility means that the optimal locations of the facilities meet the minimum sum of the distances between the population and the facilities, i.e., the formula: min∑dj(x, y).(3) Competitive Factor. Competition exists between the different facilities so that it is necessary to keep a certain distance between the facilities. If the factor of the facility scale is not considered, the proximity principle will be followed by people. Suppose the minimum distance djmin between the n facilities P1(x1, y1), P2(x2, y2),…, Pn(xn, yn) and the jth population cell Pj(xj′, yj′) exists and is given as follows:(10)djmin=mindjxj,yj=minxi−xj′2+yi−yj′2,where i = 1, 2,…, n; j = 1, 2, …, h. The parameters n and h are the number of facilities and the number of population cells, respectively.


#### 4.2.2. Objective Function

According to the above definitions, the problem should be how to locate the facilities in a certain place when the number of facilities is constant. First, the objective function is defined to solve the problem. The influence of the facilities will gradually attenuate with the distance. Thus, the objective function based on the vector data is calculated as follows:(11)F=C∑i=1n∑j=1hDjAje−kdij∑j=1hDjAjdjmin,where F is the fitness value of the objective function; C is a constant; Dj is the population density of the jth population cell; Aj is the area of the jth population cell, where j = 1,2,…,h; and k is the decay parameter. The parameter dij is equal to xi−xj2+yi−yj2, where dij is the distance between point (xi, yi) and point (xj, yj). The point (xi, yi) represents the coordinates of the ith facility, and the point (xj, yj) represents the coordinates of the jth population cell.The combinatorial optimization problem can be formulated as computing the maximum value. Therefore, the purpose of the research is to obtain the coordinates of facilitiesPg(xgbest, ygbest) under the constraint condition of the maximization of F.


#### 4.2.3. Introduction to NHBA

The NHBA treats each solution as a bat searching inD-dimensional hyperspace. Different from other optimal problems, the optimal spatial search is in a 2-dimensional geographic space in which each point includes X and Y coordinates. The concept of geographic multidimensional space is different from that of D-dimensional hyperspace. D is equal to 2n, where n is the number of facilities. Consequently, the position vector of the ith bat is Bi (xi1, yi1, xi2, yi2, …, xin, yin). Each bat flies over the search space, and its velocity vector is VIvxi1,vyi1,vxi2,vyi2,…,vxin,vyin. The bats can adjust their positions and velocities according to the current optimal value B(j).Best and global optimal value Bgbest. The updating formula of the position is expressed by(12)xit+1=xit+vxit+1,yit+1=yit+vyit+1.The velocity is updated by(13)vxit+1=vxit+xit−xgbesttfxi,vyit+1=vyit+yit−ygbesttfyi,where t is the iteration. The parameters xgbestt and ygbestt are the global optimal positions of X and Y axis directions found by all the bats in the tth iteration. The frequency f is dominated by(14)fxi=fxmin+fxmax−fxminβxi,fyi=fymin+fymax−fyminβyi,whereβxi and βyi are random vectors subjected to a uniform distribution in the range of [0, 1]. Typically, fxmin = fymin = 0, fxmax = fymax = 100, and the initial frequencies of each bat along the X and Y axis directions are uniformly distributed random numbers in the range of [fxmin, fxmax] and [fymin, fymax], respectively:(15)xinewt+1=xioldt+εxit⋅Axit,yinewt+1=yioldt+εyit⋅Ayit,wherein εxit and εyit are random numbers in the range of [−1, 1] and Axit and Ayit are the average loudness of all the bats along the X and Y axis directions in the tth iteration:(16)Axit+1=α⋅Axit,Ayit+1=α⋅Ayit.The parameterα is the adjustment coefficient and takes the appropriate value according to the specific problem:(17)Axit=∑j=1,3,52n−1rxijtn,Ayit=∑j=2,4,62nrxijtn,(18)rxit=rxi0⋅1−e−γt,ryit=ryi0⋅1−e−γt.The parametersrxi0 and ryi0 are the initial pulse rates, rxit and ryit are the tth iteration of the ith bat, and γ is the adjustment coefficient. Because the x and y coordinates pair up, j = 1, 3, 5,…, 2n−1 and j = 2, 4, 6, …, 2n along the X- and Y-axis directions, respectively.


#### 4.2.4. Implementing the Procedure of NHBA for Optimal Spatial Search

Based on the above concepts, the implementation procedure of NHBA is described as follows:(1)
Define the objective functionF(Bi) based on equation (11).(2)
Initialize the parameters for each bat: the population sizeM, position x, velocity v, loudness A, and pulse rate r.(3)
Set the arithmometer that is used to count the number of bats. GenerateM bats from the initial population with random positions and velocities, and calculate the fitness value F(Bi) of each bat of the population based on the objective function in step 1.(4)
Evaluate each bat and find the current optimal bat in each population and save it asB(j).Best.(5)
Begin to iterate from each bat in the population.(5.1)
Calculate the frequency, velocity, and position for theith bat based on the equations (12)–(14).(5.2)
Judge whether the positionsxi and yi of all the bats exceed the boundaries. If so, let xi and yi be equal to the critical values of boundaries.(5.3)
Ifrand > rxi or rand > ryi, generate a new solution for each bat based on equation (15). A random walk procedure has been used to generate a new solution for each bat to improve the variability of the possible solution, where the parameter rand is a random number and rand ∈ [-1, 1].(5.4)
Adjust the loudnessA and pulse rates r of all the bats in the tth iteration. At each iteration of the algorithm, the loudness Axt and Ayt and the pulse rates rxit and ryit are updated based on the equations (16)–(18) if the new solutions are improved. It means that these bats are moving toward the optimal solution.(5.5)
A solution is accepted if a number drawn at random is less than the loudness Ai and meets the relation F(Bi) >Bgbest, where Bgbest is the global optimum value.(5.6)
Update the frequency, velocity, position, loudness, pulse rate, and fitness values of each bat because the values of the parameters change.(6)
Calculate the Euclidean distancedij between two arbitrary bats Bi and Bj in the population, where dij=Bi−Bj=∑k=1nxik−xjk2+yik−yjk2.(7)
Implement the elimination algorithm of the niche. When the relationdij < L is met, compare the fitness values of two bats Bi and Bj. L is the niche distance. It is an empirical value and is given based on experience in the following experiment (see Tables 1 and 2). Impose the stronger penalty factor on the bat with the lower fitness value, i.e., min(F(Bi), F(Bj)) = Penalty, so that the diversity of the population is maintained. The Penalty is a very small number. For example, we let Penalty = 10−6. By doing so, it can help the population improve its average fitness value and accelerate the speed of obtaining the optimal value.(8)
The new population is sorted in descending order by the new fitness values. The bats, of which the fitness values are less than the elimination factorelm, are divided into the population V. The others are divided into the population W and memorized. The bats of population V will be eliminated and reinitialized.(9)
Judge whether it meets the termination condition. If it does not meet the termination condition, update the arithmometert = t + 1. The former W individuals produced in step 8 are treated as the new population of the next iteration. Then, go to step 5. If the termination condition is met, i.e., |Ft+1−Ft| < ε (a critical value) or the maximum iterations are exceeded, end the algorithm.(10)
Obtain the best solution, including the best positions and best fitness values.(11)
Memorize and output the results.Table 1
List of parameters whenn = 1.
Parameters
Meaning
Selected values
M
The number of bats
80
C
Constant
100,000
K
Decay coefficient
0.0004
c
1, c2
Constant
2
E
min
Minimum error
10–6
α
The adjustment coefficient
0.90
γ
Constant,γ > 0
0.96
L
The niche distance
500Table 2
List of parameters whenn > 1.
Parameters
Meaning
Values
M
The number of bats
108
C
Constant
100,000
K
Decay coefficient
0.0004
E
min
Minimum error
10–6
α
The adjustment coefficient
0.86
γ
Constant,γ>0
0.95
L
The niche distance
2000


## 4.1. Brief Description of NHBA

The standard BA has a premature convergence rate and is easily trapped in a local optimum. To overcome these problems, NHBA is proposed by means of the niche technique hybridizing BA. The idea of NHBA is as follows.First, the Euclidean distancedij is calculated between the arbitrary two bats Bati and Batj among the population. If the Euclidean distance dij, denoted as || Bati − Batj ||, is less than the niche distance L which is a threshold value and is given previously based on experience, it signifies that the bats Bi and Bj have higher similarity. The fitness values of two individuals are compared. Then, the lower fitness value is given a penalty factor in order to decrease its fitness value. Regarding the two individuals within the niche distance L, the individuals with lower fitness values will become worse and be eliminated by a large probability in the next generations. It means that there exist only elite individuals within the niche distance L. Therefore, the population diversity is maintained. Meanwhile, the two different individuals maintain a certain distance, and all the individuals are scattered over the whole constrained space. Finally, the niche hybrid bat algorithm is implemented in this way.


## 4.2. GIS-Based NHBA for Optimal Spatial Search

The GIS-based NHBA has the advantage of being applied to optimal spatial search with combined GIS techniques during a sequence of iterations. Optimal spatial search is involved in the resource or location allocation of facilities, which include hospitals, supermarkets, schools, cinemas, and fire stations. Hence, optimal spatial search plays an important role in geographic information science. This paper uses the optimal location allocation of facilities as an example. This type of problem is involved in complex nonlinear optimization scenarios, especially in multiobjective location allocation. If we use the brute-force algorithm, i.e., list all the combinatorial results to obtain the optimal value, the computational volume is so enormous that it is impossible to finish the computing work in an acceptable time. The GIS-based NHBA is an effective algorithm that can easily optimize the complex spatial search. The performance procedure is described as follows.


### 4.2.1. Define the Objective Function

The objective function is strongly associated with the population density, traffic accessibility, and competitive factors of facilities.(1) Population Density. Population density is an important factor affecting the spatial distribution of facilities. When a facility is located in a position with a large population density, it will undoubtedly increase the passenger volume. The influence of the facilities on the population rapidly attenuates with the increase in distance. Suppose the central position of the facility is P(x, y), the formula of which has an effect on a certain point Pj(xj′, yj′) in the 2-dimensional geographical space is expressed by the following equation:(8)Dj=e−kx−xj′2+y−yj′2,where k is the decay coefficient. The greater the value k, the faster the decay is, and vice versa. When the facility scale is large, its influence is also large, and when the value k is relatively low, it means there is a wider scope of influence. In addition, traffic accessibility has a greater impact on the value k. When the transportation is developed, the value k is relatively lower and the influence scope of the facilities is broader. In this paper, we suppose that every facility scale is fixed, which means the value k is constant. The facility with a good position should try to cover the maximum population density, i.e., the formula: max Dj.(2) Traffic Accessibility. Traffic accessibility is another important factor affecting the spatial distribution of facilities. Convenient traffic can attract a larger passenger volume. As a consequence, the population number will be reduced accordingly with the increasing space distance. Suppose the distance d(x, y) between an arbitrary point Pj(xj′, yj′) and the central position of a certain facility P(x, y) is described by the Euclidean distance and expressed by the following equation:(9)dx,y=x−xj′2+y−yj′2.Traffic accessibility means that the optimal locations of the facilities meet the minimum sum of the distances between the population and the facilities, i.e., the formula: min∑dj(x, y).(3) Competitive Factor. Competition exists between the different facilities so that it is necessary to keep a certain distance between the facilities. If the factor of the facility scale is not considered, the proximity principle will be followed by people. Suppose the minimum distance djmin between the n facilities P1(x1, y1), P2(x2, y2),…, Pn(xn, yn) and the jth population cell Pj(xj′, yj′) exists and is given as follows:(10)djmin=mindjxj,yj=minxi−xj′2+yi−yj′2,where i = 1, 2,…, n; j = 1, 2, …, h. The parameters n and h are the number of facilities and the number of population cells, respectively.


### 4.2.2. Objective Function

According to the above definitions, the problem should be how to locate the facilities in a certain place when the number of facilities is constant. First, the objective function is defined to solve the problem. The influence of the facilities will gradually attenuate with the distance. Thus, the objective function based on the vector data is calculated as follows:(11)F=C∑i=1n∑j=1hDjAje−kdij∑j=1hDjAjdjmin,where F is the fitness value of the objective function; C is a constant; Dj is the population density of the jth population cell; Aj is the area of the jth population cell, where j = 1,2,…,h; and k is the decay parameter. The parameter dij is equal to xi−xj2+yi−yj2, where dij is the distance between point (xi, yi) and point (xj, yj). The point (xi, yi) represents the coordinates of the ith facility, and the point (xj, yj) represents the coordinates of the jth population cell.The combinatorial optimization problem can be formulated as computing the maximum value. Therefore, the purpose of the research is to obtain the coordinates of facilitiesPg(xgbest, ygbest) under the constraint condition of the maximization of F.


### 4.2.3. Introduction to NHBA

The NHBA treats each solution as a bat searching inD-dimensional hyperspace. Different from other optimal problems, the optimal spatial search is in a 2-dimensional geographic space in which each point includes X and Y coordinates. The concept of geographic multidimensional space is different from that of D-dimensional hyperspace. D is equal to 2n, where n is the number of facilities. Consequently, the position vector of the ith bat is Bi (xi1, yi1, xi2, yi2, …, xin, yin). Each bat flies over the search space, and its velocity vector is VIvxi1,vyi1,vxi2,vyi2,…,vxin,vyin. The bats can adjust their positions and velocities according to the current optimal value B(j).Best and global optimal value Bgbest. The updating formula of the position is expressed by(12)xit+1=xit+vxit+1,yit+1=yit+vyit+1.The velocity is updated by(13)vxit+1=vxit+xit−xgbesttfxi,vyit+1=vyit+yit−ygbesttfyi,where t is the iteration. The parameters xgbestt and ygbestt are the global optimal positions of X and Y axis directions found by all the bats in the tth iteration. The frequency f is dominated by(14)fxi=fxmin+fxmax−fxminβxi,fyi=fymin+fymax−fyminβyi,whereβxi and βyi are random vectors subjected to a uniform distribution in the range of [0, 1]. Typically, fxmin = fymin = 0, fxmax = fymax = 100, and the initial frequencies of each bat along the X and Y axis directions are uniformly distributed random numbers in the range of [fxmin, fxmax] and [fymin, fymax], respectively:(15)xinewt+1=xioldt+εxit⋅Axit,yinewt+1=yioldt+εyit⋅Ayit,wherein εxit and εyit are random numbers in the range of [−1, 1] and Axit and Ayit are the average loudness of all the bats along the X and Y axis directions in the tth iteration:(16)Axit+1=α⋅Axit,Ayit+1=α⋅Ayit.The parameterα is the adjustment coefficient and takes the appropriate value according to the specific problem:(17)Axit=∑j=1,3,52n−1rxijtn,Ayit=∑j=2,4,62nrxijtn,(18)rxit=rxi0⋅1−e−γt,ryit=ryi0⋅1−e−γt.The parametersrxi0 and ryi0 are the initial pulse rates, rxit and ryit are the tth iteration of the ith bat, and γ is the adjustment coefficient. Because the x and y coordinates pair up, j = 1, 3, 5,…, 2n−1 and j = 2, 4, 6, …, 2n along the X- and Y-axis directions, respectively.


### 4.2.4. Implementing the Procedure of NHBA for Optimal Spatial Search

Based on the above concepts, the implementation procedure of NHBA is described as follows:(1)
Define the objective functionF(Bi) based on equation (11).(2)
Initialize the parameters for each bat: the population sizeM, position x, velocity v, loudness A, and pulse rate r.(3)
Set the arithmometer that is used to count the number of bats. GenerateM bats from the initial population with random positions and velocities, and calculate the fitness value F(Bi) of each bat of the population based on the objective function in step 1.(4)
Evaluate each bat and find the current optimal bat in each population and save it asB(j).Best.(5)
Begin to iterate from each bat in the population.(5.1)
Calculate the frequency, velocity, and position for theith bat based on the equations (12)–(14).(5.2)
Judge whether the positionsxi and yi of all the bats exceed the boundaries. If so, let xi and yi be equal to the critical values of boundaries.(5.3)
Ifrand > rxi or rand > ryi, generate a new solution for each bat based on equation (15). A random walk procedure has been used to generate a new solution for each bat to improve the variability of the possible solution, where the parameter rand is a random number and rand ∈ [-1, 1].(5.4)
Adjust the loudnessA and pulse rates r of all the bats in the tth iteration. At each iteration of the algorithm, the loudness Axt and Ayt and the pulse rates rxit and ryit are updated based on the equations (16)–(18) if the new solutions are improved. It means that these bats are moving toward the optimal solution.(5.5)
A solution is accepted if a number drawn at random is less than the loudness Ai and meets the relation F(Bi) >Bgbest, where Bgbest is the global optimum value.(5.6)
Update the frequency, velocity, position, loudness, pulse rate, and fitness values of each bat because the values of the parameters change.(6)
Calculate the Euclidean distancedij between two arbitrary bats Bi and Bj in the population, where dij=Bi−Bj=∑k=1nxik−xjk2+yik−yjk2.(7)
Implement the elimination algorithm of the niche. When the relationdij < L is met, compare the fitness values of two bats Bi and Bj. L is the niche distance. It is an empirical value and is given based on experience in the following experiment (see Tables 1 and 2). Impose the stronger penalty factor on the bat with the lower fitness value, i.e., min(F(Bi), F(Bj)) = Penalty, so that the diversity of the population is maintained. The Penalty is a very small number. For example, we let Penalty = 10−6. By doing so, it can help the population improve its average fitness value and accelerate the speed of obtaining the optimal value.(8)
The new population is sorted in descending order by the new fitness values. The bats, of which the fitness values are less than the elimination factorelm, are divided into the population V. The others are divided into the population W and memorized. The bats of population V will be eliminated and reinitialized.(9)
Judge whether it meets the termination condition. If it does not meet the termination condition, update the arithmometert = t + 1. The former W individuals produced in step 8 are treated as the new population of the next iteration. Then, go to step 5. If the termination condition is met, i.e., |Ft+1−Ft| < ε (a critical value) or the maximum iterations are exceeded, end the algorithm.(10)
Obtain the best solution, including the best positions and best fitness values.(11)
Memorize and output the results.Table 1
List of parameters whenn = 1.
Parameters
Meaning
Selected values
M
The number of bats
80
C
Constant
100,000
K
Decay coefficient
0.0004
c
1, c2
Constant
2
E
min
Minimum error
10–6
α
The adjustment coefficient
0.90
γ
Constant,γ > 0
0.96
L
The niche distance
500Table 2
List of parameters whenn > 1.
Parameters
Meaning
Values
M
The number of bats
108
C
Constant
100,000
K
Decay coefficient
0.0004
E
min
Minimum error
10–6
α
The adjustment coefficient
0.86
γ
Constant,γ>0
0.95
L
The niche distance
2000


## 4.2.1. Define the Objective Function

The objective function is strongly associated with the population density, traffic accessibility, and competitive factors of facilities.(1) Population Density. Population density is an important factor affecting the spatial distribution of facilities. When a facility is located in a position with a large population density, it will undoubtedly increase the passenger volume. The influence of the facilities on the population rapidly attenuates with the increase in distance. Suppose the central position of the facility is P(x, y), the formula of which has an effect on a certain point Pj(xj′, yj′) in the 2-dimensional geographical space is expressed by the following equation:(8)Dj=e−kx−xj′2+y−yj′2,where k is the decay coefficient. The greater the value k, the faster the decay is, and vice versa. When the facility scale is large, its influence is also large, and when the value k is relatively low, it means there is a wider scope of influence. In addition, traffic accessibility has a greater impact on the value k. When the transportation is developed, the value k is relatively lower and the influence scope of the facilities is broader. In this paper, we suppose that every facility scale is fixed, which means the value k is constant. The facility with a good position should try to cover the maximum population density, i.e., the formula: max Dj.(2) Traffic Accessibility. Traffic accessibility is another important factor affecting the spatial distribution of facilities. Convenient traffic can attract a larger passenger volume. As a consequence, the population number will be reduced accordingly with the increasing space distance. Suppose the distance d(x, y) between an arbitrary point Pj(xj′, yj′) and the central position of a certain facility P(x, y) is described by the Euclidean distance and expressed by the following equation:(9)dx,y=x−xj′2+y−yj′2.Traffic accessibility means that the optimal locations of the facilities meet the minimum sum of the distances between the population and the facilities, i.e., the formula: min∑dj(x, y).(3) Competitive Factor. Competition exists between the different facilities so that it is necessary to keep a certain distance between the facilities. If the factor of the facility scale is not considered, the proximity principle will be followed by people. Suppose the minimum distance djmin between the n facilities P1(x1, y1), P2(x2, y2),…, Pn(xn, yn) and the jth population cell Pj(xj′, yj′) exists and is given as follows:(10)djmin=mindjxj,yj=minxi−xj′2+yi−yj′2,where i = 1, 2,…, n; j = 1, 2, …, h. The parameters n and h are the number of facilities and the number of population cells, respectively.


## 4.2.2. Objective Function

According to the above definitions, the problem should be how to locate the facilities in a certain place when the number of facilities is constant. First, the objective function is defined to solve the problem. The influence of the facilities will gradually attenuate with the distance. Thus, the objective function based on the vector data is calculated as follows:(11)F=C∑i=1n∑j=1hDjAje−kdij∑j=1hDjAjdjmin,where F is the fitness value of the objective function; C is a constant; Dj is the population density of the jth population cell; Aj is the area of the jth population cell, where j = 1,2,…,h; and k is the decay parameter. The parameter dij is equal to xi−xj2+yi−yj2, where dij is the distance between point (xi, yi) and point (xj, yj). The point (xi, yi) represents the coordinates of the ith facility, and the point (xj, yj) represents the coordinates of the jth population cell.The combinatorial optimization problem can be formulated as computing the maximum value. Therefore, the purpose of the research is to obtain the coordinates of facilitiesPg(xgbest, ygbest) under the constraint condition of the maximization of F.


## 4.2.3. Introduction to NHBA

The NHBA treats each solution as a bat searching inD-dimensional hyperspace. Different from other optimal problems, the optimal spatial search is in a 2-dimensional geographic space in which each point includes X and Y coordinates. The concept of geographic multidimensional space is different from that of D-dimensional hyperspace. D is equal to 2n, where n is the number of facilities. Consequently, the position vector of the ith bat is Bi (xi1, yi1, xi2, yi2, …, xin, yin). Each bat flies over the search space, and its velocity vector is VIvxi1,vyi1,vxi2,vyi2,…,vxin,vyin. The bats can adjust their positions and velocities according to the current optimal value B(j).Best and global optimal value Bgbest. The updating formula of the position is expressed by(12)xit+1=xit+vxit+1,yit+1=yit+vyit+1.The velocity is updated by(13)vxit+1=vxit+xit−xgbesttfxi,vyit+1=vyit+yit−ygbesttfyi,where t is the iteration. The parameters xgbestt and ygbestt are the global optimal positions of X and Y axis directions found by all the bats in the tth iteration. The frequency f is dominated by(14)fxi=fxmin+fxmax−fxminβxi,fyi=fymin+fymax−fyminβyi,whereβxi and βyi are random vectors subjected to a uniform distribution in the range of [0, 1]. Typically, fxmin = fymin = 0, fxmax = fymax = 100, and the initial frequencies of each bat along the X and Y axis directions are uniformly distributed random numbers in the range of [fxmin, fxmax] and [fymin, fymax], respectively:(15)xinewt+1=xioldt+εxit⋅Axit,yinewt+1=yioldt+εyit⋅Ayit,wherein εxit and εyit are random numbers in the range of [−1, 1] and Axit and Ayit are the average loudness of all the bats along the X and Y axis directions in the tth iteration:(16)Axit+1=α⋅Axit,Ayit+1=α⋅Ayit.The parameterα is the adjustment coefficient and takes the appropriate value according to the specific problem:(17)Axit=∑j=1,3,52n−1rxijtn,Ayit=∑j=2,4,62nrxijtn,(18)rxit=rxi0⋅1−e−γt,ryit=ryi0⋅1−e−γt.The parametersrxi0 and ryi0 are the initial pulse rates, rxit and ryit are the tth iteration of the ith bat, and γ is the adjustment coefficient. Because the x and y coordinates pair up, j = 1, 3, 5,…, 2n−1 and j = 2, 4, 6, …, 2n along the X- and Y-axis directions, respectively.


## 4.2.4. Implementing the Procedure of NHBA for Optimal Spatial Search

Based on the above concepts, the implementation procedure of NHBA is described as follows:(1)
Define the objective functionF(Bi) based on equation (11).(2)
Initialize the parameters for each bat: the population sizeM, position x, velocity v, loudness A, and pulse rate r.(3)
Set the arithmometer that is used to count the number of bats. GenerateM bats from the initial population with random positions and velocities, and calculate the fitness value F(Bi) of each bat of the population based on the objective function in step 1.(4)
Evaluate each bat and find the current optimal bat in each population and save it asB(j).Best.(5)
Begin to iterate from each bat in the population.(5.1)
Calculate the frequency, velocity, and position for theith bat based on the equations (12)–(14).(5.2)
Judge whether the positionsxi and yi of all the bats exceed the boundaries. If so, let xi and yi be equal to the critical values of boundaries.(5.3)
Ifrand > rxi or rand > ryi, generate a new solution for each bat based on equation (15). A random walk procedure has been used to generate a new solution for each bat to improve the variability of the possible solution, where the parameter rand is a random number and rand ∈ [-1, 1].(5.4)
Adjust the loudnessA and pulse rates r of all the bats in the tth iteration. At each iteration of the algorithm, the loudness Axt and Ayt and the pulse rates rxit and ryit are updated based on the equations (16)–(18) if the new solutions are improved. It means that these bats are moving toward the optimal solution.(5.5)
A solution is accepted if a number drawn at random is less than the loudness Ai and meets the relation F(Bi) >Bgbest, where Bgbest is the global optimum value.(5.6)
Update the frequency, velocity, position, loudness, pulse rate, and fitness values of each bat because the values of the parameters change.(6)
Calculate the Euclidean distancedij between two arbitrary bats Bi and Bj in the population, where dij=Bi−Bj=∑k=1nxik−xjk2+yik−yjk2.(7)
Implement the elimination algorithm of the niche. When the relationdij < L is met, compare the fitness values of two bats Bi and Bj. L is the niche distance. It is an empirical value and is given based on experience in the following experiment (see Tables 1 and 2). Impose the stronger penalty factor on the bat with the lower fitness value, i.e., min(F(Bi), F(Bj)) = Penalty, so that the diversity of the population is maintained. The Penalty is a very small number. For example, we let Penalty = 10−6. By doing so, it can help the population improve its average fitness value and accelerate the speed of obtaining the optimal value.(8)
The new population is sorted in descending order by the new fitness values. The bats, of which the fitness values are less than the elimination factorelm, are divided into the population V. The others are divided into the population W and memorized. The bats of population V will be eliminated and reinitialized.(9)
Judge whether it meets the termination condition. If it does not meet the termination condition, update the arithmometert = t + 1. The former W individuals produced in step 8 are treated as the new population of the next iteration. Then, go to step 5. If the termination condition is met, i.e., |Ft+1−Ft| < ε (a critical value) or the maximum iterations are exceeded, end the algorithm.(10)
Obtain the best solution, including the best positions and best fitness values.(11)
Memorize and output the results.Table 1
List of parameters whenn = 1.
Parameters
Meaning
Selected values
M
The number of bats
80
C
Constant
100,000
K
Decay coefficient
0.0004
c
1, c2
Constant
2
E
min
Minimum error
10–6
α
The adjustment coefficient
0.90
γ
Constant,γ > 0
0.96
L
The niche distance
500Table 2
List of parameters whenn > 1.
Parameters
Meaning
Values
M
The number of bats
108
C
Constant
100,000
K
Decay coefficient
0.0004
E
min
Minimum error
10–6
α
The adjustment coefficient
0.86
γ
Constant,γ>0
0.95
L
The niche distance
2000


## 5. Experiment and Discussion

### 5.1. Data and Methods

The spatial search needs to combine the location data with the attribute data. Regarding the problems of locating the facilities in a certain place, GIS can provide robust support. Some important information, which includes the centroid coordinates (xi, yi) and area of each population cell and the width and height of the search space, needs to be extracted from the map of the population density distribution by means of GIS. Python, a common language tool, is used in this research. We can import and manipulate the input data and thus present the generated map.The digitalized population density map is a polygon vector map of Fangcun District, Guangzhou City, China. Each polygon represents a population cell, which is the minimal unit of the geographic area. The fields should includeOID and Population in the attribute database of GIS. The field OID represents the object identity code, which is unique. The field Population represents the population density of each population cell, which can be extracted from the attribute database of GIS. Compared with the raster data, the precision of the vector data we use in the paper is high.The map of the population density in Fangcun District, Guangzhou City, China, is shown in Figure1.Figure 1
Map of population density in Fangcun District.


### 5.2. Case for Single Facility (n = 1)

This is the simplest case whenn = 1 because there is no competition from other facilities. The parameters are listed in Table 1.The optimal value isFmax = 22.0968 after 6 iterations based on NHBA. The results are shown in Figure 2.Figure 2
Map of the population density and facility distribution (n = 1) based on NHBA in Fangcun District.We compare NHBA with BA, GA, FA, PSO, and the brute-force algorithm under the same parameters (see Table1) by our programming in C language. The comparison results are shown in Tables 3 and 4.Table 3
The comparison results of NHBA with BA, GA, FA, and PSO.
Algorithms
Optimal values
Iterations
Time complexities
BA
22.0967
7
1.22 × 106
GA
22.0965
7
1.22 × 106
FA
21.0478
95
1.66 × 107
PSO
22.0968
24
4.20 × 106
NHBA
22.0968
6
1.05 × 106Table 4
Results of the brute-force algorithm based on the spatial sampling interval.
Spatial sampling interval
Fitness values
X
Y
Time complexities
100×100
21.188921
34579.5
24920.4
8712
10×10
22.095361
34079.5
25020.4
863625
1×1
22.096775
34066.5
25038.4
8.6 × 107The optimal value of BA isFmax = 22.0967 after 7 iterations. NHBA has fewer iterations and a higher precision than BA under the same parameters. GA achieves a similar optimal value with the same iterations compared with BA. The optimal value of PSO is Fmax = 22.0968 after 24 iterations. The results are similar among PSO, BA, and NHBA. However, PSO needs more iterations. FA obtains a smaller optimal value with more iterations. It means that FA is the worst algorithm and NHBA is the best algorithm among the five algorithms.The brute-force algorithm is used in order to validate the results. The time complexities of the algorithm and the accuracies of the solutions are closely related to the spatial sampling interval. The results of the fitness values, time complexities, andx and y coordinates are shown in Table 4.The smaller the spatial sampling interval is, the greater the time complexities are, and the higher the accuracies of the solutions are. When the spatial sampling interval is 1 × 1 meters, the fitness value is 22.096775, which is roughly similar to most algorithms. However, the time complexity of the brute-force algorithm isO (8.6 × 107). BA needs 7 iterations, and the time complexity is O (1.22 × 106). NHBA only needs 6 iterations, and the time complexity is O (1.05 × 106). Consequently, NHBA achieves a faster convergence speed.The 3D perspective map regarding the fitness values and geographical positions in Fangcun District is shown in Figure3.Figure 3
The 3D perspective map of the fitness values and geographical positions in Fangcun District.The map is obtained by means of developing techniques based on Python and R languages. First, evenly generate 100∗100 original sample points in the 2D space within the Fangcun District boundary. Then, calculate the fitness values as zi coordinates corresponding to every point (xi, yi) within the boundary. Finally, create the 3D perspective map based on every point (xi, yi, zi) in R language.In Figure3, the x and y coordinates represent the position of a facility; the z coordinate represents its fitness value. We can roughly judge the optimal position in which a facility is located based on Figure 3.


### 5.3. Case for Two and More Facilities (n > 1)

When the number of facilitiesn is more than 1, it is considered a combinatorial optimization problem, which is computationally difficult. For example, if the brute-force algorithm is used, the time complexity is O (7.65 × 1058) when the spatial sampling interval is 1 meter and n = 8. Because the computing time is unacceptable, the brute-force algorithm has to be abandoned.The parameters are listed in Table2. The optimal value Fmax = 436.94 after 13 iterations (GENNHBA = 13) based on NHBA. The results are shown in Figure 4.Figure 4
Map of the population density and facilities distribution (n = 8) based on NHBA in Fangcun District.In general, the parameterm, namely, the population size or number of bats, is given based on the researcher’s experience. In this paper, the parameter m is obtained by conducting experiments.In Figure5, the larger the parameter m is, the larger the fitness value is (that means that the results are better) when m is lower than 108. However, the fitness value becomes lower when m is equal to 110. In addition, the computation time will increase with the increasing m value. As a result, the parameter m is 108 when the number of facilities is 8. The adjustment coefficients are α = 0.86 and γ = 0.95 in this experiment. The comparative results among the BA, GA, FA, PSO, and NHBA algorithms based on the same parameters are as follows.Figure 5
The relationship between the parameterm and the fitness values.The fitness values are calculated by the objective functionF(Bi) based on equation (1). The larger the fitness values are, the better the results are. In Figure 6, the fitness values of BA, GA, FA, PSO, and NHBA increase with the iterations. The fitness values of NHBA are larger than those of the other algorithms under the same iterations. The fitness values of NHBA tend to be stable under 13 iterations (GENNHBA = 13), while BA needs 15 iterations, GA needs 17 iterations, FA needs 68 iterations, and PSO needs 54 iterations. The maximum fitness values of BA, GA, FA, PSO, and NHBA are 413.42, 324.80, 363.32, 430.85, and 436.94, respectively. It is easier for NHBA to obtain the optimal values than it is for the other algorithms. This finding means that the NHBA approach has the ability to converge to optimal solutions more quickly than other algorithms.Figure 6
The comparative results among the BA, GA, FA, PSO, and NHBA algorithms with the fitness values and iterations.The time complexities are defined asO (GEN∗m∗2n∗h). Here, GEN is the number of iterations, m is the population size (m = 108), n is the number of facilities (n = 8), and h is the number of population cells (h = 1093). Consequently, the time complexities of BA, GA, FA, PSO, and NHBA are approximately O (2.83 × 107), O (3.21 × 107), O (1.28 × 108), O (1.02 × 108), and O (2.46 × 107) in this experiment, respectively (Table 5). As mentioned above, the time complexity of the brute-force algorithm is O (7.65 × 1058) under the same condition. In conclusion, the convergence of NHBA is faster than that of the other algorithms. Furthermore, the computation time of BA, GA, FA, PSO, or NHBA is significantly less than that of the brute-force algorithm.Table 5
Comparison results of BA, GA, FA, PSO, and NHBA.
Methods
Fitness values
Iterations
Time complexities
BA
413.42
15
2.83 × 107
GA
324.80
17
3.21 × 107
FA
363.32
68
1.28 × 108
PSO
430.85
54
1.02 × 108
NHBA
436.94
13
2.46 × 107


### 5.4. Experiment Based on the Expanded Shaffer Function

The NHBA has random elements, and the statistics of only one experimental result are not convincing. The expanded Shaffer function (equations (19) and (20)) is used to test the above algorithms:(19)gx,y=0.5+sin2x2+y2−0.51.0+0.001x2+y22st−100<x,y<100,(20)fx=gx1,x2+gx2,x3+gx3,x4+…+gxD−1,xD+gxD,x1.Whenx(x1, x2,…, xD) locates in the position (0, 0, …, 0), the minimum value of the expanded Shaffer function is equal to 0. After every algorithm is run 50 times in D-dimensional space (D = 50), the statistical means, the best optimal values, the worst optimal values, and the variances of experimental results are obtained and listed in Table 6.Table 6
Comparison results of NHBA and other algorithms based on the expanded Shaffer function after 50 runs (D = 50).
Methods
Means
Best optimal values
Worst optimal values
Variances
BA
21.47
19.58
22.79
460.91
GA
23.78
23.22
24.38
565.49
FA
15.28
13.25
16.98
233.50
PSO
17.59
14.60
20.06
1.30
NHBA
13.14
9.20
18.03
172.78In Table6, from the means, we can see that the mean of NHBA is 13.14, which is the smallest. The second smallest mean is achieved by FA (15.28), and the others in ascending order are as follows: PSO (17.59), BA (21.47), and GA (23.78). The order of the best optimal values is the same as that of the means. However, the smallest value of the worst optimal values is obtained by FA instead of NHBA. The smallest variance is obtained by PSO (1.30). It means that the numerical change range of PSO is the smallest. Because the means can better reflect the overall level, it means that NHBA surpasses other methods and GA is worse than other algorithms.


## 5.1. Data and Methods

The spatial search needs to combine the location data with the attribute data. Regarding the problems of locating the facilities in a certain place, GIS can provide robust support. Some important information, which includes the centroid coordinates (xi, yi) and area of each population cell and the width and height of the search space, needs to be extracted from the map of the population density distribution by means of GIS. Python, a common language tool, is used in this research. We can import and manipulate the input data and thus present the generated map.The digitalized population density map is a polygon vector map of Fangcun District, Guangzhou City, China. Each polygon represents a population cell, which is the minimal unit of the geographic area. The fields should includeOID and Population in the attribute database of GIS. The field OID represents the object identity code, which is unique. The field Population represents the population density of each population cell, which can be extracted from the attribute database of GIS. Compared with the raster data, the precision of the vector data we use in the paper is high.The map of the population density in Fangcun District, Guangzhou City, China, is shown in Figure1.Figure 1
Map of population density in Fangcun District.


## 5.2. Case for Single Facility (n = 1)

This is the simplest case whenn = 1 because there is no competition from other facilities. The parameters are listed in Table 1.The optimal value isFmax = 22.0968 after 6 iterations based on NHBA. The results are shown in Figure 2.Figure 2
Map of the population density and facility distribution (n = 1) based on NHBA in Fangcun District.We compare NHBA with BA, GA, FA, PSO, and the brute-force algorithm under the same parameters (see Table1) by our programming in C language. The comparison results are shown in Tables 3 and 4.Table 3
The comparison results of NHBA with BA, GA, FA, and PSO.
Algorithms
Optimal values
Iterations
Time complexities
BA
22.0967
7
1.22 × 106
GA
22.0965
7
1.22 × 106
FA
21.0478
95
1.66 × 107
PSO
22.0968
24
4.20 × 106
NHBA
22.0968
6
1.05 × 106Table 4
Results of the brute-force algorithm based on the spatial sampling interval.
Spatial sampling interval
Fitness values
X
Y
Time complexities
100×100
21.188921
34579.5
24920.4
8712
10×10
22.095361
34079.5
25020.4
863625
1×1
22.096775
34066.5
25038.4
8.6 × 107The optimal value of BA isFmax = 22.0967 after 7 iterations. NHBA has fewer iterations and a higher precision than BA under the same parameters. GA achieves a similar optimal value with the same iterations compared with BA. The optimal value of PSO is Fmax = 22.0968 after 24 iterations. The results are similar among PSO, BA, and NHBA. However, PSO needs more iterations. FA obtains a smaller optimal value with more iterations. It means that FA is the worst algorithm and NHBA is the best algorithm among the five algorithms.The brute-force algorithm is used in order to validate the results. The time complexities of the algorithm and the accuracies of the solutions are closely related to the spatial sampling interval. The results of the fitness values, time complexities, andx and y coordinates are shown in Table 4.The smaller the spatial sampling interval is, the greater the time complexities are, and the higher the accuracies of the solutions are. When the spatial sampling interval is 1 × 1 meters, the fitness value is 22.096775, which is roughly similar to most algorithms. However, the time complexity of the brute-force algorithm isO (8.6 × 107). BA needs 7 iterations, and the time complexity is O (1.22 × 106). NHBA only needs 6 iterations, and the time complexity is O (1.05 × 106). Consequently, NHBA achieves a faster convergence speed.The 3D perspective map regarding the fitness values and geographical positions in Fangcun District is shown in Figure3.Figure 3
The 3D perspective map of the fitness values and geographical positions in Fangcun District.The map is obtained by means of developing techniques based on Python and R languages. First, evenly generate 100∗100 original sample points in the 2D space within the Fangcun District boundary. Then, calculate the fitness values as zi coordinates corresponding to every point (xi, yi) within the boundary. Finally, create the 3D perspective map based on every point (xi, yi, zi) in R language.In Figure3, the x and y coordinates represent the position of a facility; the z coordinate represents its fitness value. We can roughly judge the optimal position in which a facility is located based on Figure 3.


## 5.3. Case for Two and More Facilities (n > 1)

When the number of facilitiesn is more than 1, it is considered a combinatorial optimization problem, which is computationally difficult. For example, if the brute-force algorithm is used, the time complexity is O (7.65 × 1058) when the spatial sampling interval is 1 meter and n = 8. Because the computing time is unacceptable, the brute-force algorithm has to be abandoned.The parameters are listed in Table2. The optimal value Fmax = 436.94 after 13 iterations (GENNHBA = 13) based on NHBA. The results are shown in Figure 4.Figure 4
Map of the population density and facilities distribution (n = 8) based on NHBA in Fangcun District.In general, the parameterm, namely, the population size or number of bats, is given based on the researcher’s experience. In this paper, the parameter m is obtained by conducting experiments.In Figure5, the larger the parameter m is, the larger the fitness value is (that means that the results are better) when m is lower than 108. However, the fitness value becomes lower when m is equal to 110. In addition, the computation time will increase with the increasing m value. As a result, the parameter m is 108 when the number of facilities is 8. The adjustment coefficients are α = 0.86 and γ = 0.95 in this experiment. The comparative results among the BA, GA, FA, PSO, and NHBA algorithms based on the same parameters are as follows.Figure 5
The relationship between the parameterm and the fitness values.The fitness values are calculated by the objective functionF(Bi) based on equation (1). The larger the fitness values are, the better the results are. In Figure 6, the fitness values of BA, GA, FA, PSO, and NHBA increase with the iterations. The fitness values of NHBA are larger than those of the other algorithms under the same iterations. The fitness values of NHBA tend to be stable under 13 iterations (GENNHBA = 13), while BA needs 15 iterations, GA needs 17 iterations, FA needs 68 iterations, and PSO needs 54 iterations. The maximum fitness values of BA, GA, FA, PSO, and NHBA are 413.42, 324.80, 363.32, 430.85, and 436.94, respectively. It is easier for NHBA to obtain the optimal values than it is for the other algorithms. This finding means that the NHBA approach has the ability to converge to optimal solutions more quickly than other algorithms.Figure 6
The comparative results among the BA, GA, FA, PSO, and NHBA algorithms with the fitness values and iterations.The time complexities are defined asO (GEN∗m∗2n∗h). Here, GEN is the number of iterations, m is the population size (m = 108), n is the number of facilities (n = 8), and h is the number of population cells (h = 1093). Consequently, the time complexities of BA, GA, FA, PSO, and NHBA are approximately O (2.83 × 107), O (3.21 × 107), O (1.28 × 108), O (1.02 × 108), and O (2.46 × 107) in this experiment, respectively (Table 5). As mentioned above, the time complexity of the brute-force algorithm is O (7.65 × 1058) under the same condition. In conclusion, the convergence of NHBA is faster than that of the other algorithms. Furthermore, the computation time of BA, GA, FA, PSO, or NHBA is significantly less than that of the brute-force algorithm.Table 5
Comparison results of BA, GA, FA, PSO, and NHBA.
Methods
Fitness values
Iterations
Time complexities
BA
413.42
15
2.83 × 107
GA
324.80
17
3.21 × 107
FA
363.32
68
1.28 × 108
PSO
430.85
54
1.02 × 108
NHBA
436.94
13
2.46 × 107


## 5.4. Experiment Based on the Expanded Shaffer Function

The NHBA has random elements, and the statistics of only one experimental result are not convincing. The expanded Shaffer function (equations (19) and (20)) is used to test the above algorithms:(19)gx,y=0.5+sin2x2+y2−0.51.0+0.001x2+y22st−100<x,y<100,(20)fx=gx1,x2+gx2,x3+gx3,x4+…+gxD−1,xD+gxD,x1.Whenx(x1, x2,…, xD) locates in the position (0, 0, …, 0), the minimum value of the expanded Shaffer function is equal to 0. After every algorithm is run 50 times in D-dimensional space (D = 50), the statistical means, the best optimal values, the worst optimal values, and the variances of experimental results are obtained and listed in Table 6.Table 6
Comparison results of NHBA and other algorithms based on the expanded Shaffer function after 50 runs (D = 50).
Methods
Means
Best optimal values
Worst optimal values
Variances
BA
21.47
19.58
22.79
460.91
GA
23.78
23.22
24.38
565.49
FA
15.28
13.25
16.98
233.50
PSO
17.59
14.60
20.06
1.30
NHBA
13.14
9.20
18.03
172.78In Table6, from the means, we can see that the mean of NHBA is 13.14, which is the smallest. The second smallest mean is achieved by FA (15.28), and the others in ascending order are as follows: PSO (17.59), BA (21.47), and GA (23.78). The order of the best optimal values is the same as that of the means. However, the smallest value of the worst optimal values is obtained by FA instead of NHBA. The smallest variance is obtained by PSO (1.30). It means that the numerical change range of PSO is the smallest. Because the means can better reflect the overall level, it means that NHBA surpasses other methods and GA is worse than other algorithms.


## 6. Conclusions

Complicated combinatorial optimization problems are involved in the optimal spatial search and are difficult to be solved by traditional algorithms such as the brute-force algorithm. In this paper, a GIS-based NHBA is proposed for solving the optimal spatial search. It is characterized as follows:(1)
The GIS-based NHBA is able to maintain the population diversity, enhance the exploration of the new search domain, and overcome the deficiency that occurs in the premature convergence to only one solution of optimal spatial search.(2)
The spatial search occurs in multidimensional geographical space. Optimal spatial search involves more complicated optimization problems. The GIS-based NHBA provides robust support for this work.(3)
The parameter population sizem is traditionally obtained by the researcher’s experience. It is derived by the experiments conducted in this paper in order to improve the precision of the results.(4)
The BA, GA, FA, PSO, and NHBA outperform the brute-force algorithm in terms of computational efficiency.(5)
The optimal value is easier to obtain by employing NHBA than it is by using BA, GA, FA, and PSO. Furthermore, the precision of NHBA is higher and the convergence of NHBA is faster than that of BA, GA, FA, and PSO under the same conditions.



---

*Source: 2843436-2020-05-11.xml*