Natural Reforestation Optimization (NRO): A Novel Optimization Algorithm Inspired by the Reforestation Process

Department of Electrical Engineering and Information Technology, Paderborn University, Paderborn 33098, Germany Department of Mechanical and Industrial Engineering, Concordia University, Montreal H3G 1M8, Canada Faculty of Systems, Telecommunications and Electronics, Universidad de Especialidades Espíritu Santo, Samborondón 092301, Ecuador Department of Electrical and Computer Engineering, National University of Singapore (NUS), Singapore 117583

Similar to other meta-heuristic algorithms, NRO algorithm optimizes a defined problem by iteratively searching for the best solution. For this, it employs a population of particles (called trees) from which new particles are produced (called seeds) which are distributed within the search space (called forest). The amount of seeds each tree produces is based on the fertility of the place it was planted i.e., how good is its current solution. To distribute the seeds, they can, on the one hand, fall from the tree without the influence of the wind (these will then land close to their parent tree) or, on the other hand, leave their parent tree due to the wind influence in which the wind direction (that is influenced by the location of the current best solution) and speed are taken into account. The reduction of the wind speed due to other trees that block its path is also considered. Once the seeds have landed, they will grow and become trees and their associated solution is evaluated. The reforestation process is then repeated.
Based on the reforestation principle, our algorithm has the following characteristics:  It prioritizes the exploration of the search space in the beginning i.e., the number of seeds influenced by the wind (which can travel a long distance) is the highest at the first iteration and reduces as the number of iteration increases. These seeds are also useful to escape from local optimal regions  It prioritizes the exploitation of the search space towards the end i.e., the number of seeds not influenced by the wind (which travel only a short distance) increases as the number of iteration increases. The area range within these seeds can land is also reduced with increasing number of iterations  For the case of the seeds influenced by the wind, the search of the new solutions is affected by the location of the current best solution (which influences the wind direction) as well as the location of other trees (if close to the trajectory of the seeds) since they can reduce the wind speed based on how good is their associated solution. Therefore, the new population of particles (seeds) have the opportunity to interact with the current population (trees) in order to search for the best solution  The seeds that are not influenced by the wind are especially useful when the solutions are close to the global solution as they can perform a more intensive search within this area (exploitation). Particles located close to the global solution are expected to be obtained as the number of iterations increases and hence the number of these seeds increases (while the ones influenced by the wind decreases) as the number of iterations rises  The better the solution obtained by a tree, the more seeds it will be assigned so that more exploration and exploitation can take place in this area The rest of the paper is arranged as follows. In Section II, the principles of the natural reforestation process, from which our algorithm is inspired, are explained. Subsequently, the algorithm formulation is detailed in section III while the method to set the values of the variables are given in section IV. In section V, the performance of the NRO algorithm is compared with other optimization algorithms to solve different optimization problems. Section VI then concludes the paper. An Appendix section is presented at the end of this document describing the benchmark functions and providing the nomenclature table (Table 3).

Natural Reforestation Process
The NRO algorithm takes into account the following elements involved in the natural reforestation process: The forest, the trees, the seeds and the wind. Here, the forest represents the range of values in which the optimization variables are bounded. The trees represent the population of solutions whose location within the forest defines the value of their optimization variables. Their height is calculated from a fitness function based on the optimization problem. In case of a minimization (maximization) problem, the lower (higher) the objective function value for a particular tree, the greater its height. The seeds represent the new solutions that are utilized to search for better solutions within the forest. They are generated from the trees (the tree that generates a seed is referred as its parent tree). The place where they are planted can depend on several factors such as the location and height of the trees as well as the wind direction and speed.
The natural reforestation process is as follows: 1) The population of trees is planted among the forest 2) The taller the tree, the more seeds it can generate 3) The seeds produced by the trees can be dispersed with two options: a) No wind influence: Some seeds can fall without the wind intervention. These will land close to their parent tree b) Wind influence: The wind can separate the seeds from their parent tree and the place where they are planted depends on their parent tree's location, the wind direction and speed. The wind speed may be reduced due to other trees blocking its path 4) Based on the soil quality of the place where the seed is planted, it will grow into a tree with a certain height (the more fertile the land, the greater its height) 5) As the forest has a limited space, the trees will compete with each other for nutrients and sunlight and the ones with the greatest heights will remain 6) The process is then repeated

Algorithm Formulation
Based on the process described in the previous section, the steps of the NRO algorithm for an optimization problem are formulated below: 1) Initialization. The following parameters are defined: a) Objective function f obj which is in terms of the n var optimization variables x(k), k = 1, 2,…, n var b) Range of values the optimization variables can have xmin,(k) < x(k) < xmax,(k), k = 1, 2,…, n var . This will constitute the size of the forest c) Population size n pop d) Total number of seeds that will be generated at each iteration N seed e) Percentage of seeds that are planted without the wind influence, called internal seeds, at the initial iteration per seed,ini and at the final iteration per seed,fin f) Multidimensional sphere size within which the internal seeds are planted (the parent tree is the center of the sphere) at the initial iteration Sph seed,ini and at the final iteration Sph seed,fin . The sphere size is measured as a percentage of the maximum distance within the forest g) Minimum distance between the traveling external seeds and the trees in order for the latter (called blocking trees) to reduce the speed of the wind that is transporting the external seeds , min wind tree dist . This distance is measured as a percentage of the maximum distance within the forest h) Maximum distance that an external seed will reach if it is generated from the highest tree of the current iteration and with the highest wind speed (without the influence of blocking trees) max, max x , k = 1, 2,…, n var , p = 1, 2,…, n pop with the condition that no particle should be within the neighborhood of another particle (to assure diversity, similar to (Chelouah and Siarry, 2000)). Therefore, for any two particles with values of the optimization variables of   x , a  b, respectively, the following condition must hold (this is represented in Fig. 1a for a problem with two optimization variables) 1 : d) Among the current population, the tree with the best solution is selected (the tree closest to the optimal result). The value of the optimization variables and the objective value of this tree are stored in Bx(k) and BObj, respectively 3) Seeds production a) The number of seeds that each tree will produce   seed p n , p = 1, 2,…, n pop is based on their height (the greater their height, the more seeds they will produce) 3 : When the algorithm starts, the percentage of internal seeds is low with a value of per seed,ini while the one of external seeds is high. As the iterations go on, the percentage of internal seeds increases reaching a maximum of per seed,fin when n it = max 2 it n , while the one of external seeds decreases. This takes place to give emphasis to the exploration at the beginning and then to the exploitation towards 1175 the end of the algorithm. The number of internal seeds   , seed int p n , p = 1, 2, …, n pop from each tree is then calculated as: c) The number of external seeds for each tree   , seed ext p n , p = 1, 2, …, n pop is: In case any of the particles goes beyond the search space, they will be brought back to the searching boundary: 4 The function rand generates a random number between zero and one drawn from the standard uniform distribution. The  symbol indicates that it can be an addition or a subtraction (the operation is randomly selected). a) The external seeds are the ones dispersed from their parent tree with the wind influence. The location where they are planted depends on the wind direction and speed, the height of their parent tree and the blocking trees. The falling seeds are assumed to experience a two dimensional -consisting of x and y axisparabolic movement which follows the kinematic equations. Here, the gravity and height of the parent tree control the movement of the seeds in the y axis (height of the seed) while their movement in the x axis (travelled distance of the seed) is affected by the wind as well as the height and location of the parent and blocking trees (the direction of the x axis is the same as the wind direction). For simplicity, the rest of this section deals with the travelled trajectory of an individual seed. This procedure can then be repeated for the other external seeds b) Once the external seed is traveling, its travelled distance dist seed , when reaching a final position of x seed,fin , k = 1, 2,…, n var within the forest, is calculated based on its initial location x seed,ini , k = 1, 2,…, n var (location of its parent tree), as: c) To estimate the value of the gravity g, it is assumed that, when the external seed is generated from the tree with the highest height in the current iteration and the wind is at maximum speed, it will be able to travel a total distance of max dir , k = 1, 2,…, n var is a n var -dimensional vector due to the n var optimization variables. It is calculated taking into account the position where the tree with the best result is located, as follows: e) The initial wind speed is randomly selected between zero (no wind) and one (maximum wind speed): f) Based on the wind direction, the position of the trees that are close to its traveling path (the blocking trees) are stored. Among the total population, the blocking trees are the ones whose minimum distance between them and the traveling seed is lower than , min wind tree dist . This is represented in Fig. 1b for the case of a problem with two optimization variables g) Once all the blocking trees are found, the traveling seeds will reach a particular blocking tree (which is located at a position of   , tree block k x k = 1, 2,…, n var ) when it has travelled a distance of: Once the seed reaches a blocking tree, its speed   wind k speed , k = 1, 2,…, n var is linearly reduced based on the height of this tree h tree,block . The greater the height of the blocking tree, the higher will be the speed reduction. If the blocking tree is the one with the maximum height at the current iteration, it will reduce the wind speed by a factor of 100%-  In case any of the particles goes beyond the search space, they will be brought back to the searching boundary: The algorithm structure previously described is summarized in Fig. 2.  Fig. 1: (a) Neighborhood constraint for the definition of the initial population for a problem with two optimization variables x(1) and x(2). The solution s2 cannot be taken into account as it is within the neighborhood of s1 while s3 can be considered as it is not within this neighborhood. (b) Selection of the blocking trees based on their shortest distance to the seed traveling path. In this example, the minimum distance from the solutions s1 and s3 (d1 and d3, respectively) are assumed to be lower than , min wind tree dist so they are considered as blocking trees while the distance d2 from s2 is assumed to be greater and therefore it is not considered as a blocking tree. Because the analyzed seed will not travel towards s4 (due to the wind direction), s4 is not considered as a blocking tree. (c) Representation of the trajectory of the external seeds. The dashed lines show the locations of the blocking trees. At the beginning the seed has a height equal to the one of its parent tree and moves with a speed equal to the one of the wind speed. As it travels, its height decreases following a parabolic trajectory. Once it reaches the location of a blocking tree, the wind speed is reduced based on its height, changing the parabolic trajectory of the seed. The seed lands in the location where it reaches a height equal to zero a) x (2) x (2) x(1) Searching field (forest) Searching field (forest)

Parameter Specifications
Based on an empirical analysis in which the algorithm was tested under several objective functions, the assigned values to the constants defined in section III Algorithm Formulation are presented in Table 1. Despite the proposed values provided in this table, they can be further improved when dealing with a specific optimization problem.

Results and Discussion
The performance of the NRO algorithm is evaluated by employing a series of benchmark functions typically employed to test optimization algorithms. These are described at the Appendix section. Furthermore, four optimization algorithms are also tested to compare their results to the ones from NRO: Particle Swarm Optimization (PSO) (Eberhart and Shi, 2000) 6 , Continuous Genetic Algorithm (CGA) (Chelouah and 6 The value of the PSO parameters are the employed in this study when the Clerc's constriction method was applied. As no stopping criteria was formulated, the ones from [16] were considered. Siarry, 2000), Enhanced Simulated Annealing (ESA) (Siarry et al., 1997) and Enhanced Continuous Tabu Search (ECTS) (Chelouah and Siarry, 1999). These are selected due to the high number of optimization problems found in the literature which have been solved with them. The authors programmed the previous algorithms by following the descriptions of the cited papers so that a deeper analysis can be made. All the simulations were run on the super computer from the National University of Singapore which is composed by clusters with RAM of 48 GB and Intel Xeon X5650 of 2.67 GHz. Their results are then summarized in Table 2 and Fig. 3.
Because of the semi-randomness presented in the optimization algorithms (e.g., the location of the initial population is variable), each of them were evaluated 100 times so that a fair analysis could be performed. Table 2 presents the obtained outcomes of these algorithms with respect to the benchmark functions, each composed by a number of n var optimization variables. Only the successful results were taken into account and their averaged values are provided. A simulated result r sim is considered successful if it is close to the known/real analytical solution of the optimization problem r real , by fulfilling the following condition: where |r sim -r real | is the error of the simulation while 1 and 2 are the coefficients to determine the success condition, both were set to 0.1. Table 2 presents the percentage of successful results (per success ), the average number of times the objective function was evaluated (nf obj ), the average simulation time (time avg ) and the average absolute error (error avg ). As previously indicated, nf obj and error avg were calculated by considering only the successful results. Here, it can be appreciated that the NRO algorithm achieves the highest percentage of success for most of the benchmark functions. In comparison to the other optimization algorithms, the NRO shows a noticeable advantage when dealing with the Bea, Pow (n var = 10, 20) and Per (n var = 10, 20) functions. In general, it can be seen that, as the number of optimization variables increases, the performance of the algorithms tend to decrease (in particular for the PSO). Nevertheless, among them, the ESA is the one who achieves 100% of success rate for the SS, SDP and Sph regardless of the employed optimization variables. For a few functions (Pow with n var = 20 and Per with n var = 5, 10, 20) the NRO algorithm achieves values of less than 80% success rate. However, with respect to the Pow (n var = 20) and Per (n var = 10, 20) functions, it still reaches a better result than the ones obtained from the benchmark optimization algorithms being overcome only by the PSO and ECTS algorithms when evaluating the Per (n var = 5) function.      Multidimensional sphere size at the initial iteration time avg Average simulation time x(k) k th optimization variable xmax,(k) (xmin,(k)) Maximum (and minimum) allowed value of the k th optimization variable x seed,ext Landing position from the external seeds x seed,fin Location of seed after travelling x seed,ini Original location of travelling seed x seed,int Landing position from the internal seeds x tree Tree location x tree,block Blocking tree location The analysis of the algorithms based on their value of nf obj is required as, for complex problems, the evaluation of the objective function tends to represent a high (if not the highest) percentage of the simulation time. So, the lower this value is, the less computational time is required. All of the evaluated optimization algorithms reveal that as the number of its optimization variable increases, so does the nf obj value. This is an expected result because by the increase of n var , the problem becomes more complex and more evaluations are needed. Among them all, the NRO algorithm achieves the lowest nf obj values for most of the analyzed functions. For the Bea, Mat, CT, SHC, DP, Pow (n var = 5, 10), SS (n var = 2), SDP (n var = 2, 5, 10, 1181 20), Sph (n var = 2, 5) and Per (n var = 2, 20) functions, the NRO algorithm reaches nf obj values lower than the ones from the benchmark algorithms with a higher (or similar) success rate. This illustrates the high computational performance of NRO. Nevertheless, for particular cases such as the SS (n var = 20) and Per (n var = 20) functions, the NRO algorithm produces the highest value of nf obj . This, however, is required in order to achieve the highest success rate among the optimization algorithm (except for the ESA algorithm which achieved a greater success rate for the case of SS with n var = 20). As expected, the average simulation time has a similar trend as nf obj .
At this point it can be noted that, when the number optimization variables increase considerably, the NRO algorithm still achieves the highest success rate by sacrificing its computational time. This is different to the benchmark algorithms that will produce a lower simulation time with a lower success rate.
Based on the average error values, it can be appreciated that in general, the algorithms achieve very low values. This means that their successful results are close to the global minimum. It is worth mentioning that the ESA algorithm produces the lowest simulation errors for most of the benchmark functions reaching even values of zero. Figure 3 presents the performance of each algorithm based on the number of evaluated optimization variables. These graphs are obtained by averaging the results from Table 2 among the benchmark functions which contain the same number of optimization variables. Here, Fig. 3a reveals that, in general, the NRO algorithm outperforms the other algorithms with respect to the success rate regardless of the number of optimization variables (except for the case of n var = 5 in which the ECTS reaches a slightly higher value). It can also be seen that the PSO algorithm experiences the highest decrease in performance with the increase of n var . Furthermore, with respect to the number of times the objective function is evaluated, Fig. 3b shows that for n var = 2, 5, 10, the NRO achieves low nf obj values, similar to the other algorithms. Nevertheless, for the case of n var = 20, a high increase of nf obj is observed for NRO. This, as previously explained, is because the optimization problem is more complex with n var = 20 and more evaluations on nf obj are required to achieve satisfactory results (this can be seen in Fig. 3a). The average simulation time shown in Fig. 3c presents a similar trend as the one from Fig. 3b due to the high influence that the number of evaluations of the objective function has on the simulation time for the analyzed benchmark functions. With respect to the simulation errors, Fig. 3d reveals that all the successful results achieve low errors. Among them, the PSO algorithm reaches the lowest errors regardless of the number of optimization variables.

Conclusion
This work proposed a new meta-heuristic optimization technique for single objective optimization problems named Natural Reforestation Optimization (NRO). The NRO, as suggested by its name, is inspired by the natural reforestation process. The characteristics of this algorithm (which enhance the performance of the searching technique) were discussed in this study, among them: The initialization procedure, the exploration and exploitation process, the interaction between the particles, the stopping criteria, among others. The performance of this algorithm was tested with many benchmark functions against four other optimization algorithms. The results revealed that in general, the NRO is capable to generate a high success rate, reaching values close to the global optimum. It was also noted that when the number of optimization variables increased considerably (20 in this study), the success rate of the NRO algorithm would still be superior to the one from the benchmark requiring, however, a higher simulation time. Therefore, the potential of this algorithm was confirmed as a strong optimization tool opening the path for more work to be done to further enhance its performance. As a future work we will aim to further develop this algorithm so that it could also be applied for multi-objective optimization problems.

Appendix
The benchmark functions used in this study, shown in The Bea function contains two variables and is evaluated in the range of xi ϵ [-4.5, 4.5], for all i = 1, 2. Its global minimum is 0 and is located in x * = (3, 0.5).

Booth Function (Boo):
The Boo function contains two variables and is evaluated in the range of xi ϵ [-10, 10], for all i = 1, 2. Its global minimum is 0 and is located in x * = (1, 3). ( 7 The searching ranges of the optimization variables are obtained from the literature.

1182
The Mat function contains two variables and is evaluated in the range of xi ϵ [-10, 10], for all i = 1, 2. Its global minimum is 0 and is located in x * = (0, 0).
The Sch_N4 function contains two variables and is evaluated in the range of xi ϵ [-100, 100], for all i = 1, 2. Its global minimum is 0.292579 and is located in x * = (0, 1.25313).

Drop-Wave Function
The DW function contains two variables and is evaluated in the range of xi ϵ [-5.12, 5.12], for all i = 1, 2. Its global minimum is -1 and is located in x * = (0, 0).

1183
The Pow function contains n variables and is evaluated in the range of xi ϵ [-4, 5], for all i = 1,…, n. Its global minimum is 0 and is located in x * = (0,…, 0).

Sum Squares Function (SS):
The SS function contains n variables and is evaluated in the range of xi ϵ [-10, 10], for all i = 1,…, n. Its global minimum is 0 and is located in x * = (0,…, 0).

Sum of Different Powers Function (SDP):
The SDP function contains n variables and is evaluated in the range of xi ϵ [-1, 1], for all i = 1,…, n. Its global minimum is 0 and is located in x * = (0,…, 0).

Sphere Function (Sph):
The Sph function contains n variables and is evaluated in the range of xi ϵ [-5.12, 5.12], for all i = 1,…, n. Its global minimum is 0 and is located in x * = (0,…, 0).