Vehicle routing problems based on Harris Hawks optimization

The vehicle routing problem (VRP) is one of the challenging problems in optimization and can be described as combinatorial optimization and NP-hard problem. Researchers have used many artificial intelligence techniques in order to try to solve this problem. Among these techniques, metaheuristic algorithms that can perform random search are the most promising because they can be used to find the right solution in the shortest possible time. Therefore, in this paper, the Harris hawks optimization (HHO) algorithm was used to attempt to solve the VRP. The algorithm was applied to 10 scenarios and the experimental results revealed that the HHO had a strong ability to check for and find the best route as compared to other metaheuristic algorithms, namely, simulated annealing and artificial bee colony optimization. The comparison was based on three criteria: minimum objective function obtained, minimum number of iterations required and satisfaction of capacity constraints. In all scenarios, the HHO showed clear superiority over the other methods.

The vehicle routing problem (VRP) is a catch-all term for a group of issues involving "vehicles" completing visits to "customers. " The basic practical difficulty of distributing products to geographically dispersed consumers utilizing a number of trucks operating out of a single goods depot gives rise to these problems [7]. The challenge is to route the cars in such a way that the overall distance traveled (or time spent, money incurred, etc.) is kept to a minimum [8].
The VRP, also described as "vehicle scheduling, " "vehicle dispatching, " or simply "delivery, " emerges frequently in real-world circumstances that are not directly connected to the delivery of commodities [9][10][11]. Collection of mail from mailboxes, pick-up of children by school buses, house-call tours by a doctor, preventive maintenance inspection tours, laundry delivery, and so on are all VRPs in which the "delivery" operation can be a collection, collection and/or delivery, or neither; and in which the "goods" and "vehicles" can take a variety of forms, some of which may not even be of a physical nature [12].
The key goal is to find the minimum-length route that starts and ends in a depot, taking into account a number of variants such as the time window, service time, and vehicle capacity among others [13]. The VRP is a combinatorial optimization and NP-hard problem [14,15] and several algorithms have been proposed in order to solve it, yet it remains an incredibly important challenge in the field of combinational optimization [16]. The only small-scale VRP technology yet discovers the ultimate global solution [17]. The expertise of stream VRP algorithms is divided into exact and heuristics algorithms [18].
Classical heuristics for the VRP can be split into constructive heuristics and improvement heuristics [19]. The origin heuristics forever, progress from a solution to a best one in its neighborhoods [20].Heuristic methods execute a comparatively limited exploration of the search space and typically require a short computation time to output solutions that have good fitness [21]. Classification heuristics for the VRP include the savings algorithm, route-first cluster-second, cluster-first route-second, and insertion heuristics. Two kinds of improvement algorithms are suitable for the VRP, namely, intra-route heuristics and inter-route heuristics [22].
Metaheuristics are also widely used in optimization problems including the VRP [23][24][25].These smart algorithms are suitable for solving hard optimization problems that require a large amount of time because they are able to reach optimal solutions in a fast time [26,27], which is largely due to their utilization of random search instead of full search when exploring the search space [28].
The major goal of developing metaheuristic algorithms is to address complicated optimization issues where previous optimization approaches have failed. These techniques are now widely regarded as some of the most practical approaches to solve a wide range of real-world issues. Metaheuristic could be used to solve any problem that can be expressed as a problem of optimization technique. Also, They are ease of implementation as well as they are efficiency and flexibility [29].
The cooperative behavior and pursuit manner of Harris' hawks in nature, known as surprise pounce, is the fundamental inspiration for HHO. Several hawks work together to attack on a victim from different directions in an attempt to catch it off guard. Based on the dynamic nature of events and the prey's escape behaviors, Harris hawks can display a variety of pursuit patterns. To design an optimization method, this work mathematically duplicates such dynamic patterns and behaviors [50,52,53]. As an outcome, the Harris hawks optimization (HHO) algorithm was used to solve the VRP in this paper.
The rest of this paper is organized as follows: In "Related works", the related works on VRP are described. Then, in "Harris Hawks optimization", the HHO algorithm is discussed in detail. This is followed by a demonstration of the proposed method of using HHO for VRP in "Proposed method: HHO for VRP". After that, in "Experimental results", the experimental results are presented and analyzed. Finally, in "Conclusion" some conclusions are made and some possible directions for future improvements are suggested.

Related works
In recent times, metaheuristic algorithms have become an important part of the methods that have been developed to solve the VRP, largely because these algorithms have the ability to find optimal solutions due to their robust local and global search capabilities. Several previous studies have used these algorithms to solve the VRP, including Szeto, et al. [54], who presented a heuristic ABC for the capacitated VRP (CVRP). To improve the original heuristic for solving the CVRP, an improved version of the ABC heuristic was developed. Computing tests revealed that the improved heuristic ABC is capable of providing significantly stronger solutions for the CVRP than the original ABC and that it also consumes less computing energy.
On the other hand, Wang, et al. [55] addressed the periodic VRP with time window and service selection issue by building a heuristic algorithm based on improved ant colony optimization (IACO) and simulated annealing (SA), which they named multiobjective simulate annealingant colony optimization (MOSA-ACO). Also using SA, Rabbouch et al. [56] proposed an empirical-type SA to solve the CVRP. Their approach works incrementally by exploiting the last component of the worst practicable option. The results showed that their method has high search precision, optimal quicker entry, and the ability to classify all optimal rates while increasing the SA algorithm's convergence speed. In another SA-based research study, Ilhan [57] proposed using the population-based SA algorithm to solve the CVRP. Three operators were used to build the process: swap, inclusion and return. In fact, the algorithm used has been evaluated on 63 occasions. The research findings demonstrated that the proposed method has the ability to evaluate the best pathways for 23 instances, along with the presence of the percentage error value and the run time.
Another type of VRP was addressed in [58], where the author suggested an optimization algorithm for a hybrid ant colony and applied it to the multi-depot VRP (MDVRP). The suggested algorithm is based on a mixture of probabilistic techniques that are based on the normal actions of ants when foraging for food as well as on specific techniques that are based on the computational concepts of annealing. The experimental results showed the proposed algorithm is superior to other methods as it is able to produce the smallest mean error. Also, in [59], the authors recommended a way to minimize the MDVRP in terms of both distance traveled and time consumed by using inferential GAs. In their method, the proposed algorithm assigns each depot to its nearest customer and allows each customer to draw one customer at a time before all customers are allocated to drive. In addition, the algorithm uses Dijkstra's shortest path to evaluate each depot's closest customer. The findings showed that the suggested algorithm was calculated from distance and time of travel and thus confirmed its efficacy.
On the other hand, Sethanan and Jamrus [60] suggested the use of both an integer linear programming formulation and a novel hybrid differential evolution algorithm involving a fuzzy logic controller genetic operator to solve the backhaul and heterogeneous fleet routing problem for multi-trip vehicles. In their work, the objective was to minimize the cumulative cost associated with the distance traveled by using the traditional method of differential evolution and differential evolution with a selected genetic operator and real-settings fuzzy logic controller.
Meanwhile, in [61], the authors suggested a mathematical model to deal with the issue of the electric VRP (EVRP), taking into account the demands of electrical charging in the operations. To solve this optimization problem, the authors built a dedicated tabu search (TS) algorithm. The results of the improvement the VRP using this approach has obtained more reliable results than other methods. In another hybrid approach, Akbar and Aurachmana [62] also used a TS algorithm, in this case with a GA, in order to refine the CVRP with time windows (CVRPTW). One of the shortcomings of their work is its lack of animation by which to visualize the outcome and thereby encourage future studies. Nevertheless, their findings indicated that the proposed hybrid algorithm not only successfully minimizes the current route, but also estimates the optimal number of homogeneous fleets.
Based on the foregoing review of related works, it can be seen that several different metaheuristic algorithms have been used to solve the VRP and its variants. The results obtained by previous research studies showed that these smart algorithms can find the best routes for vehicles to reach customers in the least possible time. The success of these algorithms is due to the precision of their random search mechanism. and to the balance they achieve between the local and global search processes. It therefore seemed logical to investigate whether HHO would be able to find optimal solutions for the VRP.

Harris Hawks optimization
The HHO algorithm is a new swarm intelligence optimization algorithm that was proposed by Heidari et al. in 2019 [50]. The algorithm has been shown to perform efficiently in the optimization domain relative to other metaheuristic algorithms. Moreover, the algorithm can be extended to solve multiple types of optimization problems successfully, exhibiting a good level of performance [63]. Centered on the effects of such an algorithm in problems of intricacy optimization, the study uses it to increase the efficiency related to the VRP [50].
The basic notion behind the improved algorithm is derived from how Harris hawks hunt for food under desert conditions which are characterized by a scarcity of food [64]. To find their prey, Harris hawks engage in a variety of complex processes of mutual hunting and coordination to find, encircle, flush out and ultimately strike a possible target or prey. The HHO algorithm mimics the behavior of these hawks that monitor and catch their prey by means of a "surprise pounce", which is also known as the "seven kills" strategy [50].In this strategy, many hawks collaboratively target a single prey, such as a rabbit that is out in the open rather than hiding under cover of vegetation for example, by using different techniques and they ultimately converge on this prey [65].
The attack can be terminated promptly by catching the shocked target in just a few seconds. However, depending on the prey's escape capability and habits, the hawks may need to undertake several, short-length, swift dives that are likely to be in the vicinity of the prey for several minutes [51] before they are successful in executing the seven kills strategy. Harris hawks display a variety of types of chasing, which depends not only on the prey's escape habits, but also on a diverse range of environmental circumstances. For instance, in a swooping chase scenario, the most successful hawk (leader) swoops at the prey to capture it, but it can miss the prey (i.e., not find the target solution) because in essence the search and attack is carried out by a single member of the group. By increasing its vulnerability, Harris hawks proceed to locate an exhausted rabbit. In the end, when only one hawk, usually the most seasoned and strongest, can easily catch the exhausted rabbit and share it with the extant members, it cannot run away from the team besiege [66].
The HHO algorithm consists of three phases: (1) discovery or exploration (global search), (2) transformation from exploration to exploitation, and (3) exploitation (local search) [64]. For the intention of finding and investigating all aspects relevant to these propositions, the discovery phase reflects the change to many locations and a new venue. In addition, this phase illustrates a new technique for studying. The transformation from discovery to exploitation phase is regarded as crucial to the success of metaheuristic algorithms. In HHO, in order to transform the specified phases, the escape energy of the rabbit prey is used. The utilization process, on the other hand, reflects the installation and use of such sites and enables the available resources to be completely implemented [67].

Exploration phase
In this algorithm, the Harris hawks perch in certain positions frequently and wait to spot a target based on two techniques that are formulated in Eq. (1): where X(t + 1) is the hawks' location vector in the next iteration t, X rabbit ( t) is the rabbit's position, X(t) is the hawks' current location vector, r 1 , r 2 , r 3 , r 4 and q are random numbers within (0,1) that are modified in each iteration, LB and UB represent the upper and lower variable borders, X rand(t) is a hawk randomly selected from the current population, and X m(t) is the mean variable position. The average position of the hawks is reached by using Eq. (2).
where the position of each hawk in iteration t is indicated by X i (t) and N refers to the total number of hawks.

Transition from exploration to exploitation
In order to represent this step, the rabbit's energy is formulated as in Eq. (3): where E denotes the fleeing energy of the rabbit, T is the cumulative number of iterations, and E 0 is the rabbit's initial energy state, and t is the current iteration.

Exploitation phase
In the final exploitation phase, the hawks employ four different types of besiege technique to capture their prey.

Soft besiege
The following rules in Eqs. (4) and (5) describe the soft besiege behavior of the Harris hawks: where X(t) is the variance between the rabbit's place vector and the specific position in iteration t during the escape process, r 5 is a random number within (0, 1), and J = 2(1 − r 5 ) defines the rabbit's random leap ability. To mimic the nature of the rabbit's movements, the J value varies arbitrarily with each iteration.

Hard besiege
For a hard or intense besiege, the existing positions of the hawks are adjusted by Eq. (6):

Soft besiege with progressive rapid dives
When a soft besiege is accompanied by rapid dives, it is assumed that the hawks decide on adopting this technique based on the following rule in Eq. (7): It is also assumed that they use the following rule in Eq. (8) to dive at the prey based on Lévy flight-based patterns: where D is the issue aspect, S is a random vector of size 1 multiplied by D, and LF is the Lévy flight function that is determined using Eq. (9): where µ and σ are random values within (0, 1), and β is the default constant that is fixed to 1.5.
The final technique that is used to upgrade the positions of the hawks that are in the process of undertaking a soft besiege with dives can then be carried out by Eq. (10): where Y and Z are obtained by using Eqs. (7) and (8).

Hard besiege with progressive rapid dives
The final technique involves the combination of a hard besiege and rapid dives, as formulated in Eq. (11): where Y and Z are obtained by using the new rules in Eqs. (12) and (13): All of these steps of HHO are described in the following algorithm pseudocode [50]:

Proposed method: HHO for VRP
In this paper, we suggest using HHO to construct a routing strategy by focusing on minimizing the number of vehicles in two steps:

Stage 1
The HHO global search functionality is used in the first step, with the HHO fitness function defined as: where |σ| is the number of routes in the σ routing strategy and c(σi) is the travel expense of the σi route ϵ σ.
In the first iteration, the rabbit randomly chooses a direction of travel. A fitness function is calculated when it hits its destination (e.g., the transport distance) and the artificial pheromones are evaluated and identified by the network. These artificial pheromones affect the next hawk's preference. In the subsequent iterations, the rabbits appear to take the shorter direction. In addition, the chance of better pathways being chosen can be improved by artificial pheromones.
The decrease in the number of vehicles is typically the complication, instead of the primary feature of the searching process, of reducing transport costs. In VRP-SPDTW, because reducing the number of vehicles and travel costs are both important targets, we concentrate on one of them in each process.

Stage 2
In the second stage, the local search of HHO minimizes the travel cost by using the routing plans that were generated in Stage 1 as the initial solution. The search in this stage continues until a certain stopping condition is met. Figure 1 shows the mechanism of the suggested method.
There are three conditions that must be taken into account when solving the routing problem: (1) the vehicle starts from the warehouse and returns to it, (2) the vehicle visits one customer only in one journey and then returns to the place of departure i.e., the warehouse) before going to another customer and (3) the customer may not request a delivery from more than one vehicle.
Furthermore, each vehicle must not exceed a specific load (capacity), which is defined as follows: where N is the number of nodes (customers) and K is the number of vehicles.
As for addressing the issue of time complexity, the HHO's computational complexity depends primarily on three processes: initialization, fitness assessment, and hawk update [43].

Experimental results
To provide a fair scientific examination, the experiments used the same work settings and situations. The proposed method was performed on a Notebook that had a Windows 10 operating system, An Intel ® CoreTM i7-6006U processor operating at 2.00 GHz Since the input variables play an essential role in the consistency of the solutions. The input parameters were set in the experiments due to the obvious results as in [65]. In order for the outcomes of the experiment to be identical, the configurations of the algorithms are equivalent. The final set of parameters is specified in Table 1.

Application of HHO to 10 VRP scenarios
First, the HHO algorithm was applied to solve the VRP using 10 scenarios, where each scenario consisted of different numbers of customers and vehicles. In the first scenario, there were eight customers and three vehicles. These numbers were increased in each subsequent scenario to 70 customers and eight vehicles in the last scenario. This was done to make the problem more complex in order to measure the efficiency of the proposed approach using HHO in reaching the optimal solution for different sized problems. To assess the performance of HHO in each scenario, three criteria were calculated: (1) the minimum objective function obtained, (2) the minimum number of iterations required, and (3) the satisfaction of capacity constraints as explained in "Stage 2" Table 2 shows the results obtained by HHO for the 10 scenarios.
It can be seen from Table 2 that, in all 10 scenarios, HHO succeeded in reaching the optimum solution with a relatively low number of iterations. In addition, HHO was able to satisfy the capacity constraints in all scenarios. Figure 2 depicts the routes taken by the vehicles in order to reach the customers according to the conditions of the VRP for each scenario. In the figure, the route of each vehicle is depicted using a different color and is the route used to reach each customer at the lowest possible cost.
As seen in Fig. 2, in all ten simulations. HHO was able to find the best solution with a very small number of iterations. Furthermore, in all instances, HHO was able to meet the capacity limits. As the number of customers increases, it takes a lot of roads that must be followed in order to reach all customers as soon as possible; therefore, the number of roads in the first simulation is different from what it is in the tenth.

Comparison of HHO against SA and ABC
In order to evaluate the efficacy of the proposed method, the outcomes of the proposed method were compared with those of two methods in the literature, namely, SA [48]   and ABC [45]. SA and ABC are effective and largely used in many recent works in optimization field [68][69][70][71][72]. On the other hand, the comparison was made with two types of algorithms; single-based algorithm represented by SA, and population-based algorithm represented by ABC. The comparison with two different types of research method is very important to judge the proposed approach. The comparison was made by using the same parameter settings and criteria outlined above. Table 3 shows the outcomes of all three methods. As can be seen from Table 3, in sim.1 with eight customers and three vehicles, all three methods obtained the same minimum objective function of 220.1634 and satisfied the capacity constraints. However, HHO did so more quickly in just 13 iterations. However, in sim. 2 with 10 customers and three cars, SA performed better than HHO and ABC with an objective function of 269.8302. However, again the least number of iterations was achieved by HHO with 28 iterations. All three methods met the capacity constraints of sim.2.  Importantly, in sim.10 that had the highest level of complexity with 70 customers and eight vehicles, the lowest objective function was also achieved by HHO with 423.0535. The proposed method also performed the best in terms of convergence speed with 186 iterations.
Furthermore, a T-test was performed to compare the effectiveness of the HHO technique. The provided methodologies are utilized to calculate the findings statistics, which are based on the simulation's minimum objective function. Table 4 shows the results of a T-test on p-values and minimum objective function by the HHO, SA, and ABC, along with a 95 percent confidence interval (alpha value = 0.05).
As shown in Table 4,the HHO is more efficient than the SA and ABC, where the total number of p-values used in all simulations is less than 0.0001.As well as, the HHO has the least standard deviation and standard error rate in all simulations expect simulation 5. These findings suggest that the HHO is effective in resolving VRP.
In sum, in all ten scenarios, all three methods, HHO, SA and ABC, were able to satisfy the capacity constraints. However, HHO was clearly superior to the other two methods in terms of the minimum objective function obtained and the minimum number of iterations required. This indicates that HHO has the ability to achieve a balance between local and global search to reach more accurate solutions in a fast time. Table 5 provides more details on the number of iterations required in each simulation to reach the optimal solution compared to the other two methods.
From Table 5, the strength of HHO in every simulation is readily apparent. The use of HHO is advantageous whether the size of the problem is minimal, moderate or large, as it is able to find an optimal solution in the least number of iterations as compared to the other compared methods. Harris hawk optimization is not only able to find an optimal solution, it can do so quickly. This is due to its architecture, which achieves a better balance between exploration and exploitation, and it can transition smoothly between these two types of search. The number of iterations, i.e., convergence speed, is an important factor in solving complex optimization problems such as the VRP in a timely and cost-effective manner. Hence these results show that HHO has great potential for solving the VRP.

Conclusion
The VRP is an NP hard real-world problem, which many researches are still trying to provide the best solutions to solve it. In this paper, we proposed a model using HHO algorithm to solve the VRP. The performance of HHO was assessed according to three criteria: the minimum objective function obtained, the minimum number of iterations required, and the satisfaction of capacity constraints. Compared with two well-known methods, namely, SA and ABC, the results showed that HHO was clearly superior in terms of both the minimum objective function obtained and the minimum number of iterations required. This is likely due to the advanced exploitative capabilities of HHO as well as its ability to achieve better integration between the local and global search processes.
On the other hand, HHO like other metaheuristic algorithms, cannot prove optimality and it cannot probably reduce the search space, as well as it is hard to guarantee the consistency of optimization outcomes acquired with the same initial condition settings. In the future, to achieve better results, HHO could be employed in a hybrid model with another metaheuristic algorithm, or its global search mechanism could be modified through the use of several methods such as a crossover operator.

Acknowledgements
This work is supported by the Gulf University for Science and Technology, Kuwait.

Author contributions
We attest that all authors contributed significantly to the creation of this manuscript, each having fulfilled criteria. We confirm that the manuscript has been read and approved by all named authors. We confirm that the order of authors listed in the manuscript has been approved by all. All authors read and approved the final manuscript.

Funding
No funding was received for this work.

Availability of data and materials
The data that support the findings of this study are available on request.

Declarations Ethics approval and consent to participate
We confirm that any aspect of the work covered in this manuscript that has involved human patients has been conducted with the ethical approval of all relevant bodies and that such approvals are acknowledged within the manuscript.

Consent for publication
We give my consent for the publication of identifiable details, which can include photograph(s) and/or videos and/or case history and/or details within the text ("Material") to be published in the above Journal and Article. I confirm that I have seen and been given the opportunity to read both the Material and the Article to be published by Journal of big data.