Bilel Elbetch
email: bilel.elbetch@univ-saida.dz

Tounsia Benzekri
email: tbenzekri@usthb.dz

Daniel Massart
email: daniel.massart@umontpellier.fr

Tewfik Sari
email: tewfik.sari@inrae.fr

The multi-Patch logistic equation

Keywords: 2010 Mathematics Subject Classification. Primary: 37N25, 92D25; Secondary: 34D23, 34D15 Fragmented habitat, Migration rate, Logistic growth, Carrying capacity, Perfect mixing

1. Introduction. Population dynamics is a wide field of mathematics, which contains many problems, for example fragmentation of population and the effect of migration in the general dynamics of population. Bibliographies can be found in the work of Levin [START_REF] Levin | Dispersion and population interactions[END_REF][START_REF] Levin | Spatial patterning and the structure of ecological communities[END_REF] and Holt [START_REF] Holt | Population dynamics in two patch environments: some anomalous consequences of an optimal habitat distribution[END_REF]. There are ecological situations that motivate the representation of space as a finite set of patches connected by migrations, for instance an archipelago with bird population and predators. It is an example of insular bio-geography. The standard question in this type of biomathematical problems, is to study the effect of migration on the general population dynamics, and the consequences of fragmentation on the persistence or extinction of the population.

Freedman and Waltman [START_REF] Freedman | Mathematical Models of Population Interactions with Dispersal I: Stabilty of two habitats with and without a predator[END_REF] were among the first to consider the case of two patches  

 dx1 dt = r 1 x 1 1 -x1 K1 + β(x 2 -x 1 ), dx2 dt = r 2 x 2 1 -x2 K2 + β(x 1 -x 2 ), (1)
where x i is the population in patch i, the parameters r i and K i are respectively the intrinsic growth rate and the carrying capacity patch i, and β is the migration rate assumed to be symmetric in both directions. They analyzed the model in the case of perfect mixing (β → +∞) and showed that the total equilibrium population x * 1 (β) + x * 2 (β) can be greater than the sum of the carrying capacities K 1 + K 2 . The 2-patch model [START_REF] Arditi | The perfect mixing paradox and the logistic equation: Verhulst vs[END_REF] was extensively studied by many authors, among them DeAngelis et al. [START_REF] Deangelis | Persistence and stability of seed-dispersel species in a patchy environment[END_REF], Holt [START_REF] Holt | Population dynamics in two patch environments: some anomalous consequences of an optimal habitat distribution[END_REF], Freedman et al. [START_REF] Freedman | Mathematical Models of Population Interactions with Dispersal II: Differential Survival in a Change of Habitat[END_REF]. It is worth noting that the two-patch model (1) had never been completely analysed until the work of Arditi et al. [START_REF] Arditi | In dispersal always beneficial to carrying capacity? New insights from the multi-patch logistic equation[END_REF] which gave the complete study of the model [START_REF] Arditi | The perfect mixing paradox and the logistic equation: Verhulst vs[END_REF], i.e the comparison of X * T with K 1 +K 2 for all values of β. One distinguishing characteristic is the depiction of three situations for the position of x * 1 (β)+x * 2 (β) compared to K 1 +K 2 . Under some conditions, fragmentation leads to a total equilibrium population strictly greater or smaller than the sum of the carrying capacities, see Fig. 1 in [START_REF] Arditi | In dispersal always beneficial to carrying capacity? New insights from the multi-patch logistic equation[END_REF].

DeAngelis and Zhang [START_REF] Deangelis | Effects of dispersal in a non-uniform environment on population dynamics and competition: a patch model approach[END_REF], DeAngelis et al. [START_REF] Deangelis | Effects of diffusion on total biomass in heterogeneous continuous and discrete-patch systems[END_REF] and Zhang et al. [START_REF] Zhang | Effects of dispersal on total biomass in a patchy, heterogeneous system: analysis and experiment[END_REF] have considered the model

dx i dt = r i x i 1 - x i K i + β(x i-1 -2x i + x i+1 ), i = 1, • • • , n (2)
where we denote x 0 = x n and x n+1 = x 1 , allowing the patches to join in a circle so that the same relationships hold between x i , x i-1 and x i+1 for all values of i. An interesting result of ( 2) is that in the case r i = K i , for i = 1, • • • , n, for any β > 0, the total population at steady state satisfies n i=1

x * i (β)

> n i=1 K i . ( 3
)
Our aim is to extend some of the results in [START_REF] Arditi | In dispersal always beneficial to carrying capacity? New insights from the multi-patch logistic equation[END_REF][START_REF] Arditi | Asymmetric dispersal in the multi-patch logistic equation[END_REF][START_REF] Deangelis | Effects of dispersal in a non-uniform environment on population dynamics and competition: a patch model approach[END_REF][START_REF] Deangelis | Effects of diffusion on total biomass in heterogeneous continuous and discrete-patch systems[END_REF] to the more general case of the n-patch model:

dx i dt = r i x i 1 - x i K i + β n j=1,j =i γ ij (x j -x i ) , i = 1, • • • , n. ( 4
)
where n is the number of the patches in the system. The parameter β represents the dispersal rate of the population; γ ij = γ ji ≥ 0 denotes the flux between patches j and i, for j = i. Note that if γ ij > 0 there is a flux of migration between patches i and j and if γ ij = 0 there is no migration. In the case where β = 0, there is no migration , and the system (4) admits (K 1 , • • • , K n ) as a non trivial equilibrium point, which furthermore is globally asymptotically stable (GAS). The problem is whether or not, the equilibrium continues to be positive and GAS for any β > 0 and whether or not, the total population at equilibrium is greater than the sum of the carrying capacities. The paper is organized as follows. In Section 2, the mathematical model of n patches, and some notation, are introduced. In Section 3, we give some results on the comparison of the total equilibrium population with the sum of the carrying capacities. In Section 4, the behavior of the model is studied when the migration rate tends to infinity. In Section 5, we consider the three-patch model and by numerical simulations we prove the existence of a new behavior of the total equilibrium population. In section 6, we consider the multi-patch logistic growth, coupled by asymmetric migration terms. In Appendix A, we give some properties of the total equilibrium population. In Appendix B, we recall some results which are used in our proofs.

In all the paper, we denote by: R + the non negative real numbers. Given a vector u ∈ R n let u T denotes its transpose. We denote by Re and |.| the real part and the module of a complex number respectively.

2. Mathematical model. We consider the model of multi-patch logistic growth, coupled by symmetric migration terms (4). This system of differential equations can be written:

dx i dt = r i x i 1 - x i K i + β n j=1 γ ij x j , i = 1, • • • , n (5)
where

γ ii = - n j=1,j =i γ ij , i = 1, • • • , n (6)
denotes the outgoing flux of patch i. We denote by Γ the matrix

Γ := (γ ij ) n×n . (7)
Its columns sum to 0 since the matrix Γ is symmetric and the diagonal elements γ ii are defined by [START_REF] Deangelis | Persistence and stability of seed-dispersel species in a patchy environment[END_REF] in such a way that each row sums to 0. The matrix

Γ 0 := Γ -diag(γ 11 , . . . , γ nn ) (8)
which is the same as Γ, except that the diagonal elements are 0, is called the connectivity matrix. It is the adjacency matrix of the weighted directed graph G, which has exactly n vertices (the patches), and there is an arrow from patch j to patch i precisely when γ ij > 0, with weight γ ij assigned to the arrow.

The model (1) studied in [START_REF] Arditi | In dispersal always beneficial to carrying capacity? New insights from the multi-patch logistic equation[END_REF][START_REF] Deangelis | Persistence and stability of seed-dispersel species in a patchy environment[END_REF][START_REF] Freedman | Mathematical Models of Population Interactions with Dispersal II: Differential Survival in a Change of Habitat[END_REF][START_REF] Freedman | Mathematical Models of Population Interactions with Dispersal I: Stabilty of two habitats with and without a predator[END_REF][START_REF] Holt | Population dynamics in two patch environments: some anomalous consequences of an optimal habitat distribution[END_REF] corresponds to the case where the migration rate γ 12 = γ 21 > 0 can be normalized to 1, by replacing β by β/γ 12 . The model [START_REF] Arditi | In dispersal always beneficial to carrying capacity? New insights from the multi-patch logistic equation[END_REF] studied in [START_REF] Deangelis | Effects of dispersal in a non-uniform environment on population dynamics and competition: a patch model approach[END_REF][START_REF] Deangelis | Effects of diffusion on total biomass in heterogeneous continuous and discrete-patch systems[END_REF][START_REF] Zhang | Effects of dispersal on total biomass in a patchy, heterogeneous system: analysis and experiment[END_REF] corresponds to the case where

γ 1n = γ n1 = γ i,i-1 = γ i-1,i = 1 for 2 ≤ i ≤ n and γ ij = 0 otherwise. ( 9
)
We have the following result:

Proposition 2.1. The domain Ω = {(x 1 , . . . , x n ) ∈ R n /x i ≥ 0, i = 1, . . . , n.
} is positively invariant for the system (4).

Proof. Assume that x j ≥ 0 for all j and there exists i such that x i = 0. We have

dx i dt = j =i γ ij x j ≥ 0.
Hence, on the boundary of Ω, the vector field associated to (4) either is tangent to the boundary of Ω, or points inward Ω. According to [START_REF] Smith | The Theory of the Chemostat: Dynamics of Microbial Competition[END_REF]Proposition B.7,page 267], no trajectory comes out of Ω. Therefore, Ω is positively invariant for (4).

Proposition 2.2. Assume that the matrix Γ := (γ ij ) n×n (or equivalently, the connectivity matrix Γ 0 ) is irreducible. The model (4) has a unique positive equilibrium point which is GAS in the positive cone R n + \ {0}.

Proof. The result follows from [START_REF] Takeuchi | Cooperative systems theory and global stability of diffusion models[END_REF].

Remark 2.3. The matrix Γ being irreducible means that the set of patches cannot be partitioned into two nonempty disjoint subsets, I and J, such that there is no migrations between a patch in subset I and a patch in subset J. The matrix Γ is assumed to be irreducible throughout the rest of the paper. Therefore species can reach any patch from any patch either directly or through other patches.

3. Comparison of the total equilibrium population with the sum of carrying capacities. In all of this work, the GAS equilibrium of the system (4), whose existence is shown in Proposition 2.2, is denoted by E * (β) = (x * 1 (β), . . . , x * n (β)). The equilibrium point E * (β) is the solution of the algebraic system:

r i x i 1 - x i K i + β n j=1,j =i γ ij (x j -x i ) = 0, i = 1, • • • , n. (10)
The sum of these equations shows that E * (β) satisfies the following equation

n i=1 r i x i 1 - x i K i = 0.
Therefore E * (β) belongs to the ellipsoid

E n-1 = x ∈ R n : Θ(x) := n i=1 r i x i 1 - x i K i = 0 . (11)
Note that this ellipsoid is independent of the migration terms β and γ ij .

Our aim is to compare the total equilibrium population

X * T (β) = x * 1 (β) + • • • + x * n (β), (12)
with the sum of carrying capacities K 1 + . . . + K n , when the rate of migration β varies from zero to infinity. Let us begin with some interesting particular cases of the system (4).

3.1. Equal growth rates. In the next proposition, we show that if the growth rates are equal in all patches, then the total equilibrium population is always smaller than the sum of the carrying capacities:

Proposition 3.1. If r 1 = • • • = r n ,
then the total equilibrium population, defined by [START_REF] Holt | Population dynamics in two patch environments: some anomalous consequences of an optimal habitat distribution[END_REF] satisfies

X * T (β) ≤ n i=1 K i , for all β ∈ [0, ∞[. Proof.
The equation of the tangent space to the ellipsoid E n-1 , defined by [START_REF] Gantmacher | The Theory of Matrices[END_REF]

, at point A = (K 1 , . . . , K n ) is given by n i=1 (x i -K i ) ∂Θ ∂x i (A) = 0, ( 13
)
where Θ is given by [START_REF] Gantmacher | The Theory of Matrices[END_REF]. Since ∂Θ ∂xi (A) = -r i , (13) can be written as follows:

n i=1 r i (x i -K i ) = 0. ( 14
)
If we take r 1 = • • • = r n in [START_REF] Levin | Dispersion and population interactions[END_REF], we get that the equation of the tangent plane to

E n-1 at the point A is n i=1 x i = n i=1 K i .
By the convexity of E n-1 , any point of E n-1 lies in the half-space defined by the inequation

n i=1 x i ≤ n i=1 K i . Therefore E * (β) satisfies n i=1 x * i (β) ≤ n i=1 K i for all β ≥ 0.
This completes the proof of the proposition.

Equal carrying capacities.

In the next proposition, we show that if the carrying capacities are equal in all patches, then the total equilibrium population does not depend on the migration rate and is always equal to the sum of the carrying capacities:

Proposition 3.2. If K 1 = • • • = K n =: K, then X * T (β) = nK for all β ≥ 0. Proof. Suppose that K 1 = • • • = K n =: K.
Then the equations [START_REF] Freedman | Mathematical Models of Population Interactions with Dispersal I: Stabilty of two habitats with and without a predator[END_REF] giving the equilibrium point are satisfied for x i = K. Therefore E * (β) = (K, . . . , K) does not depend on β and X * T (β) = nK. 3.3. Equal ratios r i /K i . In the next proposition, we show that if the ratios α i = r i /K i are equal in all patches, then the total equilibrium is always larger than the sum of the carrying capacities:

Proposition 3.3. If α 1 = . . . = α n =: α, then X * T (β) ≥ n i=1 K i , for all β ≥ 0. (15)
Moreover, if there exist i and j such

K i = K j , then X * T (β) > n i=1 K i , for all β > 0.
Proof. If α 1 = . . . = α n =: α, then, from (42) we have

X * T (β) = n i=1 K i + β α n i=1 n j=1,j =i γ ij (x * j (β) -x * i (β)) x * i (β) = n i=1 K i + β α j<i γ ij (x * j (β) -x * i (β)) x * i (β) + γ ji (x * i (β) -x * j (β)) x * j (β) = n i=1 K i + β α j<i γ ij (x * j (β) -x * i (β))x * j (β) + γ ji (x * i (β) -x * j (β))x * i (β) x * j (β)x * i (β) . Since the matrix Γ is symmetric, i.e γ ij = γ ij then X * T (β) = n i=1 K i + β α j<i γ ij (x * j (β) -x * i (β)) 2 x * j (β)x * i (β) ≥ n i=1 K i .
Equality can hold if and only if β = 0 or x * j (β) = x * i (β) for all j < i. Let us prove that if at least two patches have not the same carrying capacity then equality cannot hold for β > 0. Suppose that there exists β > 0 such that the positive equilibrium satisfies x * j (β) = x * i (β) for all j < i, then x * 1 (β) = . . . = x * n (β). Denote by x * > 0 their common value. Then all terms γ ij (x j -x i ) in [START_REF] Freedman | Mathematical Models of Population Interactions with Dispersal I: Stabilty of two habitats with and without a predator[END_REF] vanish, so that these equations become

r i x * 1 - x * K i = 0, i = 1, • • • , n.
Therefore K i = x * for all i, which gives a contradiction. Hence the equality in [START_REF] Levin | Spatial patterning and the structure of ecological communities[END_REF] can hold if and only if β = 0.

In [7, page 3092], DeAngelis et al. have shown this result in the particular case of the model ( 2), and

r i = K i , i = 1 • • • n, see (3). Case K < K n 0 J 0 J 1 J 2 r r n rn = r rn Kn = r K Case K > K n 0 J 2 J 1 J 0 r r n rn = r rn Kn = r K Figure 1
. Qualitative properties of model ( 4) when ( 16) holds. In J 0 , patchiness has a beneficial effect on total equilibrium population. This effect is detrimental in J 2 . In J 1 , the effect is beneficial for β < β 0 and detrimental for β > β 0 .

3.4. All patches but one are identical. We assume now that n -1 patches are identical, that is to say, they have the same carrying capacities and the same growth rates. We assume also that the flux of migration between the n-th patch and all these n -1 identical patches are equal. Hence we have the following conditions

K 1 = • • • = K n-1 =: K, r 1 = . . . = r n-1 =: r and γ n1 = . . . = γ nn-1 =: γ (16)
Under conditions [START_REF] Lu | Global asymptotic behavior in single-species discrete diffusion systems[END_REF], we show that the n-patch model behaves like a 2-patch model, that is, there are only three situations as depicted in Fig. 1. This figure is similar to Figure 2 of [START_REF] Arditi | In dispersal always beneficial to carrying capacity? New insights from the multi-patch logistic equation[END_REF]. This property holds regardless of the flux of migration between the identical n-1 patches, provided that the connectivity matrix Γ, defined by [START_REF] Deangelis | Effects of diffusion on total biomass in heterogeneous continuous and discrete-patch systems[END_REF], is irreducible. More precisely, define regions J 0 , J 1 and J 2 by:

If K < K n then        J 0 = (r, r n ) : r n ≥ Kn K r J 1 = (r, r n ) : Kn K r > r n > r J 2 = {(r, r n ) : r ≥ r n } (17) If K n < K then        J 0 = (r, r n ) : r n ≤ Kn K r J 1 = (r, r n ) : Kn K r < r n < r J 2 = {(r, r n ) : r ≤ r n } (18)
We have the following proposition: Proposition 3.4. Assume that the conditions (16) hold in the system (4). We consider the regions, denoted J 0 , J 1 and J 2 , depicted in Fig. 1, and defined by (17) when K < K n and by [START_REF] Othmer | A Continuum Model for Coupled Cells[END_REF]

when K > K n . Let α n = r n /K n and α = r/K. 1. If (r, r n ) ∈ J 0 then X * T (β) > (n -1)K + K n for all β > 0. 2. If (r, r n ) ∈ J 1 then X * T (β) > (n -1)K + K n for 0 < β < β 0 and X * T (β) < (n -1)K + K n for β > β 0 where β 0 is given by β 0 = r n -r γ(α -α n ) 1 1/α n + (n -1)/α . ( 19
)
3. If (r, r n ) ∈ J 2 then X * T (β) < (n -1)K + K n for all β > 0. Proof. Assume that the conditions ( 16) hold and K < K n . From item (2) in Proposition A.2 we deduce that

x * 1 (β) = • • • = x * n-1 (β) =: x * (β) < x * n (β). ( 20
)
Now, using Lemma A.1, the total equilibrium population becomes

X * T (β) = (n -1)K + K n + β   n-1 j=1 γ (x * (β) -x * n (β)) α n x * n + n-1 i=1 γ (x * n (β) -x * (β)) αx *   .
Therefore

X * T (β) = (n -1)K + K n + (n -1)βγ (x * (β) -x * n (β))(αx * (β) -α n x * n (β)) αα n x * (β)x * n (β) . (21) 1. If (r, r n ) ∈ J 0 , then α n ≥ α. Hence (αx * -α n x * n ) < α(x * -x * n ) < 0, which implies that (x * -x * n )(αx * -α n x * n ) > 0, and then X * T (β) > (n -1)K + K n for β > 0.
2. If (r, r n ) ∈ J 1 , then r n > r and α > α n . According to Equation ( 21), the equality

X * T (β) := (n -1)x * (β) + x * n (β) = (n -1)K + K n , is equivalent to β = 0 or αx * (β) -α n x *
n (β) = 0. Thus x * (β) and x * n (β) are the solutions of the linear system

αx * (β) -α n x * n (β) = 0, (n -1)x * (β) + x * n (β) = (n -1)K + K n , which admits a unique solution x * (β) = 1 α (n -1)K + K n (n -1)/α + 1/α n , x * n (β) = 1 α n (n -1)K + K n (n -1)/α + 1/α n . ( 22
) As (x * 1 (β), • • • , x * n (β)
) is an equilibrium point of the system (4), it satisfies the following system

       0 = rx * (β) 1 - x * (β) K + βγ (x * n (β) -x * (β)) , 0 = r n x * n (β) 1 - x * n (β) K n + (n -1)βγ (x * (β) -x * n (β)) , (23)
obtained from (10) by replacing x i = x * i (β) and using ( 16) and [START_REF] Takeuchi | Cooperative systems theory and global stability of diffusion models[END_REF]. Using the formulas [START_REF] Wasow | Asymptotic Expansions for Ordinary Differential Equations[END_REF], the system (23) has a unique solution β 0 given by [START_REF] Smith | The Theory of the Chemostat: Dynamics of Microbial Competition[END_REF], which is positive since r n > r and α > α n . According to Lemma A.3 we have:

dX * T dβ (0) = (n -1)γ(K -K n ) 1 r n - 1 r ,
which is positive since r < r n and K < K n . So:

• If β ∈ (0, β 0 ], then X * T (β) > (n -1)K + K n . • If β ∈ (β 0 , ∞[, then X * T (β) < (n -1)K + K n .
3. If (r, r n ) ∈ J 2 then r ≥ r n . We have:

αx * (β) -α n x * n (β) = r K x * (β) - r n K n x * n (β)
From item (2) in Proposition A.2 we deduce that x * (β)/K > 1 and 4. Perfect mixing. In this section our aim is to study the behavior of the system (4) for large migration rate, i.e when β → ∞. We need the following result: Lemma 4.1. The matrix Γ has rank n -1. Except 0 which is a simple eigenvalue of Γ, whose eigenvector is u = (1, . . . , 1) T , all other eigenvalues of Γ are negative.

x * n (β)/K n < 1. Therefore αx * (β) -α n x * n (β) > r -r n ≥ 0, which gives X * T (β) < (n -1)K + K n for β > 0. We use item 2-(b) of Proposition A.
Proof. Let s = max i=1,...,n (-γ ii ). The matrix B defined by

B = Γ + sI,
is non-negative and irreducible. Then by the Perron-Frobenius theorem [11, Theorem 2, page 53], the spectral radius of B is a simple eigenvalue, and it is the only eigenvalue of B which admits a positive eigenvector. Now observe that (1, . . . , 1) T is an eigenvector of B, and the corresponding eigenvalue is s. This proves that s = ρ(B) is a simple eigenvalue of B, and since s = ρ(B) any other eigenvalue λ of B satisfies λ < s. Therefore 0 is a simple eigenvalue of Γ = -sI + B, that is, the rank of Γ is n -1, and all other eigenvalues of Γ are negative.

We have the following result: Theorem 4.2. We have

lim β→+∞ E * (β) = n i=1 r i n i=1 r i /K i (1, . . . , 1) . Proof. Let α i = r i /K i . Denote E * (∞) = n i=1 r i n i=1 α i , . . . , n i=1 r i n i=1 α i .
Recall that the equilibrium point E * (β) is the unique positive solution of the algebraic system [START_REF] Freedman | Mathematical Models of Population Interactions with Dispersal I: Stabilty of two habitats with and without a predator[END_REF]. Therefore it is defined, in the positive cone, by the fact that it belongs to the ellipsoid E n-1 defined by [START_REF] Gantmacher | The Theory of Matrices[END_REF], and the equation F β (E * (β)) = 0, where

F β (x 1 , . . . , x n ) :=   r i β x i 1 - x i K i + n j=1,j =i γ ij (x j -x i )   i=1,••• ,n-1
.

The ellipsoid E n-1 is compact, so E * (β) has at least one limit point in E n-1 , when β goes to infinity. Besides, F β converges, uniformly on compact sets, to the linear map

R n -→ R n-1 (x 1 , . . . , x n ) T -→ Λ(x 1 , . . . , x n ) T
where Λ is the sub matrix of the matrix Γ defined by [START_REF] Deangelis | Effects of dispersal in a non-uniform environment on population dynamics and competition: a patch model approach[END_REF], obtained by dropping the last row of Γ. Hence any limit point (x 1 , . . . , x n ) of E * (β) satisfies Λ(x 1 , . . . , x n ) T = 0.

Since Λ is obtained from Γ by dropping the last row, and the rows of Γ sum to zero, the rank of Λ equals that of Γ, which is n -1 by Lemma 4.1.

Therefore, by the convexity of E n-1 , Λ(x 1 , . . . , x n ) T = 0 has exactly two solutions in E n-1 , which are the origin and E * (∞). To prove the convergence of E * (β) to E * (∞), it suffices to prove that the origin cannot be a limit point of E * (β).

Write E * (β) = (x 1 (β), . . . , x n (β)). We claim that for any β, there exists i such that x i (β) ≥ K i , which entails that E * (β) is bounded away from the origin.

Let

x i (β) = min{x 1 (β), • • • , x n (β)}.
We have

r i β x i (β) 1 - x i (β) K i + n j=1,j =i γ ij (x j -x i ) = 0 whence r i β x i (β) 1 - x i (β) K i ≤ 0
and since x i (β) cannot be negative or 0, we have

x i (β) ≥ K i .
As a corollary of the previous theorem we obtain the following result which describes the total equilibrium population for perfect mixing: Corollary 4.3. We have

X * T (+∞) = lim β→+∞ n i=1 x * i (β) = n n i=1 r i n i=1 r i /K i . ( 24
)
Proposition 4.4.

If α 1 = • • • = α n =: α, then X * T (+∞) = n i=1 K i . Proof. We use Equation (24) for α 1 = • • • = α n =: α.
Remark 4.5. Notice that the formula [START_REF] Zhang | Effects of dispersal on total biomass in a patchy, heterogeneous system: analysis and experiment[END_REF] shows that the total equilibrium population does not depend on the flux of migration γ ij . This formula was obtained in [START_REF] Deangelis | Effects of diffusion on total biomass in heterogeneous continuous and discrete-patch systems[END_REF]Theorem B.1] for the model (2) corresponding to the particular case (9) of the migration flux. It was also obtained for the 2-patch case by Freedman and Waltman [10, Theorem 3.1] and Arditi et al. [2, Formula (A.13)].

We can use the theory of singular perturbations to obtain a better understanding of the behaviour of the system in the case of perfect mixing. Theorem 4.6. Let (x 1 (t, β), . . . , x n (t, β)) be the solution of the system (4) with initial condition

(x 10 , • • • , x n0 ) satisfying x i0 ≥ 0 for i = 1 • • • n. Let Y (t) be the solution of the logistic equation dX dt = rX 1 - X nK , where r = n i=1 r i n and K = n i=1 r i n i=1 r i /K i , (25)
with initial condition Y (0) = n i=1 x i0 . Then, when β → ∞, we have

n i=1 x i (t, β) = Y (t) + o(1), uniformly for t ∈ [0, +∞) (26)
and, for any t 0 > 0, we have

x i (t, β) = Y (t) n + o(1), i = 1, . . . , n, uniformly for t ∈ [t 0 , +∞). ( 27
)
Proof. Since Γ is symmetric with real coefficients, there exists an orthogonal matrix P such that

P -1 ΓP = D
where D is the diagonal matrix whose diagonal elements are the eigenvalues

λ 1 ≤ • • • ≤ λ n-1 < λ n = 0
of Γ, see Lemma 4.1. We rewrite the system (5) in vector form:

dx dt = f (x) + βΓx, (28)
where

x = (x 1 , • • • , x n ) T and f (x) = (r 1 x 1 (1 -x 1 /K 1 ), • • • , r n x n (1 -x n /K n )) T . ( 29
)
Using the variables y = P -1 x, the system (28) becomes:

dy dt = P -1 f (P y) + βDy.
This system of n differential equations can be written

     dy i dt = g i (y) + βλ i y i , i = 1, • • • , n -1 dy n dt = g n (y), (30)
where the components g i , for i = 1, • • • , n, are defined by

P -1 f (P y) = (g 1 (y), • • • , g n (y)) T . (31)
When β → ∞, (30) is a slow-fast system, with one slow variable, y n , and n -1 fast variables,

y i for i = 1 • • • n -1.
According to Tikhonov's theorem [START_REF] Lobry | On Tykhonov's theorem for convergence of solutions of slow and fast systems[END_REF][START_REF] Tikhonov | Systems of differential equations containing small parameters in the derivatives[END_REF][START_REF] Wasow | Asymptotic Expansions for Ordinary Differential Equations[END_REF] we consider the dynamics of the fast variables in the time scale τ = βt. One obtains

dy i dτ = 1 β g i (y) + λ i y i i = 1, • • • , n -1.
In the limit β = ∞, we find the fast dynamics

dy i dτ = λ i y i , i = 1, • • • , n -1. ( 32
)
The slow manifold, which is the equilibrium point of the fast dynamics (32), is unique and is given by:

y i = 0 for i = 1 • • • n -1. ( 33
)
Since λ i < 0 for i = 1, . . . , n -1, the slow manifold is GAS for (32). So, the theorem of Tikhonov ensures that after a fast transition toward the slow manifold, the solutions of (30) are approximated by the solutions of the reduced model which is obtained by replacing y i , i = 1, • • • , n -1, by (33) into the dynamics of the slow variable y n , given by the last equation in (30). One obtains:

dy n dt = g n ((0, • • • , 0, y n ) T ). ( 34
)
Let us compute g n ((0, • • • , 0, y n ) T ) which is the last component of the vector (31), given by P -1 f P (0, • • • , 0, y n ) T . We first notice that

P (0, • • • , 0, y n ) T = 1 √ n (y n , • • • , y n ) T ,
since it is obtained by multiplying by y n the last column of P , which is simply the normalized eigenvector u/ √ n of Γ, corresponding to the eigenvalue λ n = 0, see Lemma 4.1. On the other hand, since P is orthogonal, the last row of P -1 = P T is equal to u T / √ n. Therefore, using the definition (29) of the vector f , one has

g n ((0, • • • , 0, y n ) T ) = 1 √ n n i=1 r i y n √ n 1 - y n √ nK i = ry n 1 - y n √ nK ,
where r and K are defined in (25). Hence the reduced equation ( 34) is

dy n dt = ry n 1 - y n √ nK . ( 35
)
Since (35) admits √ nK as a positive equilibrium point, which is GAS in the positive axis, the approximation given by Tikhonov's theorem holds for all t ≥ 0 for the slow variable y n and for all t ≥ t 0 > 0 for the fast variables y i , i = 1, • • • , n -1, where t 0 is as small as we want.

Using

y n = X/ √ n, where X = n i=1 x i , one obtains that n i=1 x i (t, β) is ap- proximated by a solution of the equation dX dt = rX 1 - X nK . (36)
Therefore, let Y (t) be the solution of (36) of initial condition Y (0) = n i=1 x i0 , then, when β → ∞, we have the approximation (26). To prove (27), we observe that the last column of P is u/ √ n, so for all i, x i (t, β) is y n (t, β)/ √ n, plus some linear combination of the y i (t, β) for i < n, all of which converge to zero, uniformly in t.

In the case of perfect mixing, the approximation (26) shows that the total population behaves like the unique logistic equation (25) and then, when t and β tend to ∞, the total population x i (t, β) tends toward nK = n r i / α i , where α i = r i /K i , as stated in Corollary 4.3. The approximation (27) shows that, with the exception of a thin initial boundary layer, where the population density x i (t, β) quickly jumps from its initial condition x i0 to the average Y (0)/n, each patch of the n-patch model behaves like the single logistic equation

dx dt = rx 1 - x K , where r = 1 n n i=1 r i and K = n i=1 r i n i=1 r i /K i . ( 37
)
Hence, when t and β tend to ∞, the population density x i (t, β) tends toward K = r i / α i , as stated in Theorem 4.2.

Remark 4.7. The single logistic equation (37) gives an approximation of the population density in each patch in the case of perfect mixing. The intrinsic growth rate r is the arithmetic mean of the local intrinsic growth rates r i and the carrying capacity K is the weighted harmonic mean of the local carrying capacities K i with weights the growth rates r i . We point out the similarity between our expression for the carrying capacity in the limit β → ∞, and the expression obtained in spatial homogenization, see e.g [START_REF] Yurk | Homogenization techniques for population dynamics in strongly heterogeneous landscapes[END_REF][START_REF] Othmer | A Continuum Model for Coupled Cells[END_REF][START_REF] Zaker | The Effect of Movement Behavior on Population Density in Patchy Landscapes[END_REF].

Remark 4.8. Notice that if we use the r-α formalism for the logistic equation, instead of the r-K formalism, where α = r/K is the parameter quantifying intraspecific competition, then the n-patch model ( 4) becomes

dx i dt = r i x i -α i x i 2 + β n j=1,j =i γ ij (x j -x i ) , i = 1, • • • , n.
The perfect mixing approximation (37) of each population x i becomes

dx dt = rx -αx 2 , where r = 1 n n i=1 r i and α = 1 n n i=1 α i ,
which is a single logistic equation whose intrinsic growth rate r and intraspecific competition parameter α are the arithmetic means of the local r i and α i respectively.

For more information and complements on the comparison of the r-K and r-α formalism, the reader is referred to [START_REF] Arditi | The perfect mixing paradox and the logistic equation: Verhulst vs[END_REF].

5.

Comparison of the n-patch system with a two-patch system. Our aim in this section is to show that when n ≥ 3 the behavior of the n-patch logistic equation is richer and more intricate than in the case n = 2.

5.1. Two-patch model. As noticed in Remark 4.5, the formula (24) for perfect mixing was obtained by Arditi et al. [START_REF] Arditi | In dispersal always beneficial to carrying capacity? New insights from the multi-patch logistic equation[END_REF]. These authors used the theory of singular perturbations. They showed that, if (x 1 (t, β), x 2 (t, β)) is the solution of (1), with initial condition (x 10 , x 20 ), then, when β → ∞, the total population x 1 (t, β) + x 2 (t, β) is approximated by Y (t), the solution of the logistic equation

dX dt = rX 1 - X 2K , where r = r 1 + r 2 2 and K = r 1 + r 2 r 1 /K 1 + r 2 /K 2 , (38)
with initial condition Y (0) = x 10 + x 20 . Therefore the total population behaves like the unique logistic equation given by (38). In addition, one obtains the following property: with the exception of a small initial interval, the population densities x 1 (t, β) and x 2 (t, β) are both approximated by Y (t)/2, see [2, Proposition 3]. Therefore, this approximation shows that, when t and β tend to ∞, the population density x i (t, β) tends toward r1+r2 α1+α2 , as stated in Theorem 4.2 and in addition, x i (t, β) quickly jumps from its initial condition x i0 to the average (x 10 + x 20 )/2 and then is very close to Y (t)/2. In Theorem 4.6 we extended this result to the n-patch system.

It was shown by Arditi et al. [2, Proposition 2] that only three situations can occur: the case where the total equilibrium population is always greater than the sum of carrying capacities, the case where it is always smaller, and a third case, where the effect of migration is beneficial for lower values of the migration coefficient β and detrimental for the higher values. More precisely, it was shown in [START_REF] Arditi | In dispersal always beneficial to carrying capacity? New insights from the multi-patch logistic equation[END_REF] that, if n = 2, the following trichotomy holds

• If X * T (+∞) > K 1 + K 2 then X * T (β) > K 1 + K 2 for all β > 0. • If d dβ X * T (0) > 0 and X * T (+∞) < K 1 + K 2 , then there exists β 0 > 0 such that X * T (β) > K 1 + K 2 for 0 < β < β 0 , X * T (β) < K 1 + K 2 for β > β 0 and X * T (β 0 ) = K 1 + K 2 . • If d dβ X * T (0) < 0, then X * T (β) < K 1 + K 2 for all β > 0.
Therefore, the condition X * T (β) = K 1 + K 2 holds only for β = 0 and at most for one positive value β = β 0 . The value β 0 exists if and only if d dβ X * T (0) > 0 and

X * T (+∞) < K 1 + K 2 .
In Proposition 3.4 we extended this result to the particular n-patch when all patches, but one, are identical and have the same migration rate to and from the last patch. Our proof is adapted from the proof of [2, Proposition 2]. In the next section we show that for three or more patches the total equilibrium population as a function of diffusion strength can have more types of qualitative behavior than the two-patch system. In particular there is more than one value of β for which the total equilibrium population can equal the sum of the carrying capacities.

β K 1 K 2 + K 3 X * T Figure 2.
Total equilibrium population X * T of the system (4) (n = 3) as a function of migration rate β. The parameter values are given in Table 5.2. 5.2.

β K 1 + K 2 + K 3 K 1 + K 2 + K 3 X * T X * T
β K 1 + K 2 + K 3 X * T Figure 4.
Total equilibrium population X * T of the system (4) (n = 3) as a function of migration rate β. The parameter values are given in Table 5.2. Table 1. The numerical values of parameters for the logistic growth function of the model ( 4), with n = 3, used in Fig. 2, 3, 4. All migration coefficients satisfy γ ij = 1. The derivative of the total equilibrium population at β = 0 is computed with Eq. ( 48), and the perfect mixing total equilibrium population X * T (+∞) is computed with Eq. [START_REF] Zhang | Effects of dispersal on total biomass in a patchy, heterogeneous system: analysis and experiment[END_REF]. 5.2. Three-patch model. In the numerical simulations we take n = 3 and γ ij = 1 for all i, j = 1, 2, 3. We show that we can have new behaviors of X * T (β) due to the addition of the third patch. We show that we can have the following situations, which do not exist in the two-patch model: 6. Asymmetric dispersal. Arditi et al. [START_REF] Arditi | Asymmetric dispersal in the multi-patch logistic equation[END_REF] generalized the mathematical analysis of the two-patch model [START_REF] Arditi | The perfect mixing paradox and the logistic equation: Verhulst vs[END_REF] to the asymmetric case

• We can have simultaneously dX * T dβ (0) < 0 and X * T (+∞) > K 1 + K 2 + K 3 ,
       dx 1 dt = r 1 x 1 1 - x 1 K 1 + β(γ 12 x 2 -γ 21 x 1 ), dx 2 dt = r 2 x 2 1 - x 2 K 2 + β(γ 21 x 1 -12 x 2 ),
where γ 12 = γ 21 . They showed that there are only three cases as in the symmetric dispersal case. In this section we give some results on the general multi-patch logistic growth, coupled with asymmetric migration terms :

dx i dt = r i x i 1 - x i K i + β n j=1,j =i (γ ij x j -γ ji x i ) , i = 1, • • • , n. (39)
where γ ij ≥ 0 denotes the incoming flux from patch j to patch i, for i = j. As in the symmetric dispersal case (4), the system (39) can be written in the form [START_REF] Cvetković | Stabilizing the Metzler matrices with applications to dynamical systems[END_REF] where

γ ii = - n j=1,j =i γ ji , i = 1, • • • , n (40)
denotes the outgoing flux of patch i. The connectivity matrix Γ 0 is defined by [START_REF] Deangelis | Effects of diffusion on total biomass in heterogeneous continuous and discrete-patch systems[END_REF] and the matrix Γ is defined by [START_REF] Deangelis | Effects of dispersal in a non-uniform environment on population dynamics and competition: a patch model approach[END_REF]. Its columns sum to 0 since the diagonal elements γ ii , defined by (40) in such a way that what comes out of a patch is distributed between the other n -1 patches. As in the symmetric dispersal case, the positive cone R n + is positively invariant for the system (39). We have the following result, whose proof needs the Lemma B.2, which is recalled in Appendix B. Theorem 6.1. Assume that the matrix Γ := (γ ij ) n×n (or equivalently, the connectivity matrix Γ 0 ) is irreducible. The model (39) has a unique positive equilibrium point which is GAS in the positive cone R n + \ {0}. Proof. Consider the system (39). We define the following matrix

A = βΓ + diag(r 1 , • • • , r n ),
where Γ is the matrix defined by [START_REF] Deangelis | Effects of dispersal in a non-uniform environment on population dynamics and competition: a patch model approach[END_REF]. Therefore A = (a ij ) where

a ij = βγ ij if i = j r i -β k =i γ ik if i = j (41)
Note that, the matrix (41) is the Jacobian matrix of the system (39) evaluated at x = 0. According to a result of Lu and Takeuchi [16, Corollary 1], the system (39) possesses a globally stable positive equilibrium if s(A) > 0 where s(A) is the stability modulus of the matrix A. Let us prove that s

(A) > 0. Let u = (1, • • • , 1) T .
We have

A T u = (r 1 , • • • , r n ) T ≥ λu, where λ = min{r 1 , • • • , r n } > 0.
Therefore, since A is a Metzler matrix, according to Lemma B.2, we have s

(A) = s(A T ) ≥ λ > 0.
Some of the results obtained in the symmetrical case can be extended to the asymmetric case. For instance, we have the following result, which extend Proposition 3.1 to the asymmetric case and [3, Proposition 1], obtained in the two-patch model: Proposition 6.2. If r 1 = • • • = r n , then the total equilibrium population, defined by ( 12) satisfies X * T (β) ≤ n i=1 K i , for all β ∈ [0, ∞[. Proof. The equilibrium point E * (β) is the solution of the algebraic system:

r i x i 1 - x i K i + β n j=1,j =i (γ x j -γ ji x i ) = 0, i = 1, • • • , n.
As in the symmetric dispersal case, the sum of these equations shows that E * (β) satisfies the equation [START_REF] Arditi | Asymmetric dispersal in the multi-patch logistic equation[END_REF]. Therefore E * (β) belongs to the ellipsoid defined by [START_REF] Gantmacher | The Theory of Matrices[END_REF].

The rest of the proof is the same as the proof of Prop. 3.1 and simply consists in noticing that the interior of the ellipsoid is strictly convex.

Conclusion.

The goal of this paper was to generalize to a multi-patch model the results obtained in [START_REF] Arditi | In dispersal always beneficial to carrying capacity? New insights from the multi-patch logistic equation[END_REF] for a two-patch model. The migration between patches is modeled by a symmetric Metzler matrix, called the connectivity matrix. When the connectivity matrix is irreducible, the system is shown (Prop. 2.2) to have a unique non-trivial equilibrium, which furthermore is globally asymptotically stable. One of the questions we have looked at is Question 7.1. Is it possible, depending on the migration rate, that the total equilibrium population be larger than the sum of the capacities of each patch ?

This question is of ecological importance since the answer gives the conditions under which dispersal is either beneficial or detrimental to total equilibrium population. We have answered this question in several particular cases. In the case that the patches do not differ with respect to the carrying capacity (i.e. K 1 = • • • = K n ), migration has no effect on the total equilibrium population (see Proposition 3.2). In the case when the patches do not differ with respect to the intrinsic growth rate (i.e., r 1 = • • • = r n ), the effect of migration is always detrimental (see Propositions 3.1 and 6.2). In the case that the patches do not differ with respect to the the parameter α = r/K quantifying intraspecific competition (i.e., α 1 = • • • = α n ), the effect of migration is always beneficial (see Proposition 3.3). We have also answered Question 7.1 in the particular case when all patches, but one, are identical in carrying capacity, growth rate, and migration rate to and from the "central" patch. The results are similar to those of [START_REF] Arditi | In dispersal always beneficial to carrying capacity? New insights from the multi-patch logistic equation[END_REF]: either patchiness has always a beneficial effect on the total equilibrium population, or this effect is always detrimental, or there exists a critical value β 0 of the migration rate β such that the effect is beneficial for lower values of β and detrimental for the higher values (see Proposition 3.4).

In Section 4 we looked at another particular case, that of perfect mixing, when the migration rate goes to infinity, in other words, when there is no restriction whatsoever on travel. We computed the equilibrium in this situation, and by perturbation arguments (see [START_REF] Tikhonov | Systems of differential equations containing small parameters in the derivatives[END_REF]), we proved that the dynamics in this ideal case provide a good approximation to the case when the migration rate is large.

In Section 5 we considered the case when there are only three patches. We provided numerical evidence to the fact that the answer to Question 7.1 is more subtle than in the two-patch case, in particular there may be at least two critical values where the answer changes.

In Section 6 we considered the asymmetrical dispersal case, that is, the connectivity matrix is no longer assumed to be symmetric. We proved that the globally asymptotically stable equilibrium still exists in this case (see Theorem 6.1). For the convenience of the reader, we give here the details of the proof of (47). Adding the equations of (46) we get r(n -1)x 1 -x K + r n x n 1 -

x n K n = 0,
which is the equation of an ellipse E passing through the points (0, 0), (K, 0), (0, K n ) and A = (K, K n ). The first equation of (46) represents a parabola P of equation x n = P (x), where

P (x) = x 1 - r βγ + r βγK x .
So, the solutions of system (46) are defined by the intersection points between Ellipse E and Parabola P, as it is shown in Fig. 5. Since P (K) = 1 + r βγ > 1, this intersection lies in the triangle ABC, where A = (K, K n ), B = (K, K) and C = (K n , K n ), which proves (47).

γ ij (K j -K i ) . (48)

  2 and the same procedures to study the case where K > K n .

Figure 3 .

 3 Figure 3. Total equilibrium population X * T of the system (4) (n = 3) as a function of migration rate β. The figure on the right is a zoom, near the origin, of the figure on the left. The parameter values are given in Table 5.2.

Figure r 1 r 2 r 3 K 1 K 2

 312 Figure r 1 r 2 r 3 K 1 K 2 K 3 dX * T (0) dβ X * T (+∞) 2 0.12 18 0.02 0.5 1.5 2 -79.19 4.44 > K i = 4 3 0.04 3 0.2 0.5 6 9.5 299.33 16.17 > K i = 16 4 4 0.7 0.06 5 1 4 -24.58 9.42 < K i = 10

  as shown in Fig. 2. • We can have dX * T dβ (0) > 0 and X * T (+∞) > K 1 + K 2 + K 3 and there exist values of β for which we have X * T (β) < K 1 + K 2 + K 3 , as shown in Fig. 3. • We can have dX * T dβ (0) < 0 and X * T (+∞) < K 1 + K 2 + K 3 and there exist values of β for which we have X * T (β) > K 1 + K 2 + K 3 , as shown in Fig. 4. Therefore the equality X * T (β) = K 1 + K 2 + K 3 can occur for two positive values of β, not only for a unique positive value as in the two-patch case.

Figure 5 .

 5 Figure 5. The intersection point (x * , x * n ), between Ellipse E and Parabola P, lies in the interior of triangle ABC. (a): the case K < K n . (b): the case K > K n .

Lemma A. 3 .

 3 The derivative of the total equilibrium population X * T (β) at β = 0, is given by dX *

Acknowledgments. The authors thank the Editor and the two anonymous reviewers for their constructive comments which have greatly improved this work.

Appendix A. Some properties of the total equilibrium population. In this section, we give some properties of the total equilibrium population X * T (β) = n i=1 x * i (β).

Lemma A.1. The total equilibrium population X * T satisfies the following relation:

Proof. The equilibrium point E * (β) satisfies the system

Dividing (43) by α i x * i , one obtains

.

Taking the sum of these expressions gives (42).

Proposition A.2. Consider the system (4). If (16) hold then x * 1 = . . . = x * n-1 =: x * and

Proof. Suppose that ( 16) hold, that is to say

Then the equations ( 10) are written

Now, consider the following system of algebraic equations

obtained from (44) by replacing x i = x, for i = 1 • • • n -1 and using the conditions

Proof. By differentiating the equation (42), at β = 0, we get:

, which gives (48), since

Stability modulus of a Metzler matrix. For the ease of the reader, we recall in this section the proof of a result which gives a minoration of the stability modulus of a Metzler matrix, see Lemma B.2. This Lemma is a corollary of result which gives a minoration of the spectral radius of a non negative matrix, see Lemma B.1. Recall that the spectral radius of matrix A is defined as

We have the following result [17, Lemma 8]:

Proof. If Au ≥ λu then, since A is non negative, A k u ≥ λ k u for all k. Therefore A k ≥ λ k for any matricial norm. Using the Gelfand formula ρ(A) = lim k→∞ A k 1 k , we obtain that ρ(A) ≥ λ.

Recall that the stability modulus of a matrix A is given by s(A) = max {Re(λ) : λ is an eigenvalue of A} .

A matrix A = (a ij ) is said to be Metzler if a ij ≥ 0 for i = j. We have the following result [5, Lemma 8]:

Lemma B.2. Let A be a Metzler matrix. Let u ∈ R n and λ ∈ R. If Au ≥ λu then s(A) ≥ λ.

Proof. Let A be a Metzler matrix, there exists h > 0 such that A + hI, where I is the identity matrix, is non negative. Let u and λ be such that Au ≥ λu. Since (A + hI)u ≥ (λ + h)u, using Lemma B.1, we deduce that ρ(A + hI) ≥ λ + h. According to the Perron-Frobenius Theorem [11, Theorem 3, page 66], we have s(A + hI) = ρ(A + hI). Therefore we have s(A + hI) ≥ λ + h. Using s(A + hI) = s(A) + h, we obtain s(A) ≥ λ.


Haopeng Zhang
email: haopeng.zhang@louisville.edu

Qing Hui
email: qing.hui@unl.edu

Emmanuel Moulay
email: emmanuel.moulay@univ-poitiers.fr

Patrick Coirault

Sign gradient descent method based bat searching algorithm with application to the economic load dispatch problem

Inspired by the echolocation behaviors of the bats and the swarm intelligence optimization, bat searching algorithm (BA) was developed to solve unconstrained optimization problems efficiently. However, due to the lack of the gradient term, the accuracy of the BA is not superior, and the enhancement of the algorithm is still of vital importance. The sign gradient descent method (SGD) is a first-order optimization method involving only the sign of the gradient of the function to minimize. Most importantly, the convergence and optimality issues of the SGD have been rigorously studied, which guarantees the competitive performance of SGD method. Therefore, in this paper, a combination of the BA and SGD method is proposed by integrating the SGD term into the update equation of the bats during the searching process. With the social behavior among the bats and the sign gradient descent method, the proposed algorithm shows significant improvement comparing with the original algorithm. Moreover, the convergence issue of the proposed algorithm is studied from system dynamics perspective. The numerical evaluations are provided to demonstrate the improvement of the proposed sign gradient descent method based bat searching algorithm. In the end, the economic load dispatch problem for the power system is studied as an application of the proposed BA algorithms. Based on the numerical results, the proposed BA shows superior performance.

I. INTRODUCTION

Inspired by the echolocation behavior of microbats, Bat Algorithm (BA) is a metaheuristic for global optimization which was first proposed in [START_REF] Yang | A new metaheuristic BAT-inspired algorithm[END_REF] and widely applied to solve challenging optimization problems in various areas such as robotics [START_REF] Suárez | Make robots be bats: specializing robotic swarms to the bat algorithm[END_REF], modern communication in [START_REF] Das | An efficient side lobe reduction technique considering mutual coupling effect in linear array antenna using bat algorithm[END_REF], or energy systems [START_REF] Bahmani-Firouzi | Optimal sizing of battery energy storage for micro-grid operation management using a new improved bat algorithm[END_REF], [START_REF] Ali | Optimization of power system stabilizers using bat search algorithm[END_REF], [START_REF] Adarsh | Economic dispatch using chaotic bat algorithm[END_REF]. Extensive work has been conducted to improve the performance of the original BA. In [START_REF] Meng | A novel bat algorithm with habitat selection and doppler effect in echoes for optimization[END_REF], a novel Bat searching algorithm (NBA) was proposed by mimicking bats habitat selection and their self-adaptive compensation for Doppler effect in echoes into the basic BA to improve the accuracy of the original BA. A bat algorithm based on iterative local search and stochastic inertial weight (ILSSI-WBA) was proposed in [START_REF] Gan | A new bat algorithm based on iterative local search and stochastic inertia weight[END_REF] to improve the algorithm to jump out of the local optima and to enhance the diversity and flexibility of bat population. In order to balance exploration patrick.coirault@univ-poitiers.fr behavior which expands the search area and exploitation behavior which favors a local minimum search, a new algorithm, known as directional Bat Algorithm (dBA) [START_REF] Chakri | New directional bat algorithm for continuous optimization problems[END_REF] or Enhanced Bat Algorithm (EBA) has been [START_REF] Yılmaz | A new modification approach on bat algorithm for solving optimization problems[END_REF] proposed. Cooperative bat searching algorithm was proposed in [START_REF] Zhang | Cooperative bat searching algorithm: A combined perspective from multiagent coordination and swarm intelligence[END_REF] by combining the swarm intelligence and cooperative control to improve the original BA while a couple spring forced BA [START_REF]A coupled spring forced bat searching algorithm: Design, analysis and evaluation[END_REF] was developed by guiding the bats with coupled springs' oscillation behavior and echolocation. Another main research direction is to extend the basic BA to address different types of optimization problems, such as, discrete optimization in [START_REF] Mirjalili | Binary bat algorithm[END_REF], [START_REF] Zhang | A binary cooperative bat algorithm based optimal topology design of leader-follower consensus[END_REF], multi-objective and many-objective optimization in [START_REF] Tuba | Hybridized bat algorithm for multiobjective radio frequency identification (rfid) network planning[END_REF], [START_REF] Zhang | Many objective cooperative bat searching algorithm[END_REF].

However, due to the lack of the gradient term, the accuracy of the BA is not superior, and the enhancement of the algorithm is still of vital importance. A first attempt to couple BA and SGD was proposed in [START_REF] Reddy | An automated gradient enhanced bat algorithm[END_REF] but both algorithms are totally decoupled. The sign gradient descent method (SGD) is a first-order optimization method involving only the gradient of the function to minimize. It has been first proposed in the Rprop algorithm for machine learning in [START_REF] Riedmiller | A direct adaptive method for faster backpropagation learning: the RPROP algorithm[END_REF] and its mathematical properties have been studied in [START_REF] Moulay | Properties of the sign gradient descent algorithms[END_REF]. One of the main interest of this method lies in the fact that the step size can be tuned precisely independently of the amplitude of the gradient avoiding the vanishing gradient problem explained in [START_REF] Du | Neural networks and statistical learning[END_REF]. Therefore, in this paper, a combination of the BA and SGD methods is proposed by adding the SGD term into the update equations of the bats during the searching process. With the social behavior among the bats and the sign gradient descent method, the proposed algorithm shows significant improvement compared to the original algorithm when solving the benchmark functions. Moreover, the convergence issue of the proposed algorithm is studied from system dynamics perspective. The convergence condition of the proposed algorithm is provided. The Economic Load Dispatch (ELD) problem for power system is considered and numerically solved in this paper as an application. With the numerical comparison with other BA variants, the proposed algorithm shows superior performance.

The article is organized as follows. After some preliminaries given in Section II, the SGD-BA algorithm is developed in Section III. The convergence results are provided in Section IV. Numerical evaluations are provided in Section V to illustrate the efficiency of the proposed algorithms. An application problem is studied in Section VI. Finally, a conclusion is addressed in Section VII.

II. PRELIMINARIES

Denoting y = (y 1 , . . . , y n ) T ∈ R n , the gradient of a differentiable function F : R n → R is the vector where sign refers to the sign function defined by

∇F (y) = ! ∂F (
sign(y) := # $ % $ & -1 if y < 0, 0 if y = 0, 1 if y > 0.
rand denotes a random number ranging from 0 to 1. Let x = (x 1 , . . . , x n ) ∈ R n and y = (y 1 , . . . , y n ) ∈ R n , we denote x × y = (x 1 y 1 , . . . , x n y n )

A. BA algorithm

The BA was based on the echolocation or bio-sonar characteristics of microbats, see [START_REF] Yang | A new metaheuristic BAT-inspired algorithm[END_REF].

Let x i (k) ∈ R n , v i (k) ∈ R n , f i (k) ∈ R
be the position, velocity and frequency information of bat i at iteration k, respectively. We suppose that f i belongs to set [f min , f max ]. The update rules are

f i (k) = f min + (f max -f min )β v i (k) = v i (k -1) + (p(k -1) -x i (k -1)) f i (k) (1) x i (k) = x i (k -1) + v i (k)
where k ∈ N, 1 ≤ i ≤ q where q ∈ N is the number of bats, β ∈ [0, 1] is a random vector drawn from a uniform distribution and p(k -1) is the best global solution of all bats at step k -1.

The local search is also embedded in the BA algorithm by using random walk around the current best solutions for each bat as follows

x new = x old + #A(k) (2)
where # ∈ [-1, 1] and A(k) = 〈A i (k)〉 ∈ R n is the average loudness of all the bats i at time step k. The loudness A i and pulse rate r i for bat i are governed by the following two equations

A i (k + 1) = α A i (k) r i (k + 1) = r i (k)(1 -exp (-γk)) (3)
where α and γ are constants.

Let us now recall in Algorithm 1 the BA algorithm for a cost function f : R n → R which must be minimized. And x, x, v, v are the lower bounds and upper bounds for positions and velocities, respectively.

Algorithm 1 Bat Searching Algorithm Initialize the bat's position and velocity with two uniformly distributed random vectors

x i ∼ U (x, x), v i ∼ U (v, v), respectively.
Initialize the global best position p, loudness A i (0), pulse rate r i (0), and frequency f i (0). while k < MaxIter; do while i <= q do Update the bat's frequency:

f i (k) = f min + (f max -f min )β;
Update the bat's velocity:

v i ← v i + (x i -p)f i ;
Update the bat's position:

x i (k) = x i (k -1) + v i (k); if rand > r i then

III. SGD-BA ALGORITHM

In this section, the gradient descent method is added into the update equation ( 1) of the original BA to obtain SGD-BA to improve the performance of the original bat searching algorithm. More precisely, the local search given by Equations ( 2) and ( 3) is replaced by SGD local search directly in Equation [START_REF] Yang | A new metaheuristic BAT-inspired algorithm[END_REF]. Mathematically, it reads as

f i (k) = f min + (f max -f min )β v i (k) = ωv i (k -1) + (p(k -1) -x i (k -1)) f i (k) -SGD(k -1) x i (k) = x i (k -1) + v i (k)
where ω is the weight inertia,

SGD(k -1) = d(k -1) × sign(∇(F (x i (k-1)))) and d(k) ∈ R n is the step size of the searching. Overall, the SGD-BA is described in Algorithm 2.
With the combination of the BA and SGD method, the proposed SGD-BA will search and find the optimal solution with high efficiency.

In the following, we propose two approximations of the partial derivative term.

Modification 1: for 1 ≤ i ≤ q and 1 ≤ j ≤ n:

∂F (x i (k)) ∂x ij (k) ≈ F (xi1(k),...,xij (k+1),...,xin(k))-F (xi(k)) xij (k+1)-xij (k) (4)
Modification 2: for 1 ≤ i ≤ q and 1 ≤ j ≤ n:

∂F (x i (k)) ∂x ij (k) ≈ F (X(k + 1)) -F (X(k)) x ij (k + 1) -x ij (k) (5)
with

x i (k) = (x i1 (k), . . . , x in (k)) ∈ R n .
Moreover, for the step size, we choose

d(k) = p(k -1) -x i (k -1) N of B (6)
where N of B represents the number of bats used in the algorithm. We have that lim k→∞ d(k) = 0 which is a necessary condition for the convergence of the algorithm as explained in [START_REF] Moulay | Properties of the sign gradient descent algorithms[END_REF].

Remark 3.1: Overall, Modification 1 is a better approximation of the partial derivative than Modification 1. However, SGD-BA takes more computation time when using Modification 2 than when using Modification 1 and it can be useful in practice to use Modification 1.

Algorithm 2 Sign-gradient descent based bat searching algorithm Initialize the bat's position and velocity with two uniformly distributed random vectors

x i ∼ (lbx, ubx), v i ∼ U (v, v), respectively.
Initialize the global best position p, loudness A i (0), pulse rate r i (0), and frequency f i (0). while k < MaxIter; do while i <= q do Update the bat's frequency:

f i (k) = f min + (f max -f min )β;
Update the bat's velocity:

v i ← ωv i + (x i -p)f i -SGD;
Update the bat's position:

x i (k) = x i (k -1) + v i (k); if rand >

IV. CONVERGENCE STUDY

In this section, the convergence issue of the proposed SGD-BA with deterministic parameters is studied.

Theorem 4.1: All the bats in the SGD-BA converge to the global best solution if and only if -1 < ω < 1, f > 0 and 2ω + 2 -f > 0.

Proof: For bat i, we have

f i (k) = f min + (f max -f min )β v i (k) = ωv i (k -1) + (p(k -1) -x i (k -1)) f i (k -1) -SGD(k -1) x i (k) = x i (k -1) + v i (k)
and then,

x i (k) = x i (k -1) + ωv i (k -1) + (p(k -1) -x i (k -1)) f i (k) -SGD(k -1)
by replacing v i (k -1) = x i (k -1) -x i (k -2), we have:

x i (k) = (1 + ω -f )x i (k -1) -ωx i (k -2) +p(k -1)f -SGD(k -1)
The Z-transformation is applied to the above equation:

µ i (z) = (1 + ω -f )zµ i (z) -ωz 2 µ i (z) + Z{p(k -1)f } -Z{SGD(k -1)}
where µ i (z) is the Z-transform of x i (k). Therefore,

[1 -(1 + ω -f )z -1 + ωz -2 ]µ i (z) = Z{p(k -1)f } -Z{SGD(k -1)} and µ i (z) = Z{p(k -1)f } [1 -(1 + ω -f )z -1 + ωz -2 ] - Z{SGD(k -1)} [1 -(1 + ω -f )z -1 + ωz -2 ] µ i (z) = 1 z -1 -1 f P * [1 -(1 + ω -f )z -1 + ωz -2 ] - Z{SGD(k -1)} [1 -(1 + ω -f )z -1 + ωz -2 ]
Based on the final value theorem:

lim k→∞ x i (k) = lim z→1 µ i (z)(z -1) = lim z→1 zf P * [1 -(1 + ω -f )z -1 + ωz -2 ] -lim z→1 (z -1) Z{SGD(k -1)} [1 -(1 + ω -f )z -1 + ωz -2 ]
To grantee the stability of the system, the roots of

[1 -(1 + ω -f )z -1 + ωz -2
] should be inside the unit circle, and we have

• -1 < ω < 1 • f > 0 • 2ω + 2 -f > 0 and lim k→∞ x i (k) = p.
A numerical example is provided here. Assume three bats are considered in the SGD-BA algorithm to solve the minimization problem of objective function f (x) = x 2 . Let ω = 0.8, and f i = 2, then the positions and velocities of the bats during the searching process are shown in Fig. 1 and Fig. 2, which show the convergence of the algorithm. Moreover, the best solutions found in each iteration are shown in Fig. 3, and the fitness values during the searching process are provided in Fig. 4.

V. NUMERICAL EVALUATIONS

In this section, the proposed SGD-BAs are numerically evaluated.

A. Benchmark functions

The standard benchmark functions are chosen from the 2018 IEEE Congress on Evolutionary Computation (CEC). Due to the space limitation, all the functions are not fully reviewed, for more information, one can check [START_REF] Li | Benchmark functions for the CEC'2013 special session and competition on large-scale global optimization[END_REF].

B. Numerical results

All the algorithms run 20 times with random initial positions with 100 bats. The max, min, mean and the standard deviation values are calculated and are compared, which are given in Tables I-VI. In the last line of the tables, number of the best solutions are listed. Overall, equation ( 4) being a better approximation of the partial derivative than Equation [START_REF] Ali | Optimization of power system stabilizers using bat search algorithm[END_REF], it explains why we obtain better results for SGD-BA with Equation (4) than with Equation ( 5) in Table IV. However, SGD-BA takes more computation time when using Equation (4) than when using Equation ( 5) and it can be useful in practice to use Equation [START_REF] Ali | Optimization of power system stabilizers using bat search algorithm[END_REF].

VI. APPLICATION TO ECONOMIC LOAD DISPATCH PROBLEM

To demonstrate the performance of the SGD method and the proposed BA algorithms, the Economic Load Dispatch (ELD) problem for power systems is considered and numerically solved in this section as an application. The ELD problem considers a series of power generation and power transmission facilities, each with associated operating limits, costs, and losses. Given a power demand for the entire grid, the optimal outputs for each generating station must be determined in order to minimize overall cost, see [START_REF] Park | A particle swarm optimization for economic dispatch with nonsmooth cost functions[END_REF] for details. Due to the nonlinear cost functions and the presence of both inequality and equality constraints, the ELD problem is a rigorous test that offers insight into how different configurations affect the performance of an algorithm. Mathematically, the ELD problem can be expressed as:

min : n ' i=1 F i (P i ) (7)
subject to :

P min i ≤ P i ≤ P max i , (8) n '
i=1

P i = D + P l , (9)
P l = n ' i=1 n ' i=j B i,j P i P j , (10)
F i,j (P (i)) = a i P 2 i + b i P i + c i (11)
where the system consists of n generating stations each with a load of P i . The fuel cost for each plant is given by the quadratic function F i (P i ). The load for each generating station is limited by the lower and upper bounds P min i and P max i , respectively, while the total load is P l . The transmis- sion losses are described by the B matrix. In this section, the BA, SGD-BA(4) and SGD-BA(5) are applied to solve the ELD problem ( 7)- [START_REF] Zhang | Cooperative bat searching algorithm: A combined perspective from multiagent coordination and swarm intelligence[END_REF]. 30 bats are considered in all the BA algorithms with 500 iterations. 20 executions have been conducted, and the max, min, standard derivation (SD), mean, and median values are compared and provided in Table . V. Based on the simulation results, the proposed two sign-gradient based BA significantly improve the original BA. And the second modification showed the best results.

VII. CONCLUSION

In this paper, two sign-gradient descent modifications are applied to the bat searching algorithm in order to enhance the accuracy of the BA. Due to the gradient term in the proposed algorithms, the SGD-BAs show competitive improvement compared with the original BA by solving the 2018 IEEE Congress on Evolutionary Computation's standard benchmark functions. Moreover, the convergence analysis of the proposed algorithm has been conducted, and the require- ments of the parameters are derived for convergence. In the end, the ELD problem is considered and numerically solved by the three bat algorithms. Based on the numerical results, the proposed sign-gradient BAs show better performance than the original BA and two other BA variations in the literature.

3 Fig. 1 . 3 Fig. 2 .Fig. 3 .Fig. 4 .

 313234 Fig. 1. Convergence of the SGD-BA: positions

  r i then Select a solution among the best solutions else Generate a local position around the best solution. end if rand < A i and f (x i ) < f (p) then Accept the new solutions Increase r i and reduce A i

	else
	Keep the same p
	end
	end
	end

TABLE I MAXIMUM

 I COMPARISONS BETWEEN SGD-BA(4) USES EQUATION (4)

	AND SGD-BA(5) USES EQUATION (5)
	Functions	BA	SGD-BA(4)	SGD-BA(5)
	f 1	2.2067e+11	2.2055e+11	2.2042e+11
	f 2	4.9380e+04	5.1108e+04	5.1018e+04
	f 3	2.1232e+01	2.1005e+01	2.1208e+01
	f 4	2.3045e+13	2.0412e+13	2.0674e+13
	f 5	4.6616e+07	4.0791e+07	5.0117e+07
	f 6	1.0490e+06	1.0421e+06	1.0437e+06
	f 7	4.4566e+14	2.5223e+14	2.7483e+14
	f 8	7.0595e+17	8.6122e+17	9.1962e+17
	f 9	3.9952e+09	3.6655e+09	3.8955e+09
	f 10	9.5325e+07	9.5430e+07	9.5730e+07
	f 11	1.8293e+16	1.7352e+16	1.6622e+16
	f 12	2.3174e+12	2.9851e+12	2.4451e+12
	f 13	6.5514e+12	5.2442e+12	5.7714e+12
	f 14	3.1320e+16	3.7821e+16	3.6381e+16
	f 15	3.8474e+14	1.0233e+14	1.1123e+14
	# of the best	5 out of 15	8 out of 15	2 out of 15

TABLE II MINIMUM

 II COMPARISONS BETWEEN SGD-BA(4) USES EQUATION (4)

	AND SGD-BA(5) USES EQUATION (5)
	Functions	BA	SGD-BA(4)	SGD-BA(5)
	f 1	1.8409e+11	1.8328e+11	1.8248e+11
	f 2	4.3991e+04	4.4003e+04	4.3773e+04
	f 3	2.1129e+01	2.1003e+01	2.1113e+01
	f 4	6.4297e+12	3.0125e+12	3.2055e+12
	f 5	3.4032e+07	3.9709e+07	3.7870e+07
	f 6	1.0211e+06	1.01071e+06	1.0203e+06
	f 7	3.6324e+13	3.8211e+13	4.1261e+13
	f 8	2.0124e+17	1.9841e+17	1.6494e+17
	f 9	2.6673e+09	2.6099e+09	2.6569e+09
	f 10	9.0103e+07	8.9881e+07	8.7831e+07
	f 11	2.6236e+15	2.5429e+15	2.7939e+15
	f 12	1.9961e+12	2.0186e+12	1.9886e+12
	f 13	3.1547e+11	3.0626e+11	3.2266e+11
	f 14	2.9560e+15	2.3127e+15	1.8397e+15
	f 15	2.6404e+11	2.9621e+11	3.2036e+11
	# of the best	3 out of 15	7 out of 15	7 out of 15

TABLE III STANDARD

 III DEVIATION COMPARISONS BETWEEN SGD-BA(4) USES EQUATION (4) AND SGD-BA(5) USES EQUATION[START_REF] Ali | Optimization of power system stabilizers using bat search algorithm[END_REF]

	Functions	BA	SGD-BA(4)	SGD-BA(5)
	f 1	9.1698e+09	1.2123e+10	1.0144e+10
	f 2	1.2191e+03	1.6361e+03	1.7831e+03
	f 3	2.7249e-02	2.9941e-02	2.7894e-02
	f 4	4.2594e+12	4.1901e+12	4.4449e+12
	f 5	2.7758e+06	2.6825e+06	2.7652e+06
	f 6	6.8454e+03	6.5851e+03	5.8550e+03
	f 7	1.0716e+14	6.7881e+13	6.1388e+13
	f 8	1.3153e+17	1.6193e+17	1.3103e+17
	f 9	3.1330e+08	2.8234e+08	2.7825e+08
	f 10	1.4175e+06	1.2972e+06	1.7479e+06
	f 11	3.9061e+15	3.5565e+15	3.5565e+15
	f 12	8.0107e+10	9.8551e+10	1.2875e+11
	f 13	1.6392e+12	1.1441e+12	1.3254e+12
	f 14	8.3149e+15	7.6873e+15	7.9783e+15
	f 15	8.5946e+13	1.9634e+13	2.4633e+13
	# of the best	4 out of 15	6 out of 15	5 out of 15

TABLE IV MEAN

 IV VALUE COMPARISONS BETWEEN SGD-BA(4) USES EQUATION (4)

	AND SGD-BA(5) USES EQUATION (5)
	Functions	BA	SGD-BA(4)	SGD-BA(5)
	f 1	2.0517e+11	2.0321e+11	2.0084e+11
	f 2	4.6593e+04	4.6823e+04	4.6483e+04
	f 3	2.1177e+01	2.0614e+01	2.1164e+01
	f 4	1.0324e+13	9.2241e+12	9.4264e+12
	f 5	3.9574e+07	3.6561e+07	3.8886e+07
	f 6	1.0368e+06	1.0011e+06	1.0363e+06
	f 7	1.9538e+14	9.8141e+13	1.1431e+14
	f 8	4.3026e+17	4.2181e+17	4.1478e+17
	f 9	3.4215e+09	3.3323e+09	3.2352e+09
	f 10	9.2536e+07	9.2750e+07	9.3008e+07
	f 11	6.8157e+15	7.0711e+15	6.1711e+15
	f 12	2.1201e+12	2.0853e+12	2.1583e+12
	f 13	1.8864e+12	1.7783e+12	1.4487e+12
	f 14	1.2638e+16	1.0051e+16	1.1155e+16
	f 15	3.6229e+13	1.1414e+13	1.2447e+13
	# of the best	1 out of 15	8 out of 15	6 out of 15

TABLE V COMPARISON

 V BETWEEN DIFFERENT BA ALGORITHMS FOR ELD USING 30 BATS

	BA	Max	Min	SD	Mean	Median
	BA	5.9738e+05	5.9465e+05	8.3622e+01	5.9684e+05	5.9725e+05
	SGD-BA(4) 5.2540e+05	4.7899e+05	1.5552e+04	4.9852e+05	4.9472e+05
	SGD-BA(5) 4.5768e+05	4.5695e+05	2.1379e+02	4.5727e+05	4.5729e+05

TABLE VI MEDIAN

 VI VALUE COMPARISONS BETWEEN SGD-BA(4) USES EQUATION (4) AND SGD-BA(5) USES EQUATION[START_REF] Ali | Optimization of power system stabilizers using bat search algorithm[END_REF]

	Functions	BA	SGD-BA(4)	SGD-BA(5)
	f 1	2.0723e+11	2.1201e+11	2.0070e+11
	f 2	4.6308e+04	4.6276e+04	4.6605e+04
	f 3	2.1180e+01	2.0987e+01	2.1165e+01
	f 4	8.4026e+12	8.9559e+12	8.9925e+12
	f 5	4.0037e+07	4.2661e+07	4.1156e+07
	f 6	1.1351e+06	1.0394e+06	1.0371e+06
	f 7	1.8825e+14	1.7797e+14	1.0791e+14
	f 8	4.1820e+17	3.7858e+17	3.9188e+17
	f 9	3.3901e+09	3.3451e+09	3.2382e+09
	f 10	9.2416e+07	9.1981e+07	9.2466e+07
	f 11	5.5110e+15	6.6475e+15	6.7615e+15
	f 12	2.1136e+12	2.4681e+12	2.1063e+12
	f 13	1.3413e+12	1.7013e+12	1.0074e+12
	f 14	9.8651e+15	9.9818e+16	1.0428e+16
	f 15	9.1951e+12	4.818569e+12	4.8185e+12
	# of the best	4 out of 15	4 out of 15	6 out of 15