1. INTRODUCTION
Estimation of Distribution Algorithms (EDAs) have been widely published in more than two decades. These algorithms have shown to be efficient to solve optimization problems. In addition, EDAs have been used to tackle combinatorial problems. It can be considered that there are two classifications for the EDAs, i.e., pure EDAs and hybrid EDAs.
Pure EDAs base their performance on a probability model to find new solutions from previous solutions. As pure EDAs are considered the Univariate Marginal Distribution Algorithm (UMDA) proposed by Mühlenbein (1997), the Mutual Information Maximization for Input Clustering (MIMIC) presented by De Bonet et al. (1997), the Combining Optimizers with Mutual Information Trees (COMIT) published by Baluja and Davies (1997), and the Bayesian Optimization Algorithm (BOA) proposed by Pelikan et al. (1998). These pure EDAs are some of the most frequently cited.
Hybrid EDAs base their performance on an interaction between a probability model and some another technique or method to generate offspring and/or to improve the performance of the algorithm. As hybrid EDAs is considered the Peña et al. (2004) research for solving synthetic optimizations problems. The Zhang et al. (2006) algorithm for the quadratic assignment problem. The Liu et al. (2011) study for the permutation flowshop scheduling problem. The Wang et al. (2012) method for the flexible jobshop scheduling problem. The Fang et al. (2015) algorithm for the stochastic resourceconstrained projectscheduling problem, and the Wang et al. (2016) research for the distributed permutation flowshop scheduling problems under machine breakdown. The articles mentioned above are current and representative of this group of hybrid EDAs.
Both pure EDAs and hybrid EDAs have been widely used to solve combinatorial problems. The development of EDAs does not finish in this point. Recently, new EDAs haven been published. These EDAs have a new characteristic. These utilize permutation of elementsbased representation as a solution to solve combinatorial optimization problems. That is, these algorithms propose specific probability models to integrate solutions as permutation based representation. This category is named distancebased ranking models. The proposed EDA by Ceberio et al. (2012) for flowshop scheduling problem. The study of PérezRodríguez et al. (2017) for the school bus routing problem with bus stop selection. The algorithm of PérezRodríguez and HernándezAguirre (2018) for the flexible jobshopscheduling problem with process plan flexibility, and the technique of PérezRodríguez and HernándezAguirre (2019) for the vehicle routing problem with time windows, are currently published papers belonging to this category.
From the previous classification, it can realize that recent authors have combined new and recognized methodologies to solve optimization problems such as the vehicle routing problem with time windows.
The development of new probability models continues being attractive for practitioners and researchers. The proposal of identifying new probability models, to be more efficient the performance of the EDAs, is the aim and scope for any research in this field. Therefore, as contribution of this research, the development of a new probability model is presented. The proposed EDA generates new solutions based on a radial probability model. The objective of this research is to identify new ways to improve the performance of the EDAs through application and use of new probability distributions.
It means that an atomic orbital generation function is used to produce offspring. In quantum chemistry it is established that a function, of this type, makes it possible to describe the behavior of an electron in a space occupied by an atom to which the electron belongs. However, due to the random behavior of the electron, it is more useful to describe its behavior in terms of the probability to find the electron in a specific volume of the space occupied by the atom to which the electron belongs.
In the current atom model, the electron is described in terms of a wave function ψ . The wave function is a math function that it describes the behavior of the electron in a specific space. This space is called atomic orbital.
The function ψ^{2} is proportional to the probability density of the electron in a specific point of the atomic space. If the values of ψ^{2} around of the core are considered, then it can define a contour surface (see Figure 1). This contour surface contains the volume where there exists a probability to find the electron. It permits to visualize the atomic orbital. In this research, the wave function of the hydrogen (H) is used. This function is elected because it is exactly defined in math terms. In this sense, there exists four orbitals for the hydrogen. Then, four equations describe the radial probability function.
The contribution of this research is to use such radial probability functions as probability model for the EDA scheme, and from these functions to generate offspring to solve the VRPTW. The performance of the proposed EDA, called RHEDA (Radial Hybrid Estimation of Distribution Algorithm), is compared with those recent algorithms that efficiently solve the VRPTW.
Although diverse methods and strategies have been used to solve the VRPTW, this paper contributes to the state of the art through utilization the radial probability functions of the hydrogen, for the VRPTW, as a probability model to enhance the performance of the EDA.
The results, obtained from this research, permits to conclude that using radial probability distributions is an emerging field to develop new and efficient EDAs.
2. RELATED WORK
The most studies make comparisons with other algorithms, to show the effectiveness of proposed approach. The main evolutionary algorithms used are the genetic algorithm, the particle swarm optimization algorithm, multiobjective algorithms, between others. The most representative studies such as Theophilus et al. (2021), FathollahiFard et al. (2019), Mojtahedi et al. (2021), Ruiz et al. (2019), SalavatiKhoshghalb et al. (2019), Pasha et al. (2020), Zhang et al. (2019), Trachanatzi et al. (2020), and FathollahiFard et al. (2020) are considered current research to tackle different routing problems.
Currently, other researchers show the effectiveness and efficiency of the proposed methods as a part of their manuscripts. The published methods try to get new solutions to some benchmark problems. The cited methods carried out experiments on actual examples, i.e., actual data of real environments. Some examples, for the VRPTW, of such methodologies are detailed below
Figliozzi (2010) offers an algorithm based on route construction and improvement. The algorithm mentioned can be applied sequentially to solve the vehicle routing problem (VRP) with Soft Time Windows (VRPSTW) and Hard Time Windows (VRPHTW). Computational results detail the performance of the proposed algorithm against other published solution methods using the Solomon (1987) benchmark problems.
Taş et al. (2013) study a vehicle routing problem with soft time windows and stochastic travel times. The authors propose a TS method to solve the problem mentioned. The Solomon (1987) instances were used to obtain competitive final solutions. In the first phase of the algorithm, an initial solution is constructed by Solomon’s (1987) a heuristic method. In the second phase, the initial solution is improved by a TS method based on research by Cordeau et al. (2001). In the third phase, a postoptimization method is applied to improve the solution obtained by the second phase.
Vidal et al. (2013) address efficiently a wide range of largescale VRP variants. The authors highlight that the most successful approaches to solving the VRPTW involve local search improvement procedures based on arcand nodeexchange neighborhoods. The authors indicate that the most competitive results are currently offered by the hybrid genetic algorithm (HGA) to solve VRP variants. The HGA combines methodologies such as powerful route minimization procedures, effective edge assembly cross operators, and extremely efficient local search procedures. These methods stand out in terms of simplicity and wider applicability.
PérezRodríguez and HernándezAguirre (2019) propose an estimation of distribution algorithm to solve the VRPTW. The algorithm uses the generalized Mallows distribution as a probability model to describe the distribution of the solution space.
PérezRodríguez (2020) details a hybrid approach that hybridizes an estimation of distribution algorithm and the mothflame algorithm to solve efficiently the VRPTW, and other combinatorial optimization problems.
Table 1 details the existing contributions in the related work.
3. PROBLEM STATEMENT
A complete formulation of the mathematical integer linear programming model, for the VRPTW, is develop by Toth and Vigo (2001). The aforementioned formulation is detailed below based on the seventh chapter of the Toth and Vigo’s (2001) book.
The depot is represented by the nodes 0 and n+1 . All feasible vehicle routes correspond to paths in G = (V, E) that start from node 0 and end at node n+1 . A time window is associated with nodes 0 and n+1 , i.e., $[{a}_{0},\hspace{0.17em}{b}_{0}]=[{a}_{n+1},\hspace{0.17em}{b}_{n+1}]=[E,\hspace{0.17em}L]$, where E and L represent the earliest possible departure from the depot and the latest possible arrival at the depot, respectively. Moreover, zero demands and service times are defined for these two nodes, that is, $d{e}_{0}=d{e}_{n+1}={s}_{0}={s}_{n+1}=0$. Two types of variables exist in this formulation: flow variables ${x}_{ijk}(i,j)\in A,k\in K,$ equal to 1 if arc (i, j) is used by vehicle k and 0 otherwise, and time variables ${w}_{ik},\hspace{0.17em}i\in V,$, k ∈ K, specifying the start of service at node i when serviced by vehicle k.
subject to
The objective function (1) expresses the total distance. Given $N=V\setminus \{0,n+1\}$ represents the set of customers. Constraints (2) restrict the assignment of each customer to exactly one vehicle route. Next, constraints (3)(5) characterize the flow on the path to be followed by vehicle k. Additionally, constraints (6)(9) guarantee schedule feasibility with respect to time considerations and capacity aspects, respectively. Note that for a given k, constraints (7) force w_{ik} = 0 whenever customer i is not visited by vehicle k. Finally, conditions (10) impose binary conditions on the flows variables.
An example is provided by PérezRodríguez and HernándezAguirre (2019). The data are given in Table 2.
Assuming a vehicle capacity Q = 50 , and establishing a specific route such (1, 4, 5, 2, 3, 6) , a feasible solution is depicted in Figure 2.
Based on the solution, depicted in Figure. 2, the distance travel is 1151.21. The goal is to find the best routing in order to obtain the minimum distance travel.
4. RHEDA FOR THE VRPTW
The RHEDA framework is provided below
4.1 Solution Representation
A solution for the VRPTW, in this paper, is a tour of all the customers. Each position, on the representation, contains a continuous value. Therefore, the representation is simply a sequence of continues values, called tour vector. The aforementioned continuous values, indicate the distance, in picometers, between the electron and the core. In Table 3, the representation of a tour vector is illustrated.
The representation, elected in this research, is suitable to integrate the radial probability distribution as a probability model. Anyway, the tour vectors must be decoded to represent valid routings. The decoding process is detailed as follows:
A fixed integer number is assigned for each customer. A sort on the continuous values of each tour vector is done. Assigning each continuous value to the corresponding fixed integer number that belongs to each customer. Table 4 details the previous tour vector, depicted in Table 3, and its decoding.
4.2 Fitness
The distance travel is computed for each member of the population. Using each decoding vector, a distance travel is obtained.
4.3 Probability Model
The math expressions of the radial distribution function P(r), for the hydrogen, in each atomic orbital are detailed as follows
1^{st} atomic orbital
where Z represents the atomic number of the element, i.e., the hydrogen Z =1 . a_{0} indicates the Bohr radius, and r is the distance (radius), in picometers (pm) of electron to the core.
The Bohr radius is defined as
${a}_{0}=\frac{{h}^{2}}{(4{\pi}^{2}me)}=52.9\text{pm}$, h is the Planck constant, m is the mass of the electron, and e its charge.
Figure 3 shows the radial distribution probability, for the 1^{st} atomic orbital, using the Eq. (11). As we can see, the function rapidly decays with respect to its distance to the core.
The next functions correspond to the 2^{nd}, 3^{rd}, and 4^{th} atomic orbitals for the hydrogen
With these radial distribution functions, a cumulative distribution should be built for each atomic orbital. The offspring can be generated using any cumulative radial distribution.
4.4 Sampling
The process to obtain an offspring is computed customer per customer. Firstly, a random value should be generated for each customer. Then, each random value, is interpolated in a cumulative probability distribution, previously selected, to identify which distance, between the electron and the core, should be established. Figure 4 shows an example of this process.
4.5 Replacement
The offspring should be evaluated to obtain their fitness. Finally, the replacement process used in this study is by binary tournament between the parents and the offspring.
All the stages of the proposed algorithm have been defined. All of this, within a number of the generations. In this research, 100 generations and 1000 solutions per generation, were used. These are fixed parameters.
5. RESULTS AND COMPARISON
5.1 Comparison between Atomic Orbitals
The standard benchmarking datasets, for the VRPTW, are used as input data for the mentioned comparison. The aforementioned datasets used in this research are the general and standard benchmarking dataset, i.e., Homberger and Gehring (2005) instances. These instances can be found at https://www.sintef.no/projectweb/top/vrptw/hom bergerbenchmark/200customers/
The relative percentage increase (RPI) is computed in order to compare the performance of each atomic orbital. The RPI is detailed as follows
where c_{i} is the distance travel obtained in the ith replication, and c_{*} is the best distance travel found and reported in the literature. The four atomic orbital functions, of the hydrogen, are compared in order to identify which atomic orbital is better to solve the VRPTW.
The experiments are executed in a Lanix Titan HX 4200 computer, Intel Core^{TM} i7 processor, 3.4 GHz, 8 GB of RAM, Windows 10 for 64 bits to run every instance. C++ language is used for the implementation for all the comparisons. To account for the stochastic nature of the RHEDA, we run 30 trials for each instance.
Figure 5 includes four box plots: the 1^{st} atomic orbital, the 2^{nd} atomic orbital, and so on, after running all the instances, and all the trials. As we can see, there is no practically difference between the orbitals, i.e., any orbital is suitable for solving the VRPTW instances used in the comparison.
5.2 Comparison between Recent Algorithms
In order to validate the scientific relevance of this paper, recent algorithms are proposed as a benchmark for comparison with the RHEDA scheme. The aforementioned recent algorithms are as follows

 The a cellular Genetic Algorithm (cGA) designed by Kamkar et al. (2010).

 The algorithmbased on route construction, and improvement detailed by Figliozzi (2010),

 The tabu search (TS) method presented by Taş et al. (2013),

 The hybrid genetic algorithm (HGA) proposed by Vidal et al. (2013),

 The hybrid estimation of distribution algorithm proposed by PérezRodríguez and Hernández Aguirre (2019), labeled as ‘Mallows’, and

 The hybrid approach that hybridizes an estimation of distribution algorithm and the mothflame algorithm designed by PérezRodríguez (2020), called HEDAMMF.
Three metrics are used to compare the performance of the algorithms. First, the mean absolute error (MAE)
where c_{i} is the fitness obtained in the ith replication, and c^{+} is the best fitnees found and reported in the literature for each instance used in this study.
The MAE is used to quantify the precision of the algorithm with respect to the values it should obtain
Second, the mean square error (MSE)
The MSE measures the amount of error between two data sets, that is, between the values that the algorithm returns and the values that it should obtain
Third, the relative percentage increase (RPI), Eq. (15). The RPI is used to compare two quantities while taking into account the “sizes” of the things being compared. The comparison is expressed as a ratio.
The experiments are executed in the same computer and language specification.
The Figure 6 indicates the output by the algorithms for the VRPTW using the Eq. (16). Through box and whisper charts, the dispersion of the values obtained, using the MAE metric, is appreciated. The dispersion of the values, using the MAE, shows the results over the instances and over different runs. As we can see, the algorithm labeled as ‘HEDAMMF’ obtain better results than the recent algorithms, as the RHEDA scheme. Based on the MAE metric, the RHEDA is more accurate with respect to the other algorithms.
Figure 7 details the Dunnett test; there is a statistically significant difference between all the recent algorithms and the RHEDA scheme. The RHEDA scheme outperforms all the recent algorithms for the VRPTW, and it has practically the same performance as HEDAMMF.
The Figure 8, shows the output by the algorithms for the VRPTW using the Eq. (17). The dispersion of the results is similar to the previous Figure 6. Based on the MSE metric, the RHEDA obtains the interval with the smallest error with respect to the other algorithms for the VRPTW.
Figure 9 details the Dunnett test; there is a statistically significant difference between all the recent algorithms and the RHEDA scheme. The RHEDA scheme outperforms all the recent algorithms for the VRPTW, and it has practically the same performance as HEDAMMF.
The Figure10 depicts the results obtained by the algorithms for the VRPTW using the Eq. (15). The algorithm labeled as ‘HEDAMMF’ obtain better results than the recent algorithms, and the RHEDA scheme outperforms all the previous results.
Figure 11 details the Dunnett test; there is a statistically significant difference between all the recent algorithms and the RHEDA scheme. The RHEDA scheme outperforms all the recent algorithms for the VRPTW, and it has practically the same performance as HEDAMMF.
As we can see, by using three different metrics, the performance between the RHEDA and the HEDAMMF is practically the same. The RHEDA outperforms the recent algorithms for the VRPTW. We can conclude that the radial distribution of the hydrogen is suitable and competitive against hybrid procedures to find the best routings for the VRPTW. The RHEDA does not need to be hybridized to find the best solutions for the VRPTW.
5.3 Parameters Setting
The EDA scheme considers population size, replacement (also known as generation gap), and selection strategy as key parameters. It is consistently with Grefenstette’s research (1986).

 The population size; in the current experiments, the population size ranged from 500 to 1000 solutions in increments of 500.

 The replacement; the current experiments allowed to vary the percentage of the population to be replaced during each generation between 50% and 100%, in increments of 50%.

 The selection; the experiments compared two ways to bubble, i.e., sorting by the distance travel and sorting by the number of vehicles.
A design of experiment is built to identify the best parameter of each parameter. Then parameter tuning is detailed below
Finally, the results of the parameter tuning is shown in Figure 12. There is no statistically significant difference of any of the three controlled parameters (number of generations, initial population size, and selected population size). Therefore, the parameters used are the same for all the algorithms.
Furthermore, the convergence pattern for the R1_2_7 instance, shown in Figure 13, between the RHEDA and the other algorithms is so closed. The performance of the RHEDA scheme is competitive.
Finally, the Figure 14 details the best solution, found by the RHEDA, for the RC2_2_8 instance.
6.CONCLUSIONS
This paper discusses the vehicle routing problem with time windows, which considers different routings, for each vehicle, as many logistics environments. To solve this problem, the RHEDA scheme is proposed. By means of wide and diverse numerical experiments and comparisons, the RHEDA offers a competitive performance. In addition, the RHEDA approach generates competitive solutions. Implementing these solutions can improve customer service by delivering orders on time and by avoiding delays in other process stages.
This research concludes that radial probability functions can be coupled with the EDA scheme in order to solve combinatorial optimization problems, such as VRPTW.
The computational results show that the different radial probability distributions, used for the VRPTW, with large data sets are suitable.
The results, obtained from this research, permits to conclude that using radial probability distributions is an emerging field to develop new and efficient EDAs.
A list of items as future research needs is detailed below,

•More radial probability distributions, from others elements not only the hydrogen, for the VRPTW and others combinatorial problems, should be considered in future research.

•Dynamic logistics issues should be included, in light of these results, for building new evolutionary algorithms coupled with other radial probability functions.

•Furthermore, the results of this research suggest testing the RHEDA on other related dynamic issues (traffic, heterogeneous fleet) of the routing systems.

•Since the RHEDA presents stability, it appears very suitable for implementation in software systems for practical purposes.

•Further research directions may deal with an extension of the RHEDA for building effective modules for specific users in the transportation sector. Finally, the RHEDA could be used in other types of routing systems.
We conclude that the VRPTW by means of a radial probability model has not been sufficiently considered in the literature.