# Learning Whale Optimization Algorithm for Open Vehicle Routing Problem with Loading Constraints

Authors: Nai K. Yu; Wen Jiang; Rong Hu; Bin Qian; Ling Wang
Journal: Discrete Dynamics in Nature and Society
 (2021)

---

## Abstract

This paper addresses the two-dimensional loading open vehicle routing problem with time window (2L-OVRPTW). We propose a learning whale optimization algorithm (LWOA) to minimize the total distance; an improved skyline filling algorithm (ISFA) is designed to solve the two-dimensional loading problem. In LWOA, the whale optimization algorithm is used to search the solution space and get the high-quality solution. Then, by learning and accumulating the block structure and customer location information in the high-quality solution individuals, a three-dimensional matrix is designed to guide the updating of the population. Finally, according to the problem characteristics, the local search method based on fleet and vehicle is designed and performed on the high-quality solution region. IFSA is used to optimize the optimal individual. The computational results show that the proposed algorithm can effectively solve 2L-OVRPTW.


---

## Body

## 1. Introduction

Vehicle routing problem (VRP), as a classical combinatorial optimization problem, was first proposed by Dantzig and Ramser [1]. Open vehicle routing problem (OVRP) is a variant of VRP, which widely exists in the third-party logistics distribution. This type of problem is generating a set of vehicle routes that meet vehicle load conditions, compartment volume, and customer service demand. The vehicles do not need to return to the central depot. With the rapid development of e-commerce technology, the importance of third-party logistics is more prominent, and the loading link is an essential part of the logistics field. The loading order of goods can directly affect the service efficiency of customer points. In this context, it has important practical significance to study the two-dimensional loading open vehicle routing problem with time windows (2L-OVRPTW). 2L-OVRPTW includes “routing” and “loading.” Generally speaking, we need to consider both “routing” and “loading” to obtain a satisfactory solution. In addition, due to factors such as the weight, height, and fragility of the goods, the goods cannot be stacked on each other during the loading process, such as the distribution of large household appliances and delicate home appliances. Therefore, the study of 2L-OVRPTW has an important economic value. In theory, 2L-OVRPTW comprises a vehicle routing problem and a two-dimensional loading problem, so 2L-OVRPTW is an NP-hard problem with high complexity.The research on the vehicle routing problem and the two-dimensional loading problem has been relatively mature, but the vehicle routing problem with two-dimensional loading constraints (2L-CVRP) has gradually received attention in recent years. Research is very limited. Ma et al. [2] proposed a graph-based fuzzy evolutionary algorithm to solve the two-echelon vehicle-routing problem. With the optimization goal of minimizing the cost, each offspring is generated from the parent’s allocation graph through the graph-based fuzzy allocation process, and a fuzzy local search process is designed to further improve the performance of the offspring. Iori et al. [3] proposed 2L-CVRP for the first time in 2007 and used an exact algorithm based on branch pricing to solve it, but this algorithm can only effectively solve small- and medium-sized cases. Zachariadis et al. [4] proposed a metaheuristic algorithm for 2L-CVRP. The algorithm combines the tabu search algorithm and the guiding local search principle. In addition, the paper uses storage structure record feasible loading scheme and a set of heuristic algorithms to solve the loading problem and has achieved good results. Fuellerer et al. [5] proposed a new ant colony optimization algorithm based on saving algorithm to solve 2L-CVRP. At the same time, different heuristic algorithms were used to solve the loading problem. Good results were obtained on small-scale issues. Duhamel et al. [6] proposed a GRID × ELS algorithm for 2L-CVRP, which transforms loading constraints into resource constrained project scheduling problems. The optimization framework deals with RCPSP-CVRP and finally transforms the RCPSP-CVRP solution into 2L-CVRP solution by solving special loading problems. For 2L-CVRP, Guimarans et al. [7] proposed a hybrid simulation algorithm combining Monte Carlo simulation, an iterative local search framework, and a route and bin loading heuristic method with stochastic bias to solve the problem. Wei et al. [8] proposed a variable domain search algorithm to solve the routing problem for 2L-CVRP, used skyline heuristic algorithm to solve the loading problem, and designed a data structure to record the loading feasibility information to reduce the loading time. Leung et al. [9] proposed a metaheuristic method combining tabu search and extended guided local search (EGLS). Firstly, it is proved that tabu search algorithm can solve CVRP effectively, and EGLS can avoid the algorithm from falling into local optimal. Secondly, a new heuristic loading algorithm is added to solve the loading problem, which improves the distribution cost significantly. Wei et al. [10], aiming at 2L-CVRP, proposed a simulated annealing algorithm with repeated temperature rise and fall mechanism to solve the problem and adopted a heuristic method based on open space to solve the loading problem. Based on 2L-CVRP, Rui et al. [11] established the corresponding mathematical model considering the constraints of multidepot and time window. The former is used to solve the vehicle routing problem, and the latter is used to solve the bin loading problem. Two groups of experiments are designed to verify the effectiveness of the algorithm. Wang and Zhou [12] aimed at establishing a multiobjective vehicle scheduling problem model with two-dimensional loading constraints based on comprehensively considering the constraints of time window, two-dimensional loading, and customer satisfaction. At the same time, a multiobjective ant colony optimization algorithm was proposed. In the loading stage, the two-dimensional loading strategy of improved minimum horizontal line search algorithm was adopted to improve the success rate of loading. Ma et al. [2] proposed an adaptive localized decision variable analysis approach under the decomposition-based framework to solve large-scale multiobjective optimization problems. An adaptive strategy is used to optimize the decision variables, which can adaptively balance the convergence and diversity of solutions in the target space. Leung et al. [13] established a mathematical model to minimize the transportation cost for the heterogeneous vehicle routing problem with two-dimensional loading constraints. A simulated annealing algorithm combined with heuristic local search is proposed, and a set of bin loading heuristic algorithms are designed to solve the bin loading problem. Zhang et al. [14] took the total cost minimization as the optimization objective, studied the heterogeneous vehicle routing problem with loading constraints, designed a hybrid particle swarm optimization algorithm combining artificial bee colony and artificial immune to solve the problem, and proposed several strategies to avoid the algorithm falling into local optimum. Sabar et al. [15] proposed a two-stage algorithm to solve the heterogeneous vehicle routing problem with two-dimensional bin loading constraints. In the first stage, an adaptive meme algorithm was designed to solve the routing problem. In the second stage, a hybrid algorithm based on five heuristic algorithms was used to solve the bin loading problem. Dominguez et al. [16] proposed a hybrid algorithm based on neighborhood search and loading heuristic algorithm to solve the two-dimensional vehicle routing problem with cluster backhaul and used randomization technology to guide the search process. Good results were obtained. Wang et al. [17] established a mathematical model for the logistics distribution problem with two-dimensional loading and proposed a memetic algorithm to solve the problem. The algorithm integrates the heuristic search knowledge of loading and route and improves the performance and result quality of existing algorithms. Literature research shows that the research on 2L-OVRPTW is still in its infancy.The 2L-OVRPTW can establish two different problem models. One is the 0-1 mixed integer model (MIP) and the other is the order model. For MIP, the most exact algorithms such as branch and bound and column generation mainly solve the problem. The solution space can be searched through commercial solvers (such as GUROBI or CPLEX) for small-scale problems. The optimal solution can be found reasonably, but the exact algorithm will show a geometric growth trend with the increase of the problem size for large-scale problems. For the order model, it is usually composed of customer sequence of vehicle service. The constraints of such problems are often hidden in the order model, mainly by heuristic algorithm and intelligent algorithm. The heuristic algorithm is usually based on the characteristics of the problem to construct feasible solutions, such as Johnson algorithm, palmar algorithm, and Gupta algorithm. This kind of algorithm can obtain a feasible solution quickly, but it cannot guarantee that the solution is the best solution in the global range. Intelligent algorithm is a combination of the characteristics of the problem and intelligent mechanism to solve; through simulating some mechanisms in nature to guide the algorithm search, often in a short time, we can get a satisfactory solution or approximate optimal solution of complex scheduling problems. He pointed out in [18] that the effectiveness of intelligent algorithm depends on the model characteristics of scheduling problem and the mechanism of the intelligent algorithm. The intelligent optimization algorithm only needs to spend a short time searching part of the solution space of the scheduling model to find a satisfactory solution so that it can effectively solve all kinds of complex scheduling problems.The whale optimization algorithm (WOA) is a new swarm intelligence optimization algorithm proposed by Seyedali and Andrew [19]. WOA simulates the predator-prey process of humpback whales in nature and realizes the search of high-quality solution by searching, encircling, and attacking. WOA is widely used in many fields [20–23] because of its advantages of simple operation, few parameters, easy understanding, and strong optimization performance and has achieved good results. From the current research status of WOA, reasonable selection of adaptive weights in the search process can effectively control the search direction and further improve the search ability of the algorithm. At present, WOA is mainly used to solve continuous optimization problems, but for discrete optimization problems, the real solution space needs to be discretized. WOA has a certain effect on global search, but it cannot avoid the problem that the high-quality block structure is destroyed in the process of generating new species group. Therefore, WOA is used to construct a three-dimensional matrix after the first stage search in the solution space and to search in the second stage. The three-dimensional matrix is used to learn and accumulate the block structure and its location information based on customer service order. Then, the population updating process is guided to improve the search efficiency of learning whale optimization algorithm (LWOA), for example, for the customer service order of high-quality solution [1,3,5,2,4] and [2,4,5,1,3], where [2, 4] is a block structure based on the order relationship of customer 2 before customer 4 and [2, 4] appears in both the fourth position of the first solution and the first position of the second solution. Therefore, the three-dimensional matrix can retain the block structure of the first two and the last four and also record the location information of the block structure, that is, the block structure [2, 4] can be stored in the elements with subscripts (4,2,4) and (1,2,4), respectively, in the three-dimensional matrix (the first dimension subscript in the three-dimensional matrix represents the location information, and the second and third dimensions represent the order relationship of customers).In this paper, the model and solution of 2L-OVRPTW are studied. In the model aspect, the 2L-OVRPTW model integrated by OVRPTW and BPP is established, and the optimization objective of this problem model is to minimize the driving distance. In the process of solving, through the analysis of the properties of the problem, a learning whale optimization algorithm (LWOA) is proposed to optimal vehicle routing problem, and an improved skyline method is designed to solve the loading problem. Specifically, three rules are used to initialize the population, and LOV rules are used to discretize the real solution. Secondly, in the global search, the whale optimization algorithm is used to find the optimal solution in the first stage. In the second stage, the block structure based on the order relation and the three-dimensional matrix of its position information are designed to effectively learn the solution and accumulate relevant information. Finally, the local search which combines the exchange operation, swap operation, and insert operation is introduced to carry out the efficient and detailed local search for the solution obtained from the two-stage global search, which can effectively balance the global and local search of LWOA and improve the overall search ability of the algorithm. Simulation results and algorithm comparison show that LWOA is an effective algorithm for 2L-OVRPTW.The rest of this paper is arranged as follows. In Section2, the definition of the OVRPTW is introduced and a corresponding integer programming model is also proposed. In Section 3, the LWOA algorithm is proposed and described in detail. In Section 4, experimental results and comparisons are given and discussed. The current work is summarized and the future research direction is prospected in Section 5.


## 2. Two-Dimensional Loading Open Vehicle Routing Problem with Time Window

### 2.1. Problem Description and Symbol Definition of the 2L-OVRPTW

The 2L-OVRPTW can be described as follows. In directed network graphG=V,E, V=0,1,2,…,N represents the vertex set, where vertex 0 represents central depot and V′=V\0 represents the customer set. E=i,j|i,j∈V,i≠j is the arc set, and dij is the distance from customer i to customer j. There are Nk homogeneous vehicles in the central depot for service, the capacity of each vehicle is Q, and the length and width of the carriage areL,W; the carriage area is S=L∗W. Mi is the total number of rectangular goods required by customer i, qi is the total weight of goods required by customer i, and qik is the weight of the kth item required by customer i; then, qi=∑k=1miqik. Mig is the g goods required by customer i, lig is the length of the g goods required by customer i, and wig is the width of the g goods required by customer i; then, sig=lig⋅wig. At the same time, considering the needs of customer satisfaction, the time window for customer i to receive service is set as ei,hi, ei means the earliest service time allowed, hi means the latest service time, and the service time of each customer is Ti. If the service starts within the time window of customer acceptance, the satisfaction is 1, otherwise it is 0. The time to arrive at the customer is sti. If the service arrives in advance, it needs to wait until the service start time.The objective of 2L-OVRPTW is to minimize the distance and consider the following constraints:(1)
Each vehicle starts from the central depot and does not need to return back after serving all customers(2)
When serving current customers, it is not necessary to move the goods of other customers (the principle of last in first out)(3)
The loading edge of all articles must be parallel to the carriage edge(4)
All goods in the carriage shall not be overlapped(5)
The goods loaded on each vehicle shall not exceed the length and width of the carriage and the maximum load of the vehicle(6)
Each vehicle must be serviced after the service time allowed by the customerFigure1 shows a distribution route of the 2L-OVRPTW. It can be seen from Figure 1 that three homogeneous vehicles serve 11 customers in the central depot. The service routes of vehicle 1, vehicle 2, and vehicle 3 are route1=0,1,2,3, route2=0,4,5,6,7, and route3=0,8,9,10,11. The loading scheme of vehicle 2 is shown in Figure 2. The goods are placed at the lower-left corner of the top view of the carriage (the front of the carriage is on the left). When unloading the goods, the current customer’s goods are unloaded through a straight line parallel to the length of the carriage. The thick line in the figure represents the skyline, the thick line parallel to the length of the carriage is the longitudinal skyline, and the thick line parallel to the width of the carriage is the transverse skyline; the shaded part is the redundant space (wasted space) of the vehicle.Figure 1
Distribution route of 2L-OVRPTW.Figure 2
Vehicle loading scheme.


### 2.2. Mathematical Model of the 2L-OVRPTW

(1)minf1=∑k=1Nk∑i,j∈vdijxkij,(2)s.t.∑i=1Nqiyik≤Q,k∈0,1,…,Nk,(3)∑k=1Nkyik=1,∀i∈V′,(4)∑i=1Nxijk=yik,∀j∈V′,k∈0,1,…,Nk,(5)∑j=1Nxijk=yik,∀i∈V′,k∈0,1,…,Nk,(6)0≤vkig≤W−wkig,k∈0,1,…,Nk,∀i∈V′,∀g∈Mi,(7)0≤hkig≤L−Ikig,k∈0,1,…,NK,∀i∈V′,∀l∈Mi,(8)HHii′gg′=maxyik⋅Hkig+wig,yi′k′⋅Hki′g′+wi′g′−yik⋅wig+yi′k′⋅wi′g′,∀i,i′∈V′,∀g,g′∈Mi,g≠g′,k∈0,1,…,Nk,(9)LLii′gg′=maxyik⋅Vkig+lig,yi′k′⋅Vki′g′+li′g′−yik⋅lig+yi′k′⋅li′g′,∀i,i′∈V′,∀g,g′∈Mi,g≠g′,k∈0,1,…,Nk,(10)maxHHii′gg′,LLii′gg′≥0,∀i,i′∈V′,∀g,g′∈Mi,(11)maxHHii′gg′,Hkig−Hki′g′≥0,∀i,i′∈V′,∀g,g′∈Mi,g≠g′,(12)yki⋅sti≥yki⋅ei,∀i∈V′,k∈0,1,…,Nk,(13)ei≤yki⋅sti+Ti≤hi,∀i∈V′,k∈0,1,…,Nk,(14)xijk=1,if customeriis visited before customerjby the vehiclek,0,otherwise,(15)yijk=1,if customeriis visited by the vehiclek,0,otherwise.In the above model, (1) gives the objective function of 2L-OVRPTW. Constraint (1) represents the minimization of the total distance. Constraint (2) ensures that the load of the vehicle does not exceed the maximum load of the vehicle. Constraint (3) ensures that every customer can be served. Constraints (4) and (5) ensure that every customer has and only has one vehicle to serve. Constraints (6)–(11) are the loading constraints, in which constraints (6) and (7) ensure that all items are fully loaded into the compartment. Whereas wkig and lkig, respectively, represent the g item width and length of the k vehicle serving the i customer. Hkig and Vkig, respectively, represent the horizontal and vertical coordinates of the upper-left corner of the g item of the i customer served by the k vehicle. Constraints (8)–(10) show that customers cannot overlap. Constraint (11) satisfies the principle of first in and last out. Constraints (12) and (13) show the time window constraints. Constraints (14) and (15) show the value range of decision variables.


## 2.1. Problem Description and Symbol Definition of the 2L-OVRPTW

The 2L-OVRPTW can be described as follows. In directed network graphG=V,E, V=0,1,2,…,N represents the vertex set, where vertex 0 represents central depot and V′=V\0 represents the customer set. E=i,j|i,j∈V,i≠j is the arc set, and dij is the distance from customer i to customer j. There are Nk homogeneous vehicles in the central depot for service, the capacity of each vehicle is Q, and the length and width of the carriage areL,W; the carriage area is S=L∗W. Mi is the total number of rectangular goods required by customer i, qi is the total weight of goods required by customer i, and qik is the weight of the kth item required by customer i; then, qi=∑k=1miqik. Mig is the g goods required by customer i, lig is the length of the g goods required by customer i, and wig is the width of the g goods required by customer i; then, sig=lig⋅wig. At the same time, considering the needs of customer satisfaction, the time window for customer i to receive service is set as ei,hi, ei means the earliest service time allowed, hi means the latest service time, and the service time of each customer is Ti. If the service starts within the time window of customer acceptance, the satisfaction is 1, otherwise it is 0. The time to arrive at the customer is sti. If the service arrives in advance, it needs to wait until the service start time.The objective of 2L-OVRPTW is to minimize the distance and consider the following constraints:(1)
Each vehicle starts from the central depot and does not need to return back after serving all customers(2)
When serving current customers, it is not necessary to move the goods of other customers (the principle of last in first out)(3)
The loading edge of all articles must be parallel to the carriage edge(4)
All goods in the carriage shall not be overlapped(5)
The goods loaded on each vehicle shall not exceed the length and width of the carriage and the maximum load of the vehicle(6)
Each vehicle must be serviced after the service time allowed by the customerFigure1 shows a distribution route of the 2L-OVRPTW. It can be seen from Figure 1 that three homogeneous vehicles serve 11 customers in the central depot. The service routes of vehicle 1, vehicle 2, and vehicle 3 are route1=0,1,2,3, route2=0,4,5,6,7, and route3=0,8,9,10,11. The loading scheme of vehicle 2 is shown in Figure 2. The goods are placed at the lower-left corner of the top view of the carriage (the front of the carriage is on the left). When unloading the goods, the current customer’s goods are unloaded through a straight line parallel to the length of the carriage. The thick line in the figure represents the skyline, the thick line parallel to the length of the carriage is the longitudinal skyline, and the thick line parallel to the width of the carriage is the transverse skyline; the shaded part is the redundant space (wasted space) of the vehicle.Figure 1
Distribution route of 2L-OVRPTW.Figure 2
Vehicle loading scheme.


## 2.2. Mathematical Model of the 2L-OVRPTW

(1)minf1=∑k=1Nk∑i,j∈vdijxkij,(2)s.t.∑i=1Nqiyik≤Q,k∈0,1,…,Nk,(3)∑k=1Nkyik=1,∀i∈V′,(4)∑i=1Nxijk=yik,∀j∈V′,k∈0,1,…,Nk,(5)∑j=1Nxijk=yik,∀i∈V′,k∈0,1,…,Nk,(6)0≤vkig≤W−wkig,k∈0,1,…,Nk,∀i∈V′,∀g∈Mi,(7)0≤hkig≤L−Ikig,k∈0,1,…,NK,∀i∈V′,∀l∈Mi,(8)HHii′gg′=maxyik⋅Hkig+wig,yi′k′⋅Hki′g′+wi′g′−yik⋅wig+yi′k′⋅wi′g′,∀i,i′∈V′,∀g,g′∈Mi,g≠g′,k∈0,1,…,Nk,(9)LLii′gg′=maxyik⋅Vkig+lig,yi′k′⋅Vki′g′+li′g′−yik⋅lig+yi′k′⋅li′g′,∀i,i′∈V′,∀g,g′∈Mi,g≠g′,k∈0,1,…,Nk,(10)maxHHii′gg′,LLii′gg′≥0,∀i,i′∈V′,∀g,g′∈Mi,(11)maxHHii′gg′,Hkig−Hki′g′≥0,∀i,i′∈V′,∀g,g′∈Mi,g≠g′,(12)yki⋅sti≥yki⋅ei,∀i∈V′,k∈0,1,…,Nk,(13)ei≤yki⋅sti+Ti≤hi,∀i∈V′,k∈0,1,…,Nk,(14)xijk=1,if customeriis visited before customerjby the vehiclek,0,otherwise,(15)yijk=1,if customeriis visited by the vehiclek,0,otherwise.In the above model, (1) gives the objective function of 2L-OVRPTW. Constraint (1) represents the minimization of the total distance. Constraint (2) ensures that the load of the vehicle does not exceed the maximum load of the vehicle. Constraint (3) ensures that every customer can be served. Constraints (4) and (5) ensure that every customer has and only has one vehicle to serve. Constraints (6)–(11) are the loading constraints, in which constraints (6) and (7) ensure that all items are fully loaded into the compartment. Whereas wkig and lkig, respectively, represent the g item width and length of the k vehicle serving the i customer. Hkig and Vkig, respectively, represent the horizontal and vertical coordinates of the upper-left corner of the g item of the i customer served by the k vehicle. Constraints (8)–(10) show that customers cannot overlap. Constraint (11) satisfies the principle of first in and last out. Constraints (12) and (13) show the time window constraints. Constraints (14) and (15) show the value range of decision variables.


## 3. The LWOA for Solving 2L-OVRPTW

### 3.1. Global Search

In this section, the population initialization, standard WOA, block structure and three-dimensional matrix, local search, and loading strategy are introduced. Finally, the loading flowchart and the flowchart of LWOA are given.


#### 3.1.1. Population initialization

In order to ensure the diversity and dispersion of the solutions, two individuals are generated byrule1 and rule2, and then, the rest individuals are generated by rule3. Rule1 considers the nearest neighbor principle, rule2 considers the highest satisfaction principle, and rule3 adopts the quasi-opposition strategy. Placing all unserved customers in set Z, the specific rules are as follows:Rule1 is the nearest neighbor principle.After the vehicle departs from the central depot, the distance between the central depot and each customer inZ is calculated, and the nearest customeri (if there are multiple customers, the customer with the highest satisfaction) is selected for service. The customeri from set Z was removed and the distances between customer i and other customers inZ are calculated. After selecting the nearest customer j, if the constraints are still met when continuing to serve the customer j, it is necessary to remove the customer j from set Z and calculate the distance between the customer and other customers in Z. Otherwise, the vehicle will give up serving the customer and dispatch the vehicle from the central depot to serve the customers in Z again. This cycle will last until all customers inZ are served.Rule2 is the maximum satisfaction principle.After vehicle starts from the central warehouse, we calculate the arrival time of each customer inZ from the central depot and select the customer i meeting the time window constraint to serve. If there are multiple customers whose arrival time is within the time window, the nearest customer i is selected to serve. Move customer i out of set Z and calculate the time from customer i to other customers in Z. Select customer j that meets the time window constraint. If customer j still meets the constraint, it is necessary to move customer j out of set Z and calculate the arrival times from customer Z to other customers in Z. Otherwise, the vehicle will give up serving the customer and send the vehicle from the central depot to serve the customers in Z again until all the customers in Z are served.Rule3 is quasi-antagonistic strategy.Let the population size bepopsize and use the quasi-opposition strategy to generate popsize-2 individuals with n customers (that is, the search space is n). Let πigen=πigen1,πigen2,…,πigeng,…,πigenSπigeng∈V,n+1≤S≤2n denote thei individual (distribution scheme) in the gen population. Let Pigen=Pigen1,Pigen2,…,Pigeng,…,PigennPigeng∈V′ be the customer-based discrete sort in individual πigen (i.e., the customer sequence without central depot), and let Xigen=Xigen1,Xigen2,…,Xigeng,…,XigennXigeng∈a,b be the real number sequence corresponding to the customer’s discrete sort Pigen. As shown in equation (16), the quasi-antagonistic strategy is to obtain the opposite inverse solution (quasi-antagonistic solution) Xjgen=Xjgen1,Xjgen2,…,Xjgeng,…,XjgennXjgeng∈a,b of the gen generation individual i through the calculation of equation (16) in the real number field a,b. Therefore, the population with quasi-opposite solution can be obtained when the population is randomly generated. popsize-2 individuals are selected from the two populations as the next generation population:(16)Xjgeng=randX¯,Xigeng,Xigeng>X¯,randXigeng,X¯,Xigeng≤X¯,X¯=a+b2.The population obtained byXigen through rule3 is a real number sequence over a continuous field. Therefore, this paper uses LOV (large-order value) rule to transform the individual real number sequence Xigen into discrete customer sequence Pigen.The specific steps of LOV rules are as follows:Step 1: we select the real number sequence Xigen=Xigen1,Xigen2,…,Xigeng,…,Xigenn, arrange Xigen in the descending order to get Xigen′, find the corresponding sequence number Xigen in Xigen′g, and place it in φigeng to get the intermediate sequence φigenStep 2: for discrete customers, the values in Pigen are obtained successively by the formula Pigenφigeng=gAs shown in Table1, an example of LOV application is given. If the real number sequence of an individual is Xigen=2.30,1.59,2.17,1.78,2.06,1.89,0.88, the real number sequence after descending is Xigen′=2.30,2.17,2.06,1.89,1.78,1.59,0.88. Through the above steps, the intermediate sequence φigen=1,6,2,5,3,4,7 and the customer’s discretization sequence Pigen=1,3,5,6,4,2,7 are obtained.Table 1
Solution representation of LOV.
Dimension (g)1234567Xigeng2.301.592.171.782.061.890.88φigeng1625347Pigeng1356427By calling the LOV rule, we can realize the transformation from real sequence to discrete sequence. Similarly, if the customer’s discrete sequence changes, we can adjust it by using the reverse largest order value (RLOV) rule. Therefore, the RLOV rule is used to real number the discrete customer sequence obtained byrule1 and rule2 to get the corresponding real number solution. The initialization population is composed of popsize-2 individuals generated by rule3.


#### 3.1.2. Standard WOA

Wang et al. [17] designed a whale optimization algorithm (WOA) by studying the foraging behavior of humpback whales. WOA mainly includes three search mechanisms: encircling search mechanism, spiral search mechanism, and random search mechanism.Encircling search mechanism: when finding prey (solution), individual whales communicate with each other in the whole population so that the individual whales in the population keep approaching the prey. The mathematical model of encircling search mechanism is shown in formula (17)∼(18):(17)D=CXgen∗−Xigen,(18)Xigen+1=Xgen∗−A⋅D,where gen is the current number of iterations, Xigen is the i individual in the gen generation, Xgen∗ is the position of prey in the gen generation, A and C are coefficient vectors, a is the control parameter, t is the current number of iterations, each iteration updates the population, and Tmax is the maximum number of iterations, with the increase of the number of iterations.Spiral search mechanism: after calculating the distance from the prey, the individual whale takes the prey as the center and constantly approaches it in a spiral form. The mathematical model of spiral search mechanism is shown as(19)Xigen+1=D⋅eblcos2πl+Xigen,where b is the parameter controlling the shape of the spiral, l is the random value of [−1,1], and γ is the search factor.In equation (20), the probability p is introduced to choose whether to update the position information of individual whale by using the encircling search mechanism or the spiral search mechanism:(20)Xigen+1=Xigen−A⋅D,p<γ,D⋅eblcos2πl+Xgen∗,p≥γ.Random search mechanism: when the value ofA is greater than 1, it will search outside the bounding circle. When we cannot get the effective information of the prey, we will choose a random way to search the prey to update the position information of the population, which can effectively enhance the diversity of the population. Then, the mathematical model of random search mechanism is as follows:(21)Xigen+1=Xgen′−A⋅D,(22)D=CXgen′−Xigen,where Xgen′ is a randomly selected individual in the gen generation population.


#### 3.1.3. Block Structure and Three-Dimensional Matrix

For any solution of 2L-OVRPTW, the block structure refers to the order of two consecutive adjacent customers in the solution. In the process of updating, the solution set shows certain similarities in the structure of the solution (block structure and its location information). The three-dimensional matrix proposed in the second stage of global search can effectively learn and accumulate the block structure. The location information of the solution set is obtained in the first stage so as to improve the overall effect of the algorithm. Because the same block structure in different positions will have a great impact on the target value of the solution, the first dimension of the designed three-dimensional matrix is used to record the location information of the block structure, and the other two dimensions represent two customers in the block structure.We defineBest_Pgen=Best_πgen,1,Best_πgen,2,⋯Best_πgen,g,⋯,Best_πgen,n, where Best_Pgen is the gen generation population, and Best_πgen,w=Best_π1gen,w,Best_π2gen,w,⋯Best_πggen,w,⋯,Best_πngen,w, where Best_πgen,w is the w individual and Best_πggen,w is the g customer in the w individual of gen population. Let TMl×n×ngen represent the three-dimensional matrix based on block structure and its location information in gen generation population, where TMl×n×ngenp,x,y represents the value of block structure x,y in p position and PMl×n×ngenp is the two-dimensional probability matrix normalized at p position. The specific definition is as follows:(23)ONE_TMl×n×ngen,wp,x,y1,x=Best_πpgen,w,y=Best_πp+1gen,w,0,else,p=1,2,⋯,lx,y=1,2,⋯,n,(24)TMl×n×ngenp,x,y=∑w=1popsizeONE_TMl×n×ngen,wp,x,y,p=1,2,⋯,lx,y=1,2,⋯,n,(25)TMl×n×ngenp,x=TMl×n×ngenp,x,1,TMl×n×ngenp,x,2,⋯,TMl×n×ngenp,x,n,p=1,2,⋯,lx=1,2,⋯,n,(26)TMl×n×ngenp=TMl×n×ngenp,x,1,⋯,TMl×n×ngenp,x,n=TMl×n×ngenp,1,1,⋯,TMl×n×ngenp,1,n⋮TMl×n×ngenp,n,1,⋯,TMl×n×ngenp,n,n,(27)TMl×n×ngen=TMl×n×ngen1,TMl×n×ngen2,⋯,TMl×n×ngenl,(28)PMl×n×ngenp=PMl×n×ngenp,x,1,⋯,PMl×n×ngenp,x,n=PMl×n×ngenp,1,1,⋯,PMl×n×ngenp,1,n⋮PMl×n×ngenp,n,1,⋯,PMl×n×ngenp,n,n,(29)∑x=1n∑y=1nPMl×n×ngenp,x,y=1,p∈1,2,⋯,l,where ONE_TMl×n×ngen,wp,x,y is the number of times that customer x=Best_πpgen,w and customer y=Best_πp+1gen,w appears in p position of individual w in the gen generation, that is, the number of times that block structure x,y appears in p position of individual w. For example, for individual Best_πgen,1=1,3,2,5,4,6 and individual Best_πgen,2=2,5,3,6,1,4, let p=2; then, the client block structure of individual 1 in position 2 is [3,2], and the client block structure of individual 2 in position 2 is [5,3]. Therefore, ONE_TMl×n×ngen,12,3,2=1, ONE_TMl×n×ngen,22,3,2=0, ONE_TMl×n×ngen,12,5,3=0, and ONE_TMl×n×ngen,22,5,3=1. As shown in Figure 3, the block structure information of each p position is recorded in TMl×n×ngenp.Figure 3
Three-dimensional matrix based on block structure.In the first stage of global search, the improved WOA is used to search in the whole solution space to obtain high quality. In the second stage, a three-dimensional matrix TM is constructed based on the high-quality solution obtained in the first stage. When the population is updated by using the three-dimensional matrix, the structure and location information of customer blocks in the current high-quality solution is learned and accumulated after each update and used to generate new individuals. The specific update process is as follows:Step 1: construct a three-dimensional probability matrix TMl×n×ngen by using the high-quality solution obtained by WOA and construct a list cos_sequence of remaining optional customers and a list tabu of selected customers. Perform step 2.Step 2: let p=1, and randomly generate a random number a,a∈0,1 and execute step 3.Step 3: select the matrix of position p, set the occurrence times of block structure containing customers in tabu to 0, and get the probability matrix of position p after TMl×n×ngenp normalization as PMl×n×ngenp; sum=∑x=1n∑y=1nPMl×n×ngenp,x,y; if sum>0, perform step 4; otherwise, perform step 5.Step 4: if there is a selectable block structure, the roulette method is used to select the customer block structure for PMl×n×ngenp. Let p=p+2; perform step 6.Step 5: if there is no optional block structure in the current position, randomly select a customer in cos_sequence and place it in this position; let p=p+1, and then, go to step 6.Step 6: delete the selected customers from cos_sequence and add the selected customers to tabu.Step 7: if p>num_customer, output the generated new individual tabu; if p=num_customer, return to step 5; otherwise, return to step 3.The pseudocode of the LWOA algorithm is as follows.Algorithm 1: Learning whale optimization algorithm.
(1)
Input:(2)
Let the number of customers beN.Mn represents the total number of items required by the customer n.gen be the current number of iterations; the maximum number of iterations is genmax.πigen is the individual i of generation gen, and the population size is popsize.γ is the search factor.The update times of three-dimensional matrix ismax_l.The vehicle number isk,k∈K.Ing indicates the item g of the customer n(3)
LetA=2ar1−a,C=2r2,r1,r2∈random0,1, a=2−2∗gen/genmax, b=1,l∈random−1,1, p∈random0,1(4)
begin:(5)
Population initialization: using rule1 and rule2 in Section 3.1.1 to generate an individual and using rule3 to generate popsize−2 individuals.(6)
for each gen∈genmaxdo(7)
for each i∈popsizedo(8)
ifA>1then(9)
updateπigen by equation (21)(10)
else(11)
updateπigen by equation (20)(12)
end if(13)
end for(14)
end for(15)
for each gen∈max_ldo(16)
Choosing the best individualsπgen to construct the three-dimensional matrix(17)
Update the population by the three-dimensional matrix(18)
end for(19)
Choosing the best individual asπbest∗(20)
Loading initialization: Initialize the loading sequence of all customers using rules 1 and 2 in 3.4.1.(21)
for each k∈πbest∗do(22)
for each n∈kdo(23)
for each Ing∈ndo(24)
Loading the customer’s goods into the carriage through the loading strategy in Section3.4.2(25)
end for(26)
end for(27)
end for(28)
Output:(29)
returnπbest∗The values range ofγ, popsize, and max_l are analyzed in detail in Section 4.2.


### 3.2. Local Search

According to the characteristics of the 2L-OVRPTW problem, the solution space is huge, and the high-quality solutions are usually distributed in each local region of the solution space. Therefore, it is necessary to search the high-quality solutions from the global search; VNS includes three neighborhood operations: interchange operation, insert operation, and swap operation. Interchange operation is used between vehicles, and insert operation and swap operation are used inside vehicles. In this way, local search can be carried out more carefully to update the current optimal solution. The variable neighborhood local search execution pseudocode is as follows.Algorithm 2: Variable neighborhood local search.
(1)
Input:(2)
Choose the best individualπi and convert it to the permutation π0 according to the LOV rule. Perform a neighborhood operation on π0 to get π(3)
begin:(4)
Set loop = 1, max-loop = 40,n = 1, max-n = 20(5)
whilen < max-ndo(6)
Random selectm-k, m-k∈ {1, 2, 3}(7)
Randomly selectk and l, where k≠ l(8)
ifm-k = 1, thenπ′=insertπ,k,l(9)
else ifm-k = 2, thenπ′=swapπ,k,l(10)
else if m-k = 3, thenπ′=inverseπ,k,l(11)
end if(12)
end if(13)
end if(14)
ifπ′≺π, thenπ=π′(15)
elsen = n + 1(16)
end if(17)
end while(18)
π′=interchangeπ,k,l(19)
ifπ′≺π,then π=π′(20)
end if(21)
until loop = max-loop(22)
ifπ≺π0, then π0=π(23)
end if(24)
Output:(25)
Convertπ0 back toπi according to the RLOV rule


### 3.3. Complexity Analysis of the LWOA

Let the population size be popsize, the number of running iterations of global search be genglobal, and the number of iterations of each generation of local search be genlocal. The number of times that local search actually uses neighborhood operators in each generation of search is Z(Z = 20 ∗ 40). The algorithm complexity of LWOA TLWOA is composed of global search complexity and local search complexity, respectively. In each generation, the global search complexity consists of population initialization complexity O(popsize ∗ N, where N is the number of customers), discretization complexity (2 ∗ popsize ∗ N!), three-dimensional probability matrix update complexity O(N2 ∗ (N − 1)), complexity O(popsize∗ (1 + N2 ∗ (N − 1))) of generating new population based on three-dimensional probability matrix, and population evaluation complexity (popsize ∗ N); the result is O(popsize ∗ (1 + N2 ∗ (N − 1))). The complexity of local search is O(genlocal ∗ Z ∗ (N)) (the worst-case scenario).Through the above analysis,TLWOA = O(genglobal ∗ (popsize ∗ (1 + N2 ∗ (N − 1)) + genlocal ∗ Z ∗ (N))).


### 3.4. Loading Strategy

For the loading problem based on the principle of “first in, last out” in the customer service order, this paper adopts the loading strategy of the skyline filling method to check every feasible solution generated. Only when the loading constraint is satisfied can the solution be retained. This section mainly includes initialization of loading order and loading mechanisms.


#### 3.4.1. Initialize Loading Sequence

Each time we get a solution, we need to make sure that all the customers’ goods can be loaded into the car in the reverse order of the customers’ order. However, the loading order of multiple goods of each customer is not unique. Therefore, this paper proposes two rules for the initial loading order of customers’ goods.Rule1: sort items according to their lengthRule2: sort items by area from small to largeFor example, if the vehicle service sequence is customer 1, customer 2, customer 3, and customer 4, customer 4 of the final service usesrule1 for initial loading and the other three customers use rule2 for initial loading.


#### 3.4.2. Loading Method

The loading process of articles is based on the skyline. Initially, the horizontal skyline and the vertical skylines are the width of the carriage and the length of both sides of the vehicle, respectively. The loading position of articles is determined by the score of articles placed on each section of the horizontal skyline. (1) When the width is equal to the horizontal skyline, if the length is equal to either side of the adjacent vertical skyline, the score is 5. If the length is lower than the longitudinal skyline of the lowest adjacent side, the score is 4. If the length is between the two sides, the score is 3. Otherwise, the score is 2. (2) If the width is less than the transverse skyline, the score is 4 if the length is equal to either side of the longitudinal skyline of the two adjacent sides. If the length is lower than the longitudinal skyline of the lowest adjacent side, the score is 3. Otherwise, the score is 2. (3) If the width is greater than the transverse skyline, the score is 1. If all the horizontal skylines are less than the width of the object, it is necessary to fill the skyline before loading the object. In addition, when placing the object, it is also necessary to consider whether it exceeds the length and width limits of the vehicle. Loading the object into the vehicle and updating the skyline at the position with the highest score.


### 3.5. Flowchart of LWOA

This section gives the flowchart of LWOA. In the whole operation process, the individual needs to be boxed after each update of the individual, and the individual meeting the boxing constraints will be retained. The overall schematic diagram of the LWOA algorithm and the detailed flowchart of local search are shown in Figures4 and 5.Figure 4
LWOA flowchart.Figure 5
Local search flowchart.


## 3.1. Global Search

In this section, the population initialization, standard WOA, block structure and three-dimensional matrix, local search, and loading strategy are introduced. Finally, the loading flowchart and the flowchart of LWOA are given.


### 3.1.1. Population initialization

In order to ensure the diversity and dispersion of the solutions, two individuals are generated byrule1 and rule2, and then, the rest individuals are generated by rule3. Rule1 considers the nearest neighbor principle, rule2 considers the highest satisfaction principle, and rule3 adopts the quasi-opposition strategy. Placing all unserved customers in set Z, the specific rules are as follows:Rule1 is the nearest neighbor principle.After the vehicle departs from the central depot, the distance between the central depot and each customer inZ is calculated, and the nearest customeri (if there are multiple customers, the customer with the highest satisfaction) is selected for service. The customeri from set Z was removed and the distances between customer i and other customers inZ are calculated. After selecting the nearest customer j, if the constraints are still met when continuing to serve the customer j, it is necessary to remove the customer j from set Z and calculate the distance between the customer and other customers in Z. Otherwise, the vehicle will give up serving the customer and dispatch the vehicle from the central depot to serve the customers in Z again. This cycle will last until all customers inZ are served.Rule2 is the maximum satisfaction principle.After vehicle starts from the central warehouse, we calculate the arrival time of each customer inZ from the central depot and select the customer i meeting the time window constraint to serve. If there are multiple customers whose arrival time is within the time window, the nearest customer i is selected to serve. Move customer i out of set Z and calculate the time from customer i to other customers in Z. Select customer j that meets the time window constraint. If customer j still meets the constraint, it is necessary to move customer j out of set Z and calculate the arrival times from customer Z to other customers in Z. Otherwise, the vehicle will give up serving the customer and send the vehicle from the central depot to serve the customers in Z again until all the customers in Z are served.Rule3 is quasi-antagonistic strategy.Let the population size bepopsize and use the quasi-opposition strategy to generate popsize-2 individuals with n customers (that is, the search space is n). Let πigen=πigen1,πigen2,…,πigeng,…,πigenSπigeng∈V,n+1≤S≤2n denote thei individual (distribution scheme) in the gen population. Let Pigen=Pigen1,Pigen2,…,Pigeng,…,PigennPigeng∈V′ be the customer-based discrete sort in individual πigen (i.e., the customer sequence without central depot), and let Xigen=Xigen1,Xigen2,…,Xigeng,…,XigennXigeng∈a,b be the real number sequence corresponding to the customer’s discrete sort Pigen. As shown in equation (16), the quasi-antagonistic strategy is to obtain the opposite inverse solution (quasi-antagonistic solution) Xjgen=Xjgen1,Xjgen2,…,Xjgeng,…,XjgennXjgeng∈a,b of the gen generation individual i through the calculation of equation (16) in the real number field a,b. Therefore, the population with quasi-opposite solution can be obtained when the population is randomly generated. popsize-2 individuals are selected from the two populations as the next generation population:(16)Xjgeng=randX¯,Xigeng,Xigeng>X¯,randXigeng,X¯,Xigeng≤X¯,X¯=a+b2.The population obtained byXigen through rule3 is a real number sequence over a continuous field. Therefore, this paper uses LOV (large-order value) rule to transform the individual real number sequence Xigen into discrete customer sequence Pigen.The specific steps of LOV rules are as follows:Step 1: we select the real number sequence Xigen=Xigen1,Xigen2,…,Xigeng,…,Xigenn, arrange Xigen in the descending order to get Xigen′, find the corresponding sequence number Xigen in Xigen′g, and place it in φigeng to get the intermediate sequence φigenStep 2: for discrete customers, the values in Pigen are obtained successively by the formula Pigenφigeng=gAs shown in Table1, an example of LOV application is given. If the real number sequence of an individual is Xigen=2.30,1.59,2.17,1.78,2.06,1.89,0.88, the real number sequence after descending is Xigen′=2.30,2.17,2.06,1.89,1.78,1.59,0.88. Through the above steps, the intermediate sequence φigen=1,6,2,5,3,4,7 and the customer’s discretization sequence Pigen=1,3,5,6,4,2,7 are obtained.Table 1
Solution representation of LOV.
Dimension (g)1234567Xigeng2.301.592.171.782.061.890.88φigeng1625347Pigeng1356427By calling the LOV rule, we can realize the transformation from real sequence to discrete sequence. Similarly, if the customer’s discrete sequence changes, we can adjust it by using the reverse largest order value (RLOV) rule. Therefore, the RLOV rule is used to real number the discrete customer sequence obtained byrule1 and rule2 to get the corresponding real number solution. The initialization population is composed of popsize-2 individuals generated by rule3.


### 3.1.2. Standard WOA

Wang et al. [17] designed a whale optimization algorithm (WOA) by studying the foraging behavior of humpback whales. WOA mainly includes three search mechanisms: encircling search mechanism, spiral search mechanism, and random search mechanism.Encircling search mechanism: when finding prey (solution), individual whales communicate with each other in the whole population so that the individual whales in the population keep approaching the prey. The mathematical model of encircling search mechanism is shown in formula (17)∼(18):(17)D=CXgen∗−Xigen,(18)Xigen+1=Xgen∗−A⋅D,where gen is the current number of iterations, Xigen is the i individual in the gen generation, Xgen∗ is the position of prey in the gen generation, A and C are coefficient vectors, a is the control parameter, t is the current number of iterations, each iteration updates the population, and Tmax is the maximum number of iterations, with the increase of the number of iterations.Spiral search mechanism: after calculating the distance from the prey, the individual whale takes the prey as the center and constantly approaches it in a spiral form. The mathematical model of spiral search mechanism is shown as(19)Xigen+1=D⋅eblcos2πl+Xigen,where b is the parameter controlling the shape of the spiral, l is the random value of [−1,1], and γ is the search factor.In equation (20), the probability p is introduced to choose whether to update the position information of individual whale by using the encircling search mechanism or the spiral search mechanism:(20)Xigen+1=Xigen−A⋅D,p<γ,D⋅eblcos2πl+Xgen∗,p≥γ.Random search mechanism: when the value ofA is greater than 1, it will search outside the bounding circle. When we cannot get the effective information of the prey, we will choose a random way to search the prey to update the position information of the population, which can effectively enhance the diversity of the population. Then, the mathematical model of random search mechanism is as follows:(21)Xigen+1=Xgen′−A⋅D,(22)D=CXgen′−Xigen,where Xgen′ is a randomly selected individual in the gen generation population.


### 3.1.3. Block Structure and Three-Dimensional Matrix

For any solution of 2L-OVRPTW, the block structure refers to the order of two consecutive adjacent customers in the solution. In the process of updating, the solution set shows certain similarities in the structure of the solution (block structure and its location information). The three-dimensional matrix proposed in the second stage of global search can effectively learn and accumulate the block structure. The location information of the solution set is obtained in the first stage so as to improve the overall effect of the algorithm. Because the same block structure in different positions will have a great impact on the target value of the solution, the first dimension of the designed three-dimensional matrix is used to record the location information of the block structure, and the other two dimensions represent two customers in the block structure.We defineBest_Pgen=Best_πgen,1,Best_πgen,2,⋯Best_πgen,g,⋯,Best_πgen,n, where Best_Pgen is the gen generation population, and Best_πgen,w=Best_π1gen,w,Best_π2gen,w,⋯Best_πggen,w,⋯,Best_πngen,w, where Best_πgen,w is the w individual and Best_πggen,w is the g customer in the w individual of gen population. Let TMl×n×ngen represent the three-dimensional matrix based on block structure and its location information in gen generation population, where TMl×n×ngenp,x,y represents the value of block structure x,y in p position and PMl×n×ngenp is the two-dimensional probability matrix normalized at p position. The specific definition is as follows:(23)ONE_TMl×n×ngen,wp,x,y1,x=Best_πpgen,w,y=Best_πp+1gen,w,0,else,p=1,2,⋯,lx,y=1,2,⋯,n,(24)TMl×n×ngenp,x,y=∑w=1popsizeONE_TMl×n×ngen,wp,x,y,p=1,2,⋯,lx,y=1,2,⋯,n,(25)TMl×n×ngenp,x=TMl×n×ngenp,x,1,TMl×n×ngenp,x,2,⋯,TMl×n×ngenp,x,n,p=1,2,⋯,lx=1,2,⋯,n,(26)TMl×n×ngenp=TMl×n×ngenp,x,1,⋯,TMl×n×ngenp,x,n=TMl×n×ngenp,1,1,⋯,TMl×n×ngenp,1,n⋮TMl×n×ngenp,n,1,⋯,TMl×n×ngenp,n,n,(27)TMl×n×ngen=TMl×n×ngen1,TMl×n×ngen2,⋯,TMl×n×ngenl,(28)PMl×n×ngenp=PMl×n×ngenp,x,1,⋯,PMl×n×ngenp,x,n=PMl×n×ngenp,1,1,⋯,PMl×n×ngenp,1,n⋮PMl×n×ngenp,n,1,⋯,PMl×n×ngenp,n,n,(29)∑x=1n∑y=1nPMl×n×ngenp,x,y=1,p∈1,2,⋯,l,where ONE_TMl×n×ngen,wp,x,y is the number of times that customer x=Best_πpgen,w and customer y=Best_πp+1gen,w appears in p position of individual w in the gen generation, that is, the number of times that block structure x,y appears in p position of individual w. For example, for individual Best_πgen,1=1,3,2,5,4,6 and individual Best_πgen,2=2,5,3,6,1,4, let p=2; then, the client block structure of individual 1 in position 2 is [3,2], and the client block structure of individual 2 in position 2 is [5,3]. Therefore, ONE_TMl×n×ngen,12,3,2=1, ONE_TMl×n×ngen,22,3,2=0, ONE_TMl×n×ngen,12,5,3=0, and ONE_TMl×n×ngen,22,5,3=1. As shown in Figure 3, the block structure information of each p position is recorded in TMl×n×ngenp.Figure 3
Three-dimensional matrix based on block structure.In the first stage of global search, the improved WOA is used to search in the whole solution space to obtain high quality. In the second stage, a three-dimensional matrix TM is constructed based on the high-quality solution obtained in the first stage. When the population is updated by using the three-dimensional matrix, the structure and location information of customer blocks in the current high-quality solution is learned and accumulated after each update and used to generate new individuals. The specific update process is as follows:Step 1: construct a three-dimensional probability matrix TMl×n×ngen by using the high-quality solution obtained by WOA and construct a list cos_sequence of remaining optional customers and a list tabu of selected customers. Perform step 2.Step 2: let p=1, and randomly generate a random number a,a∈0,1 and execute step 3.Step 3: select the matrix of position p, set the occurrence times of block structure containing customers in tabu to 0, and get the probability matrix of position p after TMl×n×ngenp normalization as PMl×n×ngenp; sum=∑x=1n∑y=1nPMl×n×ngenp,x,y; if sum>0, perform step 4; otherwise, perform step 5.Step 4: if there is a selectable block structure, the roulette method is used to select the customer block structure for PMl×n×ngenp. Let p=p+2; perform step 6.Step 5: if there is no optional block structure in the current position, randomly select a customer in cos_sequence and place it in this position; let p=p+1, and then, go to step 6.Step 6: delete the selected customers from cos_sequence and add the selected customers to tabu.Step 7: if p>num_customer, output the generated new individual tabu; if p=num_customer, return to step 5; otherwise, return to step 3.The pseudocode of the LWOA algorithm is as follows.Algorithm 1: Learning whale optimization algorithm.
(1)
Input:(2)
Let the number of customers beN.Mn represents the total number of items required by the customer n.gen be the current number of iterations; the maximum number of iterations is genmax.πigen is the individual i of generation gen, and the population size is popsize.γ is the search factor.The update times of three-dimensional matrix ismax_l.The vehicle number isk,k∈K.Ing indicates the item g of the customer n(3)
LetA=2ar1−a,C=2r2,r1,r2∈random0,1, a=2−2∗gen/genmax, b=1,l∈random−1,1, p∈random0,1(4)
begin:(5)
Population initialization: using rule1 and rule2 in Section 3.1.1 to generate an individual and using rule3 to generate popsize−2 individuals.(6)
for each gen∈genmaxdo(7)
for each i∈popsizedo(8)
ifA>1then(9)
updateπigen by equation (21)(10)
else(11)
updateπigen by equation (20)(12)
end if(13)
end for(14)
end for(15)
for each gen∈max_ldo(16)
Choosing the best individualsπgen to construct the three-dimensional matrix(17)
Update the population by the three-dimensional matrix(18)
end for(19)
Choosing the best individual asπbest∗(20)
Loading initialization: Initialize the loading sequence of all customers using rules 1 and 2 in 3.4.1.(21)
for each k∈πbest∗do(22)
for each n∈kdo(23)
for each Ing∈ndo(24)
Loading the customer’s goods into the carriage through the loading strategy in Section3.4.2(25)
end for(26)
end for(27)
end for(28)
Output:(29)
returnπbest∗The values range ofγ, popsize, and max_l are analyzed in detail in Section 4.2.


## 3.1.1. Population initialization

In order to ensure the diversity and dispersion of the solutions, two individuals are generated byrule1 and rule2, and then, the rest individuals are generated by rule3. Rule1 considers the nearest neighbor principle, rule2 considers the highest satisfaction principle, and rule3 adopts the quasi-opposition strategy. Placing all unserved customers in set Z, the specific rules are as follows:Rule1 is the nearest neighbor principle.After the vehicle departs from the central depot, the distance between the central depot and each customer inZ is calculated, and the nearest customeri (if there are multiple customers, the customer with the highest satisfaction) is selected for service. The customeri from set Z was removed and the distances between customer i and other customers inZ are calculated. After selecting the nearest customer j, if the constraints are still met when continuing to serve the customer j, it is necessary to remove the customer j from set Z and calculate the distance between the customer and other customers in Z. Otherwise, the vehicle will give up serving the customer and dispatch the vehicle from the central depot to serve the customers in Z again. This cycle will last until all customers inZ are served.Rule2 is the maximum satisfaction principle.After vehicle starts from the central warehouse, we calculate the arrival time of each customer inZ from the central depot and select the customer i meeting the time window constraint to serve. If there are multiple customers whose arrival time is within the time window, the nearest customer i is selected to serve. Move customer i out of set Z and calculate the time from customer i to other customers in Z. Select customer j that meets the time window constraint. If customer j still meets the constraint, it is necessary to move customer j out of set Z and calculate the arrival times from customer Z to other customers in Z. Otherwise, the vehicle will give up serving the customer and send the vehicle from the central depot to serve the customers in Z again until all the customers in Z are served.Rule3 is quasi-antagonistic strategy.Let the population size bepopsize and use the quasi-opposition strategy to generate popsize-2 individuals with n customers (that is, the search space is n). Let πigen=πigen1,πigen2,…,πigeng,…,πigenSπigeng∈V,n+1≤S≤2n denote thei individual (distribution scheme) in the gen population. Let Pigen=Pigen1,Pigen2,…,Pigeng,…,PigennPigeng∈V′ be the customer-based discrete sort in individual πigen (i.e., the customer sequence without central depot), and let Xigen=Xigen1,Xigen2,…,Xigeng,…,XigennXigeng∈a,b be the real number sequence corresponding to the customer’s discrete sort Pigen. As shown in equation (16), the quasi-antagonistic strategy is to obtain the opposite inverse solution (quasi-antagonistic solution) Xjgen=Xjgen1,Xjgen2,…,Xjgeng,…,XjgennXjgeng∈a,b of the gen generation individual i through the calculation of equation (16) in the real number field a,b. Therefore, the population with quasi-opposite solution can be obtained when the population is randomly generated. popsize-2 individuals are selected from the two populations as the next generation population:(16)Xjgeng=randX¯,Xigeng,Xigeng>X¯,randXigeng,X¯,Xigeng≤X¯,X¯=a+b2.The population obtained byXigen through rule3 is a real number sequence over a continuous field. Therefore, this paper uses LOV (large-order value) rule to transform the individual real number sequence Xigen into discrete customer sequence Pigen.The specific steps of LOV rules are as follows:Step 1: we select the real number sequence Xigen=Xigen1,Xigen2,…,Xigeng,…,Xigenn, arrange Xigen in the descending order to get Xigen′, find the corresponding sequence number Xigen in Xigen′g, and place it in φigeng to get the intermediate sequence φigenStep 2: for discrete customers, the values in Pigen are obtained successively by the formula Pigenφigeng=gAs shown in Table1, an example of LOV application is given. If the real number sequence of an individual is Xigen=2.30,1.59,2.17,1.78,2.06,1.89,0.88, the real number sequence after descending is Xigen′=2.30,2.17,2.06,1.89,1.78,1.59,0.88. Through the above steps, the intermediate sequence φigen=1,6,2,5,3,4,7 and the customer’s discretization sequence Pigen=1,3,5,6,4,2,7 are obtained.Table 1
Solution representation of LOV.
Dimension (g)1234567Xigeng2.301.592.171.782.061.890.88φigeng1625347Pigeng1356427By calling the LOV rule, we can realize the transformation from real sequence to discrete sequence. Similarly, if the customer’s discrete sequence changes, we can adjust it by using the reverse largest order value (RLOV) rule. Therefore, the RLOV rule is used to real number the discrete customer sequence obtained byrule1 and rule2 to get the corresponding real number solution. The initialization population is composed of popsize-2 individuals generated by rule3.


## 3.1.2. Standard WOA

Wang et al. [17] designed a whale optimization algorithm (WOA) by studying the foraging behavior of humpback whales. WOA mainly includes three search mechanisms: encircling search mechanism, spiral search mechanism, and random search mechanism.Encircling search mechanism: when finding prey (solution), individual whales communicate with each other in the whole population so that the individual whales in the population keep approaching the prey. The mathematical model of encircling search mechanism is shown in formula (17)∼(18):(17)D=CXgen∗−Xigen,(18)Xigen+1=Xgen∗−A⋅D,where gen is the current number of iterations, Xigen is the i individual in the gen generation, Xgen∗ is the position of prey in the gen generation, A and C are coefficient vectors, a is the control parameter, t is the current number of iterations, each iteration updates the population, and Tmax is the maximum number of iterations, with the increase of the number of iterations.Spiral search mechanism: after calculating the distance from the prey, the individual whale takes the prey as the center and constantly approaches it in a spiral form. The mathematical model of spiral search mechanism is shown as(19)Xigen+1=D⋅eblcos2πl+Xigen,where b is the parameter controlling the shape of the spiral, l is the random value of [−1,1], and γ is the search factor.In equation (20), the probability p is introduced to choose whether to update the position information of individual whale by using the encircling search mechanism or the spiral search mechanism:(20)Xigen+1=Xigen−A⋅D,p<γ,D⋅eblcos2πl+Xgen∗,p≥γ.Random search mechanism: when the value ofA is greater than 1, it will search outside the bounding circle. When we cannot get the effective information of the prey, we will choose a random way to search the prey to update the position information of the population, which can effectively enhance the diversity of the population. Then, the mathematical model of random search mechanism is as follows:(21)Xigen+1=Xgen′−A⋅D,(22)D=CXgen′−Xigen,where Xgen′ is a randomly selected individual in the gen generation population.


## 3.1.3. Block Structure and Three-Dimensional Matrix

For any solution of 2L-OVRPTW, the block structure refers to the order of two consecutive adjacent customers in the solution. In the process of updating, the solution set shows certain similarities in the structure of the solution (block structure and its location information). The three-dimensional matrix proposed in the second stage of global search can effectively learn and accumulate the block structure. The location information of the solution set is obtained in the first stage so as to improve the overall effect of the algorithm. Because the same block structure in different positions will have a great impact on the target value of the solution, the first dimension of the designed three-dimensional matrix is used to record the location information of the block structure, and the other two dimensions represent two customers in the block structure.We defineBest_Pgen=Best_πgen,1,Best_πgen,2,⋯Best_πgen,g,⋯,Best_πgen,n, where Best_Pgen is the gen generation population, and Best_πgen,w=Best_π1gen,w,Best_π2gen,w,⋯Best_πggen,w,⋯,Best_πngen,w, where Best_πgen,w is the w individual and Best_πggen,w is the g customer in the w individual of gen population. Let TMl×n×ngen represent the three-dimensional matrix based on block structure and its location information in gen generation population, where TMl×n×ngenp,x,y represents the value of block structure x,y in p position and PMl×n×ngenp is the two-dimensional probability matrix normalized at p position. The specific definition is as follows:(23)ONE_TMl×n×ngen,wp,x,y1,x=Best_πpgen,w,y=Best_πp+1gen,w,0,else,p=1,2,⋯,lx,y=1,2,⋯,n,(24)TMl×n×ngenp,x,y=∑w=1popsizeONE_TMl×n×ngen,wp,x,y,p=1,2,⋯,lx,y=1,2,⋯,n,(25)TMl×n×ngenp,x=TMl×n×ngenp,x,1,TMl×n×ngenp,x,2,⋯,TMl×n×ngenp,x,n,p=1,2,⋯,lx=1,2,⋯,n,(26)TMl×n×ngenp=TMl×n×ngenp,x,1,⋯,TMl×n×ngenp,x,n=TMl×n×ngenp,1,1,⋯,TMl×n×ngenp,1,n⋮TMl×n×ngenp,n,1,⋯,TMl×n×ngenp,n,n,(27)TMl×n×ngen=TMl×n×ngen1,TMl×n×ngen2,⋯,TMl×n×ngenl,(28)PMl×n×ngenp=PMl×n×ngenp,x,1,⋯,PMl×n×ngenp,x,n=PMl×n×ngenp,1,1,⋯,PMl×n×ngenp,1,n⋮PMl×n×ngenp,n,1,⋯,PMl×n×ngenp,n,n,(29)∑x=1n∑y=1nPMl×n×ngenp,x,y=1,p∈1,2,⋯,l,where ONE_TMl×n×ngen,wp,x,y is the number of times that customer x=Best_πpgen,w and customer y=Best_πp+1gen,w appears in p position of individual w in the gen generation, that is, the number of times that block structure x,y appears in p position of individual w. For example, for individual Best_πgen,1=1,3,2,5,4,6 and individual Best_πgen,2=2,5,3,6,1,4, let p=2; then, the client block structure of individual 1 in position 2 is [3,2], and the client block structure of individual 2 in position 2 is [5,3]. Therefore, ONE_TMl×n×ngen,12,3,2=1, ONE_TMl×n×ngen,22,3,2=0, ONE_TMl×n×ngen,12,5,3=0, and ONE_TMl×n×ngen,22,5,3=1. As shown in Figure 3, the block structure information of each p position is recorded in TMl×n×ngenp.Figure 3
Three-dimensional matrix based on block structure.In the first stage of global search, the improved WOA is used to search in the whole solution space to obtain high quality. In the second stage, a three-dimensional matrix TM is constructed based on the high-quality solution obtained in the first stage. When the population is updated by using the three-dimensional matrix, the structure and location information of customer blocks in the current high-quality solution is learned and accumulated after each update and used to generate new individuals. The specific update process is as follows:Step 1: construct a three-dimensional probability matrix TMl×n×ngen by using the high-quality solution obtained by WOA and construct a list cos_sequence of remaining optional customers and a list tabu of selected customers. Perform step 2.Step 2: let p=1, and randomly generate a random number a,a∈0,1 and execute step 3.Step 3: select the matrix of position p, set the occurrence times of block structure containing customers in tabu to 0, and get the probability matrix of position p after TMl×n×ngenp normalization as PMl×n×ngenp; sum=∑x=1n∑y=1nPMl×n×ngenp,x,y; if sum>0, perform step 4; otherwise, perform step 5.Step 4: if there is a selectable block structure, the roulette method is used to select the customer block structure for PMl×n×ngenp. Let p=p+2; perform step 6.Step 5: if there is no optional block structure in the current position, randomly select a customer in cos_sequence and place it in this position; let p=p+1, and then, go to step 6.Step 6: delete the selected customers from cos_sequence and add the selected customers to tabu.Step 7: if p>num_customer, output the generated new individual tabu; if p=num_customer, return to step 5; otherwise, return to step 3.The pseudocode of the LWOA algorithm is as follows.Algorithm 1: Learning whale optimization algorithm.
(1)
Input:(2)
Let the number of customers beN.Mn represents the total number of items required by the customer n.gen be the current number of iterations; the maximum number of iterations is genmax.πigen is the individual i of generation gen, and the population size is popsize.γ is the search factor.The update times of three-dimensional matrix ismax_l.The vehicle number isk,k∈K.Ing indicates the item g of the customer n(3)
LetA=2ar1−a,C=2r2,r1,r2∈random0,1, a=2−2∗gen/genmax, b=1,l∈random−1,1, p∈random0,1(4)
begin:(5)
Population initialization: using rule1 and rule2 in Section 3.1.1 to generate an individual and using rule3 to generate popsize−2 individuals.(6)
for each gen∈genmaxdo(7)
for each i∈popsizedo(8)
ifA>1then(9)
updateπigen by equation (21)(10)
else(11)
updateπigen by equation (20)(12)
end if(13)
end for(14)
end for(15)
for each gen∈max_ldo(16)
Choosing the best individualsπgen to construct the three-dimensional matrix(17)
Update the population by the three-dimensional matrix(18)
end for(19)
Choosing the best individual asπbest∗(20)
Loading initialization: Initialize the loading sequence of all customers using rules 1 and 2 in 3.4.1.(21)
for each k∈πbest∗do(22)
for each n∈kdo(23)
for each Ing∈ndo(24)
Loading the customer’s goods into the carriage through the loading strategy in Section3.4.2(25)
end for(26)
end for(27)
end for(28)
Output:(29)
returnπbest∗The values range ofγ, popsize, and max_l are analyzed in detail in Section 4.2.


## 3.2. Local Search

According to the characteristics of the 2L-OVRPTW problem, the solution space is huge, and the high-quality solutions are usually distributed in each local region of the solution space. Therefore, it is necessary to search the high-quality solutions from the global search; VNS includes three neighborhood operations: interchange operation, insert operation, and swap operation. Interchange operation is used between vehicles, and insert operation and swap operation are used inside vehicles. In this way, local search can be carried out more carefully to update the current optimal solution. The variable neighborhood local search execution pseudocode is as follows.Algorithm 2: Variable neighborhood local search.
(1)
Input:(2)
Choose the best individualπi and convert it to the permutation π0 according to the LOV rule. Perform a neighborhood operation on π0 to get π(3)
begin:(4)
Set loop = 1, max-loop = 40,n = 1, max-n = 20(5)
whilen < max-ndo(6)
Random selectm-k, m-k∈ {1, 2, 3}(7)
Randomly selectk and l, where k≠ l(8)
ifm-k = 1, thenπ′=insertπ,k,l(9)
else ifm-k = 2, thenπ′=swapπ,k,l(10)
else if m-k = 3, thenπ′=inverseπ,k,l(11)
end if(12)
end if(13)
end if(14)
ifπ′≺π, thenπ=π′(15)
elsen = n + 1(16)
end if(17)
end while(18)
π′=interchangeπ,k,l(19)
ifπ′≺π,then π=π′(20)
end if(21)
until loop = max-loop(22)
ifπ≺π0, then π0=π(23)
end if(24)
Output:(25)
Convertπ0 back toπi according to the RLOV rule


## 3.3. Complexity Analysis of the LWOA

Let the population size be popsize, the number of running iterations of global search be genglobal, and the number of iterations of each generation of local search be genlocal. The number of times that local search actually uses neighborhood operators in each generation of search is Z(Z = 20 ∗ 40). The algorithm complexity of LWOA TLWOA is composed of global search complexity and local search complexity, respectively. In each generation, the global search complexity consists of population initialization complexity O(popsize ∗ N, where N is the number of customers), discretization complexity (2 ∗ popsize ∗ N!), three-dimensional probability matrix update complexity O(N2 ∗ (N − 1)), complexity O(popsize∗ (1 + N2 ∗ (N − 1))) of generating new population based on three-dimensional probability matrix, and population evaluation complexity (popsize ∗ N); the result is O(popsize ∗ (1 + N2 ∗ (N − 1))). The complexity of local search is O(genlocal ∗ Z ∗ (N)) (the worst-case scenario).Through the above analysis,TLWOA = O(genglobal ∗ (popsize ∗ (1 + N2 ∗ (N − 1)) + genlocal ∗ Z ∗ (N))).


## 3.4. Loading Strategy

For the loading problem based on the principle of “first in, last out” in the customer service order, this paper adopts the loading strategy of the skyline filling method to check every feasible solution generated. Only when the loading constraint is satisfied can the solution be retained. This section mainly includes initialization of loading order and loading mechanisms.


### 3.4.1. Initialize Loading Sequence

Each time we get a solution, we need to make sure that all the customers’ goods can be loaded into the car in the reverse order of the customers’ order. However, the loading order of multiple goods of each customer is not unique. Therefore, this paper proposes two rules for the initial loading order of customers’ goods.Rule1: sort items according to their lengthRule2: sort items by area from small to largeFor example, if the vehicle service sequence is customer 1, customer 2, customer 3, and customer 4, customer 4 of the final service usesrule1 for initial loading and the other three customers use rule2 for initial loading.


### 3.4.2. Loading Method

The loading process of articles is based on the skyline. Initially, the horizontal skyline and the vertical skylines are the width of the carriage and the length of both sides of the vehicle, respectively. The loading position of articles is determined by the score of articles placed on each section of the horizontal skyline. (1) When the width is equal to the horizontal skyline, if the length is equal to either side of the adjacent vertical skyline, the score is 5. If the length is lower than the longitudinal skyline of the lowest adjacent side, the score is 4. If the length is between the two sides, the score is 3. Otherwise, the score is 2. (2) If the width is less than the transverse skyline, the score is 4 if the length is equal to either side of the longitudinal skyline of the two adjacent sides. If the length is lower than the longitudinal skyline of the lowest adjacent side, the score is 3. Otherwise, the score is 2. (3) If the width is greater than the transverse skyline, the score is 1. If all the horizontal skylines are less than the width of the object, it is necessary to fill the skyline before loading the object. In addition, when placing the object, it is also necessary to consider whether it exceeds the length and width limits of the vehicle. Loading the object into the vehicle and updating the skyline at the position with the highest score.


## 3.4.1. Initialize Loading Sequence

Each time we get a solution, we need to make sure that all the customers’ goods can be loaded into the car in the reverse order of the customers’ order. However, the loading order of multiple goods of each customer is not unique. Therefore, this paper proposes two rules for the initial loading order of customers’ goods.Rule1: sort items according to their lengthRule2: sort items by area from small to largeFor example, if the vehicle service sequence is customer 1, customer 2, customer 3, and customer 4, customer 4 of the final service usesrule1 for initial loading and the other three customers use rule2 for initial loading.


## 3.4.2. Loading Method

The loading process of articles is based on the skyline. Initially, the horizontal skyline and the vertical skylines are the width of the carriage and the length of both sides of the vehicle, respectively. The loading position of articles is determined by the score of articles placed on each section of the horizontal skyline. (1) When the width is equal to the horizontal skyline, if the length is equal to either side of the adjacent vertical skyline, the score is 5. If the length is lower than the longitudinal skyline of the lowest adjacent side, the score is 4. If the length is between the two sides, the score is 3. Otherwise, the score is 2. (2) If the width is less than the transverse skyline, the score is 4 if the length is equal to either side of the longitudinal skyline of the two adjacent sides. If the length is lower than the longitudinal skyline of the lowest adjacent side, the score is 3. Otherwise, the score is 2. (3) If the width is greater than the transverse skyline, the score is 1. If all the horizontal skylines are less than the width of the object, it is necessary to fill the skyline before loading the object. In addition, when placing the object, it is also necessary to consider whether it exceeds the length and width limits of the vehicle. Loading the object into the vehicle and updating the skyline at the position with the highest score.


## 3.5. Flowchart of LWOA

This section gives the flowchart of LWOA. In the whole operation process, the individual needs to be boxed after each update of the individual, and the individual meeting the boxing constraints will be retained. The overall schematic diagram of the LWOA algorithm and the detailed flowchart of local search are shown in Figures4 and 5.Figure 4
LWOA flowchart.Figure 5
Local search flowchart.


## 4. Simulation Experiment and Comparative Analysis

The main content of this section includes the design of customer goods, performance index and parameter setting, and simulation experiment comparison and analysis. All experiments are run on Windows 10 platform, which has 3.0 GHz CPU and 8 GB RAM, and are on a single thread. The algorithm in this paper is implemented by Python 3.7.


### 4.1. Experimental Setup

Since there is no standard testing benchmark suitable for 2L-OVRPTW, this paper uses the dataset and generates customer items according to [24], as shown in Table 2. In the homogeneous fleet, the length and width (L, W) of the car are (40, 20). Category 1 corresponds to the pure CVRP instance, and no boxing constraint is imposed on it. Therefore, this paper uses categories 2–5 to generate rectangular items needed by customers. As shown in Table 2, select category 4 for the dataset with less than 25 customers, that is, randomly select the number of items required by customers in [1, 4] and then assign one of three possible shapes to each item with the same probability. The three shapes are vertical (relative height is greater than relative width), uniform (relative height and width are generated in the same interval), or horizontal (relative height is less than relative width). Similarly, for datasets with more than 25 customers, category 3 is used.Table 2
Item generation category.
CategorymiVerticalRegularLevelRelatively longRelatively wideRelatively longRelatively wideRelatively longRelatively wide2[1, 2][0.4L, 0.9L][0.1W, 0.2W][0.2L, 0.5L][0.2W, 0.5W][0.1L, 0.2L][0.4W, 0.9W]3[1, 3][0.3L, 0.8L][0.1W, 0.2W][0.2L, 0.4L][0.2W, 0.4W][0.1L, 0.2L][0.3W, 0.8W]4[1, 4][0.2L, 0.7L][0.1W, 0.2W][0.1L, 0.4L][0.1W, 0.4W][0.1L, 0.2L][0.2W, 0.7W]5[1, 5][0.1L, 0.6L][0.1W, 0.2W][0.1L, 0.3L][0.1W, 0.3W][0.1L, 0.2L][0.1W, 0.6W]


### 4.2. Parameter Setting

The parameters involved in LWOA include population size, search factor, and learning times of the three-dimensional matrix. This paper selects problem C201 with 25 customers for experimental analysis and uses design of experiment (DOE) to explore the influence of parameters on LWOA performance. Four level values are selected for each parameter, as shown in Table3, and the orthogonal test table of the L16 (43) scale is established. In each group of parameters, the algorithm runs independently 20 times, intending to minimize the mileage. The experimental results are as follows. The experimental results are as follows.Table 3
Parameter level.
ParameterLevel1234popsize30507090γ0.20.40.60.8max_l20406080It can be seen from Tables3–5 that different parameter combinations have a great influence on the algorithm. Therefore, to maximize the algorithm’s performance, this paper sets the parameter combination based on the experimental analysis. That is, the later experimental part will also be based on this parameter combination.Table 4
Results of orthogonal experiment.
ParameterLevelAVEpopsizeγmax_l1143395.9882134404.8563122412.8294111411.8335241379.2916232393.7957224383.1898213382.6939341387.11210333374.38411322359.40912314379.46413443377.91614432369.74615421380.63216414376.294Table 5
Average response value of parameters.
Levelpopsizeγmax_l1406.376390.176389.7172384.742384.015383.9453375.092385.695382.7454377.569385.077387.763Range31.2846.1616.972Grade123


### 4.3. Comparison and Analysis of Simulation Experiments

In order to verify the effectiveness of LWOA in solving 2L-OVRPTW, this paper compares LWOA with WOA, GA, and EDA. The objective of the optimization is to minimize the driving mileage for 20 independent experiments. The three performance indexes of optimal value (BST), the worst (WST) and the average value (AVG) are selected for comparative analysis, and the optimal results are roughened. The experimental results are shown in Table6. It can be seen from Table 6 that the three indexes of LWOA in most cases are significantly better than other comparison algorithms. This shows that the proposed three-dimensional matrix can effectively learn and accumulate block structure information based on customer order relationship, improve the overall search effect of the algorithm, and verify the effectiveness of LWOA in solving 2L-OVRPTW.Table 6
Comparison of LWOA, WOA, GA, and EDA.
ProblemLWOAWOAGAEDABSTWSTAVGBSTWSTAVGBSTWSTAVGBSTWSTAVGC101242.01311.58279.80245.30327.05288.46259.00333.53294.73243.28323.21288.33C102247.62317.52284.29252.38333.87295.02261.09321.98292.14251.17327.60287.46C103224.83315.95279.87252.77326.05283.24252.57354.41290.81246.63341.71286.53C104245.98308.16280.61268.13338.33296.82258.78314.58299.73245.44330.65292.10C105242.79301.49280.94234.54311.92277.52254.55333.01289.48253.74326.57283.90C106242.15317.8283.21266.43345.98297.04230.55332.64288.16254.22329.45295.28C107248.19306.76279.81261.62310.61283.85248.30316.79282.72255.93328.17294.70C108232.96321.57285.42260.09329.08292.76241.78326.56285.77242.38314.09279.75C109224.48318.68279.85253.77338.80287.91236.75330.6294.60247.37339.38288.42C201328.36406.8373.00357.83420.27389.59343.67413.31384.74342.78420.57382.82C202348.97410.72385.56352.52420.19389.23332.24425.89386.71353.85446.74392.55C203353.11425.83387.58358.11435.19395.75356.26451.24388.47357.60428.34383.84C204337.79415.66379.23354.86431.27384.55294.66448.28384.45337.60456.18399.28C205347.33416.55382.22361.29417.95394.90317.31429.01383.39356.92441.02387.02C206337.47409.87388.76364.22423.06395.74365.30433.86391.15363.18436.91395.05C207336.32418.19382.28359.74420.08388.47340.18422.48385.61337.76431.17391.35C208339.42427.92376.57346.28425.14381.89340.07446.34386.12340.42438.14394.69For the objective function values (BST, WST, and AVG) as test samples, nonparametric paired sample test is carried out on the operation results of each algorithm with 95% confidence. The results are shown in Table7.Table 7
Results of statistical analysis.
Paired testThe value ofPBSTWSTAVGLWOA and WOA0.0008460.0005030.000420LWOA and GA0.2273850.0002930.000293LWOA and EDA0.0005030.0004210.001184


## 4.1. Experimental Setup

Since there is no standard testing benchmark suitable for 2L-OVRPTW, this paper uses the dataset and generates customer items according to [24], as shown in Table 2. In the homogeneous fleet, the length and width (L, W) of the car are (40, 20). Category 1 corresponds to the pure CVRP instance, and no boxing constraint is imposed on it. Therefore, this paper uses categories 2–5 to generate rectangular items needed by customers. As shown in Table 2, select category 4 for the dataset with less than 25 customers, that is, randomly select the number of items required by customers in [1, 4] and then assign one of three possible shapes to each item with the same probability. The three shapes are vertical (relative height is greater than relative width), uniform (relative height and width are generated in the same interval), or horizontal (relative height is less than relative width). Similarly, for datasets with more than 25 customers, category 3 is used.Table 2
Item generation category.
CategorymiVerticalRegularLevelRelatively longRelatively wideRelatively longRelatively wideRelatively longRelatively wide2[1, 2][0.4L, 0.9L][0.1W, 0.2W][0.2L, 0.5L][0.2W, 0.5W][0.1L, 0.2L][0.4W, 0.9W]3[1, 3][0.3L, 0.8L][0.1W, 0.2W][0.2L, 0.4L][0.2W, 0.4W][0.1L, 0.2L][0.3W, 0.8W]4[1, 4][0.2L, 0.7L][0.1W, 0.2W][0.1L, 0.4L][0.1W, 0.4W][0.1L, 0.2L][0.2W, 0.7W]5[1, 5][0.1L, 0.6L][0.1W, 0.2W][0.1L, 0.3L][0.1W, 0.3W][0.1L, 0.2L][0.1W, 0.6W]


## 4.2. Parameter Setting

The parameters involved in LWOA include population size, search factor, and learning times of the three-dimensional matrix. This paper selects problem C201 with 25 customers for experimental analysis and uses design of experiment (DOE) to explore the influence of parameters on LWOA performance. Four level values are selected for each parameter, as shown in Table3, and the orthogonal test table of the L16 (43) scale is established. In each group of parameters, the algorithm runs independently 20 times, intending to minimize the mileage. The experimental results are as follows. The experimental results are as follows.Table 3
Parameter level.
ParameterLevel1234popsize30507090γ0.20.40.60.8max_l20406080It can be seen from Tables3–5 that different parameter combinations have a great influence on the algorithm. Therefore, to maximize the algorithm’s performance, this paper sets the parameter combination based on the experimental analysis. That is, the later experimental part will also be based on this parameter combination.Table 4
Results of orthogonal experiment.
ParameterLevelAVEpopsizeγmax_l1143395.9882134404.8563122412.8294111411.8335241379.2916232393.7957224383.1898213382.6939341387.11210333374.38411322359.40912314379.46413443377.91614432369.74615421380.63216414376.294Table 5
Average response value of parameters.
Levelpopsizeγmax_l1406.376390.176389.7172384.742384.015383.9453375.092385.695382.7454377.569385.077387.763Range31.2846.1616.972Grade123


## 4.3. Comparison and Analysis of Simulation Experiments

In order to verify the effectiveness of LWOA in solving 2L-OVRPTW, this paper compares LWOA with WOA, GA, and EDA. The objective of the optimization is to minimize the driving mileage for 20 independent experiments. The three performance indexes of optimal value (BST), the worst (WST) and the average value (AVG) are selected for comparative analysis, and the optimal results are roughened. The experimental results are shown in Table6. It can be seen from Table 6 that the three indexes of LWOA in most cases are significantly better than other comparison algorithms. This shows that the proposed three-dimensional matrix can effectively learn and accumulate block structure information based on customer order relationship, improve the overall search effect of the algorithm, and verify the effectiveness of LWOA in solving 2L-OVRPTW.Table 6
Comparison of LWOA, WOA, GA, and EDA.
ProblemLWOAWOAGAEDABSTWSTAVGBSTWSTAVGBSTWSTAVGBSTWSTAVGC101242.01311.58279.80245.30327.05288.46259.00333.53294.73243.28323.21288.33C102247.62317.52284.29252.38333.87295.02261.09321.98292.14251.17327.60287.46C103224.83315.95279.87252.77326.05283.24252.57354.41290.81246.63341.71286.53C104245.98308.16280.61268.13338.33296.82258.78314.58299.73245.44330.65292.10C105242.79301.49280.94234.54311.92277.52254.55333.01289.48253.74326.57283.90C106242.15317.8283.21266.43345.98297.04230.55332.64288.16254.22329.45295.28C107248.19306.76279.81261.62310.61283.85248.30316.79282.72255.93328.17294.70C108232.96321.57285.42260.09329.08292.76241.78326.56285.77242.38314.09279.75C109224.48318.68279.85253.77338.80287.91236.75330.6294.60247.37339.38288.42C201328.36406.8373.00357.83420.27389.59343.67413.31384.74342.78420.57382.82C202348.97410.72385.56352.52420.19389.23332.24425.89386.71353.85446.74392.55C203353.11425.83387.58358.11435.19395.75356.26451.24388.47357.60428.34383.84C204337.79415.66379.23354.86431.27384.55294.66448.28384.45337.60456.18399.28C205347.33416.55382.22361.29417.95394.90317.31429.01383.39356.92441.02387.02C206337.47409.87388.76364.22423.06395.74365.30433.86391.15363.18436.91395.05C207336.32418.19382.28359.74420.08388.47340.18422.48385.61337.76431.17391.35C208339.42427.92376.57346.28425.14381.89340.07446.34386.12340.42438.14394.69For the objective function values (BST, WST, and AVG) as test samples, nonparametric paired sample test is carried out on the operation results of each algorithm with 95% confidence. The results are shown in Table7.Table 7
Results of statistical analysis.
Paired testThe value ofPBSTWSTAVGLWOA and WOA0.0008460.0005030.000420LWOA and GA0.2273850.0002930.000293LWOA and EDA0.0005030.0004210.001184


## 5. Conclusions

In this paper, a learning whale optimization algorithm (LWOA) combined with improved skyline filling algorithm (ISFA) is proposed to solve 2L-OVRPTW. The specific contributions are as follows. (1) Establish the 2L-OVRPTW model with the aim of minimizing the total distance. (2) In the route optimization stage, after the whale optimization algorithm is used to update the population, a three-dimensional matrix update strategy based on customer order relationship is designed to learn and accumulate the block structure and customer location information in the high-quality solution, which can effectively guide the population update. According to the characteristics of the problem, the variable neighborhood local search between and within vehicles is designed to carry out the detailed and efficient local search for the high-quality solution individuals after global search. (3) In the loading optimization stage, to improve the feasibility of the post solution of path optimization, ISFA is used to optimize the loading process, and the customer’s goods are effectively loaded into the carriage through the scoring mechanism based on the skyline.Future work will focus on the multiobjective problem with loading and time window constraints, refine the nature of the problem, and design an efficient algorithm to solve it so that the research is more in line with the actual distribution of logistics transportation.



---

*Source: 8016356-2021-12-26.xml*