# On the Theoretical Analysis of the Plant Propagation Algorithms

Authors: Muhammad Sulaiman; Abdellah Salhi; Asfandyar Khan; Shakoor Muhammad; Wali Khan
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

Plant Propagation Algorithms (PPA) are powerful and flexible solvers for optimisation problems. They are nature-inspired heuristics which can be applied to any optimisation/search problem. There is a growing body of research, mainly experimental, on PPA in the literature. Little, however, has been done on the theoretical front. Given the prominence this algorithm is gaining in terms of performance on benchmark problems as well as practical ones, some theoretical insight into its convergence is needed. The current paper is aimed at fulfilling this by providing a sketch for a global convergence analysis.


---

## Body

## 1. Introduction

The theoretical analysis of stochastic algorithms for global optimisation is not new and can be found in a number of sources such as [1–5]. The majority of the algorithms considered use random search one way or another to find the optimum solution [6–13]. Here, we consider the algorithmic scheme of the Plant Propagation Algorithm for continuous optimisation or PPA-C [14] and theoretically investigate its global convergence to the optimum solution. The optimisation problems of concern are continuous and defined in finite n-dimensional domains.The basic version of PPA [15] models the propagation of strawberry plants. The scheme uses short runners for exploitation or local search refinement while long runners are used for diversification and exploration of the search space. Since the propagation of strawberries is due to seeds as well as runners, a Seed-based Plant Propagation Algorithm (SbPPA) has also been introduced in [16]. Both PPA-C and SbPPA have been shown to be efficient on continuous unconstrained and constrained optimisation problems; statistical convergence analyses of PPA-C and SbPPA can be found in [9–11, 14–16].PPA-C [14, 15] consists of two steps:(1)
Initialization: a population of parent plants is generated randomly.(2)
Propagation: a new population is created from persistent parents (strawberry plants) and their children (new strawberry plants at the end of runners, i.e., a distance away from parent plants).LetS denote the search space such that S⊂Rn, where n is its dimension. By an iteration of PPA-C we mean a new generation of child plants produced by parent plants. These child plants are the result of either short or long runners [14, 16]. This is the basic setup that we consider to sketch a proof of convergence to the global optimum of a given continuous optimisation problem.The paper is organised as follows. Section2 presents the terminology used in the analysis of PPA-C. Section 3 analyses a population of plants. Section 3.1 describes the convergence analysis of PPA-C. Section 4 is the conclusion.


## 2. Terminology and Notation

We consider single objective minimization problems [17]. Xoptimal∈S such that f(Xoptimal)≤f(X) for all X∈S, where the objective function is defined as f:S⊂Rn→R, denotes the best spot for a plant in the search space. X is an n-dimensional position vector.The population at thegth iteration is denoted by popg={X1,g,X2,g,…,XNP,g}, where NP is the population size. The coordinates of runners, or more precisely their endpoints, are denoted by Xi=(xi1,xi2,…,xin)T, where n is the space dimension of the given problem.


### 2.1. Search Equations and Evaluation of New Plants

Variants of PPA can be found in [14–16]. In this paper we analyse PPA-C as Algorithm 1 of [14].In order to send a short or long runner,Xi′ is generated [14, 19–21], as in (1a), (1b), and (1c)(1a)xi,j′=xi,j+βjxi,jifrandj≤Pm,r≤NP,(1b)xi,j′=xi,j+xl,j-xk,jβjifrandj≤Pm,r≤NP(1c)xi,j′=xi,j+xi,j-xk,jβjifrandj≤PmorINj<4,r>NP,where NP is the population size, r is the Monte Carlo trial run counter, Pm is the modification probability, and randj∈(01) is a randomly generated number for each jth entry, j=1,2,…,n. The indices i,l,k=1,2,…,NP are mutually exclusive; that is, i≠k≠l. Another version of PPA called SbPPA [16] which is inspired by propagation via seeds implements the following search equation instead: (2)xi,j∗=xi,j+Lixi,j-θjifPR≤0.8,θj∈ajbji=1,2,…,NP;j=1,2,…,nxi,jOtherwise,where Li is a step drawn from the Lévy distribution [22] and θj is a random coordinate within the search space. Equations (1a), (1b), (1c), and (2) perturb the current solution, the results of which can be seen in Figures 1(a) and 1(b), respectively.Figure 1
Overall performance of (1a), (1b), (1c), and (2) on a design optimisation problem given in Appendix [16].
(a)
Perturbations by (1a), (1b), and (1c) (b)
Perturbations by (2)


### 2.2. A Case Study

Letc be the class of runners sent by the ith parent plant and stored in F. Each runner in class c is decomposed into two vectors τc and ωc, where τc denotes the vector of indices which are perturbed with respect to the current position of the plant, while ωc represents the vector of corresponding indices of the unperturbed coordinates with respect to the current position of plants. This can be represented as (3)τc∪ωc=1,2,…,n.To clarify this idea, let us take an example [17] of a newly generated runner by the ith plant as (4)Xi′=x1′,x2′,x3′,x4′,x5′T,such that (5)τc=1,2,4,ωc=3,5;then we can write (6)Xi′=Xτc′,Xωc′,where (7)Xτc′=x1′,x2′,0,x4′,0,Xωc′=0,0,x3′,0,x5′.The dot product of these vectors is zero, which shows that they are mutually orthogonal. Mathematically, this can be written as(8)Xτc′•Xωc′=0.LetVτc and Vωc denote two vector spaces such that (9)Vτc=containing  vectors  having  dimensions  as  inτc,Vωc=containing  vectors  having  dimensions  as  inωc.V τ c and Vωc are subspaces of Rn. This implies that Xτc′∈Vτc and Xωc′∈Vωc.A scalar objective functionf defined over Xi′ can be represented as (10)fXi′=fx1′,x2′,x3′,x4′,x5′,where x1′=Xτ1, x2′=Xτ2, x4′=Xτ3 and x3′=Xω1, x5′=Xω2, (11)⟹fXτ1,Xτ2,Xω1,Xτ3,Xω2=fXτc,Xωc,where (11) represents an objective value corresponding to a new runner in position Xi′. Similarly, different runners are produced to correspond to different classes c and evaluated by the same procedure. This procedure can be generalized for n-dimensional problems [1–3].


## 2.1. Search Equations and Evaluation of New Plants

Variants of PPA can be found in [14–16]. In this paper we analyse PPA-C as Algorithm 1 of [14].In order to send a short or long runner,Xi′ is generated [14, 19–21], as in (1a), (1b), and (1c)(1a)xi,j′=xi,j+βjxi,jifrandj≤Pm,r≤NP,(1b)xi,j′=xi,j+xl,j-xk,jβjifrandj≤Pm,r≤NP(1c)xi,j′=xi,j+xi,j-xk,jβjifrandj≤PmorINj<4,r>NP,where NP is the population size, r is the Monte Carlo trial run counter, Pm is the modification probability, and randj∈(01) is a randomly generated number for each jth entry, j=1,2,…,n. The indices i,l,k=1,2,…,NP are mutually exclusive; that is, i≠k≠l. Another version of PPA called SbPPA [16] which is inspired by propagation via seeds implements the following search equation instead: (2)xi,j∗=xi,j+Lixi,j-θjifPR≤0.8,θj∈ajbji=1,2,…,NP;j=1,2,…,nxi,jOtherwise,where Li is a step drawn from the Lévy distribution [22] and θj is a random coordinate within the search space. Equations (1a), (1b), (1c), and (2) perturb the current solution, the results of which can be seen in Figures 1(a) and 1(b), respectively.Figure 1
Overall performance of (1a), (1b), (1c), and (2) on a design optimisation problem given in Appendix [16].
(a)
Perturbations by (1a), (1b), and (1c) (b)
Perturbations by (2)


## 2.2. A Case Study

Letc be the class of runners sent by the ith parent plant and stored in F. Each runner in class c is decomposed into two vectors τc and ωc, where τc denotes the vector of indices which are perturbed with respect to the current position of the plant, while ωc represents the vector of corresponding indices of the unperturbed coordinates with respect to the current position of plants. This can be represented as (3)τc∪ωc=1,2,…,n.To clarify this idea, let us take an example [17] of a newly generated runner by the ith plant as (4)Xi′=x1′,x2′,x3′,x4′,x5′T,such that (5)τc=1,2,4,ωc=3,5;then we can write (6)Xi′=Xτc′,Xωc′,where (7)Xτc′=x1′,x2′,0,x4′,0,Xωc′=0,0,x3′,0,x5′.The dot product of these vectors is zero, which shows that they are mutually orthogonal. Mathematically, this can be written as(8)Xτc′•Xωc′=0.LetVτc and Vωc denote two vector spaces such that (9)Vτc=containing  vectors  having  dimensions  as  inτc,Vωc=containing  vectors  having  dimensions  as  inωc.V τ c and Vωc are subspaces of Rn. This implies that Xτc′∈Vτc and Xωc′∈Vωc.A scalar objective functionf defined over Xi′ can be represented as (10)fXi′=fx1′,x2′,x3′,x4′,x5′,where x1′=Xτ1, x2′=Xτ2, x4′=Xτ3 and x3′=Xω1, x5′=Xω2, (11)⟹fXτ1,Xτ2,Xω1,Xτ3,Xω2=fXτc,Xωc,where (11) represents an objective value corresponding to a new runner in position Xi′. Similarly, different runners are produced to correspond to different classes c and evaluated by the same procedure. This procedure can be generalized for n-dimensional problems [1–3].


## 3. Graphical and Theoretical Analysis of a Population of Plants

Algorithm1 states that the jth coordinate of an ith parent plant is perturbed with probability Pm and it remains unchanged with probability 1-Pm. Thus there are 2n possible runners to be generated for each ith parent plant using (1a), (1b), and (1c), where n is the space dimension of the given problem.Algorithm 1:PPA for constrained optimisation [14].
(1)  r← Counter for trial runs; NP← Population size (2) (3)  F← Population of runners (4) (5)  for  r=1 : 100 do (6) (7)if  r≤NP  then (8) (9)Create a random population of plants pop ={Xi∣i=1,2,…,NP}, and gather the best solution from each run. (10) (11)end if (12) (13)while  r>NP  do (14) (15)Use population popg formed by gathering all the best solutions of previous runs. Calculate INj value for each column j of popg. (16) (17)end while (18) (19)Evaluate the population pop. (20) (21)Assume number of runners to be nr=3, (22) (23)while   (the stopping criteria is not satisfied) do (24) (25)for  i=1 to NP  do (26) (27)for  h=1 to nrdo (28) (29)if  r≤NPthen (30) (31)if  rand ≤Pm  then (32) (33)Generate a new solution X′1 according to Equation (1a); (34) (35)Evaluate it and store it in F; (36) (37)end if (38) (39)if  rand ≤Pmthen (40) (41)Generate a new solution X′2 according to Equation (1b); (42) (43)Evaluate it and store it in F; (44) (45)end if (46) (47)else (48) (49)for  j=1:ndo (50) (51)if  (INj<4)  or  (rand ≤Pm)  then (52) (53)  Update the jth entry of Xi,i=1,2,…,NP, by using Equation (1c); (54) (55)end if (56) (57)Evaluate new solution X′3 and store it in F; (58) (59)end for (60) (61)end if (62) (63)end for (64) (65)end for (66) (67)Append F to current population; (68) (69)Sort the population in ascending order of objective values; (70) (71)Update current best; (72) (73)end while (74) (75)Return: Updated population and global best solution. (76) (77)  end for (78)Let, at any generationg, the random population be represented as pop={X1,g,X2,g,…,XNP,g}, where NP denotes the population size. To create a runner by using (1a), (1b), and (1c), for next generation g+1, PPA uses a population of parent plants at generation g for this purpose. It is not required to know about any other runner in generation g+1. This shows that all runners created at generation g+1 are statistically mutually independent. Furthermore, the initial population is random and all parent plants do not depend on each other. Thus, by induction, the runners at any further generations are mutually independent.From (1a), (1b), and (1c), a runner Xi′ may be formed by itself or by choosing three different coordinates from current population. In case of using (1a) [19], there are NP possibilities to send (long or short) runners as in Figure 1. On the other hand, by using (1b) or (1c), NP-1 vectors are used to send a new runner. Thus, in later cases there are P3NP-1 possibilities to send new runners. In (1a)–(1c), different possibilities are of the form (12)Xg′1=Xi,g+β·Xi,g,Xg′2=Xi,g+β·Xj,g,Xg′3=Xi,g+β·Xk,g,where Xi,g, Xj,g, and Xk,g are calculated according to (1a), (1b), and (1c), in which g denotes the current generation and β is an n-dimensional random vector within interval [-11]. The probability density functions (PDFs) [17, 23] of these new vectors Xi,g, β·Xi,g, β·Xj,g, and β·Xk,g can be written as pXi,gxi,g, pXi,gxi,g/β/β, pXj,gxj,g/β/β, pXk,gxk,g/β/β, respectively.The PDFs of new runners created with (1a)–(1c) are given in (15).Definition 1 (convolution operation⋆ [24]).
Letf(x) and g(x) be Laplace transformable piecewise continuous functions defined on 0∞. The convolution product of these two functions is again a function of x defined as (13)
f
⋆
g
x
=
∫
0
x
f
ξ
g
x
-
ξ
d
ξ
,
(14)
p
X
g
′
1
x
g
′
1
=
1
β
p
X
i
,
g
x
g
′
1
⋆
p
X
i
,
g
x
g
′
1
β
,
p
X
g
′
2
x
g
′
2
=
1
β
p
X
i
,
g
x
g
′
2
⋆
p
X
j
,
g
x
g
′
2
β
,
p
X
g
′
3
x
g
′
3
=
1
β
p
X
i
,
g
x
g
′
3
⋆
p
X
k
,
g
x
g
′
3
β
,
⇓
(15)
p
X
i
,
g
′
1
x
i
,
g
′
1
=
1
β
i
,
1
N
P
∑
∑
p
X
i
,
g
x
i
,
g
′
1
⋆
p
X
i
,
g
x
i
,
g
′
1
β
i
,
1
,
p
X
i
,
g
′
2
x
i
,
g
′
2
=
1
β
i
,
2
P
3
N
P
-
1
∑
∑
i
≠
j
p
X
i
,
g
x
i
,
g
′
2
⋆
p
X
j
,
g
x
i
,
g
′
2
β
i
,
2
,
p
X
i
,
g
′
3
x
i
,
g
′
3
=
1
β
i
,
3
P
3
N
P
-
1
∑
∑
i
≠
k
p
X
i
,
g
x
i
,
g
′
3
⋆
p
X
k
,
g
x
i
,
g
′
3
β
i
,
3
.LetXg+1 be any plant in generation g+1, Xg a plant position in current generation, and Xg′h the runner produced by the ith plant Xi,g at the generation g. Then the joint PDF of the parent plant in the next generation Xg+1 based on the parent Xg and runner Xg′h, where h=1,2,3, is given by (16)pXg+1,Xg,Xg′hxg+1,xg,xg′h=pXg+1∣Xg,Xg′hxg+1∣xg,xg′h·pXg′h∣Xgxg′h∣xg·pXgxg·pXg′hxg′h.Note that a runner is selected for the next population only if its rank is less than or equal toNP, the population size. Its objective value is less than the maximum objective value (in case of minimization problem) in the current population. Note also that instead of greedy selection we sort the population and eliminate those plants whose rank is higher than NP. The model for this selection mechanism can be represented as follows:(17)pXg+1∣Xg,Xg′hxg+1∣xg,xg′h=δxg+1-xg′h,∀Xg′h∈F;∀Xg∈pop;fXg≥fXg′h.


### 3.1. Convergence Analysis of PPA-C

For illustration purposes, we have implemented a combined version of (1a), (1b), (1c), and (2). This version of PPA-C called H-PPA-SbPPA is a hybridisation of PPA and SbPPA [10]. We have plotted the position of plants in populations through solving the Branin and Matyas test functions (see Figures 2, 3, 4, and 5). It is obvious from Figures 2 and 4 that (1a), (1b), and (1c) have generated short runners which exploit the search space locally. On the other hand, in Figures 3 and 5, (2) has generated a diverse range of solutions which are spread over the whole search space. This equation helps the algorithm escape from local minima and to explore the solution space better and hence the global search qualities of this algorithm. Mathematically this can be shown as follows.Figure 2
The exploitation capability of PPA while solving Branin test function.Figure 3
The exploration capability of SbPPA while solving Branin test function.Figure 4
A scatter plot of plants produced by PPA when optimising Matyas function [15].Figure 5
A scatter plot of plants produced by SbPPA when optimising Matyas function [16, 18].LetS denote the search space containing the solution of a given optimisation problem defined as (18)min⁡fX∣X∈S,where f(X) is the objective function. Then the optimal solution set [25] can be represented as (19)S∗=Xoptimal∣fXoptimal=minfX,∀X∈S,where Xoptimal is the optimum solution. The region of attraction [25] of the solution set S∗ is defined as (20)Sξ∗=X∣fXbest-fX<ξ,where ξ is a small positive real number and Xbest is the current best solution.In PPA-C, each parent plant producesXg′h runners (solutions), where h=1,2,3. The probability that at generation g a subset of the temporary population F′⊂F, containing solutions which are not good enough to be retained in the next generation by the selection model as in (17), is given as (21)pFg′∩Sξ∗=⌀≤1-ξg,where ξg is a small positive real number. Obviously, in all previous generations g-1, some of the solutions died and some succeeded to survive into the next generation. This shows that in previous generations we have some solutions Xg′h which do not belong to the region of attraction Sξ∗. (22)∏i=1g-1pFi′∩Sξ∗=⌀≤∏j=1g-11-ξj.At the end of each generation, the temporary populationF is appended to the main population popg. Then all individuals are sorted with respect to their objective values. The individuals with higher rank than size of population are omitted. Thus the probability that a generation g does not contain an optimum is given below (23)ppopg∩Sξ∗=⌀≈∏i=1g-1pFi′∩Sξ∗=⌀≤∏j=1g-11-ξj.After sorting the final population at the end of each generation, the probability that the optimum may exist in a subpopulationF′ (population of weak or dead runners/solutions) is less than that of the population popg. This can be represented as(24)limg→∞⁡ppopg∩Sξ∗≠⌀≥limg→∞⁡pFg′∩Sξ∗≠⌀=1-limg→∞⁡pFg′∩Sξ∗=⌀=1-∏g=1+∞1-ξg.Following [25–27], the right hand side term of inequality (24) is zero if the series ∑g=1+∞(ξg) diverges. The convergence of PPA follows and can be summarised in the theorem below.Theorem 2.
PPA converges to the global optimum with probability 1 if it is left to run for a reasonable amount of time [1–3]; in other words,(25)limg→∞⁡ppopg∩Sξ∗≠⌀=1.Remark 3.
Every populationpopg includes some solutions which are in set Sξ∗.Remark 4.
The above remark is due to the exploitation characteristic of PPA-C.Remark 5.
Remark3 implies that Xbest is always improving or changing its position until the optimum is reached.Remark 6.
X b e s t converges approximately to Xoptimal, as generation g grows.


## 3.1. Convergence Analysis of PPA-C

For illustration purposes, we have implemented a combined version of (1a), (1b), (1c), and (2). This version of PPA-C called H-PPA-SbPPA is a hybridisation of PPA and SbPPA [10]. We have plotted the position of plants in populations through solving the Branin and Matyas test functions (see Figures 2, 3, 4, and 5). It is obvious from Figures 2 and 4 that (1a), (1b), and (1c) have generated short runners which exploit the search space locally. On the other hand, in Figures 3 and 5, (2) has generated a diverse range of solutions which are spread over the whole search space. This equation helps the algorithm escape from local minima and to explore the solution space better and hence the global search qualities of this algorithm. Mathematically this can be shown as follows.Figure 2
The exploitation capability of PPA while solving Branin test function.Figure 3
The exploration capability of SbPPA while solving Branin test function.Figure 4
A scatter plot of plants produced by PPA when optimising Matyas function [15].Figure 5
A scatter plot of plants produced by SbPPA when optimising Matyas function [16, 18].LetS denote the search space containing the solution of a given optimisation problem defined as (18)min⁡fX∣X∈S,where f(X) is the objective function. Then the optimal solution set [25] can be represented as (19)S∗=Xoptimal∣fXoptimal=minfX,∀X∈S,where Xoptimal is the optimum solution. The region of attraction [25] of the solution set S∗ is defined as (20)Sξ∗=X∣fXbest-fX<ξ,where ξ is a small positive real number and Xbest is the current best solution.In PPA-C, each parent plant producesXg′h runners (solutions), where h=1,2,3. The probability that at generation g a subset of the temporary population F′⊂F, containing solutions which are not good enough to be retained in the next generation by the selection model as in (17), is given as (21)pFg′∩Sξ∗=⌀≤1-ξg,where ξg is a small positive real number. Obviously, in all previous generations g-1, some of the solutions died and some succeeded to survive into the next generation. This shows that in previous generations we have some solutions Xg′h which do not belong to the region of attraction Sξ∗. (22)∏i=1g-1pFi′∩Sξ∗=⌀≤∏j=1g-11-ξj.At the end of each generation, the temporary populationF is appended to the main population popg. Then all individuals are sorted with respect to their objective values. The individuals with higher rank than size of population are omitted. Thus the probability that a generation g does not contain an optimum is given below (23)ppopg∩Sξ∗=⌀≈∏i=1g-1pFi′∩Sξ∗=⌀≤∏j=1g-11-ξj.After sorting the final population at the end of each generation, the probability that the optimum may exist in a subpopulationF′ (population of weak or dead runners/solutions) is less than that of the population popg. This can be represented as(24)limg→∞⁡ppopg∩Sξ∗≠⌀≥limg→∞⁡pFg′∩Sξ∗≠⌀=1-limg→∞⁡pFg′∩Sξ∗=⌀=1-∏g=1+∞1-ξg.Following [25–27], the right hand side term of inequality (24) is zero if the series ∑g=1+∞(ξg) diverges. The convergence of PPA follows and can be summarised in the theorem below.Theorem 2.
PPA converges to the global optimum with probability 1 if it is left to run for a reasonable amount of time [1–3]; in other words,(25)limg→∞⁡ppopg∩Sξ∗≠⌀=1.Remark 3.
Every populationpopg includes some solutions which are in set Sξ∗.Remark 4.
The above remark is due to the exploitation characteristic of PPA-C.Remark 5.
Remark3 implies that Xbest is always improving or changing its position until the optimum is reached.Remark 6.
X b e s t converges approximately to Xoptimal, as generation g grows.


## 4. Conclusion

The Plant Propagation Algorithm (PPA) and its variants for continuous optimisation problems are getting notoriety as flexible and powerful solvers. PPA is a heuristic inspired by the way plants and in particular the strawberry plant propagate. It is also referred to as the Strawberry Algorithm. While there is a growing body of experimental and computational works that show its good behaviour and performance against well-established algorithms and heuristics, there is very little if any in terms of theoretical investigation. This gap, therefore, needs to be filled. The aim of course, in analysing the convergence of any algorithm (in this case PPA-C), is to give confidence to the potential users that the solutions that it returns are of good quality. The convergence analysis put forward in this paper relies on the exploitation and exploration characteristics of the algorithm. Since it does not get stuck in local optima and explores thoroughly the search space it is only a matter of time before the global optimum is discovered. The approach is probabilistic in nature and ascertains that the global optimum will be found with probability 1 provided the algorithm is run reasonably long enough. The argument for this to hold is that at each iteration new and better solutions are generated which means that, in the limit, the global optimum is reached. Questions still remain concerning what is considered a reasonable amount of time. Bounds on the time it will take to converge are being developed and results will be presented in a follow-up paper.



---

*Source: 6357935-2018-03-21.xml*