# Stability Analysis of the Bat Algorithm Described as a Stochastic Discrete-Time State-Space System

Authors: Janusz Piotr Paplinski
Journal: Complexity
 (2018)

---

## Abstract

The main problem with the soft-computing algorithms is a determination of their parameters. The tuning rules are very general and need experiments during a trial and error method. The equations describing the bat algorithm have the form of difference equations, and the algorithm can be treated as a stochastic discrete-time system. The behaviour of this system depends on its dynamic and preservation stability conditions. The paper presents the stability analysis of the bat algorithm described as a stochastic discrete-time state-space system. The observability and controllability analyses were made in order to verify the correctness of the model describing the dynamic of BA. Sufficient conditions for stability are derived based on the Lyapunov stability theory. They indicate the recommended areas of the location of the parameters. The analysis of the position of eigenvalues of the state matrix shows how the different values of parameters affect the behaviour of the algorithm. They indicate the recommended area of the location of the parameters. Simulation results confirm the theory-based analysis.


---

## Body

## 1. Introduction

In recent years, the nature-inspired metaheuristic algorithms for optimization problems become very popular. Though these algorithms do not guarantee the optimal solution, they generally have a tendency to find a good solution and become powerful methods for solving many difficult optimization problems [1–3]. The heuristic methods are based on the many different mechanisms occurring in nature. The genetic algorithms [4] are based on the biological fundamentals, tabu search is based on the social behaviour [5], and ant colony optimization [6, 7] or particle swarm optimization (PSO) [8] is based on the swarm behaviour. The bat algorithm (BA) proposed by Yang [9] belongs to the last. The bats use some type of sonar, called echolocation, to detect prey or to avoid obstacles. The echolocation guides their search and allows discrimination of different types of insects, even in complete darkness.There are some much powerful modifications of BA, for instance, BA based on differential operator and Lévy flight trajectory (DLBA) proposed by Xie et al. [10], the improved bat algorithm (IBA) proposed by Yilmaz and Küçüksille [11], and enhanced bat algorithm proposed also by Yilmaz and Küçüksille [12], but the simple BA is a base and is more popular than other modifications. For this reason, the analysis for simple BA is made in the paper.The main problem with the soft-computing algorithms is a determination of their parameters. The tuning rules are very general and need experiments during a trial and error method. The main idea in this process is to balance the running algorithm between exploration and exploitation. In the case of too little exploration and intensive exploitation, the algorithm can converge to a local optimum. Otherwise, too much exploration and too little exploitation can give the algorithm with a very small convergence [13–15].The behaviour of the algorithm and ability to converge to the global optimum depend on its dynamic, which is described by difference equations. This behaviour depends on the stability works of algorithms especially. There are some stability analyses of PSO algorithm based on the location of the roots [16–18] and the Lyapunov function [19–21]. We made a similar study for BA in [22].The present paper is the extended version of the conference paper [22] published in INnovations in Intelligent SysTems and Applications (INISTA) 2017. The way of the approach presented in the paper is largely new and much more general. New in the paper is an extension of the dynamic description of the bat algorithm to the third order by taking into account all variables used by individuals in the population. Thereby, the description of BA becomes more complete and general. Then the use of a description in the form of the state-space and checking the controllability and observability allow reducing the order of dynamic. The method of stability analysis proposed in [22] is based on the location of the roots of the difference equation which is appropriate for linear time-invariant systems. This method was used after omitting the randomness of parameters and treatment of the algorithm as stationary. In the present paper, the new approach of the stability analysis of the dynamics of BA by using the Lyapunov stability theorem and Sylvester’s criterion is made. It provides to obtain the desired range of parameters providing stability work of the algorithm. The Lyapunov theory defines the necessary and sufficient conditions for absolute stability of nonlinear and/or time-varying systems. Any simplifications are not needed during stability analysis, which is important because the obtained results are more reliable than the results from the previous work [22]. Both methods give the similar solution, then the present paper confirms the correctness of simplicity used in preceding work [22]. As an illustrative example, the same four benchmark functions were used, but the graphs were made for not presented earlier Griewank function.The paper is organized as follows. In Section2, the BA algorithm is described. The next section details the Lyapunov stability theory. After this, the dynamic of BA is described and analysed with a special focus on the stability condition. The example experiments and discussion are presented in the last section.


## 2. Bat Algorithm

The bats have fascinating abilities such as finding prey and discriminating different types of insects even in complete darkness. Bats use echolocation by emitting high-frequency audio signals and receiving a reflection of those. The time delay between emission and detection of the echo and its variation of loudness allow bats to recognize surroundings.The metaheuristic BA uses some simplicity and idealized rules:(1)
All bats use echolocation to appoint a distance and direction to the food. They also can recognize the difference between food/prey and background barriers
(2)
Thei-th bat is at position xi and flies randomly with a velocity vi. It emits an audio signal with a variable frequency between fmin,fmax, a varying wavelength λi, and loudness Ai to search for food. It can automatically adjust the wavelength (or frequency) of its emitted pulses and adjust the rate of pulse emission r∈0,1, depending on the proximity of its target
(3)
The loudness can vary in many ways. We assume that the loudness varies from a large (positive)L0 to the smallest constant value LminEach of the artificial bats in thek-th step has a position vector xik, velocity vector vik, and frequency vector fik which is updated during iterations by using the below relations, from (1) to (3). The position vector of the bat represents some specific solution of the optimization problem. Every bat emits an audio signal with a randomly assigned frequency fik, which is drawn uniformly from the range fmin,fmax:
(1)fik=fmin+fmax–fminβk,where βk∈0,1 is a random vector with a uniform distribution. The velocity of the i-th bat in the k-th step vik depends on the position of the current global best solution achieved so far xbsfk:
(2)vik=vik−1+xik−1−xbsfkfik.The new position of the bat and thus the new solution of the problem follow from his earlier position and velocity:(3)xik=xik−1+vik.The local search procedure is also used. A new solution for a bat is generated locally using current best solution and local random walk:(4)xnew=xold+εLk,where ε∈−1,1 is a random number with a uniform distribution, while Lk is the average loudness of all bats at the k-th time step.We can consider BA as a balanced combination of exploration, realized by an algorithm similar to the standard particle swarm optimization and exploitation realized by an intensive local search. The balance between these techniques is controlled by the loudnessL and emission rate r, updated as follows:
(5)Lik+1=αLik,(6)rik+1=ri01–exp–γk,where the coefficients α and γ are constants. In the simplification case, α=γ is often used. We can consider the parameter α as similar to the cooling factor in simulated annealing. The loudness and the pulse emission rate are updated only when the new solution is improved. It means when the bat is moving toward the best solution.The operation of BA can be described as follows. In the beginning, the metaheuristic BA initializes a population of bats, by assigning to individuals values of its parameters. They are most commonly defined randomly. Every bat will move from initial solutions toward the global best solution with each iteration using the position of the current global best solution attained so far. If any bat finds a better solution after moving, the best so far solution, pulse emission and loudness are updated. This process is repeated continuously until the termination criteria are satisfied. The best so far solution achieved is considered the final best solution. The pseudocode of BA is presented in Pseudocode1.Pseudocode 1: The pseudocode of BA.
1. Initialize the bat population xi0i=1,2…,nandvi0
2. Initialize pulse frequency fi0, pulse rates ri0, and the loudness Li0
3. While (the stop condition is not fulfilled)
4.  Generate new solutions by adjusting frequency,
5.  Updating velocities and locations (eq. (2,3))
6.  If rand>rik
7.   Select a solution among the best solutions
8.   Generate a local solution (eq. (4))
9.  End if
10.  Generate a new solution by flying randomly
11.  If rand<Lik&fxik<fxbsfk
12.   Accept the new solutions
13.   Reduce Lik+1 and Increase rik+1 (eq. (5,6))
14.  End if
15.  Rank bats and find current best xbsfk
16. End while


## 3. Lyapunov Stability Theory

The basic theorems for system stability are the Lyapunov stability theorems [23–25]. The second Lyapunov criterion (direct method) allows proving the local and global stabilities of an equilibrium point using proper scalar functions, called Lyapunov functions, defined in the state-space. This criterion refers to particular “positive definite” or “positive semidefinite” scalar functions, which often have the meaning of “energy functions.” By looking at how this energy-like function changes over time, we might conclude that a system is stable or asymptotically stable without solving the nonlinear differential equation. The necessary and sufficient condition for the Lyapunov function, described by the matrix, to be positive definite is determined by Sylvester’s criterion.Theorem 1 (see [26]).
Consider the equilibrium pointx=0 of the stochastic discrete-time system, defined by the state-space equation:
(7)xk+1=Akxk,where xk is a state vector at time k and Ak∈Rn×n is a nonsingular matrix with stochastic values. The equilibrium point is asymptotically stable if there is a nonnegative scalar Lyapunov function Vxk defined as
(8)Vxk=xkTPxk,where P is a positive definite symmetric matrix, with V0=0, which satisfies that the expected value of changes of the Lyapunov function E∆V is greater than zero:
(9)EΔV=EVxk+1−Vxk<0.
We can write it as(10)EΔV=Exk+1TPxk+1−xkTPxk<0,which after simple transformation, using (7) and (10), leads to a useful formula:
(11)EΔV=ExkTATPA−Pxk<0.Remark 1.
The stochastic discrete-time system (7) is asymptotically stable if and only if for any positive definite matrix Q there exists a positive definite symmetric matrix P that satisfies the Lyapunov equation [27, 28]:
(12)EATPA−P=−Q.
The symmetric matrixP is positive definite if it fulfils Sylvester’s criterion.Theorem 2 (see [29]).
The necessary and sufficient condition for the matrix to be positive definite is that the determinants of all the successive principal minors of the matrix are positive.
For the symmetric matrixP defined as
(13)P=p11p12⋯p1np21p22⋯⋯⋯⋯⋯⋯pn1pn2⋯pnn,where P=PT; the successive principal minors must be positive:
(14)p11>0,p11p12p21p22>0,…,p11p12⋯p1np21p22⋯⋯⋯⋯⋯⋯pn1pn2⋯pnn>0.


## 4. Analysis of the Bat Algorithm as a Dynamic System

By assuming the velocity of the batvik, its position xik and loudness Lik, as state variables ξk=vikxikLikT and the position of the current global best solution achieved so far xbsfk as an input uk, the dynamics of the BA described by (2), (3), and (5) can be presented in the state-space form:
(15)vik+1xik+1Lik+1=1fik011+fikε00avikxikLik+−fik−fik0xbsfk,yik=010vikxikLik,where the locus of the i-th bat is treated as an output yik.Equation (6) describing the emission rate affects only on the control of the algorithm in step 6 of its pseudocode (Pseudocode 1). For some random iteration, the local search around the best individual is made. It has no influence on the trajectory of individuals and was omitted in the description of the dynamics of BA in relations (15).Each bat, thei-th dimension of (15), updates independently from the others; thus, without losing the generality, the analysis of the algorithm can be reduced to the one-dimensional case. Therefore, consequently to the general form of the state-space:
(16)ξk+1=Akξk+Buk,yk=Cξk+D,and a one-dimensional case of (15), we obtain the state matrix A, the input matrix B, the output matrix C, and the feedforward matrix D:
(17)Ak=1fk011+fkε00a,B=−fk−fk0,C=010,D=0.For the dynamic system described by the state-space form, the analyses of its observability and controllability are essential. Without losing the generality of this analysis, we can assume the constant value of frequencyfk equals its expected value fk=Efk=f.Theorem 3 (see [30]).
The system is completely observable if any initial state vectorxt0 can be reconstructed by examining the output of the system yt over the finite period of time from t0 to tf. Then the system is completely observable if and only if the set of vectors CCA…CAn−1T is literary independent, i.e., rankCCA…CAn−1T=n, where n is a size of state-space vector ξk∈Rn×1.
The rank of the matrix of observability for BA dynamics (15) is equal:
(18)rankCCACA2T=rank01011+fε1+f2+f1+2fε1+f+a=3,which means the system is observable.Theorem 4 (see [23]).
The system is completely controllable if there exists a control signalut defined over a finite interval 0≤t≤tf, which can force system states xtf to any desired value. Then the system is completely controllable if and only if the set of vectors BAB…An−1B is literary independent, i.e., rankBAB…An−1B=n.
The rank of the matrix of controllability for BA dynamics described by (15) is equal:
(19)rankBABA2B=−f−f2−f−f3−3f2−f−f−f2−2f−f3−4f2−3f000=2.
The obtained rank of matrixA equals 2 and is less than the size of the state-space vector n=3. The system (15) is uncontrollable. It results that the order of the system is reduced from third to second. We can see it as simplifying the common expression z−a from the numerator and denominator while calculating the transfer function describing this system:
(20)Gz=CzI−A−1B+D=−fzz−az2−2+fz+1z−a=−fzz2−2+fz+1.
The uncontrollable state-variable is the loudnessLk, but if we assume α∈0,1, then Lk is decreasing during iteration and can be treated by the system as some disturbance. After leave out the earlier simplification and take into account the variability of the value of frequency fk and present the state-space form as
(21)vik+1xik+1=1fk11+fkvikxik+−fk−fkxbsf+0εLik,yik=01vikxik.
The description of the system can be modified to the autonomous form by taking new state variables and neglecting disturbances with the equilibrium point atξk=00T:
(22)x~k=xk−xbsf,then state-space take the form
(23)vik+1x~ik+1=1fk11+fkvikx~ik,yik=01vikx~ik.
According to Theorem1, for any positive definite matrix Q, there must exist symmetric positive definite matrix P fulfilling the Lyapunov function (8). We can define matrix Q as
(24)Q=c100c2,with c1,c2>0. The symmetric matrix P equals
(25)P=p1p2p2p3.
Substituting the matricesP and Q into relation (12), we obtain
(26)Ep1+2p2+p3fkp1+1+2fkp2+1+fkp3fkp1+1+2fkp2+1+fkp3fk2p1+2fk1+fkp2+1+fk2p3=−c100−c2,which leads to the system of relations:
(27)p1+2p2+p3=−c1,Efkp1+1+2Efkp2+1+Efkp3=0,Efk2p1+2Efk1+Efkp2+1+Efk2p3=−c2.
Assumingf=Efk, we can easily calculate that
(28)p2=f2p1−c1,p3=−fp1+−1+fc1.
MatrixP must be positive definite. According to Sylvester’s criterion in Theorem 2, the successive principal minors must be positive, which leads to equations
(29)p1>0,p1p3−p22>0,where
(30)p1p3−p22=−f24+fp12+f22+f−1c1p1−f24c12>0.
Considering the quadratic polynomialap12+bp1+c, we obtain
(31)a=−f24+f,b=f22+f−1c1,c=−f24c12<0,and well-known relations
(32)Δ=b2−4ac,λ1,2=−b±Δ2a.
We will look for a range of parameterf for which the polynomial (30) is greater than zero. We can divide the range of value of frequency f into some subsets.
(I)
a>0
This is met for the mean value of frequencyf∈−4,0. Assuming the big enough value of p11>λ1,2 and p11>0, we obtain that BA is stable in this range of f.
(II)
a<0
The conditionp11>0 can be satisfied only when Δ>0, which is met for f<0.5. We need to examine two conditions:
(i)
f<−4
Because ofa<0 and c<0, b2>Δ, and consequently, λ1,2<0. The condition λ1>p11>λ2>0 indicate unstable dynamics of BA.
(ii)
0<f<0.5
For this range off, the value of b<0 and consequently λ1,2<0 which indicate unstable dynamics of BA.
Lyapunov theory leads to the conclusion that the dynamic of BA is stable for the frequency belonging to the rangef∈−4,0.Theorem 5.
A linear discrete-time system described by the state (16) is asymptotically stable if and only if all eigenvalues of A have magnitude smaller than one. For eigenvalues lying on the unit circle, the system is on the stability border.
The eigenvalues ofA are calculated from the characteristic equation defined as
(33)detzI−A=0,where I=1001, or they are defined as the roots of the denominator of the transfer function (20) equivalently, which leads to the equation
(34)z2−2+fz+1=0.
The roots are equal:(35)z1,2=2+f±f2+4f2.
Figure1 presents the eigenvalues of matrix A on the z-plane. For f∈−4,0, the eigenvalues are complex, and their absolute value equals one ∣z1,2∣=1. In this case, they lie on the stability border, and the algorithm behaves as an undamped oscillator. The samples of the responses of the algorithm are presented in Figures 2(b) and 2(c). The response of the algorithm is periodic and oscillates around the best individuals. For f<−2, the ringing occurs, and the unit changes of the position of individuals are big. It is visible in Figure 2(b) as a high oscillation. For f∈−2,0 and especially for f→0, the changes are smaller, and the algorithm systematically scans the solution space. It is visible in Figure 2(c) as intermediate values between the greatest and the smallest values of the response.
Forf<−4 and f>0, the eigenvalues have only a real part, and at least one of them lies outside the unit circle ∣z1,2∣>1; the algorithm is unstable. The samples of the response of the algorithm for f=−4.2 and f=0.2 are presented in Figures 2(a) and 2(d), respectively. For the frequency f smaller than f<−4, the response has an oscillatory character with amplitude exponential growing; the faster, the farther from the border frequency f=−4. For the positive frequency f>0, the response is aperiodic and growing exponentially, the faster, the farther from the f>0. For already shown points, lying close to the limit values, the value of the amplitude is near 103 after only 15 iterations. This causes that the individuals often are going beyond search space, in the case of constrained optimization. The procedure of repair infeasible individuals becomes important. This procedure is not predefined in the algorithm and strongly depends on the designer of the algorithm. The simple replacing by the limit value is often used. It is always a type of heuristic, and it causes that the echolocation does not work correctly. The procedure of mutation of the best individuals dominates in that algorithm.Figure 1
The eigenvalues of the matrixA on the z-plane.Figure 2
The responses of the algorithm with eigenvalues from Figure1 to initial conditions.
(a)
f=−4.2
(b)
f=−3.2
(c)
f=−0.8
(d)
f=0.2


## 5. The Example Experiments and Discussion

The size of the population of BA was equalN=40; the largest number of function evaluations for each call of BA has been set as 50·103. Step 7 of the code of BA in Pseudocode 1 (“Select a solution among the best solutions”) is not predefined and can be used any, like a roulette wheel, stochastic universal sampling or other. The first one of the above-mentioned methods was used during the experiments. Four benchmark functions, presented in Table 1, were used in order to check and analyse the efficiency of the BA. They are used as a quality function J during searching for global minimum value by BA. The BA was run 50 times for each benchmark function. The limitation of search space was presented in Table 1. This table also includes the values of parameters, the loudness L and emission rate r, used during experiments.Table 1
Benchmark functions used for calculating the quality functionJ and analysing the efficiency of BA with specified parameters (the loudness L and emission rate r) used during experiments.
Function
Search range
Type of function
min
LoudnessL
Emission rater
Sphere
[−5.12–5.12]
Unimod.
0
0.01
0.1
Griewank
[−600–600]
Multim.
0
4
0.1
Ackley
[−32,768–32,768]
Multim.
0
0.7
0.6
Schwefel
[−500–500]
Multim.
0
1.1
0.1The range of variability of frequencyf was divided into three types C presented in Table 2. The frequency was determined by (1) using the lower and upper bounds appropriately for each range and value of parameter B. In the second type of constraints, the mean value of frequency is equal to zero, and the only range of variability of frequency f is changed which influences the variance of the distribution function. The first and third types of constraints lie fully on the stable or unstable region, and both the mean value and variance are changing.Table 2
The type of constraints of the frequencyf as a function of the absolute value B (B>0).
Type of constraintsC
Lower boundfmin
Upper boundfmax
I
−B
0
II
−B
B
III
0
BThe main influences of frequencyf on the stability of algorithms can be seen as a number of individuals out of the boundaries of the seeking area. The problem of crossing the border can exist in both, for the location and speed. The per cent of new individuals with location and velocity outside the permitted area for the Griewank function is presented in Figures 3(a) and 3(b). Similar figures for the Schwefel and the Ackley functions were presented in [22]. A low absolute value of B gives generally less per cent of new individuals crossing the bounds. The number of individuals crossing the bounds for the stable area is smaller even than for the frequency with a mean value equal to zero, i.e., for the type II of the constraints of frequency. In this second case, the variation of the sign of the frequency f results in balancing of positive and negative values of frequency and consequently the lower number of new individuals with speed and position out of the permissible value. Figure 4 presents the mean value of quality Jmin as a function of the lower fmin and upper fmax bounds of frequency. The best results are obtained for the negative frequency fk∈−2,0, which fulfils the stability condition f∈–4,0. The best mean value of the frequency equals fmean=−1. Table 3 presents, obtained during experiments, the optimal mean values of frequency and the best mean normalized value of the fitness function for different bounds of frequency. The mean fitness functions Jmean for all types of constraints C of frequency, from Table 2, are normalized using the best of them Jmeanbest as a normalizing factor:
(36)Jmean¯C=JmeanCJmeanbest.Figure 3
The per cent of new individuals with (a) location and (b) velocity, outside the permitted area for the Griewank function.
(a)
(b)Figure 4
The mean value of the qualityJmean of the algorithm as a function of the absolute value of the lower and upper bounds of frequency for Griewank function.Table 3
The mean normalized value of the performanceJmean¯ for the benchmark function for different bounds of frequency and the best found mean values of the frequency.
Frequency symmetry
The best mean value of frequencyfmean
Negative
Symmetric
Positive
−B0
−BB
0B
Sphere
1.00
1.06E+05
1.16E+06
−4
Griewank
1.00
1.50
2.26
−1
Ackley
1.00
2.33
3.74
−1
Schwefel
1.00
2.67E+08
3.96E+09
−2.5The best solutions are bolded in the table. The Sphere and the Schwefel functions are very sensitive to the value of frequency and the stability of BA. The Ackley and the Griewank functions have the best solution for negative frequency, but the dominance of these frequencies is not significant.


## 6. Conclusions

In the paper, we described the dynamic of BA by the state-space form. The analyses of observability and controllability allow reducing the dynamics of BA from third to second order. The Lyapunov stability theory and Sylvester’s criterion applied to the stochastic dynamic of BA determine the condition for asymptotic stability and convergence to the equilibrium point. We also made the analysis of the location of eigenvalues of the state transition matrix and its influence on the response of the algorithm. Presented in the paper results can be used during the process of design and tuning the BA. As an illustrative example, four benchmark functions were used, with particular emphasis on the behaviour of BA used for Griewank function.



---

*Source: 9837462-2018-12-03.xml*