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

The main problemwith 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.


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][2][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][14][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][17][18] and the Lyapunov function [19][20][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 Sys-Tems 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 Section 2, 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.

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) The i-th bat is at position x i and flies randomly with a velocity v i .It emits an audio signal with a variable frequency between f min , f max , a varying wavelength λ i , and loudness A i 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) L 0 to the smallest constant value L min Each of the artificial bats in the k-th step has a position vector x k i , velocity vector v k i , and frequency vector f k i 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 f k i , which is drawn uniformly from the range f min , f max : where β k ∈ 0, 1 is a random vector with a uniform distribution.The velocity of the i-th bat in the k-th step v k i depends on the position of the current global best solution achieved so far x k bsf : The new position of the bat and thus the new solution of the problem follow from his earlier position and velocity: The local search procedure is also used.A new solution for a bat is generated locally using current best solution and local random walk: where ε ∈ −1, 1 is a random number with a uniform distribution, while L k 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 loudness L and emission rate r, updated as follows: 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 Pseudocode 1.

Lyapunov Stability Theory
The basic theorems for system stability are the Lyapunov stability theorems [23][24][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 point x = 0 of the stochastic discrete-time system, defined by the statespace equation: where x k is a state vector at time k and A k ∈ R n×n is a nonsingular matrix with stochastic values.The equilibrium point is asymptotically stable if there is a nonnegative scalar Lyapunov function V x k defined as where P is a positive definite symmetric matrix, with V 0 = 0, which satisfies that the expected value of changes of the Lyapunov function E ΔV is greater than zero: We can write it as which after simple transformation, using ( 7) and (10), leads to a useful formula: 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]: The symmetric matrix P 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.Updating velocities and locations (eq.( 2,3)) 6.
If rand > r k i 7.
Select a solution among the best solutions 8.
Generate a new solution by flying randomly 11. If Accept the new solutions 13.
Reduce By assuming the velocity of the bat v k i , its position x k i and loudness L k i , as state variables T and the position of the current global best solution achieved so far x k bsf as an input u k , the dynamics of the BA described by ( 2), (3), and ( 5) can be presented in the state-space form: where the locus of the i-th bat is treated as an output y k i .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, the i-th dimension of (15), updates independently from the others; thus, without losing the generality, the analysis of the algorithm can be reduced to the onedimensional case.Therefore, consequently to the general form of the state-space:

16
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: 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 frequency f k equals its expected value f k = E f k = f .Theorem 3 (see [30]).The system is completely observable if any initial state vector x t 0 can be reconstructed by examining the output of the system y t over the finite period of time from t 0 to t f .Then the system is completely observable if and The rank of the matrix of observability for BA dynamics (15) is equal: which means the system is observable.
Theorem 4 (see [23]).The system is completely controllable if there exists a control signal u t defined over a finite interval 0 ≤ t ≤ t f , which can force system states x t f to any desired value.Then the system is completely controllable if and only if the set of vectors B AB … A n−1 B is literary independent, i.e., rank B AB … A n−1 B = n.The rank of the matrix of controllability for BA dynamics described by ( 15) is equal: 4 Complexity The obtained rank of matrix A 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: The uncontrollable state-variable is the loudness L k , but if we assume α ∈ 0, 1 , then L k 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 f k and present the statespace form as 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 = 0 0 T : then state-space take the form According to Theorem 1, 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 with c 1 , c 2 > 0. The symmetric matrix P equals Substituting the matrices P and Q into relation (12), we obtain which leads to the system of relations: Assuming f = E f k , we can easily calculate that Matrix P must be positive definite.According to Sylvester's criterion in Theorem 2, the successive principal minors must be positive, which leads to equations where Considering the quadratic polynomial ap 1 2 + bp 1 + c, we obtain

31
and well-known relations We will look for a range of parameter f 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 frequency f ∈ −4, 0 .Assuming the big enough value of p 11 > λ 1,2 and p 11 > 0, we obtain that BA is stable in this range of f .
(II) a < 0 The condition p 11 > 0 can be satisfied only when Δ > 0, which is met for f < 0 5.We need to examine two conditions: (ii) 0 < f < 0 5 For this range of f , 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 range f ∈ −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 of A are calculated from the characteristic equation defined as where I = 1 0 0 1 , or they are defined as the roots of the denominator of the transfer function (20) equivalently, which leads to the equation The roots are equal: Figure 1 presents the eigenvalues of matrix A on the zplane.For f ∈ −4, 0 , the eigenvalues are complex, and their absolute value equals one |z 1,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.
For f <−4 and f > 0, the eigenvalues have only a real part, and at least one of them lies outside the unit circle |z 1,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 10 3 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.6 Complexity

The Example Experiments and Discussion
The size of the population of BA was equal N = 40; the largest number of function evaluations for each call of BA has been set as 50•10 The range of variability of frequency f 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 on the stable or unstable region, and both the mean value and variance are changing.
The main influences of frequency f 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 J min as a function of the lower  7 Complexity f min and upper f max bounds of frequency.The best results are obtained for the negative frequency f k ∈ −2, 0 , which fulfils the stability condition f ∈ -4, 0 .The best mean value of the frequency equals f mean = −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 J mean for all types of constraints C of frequency, from Table 2, are normalized using the best of them J mean best as a normalizing factor: The 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.

Conclusions
In the paper, we described the dynamic of BA by the statespace form.The analyses of observability and controllability     8 Complexity 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.

Figure 1 :
Figure 1: The eigenvalues of the matrix A on the z-plane.

Figure 2 :
Figure 2: The responses of the algorithm with eigenvalues from Figure 1 to initial conditions.

Figure 3 :
Figure 3: The per cent of new individuals with (a) location and (b) velocity, outside the permitted area for the Griewank function.

Figure 4 :
Figure 4: The mean value of the quality J mean of the algorithm as a function of the absolute value of the lower and upper bounds of frequency for Griewank function.
3. 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 abovementioned methods was used during the experiments.Four benchmark functions, presented in Table1, 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 Table1.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 function J and analysing the efficiency of BA with specified parameters (the loudness L and emission rate r) used during experiments.

Table 2 :
The type of constraints of the frequency f as a function of the absolute value B (B > 0).

Table 3 :
The mean normalized value of the performance J mean for the benchmark function for different bounds of frequency and the best found mean values of the frequency.