Groundwater Flow Algorithm: A Novel Hydro-Geology Based Optimization Algorithm

A novel meta-heuristic nature-inspired optimization algorithm known as Groundwater Flow Algorithm (GWFA) is proposed in this paper. GWFA is inspired by the movement of groundwater from recharge areas to discharge areas. It follows a position update procedure guided by Darcy’s law which provides a mathematical framework of groundwater flow. The proposed optimization algorithm has been evaluated on 23 benchmark functions. The significance of the results is statistically validated using the Wilcoxon rank-sum, Friedman, and Kruskal-Walis tests. To prove the robustness of the algorithm, it has been further applied on several standard engineering problems. From these exhaustive experiments, it has been observed that the proposed GWFA can outperform many state-of-the-art optimization algorithms. Source code of this work is available at: https://github.com/Ritam-Guha/GWFA.


I. INTRODUCTION
Optimization deals with the process of searching for the best solution in a given scenario. An algorithm employed to find such a feasible solution is called an optimization algorithm [1], [2], [3]. The set of optimization algorithms can be broadly classified into two different categories: deterministic and stochastic. Deterministic optimization algorithms [4], [5] utilize analytical properties of data to search for the global optimum solution with theoretical guarantees. These procedures are used when it is absolutely necessary to find the global optimum solution. On the other hand, stochastic algorithms try to find near-optimal results through some approximations. A class of stochastic algorithms is called heuristic [6], [7]. Although heuristics do not guarantee finding the best feasible solution, they surely find a near-optimal solution within a reasonable time duration. In the case of NP-hard problems, it is almost impossible to find the best solution within a finite amount of time. That The associate editor coordinating the review of this manuscript and approving it for publication was Mehul S. Raval .
is why heuristic algorithms have become a very popular and widely-explored research topic in recent times to solve various time-constrained optimization problems.
Heuristic algorithms consider a particular problem as a black box with a set of inputs and outputs. The inputs are the variables of the problem and the outputs are the optimum solutions and their corresponding values. A heuristic search starts with creating a set of random inputs as the candidate solutions to the problem. The search is continued by evaluating each solution, noting the objectives and values, and modifying the solutions based on their current outputs. These steps are repeated until a termination criterion is met, which can be either a threshold on the maximum number of iterations or the maximum number of function evaluations. Regardless of the specific structure, at each stage, such algorithms require a way to compare two solutions to decide which one is better. An objective function or fitness function is used to evaluate the merit of each solution.
Although heuristics are very effective in solving real-life challenging problems, there are many difficulties when solving optimization problems [8]. Besides, optimization VOLUME 10, 2022 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ problems have diverse characteristics which are not similar at all. Therefore, it is a challenge for the researchers to design optimization procedures that are problem-independent in nature. Heuristic algorithms typically use problem-specific properties to find better approximations. To overcome such difficulties, a new class of heuristics has been created, which is known as meta-heuristics. The word 'meta-heuristic' has been first coined by Fred Glover [9]. It is a high-level problem-independent framework with some specific strategies used to devise optimization algorithms. In the last two decades, meta-heuristic algorithms have become quite popular in the literature [10] due to their (i) simple concepts and structures, (ii) almost derivation-free mechanisms, (iii) ability to handle local optima, (iv) flexibility, and (v) easy and effective hardware implementation. The new question that has surfaced in the field of meta-heuristic is what is the need for new meta-heuristic algorithms when there are so many existing algorithms that can solve the problems efficiently. A lot of researchers have also argued [11], [12], [13] that most of these algorithms are so similar in nature that their ''novelty'' is metaphorical. Some key points to note in this regard are: 1) Meta-heuristic algorithms are stochastic methods. So, these methods use the power of randomness to guide the solutions in a better direction. When the methods are not deterministic, there is no way to state that one algorithm will always perform better than the other algorithms, which eventually supports the same old fact as stated in the ''No Free Lunch'' theorem. It basically says that we cannot conclude that one algorithm is better than the other for all possible problems. So, it is a necessity to work on different approaches suitable for different problem domains. 2) Due to the stochasticity, meta-heuristic algorithms can never ensure finding the best solution in the search space. However, the empirical results show that they can reach a near-optimal solution in a reasonable time.
It is an effort to trade off the best solution with the time and computational needs. 3) Inspiration from nature is an important part of many meta-heuristic algorithms. Most of the ideas in the domain of automation are from the perspective of how human beings will do the same job, and then replicate the procedure using computational devices. The same concept applies to nature-inspired optimization as well. The researchers are looking for the natural phenomena which result in different optimizing procedures. The most important characteristic [14] of a good meta-heuristic algorithm is to find a good balance in its exploration and exploitation capabilities. With exploration or diversification property, the algorithm aims to find out the 'promising' areas where the global optima may lie, whereas the exploitation or intensification property is responsible for searching the areas already discovered in a more precise manner in order to pinpoint the solution. Exploration preserves the diversity and stops the algorithm to converge prematurely to some local optima. Exploitation helps the algorithm converge. Interestingly, many such meta-heuristic algorithms are intuitively inspired by nature. Simple natural elements have proved their skills to perform optimization in the most trivial way possible. As a result of this, researchers have developed various nature-inspired metaheuristic optimization algorithms such as Genetic Algorithm (GA) [15] inspired by genetic crossover and mutation, Particle Swarm Optimization (PSO) [16] based on flocking of birds, Ant Colony Optimization (ACO) [17] inspired from foraging behavior of ants, etc. This has motivated us to carefully observe natural occurrences and propose a novel meta-heuristic developed based on the flow of groundwater from recharge areas to discharge areas.
The contributions of this paper are: 1) Developed a novel meta-heuristic nature-inspired optimization algorithm called Groundwater Flow Algorithm (GWFA) inspired by the flow of groundwater. 2) Evaluated GWFA against 23 standard benchmark functions consisting of 7 -uni-modal, 6 multi-modal, and 10 fixed-dimensional multi-modal functions. 3) Compared GWFA with some classic and state-of-theart optimization algorithms used in the literature. 4) Explained statistical significance, convergence, and overall superiority of the algorithm in terms of obtained results. 5) Applied the algorithm over 5 engineering application problems to prove its robustness. The rest of the paper is organized as follows: Section II provides a brief overview of the research work being carried out in the domain of nature-inspired optimization, Section III presents the proposed algorithm in detail with proper explanation of its ability to handle exploration-exploitation balance, Section IV provides various test results obtained for the proposed algorithm and their comparison with some classic as well as recently proposed optimization algorithms, Section V finally concludes the manuscript and provides future scope of improvement.

II. RELATED WORK
As a modern trend in optimization, researchers across the globe are proposing various meta-heuristic algorithms to tackle different avenues of optimization problems. Metaheuristic algorithms can be divided into different categories based on varied criteria: single solution-based and population-based [18], nature-inspired and non-natureinspired [19], metaphor-based and non-metaphor based [20]. From the 'inspiration' point of view, these algorithms can roughly be divided into four categories [21]: Evolutionary, Swarm inspired, Physics-based, and Human related.
• Evolutionary algorithms are basically inspired from biology [22]. They utilize crossover and mutation operators for the evolution of the initial population, usually selected in a random fashion, over the iterations and eliminate the poor solutions, thereby ensuring the improved solution. GA is a well-known method of this category which follows Darwin's theory of evolution. Co-evolving algorithm [23], Cultural algorithm [24], Genetic programming [25], Grammatical evolution [26], Bio-geography based optimizer [27], Stochastic fractal search [28] etc. are some well-known evolutionary algorithms.
• Swarm inspired algorithms, on the other hand, imitate individual and social behavior of swarms, herds, schools, teams, or any group of animals [29]. Every individual has their own specific behavior, but the behavior of the collection of individuals gives us a way to solve complex optimization problems. A popular algorithm belonging to this category is PSO, devised by following the behavior of a flock of birds. Another important algorithm of this category is ACO, designed by mimicking the foraging process of the ant species. Several popular and recent algorithms of this category are: Shuffled frog leaping algorithm [30], Bacterial foraging [31], Artificial bee colony [32], Firefly algorithm [33], Grey Wolf optimizer (GWO) [34], Ant Lion optimizer (ALO) [35], Whale optimization algorithm [36], Grasshopper optimization algorithm (GOA) [37], Squirrel search algorithm [38], Harris Hawks optimization (HHO) [39] etc.
• Physics-based algorithms are inspired following the rules used to govern a physical process [40]. Some inspiring physical processes are based on metallurgy, physics, chemistry, complex dynamic systems, and even music. The oldest one belonging to this category is Simulated Annealing (SA) [41], formulated by following the annealing [42] process of metals in metallurgy and materials sciences. Another popular method of this category is the Gravitational search algorithm (GSA) [43], designed by following gravity and mass interaction. Besides, there are other methods belonging to this category which include Self-propelled particles [44], Harmony search (HS) algorithm [45], Black hole optimization [46], Sine Cosine algorithm [47], Multi-verse optimizer [48], Find-Fix-Finish-Exploit-Analyze [49], Hydrogeological cycle algorithm [50], etc. Rubio et al. has also provided an overview of recently developed water-related metaheuristics in [51].
• Human related algorithms look for the global optima by following human behavior [52]. Teaching-Learning-Based optimization [53] is a popular one belonging to this category which is formulated by following the improving procedure of class grade. Apart from this, some other important methods of this category include Society and civilization [54], League championship algorithm [55], Fireworks algorithm [56], Tug of war optimization [57], Volleyball Premier League algorithm [58]. The presence of such a significant number of meta-heuristic algorithms clearly raises the question of the need for another meta-heuristic algorithm. However, as indicated by No Free Lunch [59] theorem for optimization, there cannot be any single algorithm that will be equally applicable for all the optimization problems desiring optimal solutions. With each new algorithm following any regular or natural phenomenon, researchers primarily aim to provide some new facet to the algorithm where both exploration and exploitation will have a superior trade-off, thereby trying to get away from the local optima and eventually compass to the global optima. Nevertheless, accomplishing these objectives is not straightforward, hence motivating researchers to propose a new algorithm that can be applied to different problem domains. It is to be noted that many state-of-the-art feature selection algorithms follow the same architecture as PSO. Krill Herd (KH) algorithm [60] and GWFA also follow that architecture, but they are not the same algorithm. The commonalities are visible because all these algorithms start by initializing candidate solutions, updating them over the iterations, and finally providing a result. However, this is considered to be a very good architecture and it has proved its effectiveness in multiple scenarios. That is why GWFA has been implemented following that architecture, however, it has many traits which are different from both PSO and KH.

III. GROUNDWATER FLOW ALGORITHM
In this section, a detailed outline of the proposed algorithm is provided which is followed by the mathematical description of the optimization procedure.

A. HYDRO-GEOLOGICAL PERSPECTIVE/INSPIRATION
Three-fourths of the Earth's surface is covered with water. Most of this water can be seen with bare eyes when people go to beaches or stand at the side of rivers or ponds. But still, there is water that cannot be seen unless a hole is dug up in the ground deep enough to reach the water tables. This water is known as groundwater which accounts for 1.7% of all of Earth's water and 30.1% of the freshwater [61]. The depth of ground water varies from place to place. In some places (like marsh), water may occur very close to the ground surface whereas, in some other places, it may appear hundreds of feet below the ground level (as in deserts). Water nearing the ground surface may be a few hours old, at moderate depth, it may be a few hundred years old and the water nearing the earth's core may be some thousand years old flowing from one place to another underground. The level of water underground is known as the water table.
This water occupies the pores in the soil and moves from one place to another. But not 100% of the pores are occupied by water, some space may be occupied by air as well. Especially the pores of the soil closer to the Earth's surface are significantly occupied by air. This part of the soil is called the unsaturated region whereas the part where all the pores are occupied by water is known as the saturated region. In the final stage of the water cycle, i.e. precipitation, water falls on the Earth's surface as rain, snow, hail, etc. The precipitated water either flows over the surface following the slope and reaches a stream, or it infiltrates the ground and enters into the unsaturated region. Some of this water gets used by plants or evaporates due to heat but the rest of the water goes deeper into the ground and recharges the groundwater. The area where the precipitated water flows past the unsaturated zone to the water table is called the Recharge area (RA) while the areas where the groundwater flows to (streams, lakes, etc.) are called Discharge Areas (DAs). Groundwater flow is known as the flow of groundwater from RAs to DAs. The concept can be made clear by looking at Figure 1. The reason for this flow is typically the gravitational force. Water at RAs has higher potential energy due to their high elevations. So, they move to DAs with lower elevation by overcoming the internal frictions (determined by viscosity). Due to the gravitational pull, groundwater slowly flows through layers of rocks, soil, and sand which are known as aquifers (originated from two Latin words 'aqua' or water and 'ferre' or carrier). The flow of groundwater has always been of great interest to mankind apparent from many Biblical references as well. Finally, in 1855 and 1856, a French engineer named Henry Darcy performed various experiments [62] and was able to mathematically define the flow of groundwater which has been named the Darcy's Law. Darcy defined the flow of water through sand beds during his experimentation but later it has been generalized for other fluids and porous mediums as well.
In this work, the flow of groundwater is mathematically modeled to solve optimization problems. Darcy's law helps in the movement of searching agents efficiently while performing optimization.

B. MATHEMATICAL MODEL
The most important aspect of the proposed optimization approach is the position update policy which is guided by Darcy's law. In this section, the mathematical definitions of Darcy's law are discussed.
During experimentation, Darcy discovered that the velocity with which groundwater flows is dependent on two major factors: height difference and gap in positions. Darcy's experimental setup has been presented in Figure 2. Consider that groundwater is flowing from point A to point B as shown in Figure 2. The velocity of this flow is directly proportional to the height difference ( h) between the two points but inversely proportional to the length of the gap (L) between them. The ratio of these two factors (i) is termed the Hydraulic Gradient (HG).
HG clearly suggests that when the height gap between the points is more and the length between them is shorter, the water flows more rapidly. The velocity is directly proportional to the HG. So, the velocity term can be expressed as: where k is the constant of proportionality also known as the coefficient of permeability. This velocity is known as the discharge velocity. It is an idealistic term that considers water flows through the entire soil between the two points but it is not realistic as water only flows through the pores. So, there is another term used for defining velocity, called seepage velocity, which is the velocity of the water flowing through the pores of the soil. It can be represented as: where φ is called the porosity of the soil. Porosity is the ratio of the combined volume of pores to the entire volume of soil. Seepage velocity is ultimately the velocity with which the groundwater flows from one place to another. Depending on the value of the porosity (property of the soil), the seepage velocity may vary. The proposed model uses seepage velocity to update the position of the candidate solutions.

C. OPTIMIZATION APPROACH
In any population-based meta-heuristic optimization algorithm, it starts with some candidate solutions which act as the search agents. These search agents traverse various paths in the search space in pursuit of better solutions. GWFA uses some groundwater sources (GWSs) as its search agents from where water flows in different directions to search for the global optimum solution. The optimum solutions are represented by the DAs where the GWSs are guided. Hence, the objective here is to look for the most optimal DA by flowing groundwater from the RAs. It is to be noted that in this section GWSs, agents, and candidate solutions are used to represent the same thing. In a broader sense, we can say that each GWS is a search agent which carries a candidate solution. The algorithm consists of three separate phases: initialization, position update, and exploitation of DAs. A flowchart describing the entire algorithm is presented in Figure 3.

1) INITIALIZATION
The GWSs act as the search agents in the proposed optimization approach traversing various paths in search of a solution to the problem. At any stage of the algorithm, the solution represented by any GWS is called a candidate solution. At the first stage of the algorithm, a certain number of GWSs are initialized with random candidate solutions. Suppose n GWSs are initialized and the function under consideration describing the optimization objective (say f ) has D-dimensional design space. Then the GWSs will be represented as: where x ij represents the value of i th candidate solution in j th dimension, ub j (f ) and lb j (f ) are the upper bound and lower bound of design space of function f respectively and rand j (1) is a random number restricted in [0, 1]. In this way, n GWS positions are initialized which act as the initial solutions to the optimization problem. For the rest of the discussion, x i (t) represents the i th GWS in the population in t th iteration of the algorithm.
Each candidate solution is guided by a velocity term. The velocity of every individual GWS is also initialized with 0.
Here v i (0) represents the initial velocity of i th GWS.

2) POSITION UPDATE
This is the most crucial part of the proposed algorithm. After getting the initial GWSs, their positions are updated based on groundwater flow rules and abiding by Darcy's law. At the beginning of the position update policy, the quality of all the solutions is calculated using the objective function. It is to be noted that the objective function in an optimization problem is the function whose value needs to be minimized (minimization problems) or maximized (maximization problem). The objective function provides a fitness measure for the solutions. The lower (or higher) the function value is for the minimization (or maximization) problem, the fitter the solution is. After getting quality measures for each and every GWS, the fittest p% of them are assigned as DAs. The rest of the GWSs are then considered to form RA. The objective for the GWSs in RA is to approach one of these appropriate DAs, thereby improving the fitness of the candidate solutions.
Apart from the DAs, if there are GWSs in the RA which is fitter than the others then the less fit GWSs will also get guided toward the fitter GWSs. Intuitively it can be considered that fitter GWSs in the RA are at lower elevations in RA. Hence, water from higher GWSs keeps flowing toward these GWSs. The movement of water will be guided by three important factors: 1) Previous flow of the water.
3) Flow towards the fitter GWSs in the RA. For each of the GWSs in the RA, one DA is selected randomly as shown in Equation 6. But instead of defining a flow to every fitter GWS in RA, an average of the lower GWSs is constructed to reduce computation. The solution representing GWS is then guided towards the selected DA and local averaged GWS (LA) using Darcy's law of groundwater flow. From the discussion provided in Section III-B, it can be observed that groundwater is mainly guided by two important terms: height difference ( h) and length of gap (L). In the proposed optimization approach, h is calculated as the difference in the current candidate solutions and L represents the ratio of the maximum difference in the functional values in the population and the difference in their current functional values. The mathematical representations of these terms are shown below: where DA i (t) and LA i (t) represent the selected DA and the local averaged GWS for i th GWS.
Here the subscripts d and l correspond to the DA and LA counterparts respectively.
HG for each GWS is then calculated as the ratio of h and L terms. Darcy's law states that the discharge velocity is directly proportional to the HG. Therefore, by multiplying the HG with a constant of proportionality (k), discharge velocity (vd) for every individual is computed. k is also known as the coefficient of permeability.
The discharge velocities directed towards the DA and LA are combined using Equation 16. The weightage provided to vd id (t) is α while vd il (t) gets the weightage of 1 − α where α is in [0, 1]. The value of α is determined through experimentation provided in later sections. α controls the level of importance of guidance provided by DA and LA.
Discharge velocity is velocity of the water when there is no obstruction along the flow. In general, there are porous mediums known as aquifers in the path of the groundwater. So, the water only flows through the pores of the aquifers. Hence, a new term is introduced which is known as the seepage velocity (vs). Seepage velocity is determined by dividing the discharge velocity by the porosity of the intermediate medium. The expression for seepage velocity can be represented as: where φ is the porosity of the medium. φ is a very special term in the algorithm because it helps to control the exploration-exploitation balance in the process. When the position of a GWS gets updated, all the dimensions of the solution are not updated. The number of dimensions to be updated (d) is guided by a factor called the control factor (CF) which is changed over the iterations as: where T is the maximum number of iterations. The number of dimensions to be updated is selected by randomly choosing an integer d in the range [CF(t), D]. As evident from expression 18, the value of CF(t) decreases as the iterations progress. Therefore, in the beginning, there are more chances of a high number of dimensions getting updated (i.e. ensuring exploration) compared to the later stages (i.e. ensuring exploitation). The value of φ can be expressed as: Finally, the velocity of any individual GWS can be calculated by taking the combination of its previous velocity (learning from past experience) and its seepage velocity (learning from recent experience). The contributions of both factors to the velocity is the same. The velocity can be represented as: Here β = γ = 0.5. The value of the velocity is used to update the positions of the GWSs (i.e. candidate solutions). For each individual, d random dimensions are selected from the entire set. The values of these positions are modified according to the velocity: After the update, it may so happen that the GWSs get converged. At this situation, every GWS will represent the same candidate solution. It can be identified when all the GWSs are having same fitness measure. If such a situation occurs, the algorithm stores the best GWS found till now, discard all the solutions and re-initialize every GWS in the same way as discussed in the Initialization section.

3) DA EXPLOITATION
In this procedure, the GWSs flow towards DAs. So, if the DAs are not changed over iterations, the solutions will automatically get converged. In order to avoid that scenario, the DAs perform exploitation by checking their neighborhoods. Some dimensions of every DA are selected and their values are changed to find a neighborhood solution. In this way, for each DA, 5 neighborhood solutions are constructed and the one with the best fitness measure is selected to replace the older and less fit DA. If there is no fitter neighbor, DA is not changed. If at a certain point, the solutions get converged, the DAs do not participate in the exploitation and they also get re-initialized.

a: EXPLORATION-EXPLOITATION BALANCE
One of the major concerns in any optimization procedure is to maintain a good balance between the exploration and the exploitation of the search. It is a dilemma requiring the process to make a choice between searching an unexplored area of the search space and searching the same area with deeper insights. Hence, one needs to make an appropriate decision at this point to proceed. GWFA tries to control the trade-off in a very effective way. GWFA uses not one but multiple DAs which are treated as the global best of the current iteration which are followed by the other GWSs. For each GWS, a DA is selected randomly. It gives rise to the exploration of the search space because more avenues are to be looked for. Had it been only one DA, all the other GWSs would have followed its path, and finally the searching procedure would have found limited or only a small region of the search space. The next important term which helps in achieving the balance is the CF. As discussed before, CF supervises how many dimensions of the GWSs get updated by the velocity. At the beginning of the algorithm, there are chances of many dimensions getting modified but during later stages, the number of dimensions gets updated is restricted to a lower value. More updates of dimensions lead to exploration because it helps to inspect new areas of the search space. On the other hand, when the lesser number of dimensions get updated, only nearby neighborhoods are checked, thus enhancing the exploitation capability of the algorithm. In the later stages of the algorithm, DAs participate in exhaustive exploitation to intensify the exploitation in the search space. Therefore, as evident from the discussion, GWFA has the ability to maintain a good trade-off between the exploration and the exploitation phases. Overall TC: = O(ndI ) From the discussion, it can be seen that GWFA requires polynomial time for computation. This proves that it is a time-efficient algorithm and is applicable to optimization.

D. NOVELTIES OF GWFA
A literature survey reveals that nature-inspired algorithms have gained enormous popularity over the last decade. However, there is a reason why researchers find inspiration and come up with new optimization algorithms. First of all, as suggested in [63], researchers often use ''No Free Lunch'' theorem to justify the need for new optimization algorithms as no single optimization algorithm can work efficiently for every optimization problem. Thus, novel algorithms diversify the usage of optimization for different kinds of problems. Although most of the nature-inspired optimization algorithms follow a common framework mentioned in [64], they are derived from separate natural phenomena and use a different method to adjust the exploration-exploitation trade-offs. Secondly, these algorithms have proved their effectiveness in various scenarios. When an idea from a different field of research is transferred to another field of application, it actually tries to support interdisciplinary activities, This eventually widens the scope of the research, as many people from varied domains come up with their very own and unique ideas to fit such algorithms into different research problems. The novelties of the proposed algorithm can be expressed using the following points: • Hydro-geological Perspective: GWFA is inspired by the natural flow of groundwater from RAs to DAs. The motion is guided by Darcy's law of groundwater velocity. This hydro-geological perspective behind the formation of GWFA makes it very intuitive and innovative.
• Distributed Elitism: Most of the nature-inspired algorithms use a single elite candidate solution to guide the other solutions. However, in GWFA, every RA (solution) is guided by a randomly selected DA (elite solution) and other RAs which are better than the current RA. This step clearly expands the search area of the individual RAs and helps them circumvent local optima. When the candidate solutions keep on following a single elite solution, their probability of premature convergence becomes very high.

• Quality-based Exploration-Exploitation Balance:
The key to a successful search over the large search space is the trade-off between exploration and exploitation of the searching algorithm. Usually, this is achieved by adding a local search procedure to enhance exploitation or a global searching technique to bring in exploration. However, GWFA achieves this by performing two simple operations. First, for RAs, the CF is used to support exploration in the first phases of the algorithm and exploitation in the later stages. Secondly, the algorithm uses a very simple local search technique to exploit the DAs to find better solutions. This kind of updated policy is inspired by a very natural hypothesis. The scope of improvement for a level 50 agent is not as much as the scope of improvement for a level 5 agent. Hence, instead of trying to improve the best solutions, GWFA tries to improve the poorer solutions through exploration in the first phase and exploitation in later phases and the for i=1 to n do 5: Calculate fitness value of each GWS as fit(i) = f (x i (t)). 6: end for 7: Sort the GWSs according to increasing (or decreasing) order of fitness values for minimization (or maximization) problem. 8: if fit(1) == fit(n) then 9: Store the best GWS, re-initialize all GWSs. 10: t ← t + 1.

12:
end if 13: Select the top p% (or j = floor( n×p 100 )) of the GWSs as DAs.
14: Create the list of dimensions to be updated list i . 25: 26: else 27: Create 5 neighbors from x i (t) through perturbation 28: if Fitter neighbor found then 29: Replace the fittest neighbor with x i (t) 30: end if 31: end if 32: end for 33: end while best solutions are improved by a simple local search, So, we can see that more efforts are provided to improve the poorer solutions as compared to improving the best solutions. All these factors make GWFA very unique and better than many algorithms in the field. The results provided in Section IV clearly demonstrate this superiority over other well-established algorithms.

IV. EXPERIMENTATION
This section includes the performance of the proposed optimization algorithm called GWFA over some standard benchmark functions and a comparative study with some other state-of-the-art methods. The section is further divided into six sub-sections. The first sub-section contains the demonstration of the benchmark functions. The second sub-section comprises parameter tuning experimentation for GWFA. The third sub-section reports the obtained results and a detailed comparative study with other optimization methods. The next two sub-sections consist of testing outcomes for convergence and statistical significance respectively. The final sub-section includes the outcomes of the proposed algorithm when evaluated on five standard engineering applications.

A. BENCHMARK FUNCTION DESCRIPTION
In order to evaluate the proposed approach, 23 benchmark functions stated in [65] are selected. All the functions are exactly the same, except we have used 'F', instead 'f ' to denote a function (e.g. we have used F 1 , instead of f 1 ). The 23 benchmark functions can be broadly classified into three categories: 7 uni-modal functions (F 1 to F 7 ), 6 multimodal functions (F 8 to F 13 ), and 10 multi-modal functions with fixed dimensions (F 14 to F 23 ). Please refer to [65] for the detailed description of the uni-modal, multi-modal, and multi-modal functions with fixed dimensions. The uni-modal functions have only one optimal value, and obtaining this value helps estimate the exploitation ability of any algorithm. The multi-modal functions have more than one optimal value, which essentially tests the exploration ability of the algorithm. The final category is similar to the multi-modal functions but with small dimensions.

B. PARAMETER TUNING
The proposed algorithm has two control parameters namely α and k. This section mainly focuses on finding the most optimal parameter combination. For this purpose, two functions from each category are selected: F 1 , F 6 from uni-modal functions, F 9 , F 11 from multi-modal functions, and F 20 , F22 from multi-modal functions with fixed dimensions. The values considered for each parameter are as follows: α = [0.5, 0.55, 0.6, 0.65, 0.7] and k = [0.4, 0.5, 0.6, 0.7]. As a result, a total of 20 (5 * 4 = 20) combinations are checked. During the testing, the number of runs, GWSs, and iterations remains fixed as 30, 50 and 1000 respectively.
The best values obtained for each function are recorded and summed for every parameter combination. The best parameter combination is the one that has produced the best sum value. A visual interpretation of the obtained results is provided in Figure 4. From the experimentation, the best values of α and k are found to be 0.65 and 0.7 respectively. The final values for all the parameters are provided in Table 1.  Even though it is always a good idea to tune the parameters of an algorithm for a particular application, there are certain intuitions that can be used to tune them in an appropriate way, thereby reducing the effort. The parameter α determines the importance of the DA and the LA in the movement of the current RA. If it is known that the function is uni-modal, it would be better to provide more importance to the DA and the value of α should be increased. On the other hand, if there is no information about the functions in hand, its value should be closer to 0.5 (equal importance to the DA and the LA). The value of p denotes the percentage of DAs in the population. This number is very useful in guiding candidates in multiple directions. It helps when we deal with multi-modal functions. Hence, this value should be increased when it is known that the underlying function is a multimodal function, else it should be reduced. The value of 20 is a reasonable choice to maintain the appropriate diversity and is recommended for use when the user has no idea about the test function. Another important parameter in our algorithm is k. It basically converts the hydraulic gradients to velocity information and is a scaling mechanism. We have found the value 0.7 through experimentation and readers can use this value to find the scaling required for the application function.

C. RESULTS AND COMPARISON
The proposed algorithm is evaluated over 23 benchmark functions mentioned in Section IV-A. For proper experimentation, the algorithm is run 30 times, and the best, average, and standard deviation of the obtained values are stored. The result obtained for GWFA is compared with other state-of-theart methods that include GSA, PSO, Gravitational Particle Swarm (GPS), PSOGSA, Grey Wolf Optimizer (GWO), Coyote Optimization Algorithm (COA) and Equilibrium Optimizer (EO). To maintain a fair comparison environment, population size and the number of iterations for each algorithm is taken as 50 and 1000 respectively (same as GWFA). The best value, average value, and standard deviation are also calculated for every method. All the compared methods have some manually tuned parameters. The values of those parameters for each algorithm are tabulated in Table 2. The detailed comparative study is separately carried out for all three categories of functions. Tables 3, 4 and 5 contain the results for uni-modal functions, multi-modal functions and multi-modal functions with fixed dimensions respectively. A method is said to have the best result for a particular function if it has obtained the lowest best value. If multiple methods have the same best values, average values are checked and if still some methods are indistinguishable in terms of average value, finally standard deviation is checked to decide the best one.
The uni-modal functions mainly test the exploitation ability of any algorithm. As it is evident from Table 3, the proposed GWFA has outperformed all the other algorithms in six (F 1 − F 4 , F 6 , F 7 ) out of seven uni-modal functions. Not only that, the proposed algorithm has also surpassed all the other algorithms in terms of average value and standard deviation for these six functions, which proves the superiority of the algorithm. Though the proposed algorithm has not achieved the best outcome in F 5 , the result is indeed competitive with others. Hence, based on the properties of the uni-modal functions, it can be stated that GWFA is capable of handling exploitation in a competent way. VOLUME 10, 2022 Multi-modal functions are predominantly used for analyzing the exploration property due to its large number of local optima. Tables 4 show the outstanding performance of the proposed GWFA. GWFA has achieved the best result in three (F 8 , F 9 and F 11 ) out of six functions not only in terms of best value but also in terms of average value and standard deviation. GWFA achieves the second best result in F 10 and the result difference is marginal. Obtained outcomes in the other two functions (F 12 , F 13 ) are very close to the global optima. The performance in multi-modal functions clearly reveals the strong exploration ability of GWFA. Apart from GWFA, EO is the next best performer with 5 best results (2 multi-modal and 3 fixed-dimensional multi-modal functions) and SOGWO is the third best performer with 4 best values (3 multi-modal and 1 fixed-dimensional multi-modal functions). Table 5 shows the result of GWFA in multi-modal functions with fixed dimensions. In this category, GWFA has achieved impressive results as well. It has outperformed others in three (F 17 , F 20 , F 22 ) out 10 functions. If we consider only the best values, GWFA is able to achieve the best result for 7 out of the 10 functions. In the case of F 16 and F 18 , the difference between the best result and GWFA result is only in terms of standard deviation. So, it is a clear indication that although GWFA has not achieved the best results in certain situations, it has given a very close competition to the best method in such situations.
In total, the proposed GWFA has achieved the best results for 12 out of 23 (approx 52%) benchmark functions. The outcomes in other functions are also very close to global optima and competitive with other algorithms. Excellent performance on both uni-modal and multi-modal functions also establishes the strong exploration and exploitation abilities of GWFA.

D. CONVERGENCE TESTING
For any optimization algorithm, convergence is one of the most important criteria for comparison. A good optimization approach should be able to provide fast but steady convergence, which means that the approach should provide fast reduction in the objective score, but it should not oscillate too much. That is why a convergence test is necessary for any optimization algorithm to prove its efficiency. This section describes the convergence capability of the proposed GWFA. Figures 5, 6 and 7 demonstrate the convergence graphs for uni-modal functions, multi-modal functions, and multi-modal functions with fixed dimensions respectively. There are two sub-parts of each convergence graph: parameter space and objective space. Parameter space provides the visual representation of the functions in 2D, while the objective space contains the convergence curves. The X-axis in objective space indicates the iteration number, whereas the Y-axis depicts the global best value among all populations obtained till the current iteration.
The number of GWSs is kept at 50 and 500 iterations are used to test the convergence capability of GWFA.
For comparison of convergence, we have selected the top three algorithms according to Section IV-C namely GWFA, EO, and SOGWO. The convergence experiments clearly indicate that GWFA has the ability to provide steady convergence. However, it is also clear that EO and SOGWO are good performers, in terms of convergence. GWFA has provided marginally better convergence compared to both its competitors over the uni-modal functions. The trend remains almost the same in the case of multi-modal functions, except for function F 8 , where SOGWO provides better convergence. But in the case of fixed-dimensional multi-modal functions, EO displays the best convergence capabilities in most cases.

E. APPLICATION OF GWFA IN ENGINEERING PROBLEMS
Engineering problems are some real-world problems that involve designing and building systems and/or products. It is basically a decision-making process that contains complex objective functions and a large number of decision variables. Generally, meta-heuristic algorithms perform better than any traditional algorithm due to the proper convergence of the former to an optimal solution. Besides, meta-heuristic algorithms also have the ability to handle non-convex and non-differentiable functions. Mainly complex engineering problems involve a large number of design variables. The influence of those variables on the objective function, which is to be optimized, turns out to be very troublesome and non-linear in nature. Five classical engineering design problems namely, spring, gear train, welded beam, pressure vessel and closed coil helical spring design are considered in this study. The problems are explained in this section. These problems have many local optima, but the task of the algorithm is to find the global optimum. Hence, an efficient optimization algorithm is required to solve the problems.

1) TENSION/COMPRESSION SPRING DESIGN PROBLEM
In this problem, the main objective is to minimize the weight of the coil involving three decision parameters -wire diameter d, mean coil diameter D, and the number of active coils N along with four inequality constraints. The objective function is given in equation 23. Any solution to this problem can be represented as: The visual description of this problem is provided in Figure 8.

2) GEAR TRAIN DESIGN PROBLEM
The aim of this design problem is to minimize the cost of the gear ratio of the gear train. The four decision parameters present in this problem are T a , T b , T d and T f . There are no inequality constraints present in this problem. A pictorial representation of the gear train design problem is given in Figure 9. The gear ratio is formulated as     On the other hand, the objective function is demonstrated in equation 24. Representation of a solution can be given as:

3) WELDED BEAM DESIGN PROBLEM
The welded beam design problem is mainly a minimization problem consisting of four variables. The variables are weld thickness (h), length of the bar attached to the weld (l), bar's height (t), and bar's thickness (b). The available constraints for this problem include bending stress (θ), bean deflection (δ), shear stress (τ ), buckling load (P c ), and other side constraints. The problem is displayed in Figure 10.

4) PRESSURE DESIGN VESSEL PROBLEM
The main objective of this problem is to minimize the manufacturing, welding, and material cost of the pressure vessel. Four variables associated with this problem are -the thickness of shell (T S ), the thickness of the head (T h ) which are discrete decision variables, the inner radius (R), and length of the cylindrical section of the vessel (L) those are continuous decision variables. An image is provided in Figure 11 for better visualization of the problem. Any population is formulated as:

5) CLOSED COIL HELICAL SPRING DESIGN PROBLEM
The main aim of this design constraint problem is to decrease the volume of closed coil helical springs. A helical spring consists of closed coiled wire having the shape of a helix and is intended for the tensile and compressive load. The two variables related to this problem are as follows -coil diameter(D) and wire diameter(d). There is another parameter, the number of coils (n), which can be set beforehand. The pictorial description of the problem is given in Figure 12. Any population of this problem can be devised as  The objective function is formulated in Equation 27.
All the experiments are carried out using the optimal values of the two parameters -α and k. The number of iterations is kept at 1000 and the number of runs is 30. The detailed result obtained on the five engineering design problems is given in Table 6. The best value, average value, and standard deviation of all runs are also provided in Table 6. The outcomes are compared with six other state-of-the-art methods -PSOGSA, GPS, PSO, GWO, SOGWO, and EO. It is to be noted that the proposed GWFA algorithm is able to achieve the lowest value for three out of five engineering problems. For the rest two functions, GWFA has achieved results very close to the global optimum. Thus, this test clearly indicates that GWFA is highly robust in nature and is applicable to a large variety of mathematical optimization problems.

F. STATISTICAL ANALYSIS
For statistical analysis of the results of the proposed method, we have relied on three methods namely, Wilcoxon ranksum test, Friedman test, and Kruskal-Walli's test. Broadly, in any statistical test, we compare samples from two or more groups to determine how the independent and the dependent variable relate to each other. These statistical tests consider a null hypothesis of no relationship between the groups. After that, they find whether the observed data fall outside of the range of values predicted by the null hypothesis. The p-value (probability value) determines how probable it is that you would perceive the change defined by the test statistic if the null hypothesis of no relationship were true. If the value of the test statistic exceeds the statistic calculated from the null hypothesis, then it can be inferred that there is a statistically significant relationship between the predictor and outcome variables. If the value of the test statistic is less than the one calculated from the null hypothesis, then it can be inferred that no statistically significant relationship between the predictor and outcome variables exists. In general, statistical tests are divided into two types, namely, parametric and nonparametric tests. Non-parametric tests do not make as many assumptions about the data and are useful when one or more of the common statistical assumptions are not satisfied. So, non-parametric tests are considered in our work.  Firstly, a non-parametric statistical test, known as Wilcoxon rank-sum test [66], has been performed. This is done in order to check whether the results of an algorithm are statistically different from other algorithms [67]. The null hypothesis states that if two sets of results are from the same distribution, any difference in the two mean ranks comes only from the sampling error. If the distributions of the two results are statistically different, then the generated p-value from the test statistics will be < 0.05 (level of significance), as we have performed the test at 0.05% significance level, resulting in the rejection of the null hypothesis. To obtain the p value, we have considered two lists obtained from the proposed algorithm and the comparing methods for each function. The test is performed for 23 benchmark functions and 5 engineering applications. The final results of the Wilcoxon test are provided in Table 7.
Kruskal-Walli's test [68] is also performed to prove the consistency of the proposed GWFA. This is a non-parametric method that compares two or more groups of results obtained. To obtain the p value, we have considered two lists obtained VOLUME 10, 2022  from the proposed algorithm and the comparing methods for each function. This test is conducted at 0.05% significance level for 23 benchmark functions and 5 engineering applications and the obtained results, reported in Table 8, prove the significance of the proposed method.
Friedman test [69] is also performed for determining the statistical significance of the results obtained by GWFA. Friedman test is a non-parametric statistical test that provides the measure of the differences between multiple methods. It is normally used for answering the following kind of question: in a set of n number of results (where n ≥ 2), does at least two of the sets represent results with different median values? The null hypothesis considered here is no significant difference in the results of the methods at 0.05% significance level. To obtain the p value, we have considered two lists obtained from the proposed algorithm and all the comparing methods for each function. Table 9 contains the p values obtained from 23 benchmark functions and 5 engineering applications. From the test results provided in Table 9, all the obtained p-values are < 0.05. This clearly proves the presence of at least one set of significant results.

V. CONCLUSION
In this paper, we have proposed an optimization algorithm based on a hydro-geological phenomenon known as groundwater flow. The algorithm, named GWFA, is inspired by the movement of groundwater from recharge areas to discharge areas following Darcy's law. The Exploration and exploitation abilities of the algorithm are carefully taken care of through the usage of the control factor. From the experimental outcomes and associated discussion, we have shown the applicability of GWFA for mathematical optimization. The algorithm has been tested over 23 benchmark functions and 5 classical engineering problems. Some classical and recently published algorithms are compared with the proposed algorithm. GWFA has been able to outperform the other methods over 12 out of the 23 functions. Besides, significant tests are conducted to ensure the statistical significance of the proposed algorithm. Therefore, by analyzing the results and related discussion, it can be concluded that GWFA is comparable to most state-of-the-art algorithms and applicable to a variety of optimization problems. Every optimization algorithm suffers from certain limitations and GWFA is no exception to that. GWFA has a slower convergence speed due to the presence of multiple DAs and guidance directions. GWFA also uses a very simplistic local search for the DA exploitation procedure, however, it helps to keep the computational requirements to a minimum. In the future, we want to modify the model further by hybridization with some other meta-heuristic methods and apply the new model in other areas like hyperparameter optimization of neural networks, fast source estimation, feature selection, image contrast enhancement, etc.
RITAM GUHA received the undergraduate degree from the Department of Computer Science and Engineering, Jadavpur University, Kolkata, India, in 2020. He is currently pursuing the Ph.D. degree with Michigan State University, MI, USA. He is currently a Research Assistant with Hemlock Semiconductor Corporation, Hemlock, MI, USA. His research interests include intersection of evolutionary computation, AutoML, and interpretable AI.
SOULIB GHOSH received the undergraduate degree from the Department of Computer Science and Engineering, Jadavpur University, Kolkata, India, in 2020. He is currently working with the Algorithm and Optimization Team under the Samsung Sensor Team. His primary work is to develop a camera sensor algorithm for mobile cameras. His main research interests include image processing and machine learning.
KUSHAL KANTI GHOSH received the undergraduate degree from the Department of Computer Science and Engineering, Jadavpur University, Kolkata, India, in 2020. He is currently a Software Engineer with Microsoft, India.