Multi-Objective Material Logistics Planning with Discrete Split Deliveries Using a Hybrid NSGA-II Algorithm

Abstract

To schedule material supply intelligently and meet the production demand, studies concerning the material logistics planning problem are essential. In this paper, we consider the problem based on the scenario that more than one vehicle may visit each station in batches. The primary objective is to satisfy the demands in the time windows, followed by logistics planning with the minimum vehicles and travel time as the optimization objective. We construct a multi-objective mixed-integer programming model for the scenario of discrete material supply in workshops. First, we propose a hybrid heuristic algorithm combining NSGA-II and variable neighborhood search. This proposed algorithm combines the global search capability of NSGA-II and the strong local search capability, which can balance intensification and diversification well. Second, to maintain the diversity of population, we design the population diversity strategy and various neighborhood operators. We verify the effectiveness of the hybrid algorithm by comparing with other algorithms. To test the validity of the proposed problem, we have carried out research and application in a construction machinery enterprise.

1. Introduction

Production logistics is a variant application of logistics in the manufacturing system, which has a direct impact on the production cost and delivery time of products. Making reasonable material supply planning under the known material plan is an important part of planning and scheduling research. A reasonable and effective logistics system can not only improve the efficiency of distribution, but also support the continuous production and stability of the assembly line.

Materials are stored in the warehouse in batches. If the batch is the transportation unit, the unpacking and packaging processes of materials do not need extra manpower. Reasonable split deliveries by batch help derive the most use out of a vehicle’s loading capacity and reduces the number of vehicles and delivery time. Just-in-time delivery is an effective control method to avoid production interruption caused by material shortage.

Based on the premise of multiple batches of various materials, each station may be visited by more than one vehicle in batches. This kind of problem belongs to the Split Delivery Vehicle Routing Problem (SDVRP), which was first studied in 1989 [1]. The existing studies on SDVRP follow combinatorial splitting, based on the arbitrary splitting of node demands. However, in production logistics, most materials are distributed in batches [2,3]. Therefore, the batch can be considered as the minimum distribution unit that is not allowed to be split.

Based on the analysis of the mixed-flow assembly line, we explore the vehicle routing problem with time windows and discrete split deliveries by batch (VRPTWDB) with meeting the load constraint, time windows constraint, material demands, and other conditions. As far as we know, few studies consider the multi-objective optimization on SDVRP in the mixed-flow assembly line scenarios.

The primary objective of VRPTWDB is to meet the demands within the time constraint, followed by logistics planning with the minimum vehicles and total travel time as the optimization objective. Therefore, we construct a multi-objective mixed integer programming model for the scenario of the discrete material supply in workshops.

To solve the NP-hard problem, we propose the NSGA-II hybrid algorithm with variable neighborhood descent (VND) (HNSGA-II-VND). The proposed algorithm combines the global search capability of NSGA-II and the strong local search capability of VND, which can balance intensification and diversification. The proposed algorithm provides a series of Pareto optimal solutions with the minimum vehicles and total travel time. The proposed algorithm uses the Pareto optimization concept to solve all the objectives simultaneously. To maintain the diversity of population, we design the population diversity strategy in the hybrid algorithm. To reduce computing time, the adaptive expectation factors and various neighborhood operators are designed. The effectiveness of the hybrid algorithm is verified by comparison with other algorithms. To test the validity of the proposed problem, we have carried out application in a construction machinery enterprise.

The current research aims to solve this multi-objective problem with waiting time optimization by an efficient MILP model and the hybrid algorithm HNSGA-II-VND. The major contributions of this paper are as follows:

• We study the multi-objective VRPTWDB problem in the mixed-flow assembly line scenario. Based on the highly mixed demand environment, the relationship between station visit sequence and multiple optimization objectives (minimize the routes, transportation time and waiting time) is comprehensively considered.
• The decision optimization mathematical model is constructed to quantitatively represent the material logistics planning problem of highly mixed production.
• An effective multi-objective HNSGA-II-VND algorithm is proposed to solve the problem. The optimization potential of the solution is further explored by crossover operator, VND local search and the population diversity mechanism. Batch combination operation and batches random binding operation are designed.

The rest of the paper is organized as follows. Section 2 introduces the relevant existing studies. Section 3 constructs the MILP model. Section 4 elaborates the hybrid algorithm in detail. The analysis of results and the experimental evaluation are illustrated in Section 5 and Section 6. Section 7 presents conclusions and future directions.

2. Literature Review

2.1. The Distribution Logistics

The logistics system of the assembly workshop needs to “deliver the right materials in the right quantity to the right station at the right time under the right conditions” [4]. Vehicle scheduling is the key of the logistics system, which is generally simplified to scheduling strategy selection problem and vehicle routing problem (VRP) [5].

In the material transportation research of manufacturing systems, the optimization of vehicle scheduling strategies takes up a large research proportion. Considering the limitations of actual physical factors in the workshop, some scholars study the dynamic scheduling of vehicles to avoid collisions [6,7]. Many heuristic algorithms have been designed to solve this problem [8,9,10].

The distribution logistics is an important application of VRP in manufacturing systems. Planning routes reasonably is the key to realizing efficient transportation. The optimization goal of logistics planning is to arrange vehicles to complete the tasks with minimum costs and travel distances when the material demands plan is known [11]. Heuristic algorithms are designed to solve the logistics planning problem, which mainly focuses on the delivery time, delivery quantity, and transportation route [12,13,14].

2.2. Discrete Split Delivery Vehicle Routing Problem

In actual transportation, discrete splitting with the smallest non-detachable unit is common. Few research considers the SDVRP in discrete splitting. Nakao and Nagamochi [15] propose the concept of discrete-type split delivery and define it as Discrete Split Delivery Vehicle Routing Problem (DSDVRP). Salani and Vacca [16] formulate the DSDVRP with time window, assuming all feasible orders are known in advance.

Gulczynski et al. [17] consider SDVRP with minimum delivery amounts are allowed if a minimum fraction of demand is serviced by a vehicle. Many researchers continue to study the problem from the perspective of minimum delivery amounts [18,19,20].

Archetti et al. [21] allow the splitting of the demand of a customer only for different commodities. They present a branch-price-and-cut algorithm to solve the SDVRP. Chen et al. [3] propose a novel and efficient approach to solve the SDVRP using split strategy. Xia et al. [22] discuss the problem of discrete split only through backpacks from a green and low-carbon perspective. They use the hierarchical approach to construct a double-objective mathematical model. Xia and Fu study the DSDVRP from the perspective of backpack splitting and order splitting [23,24,25]. They improve tabu search with dynamic tabu lists and multi-neighborhood operators.

SDVRP with time windows means that each customer can only be accessed during a certain period of time. If the actual visit time of the vehicle to the customer is not within the required time period, it will be punished. Therefore, they can be divided into hard time window and soft time window accordingly.

The soft time window model is more flexible, which is beneficial to enhance the flexibility of the distribution system [23]. More feasible solutions can be obtained, and the range of solutions can be enlarged by soft time window [26,27,28].

In fact, the application of hard time windows of SDVRP is necessary in emergency scenarios [11,29,30]. The application of hard time window ensures the on-time delivery of vehicles without delaying the use of products [31,32]. Salani and Vacca [16] proposed the hard time window constraint of DSDVRP.

In view of the background of mixed flow assembly workshop production, only hard time window is appropriate. If the vehicle is earlier than the earliest acceptance time of the station, the vehicle must wait until that time for service. However, deliveries outside the time window can lead to assembly line shutdowns.

In the study of SDVRP with soft time window, vehicles are allowed to arrive at customers earlier and later, so the factor of waiting time is taken into account. However, the wait time and delay time are usually considered together. In the study of single objective problems, the minimization of vehicle cost and total traveling distance are generally considered together with waiting time and delay time, which are generally set with certain penalty values, respectively [23,30]. However, in multi-objective problems, the sum of waiting time and delay time is optimized as a single objective [27,28,33].

In the research of VRP with hard time window, most studies only consider the minimum cost and total traveling distances. However, waiting time is also a part of the total duration of vehicles. Reducing waiting time is also a way to reduce the total duration [31,32]. Based on the background of mixed-flow assembly workshop, materials need to closely abide by the time constraint, so we consider the optimization of waiting time in the proposed problem.

2.3. Multi-Objective SDVRP

Research into multi-objective SDVRP is becoming popular. Belfiore and Yoshizaki [34] designed a scatter search algorithm for heterogeneous fleet VRP with time windows and split deliveries to find good solutions, including swapping the same route, demand reallocation, route elimination and combination, and route addition. They consider the problem because it occurs in a major Brazilian retail group. Then, they [29] applied this method to solve the fleet size and mixed VRP with time windows and split deliveries. Their research extends to general instances for verification. However, the scatter search algorithm does not elaborate on the processing strategy of the two objectives. Chan et al. [26] study dynamic scheduling of oil tankers with the splitting of cargo at pickup and delivery locations. They developed a mathematical model aimed at minimizing total costs and maximizing service level. A modified multi-objective ant colony algorithm is proposed to solve the problem. Based on the concept that the Pareto optimal solution set optimizes two objectives simultaneously, Chiang and Hsu [11] used a multi-objective evolutionary algorithm combining enhanced crossover and mutation operators to solve VRPTW.

Based on the requirements of low-carbon logistics, Xia et al. [22] propose the logistics VRP with split deliveries by backpack. They used the hierarchical approach to construct a double-objective mathematical model to minimize the number of vehicles used and the travel time. Gupta et al. [35] study SDVRP with fuzzy time-distances characteristics to minimize fuel emissions. Fast comparisons are made based on fuzzy rules. A discrete fuzzy hybrid genetic algorithm is developed to solve multi-objective SDVRP. For solving the problem involved in classified waste disposal under an IoT environment, Cao B et al. [36] propose an improved SDVRP multi-objective model and a modified ant colony optimization (MACO) algorithm. Hasani Goodarzi, A [27] address a VRP with cross-docking that considers vehicle scheduling, splitting pick-up and delivery with time-windows at supplier and retailer locations, while optimizing two conflicting objectives. Shahabi-Shahmiri, R et al. [28] propose a multi-objective mixed-integer programming model for the scheduling and routing of heterogeneous vehicles carrying perishable goods across multiple cross-docking systems.

Research into multi-objective SDVRP is becoming popular. However, there are a few multi-objective optimizations for DSDVRP. Through multi-objective research, decision makers can make more favorable choices.

For DSDVRP problems without any batch combination process, each batch that can be split at each station is regarded as a single node. More specifically, the distance between nodes at the same station is regarded as 0 [22,23,24,37]. Then, the problem can be transformed into general VRPs and solved by existing algorithms [38].

However, Qiu. M et al. [25] are aware of such problems. To avoid the same customer being visited more than once, increasing the route distance compared to a normal visit, a batch combination operation is designed. A single batch as an operator results in complex and time-consuming calculations. The authors designed a neighborhood move object instead of customer or batch as a moving object by randomly combining batches. The split shift operator proposed by Nagy et al. [39] is very suitable for neighborhood exploration of this move object.

A comparative analysis of the existing literature is shown in

Table 1. Comparative analysis with existing literature.

Properties	Multi-Objective	Linear Objective Function	Min Vehicles	Min the Traveling Distances	Min Waiting Times	Discrete Split Deliveries	Split Deliveries	Hard Time Windows	Soft Time Windows
Nakao and Nagamochi [15]				√		√
Salani and Vacca [16]		√		√		√		√
Gulczynski et al. [17]		√		√			√
Archetti et al. [21]		√		√			√
Chen et al. [3]			√	√			√
Xia et al. [22]	√	√	√	√		√
Xia et al. [23]		√	√	√	√	√			√
Belfiore and Yoshizaki [29]		√	√	√			√	√
Chan et al. [26]	√	√		√			√		√
Chiang and Hsu [11]	√		√	√				√
Gupta et al. [35]	√			√			√
Cao B et al. [36]	√	√	√	√			√
Hasani Goodarzi, A [27]	√	√		√	√		√		√
Shahabi-Shahmiri, R et al. [28]	√	√		√	√		√		√
Qiu, M et al. [25]		√		√		√
Nagy et al. [39]		√		√			√
Xia et al. [37]	√	√	√	√		√
Lim, H et al. [38]		√		√			√
Belhaiza, S et al. [30]		√	√	√	√				√
Pan, J et al. [40]		√		√			√
Cordeau, J.F et al. [31]				√	√			√
Tricoire, F et al. [32]		√		√	√			√
Rocha, D et al. [41]	√	√		√
Ombuki, B et al. [42]	√		√	√				√
Xu, H et al. [33]	√		√	√	√				√
This paper	√	√	√	√	√	√		√

. The literature review shows that the SDVRP problem has been investigated widely and various heuristic algorithms have been designed to solve SDVRP. The VRPTWDB problem is a relevant problem in the industry, which has been gradually studied by academics in recent years. The type of problems addressed, although similar, include different constraints or objectives.

We study the VRPTWDB problem from a multi-objective perspective, based on the practical assembly workshop scene. Our work includes a new constraint on the hard time windows. This perspective has its own unique importance that no one has considered before. The study of this problem has important guiding significance for workshop logistics planning.

3. Mathematical Formulation

3.1. Problem Description

The mixed-flow assembly line is a very important part of the manufacturing industry. It can produce various products of the same category within a certain production cycle. On the mixed-flow assembly line, each workstation performs different tasks to assemble given various materials to produce multiple specific products.

According to the logic of delivery time and quantity, the planning layer summarizes the demands to generate the distribution planning. Therefore, when analyzing logistics optimization from the perspective of the logistics operation layer, the key is on the distribution logistics planning. The purpose is to ensure that each station receives the required materials within the time windows. More consideration should be given to various complicated time constraints in the production and distribution process to improve the reliability of distribution. Under the premise of not causing material shortage, the material of the same station can be relatively flexible, allowing it to be split and distributed. Reasonable split delivery by batch is beneficial to fully utilize the vehicle’s loading capacity and reduce the number of vehicles and delivery time.

The problem is stated as follows: the distribution center has K vehicles to dispatch, which cannot depart from the center at the same time. Each station has a demand time window and demands. Split deliveries by batch are allowed without interruptions to the production. The same station can be visited by multiple vehicles. The goal is to plan a reasonable transportation route and arrange the appropriate vehicle departure sequence to complete the delivery task to minimize the total cost and total delivery time.

Moreover, we make the following assumptions:

• The workshop roads network is an undirected network without capacity limitation.
• The vehicles will not stop due to failure or traffic jams during transportation.
• Assembly planning and material distribution lists are known. The materials are qualified and meet the requirements of homogeneity. Products will not be changed during the production process.
• A certain number of materials have already been stored in the line-side before starting production. The material is consumed evenly during production.
• Materials are converted to uniform equivalent according to weight and volume.
• The speed and rated load of the vehicles are the same.

3.2. The MILP Model

Let $\mathrm{G}=(V,E$) be a graph. Let$\mathrm{V}=\left\{0\right\}{{\displaystyle \cup }}^{}\mathrm{N}$ be the set of distribution center and workstations. Let$E=\left(\mathrm{V}\times V\right)$ be the set of arcs: there are travel times $t_{ij}$ associated with each arc. Let $\mathrm{k}\in K$ be the set of the vehicles. It also refers specifically to a distribution route. The demands of station $i$ are divided into $\mathrm{r}$ batches of materials, so the material in batch $r$ of station $i$ is described as $i^r$. Combined with the classical MILP model [23], the MILP model of the proposed problem is constructed.

The following list provides a description of the parameters and variables involved in the proposed model, as shown in Abbreviation.

Under JIT production mode, the first goal of the distribution is to supply the material on time and accurately to meet the requirements of the workshop. The second goal is to promote the orderly distribution and optimize vehicle loading capacity.

As minimizing the routes is on the tactical level, objective (1) should be taken as the first goal of planning and given a larger weight. When analyzing logistics optimization from the perspective of the logistics operation layer, the purpose is to ensure that each station receives the required materials within time windows and other objectives.

The objective Function (1) is to minimize the routes by using as few vehicles as possible to reduce fixed costs. The objective Function (2) is to minimize the sum of the transportation time, which is related to the feasibility of vehicle distribution. The objective Function (3) is to minimize the waiting times to improve service satisfaction.

$$Min z_1={\displaystyle \sum }_{j=1}^N{x}_{oj}^k$$

(1)

$$Min z_2={\displaystyle \sum }_{k=1}^K{\displaystyle \sum }_{i=1}^N{\displaystyle \sum }_{j=1}^N{x}_{ij}^k{t}_{ij}$$

(2)

$$Min z_3={\displaystyle \sum }_{i=1}^N{\displaystyle \sum }_{k=1}^K{w}_{ik}$$

(3)

Subject to:

$${\displaystyle \sum }_{j=1}^N{x}_{0j}^k\le 1,\forall k\in K$$

(4)

$${\displaystyle \sum }_{j=0}^V{x}_{ij}^k={\displaystyle \sum }_{r=1}^{R_i}{y}_{ir}^k,\forall k\in K,\forall i\in V$$

(5)

$${\displaystyle \sum }_{i=0}^V{x}_{ip}^k={\displaystyle \sum }_{j=0}^V{x}_{pj}^k,\forall p\in N,\forall k\in K$$

(6)

$${\displaystyle \sum }_{i=0}^V{\displaystyle \sum }_{k=1}^K{x}_{ij}^k\ge 1,\forall j\in N$$

(7)

$${\displaystyle \sum }_{k=1}^K{y}_{ir}^k=1,\forall r\in \mathrm{R},\forall k\in K$$

(8)

$${\displaystyle \sum }_{k=1}^K{\displaystyle \sum }_{r=1}^R{d}_i^r{y}_{ir}^k=d_i,\forall i\in N$$

(9)

$${\displaystyle \sum }_{i=1}^N{\displaystyle \sum }_{r=1}^R{d}_i^r{y}_{ir}^k\le Q,\forall k\in K$$

(10)

$$T_{ik}+w_{ik}+s_i+t_{ij}-T_{jk}\le \mathsf{\Omega }\left(1-x_{ij}^k\right),\forall i\in N,\forall j\in V,\forall k\in K$$

(11)

$$T_{ik}+w_{ik}+s_i+t_{ij}-T_{jk}\ge \mathsf{\Omega }\left(x_{ij}^k-1\right),\forall i\in N,\forall j\in V,\forall k\in K$$

(12)

$$T_{ik}+w_{ik}\ge e_i{\displaystyle \sum }_{j=0}^V{x}_{iJ}^k,\forall i\in N,\forall k\in K$$

(13)

$$T_{ik}\le l_i{\displaystyle \sum }_{j=0}^V{x}_{iJ}^k,\forall i\in N,\forall k\in K$$

(14)

$$w_{ik}=\{\begin{array}{c}{e}_i-T_{ik}, T_{ik}<e_{ik}\\ 0, T_{ik}\ge e_{ik}\end{array},\forall i\in N,\forall k\in K$$

(15)

$$x_{ij}^k\in \left\{0,1\right\},\forall k\in K,\forall \left(i,j\right)\in E$$

(16)

$$y_{ir}^k\in \left\{0,1\right\},\forall i\in N,\forall k\in K,\forall r\in R$$

(17)

$$T_{ik}\ge 0,\forall i\in N,\forall k\in K$$

(18)

$$w_{ik}\ge 0,\forall i\in N,\forall k\in K$$

(19)

Constraint (4) states that each vehicle used for distribution starts and ends at the warehouse. Constraint (5)–(7) represent the flow conservation law and that each station is serviced by at least one vehicle. Constraint (8) represents that each batch can only be distributed by one vehicle, meaning that each batch cannot be split. Constraint (9) indicates that the material requirements of each station need to be met. Constraint (10) enforces the maximum loading capacity constraint of the vehicle. Constraints (11) and (12) indicate that the arrival time at station$j$ is equal to the arrival time at station $i$ plus the waiting time and service time plus the travel time, only if this arc is assigned to vehicle$k$. Constraints (13) and (14) illustrate that the arrival time at station$i$ plus the waiting time must meet the time window of station $i$. The waiting time of station i is given in Constraint (15). Constraints (16)–(19) state the domain for the decision variables.

4. Hybrid Heuristic Algorithm HNSGA-II-VND

The demand discrete separable problem is a special form of SDVRP. Therefore, VRPTWDB is an NP-Hard problem [16]. Designing good meta-heuristic algorithms is the most effective solution to solve the NP-Hard problem.

The multi-objective algorithm can find a set of Pareto optimal solutions for the tradeoff between non-parallel multiple objectives, which can meet the requirements of minimum travel time and the minimum number of vehicles. Individual expressions and local searches determine the performance of the proposed algorithms. The NSGA-II algorithm has been implemented and proved to be a very effective heuristic algorithm for solving general multi-objective problems [43]. Variable Neighborhood Search (VNS) is a general local search method. Neighborhood change is applied to avoid local optimality during the descent and exploration stages.

To better solve the NP-Hard problem, the hybrid swarm intelligence algorithm with local search steps was proposed [44]. Considering that the NSGA-II and VND algorithm have good performance in solving the large-scale VRP problems, we propose the hybrid algorithm HNSGA-II-VND. The HNSGA-II-VND algorithm combines the global search capability of NSGA-II and the strong local search capability of VND, which can balance intensification and diversification.

4.1. The Flowchart of the HNSGA-II-VND

To facilitate the understanding of the proposed algorithm from the overall level, the flowchart of the HNSGA-II-VND is given in Figure 1. The following sections detail various aspects of the proposed algorithm.

The termination criteria are set: The proposed algorithm stops when the number of iterations reaches the preset maximum value, or the best solution is not improved within a certain number of iterations. If the proposed algorithm reaches the termination criteria, the Pareto optimal solution is printed out.

4.2. Solution Representation and Initialization

4.2.1. The Expression of Solutions

In the description of the solution, the material in batch $r$ of station $i$ is described as $i^r$.

The solutions can be represented by an arrangement of the distribution center 0 with the stations, in which the two stations nearest 0 and the middle part form a route. Figure 2 shows a graphical representation of the expression of solutions. The first route means that the vehicle starts from distribution center 0, arrives at station 1 (for 2 batches), 2 (for 3 batches) and 7 (for 1 batches) for unloading, and finally returns to 0.

4.2.2. Initial Solution

The hybrid algorithm is sensitive to the initial solution. A good initial solution can make the proposed algorithm converge to the optimal solutions quickly. The initial solution used by the proposed algorithm does not split the demands of all stations and generates the routes that serve all stations.

For the initial generation phase, we used the greedy insertion algorithm to generate the initial solution quickly [40]. The application of time difference ensures the fast insertion and feasibility of the station. During the initial generation phase, the stations assigned to the routes are those that result in the least increase in total route duration. Because this idea might produce a set of routes that do not serve all stations, unserved stations will be inserted into the routes. To generate a better feasible solution, the HNSGA-II-VND uses the insertion detection method for initialization.

The earliest completion time of vehicle k at station j is (vehicle passing through arc (i, j)): $E{F}_j=\max \left\{e_j+s_j,E{F}_i+t_{ij}+s_j\right\}$.

The latest start time of vehicle k at station j is: $L{S}_j=\min \left\{l_j,L{S}_k-t_{jk}-s_j\right\}$.

Theorem 1. 

Let j and u be two consecutive points in the route, and insert h between j and u under the necessary and sufficient conditions:

$$\{\begin{array}{c}E{F}_h-s_h<l_h\hfill \\ T{D}_{ju}=L{S}_u-E{F}_j\ge t_{jh}+s_h+t_{hu}\hfill \end{array}$$

The steps to plan to insert h into the ${\mathit{S}}_{\mathit{k}}$ route (${\mathit{S}}_{\mathit{k}}$ is the current viable sub-route)

Step 1: Is it overloaded after inserting h into the $S_k$ route? If it is overloaded, then k + 1, and Step 1 is repeated; otherwise, proceed to Step 2.

Step 2: Calculate the $E{F}_i$ and $L{S}_i$ of each station in $S_k$.

Step 3: Calculate the $T{D}_{ju}$ and $E{F}_h$ of all possible insertion positions of h and calculate the increment of the objective function $f_1$ after the insertion satisfies Theorem 1.

Step 4: Insert h into the position with the smallest increment in the objective function $f_1$.

Step 5: Select the next waiting insertion point to recirculate Step 1.

4.2.3. Batch Combination Operation

Discrete batches of materials are regarded as the operation objects, so batches of the same station are visited many times in the transformation, which increases the travel time cost.

Batches for the same station on the same route should be combined to avoid unnecessary route costs. The operation moves batches of the same station into one visit. After each neighborhood operation, only the first 0 of each neighboring 0 is retained.

No other points will be inserted between two batches for the same station, because this process would not make any improvement and would be a waste of computing time.

4.2.4. Batches Random Binding Operation

The number of batches far exceeds the number of stations. If the batch unit is regarded as the operation object, it will lead to a computation-intensive and time-consuming optimization process.

We design a single neighborhood moving object consisting of multiple batches at the same station [25]. The object creation operation is used for reference in the algorithm. Random binding operation of batches at the same station in the same route is considered for the generation of the actual operation objects in neighborhood operations, as shown in Figure 3. The specific operation mode is as follows:

• Select the feasible batches randomly from the current solution.
• Check whether there exist other batches that have the same station origin in the route.
• Select several consecutive batches randomly for binding as the actual operation objects for the next neighborhood operations.

4.3. Minimizing the Route Duration

After obtaining the solutions, the starting time of the vehicles departing from the warehouse can be delayed without violating the time window required by the stations, thus reducing the minimum waiting time and route duration of the solutions.

Cordeau et al. [31] apply the forward time slack to the multi-warehouse VRPTW case study. The quality of the solutions is greatly improved. Tricoire et al. [32] propose an accurate algorithm that traverses all solutions and tracks the best solution to minimize the route duration in multi-time windows team directional motion problems. Belhaiza et al. [30] propose a minimum backward time slack algorithm when studying the VRPs with multiple time windows. The proposed algorithm records the minimum waiting time and minimum delay time in the routing generation process and adjusts the arrival and departure times inversely.

We improve the method of minimizing the route duration to the proposed problem. Therefore, the forward time slack is calculated with station i as the unit (the batch unit is not differentiated in the route). The warehouse time window constraint is ignored.

Let ${\mathsf{\theta }}_i$ be the forward time slack of station i. ${\mathrm{b}}_i$ indicates the service start time at station i. ${\mathrm{s}}_i$ indicates the service time of station i. $d_c$ indicates the distribution center.

While the total waiting time in the route ${{\displaystyle \sum }}_{0<p<M}{w}_{pk}<{\mathsf{\theta }}_i,\forall k\in K$, the departure time from the distribution center is delayed to ${\mathsf{\theta }}_i$. This will lead to the departure time of the last station being delayed to ${\mathsf{\theta }}_i-{{\displaystyle \sum }}_{0<p<M}{w}_{pk},\forall k\in K$. Therefore, to obtain the minimum waiting time and the route duration of a certain route without increasing the time window constraint, the delayed departure time of vehicle K from the distribution center can be computed by Equation (21):

$${\mathsf{\theta }}_i=\underset{i\le j<d_c}{\min }\left\{{\displaystyle \sum }_{i<p\le j}{w}_{pk}+\max \left\{l_j-b_j,0\right\}\right\},\forall k\in K$$

(20)

$${\mathsf{\theta }}_0=min\left\{{\mathsf{\theta }}_0,{\displaystyle \sum }_{0<p<M}{w}_{pk}\right\},\forall k\in K$$

(21)

4.4. Selection Operator

The function of the selection operator is to update the population by selecting the current exceptional individuals. The introduction of Pareto in individual evaluation means that the non-dominant solutions are considered as the individual with relatively high fitness value. A set of solutions that rank equally defines the Pareto frontier. In addition, the elite mechanism and crowding distance are used to ensure the validity of selection [41,45]. The crowding distance defines the distance between the solution and its optimal neighbor in the pareto frontier.

Selection operator

Step 1: Input the current population $\mathrm{O}$ and previous generation $\mathrm{Q}$.

Step 2: Let $\mathrm{R}=\mathrm{O}{{\displaystyle \cup }}^{}\mathrm{Q}$ be the new population. Ranking all individuals of R according to the crowding distance and Pareto optimization.

Step 3: Select the best individuals to form the next generation population.

The elite solution is established to ensure that the best individuals are passed on to the next generation. Copying the current best solution from the previous generation to the next means that the best solution produced by all the best chromosomes will never deteriorate from one generation to the next. Although the previous optimal individual is passed on unchanged to the next generation, it is forced to compete with the new individual with a higher fitness value.

4.5. The Population Diversity

In the algorithm, the offspring must retain the non-dominant solution of the parent, which also means that the number of non-dominant solutions will not decrease. In addition, the binary championship selection operator tends to select the non-dominant solution to join the mating pool, which makes it easier for the non-dominant solution to be frequently added to the mating pool. It is possible that the two parents selected during the crossover operation are the same, which will lead to local optimization.

To maintain the diversity of the population, the population diversity strategy is added before sorting elite solutions. The principle of the strategy is to mutate the repeated solution in the whole large population after the merging of the offspring and the parent. The purpose of mutating repeated solutions into new solutions is to avoid excessive repeated solutions in the next generation, maintain the diversity of the population and expand the scope of the search.

4.6. Crossover Operator

The crossover operator is an important step in the hybrid algorithm. In this paper, the Best Cost Route Crossover (BCRC) operator is used in combination with NSGA-II multi-objective algorithm to meet the feasibility constraints while minimizing the number of vehicles and total travel time. Figure 4 illustrates a small case of the BCRC operator [42].

Note that when inserting the deleted nodes, the inserted node is first selected randomly. For example, when C1 is created, the insertion order of 7 and 4 is arbitrary. Based on satisfying the constraints, the best insertion position is the position with the smallest increment in total travel time. If no viable insertion points are found, a new route needs to be created, as shown by C1 in Figure 4.

Batches need to be checked for missing supplies after the neighborhood operation. If a batch is un-routed, the missing batch is reinserted into a new route. During the neighborhood operation, if a removed batch cannot be reinserted into any existing route, then a new route is created for this batch.

4.7. Multi-Mode Mutation Operator

The mutation operators avoid premature convergence of the HNSGA-II-VND algorithm. The mutation operators explore more neighborhood space by generating different individuals to increase the genetic diversity of the population. The multi-mode mutation pattern is adopted to increase the solution space by supplementing local search information. This section defines that the probability of selecting each operator is the same, meaning that one of the mutation operators is chosen at random [45].

Since the mutation operators may damage the chromosome structure and change the service sequence in the route, the probability of mutation in each chromosome is very low. The mutation operators randomly conduct stations exchange without considering the constraints on vehicle capacity and route time.

Mutation operator

Input: Parent population Q

Output: Offspring population O

Create O from Q by mutation routine:

While (i ≤ N)

Calculate individual mutation-probability $P_m$

If: (random() < $P_m$):

Select the mth operator

O ← Apply the mutation operator to the offspring

else

end if

i = i+1

End while

Output

The mth mutation operator is randomly selected from the defined neighborhood structures set, even though exploring the current neighborhood method increases the cost of the solution.

To define the neighborhood structures set, an appropriate balance must be struck between breaking the current solution and maintaining parts of it. This section defines neighborhoods of different sizes including Relocate Operator, Node Exchange Operator, 2-Opt, Inter-Route Swap Operator, Inter-Route Cross Operator.

4.8. Adaptive Probability

With the evolution of the population, an adaptive probability adjustment strategy is adopted to protect the optimum solution and accelerate the population evolution rate. The adaptive probability mechanism can adjust the mutation and crossover probability according to the iterative process. At the beginning of the search, large crossover probability and mutation probability are adopted to expand the search space. As the search progresses, the quality of the solutions improve, leading to an increase in fitness. The crossover probability and mutation probability should be reduced to avoid destroying the sequence of non-dominant solutions.

$P_C$ is the adaptive crossover probability. $P_1$ is the initial crossover probability. $g$ is the current iteration number. $g_{max}$ is the maximum iteration number. $P_m$ is the adaptive mutation probability. $P_2$ is the initial mutation probability.

$$P_C=\{\begin{array}{c}{P}_1, g<0.5g_{max}\hfill \\ P_1\frac{g_{max}-0.5g}{g_{max}}, 0.5g_{max}\le g\le g_{max}\hfill \end{array}$$

(22)

$$P_m=\{\begin{array}{c}{P}_2, g<0.5g_{max}\hfill \\ P_2\frac{g_{max}-0.5g}{g_{max}}, 0.5g_{max}\le g\le g_{max}\hfill \end{array}$$

(23)

4.9. The VND Algorithm

As an effective algorithm based on local search, VND has been widely used in the research of VRPs and other variants [46]. VND can use different neighborhood operators to explore systematically and compare their local optimal solutions, to effectively update the global optimal solution step by step. VND further explores the larger neighborhood space of the current solutions. The main idea of VND is that a locally optimal solution in one neighborhood is not necessarily the locally optimal solution in another neighborhood [33,44].

Based on the proposed problem, we made some improvements for local search structures that consider batch combination operation to increase the probability of disposing of the local optimum. Each non-dominant solution can be modified with local search structures. Neighborhood operators with low complexity are frequently used due to their high ranking.

The improved VND algorithm

Step 1.  Initial solution. Input the initial solution x0, x = x0.

Step 2.  Define a group of neighborhood structures  $M_m\left(m=1,2, \dots m_{max}\right)$$,\mathrm{set}m=1$ .

Step 3.  If $m\le m_{max}$ , repeat Step 4; otherwise, perform Step 5.

Step 4.  Apply the local search structures to update the solutions.

Step 4.1.  Apply random binding operation to create corresponding operation objects.

Step 4.2.  Find the local optimal solution x* of x $\mathrm{in}\mathrm{the}{M}_m$ neighborhood space.

Step 4.3.  Update solution x* according to batch combination operation and route elimination operation at the same station.

Step 4.4.  Determine the dominant relationship between x* and x. If x* dominates x, then x = x*, $m=1$.$\mathrm{Otherwise},m=m+1$.

Step 4.5.  Repeat Step 3–4.

Step 5.  Output the optimal solution x.

Neighborhood search contains the idea of population evolution to improve the optimization ability of the HNSGA-II-VND algorithm. Reasonable neighborhood structures can improve the quality of the solutions and improve the efficiency of the algorithm. The neighborhood structures adopted in this stage must conform to the constraints and be able to expand the search space.

Multiple classic neighborhood structures are adopted and divided into intra-route and inter-route operations. In each neighborhood operation, different routes $R_1$ and $R_2$ are selected as operation routes. Each operation object is randomly selected for the transformation operation, except distribution center 0. The first accept strategy is applied to reduce computing time. The strategy stops the search of the current operation and starts the next local search process after the improved solution is found for the first time.

The multi-neighborhood operations have different effects on the batch combination operation proposed above. Intra-line neighborhood operations have no effect on this operation but may reduce the total travel time. The relocate operator may eliminate the redundant routes. In the process of optimization, the routes number is constantly reduced to achieve the minimum number of vehicles.

Different intra-line operations are applied to each route, as shown in Figure 5 and Figure 6. The operations are applied in order of increasing complexity [23].

• Relocate Operator;
• Node Exchange Operator;
• When the intra-line operations are not improved, the interline operations are continued to be applied to a group of routes;
• Inter-route Relocate Operator;
• Inter-route Swap Operator;
• Inter-route Cross Operator;
• 2-opt*.

5. Computational Results and Discussion

This section summarizes the results of the experiments conducted to evaluate the performance of the proposed hybrid algorithm. All numerical experiments are performed in Python on a computer with the following specifications: Intel Core (TM) I5-7500 CPU (3.40 GHz) 8 GB RAM.

First, we describe the test environment used for testing and the parameters of the HNSGA-II-VND algorithm that are optimized by the Taguchi method. Second, we test multiple strategies included in the proposed algorithm to confirm the contribution. Third, to compare the different algorithms, we discuss the experimental results of various scale testing problems.

5.1. Testing Problems

Since the proposed problem is novel and existing research has no standard test instances for such a problem, we generated the test instances based on the classical VRPTW benchmark examples to assess the effect of the improved algorithm.

A total of 20 examples including R1 and RC1 were selected for testing. Each example was given the maximum vehicle loading capacity, and location, service time and time window for each node. The demands of each station were created according to the characteristics of the proposed problem and randomly split into 1–4 batches, indicating the batches number of different materials [23]. To study the influence of different scale instances on the hybrid algorithm, 25 nodes and 50 nodes were selected to form other part single examples.

5.2. Performance Metrics

To evaluate the behavior of different algorithms for solving the proposed multi-objective problem, two performance metrics were utilized [47].

• Generational Distance (GD). This metric is used to represent the distance between the current Pareto frontier (PF) and the real Pareto frontier (PF*).
• Inverted Generational Distance (IGD). This refers to the average distances of the nearest solutions in the Pareto solution set, considering both distribution and convergence at the same time [12,47].

These objective functions are different, so these metrics are normalized to facilitate comparison. It should be noted that the real Pareto frontier PF* is probably not available. Therefore, the obtained solutions are sorted in a non-dominated order to obtain approximate PF* in each instance.

5.3. Parameters Optimization

The parameters are very important to the effectiveness and efficiency of the algorithms. This section designs a Taguchi experiment design method to choose the best combination of experimental parameters to gain empirical insight into the influence of the parameters [45,48,49,50].

The key parameters of the hybrid algorithm: population size S, crossover probability $P_C$, mutation probability $P_m$. The orthogonal array L16(4^3) is designed with four levels for each parameter combination, as shown in

Table 2. The orthogonal array L16(4^3).

Experiments	S	${\mathit{P}}_{\mathit{C}}$	${\mathit{P}}_{\mathit{m}}$	IGD
1	100	0.6	0.1	0.1351
2	100	0.7	0.2	0.1412
3	100	0.8	0.3	0.1386
4	100	0.9	0.4	0.1165
5	120	0.6	0.2	0.1206
6	120	0.7	0.1	0.1504
7	120	0.8	0.4	0.0963
8	120	0.9	0.3	0.1173
9	140	0.6	0.3	0.1203
10	140	0.7	0.4	0.1239
11	140	0.8	0.1	0.1307
12	140	0.9	0.2	0.1185
13	160	0.6	0.4	0.1060
14	160	0.7	0.3	0.1273
15	160	0.8	0.2	0.0749
16	160	0.9	0.1	0.1461

.

Multiscale cases are selected from the test instances. Each parameter combination is tested 20 times. The average IGD values are collected as the response variables, as shown in Table 2. Figure 7 shows the main effects of the three parameters of the IGD index. According to the comparison of the parameters on algorithm performance, the variation trend of each parameter is statistically analyzed in Figure 7. Based on the results, it is recommended that the optimal parameter combination of the current algorithm settings be: $\mathrm{S}=160,P_C=0.8,P_m=0.4$. The Taguchi method has been employed on testing other part instances, and this parameter setting also works well. These values are used to conduct comprehensive experiments.

5.4. The Effectiveness of Strategies

This section evaluates the contribution of the proposed improved strategies to the hybrid algorithm. These improved strategies comprise the diversity of population and VND local search. The population diversity strategy is to avoid excessive repeated solutions in the next generation, maintain the diversity of the population and expand the scope of search. This strategy explores the global domain. The principle of the strategy is to mutate the repeated solution in the whole large population.

As an effective algorithm based on local search, VND further explores the larger neighborhood space of the current solutions. The main idea of VND is that a locally optimal solution in one neighborhood is not necessarily the locally optimal solution in another neighborhood.

There are three variants: NSGA-II-1, NSGA-II-2, NSGA-II-3. NSGA-II-1 represents the HNSGA-II-VND without both strategies. NSGA-II-2 represents the HNSGA-II-VND without the diversity of population. NSGA-II-3 indicates the HNSGA-II-VND without VND local search.

To make a fair comparison, the operation of all the compared algorithms is the same as that of the proposed hybrid algorithm. Each algorithm runs 20 experiments on the selected test instances. The average values of GD and IGD obtained by all the comparing algorithms are shown in

Table 3. The mean value of the metrics between HNSGA-II-VND and its variants in R1.

Instances	HNSGA-II-VND		NSGA-II-1		NSGA-II-2		HNSGA-II-3
	GD	IGD	GD	IGD	GD	IGD	GD	IGD
R101-25	0.0113	0.0135	0.0499	0.0814	0.0430	0.0717	0.0190	0.0162
R102-25	0.0456	0.0214	0.0958	0.1137	0.0752	0.0889	0.0653	0.0279
R103-25	0.0441	0.0336	0.2250	0.2015	0.1380	0.1837	0.0426	0.1069
R104-25	0.0329	0.0093	0.5659	0.2405	0.3589	0.2023	0.1538	0.0679
R105-25	0.0166	0.0184	0.1775	0.0990	0.1161	0.0816	0.1119	0.0326
R106-25	0.0189	0.0113	0.2438	0.1876	0.0886	0.1555	0.0392	0.0435
R107-25	0.0771	0.0512	0.2965	0.2173	0.0740	0.1320	0.2009	0.0687
R108-25	0.0127	0.0025	0.4729	0.2444	0.2884	0.2002	0.0423	0.0350
R109-25	0.1022	0.1649	0.5275	0.2707	0.3705	0.2012	0.2099	0.1884
R110-25	0.0169	0.0056	0.6580	0.2193	0.4120	0.1373	0.1755	0.0152
R111-25	0.6463	0.0553	1.2938	0.3930	1.0264	0.2963	0.7450	0.0924
R112-25	0.0021	0.0007	0.5866	0.1719	0.2450	0.0817	0.2122	0.0707
R101-50	0.0344	0.0187	0.0881	0.0890	0.0607	0.0789	0.0522	0.0257
R102-50	0.0324	0.0214	0.1127	0.0912	0.0884	0.0682	0.0753	0.0305
R103-50	0.0822	0.0297	0.1538	0.0932	0.0874	0.0815	0.1020	0.0463
R104-50	0.1374	0.0399	0.4401	0.1855	0.2893	0.1383	0.2899	0.0895
R105-50	0.0852	0.0392	0.4789	0.1674	0.3425	0.0848	0.1797	0.0468
R106-50	0.0537	0.0830	0.2724	0.1736	0.1288	0.1471	0.0961	0.1217
R107-50	0.0967	0.0855	0.4939	0.1996	0.3285	0.1678	0.1761	0.1080
R108-50	0.0038	0.0216	0.6555	0.2019	0.2879	0.1113	0.3219	0.0641
R109-50	0.0248	0.0933	0.5487	0.1994	0.3656	0.1502	0.2125	0.1125
R110-50	0.1680	0.0087	0.6565	0.2188	0.5965	0.1691	0.4785	0.0643
R111-50	0.1121	0.0069	0.6932	0.2311	0.6059	0.1627	0.3747	0.0581
R112-50	0.0327	0.0109	0.8218	0.2739	0.3117	0.0742	0.1654	0.0551
R101-100	0.1032	0.0400	0.1862	0.0839	0.1421	0.0721	0.1515	0.0548
R102-100	0.0629	0.0374	0.1537	0.0728	0.0970	0.0696	0.1383	0.0500
R103-100	0.1575	0.0484	0.2382	0.1175	0.1582	0.1045	0.1941	0.0764
R104-100	0.1670	0.0743	0.4975	0.1856	0.2883	0.1242	0.2346	0.0937
R105-100	0.2265	0.0507	0.6881	0.1818	0.5152	0.1358	0.6255	0.1402
R106-100	0.1921	0.0962	0.6595	0.2502	0.3423	0.1604	0.3839	0.1398
R107-100	0.0693	0.0934	0.7846	0.2782	0.3639	0.1602	0.3472	0.1266
R108-100	0.1435	0.0192	1.0718	0.3001	0.7234	0.2199	0.5115	0.1031
R109-100	0.1752	0.0692	0.7903	0.2802	0.3935	0.1407	0.4905	0.1406
R110-100	0.0537	0.0931	0.7328	0.2524	0.3659	0.1453	0.3237	0.1315
R111-100	0.0338	0.0753	0.7218	0.2596	0.4082	0.1720	0.4620	0.1286
R112-100	0.2117	0.0645	1.0687	0.3407	0.7864	0.2621	0.5972	0.1834
mean	0.0968	0.0447	0.5056	0.1991	0.3143	0.1398	0.2500	0.0821
hit rate	36/36	36/36	0/36	0/36	0/36	0/36	0/36	0/36

and

Table 4. The mean value of the metrics between HNSGA-II-VND and its variants in RC1.

Instances	HNSGA-II-VND		NSGA-II-1		NSGA-II-2		HNSGA-II-3
	GD	IGD	GD	IGD	GD	IGD	GD	IGD
RC101-25	0.0637	0.0547	0.2484	0.1575	0.0827	0.1416	0.0792	0.0684
RC102-25	0.0010	0.0003	0.8895	0.2965	0.3354	0.1118	0.0892	0.0085
RC103-25	0.0034	0.0144	0.5100	0.2792	0.1319	0.1425	0.0167	0.0212
RC104-25	0.0111	0.0107	0.8411	0.2841	0.3595	0.0861	0.5521	0.0876
RC105-25	0.0174	0.0127	0.3038	0.1445	0.2204	0.1048	0.0199	0.0255
RC106-25	0.0013	0.0004	0.3159	0.1053	0.1938	0.0646	0.2828	0.0943
RC107-25	0.2000	0.0667	0.3464	0.1155	0.2000	0.0667	0.2000	0.0667
RC108-25	0.2829	0.0943	0.6840	0.2280	0.5160	0.1720	0.3374	0.1125
RC101-50	0.0464	0.0541	0.5669	0.2508	0.2626	0.1532	0.3777	0.1136
RC102-50	0.1309	0.0227	1.1434	0.3802	0.9077	0.3034	0.8949	0.2938
RC103-50	0.0534	0.0227	0.4217	0.1804	0.0959	0.0942	0.1053	0.0335
RC104-50	0.0488	0.0294	0.3626	0.1736	0.2355	0.1693	0.1270	0.0485
RC105-50	0.0964	0.0483	0.3358	0.1068	0.2240	0.1042	0.1182	0.0529
RC106-50	0.0120	0.0040	0.7316	0.2439	0.5774	0.1925	0.2126	0.0709
RC107-50	0.1475	0.0492	0.7279	0.2426	0.5411	0.1804	0.4680	0.1560
RC108-50	0.0364	0.0121	0.6627	0.2209	0.4011	0.1337	0.3710	0.1237
RC101-100	0.1360	0.0745	0.4757	0.2305	0.1946	0.1236	0.2557	0.1079
RC102-100	0.0804	0.1308	0.5523	0.2671	0.3229	0.1866	0.3610	0.1589
RC103-100	0.1963	0.0561	0.9184	0.3159	0.5057	0.1740	0.6109	0.2055
RC104-100	0.0519	0.1454	0.7948	0.3664	0.4012	0.2287	0.2952	0.1681
RC105-100	0.2635	0.1297	0.6703	0.2860	0.4282	0.2042	0.5407	0.2297
RC106-100	0.1495	0.0498	0.7725	0.2368	0.6021	0.1888	0.5061	0.1217
RC107-100	0.0194	0.1383	0.7117	0.3242	0.5267	0.2645	0.3816	0.2368
RC108-100	0.0762	0.0908	0.7588	0.2966	0.5431	0.2665	0.5287	0.2602
mean	0.0886	0.0547	0.6144	0.2389	0.3671	0.1608	0.3222	0.1194
hit rate	24/24	24/24	0/24	0/24	1/24	1/24	1/24	1/24

.

Table 3 and Table 4 show that HNSGA-II-VND is superior to other algorithms in overall performance. From the perspective index values in all instances, the average value of HNSGA-II-VND index is lower than that of the HNSGA-II-VND variants. It can be seen from the results that the variants have no absolute advantage in comparison, except NSGA-II-1. Apparently, NSGA-II-3 has better results than NSGA-II-2 because the population diversity increases the possibility of exploring the solution space.

Figure 8 shows the comparison interval diagram between HNSGA-II-VND and its variants to analyze the statistically significant difference in the results. At the 95% CIs, HNSGA-II-VND has the lowest interval and the smallest distribution range. This shows that the HNSGA-II-VND algorithm has a more stable effect than its other variants. The analysis shows that the proposed strategies can contribute to the hybrid HNSGA-II-VND algorithm to solve the proposed problem.

5.5. Comparison with the Model

To evaluate the accuracy of the model, this section analyzes the results of solving the mathematical model by the CPLEX program. Based on the weighted coefficient method, the solutions in different directions can be obtained by changing the weighted coefficients of the time objectives. As the minimization of the routes is on the tactical level, objective (1) should be taken as the first goal of planning and given a larger weight. The transport time and waiting times as the second level of objectives. To calculate the results in an effective time, several groups of small-scale instances are randomly combined for testing.

The results of solving the multi-objective problem under the combination of weighted coefficients are counted in

Table 5. Results of MILP model.

No.	Weight Coefficient		R1-25			R2-25			R3-25
	W1	W2	z1	z2	z3	z1	z2	z3	z1	z2	z3
1	1	7	8	788.68	9.48	8	768.8	13.1	7	754	15
2	1	6	8	768.96	12.34
3	1	5
4	1	4	8	735.68	20.44
5	1	3
6	1	2	8	718.54	28.23				7	729	24
7	1	1	8	655.47	74.53	8	683.62	72.38	7	649.31	95.45
8	2	1	8	643.34	91.32	8	647.25	128.5	7	629.3	119.43
9	3	1	8	624.76	141.73	8	633.96	137.12	7	620.2	142.3
10	4	1	8	618.08	164.65	8	628.6	154.62
11	5	1
12	6	1				8	622.77	187.37	7	608.3	197.3
13	7	1
CPU/s			987			1056			1206
gap/%			0			0			0

. The number and distribution of Pareto solutions of different instances are different. The smaller the number of solutions, the smaller the conflict. Obviously, more than one Pareto solution is obtained in each instance, so it can be seen that there are conflicts between the objectives. For the same instances, the optimal Pareto solutions obtained by the hybrid algorithm for several times are shown in

Table 6. Results of the HNSGA-II-VND.

No.	R1-25			R2-25			R3-25
	z1	z2	z3	z1	z2	z3	z1	z2	z3
1	8	788.68	9.48	8	768.8	13.1	7	753	20.37
2	8	643.34	91.32	8	761.81	22.47	7	729	24
3	8	638.25	102.83	8	759.9	24.79	7	714.82	37.81
4	8	683.87	57.79	8	756.55	29.85	7	710.57	42.66
5	8	701.25	56.14	8	746.44	34.07	7	678.34	88.03
6	8	755.4	17.58	8	706.71	62.43	7	677.85	91.41
7	8	674.85	69.02	8	733.83	41.92	7	649.31	95.45
8	8	707.91	45.99	8	710.26	54.92	8	776.48	8.04
9	8	655.47	74.53	8	683.62	72.38	8	746.02	10.77
10	8	768.96	12.34	8	679.19	86.26	8	732.16	18.68
11	8	681.46	62.69	8	675.82	94.89	8	729.81	19.32
12	9	835.07	3.1	8	647.25	128.5	8	722.41	23.16
13	9	826.48	6.38	8	654.52	112.93	8	707.31	35.51
14	9	750.64	17.92	8	649.53	115.92	8	701.36	48.69
15							9	809.17	1.36
CPU/s	121			143			157

. A few combinations of weight coefficients cannot obtain all Pareto frontier solutions.

The Pareto solutions obtained by the algorithm are more extensive. The proposed algorithm is mainly oriented to large-scale instances and can find effective solutions in the valid computation time. The construction of mathematical models is also very important. This is because the optimal solution obtained by mathematical models can be used as a reference standard for developing meta-heuristic algorithms.

5.6. Comparison Results of Multiple Algorithms

This section further studies the effectiveness of the hybrid algorithm after the above settings in different scale cases. The hybrid algorithm and the multi-objective algorithms proposed in other literature are evaluated and analyzed.

The comparison algorithms comprise: GA [42], NSGA-II [33] and KBEA [11].

We use these comparison algorithms to solve multi-objective VRPTW. A few modifications are made to these algorithms to solve the proposed problem. All algorithms use the same coding method. To be fair, GA and NSGA-II adopt the same crossover and mutation search operations as HNSGA-II-VND. Similar to the parameter optimization of the proposed hybrid algorithm (Section 5.3), the effects of parameters on performance in these algorithms are studied. The results show that the parameters settings of NSGA-II are the same as HNSGA-II-VND. The stop criteria are set to be the same.

Table 7. The index values of all comparison algorithms in R cases.

Instances	HNSGA-II-VND		GA		NSGA-II		KBEA
	GD	IGD	GD	IGD	GD	IGD	GD	IGD
R101-25	0.0140	0.0160	0.0220	0.0271	0.0398	0.0669	0.0140	0.0160
R102-25	0.0157	0.0216	0.0853	0.0270	0.1664	0.1414	0.0334	0.0533
R103-25	0.0496	0.0369	0.0721	0.0925	0.2204	0.2166	0.0966	0.0645
R104-25	0.0255	0.0071	0.3359	0.0713	1.0159	0.3339	0.2281	0.0453
R105-25	0.0168	0.0293	0.1432	0.0652	0.2153	0.1315	0.0918	0.0517
R106-25	0.0170	0.0115	0.0222	0.0482	0.1752	0.1626	0.0189	0.0481
R107-25	0.0749	0.0506	0.2047	0.0662	0.3013	0.1938	0.1679	0.0934
R108-25	0.0122	0.0023	0.0752	0.0325	0.2864	0.2102	0.0982	0.0552
R109-25	0.0604	0.1139	0.2772	0.1432	0.3982	0.2463	0.2599	0.2211
R110-25	0.0341	0.0114	0.3187	0.1053	0.5471	0.1824	0.3384	0.1005
R111-25	0.0309	0.0515	0.2440	0.1454	0.2415	0.1480	0.0309	0.0515
R112-25	0.0000	0.0000	0.4345	0.1448	0.7527	0.2509	0.0109	0.0036
R101-50	0.0690	0.0331	0.0707	0.0354	0.1205	0.0713	0.0892	0.0335
R102-50	0.0296	0.0232	0.0666	0.0297	0.0951	0.0633	0.0577	0.0277
R103-50	0.0718	0.0349	0.1219	0.0524	0.1006	0.0937	0.0781	0.0389
R104-50	0.1457	0.0429	0.2104	0.0840	0.5570	0.2202	0.4041	0.1214
R105-50	0.0722	0.0539	0.2898	0.0868	0.3382	0.1234	0.2338	0.0814
R106-50	0.0530	0.0480	0.1492	0.0758	0.3671	0.1735	0.1252	0.0866
R107-50	0.0628	0.0968	0.4216	0.1736	0.3954	0.2299	0.2153	0.1431
R108-50	0.0100	0.1140	0.4498	0.2139	0.3921	0.2076	0.1498	0.1347
R109-50	0.1140	0.1094	0.4204	0.1421	0.5213	0.1966	0.2740	0.1424
R110-50	0.1897	0.0156	0.5984	0.1470	0.8582	0.2464	0.5841	0.1507
R111-50	0.0256	0.1416	0.3428	0.2401	0.5995	0.2863	0.4181	0.2269
R112-50	0.0923	0.0308	0.7899	0.2305	0.9758	0.3116	0.7564	0.2029
R101-100	0.0747	0.0372	0.0753	0.0453	0.1065	0.0633	0.0755	0.0400
R102-100	0.0767	0.0388	0.2151	0.0615	0.0776	0.0719	0.1341	0.0481
R103-100	0.0915	0.0358	0.2425	0.0837	0.2290	0.1067	0.1304	0.0421
R104-100	0.1360	0.1007	0.2684	0.1300	0.3717	0.1960	0.2604	0.1405
R105-100	0.2215	0.1125	0.4521	0.1499	0.4777	0.1841	0.4055	0.1333
R106-100	0.1725	0.1124	0.5564	0.2181	0.4689	0.2270	0.4433	0.1566
R107-100	0.3355	0.0971	0.4663	0.1602	0.7154	0.2430	0.6200	0.2037
R108-100	0.1118	0.0053	0.7773	0.1416	0.6415	0.1975	0.5568	0.0978
R109-100	0.2809	0.1128	0.5296	0.1990	0.9110	0.2916	0.4902	0.1884
R110-100	0.0960	0.1015	0.6138	0.2209	0.5877	0.2122	0.4530	0.1828
R111-100	0.1202	0.0924	0.3865	0.1512	0.5356	0.1870	0.4649	0.1544
R112-100	0.3673	0.1072	0.9602	0.3018	0.9711	0.3237	0.9512	0.2994
Mean	0.0937	0.0569	0.3253	0.1206	0.4382	0.1892	0.2714	0.1080
hit rate	36/36	36/36	0/36	0/36	0/36	0/36	2/36	2/36

and

Table 8. The index values of all comparison algorithms in RC cases.

Instances	HNSGA-II-VND		GA		NSGA-II		KBEA
	GD	IGD	GD	IGD	GD	IGD	GD	IGD
RC101-25	0.0251	0.0465	0.0867	0.0581	0.0665	0.1447	0.0345	0.0570
RC102-25	0.0002	0.0001	0.3042	0.0694	0.2169	0.0723	0.0974	0.0078
RC103-25	0	0.0141	0.0211	0.0161	0.3135	0.1891	0.0034	0.0224
RC104-25	0.0070	0.0084	0.3882	0.0711	0.6056	0.1761	0.0070	0.0084
RC105-25	0.0259	0.0154	0.0426	0.0201	0.1660	0.1367	0.0306	0.0299
RC106-25	0	0	0.2828	0.0943	0.0400	0.0133	0	0
RC107-25	0.0072	0.0024	0.5464	0.1821	0.0973	0.0324	0.0400	0.0133
RC108-25	0	0	0.2828	0.0943	0.6828	0.2276	0	0
RC101-50	0.1242	0.0620	0.4029	0.1384	0.6463	0.2460	0.4925	0.1702
RC102-50	0.1540	0.0292	0.5117	0.1685	0.6487	0.1935	0.4666	0.1497
RC103-50	0.0082	0.0248	0.1245	0.0776	0.5435	0.2465	0.1142	0.0494
RC104-50	0.0509	0.0310	0.1279	0.0496	0.2957	0.2542	0.0975	0.0679
RC105-50	0.0328	0.0221	0.0643	0.0332	0.2464	0.0743	0.0917	0.0431
RC106-50	0.0280	0.0093	0.2102	0.0701	0.6474	0.2158	0.3797	0.1266
RC107-50	0.2058	0.0686	0.7990	0.2663	0.4024	0.1341	0.3669	0.1223
RC108-50	0.0545	0.0182	0.4046	0.1349	0.7067	0.2356	0.1590	0.0530
RC101-100	0.1435	0.0738	0.2280	0.1115	0.4804	0.2034	0.3176	0.1118
RC102-100	0.1308	0.1027	0.6860	0.2266	0.4482	0.1558	0.3308	0.1228
RC103-100	0.1859	0.0514	0.5359	0.1537	0.7024	0.2161	0.6031	0.1773
RC104-100	0.1230	0.0639	0.6358	0.1967	0.3454	0.1437	0.4017	0.1491
RC105-100	0.1464	0.1056	0.5454	0.2500	0.5432	0.2243	0.2921	0.2279
RC106-100	0.1648	0.0549	0.5193	0.1471	0.4905	0.1635	0.4022	0.1069
RC107-100	0.0235	0.1401	0.4938	0.2533	0.5767	0.2998	0.3743	0.2278
RC108-100	0.0762	0.0907	0.7178	0.2724	0.5485	0.2806	0.3593	0.1689
mean	0.0716	0.0431	0.3734	0.1315	0.4359	0.1783	0.2652	0.1014
hit rate	24/24	24/24	0/24	0/24	0/24	0/24	3/24	3/24

summarize the results of the metrics obtained by executing all algorithms 20 times in all test instances. “Mean” represents the mean value. The last line gives the hit rate of the proposed algorithm compared with other algorithms.

In Table 7 and Table 8, we compare the proposed algorithm with other algorithms. We find that HNSGA-II-VND achieves the optimal results in the average index values of all instances, which means that HNSGA-II-VND has the best overall performance. HNSGA-II-VND obtained the minimum mean of the index values of all instances, indicating that the robustness of HNSGA-II-VND is superior to other algorithms. The metrics results confirm that the hybrid algorithm achieves a more efficient frontier of non-dominated solutions. It is normal that the same solutions can be found in a very small number of small cases.

Due to the strong randomness, Figure 9 shows the comparison interval diagram of indicators between all comparison algorithms to analyze the statistically significant difference in the results. The confidence level of all questions is set at 95%. The index mean of HNSGA-II-VND is the lowest and the range is relatively small, which indicates that HNSGA-II-VND is superior to other multi-objective algorithms. This is because HNSGA-II-VND is a hybrid algorithm, in which VND is inserted into NSGA-II to improve its local search capability.

To display the performance of these algorithms more intuitively, Figure 10 shows the Pareto frontier corresponding to the best IGD index obtained by different algorithms in the randomly selected instance. It is obvious that compared with other multi-objective algorithms, the HNSGA-II-VND algorithm can obtain better convergence and distribution of non-dominated solutions. These statistical results indicate that the HNSGA-II-VND algorithm can stably obtain better solutions for the current problems.

In summary, the experimental results show that the proposed HNSGA-II-VND algorithm is more effective and robust than other algorithms in solving the proposed multi-objective problems. The proposed algorithm performs better both in global and local search abilities, which can be used to solve this problem effectively to obtain reasonable planning routes.

6. Industrial Case Study

To test the usefulness of the proposed problem, we implement case studies with the real data originating from an automobile manufacturing enterprise. To verify the practical application effect of the hybrid algorithm, this section takes the actual data of the enterprise, for example, calculation. There are 32 station demand stations on the assembly line of the final assembly shop. The information required for material distribution tasks is obtained according to the reasonable production plan. According to BOM (Bill of material) and material properties, the vehicle delivery tasks are planned. The logistics distribution diagram of the automobile assembly workshop is shown in Figure 11.

The demands and time window of the corresponding workstations are calculated in advance. According to the actual situation of the workshop, the specified time window length is 15 min, and the unloading time is 2 min. The material distribution center arranges multiple delivery vehicles with a deadweight of 100 units to be responsible for distribution.

To make a better comparison, the results of routes optimization are compared between the empirical distribution mode, the traditional distribution mode with non-separable demand of stations, and the proposed discrete split distribution mode. The experiential distribution mode of vehicles is formulated by distribution dispatchers based on experience and material demand of each station.

In the experimental calculation, this section adopts the same parameter configuration designed by the Taguchi experiment. The experimental results obtained after 20 random experiments are shown in

Table 9. The best result of case study.

No.	Empirical Distribution			Traditional Distribution			Split Distribution
	Z1	Z2	Z3	Z1	Z2	Z3	Z1	Z2	Z3
1	13	1312	191	11	1069.88	41.82	10	1075.23	23.39
2				11	1074.35	9	10	1093.03	17.39
3				11	1092.99	0	10	1059.57	67.84
4				11	1080.36	6.06	10	1108.08	0
5				10	1151.44	0.37	10	1082.32	18.66
6				10	1147.77	3.84	10	1072.92	38.39
7				10	1056.75	88.82	10	1061.88	52.84
8				10	1083.54	47.33	10	1094.19	8
9				10	1093.1	41.27	10	1067.88	44.84
mean	13	1312	191	10.44	1094.46	26.50	10	1079.46	30.15
Improvement (100%)				+19.66	+16.58	+86.13	+23.08	+17.72	+84.21

. The comparison results between the best Pareto frontier solutions provided by the hybrid algorithm for the current production task are statistically analyzed in Figure 12.

The results show that the running solutions of the hybrid algorithm satisfy the constraints of vehicle capacity and time window. The hybrid algorithms offer a choice of nine high-quality non-dominant solutions. These solutions are good results of routes optimization, which can be used as the basis of material distribution routes selection. They have certain practical application values. The solutions given by the proposed algorithm are better than the solution implemented in the workshop based on the given objectives, which saves three vehicles and greatly saves the fixed cost.

It can be clearly seen that the results of the proposed discrete split distribution mode are better than those of the traditional distribution mode. This is because when completing the same material distribution tasks, the traditional distribution mode is difficult to ensure the utilization rate of vehicles. In most cases, the loading rates of distribution vehicles are low. Delivery vehicles in this mode also travel longer distances.

The proposed method can be applied to the actual system to solve the problem effectively. Compared with the factory implementation scheme, the schemes obtained from the hybrid HNSGA-II-VND algorithm can provide the reasonable vehicle routes and meet the production demand.

7. Conclusions

The effective operation of logistics systems has a direct impact on the production cost and delivery cycle of products. Based on the analysis of the mixed-flow assembly line, vehicle routing planning and scheduling with time windows and discrete split deliveries by batch of the demands are studied in this paper. The MILP model and the effective hybrid algorithm are proposed to solve the problem. This problem has wide application value in the manufacturing system. Combined with the material requirement characteristics of the assembly workshop, the corresponding multi-objective MILP model is constructed to describe the problem completely. The hybrid HNSGA-II-VND algorithm combines the global search capability of NSGA-II and the strong local search capability of VND, which can balance intensification and diversification well. To maintain the diversity of the population, we designed a strategy to remove repetition in the hybrid algorithm. To reduce the computing time, the adaptive expectation factors and multiple neighborhood structures were designed, especially including the same station merging operation. Finally, the effectiveness of the hybrid algorithm was verified by testing and comparison with other algorithms. The multiple strategies included in the proposed algorithm were tested to confirm the contribution. Then, the effectiveness of the proposed hybrid algorithm for solving multiple scale examples was verified. The proposed research can be effectively applied to the distribution system of the workshop to solve the problem effectively. The schemes solved by the hybrid algorithm can plan the reasonable vehicle routes to meet the production demand.

In future research, it will be necessary to further study the problems not involved. In the process of materials supply in the workshop, the solution of fuzzy uncertainty needs to be further studied, such as the temporary change in materials, and the uncertain time window. In addition, in other scenarios, there may be more than one warehouse to supply materials. Multi-warehouses constraint should be considered to divide multi-distribution areas and allocate vehicles. Furthermore, to better verify the effectiveness of the proposed algorithm, it should be compared with state-of-the-art methods.