Chemical reaction inspired approach for routing problems with hard time constraints

ABSTRACT

One of the most common issues in logistics and transportation planning is the vehicle routing problem. The objective is to find minimum cost-vehicle routes for serving a set of dispersed customers with deterministic demands while satisfying some constraints. In this paper, we focus on two main constraints, which are vehicle capacity and time windows for customers. We propose to tackle this NP-hard problem with a bio-inspired approach based on a hybridization of Chemical Reaction Optimization and local search method. The solving process of our approach mimics the chemical reactions between molecules to attain their stability with lowest energy. For more diversification of the search space, we involved in the initialization stage a local search method based on greedy adaptive procedure. The performances of our hybrid approach are assessed on the base of various benchmark instances. Simulation results show the good performances and the efficiency of the proposed approach.

KEYWORDS:

• Logistics and transportation
• routing problem
• bio-inspired approach
• greedy adaptive procedure
• chemical reaction optimization

1. Introduction

The vehicle routing problem (VRP) is the generic model of a whole class of combinatorial optimization problems that has theoretical and practical relevance to scientific communities and industries, especially logistics, transportation and telecommunications. In this model, a set of routes have to be determined for a fleet of vehicles based at one or several depots in order to deliver geographically scattered customers with known demands. This set involves minimum-cost vehicle routes originating and terminating at a depot. Determining the optimal solution is an NP-hard problem in combinatorial optimization. Hence, the size of instances that can be optimally solved is restricted. Therefore, industrial operators tend to use heuristics due to the size of real-world problems and the frequency of their emergence.

The VRP has many obvious applications in industry mainly in the fields of transportation, distribution, and logistics. Each real application has its specific constraints and characteristics. Consequently, several variations and specializations of the vehicle routing problem exist as stated in Han & Wang (2018). The problem was classified by (Irnich et al., 2014) according to:

• Network structure: node routing, arc routing.

• Type of transportation requests: delivery and collection, simple visits and vehicle scheduling, point-to-point transportation, repeated supply, non-split and split services, combined shipment and multimodal services, dynamic and stochastic routing.

• Route constraints: vehicle capacity, route length, multiple use of vehicles, time windows and scheduling aspects.

• Inter-route constraints: balancing constraints, inter-route resource constraints, synchronization constraints (task synchronization, operation synchronization, movement synchronization, load synchronization, resource synchronization).

• Fleet composition and location: homogeneous/heterogeneous or mixed fleet VRP

• Optimization objectives: single optimization objective, hierarchical objectives, multi-criteria optimization

The most studied variant is the basic version of VRP, the Capacitated Vehicle Routing Problem (CVRP). In the CVRP, the transportation request consists of a distribution from a single depot to a set of customers with deterministic demand. The vehicle fleet is initially available at the depot and assumed homogeneous meaning that all vehicles have the same capacity and operate at similar costs. The objective is to find minimum-cost routes to cover all customers while satisfying the capacity constraint for each route. We also notice that the vehicle routes must start and terminate at the depot and that no customer has visited more than once. In Toth & Vigo (2014), an overview of early and recent solving methods for the CVRP is presented. This overview covers classical exact methods (Poggi & Uchoa, 2014; Semet et al., 2014) as well as heuristics approaches (Laporte et al., 2014).

The second important variant of the VRP is the Vehicle Routing Problem with Time Windows (VRPTW). In addition to vehicle capacity constraint, each delivery location has a time window that must be respected. Each time window defines the early and late date of service. When time windows have to be rigorously respected, the problem is known as the VRP with Hard Time Windows (VRPHTW). However, if these windows may be violated by adding penalties in the objective function, the problem is called VRP with Soft Time Windows (VRPSTW). In literature, different approaches have been proposed to solve the two variants.

In this paper, our interest is focused on the VRP with Hard Time Windows (VRPHTW) for which we propose a new hybrid stochastic approach based on two metaheuristics: Greedy Randomized Adaptive Search Procedure (GRASP) and Chemical Reaction Optimization (CRO). The First method is based on local search and has been successfully applied to a big number of classical and real-world combinatorial optimization problems in different fields of science, networking, business and technology (Resende & Ribeiro, 2008). CRO is a nature-inspired and a population-based metaheuristic that was first proposed in Lam & Li, 2010; Siddique & Adeli (2017) where its applicability to solve some NP-hard problems in operation research is demonstrated. Among these problems, we mention Quadratic Assignment Problem, Resource-Constrained Project Scheduling Problem and a practical channel assignment problem in wireless networks. Then, CRO has been successfully applied to tackle diverse real-world problems in both discrete and continuous domains such as Task Scheduling in Grid Computing (Xu et al., 2011), continuous optimization problems with a large set of standard continuous benchmark functions (Lam et al., 2012), cloud job scheduling (Zain & Yousif, 2020), and recently for workflow scheduling in hybrid cloud-fog environment (Ramesh et al., 2024). These applications have shown the splendid performance of CRO method. Besides its good performance in solving NP-hard problems, and to the best of our knowledge CRO method has not been adopted before for tackling VRPs.

Therefore, we chose to perform the optimization process according to CRO scheme. Moreover, to improve the exploration of the search space through intensification and diversification strategies, we embedded in our approach to the local search based method GRASP.

We notice that this paper extends our works in Boudali & Ragmoun (2022). We focus on the proposed diversification and intensification strategies during the exploration of the search space. The involved bio-inspired operators are detailed and illustrated through examples to explain how the exploration of promising regions is intensified and how the local optima are escaped. We also emphasize on parameter configuration step through a set of tests. Furthermore, we provide a comparative study of the best-obtained solutions with the proposed approach and the best-known ones in literature.

The remainder of this paper is organized as follows. Section 2 presents a mathematical formulation of the routing problem with hard time windows. Then, Section 3 & 4 exposes related works and discusses the most relevant existing approaches in order to justify our contribution. Section 5 describes the CRO process. Afterwards, we focus on our hybrid approach in Section 6. We present in Section 7, the performances of the proposed approach through simulations based on benchmark problems. We conclude this paper and give some potential future works in Section 8.

2. Mathematical formulation

In this section, we provide a mathematical formulation of the VRPHTW (Toth & Vigo, 2014). We denote by:

• G (N, A): non oriented graph;

• N: set of nodes {0,.., n} including n customers and the depot 0;

• A: set of arc (i, j) between two customers i and j;

• V: set of m identical vehicles;

• Q: capacity of each vehicle;

• q_i: the demand of customer i;

• c_ij: the cost for arc (i, j);

• e_i, l_i: the earlier and latest date of service for customer i;

• a_i^v: the arriving date of vehicle v at customer i;

• s_i: service time of customer i;

We also consider the following decision variable: $x_{\mathrm{ij}}^v=\{\begin{array}{c}1\mathrm{if}\mathrm{arc}(i,j)\mathrm{is}\mathrm{travelled}\mathrm{by}\mathrm{vehicle}v\\ 0\mathrm{otherwise}\end{array}$The objective function corresponds to the minimization of total travelled distances by vehicles. Formally, it is defined by equation (1) (1) $\mathrm{Minimize}\sum _{v\in V}\sum _{i\in N}\sum _{j\in N}{c}_{\mathrm{ij}}{x}_{\mathrm{ij}}^v$(1) Subject to: (2) $\sum _{v\in V}\sum _{j\in N}{x}_{\mathrm{ij}}^v=1\mathrm{\forall i}\in N$(2) (3) $\sum _{i\in C}{q}_i\sum _{j\in N}{x}_{\mathrm{ij}}^v\le Q\mathrm{\forall v}\in V$(3) (4) $\sum _{j\in N}{x}_{0j}^v=1\mathrm{\forall v}\in V$(4) (5) $\sum _{i\in N}{x}_{\mathrm{ik}}^v-\sum _{j\in N}{x}_{\mathrm{kj}}^v=0\mathrm{\forall k}\in N,\mathrm{\forall v}\in V$(5) (6) $\sum _{i\in N}{x}_{i0}^v=1\mathrm{\forall v}\in V$(6) (7) $x_{\mathrm{ij}}^v(a_i^v+s_i+c_{\mathrm{ij}})\le a_j^v\mathrm{\forall }(i,j)\in A,\mathrm{\forall v}\in V$(7) (8) $e_i\le a_i^v\le l_i\mathrm{\forall }i\in N,\mathrm{\forall v}\in V$(8)

Equation (2) ensures that each customer is visited by exactly one vehicle, while Equation (3) guarantees the respect of vehicle capacity constraint. Each vehicle route must start and finish at the depot as stated in Equations (4), (5) and (6). Hard time windows constraints are defined by Equation (7) and Equation (8). In fact, a vehicle must wait when he arrived before e_i at customer i, but he is not accepted when he arrived after l_i.

3. Related works

In Desaulniers et al. (2014), authors present a literature review of exact methods, heuristics and metaheuristics that were proposed for the VRPTW. Among exact approaches, we find subset-row inequalities (Jepsen et al., 2008), generalized k-path inequality (Desaulniers et al., 2008) and route relaxation and pricing strategies (Baldacci et al., 2011). These strategies are based on set partitioning formulation of the problem. For further details on these works, we can refer to a review on exact methods for the VRPTW proposed by (Baldacci et al., 2012; Kallehauge, 2008). In the past decades, research on heuristic methods for the VRPTW has mainly focused on metaheuristics (Swan et al., 2022). In fact, a little interest has been conducted on heuristic techniques that are organized into three categories: constructive methods (saving methods, insertion methods), improvement methods (β-Opt) and two-phases methods (Cluster-first Route-second/Route-first, Cluster second) (Desaulniers et al., 2014). As constructive algorithm, we mention Clarke and Wright saving algorithm that has been recently adapted for a real case of VRPTW arising in urban delivery (Cao et al., 2020).

Furthermore, most of metaheuristics have been adapted to the case of VRPTW (Dixit et al., 2019): Genetic Algorithms (GA), Ant Colony Optimization (ACO), Tabu Search (TS), Artificial Bee Colony algorithm (ABC), Particle Swarm Optimization (PSO), etc. These metaheuristics have been successfully applied separately or in a combined way. Hybrid approaches involve the hill-climbing algorithm with the genetic algorithm, which effectively improved the speed of VRPTW (Fan & Feng, 2018). In addition, a multi-type ant system (MTAS) algorithm was developed based on the multi-ant colony algorithm with local search process (Deng et al., 2018). In (Osaba et al., 2017), authors proposed for VRPTW an evolutionary discrete bat algorithm (EDBA) by randomly inserting multiple heuristic operators. A hybridization of Tabu Search and Neural Networks was also implemented for VRPTW on Google Maps (Kirci, 2016). Other improved population-based approaches have been also developed for the VRPTW, such as: improved population-based incremental learning algorithm (IPBIL) (Xie et al., 2016), improved artificial bee colony algorithm IABC (Yao et al., 2017), dynamic hybrid ant colony optimization algorithm (DHACO) based on fusion strategy and cloud association rules (Ge et al., 2015). Another hybrid approach for the VRPTW, was proposed by (Saksuriya & Likasiri, 2022) based on Local Search (LS), Ruin and Recreate procedure (R&R) and Particle Swarm Optimization (PSO). A recent contribution was proposed by (Dornemann, 2023) for the capacitated VRPTW, by using a supervised deep learning approach. The proposed model employs a graph convolutional neural network that assists in finding solutions by providing a probability distribution to guide a tree search.

4. Discussion

As stated in (Toth & Vigo, 2014), exact approaches quickly showed their limitations in solving complex and large instances of vehicle routing problems. The emergence of specific heuristics has contributed to obtaining fast good solutions in a short amount of time, which is of high importance in real-world cases. However, the main disadvantages of heuristic techniques are the lack of flexibility (adding/changing decision variables, goals, constraints, etc.), and their inability to produce optimal solutions.

With the advent of metaheuristics during the last decades, we noticed huge progress in producing efficiently high-quality solutions for various instances with large sizes (Swan et al., 2022). Another benefit of metaheuristics are their flexibility to integrate new constraints, goals, decision variables, which is an important feature for solving real-world problems. Bio-inspired metaheuristics, like ant colony optimization, genetic algorithms, particle swarm, harmony search, etc. have been widely and successfully applied to the VRPTW in a separate or combined way as discussed above.

However, another important bio-inspired metaheuristic, CRO, has effectively contributed to solving various real-world optimization problems, such as Task Scheduling in Grid Computing (Xu et al., 2011), cloud job scheduling (Zain & Yousif, 2020), and recently used for workflow scheduling in hybrid cloud-fog environment (Ramesh et al., 2024). To the best of our knowledge, this method has not been explored in literature for the VRPHTW. The contribution of this work is to propose a new CRO-based approach for our problem by integrating a local search that promotes the diversification search. GRASP is selected for the local search as a commonly applied method in industrial field. This method provides efficient neighbourhood search techniques for solving various combinatorial optimization problems.

5. Nature inspired metaheuristic: CRO

Nature-inspired meta-heuristics have dominated the scientific literature of optimization in the last decades. In this context, CRO is a bio-inspired metaheuristic which is based on population of solutions. CRO has been designed on the basis of interactions between molecules in chemical reactions (Lam et al., 2013; Lam & Li, 2010; Lam & Li, 2012; Siddique & Adeli, 2017). During chemical reaction process, molecules continuously change in order to attain a stable state with the lowest free energy. Each molecule holds a molecular structure and two kinds of energies: potential energy (PE) and kinetic energy (KE). We notice that PE corresponds to solution fitness value while KE quantifies the tolerance of less stable structure. Thus, KE specifies the solution’s ability to escape local optima: large value of KE indicates that molecular structure could have a rigorous change.

Moreover, CRO relies on the conservation of energy law, which states that energy in a system is conserved and can neither be created nor destroyed. According to this law, energy is transmitted among molecules or within a molecule. To sustain kinetic energy of molecules, a central buffer is used. As the process progresses, molecules tend to have less KE.

For more explanation of the chemical reaction process, we consider a closed container with a set of molecules. During this process, molecules experience a series of elementary reactions that may alter their structures and related energies. We distinguish two sorts of elementary reactions: uni-molecular and inter-molecular collisions. In the first category, only one molecule is involved to undergo decomposition or on-wall ineffective collision. While in the second category, more than one molecule are concerned with inter-molecular ineffective collision or synthesis.

In the two ineffective collisions, molecules are modified with new structures of their surroundings in the search space, which is known as the potential energy surface (PES). Thus, these two collisions correspond to intensification operators. Conversely, synthesis and decomposition result in new molecules with different molecular structures. Consequently, these operators provide the ability to explore new regions of the PES and so, a diversification of the search. The energies are then, redistributed (PE and KE) among the different molecules in the container. Through the sequence of reactions, different regions of the PES are explored in an intelligent way.

The CRO algorithm is described by three stages (Lam et al., 2012; Lam & Li, 2010), (Zain & Yousif, 2020): initialization, optimization and final stage. In the first stage, some control parameters are defined and the initial population of solutions is generated. Afterwards, an optimization stage is iteratively performed through a sequence of reactions until a stopping criterion is encountered. At each iteration, the algorithm first decides whether a uni-molecular or inter-molecular collision will be performed according to molecular collision parameter. Afterwards, the algorithm checks the decomposition criterion in the case of uni-molecular collision and assesses the synthesis criterion in the case of inter-molecular collision. Finally, the solution with the lowest energy is provided as the output of the process.

6. The proposed chemical reaction-inspired approach

The proposed approach for the vehicle routing problem with hard time windows is defined by a preliminary step during which a parameter setting and population initialization are performed. Then, an optimization process is launched according to a CRO scheme. The following subsections describe the proposed solution encoding, the preliminary step as well as the optimization process.

6.1. Solution encoding

A Molecule or problem solution is a set of feasible routes that begin and end at the depot. Therefore, the molecular structure of a solution M is defined by an ordered sequence of customer identifiers as illustrated in Figure 1. In the presented solution structure, two routes are defined. The limits of a route in this sequence are detected when capacity and/or time windows constraints are no longer satisfied. Notice that the depot identifier corresponds to 0.

Figure 1. Solution encoding.

6.2. Parameter initialization

The proposed nature-inspired method begins with a preliminary stage during which an initialization of CRO parameters is performed. The CRO parameters are described as follows (Boudali & Ragmoun, 2022):

• InitialKE: initial KE for each initial solution;

• Buffer: central energy buffer, initially empty with Buffer = 0;

• PopSize: initial size of solution population;

• KELoss: upper limit of KE loss in on-wall ineffective collisions.

• MoleColl: faction of all elementary reactions corresponding to inter-molecular reactions.

• α: specifies the decomposition criterion. It corresponds to maximum time that a selected molecule is allowed to stay in a stable state;

• β: specifies the synthesis criterion. It corresponds to the least kinetic energy a molecule should possess to undergo a synthesis.

Once these parameters are defined, it is crucial to generate an initial population with size PopSize. We notice that the population size varies during the optimization process due to decomposition and synthesis operators.

6.3. Population initialization

Since the performances of the optimization process highly depend on the quality of initial population, we embedded in our approach random procedure with GRASP for generating initial population. The greedy randomized search method has been widely and successfully applied to classical and real-world optimization problems in different areas of business, science and technology (Resende & Ribeiro, 2008). As shown in Figure 2, two main steps define the iterative scheme of the procedure: construction and local search.

Figure 2. Greedy randomized adaptive search procedure for population initialization.

The construction step generates a greedy randomized solution by progressively inserting elements to the solution set and preserving feasibility. Inserted elements are selected from a list of candidates known as Restricted Candidate List (RCL). In this list, solutions are ranked according to the quality of solution they will achieve. Given the obtained feasible solution, a local search is applied in its neighbourhood for improvement until a local minimum is encountered. Then, the best solution is retained as a result. At the end of this procedure, PopSize solutions are generated and memorized in Pop_init. Each solution (or molecule) is described by some properties:

• PE: the fitness in the objective function;

• KE: Kinetic energy initialized with InitialKE;

• NumHit: number of hits, initially 0;

• MinHit: number of hit when the molecule has the best solution;

• MinVal: best encountered value of PE.

6.4. Optimization mechanism

Once parameter setting is determined and initial population is generated, the optimization process is launched as illustrated in Figure 3. This process is iteratively performed until stopping criteria are satisfied (Maximum iterations without improvements, fixed CPU time, etc.). In order to explain the bio-inspired optimization process, we assume an initial set of molecules (or solutions) in a closed container. These molecules with excessive energy move and collide either on the container walls or on each other. Through these collisions, molecules undergo a series of elementary reactions that may change their molecular structures and energies until they reach a stable state with the lowest free energy. In fact, at each iteration of this process, a random value in [0,1] is generated to decide uni-molecular or inter-molecular collision. Accordingly, a selection of one or two molecule(s) is performed to undergo a specific elementary reaction (see Table 1): decomposition, on-wall ineffective collision, synthesis or inter-molecular ineffective collision. The fitness value of new molecule(s) is computed and the energy conservation law is checked before the update of the population. In the following section, we provide details about the properties and the different steps of each involved elementary reaction in the optimization process.

Figure 3. Flow diagram of the optimization process.

Table 1. Elementary reactions according to molecularity.

	Uni-molecular reactions	Inter-molecular reactions
Fixed population size	On-wall ineffective collision	Inter-molecular ineffective collision
Variable population size	Decomposition	Synthesis

6.5. Properties of elementary reactions

The different elementary reactions may be classified according to the number of involved molecules. As shown in Table 1, the population size is variable with the application of decomposition and synthesis throughout the optimization process. So, in one hand we find operators involving only one selected molecule to generate a new neighboring one or two new molecules. In the other hand, we find operators implying two molecules to result in a combined molecule or new neighbouring ones. The success of reactions through the checking of energy conservation law, decides the substitution of molecules in the population.

Moreover, as illustrated in Table 2, some reactions act as intensification operators since they produce neighbouring molecules. While other reactions are considered as diversification operators given their explorative capabilities of other regions of the search spaces by generating new molecular structures.

Table 2. Elementary reactions according to intensification and diversification.

Elementary reactions	Intensification	Diversification
On-wall ineffective collision	X
Decomposition		X
Inter-molecular ineffective collision	X
Synthesis		X

In the remainder of this section, we detail the different steps of each elementary reaction in our optimization process.

6.5.1. On-wall ineffective collision

In this reaction, the selected molecule M hits a wall and bounces back. This collision causes minor changes to the molecular structure of the molecule. Thus, we get a new molecule M’ from the initial one by performing a neighbourhood search operator such as pair-exchange or 2-opt as illustrated in Figure 4. The fitness value PE_M’ is calculated and the energy conservation law is verified: PE_M + KE_M ≥ PE_M’. After, KE_M’ is computed and the central energy Buffer is updated. Solution M is then updated with the new molecular structure of M’ and the corresponding energies PE_M’ and KE_M’. The pseudo-code of this elementary reaction is illustrated in Algorithm 1.

Figure 4. Illustration of on-wall ineffective collision.

Table

Algorithm 1. On-wall ineffective collision
Input: M, Buffer, objective function f, KELoss
Output: M, Buffer
generate M’ = Neighbourhood(M)
2. PEM = f(M)
3. If PEM + KEM≥PEM’ Then
4.   generate q in [KELoss, 1]
5.   KEw'= (PEw + KEw – PEw’) × q
6.   Buffer = Buffer + (PEw + KEw – PEw’) × (1 – q)
7.   M = M’

6.5.2. Inter-molecular ineffective collision

In this collision, two selected molecules M1, M2 hit each other and then bounce away. This reaction causes small changes in their structures. As shown in Figure 5, we get two new molecules M1’, M2’ from the existing ones by randomly swapping two values of their molecular structures. We can apply the same operator as on-wall ineffective collision to generate the new molecules from the existing ones. For each new molecule, PE is calculated and the energy conservation law is checked. Then, KE of the transformed molecules share the remaining energy in the sub-system and the population is updated with the new molecules. This collision operator is detailed in Algorithm 2.

Figure 5. Illustration of inter-molecular ineffective collision.

Table

Algorithm 2. Inter-molecular ineffective collision
Input: M1, M2, objective function f
Output: M1, M2
M1’ =  Neighbourhood(M1)
2. M2’ = Neighbourhood (M2)
3. PEM1’ = f(M1’)
4. PEM2’ = f(M2’)
5. Diff = (PEM1 + PEM2 + PEM1’ + KEM1 + KEM2) – (PEM1’ – PEM2’)
6. If Diff > 0
7.   generate p in [0, 1]
8.   KEM1’ = Diff * p
9.   KEM2’ = Diff * (1 – p)
10   M1 = M1’
11.   M2 = M2’

6.5.3. Decomposition

In decomposition, a selected molecule M hits the wall, which causes a significant change to the molecular structure. As a result, the molecule M decomposes into two new molecules M1 and M2. In Figure 6, we illustrate an example of decomposition operator. The two new molecular structures are obtained from the original structure by randomly generating two indexes ind1, ind2 in [−n, +n], with n the size of a solution (i.e. the number of customers). These indexes specify the direction and the number of cyclic swaps to be applied on molecular structure of M to obtain M1 and M2. In the example of Figure 6, we assume ind1 = 2 and ind2 = −1. Thus, molecular structure of M1 is obtained from M by performing 4 cyclic permutations in the left-right direction. While molecular structure of M2 is obtained by 3 cyclic permutations in the opposite direction. The fitness value (PE) of each molecular structure is calculated and the conservation energy law is assessed. In case of success, the population is updated with the new molecules M1, M2 instead of M. The different steps of this elementary reaction are presented in Algorithm 3.

Figure 6. Illustration of decomposition.

Table

Algorithm 3. Decomposition
Input: M, Buffer, objective function f, Success
Ouput: M1, M2, Buffer
1. M1 = Neighbourhood(M)
2. M2 = Neighbourhood(M)
3. PEM1 = f(M1)
4. PEM2 = f(M2)
5. Diff = PEM + KEM – PEM1 – PEM2
6. If Diff ≥ 0
7.   Success = True
8.   generate k in [0, 1]
9.   KEM1 = Diff * k
10   KEM2 = Diff * (1 – k)
11. Else If Diff + Buffer ≥ 0
12.   Success = True
13.   generate m1, m2, m3 et m4 in [0, 1]
14.   KEM1 = (Diff + Buffer) * m1* m2
15.   KEM2 = (Diff + Buffer – KEM1) * m3* m4
16. Else
17.   Success = False

6.5.4. Synthesis

In synthesis, two randomly selected molecules collide with each other and fuse to generate a new molecule. As shown in Figure 7, two molecules M1 and M2 hit against each other to produce new molecule M’. In this example, two distinct values are randomly picked from the structure of M1. These values are placed at the same position in the new molecular structure of M’. Then, the remainder values are filled in M’ according to their appearance in structure of M2. Accordingly, each customer identifier appears only once in the molecular structure of M’. Afterwards, the energy conservation law is checked to decide the success of synthesis collision, as shown in Algorithm 4.

Figure 7. Illustration of synthesis.

Table

Algorithm 4. Synthesis
Input: M1, M2
Output: M’, Success
1. generate M’ from M1, M2
2. PEM’ = f (M’)
3. If PEM1 + PEM2 + KEM1 + KEM2 ≥ PEM’
4.   Success = True
5. Else
6.   Success = False

7. Evaluation and experimentation

In order to assess our hybrid approach for the VRPHTW, we implemented it on jMetal framework by using a Desktop Intel® Core(TM) i5-1135G7 CPU @ 2.4 GHz with 16.0 Go of RAM. Afterwards, we considered benchmark instances of different types and sizes. In this section, we first present the characteristics of these instances (number of customers, geographical dispersion, vehicle capacity, time windows, etc.). Then, we provide an empirical study to define the parameter configuration of our approach. Finally, we present a comparative study of our best results with the best-known solutions in literature based on benchmark instances (Solomon’s instances).

7.1. Benchmark instances

In our experimental study, we consider 18 instances from the predefined 56 Euclidien instances of Solomon (1987) as described in Table 3. These problems are characterized by one depot, a homogeneous fleet of vehicles, a number of customers up to 100 with geographic coordinates and spatio-temporal constraints. These instances are organized according to three categories:

• class C (Cluster): holds instances with geographically clustered customers;

• class R (Random): embraces instances with randomly dispersed customers;

• class (RC): holds instances with combined characteristics;

Table 3. Characteristics of considered benchmark instances.

Instances	Number of customers	Vehicle capacity	Distribution of customers	Time windows
R101	25, 50 and 100	200	Random	Narrow
R102	25, 50 and 100	200	Random	Narrow
R103	25, 50 and 100	200	Random	Narrow
R201	25, 50 and 100	1000	Random	Wide
R202	25, 50 and 100	1000	Random	Wide
R203	25, 50 and 100	1000	Random	Wide
C101	25, 50 and 100	200	Cluseterd	Narrow
C102	25, 50 and 100	200	Cluseterd	Narrow
C103	25, 50 and 100	200	Cluseterd	Narrow
C201	25, 50 and 100	700	Cluseterd	Wide
C202	25, 50 and 100	700	Cluseterd	Wide
C203	25, 50 and 100	700	Cluseterd	Wide
RC101	25, 50 and 100	200	Random/Cluseterd	Narrow
RC102	25, 50 and 100	200	Random/Cluseterd	Narrow
RC103	25, 50 and 100	200	Random/Cluseterd	Narrow
RC201	25, 50 and 100	1000	Random/Cluseterd	Wide
RC202	25, 50 and 100	1000	Random/Cluseterd	Wide
RC203	25, 50 and 100	1000	Random/Cluseterd	Wide

Moreover, in each class we find subclasses as follows:

• Instances of subclasses R1, C1 and RC1 have low vehicle capacity, reducing the number of clients per vehicle;

• Instances of subclasses R2, C2 and RC2 have high-capacity vehicles, allowing a large number of customers per vehicle.

Let us notice that the service time is 10 for classes R, RC and 90 for class C.

7.2 Simulation results

As the first step of experimental study, we performed some preliminary tests to determine the optimal parameter configuration of our approach. These parameters are: maximum number of iteration NI, the rate KELoss, MoleColl, InitialKE, α and β. Tables 4–6 illustrate parts of this preliminary study by varying the problem size 25, 50 and 100 customers, respectively, and by using parameter values inspired by literature of CRO. According to this tuning process, we noticed that the performance of our approach is getting better when the maximum number of iteration NI increases, the values of KELoss and α rise with a decreasing quantity of MoleColl and β. Following this study, we retained these parameters values: NI = 18,000, InitialKE = 10; KELoss = 0.9; MoleColl = 0.1; α =  100; β =  0.1. Once this preliminary step is finalized, the generation of initial population is performed according to random procedure and GRASP method. This population holds initial molecules that will undergo the presented optimization process.

Table 4. Part of parameter configuration for instances with 25 customers.

	NI = 1000 KELoss = 0.1 MoleColl = 0.7 α = 25 β = 0.5	NI = 15,000 KELoss = 0.5 MoleColl = 0.5 α = 50 β = 0.3	NI = 18,000 KELoss = 0.9 MoleColl = 0.1 α = 100 β = 0.1
R101	736.99	635.50	618.32
R102	682.08	573.69	549.35
R103	649.74	512.36	461.55
R201	666.86	493.52	465.18
R202	626.08	480.09	429.35
R203	607.10	441.07	425.51
C101	802.35	628.36	438.69
C102	667.98	524.56	434.30
C103	659.79	512.05	392.31
C201	794.12	574.24	457.53
C202	656.90	566.59	431.59
C203	655.37	528.11	387.33
RC101	1391.82	1126.56	984.65
RC102	1176.23	1058.01	880.66
RC103	1012.57	962.20	735.56
RC201	1197.17	1032.88	778.35
RC202	1073.30	872.61	708.11
RC203	954.23	862.35	643.88

Table 5. Part of parameter configuration for instances with 50 customers.

	NI = 1000 KELoss = 0.1 MoleColl = 0.7 α =  25 β =  0.5	NI = 15,000 KELoss = 0.5 MoleColl = 0.5 α =  50 β =  0.3	NI = 18,000 KELoss = 0.9 MoleColl = 0.1 α =  100 β =  0.1
R101	1347.89	1178.14	1068.42
R102	1181.89	1026.97	925.75
R103	1316.90	1053.95	791.90
R201	1285.25	1006.19	829.59
R202	1047.12	917.40	751.47
R203	883.08	729.01	662.14
C101	802.35	628.36	438.69
C102	667.98	524.56	434.30
C103	659.79	512.05	392.31
C201	794.12	574.24	457.53
C202	656.90	566.59	431.59
C203	655.37	528.11	387.33
RC101	1391.82	1126.56	984.65
RC102	1176.23	1058.01	880.66
RC103	1012.57	962.20	735.56
RC201	1197.17	1032.88	778.35
RC202	1073.30	872.61	708.11
RC203	954.23	862.35	643.88

Table 6. Part of parameter configuration for instances with 100 customers.

	NI = 1000 KELoss = 0.1 MoleColl = 0.7 α = 25 β = 0.5	NI = 15,000 KELoss = 0.5 MoleColl = 0.5 α = 50 β = 0.3	NI = 18,000 KELoss = 0.9 MoleColl = 0.1 α = 100 β = 0.1
R101	2292.05	2083.45	1721.11
R102	1980.61	1723.44	1545.36
R103	2416.67	1516.78	1268.39
R201	1944.10	1838.65	1296.00
R202	1976.17	1418.89	1191.39
R203	1713.91	1277.87	1027.29
C101	2164.96	1528.12	1197.03
C102	2045.35	1538.57	1112.34
C103	1938.25	1335.44	1029.15
C201	1982.80	1372.48	988.50
C202	1747.77	1116.16	910.23
C203	1785.02	1232.58	816.52
RC101	2723.29	2141.78	1827.18
RC102	2594.57	1887.86	1618.16
RC103	2686.07	1643.25	1423.28
RC201	2359.41	1878.72	1516.65
RC202	2148.55	1540.36	1335.08
RC203	2000.21	1365.97	1105.34

In Figures 8–10, we present the obtained results with the adopted configuration by considering instances with 25, 50 and 100 customers, respectively. In each figure, we illustrate the best-found solution by the proposed Chemical Reaction-inspired Approach (GRASP-CRO) and the best known solution from literature. We notice that for each instance with specific size, we performed 10 runs of the algorithm and we retained the best found solution.

Figure 8. Simulation results for instances with 25 customers.

Figure 9. Simulation results for instances with 50 customers.

Figure 10. Simulation results for instances with 100 customers.

According to Figure 8, the shape of the curves are very similar. In fact, the obtained results with the hybrid approach are very close to the best-known solutions for most instances with 25 instances, especially for instances of class C and RC. Thus, our approach can reach optimality for small instances with up to 25 customers that may be geographically clustered or randomly dispersed.

In Figure 9, we notice similar shapes of curves with a decreasing deviation from instances of class R to class RC. The obtained results for the three classes with 50 customers are similarly promising and even competitive for problems of class RC.

In Figure 10, we show simulation results for large problems with 100 customers. As illustrated, the deviation of the obtained solutions from the best-known ones is very negligible for instances of class R and RC.

As we considered the best-known solutions from literature as reference points in the search space, we chose the percent deviation measure to assess the quality of our best found solutions with the hybrid approach. This measure ensures a relative comparison of the obtained solutions to reference ones by quantifying the differences between them as shown in Equation 9: (9) $\mathrm{Deviation}=(X_{\mathrm{Best}}-X^{\ast })/X^{\ast }\ast 100$(9) where

• Deviation is the percent deviation;

• X_Best is the cost of the best-obtained solution with our hybrid approach;

• X* is the cost of the best-known solution in literature (reference solution);

In Tables 7–9, we provide the percent deviation for instances with 25, 50 and 100 respectively. This measure is computed for each instance based on the costs of the best-found solution X_Best using GRASP-CRO and the best-known solution X* in literature.

Table 7. Percent deviations from the best known solutions for instances with 25 customers.

Instances	R101	R102	R103	R201	R202	R203	C101	C102	C103	C201	C202	C203	RC101	RC102	RC103	RC201	RC202	RC203
Deviations	0.19	0.33	1.52	0.23	3.61	7.25	0.26	0.23	0.23	0.39	0.39	2.04	0.41	0.27	0.39	0.28	0.25	0.47

Table 8. Percent deviations from the best known solutions for instances with 50 customers.

Instances	R101	R102	R103	R201	R202	R203	C101	C102	C103	C201	C202	C203	RC101	RC102	RC103	RC201	RC202	RC203
Deviations	2.33	14.84	11.39	20.48	21.61	9.39	21.05	15.88	6.57	9.71	7.85	4.88	4.3	3.17	−4.45	0.41	1.79	3.16

Table 9. Percent deviations from the best known solutions for instances with 100 customers.

Instances	R101	R102	R103	R201	R202	R203	C101	C102	C103	C201	C202	C203	RC101	RC102	RC103	RC201	RC202	RC203
Deviations	4.21	3.47	−2.22	1.74	−4.54	5.77	37.07	24.56	24.7	17.18	17.62	18.82	3.91	0.66	5.43	1.85	−6.56	2.32

As shown in Table 7, the deviations are very close to 0% for most instances with 25 customers, particularly for instances of class C and RC. These results show that our approach can attain optimality for small instances with 25 customers geographically clustered and/or randomly scattered.

In Table 8, we present the percent deviations for instances with 50 customers. Our approach has produced good solutions for most instances of class C (C103 –> C203) with a deviation less than 10%. Furthermore, the proposed approach attained promising results for instances of class RC since deviations are less than 5%. We notice that the obtained solution for instance RC103 is better than the best-known one with a deviation of −4.45%, which shows the competitiveness of the proposed approach. These results also emphasize the good performance of our approach for instances with 50 customers geographically clustered and randomly dispersed. However, the deviations are higher for instances with only a random distribution of customers (class R).

In the case of large instances with 100 customers, percent deviations are less or equal to 5% for classes R and RC as shown in Table 9. Moreover, competitive solutions with lower costs than the best-known ones have been found (instances R103, R202 and RC202), which explains the negative percent deviations. Nevertheless, this percentage is relatively higher for problems of class C, since an average percent deviation of 23% is obtained for class C. These results show the good performances of the proposed approach for large instances with up to 100 customers randomly scattered.

8. Conclusion

In this paper, we proposed a hybrid approach to deal with the Vehicle Routing Problem with Hard Time Windows. The optimization process of our approach is based on a bio-inspired metaheuristic, which is the CRO. We integrated the GRASP in our approach to provide more diversification of the search space. After presenting the adopted solution encoding, the preliminary stage as well as the optimization process of our approach GRASP-CRO, we performed an experimental study for its evaluation. This study was conducted on the base of benchmark problems for the VRPHTW. During this study, we considered different types of instances (R, C and RC) with different sizes and various characteristics in order to define the best configuration of our approach. The obtained results have shown the competitive and promising performances of the proposed approach for problems with 25, 50 and 100 customers, since the overall deviation is very low and competitive solutions to the best-known ones have been found. These good performances are promising to apply our approach to real-world routing problems in the fields of networking and transportation.

Disclosure statement

No potential conflict of interest was reported by the author(s).