# Corrigendum to “The Fourth-Party Logistics Routing Problem Using Ant Colony System-Improved Grey Wolf Optimization”

Authors: Fuqiang Lu; Wenjing Feng; Mengying Gao; Hualing Bi; Suxin Wang
Journal: Journal of Advanced Transportation
 (2022)

---

## Body

## 1. Introduction

Fourth-party logistics (4PL) provide customer with satisfactory supply chain solutions by integrating their own and other resources and capabilities. Fourth-party logistics routing problem (4PLRP) is a basic and important research issue in the process of supply chain integration. For various operation modes of 4PL, scholars have studied the routing problems for more than a decade by using a variety of optimization methods.Aiming to maximize the degree of customers’ satisfaction, benefit third-party logistics providers, and minimize transport costs simultaneously, Huang et al. [1] proposed a mathematical model of the point-to-point single task path optimization in fourth-party logistic with soft time window (4PLRPSTW). A harmony search (HS) algorithm was designed to solve this problem. In order to find a route of the minimum cost with constraints under uncertain environments, Huang et al. [2] established a mathematical model with fuzzy duration time (4PLRPF), and a two-step fuzzy simulation genetic algorithm was designed to solve the problem. On the basic of the reliability theory and multigraph, Li et al. [3] proposed a chance-constrained programming model of 4PL, aiming at minimum cost and time constraints. A messy genetic algorithm with double arrays encoding was designed to solve the problem. To handle the 4PL routing optimization under emergency conditions that lack historical data, Huang et al. [4] proposed an uncertain programming model (UPM) with uncertain delivery time. The UPM was based on uncertainty theory (UT), and the solution can satisfy the belief degree constrain effectively. Lu et al. [5] proposed an uncertain delivery time control model for 4PLRP, which considers the selection of 3PL suppliers and transportation routes, delivery time, and transportation cost. An improved genetic algorithm (I-GA) was developed to solve this problem. Aiming at dynamic time planning in the logistics transportation network, Cui et al. [6] proposed a 4PL model considering the transportation, stay, and transit cost with the time constraint, and the dynamic time is updated at the transit nodes. A two-phase solution method based on the ant colony optimization (TACO) was established to solve the problem. In many previous studies, it was considered that there was only one 3PL between nodes to undertake transportation tasks. However, it is more in line with actual operation to consider that more 3PLs between nodes jointly undertake transportation tasks.Traditional expected utility theory (EUT) and prospect theory (PT) [7] in behavioral economics are two important theories currently studying people’s behavioral decision-making patterns. However, for EUT, it is difficult to make a reasonable explanation for the phenomena such as Allais paradox [8]. PT uses the amount of change in wealth as a reference for people to make decisions, but converts the objective probabilities into subjective probabilities with greater randomness. Wang and Kong [9] extended the two and proposed a new consumer behavior decision-making mechanism: proportional utility theory. Ratio utility theory (RUT) holds that the change value and the absolute value both affect the behavior of the decision maker [10]. It has been proved based on the panel data of 1997–2010 [11] and has been successfully applied to many fields [12, 13]. As a result, the value function of RUT has more advantages in representing customers’ risk aversion in 4PLRP.In 1991, M. Dorigo et al. proposed the first ACO algorithm-ant system (AS) and successfully used it to solve the TSP problem [14–16]. Experimental results show that the AS algorithm has strong robustness and the ability to find better solutions, but it also has some defects, such as slow convergence speed and stagnation. In order to improve these problems, M. Dorigo [17] proposed the ant colony system (ACS). This algorithm is a simple and mature ant colony algorithm framework, which has the advantages of parallel computing of group biomimetic algorithms. As a probabilistic algorithm, it is usually used to find the best path on the graph. For a complex combination optimization problem such as 4PLPR that simultaneously performs path optimization and 3PL vendor selection, if the ant’s autonomous optimization is used to obtain the initial solution, it will fully reflect the advantages of the ant colony algorithm and solve the problem efficiently.In 2014, Mirjalili et al. [18] proposed a novel intelligent algorithm, grey wolf optimization (GWO) algorithm, based on the cooperative hunting method of the grey wolf group intelligence. This algorithm fully simulates the pyramid-like social hierarchy of grey wolves and the mechanism of communication and sharing among grey wolves. In the GWO algorithm, the initial position, convergence factor, and the position update formula of the wolves will have a certain impact on the performance of the algorithm. In order to improve the optimization accuracy and search efficiency of the GWO algorithm, researchers have done a lot of improvement work. To strengthen the algorithm’s local search capability and speed up the convergence rate, Saremi et al. [19] improved the GWO algorithm by introducing a dynamic evolution population operator, and the efficiency of the algorithm was proved by solving function optimization problems. Heidari and Pahlavani [20] introduced Levy flight and greedy selection strategies into the GWO algorithm and improved the hunting stage. Experimental results and statistical tests show that LGWO’s performance is significantly better than other algorithms. In order to overcome the premature convergence problem of the GWO algorithm, Wang et al. [21] combined the basic GWO with Gaussian distribution estimation, called the GEDGWO algorithm, and proposed a poor solution repair based on the Gaussian distribution method to modify the morbid distribution of the population. This algorithm applied the Gaussian disturbance random walk method to enhance the local exploration ability and was used to solve the problem of multi-UAV multitarget tracking path planning. In [22], in order to improve the search ability of the grey wolf, a GWO algorithm RW-GWO based on random walk was proposed. The algorithm was subjected to the nonparametric test, Wilcoxon test, and performance index analysis, and the improved algorithm showed effectiveness. In [23], a multiobjective GWO algorithm based on decomposition was proposed. By defining a neighborhood relationship between the scaled subproblems of the decomposition multiobjective problem, the Pareto solution was cooperatively approximated. This algorithm was comparable to the performance of the six most advanced biological heuristic techniques. The comparison shows that it has high performance on the famous benchmark problem and two real-world engineering problems. Scholar Kamboj [24] designed a mixed particle swarm optimization algorithm of grey wolf, which combines grey wolf optimization algorithm and particle swarm optimization algorithm to solve the unit commitment problem (UCP). GWO has the advantages of fewer parameters and fast convergence speed and is widely used in engineering fields. Through the improvement of the algorithm convergence mechanism and weight update mechanism, GWO will highlight its advantages in the solution of choosing multi-3PL in 4PLRP.The uncertain environment always leads to delays in transport time and increase in transport cost, and the related research of 4PL routing optimization based on risk preference is still relatively rare. In this study, a 4PLRP model based on RUT is established to indicate the decision problem when customer is risk aversion. This study applied an improved grey wolf optimizer (IGWO) [25] to solve the proposed model. The IGWO attempts to improve the convergence factor and the proportional weight to solve the model. IGWO adopts the convergence factor based on the exponential law change, balances the global search and local search ability of the algorithm, and introduces the proportional weight based on the step of Euclidean distance to update the grey wolf’s position and accelerate the convergence speed of the algorithm. Among them, by considering the dependence of GWO on the initial population, the principle of the ant colony system (ACS) is employed to provide a good initial population for IGWO. Therefore, the designed algorithm is named as ant colony system-improved grey wolf optimizer (ACS-IGWO). Finally, the validity of the model is verified by the experiments, which further emphasizes the necessity of considering the customer’s behavior characteristics in 4PLRP.The contributions of our work are manyfold, which are listed as follows.(1)
In the transportation scheme provided to customers, the situation of multiple 3PL cooperative transportation between the two nodes will be considered to avoid tardiness or overspending caused by insufficient single 3PL carrying capacity(2)
Based on the ratio utility theory, this study establishes the utility function and 4PLRP model from the perspective of customers’ risk appetite for time and cost(3)
First, the grey wolf algorithm is improved by using a nonlinear convergence factor, so that the algorithm expands the search range at the early stage and finds more high-quality solutions. Second, the introduction of proportional-based dynamic weighting strategy improves the iteration method of the next generation of individuals, making the algorithm has powerful adaptability and faster convergence speed.The rest of this study is organized as follows. Section2 is problem description and parameter definition. Section 3 describes our 4PLRP mathematical model based on ratio utility theory. In Section 4, the mechanism of population initialization based on ACS is introduced and the improve strategy of GWO is proposed. Section 5 presents the model parameters analyses, validates the advantages of cooperation of multiple 3PLs, and evaluates the performance of ACS-IGWO heuristic.


## 2. Problem Description

The 4PLRP can be defined as selecting routes and third-party logistics (3PL) for a 4PL to optimize the supply chain. A key problem that the 4PL should consider is how to optimize the transport route from the source to the destination to ensure the maximization of utility and the satisfaction of relevant constraints. Assume that a 4PL undertakes a task design for transporting supplies. Information of the current transportation network, 3PLs, due date, and total cost are obtained. And one or more 3PLs can be selected between nodes, and the transport route must meet the 3PL provider’s transportation capacity and reputation constraints at the same time. To describe the problem more concisely and clearly, an undirected multigraph is used to describe the transportation network, as shown in Figure1.Figure 1
Multigraph of 4PLRP.In Figure1, the transportation network of 4PLRP is represented by a multigraph G (V, E), where V is the set of nodes and E is the set of edges. Node v1 represents the source node vs and node v7 represents the destination node ve, where other nodes are the transit nodes. Each node has cost and time properties, that is, the transit cost and time should be considered when the task transits at this node. Each edge represents a 3PL provider who wants to undertake a portion of the task between pairs of nodes. Each edge also has cost and time properties, which are the cost and time needed for the 3PL provider to complete the transportation service between two nodes. In such a graph, there may be more than one edge (3PL provider) between pairs of nodes to undertake the task. For a given task from the source to the destination, the 4PL needs to select a feasible route, including nodes and edges. To formulate the problem, parameters are defined in Table 1.Table 1
Parameter definition of the model.
ParameterParameter descriptionnThe number of city nodes in the transportation network also indicates the label of end pointrijThe number of edges (3PLs) between nodevi and node vjeijkThe kth edge between node vi and node vj, where i,j∈1,2,…,n, and k is the index number of edgesTiThe time of nodevi, including the time of processing, inventory, loading, and unloadingCiThe cost of nodevi, including the cost of processing, inventory, loading, and unloadingTijkThe time of the 3PL provider in edgeeijkCijkThe cost of the 3PL provider in edgeeijkTRThe transport time of the routeCRThe transport cost of the routeTijkSThe minimum time of the 3PL provider in edgeeijkTijkLThe maximum time of the 3PL provider in edgeeijkTijSThe shortest transport time between nodevi and node vjTijLThe longest transport time between nodevi and node vjTSThe minimum transport time of the routeTLThe maximum transport time of the routeCijkSThe minimum cost of the 3PL provider in edgeeijkCijkLThe maximum cost of the 3PL provider in edgeeijkCijSThe maximum transport cost between nodevi and node vjCijLThe maximum transport cost between nodevi and node vjCSThe minimum transport cost of the routeCLThe maximum transport cost of the routeT0The due date required by customerC0The total cost required by customerφ1Weight of the timeφ2Weight of the costQijkTransport capacity of the 3PL provider in edgeeijkDijkReputation of the 3PL provider in edgeeijkQThe required transport capacity from the taskΔQiThe amount of task that transits at the nodeptThe delay probability of the transportation taskspcThe probability of overspending of the transportation tasksDThe required reputation required by customerRThe route from nodevs to node vehijkThe proportion of transport tasks undertaken by edgeeijkxijWhether the 3PL provider is selected between nodevi and node vj to undertake the transportation task, which is a 0-1 variablexijkWhether the edgeeijk is selected to undertake the transportation task, which is a 0-1 variableyiWhether to select nodevi as a transit node, which is a 0-1 variable


## 3. Mathematical Model

### 3.1. Description of the Distribution Mechanism

3PL providers are responsible for the transportation task between the two nodes. A reasonable task distribution for the total task volume is needed, and the certain amount of tasks assigned to the 3PL under the condition that the capacity of each selected 3PL is satisfied, which is shown in Figure2. Assume that three 3PLs are selected between nodes v2 and v5, and the task distribution ratios of the three providers are h251=0.4, h252=0.4, and h253=0.2. In order to obtain a better performance of the task assignment and the optimal transport plan, this study distributes the total task randomly according to the number of 3PLs and the transport capacity conditions. In addition, as can be seen from Figure 2, 80% of the goods at node v5 are transited, so the transit time and the transit cost at node v5 are, respectively, 80% T5 and 80% C5.Figure 2
Schematic diagram of the distribution mechanism.


### 3.2. Mathematical Model Based on Ratio Utility Theory

People often focus on the change value and the absolute value when making decisions in an uncertain environment. Therefore, this study uses the ratio of the two to define the utility of the consumers, and the three hypotheses [15] of RUT are as follows:(1)
Consumer will maintain his wealth stateG0 until there exists an option ΔG. G0=G1,G2,G3,…,Gi and ΔG=ΔG1,ΔG2,ΔG3,…,ΔGi, where G denotes the quantity of wealth and “i” indicates the category of wealth.(2)
When option exists, the ratio of change valueΔG to the final value G0+ΔG indicates consumer utility, named RUT as follows: RU=∑iΔGi/Gi+ΔGi; if RU > 0, consumer has a willing to receive the option ΔG; otherwise, consumer will reject the option.(3)
The sensitivity of RUT is limitedIf multiple 3PLs undertake transportation tasks, the shortest transport time and the longest transport time should be calculated by equations (1) and (2), respectively. The minimum and maximum transport costs between the two nodes are calculated by equations (3) and (4), respectively.(1)TijS=Maxxij1Tij1S,xij2Tij2S,…,xijkTijkS,…,xijrijTijrijS.(2)TijL=Maxxij1Tij1L,xij2Tij2L,…,xijkTijkL,…,xijrijTijrijL.(3)CijS=∑i=1n∑j=1n∑ijkrijhijkQQijkCijkSxijk.(4)CijL=∑i=1n∑j=1n∑k=1rijhijkQQijkCijkLxijk.TR, CR, TS, TL, CS, and CL are as follows:(5)TR=∑i=1n∑j=1ntijxij+∑i=1nΔQiQTiyi,CR=∑i=1n∑j=1n∑k=1rijhijkQQijkcijkxijk+∑i=1nΔQiQCiyi,TS=∑i=1n∑j=1nTijSxij+∑i=1nΔQiQTiyi,TL=∑i=1n∑j=1nTijLxij+∑i=1nΔQiQTiyi,CS=∑i=1n∑j=1n∑k=1rijhijkQQijkCijkSxijk+∑i=1nΔQiQCiyi,CL=∑i=1n∑i=1n∑k=1rijhijkQQijkCijkLxijk+∑i=1nΔQiQCiyi.Tij∼UTijS,TijL and Cijk∼UCijkS,CijkL, so TR∼UTS,TL and CR∼UCS,CL. Then, the tardiness and overspending probability of the transport task are as follows:(6)pt=0,T0>TL,TL−T0TL−TS,TS<T0<TL1,T0<TS,,pc=0,C0>CL,CL−C0CL−CS,CS<C0<CL,1,C0<CS.This study proposes a decision-making method based on RUT for customer who is risk aversion on tardiness and overspending. The method gives the corresponding RUT value, and in order to make it more reasonable, it refers to the form of the prospect theoretical value function and introduces the risk attitude coefficientυ and the loss aversion coefficient λ. This method calculates and obtains the utility value of the transport time and transport cost of the distribution plan. The customer’s required due date T0 and cost C0 are set as reference points, respectively. Therefore, the utility functions of T(R) and C(R) are shown as follows:(7)Ut=−λptTR−T0TRυ,T0−TR≤0,1−ptT0−TRTRυ,T0−TR>0,(8)Uc=−λpcCR−C0CRυ,C0−CR≤0,1−pcC0−CRCRυ,C0−CR>0,Based on RUT, the model considering a risk aversion customer is as follows:(9)maxω1Ut+ω2Uc,(10)s.t.∑i=1n∑j=1n∑k=1nhijk=1,vi,vj∈R,(11)hijkQ≤Qijk,vi,vj∈R,(12)0≤hijk≤1,vi,vj∈R,(13)xijk=1,0≤hijk≤1,0,hijk=0,(14)xijk=1,vi,vj∈R,0,otherwise,(15)yi=1,vi∈R,0,otherwise,(16)Dijk≥Dxijk,(17)y1=yn=1,(18)R=vs,…,vi,1h1/ij,…,kh1/ij,…,rijh1/ij,vj,…,ve.Equation (9) is the objective function; it means to maximum the total utility, including the utility of transport time and the utility of transport cost; equation (10) represents that for the task distribution ratio of the 3PL provider, the sum of them is 1; equation (11) represents the capability constraint of 3PL; equation (12) represents that the task distribution ratio of each 3PL which undertakes the transportation task between the two nodes is in [0, 1]; equation (13) represents that whether 3PL undertakes the transportation task between the two nodes; equation (14) represents that whether choose 3PL between two nodes; equation (15) represents that whether the transport route passes through the node; equation (16) represents that the reputation of 3PL must meet customer’s needs; equation (17) represents that both the start node and the destination node must be selected; equation (18) is used to ensure the selected route is legal, which means that nodes and edges of the route must be interconnected one by one, and must start from the source and end at the destination.


## 3.1. Description of the Distribution Mechanism

3PL providers are responsible for the transportation task between the two nodes. A reasonable task distribution for the total task volume is needed, and the certain amount of tasks assigned to the 3PL under the condition that the capacity of each selected 3PL is satisfied, which is shown in Figure2. Assume that three 3PLs are selected between nodes v2 and v5, and the task distribution ratios of the three providers are h251=0.4, h252=0.4, and h253=0.2. In order to obtain a better performance of the task assignment and the optimal transport plan, this study distributes the total task randomly according to the number of 3PLs and the transport capacity conditions. In addition, as can be seen from Figure 2, 80% of the goods at node v5 are transited, so the transit time and the transit cost at node v5 are, respectively, 80% T5 and 80% C5.Figure 2
Schematic diagram of the distribution mechanism.


## 3.2. Mathematical Model Based on Ratio Utility Theory

People often focus on the change value and the absolute value when making decisions in an uncertain environment. Therefore, this study uses the ratio of the two to define the utility of the consumers, and the three hypotheses [15] of RUT are as follows:(1)
Consumer will maintain his wealth stateG0 until there exists an option ΔG. G0=G1,G2,G3,…,Gi and ΔG=ΔG1,ΔG2,ΔG3,…,ΔGi, where G denotes the quantity of wealth and “i” indicates the category of wealth.(2)
When option exists, the ratio of change valueΔG to the final value G0+ΔG indicates consumer utility, named RUT as follows: RU=∑iΔGi/Gi+ΔGi; if RU > 0, consumer has a willing to receive the option ΔG; otherwise, consumer will reject the option.(3)
The sensitivity of RUT is limitedIf multiple 3PLs undertake transportation tasks, the shortest transport time and the longest transport time should be calculated by equations (1) and (2), respectively. The minimum and maximum transport costs between the two nodes are calculated by equations (3) and (4), respectively.(1)TijS=Maxxij1Tij1S,xij2Tij2S,…,xijkTijkS,…,xijrijTijrijS.(2)TijL=Maxxij1Tij1L,xij2Tij2L,…,xijkTijkL,…,xijrijTijrijL.(3)CijS=∑i=1n∑j=1n∑ijkrijhijkQQijkCijkSxijk.(4)CijL=∑i=1n∑j=1n∑k=1rijhijkQQijkCijkLxijk.TR, CR, TS, TL, CS, and CL are as follows:(5)TR=∑i=1n∑j=1ntijxij+∑i=1nΔQiQTiyi,CR=∑i=1n∑j=1n∑k=1rijhijkQQijkcijkxijk+∑i=1nΔQiQCiyi,TS=∑i=1n∑j=1nTijSxij+∑i=1nΔQiQTiyi,TL=∑i=1n∑j=1nTijLxij+∑i=1nΔQiQTiyi,CS=∑i=1n∑j=1n∑k=1rijhijkQQijkCijkSxijk+∑i=1nΔQiQCiyi,CL=∑i=1n∑i=1n∑k=1rijhijkQQijkCijkLxijk+∑i=1nΔQiQCiyi.Tij∼UTijS,TijL and Cijk∼UCijkS,CijkL, so TR∼UTS,TL and CR∼UCS,CL. Then, the tardiness and overspending probability of the transport task are as follows:(6)pt=0,T0>TL,TL−T0TL−TS,TS<T0<TL1,T0<TS,,pc=0,C0>CL,CL−C0CL−CS,CS<C0<CL,1,C0<CS.This study proposes a decision-making method based on RUT for customer who is risk aversion on tardiness and overspending. The method gives the corresponding RUT value, and in order to make it more reasonable, it refers to the form of the prospect theoretical value function and introduces the risk attitude coefficientυ and the loss aversion coefficient λ. This method calculates and obtains the utility value of the transport time and transport cost of the distribution plan. The customer’s required due date T0 and cost C0 are set as reference points, respectively. Therefore, the utility functions of T(R) and C(R) are shown as follows:(7)Ut=−λptTR−T0TRυ,T0−TR≤0,1−ptT0−TRTRυ,T0−TR>0,(8)Uc=−λpcCR−C0CRυ,C0−CR≤0,1−pcC0−CRCRυ,C0−CR>0,Based on RUT, the model considering a risk aversion customer is as follows:(9)maxω1Ut+ω2Uc,(10)s.t.∑i=1n∑j=1n∑k=1nhijk=1,vi,vj∈R,(11)hijkQ≤Qijk,vi,vj∈R,(12)0≤hijk≤1,vi,vj∈R,(13)xijk=1,0≤hijk≤1,0,hijk=0,(14)xijk=1,vi,vj∈R,0,otherwise,(15)yi=1,vi∈R,0,otherwise,(16)Dijk≥Dxijk,(17)y1=yn=1,(18)R=vs,…,vi,1h1/ij,…,kh1/ij,…,rijh1/ij,vj,…,ve.Equation (9) is the objective function; it means to maximum the total utility, including the utility of transport time and the utility of transport cost; equation (10) represents that for the task distribution ratio of the 3PL provider, the sum of them is 1; equation (11) represents the capability constraint of 3PL; equation (12) represents that the task distribution ratio of each 3PL which undertakes the transportation task between the two nodes is in [0, 1]; equation (13) represents that whether 3PL undertakes the transportation task between the two nodes; equation (14) represents that whether choose 3PL between two nodes; equation (15) represents that whether the transport route passes through the node; equation (16) represents that the reputation of 3PL must meet customer’s needs; equation (17) represents that both the start node and the destination node must be selected; equation (18) is used to ensure the selected route is legal, which means that nodes and edges of the route must be interconnected one by one, and must start from the source and end at the destination.


## 4. Algorithm Design

Since the selection of multiple 3PL 4PLRP is a complex combinatorial optimization problem and is also an NP-hard problem, therefore, the grey wolf optimizer (GWO) is applied to solve the model. GWO has a simple structure, requires few parameters to be set, is easy to experiment with coding, and is widely used in terms of attribute simplicity and feature selection. However, GWO has disadvantages of low solution accuracy and slow convergence speed. Therefore, this study attempts to revise the convergence factor and the proportional weight of GWO. The ACS algorithm is used to generate the initial solution, and the GWO algorithm is used to solve the problem, which not only reflects the advantages of the algorithm but also solves the problem. Therefore, an ACS-IGWO is designed.


### 4.1. Population Initialization

The GWO has a certain dependence on the initial population. In order to find the optimal route, the principle of ACS was used to provide the initial population for GWO.


#### 4.1.1. Coding Method

This study adopts the integer number encoding method. First, the multigraph is represented by an adjacency list shown in Table2. If there is no edge between two nodes, the value of corresponding element in the list will be set to 0; otherwise, the element value will be set to the number of edges between the two nodes.Table 2
Logistics network adjacency list (7 nodes).
Node123456710340000200222003000304040000323500000036000000270000000


#### 4.1.2. Transfer Strategy

The ant determines the transfer direction according to the pseudorandom scale rule of the following equation:(19)j=argmaxj∈JsiτijkNgαηijkβ,q≤q0,Eq.20,otherwise.where Ng is the current number of iterations, τijkNg is the pheromone concentration on eijk in iteration Ng, ηijk is its route heuristic information, Jsi is the set of nodes that ant s can visit after visiting node vi, q is a random number (q∼U0,1), and q0 is a fixed algorithm parameter, q0∈0,1; when q>q0, the ant will determine the node to be transferred according to the following equation:(20)Pijk=argτijkNgαηijkβ∑j∈JsiτijkNgαηijkNgβ,j∈Jsi,0,otherwise,where α and β are the algorithm parameters, representing the importance of pheromone and route heuristic information, respectively. ηijk=1/Tijk+1/Cijk+1/Ti+1/Ci is the heuristic information, and it will provide in advance. Before and after the node, if a same one 3PL provider is responsible for the transportation task, the time and the cost of the node will not be calculated, that is, ηijk=1/Tijk+1/Cijk.


#### 4.1.3. Construction of Feasible Solution

The feasible solution of ACS is gradually generated by ants in multigraph. At the initial moment of the iteration, the ant is placed at the start nodevs, and the ant performs a proportional selection according to the transition probability, thereby determining the next direction of the transfer. After the ant moves to the next node, the transfer is performed according to the transfer strategy, and such step is repeated until reaching the destination node ve. And R, which is the optimal route from the start to the destination, is obtained. In order to prevent the situation where the ant cannot reach the destination node, it is assumed that if the next node transferred by the ant is a node that the ant has traveled before, the connected edge with the node is not considered within the feasible transfer direction.


#### 4.1.4. Global Pheromone Update

In ACS, the pheromone of the optimal ant route was updated in each cycle globally.(21)τijkbsNg=1−ρτijkbsNg+ρΔτijkbsNg,where Δτijkbs=Q/TRbs+CRbs is the pheromone increment of the current optimal route. TRE and CRbs are the transport time and transport cost of the current optimal route. ρ and Q are the algorithm parameters. ρ is the pheromone volatilization coefficient, ρ∈0,1.Q is the certain amount of pheromone concentration.


#### 4.1.5. Local Pheromone Update

After all the ants complete a transfer, perform a local update of the pheromone according to the following equation, whereτ0 is the initial pheromone concentration.(22)τijkNg=1−ρτijkNg+ρτ0.


### 4.2. Grey Wolf Optimization

Grey wolves are top carnivores, located at the top of the food chain. Their lifestyles are mostly group-based with an average of 5 to 12 wolves. The characteristics of the GWO can be described in the following three aspects:


#### 4.2.1. Social Hierarchy

The grey wolf group has a strict hierarchy,α is the highest level of grey wolf, β is a subordinate of α, δ obeys the command of α and β, and the grey wolf with the lowest rank is called ω.


#### 4.2.2. Search and Surround Prey

Suppose that the number of grey wolves isN, the position of the ith grey wolf is Xi, the optimal solution of the group is α, the suboptimal is β, the third is δ, and the other individuals are ω. Then the grey wolf searches and surrounds. The mathematical model of prey behavior is described as follows:(23)D=C⋅Xpt−Xt,Xt+1=Xpt−Α⋅D.Among them,t represents the current number of iterations, Xpt is the position vector of the prey, Xt is the position vector of the grey wolf individual, D is the distance between the individual and the prey, and A and C are the coefficient vectors.(24)a=2−2ttmax,(25)A=2⋅a⋅r2−a,(26)C=2⋅r1,r1,r2 are the random numbers, r1,r2∈ [0, 1]. tmax is the maximum number of iterations, and a is the convergence factor, decreasing linearly from 2 to 0 as the number of iterations.


#### 4.2.3. Hunting

During the hunting process of the grey wolf group, whenα perceives the position of the prey, α will combine with other wolves at the leadership level to command the entire group, guide the wolves to close the prey from all directions, and further achieve the ultimate predation. The location of grey wolves ω in the population is determined by the location of α, β, and δ. Xt represents the position of ω. Xα, Xβ, and Xδ are the positions of α, β, and δ, respectively. C1, C2, and C3 represent random perturbations for α, β, and δ, respectively. X1, X2, and X3 represent the positions guided by α, β, and δ, respectively. The updated location X(t + 1) represents the new position of ω:(27)Dα=C1⋅Xαt−Xt,X1=Xαt−A1Dα,Dβ=C2⋅Xβt−Xt,X2=Xβt−A2Dβ,Dδ=C3⋅Xδt−Xt,X3=Xδt−A3Dδ,Xt+1=X1+X2+X33.


### 4.3. Nonlinear Convergence Factor Based on Exponential Function

IfA>1, the algorithm will look for more possible solutions globally, and the scope of the search is further expanded, which search globally, and the convergence speed is fast; when A>1, the algorithm will shrink the search range and seek the current range more locally, and the convergence speed is slow. Therefore, the value of A is related to the global search and local search capabilities of GWO. It can be known from equation (25) that A changes with the convergence factor a, and the convergence factor a decreases linearly from 2 to 0 as the number of iterations. But the algorithm is not linear in the process of constant convergence. The decreasing convergence factor does not fully reflect the actual optimization search process. Therefore, this study adopts a convergence factor based on the change of exponential law [33], and its modified expression is equation (28). The comparison before and after the improvement is shown in Figure 3.(28)a=2−21e−1×et/tmax−1.Figure 3
The change of the convergence factor.It can be seen from Figure3 that the original convergence factor a is linearly decreasing and decreases at the same rate during the iterative process, and the improved convergence factor a is a curve based on the exponential law, which is reduced in the initial stage of the iteration. In the initial stage of the iteration, a decreases slowly, so that a keeps a large value for a long time, and A keeps a larger value for a longer period to improve the searching efficiency. At the end of the iteration, a decreases faster, so that a keeps a smaller value for a longer period, and A keeps a smaller value for a longer period to improve the search accuracy. Thus, the global search and local search capabilities of the algorithm are balanced.


### 4.4. Proportional Weight Based on Step of Euclidean Distance

(29)W1=X1X1+X2+X3,(30)W2=X2X1+X2+X3,(31)W3=X3X1+X2+X3,(32)Xt+1=X1×W1+X2×W2+X3×W33,where W1, W2, and W3 represent the learning rates for α, β, and δ, respectively. The proportional weight is used to update the grey wolf position to speed up the convergence speed of the algorithm.The transportation plan update method represented byω is shown in Figure 4. The transportation plan represented by αω in Figure 5 is the result that ω learns from α wolf. ω learns β and δ in the same way, and the most effective scheme in ω, βω, and δω is called ωmax. Compare ωmax with X(t + 1); if ωmax≥Xt+1, update the transportation plan represented by the wolf to the transportation plan; if ωmax≤Xt+1, rerandomly distribute 3PLs between nodes to undertake the proportion of transportation tasks 10 times, and the maximum utility value is ω10. If ω10≥Xt+1, the transportation plan is updated; otherwise, the transportation plan waits for the next iteration to update.Figure 4
The update mode of transportation plan.Figure 5
Optimal route scheme for example I.


### 4.5. Procedures of ACS-IGWO

Step 1: set the maximum number of iterationstmax, A, and C. A and C are randomly generated.Step 2: the population is initialized withNP individuals by the principle of ACSStep 3: calculate the utility value of all grey wolves in the population and sort according to the utility value. Then, select the top three best wolves and record their positions asXα, Xβ, and Xδ.Step 4: update the positions of other grey wolves in the population by equations (29)–(32)Step 5: updateα by equation (28), and update A and C by equations (25) and (26), respectively.Step 6: iftmax is reached, the calculation is stopped, output the current optimal position Xα and its utility value; otherwise, go to step 3.


## 4.1. Population Initialization

The GWO has a certain dependence on the initial population. In order to find the optimal route, the principle of ACS was used to provide the initial population for GWO.


### 4.1.1. Coding Method

This study adopts the integer number encoding method. First, the multigraph is represented by an adjacency list shown in Table2. If there is no edge between two nodes, the value of corresponding element in the list will be set to 0; otherwise, the element value will be set to the number of edges between the two nodes.Table 2
Logistics network adjacency list (7 nodes).
Node123456710340000200222003000304040000323500000036000000270000000


### 4.1.2. Transfer Strategy

The ant determines the transfer direction according to the pseudorandom scale rule of the following equation:(19)j=argmaxj∈JsiτijkNgαηijkβ,q≤q0,Eq.20,otherwise.where Ng is the current number of iterations, τijkNg is the pheromone concentration on eijk in iteration Ng, ηijk is its route heuristic information, Jsi is the set of nodes that ant s can visit after visiting node vi, q is a random number (q∼U0,1), and q0 is a fixed algorithm parameter, q0∈0,1; when q>q0, the ant will determine the node to be transferred according to the following equation:(20)Pijk=argτijkNgαηijkβ∑j∈JsiτijkNgαηijkNgβ,j∈Jsi,0,otherwise,where α and β are the algorithm parameters, representing the importance of pheromone and route heuristic information, respectively. ηijk=1/Tijk+1/Cijk+1/Ti+1/Ci is the heuristic information, and it will provide in advance. Before and after the node, if a same one 3PL provider is responsible for the transportation task, the time and the cost of the node will not be calculated, that is, ηijk=1/Tijk+1/Cijk.


### 4.1.3. Construction of Feasible Solution

The feasible solution of ACS is gradually generated by ants in multigraph. At the initial moment of the iteration, the ant is placed at the start nodevs, and the ant performs a proportional selection according to the transition probability, thereby determining the next direction of the transfer. After the ant moves to the next node, the transfer is performed according to the transfer strategy, and such step is repeated until reaching the destination node ve. And R, which is the optimal route from the start to the destination, is obtained. In order to prevent the situation where the ant cannot reach the destination node, it is assumed that if the next node transferred by the ant is a node that the ant has traveled before, the connected edge with the node is not considered within the feasible transfer direction.


### 4.1.4. Global Pheromone Update

In ACS, the pheromone of the optimal ant route was updated in each cycle globally.(21)τijkbsNg=1−ρτijkbsNg+ρΔτijkbsNg,where Δτijkbs=Q/TRbs+CRbs is the pheromone increment of the current optimal route. TRE and CRbs are the transport time and transport cost of the current optimal route. ρ and Q are the algorithm parameters. ρ is the pheromone volatilization coefficient, ρ∈0,1.Q is the certain amount of pheromone concentration.


### 4.1.5. Local Pheromone Update

After all the ants complete a transfer, perform a local update of the pheromone according to the following equation, whereτ0 is the initial pheromone concentration.(22)τijkNg=1−ρτijkNg+ρτ0.


## 4.1.1. Coding Method

This study adopts the integer number encoding method. First, the multigraph is represented by an adjacency list shown in Table2. If there is no edge between two nodes, the value of corresponding element in the list will be set to 0; otherwise, the element value will be set to the number of edges between the two nodes.Table 2
Logistics network adjacency list (7 nodes).
Node123456710340000200222003000304040000323500000036000000270000000


## 4.1.2. Transfer Strategy

The ant determines the transfer direction according to the pseudorandom scale rule of the following equation:(19)j=argmaxj∈JsiτijkNgαηijkβ,q≤q0,Eq.20,otherwise.where Ng is the current number of iterations, τijkNg is the pheromone concentration on eijk in iteration Ng, ηijk is its route heuristic information, Jsi is the set of nodes that ant s can visit after visiting node vi, q is a random number (q∼U0,1), and q0 is a fixed algorithm parameter, q0∈0,1; when q>q0, the ant will determine the node to be transferred according to the following equation:(20)Pijk=argτijkNgαηijkβ∑j∈JsiτijkNgαηijkNgβ,j∈Jsi,0,otherwise,where α and β are the algorithm parameters, representing the importance of pheromone and route heuristic information, respectively. ηijk=1/Tijk+1/Cijk+1/Ti+1/Ci is the heuristic information, and it will provide in advance. Before and after the node, if a same one 3PL provider is responsible for the transportation task, the time and the cost of the node will not be calculated, that is, ηijk=1/Tijk+1/Cijk.


## 4.1.3. Construction of Feasible Solution

The feasible solution of ACS is gradually generated by ants in multigraph. At the initial moment of the iteration, the ant is placed at the start nodevs, and the ant performs a proportional selection according to the transition probability, thereby determining the next direction of the transfer. After the ant moves to the next node, the transfer is performed according to the transfer strategy, and such step is repeated until reaching the destination node ve. And R, which is the optimal route from the start to the destination, is obtained. In order to prevent the situation where the ant cannot reach the destination node, it is assumed that if the next node transferred by the ant is a node that the ant has traveled before, the connected edge with the node is not considered within the feasible transfer direction.


## 4.1.4. Global Pheromone Update

In ACS, the pheromone of the optimal ant route was updated in each cycle globally.(21)τijkbsNg=1−ρτijkbsNg+ρΔτijkbsNg,where Δτijkbs=Q/TRbs+CRbs is the pheromone increment of the current optimal route. TRE and CRbs are the transport time and transport cost of the current optimal route. ρ and Q are the algorithm parameters. ρ is the pheromone volatilization coefficient, ρ∈0,1.Q is the certain amount of pheromone concentration.


## 4.1.5. Local Pheromone Update

After all the ants complete a transfer, perform a local update of the pheromone according to the following equation, whereτ0 is the initial pheromone concentration.(22)τijkNg=1−ρτijkNg+ρτ0.


## 4.2. Grey Wolf Optimization

Grey wolves are top carnivores, located at the top of the food chain. Their lifestyles are mostly group-based with an average of 5 to 12 wolves. The characteristics of the GWO can be described in the following three aspects:


### 4.2.1. Social Hierarchy

The grey wolf group has a strict hierarchy,α is the highest level of grey wolf, β is a subordinate of α, δ obeys the command of α and β, and the grey wolf with the lowest rank is called ω.


### 4.2.2. Search and Surround Prey

Suppose that the number of grey wolves isN, the position of the ith grey wolf is Xi, the optimal solution of the group is α, the suboptimal is β, the third is δ, and the other individuals are ω. Then the grey wolf searches and surrounds. The mathematical model of prey behavior is described as follows:(23)D=C⋅Xpt−Xt,Xt+1=Xpt−Α⋅D.Among them,t represents the current number of iterations, Xpt is the position vector of the prey, Xt is the position vector of the grey wolf individual, D is the distance between the individual and the prey, and A and C are the coefficient vectors.(24)a=2−2ttmax,(25)A=2⋅a⋅r2−a,(26)C=2⋅r1,r1,r2 are the random numbers, r1,r2∈ [0, 1]. tmax is the maximum number of iterations, and a is the convergence factor, decreasing linearly from 2 to 0 as the number of iterations.


### 4.2.3. Hunting

During the hunting process of the grey wolf group, whenα perceives the position of the prey, α will combine with other wolves at the leadership level to command the entire group, guide the wolves to close the prey from all directions, and further achieve the ultimate predation. The location of grey wolves ω in the population is determined by the location of α, β, and δ. Xt represents the position of ω. Xα, Xβ, and Xδ are the positions of α, β, and δ, respectively. C1, C2, and C3 represent random perturbations for α, β, and δ, respectively. X1, X2, and X3 represent the positions guided by α, β, and δ, respectively. The updated location X(t + 1) represents the new position of ω:(27)Dα=C1⋅Xαt−Xt,X1=Xαt−A1Dα,Dβ=C2⋅Xβt−Xt,X2=Xβt−A2Dβ,Dδ=C3⋅Xδt−Xt,X3=Xδt−A3Dδ,Xt+1=X1+X2+X33.


## 4.2.1. Social Hierarchy

The grey wolf group has a strict hierarchy,α is the highest level of grey wolf, β is a subordinate of α, δ obeys the command of α and β, and the grey wolf with the lowest rank is called ω.


## 4.2.2. Search and Surround Prey

Suppose that the number of grey wolves isN, the position of the ith grey wolf is Xi, the optimal solution of the group is α, the suboptimal is β, the third is δ, and the other individuals are ω. Then the grey wolf searches and surrounds. The mathematical model of prey behavior is described as follows:(23)D=C⋅Xpt−Xt,Xt+1=Xpt−Α⋅D.Among them,t represents the current number of iterations, Xpt is the position vector of the prey, Xt is the position vector of the grey wolf individual, D is the distance between the individual and the prey, and A and C are the coefficient vectors.(24)a=2−2ttmax,(25)A=2⋅a⋅r2−a,(26)C=2⋅r1,r1,r2 are the random numbers, r1,r2∈ [0, 1]. tmax is the maximum number of iterations, and a is the convergence factor, decreasing linearly from 2 to 0 as the number of iterations.


## 4.2.3. Hunting

During the hunting process of the grey wolf group, whenα perceives the position of the prey, α will combine with other wolves at the leadership level to command the entire group, guide the wolves to close the prey from all directions, and further achieve the ultimate predation. The location of grey wolves ω in the population is determined by the location of α, β, and δ. Xt represents the position of ω. Xα, Xβ, and Xδ are the positions of α, β, and δ, respectively. C1, C2, and C3 represent random perturbations for α, β, and δ, respectively. X1, X2, and X3 represent the positions guided by α, β, and δ, respectively. The updated location X(t + 1) represents the new position of ω:(27)Dα=C1⋅Xαt−Xt,X1=Xαt−A1Dα,Dβ=C2⋅Xβt−Xt,X2=Xβt−A2Dβ,Dδ=C3⋅Xδt−Xt,X3=Xδt−A3Dδ,Xt+1=X1+X2+X33.


## 4.3. Nonlinear Convergence Factor Based on Exponential Function

IfA>1, the algorithm will look for more possible solutions globally, and the scope of the search is further expanded, which search globally, and the convergence speed is fast; when A>1, the algorithm will shrink the search range and seek the current range more locally, and the convergence speed is slow. Therefore, the value of A is related to the global search and local search capabilities of GWO. It can be known from equation (25) that A changes with the convergence factor a, and the convergence factor a decreases linearly from 2 to 0 as the number of iterations. But the algorithm is not linear in the process of constant convergence. The decreasing convergence factor does not fully reflect the actual optimization search process. Therefore, this study adopts a convergence factor based on the change of exponential law [33], and its modified expression is equation (28). The comparison before and after the improvement is shown in Figure 3.(28)a=2−21e−1×et/tmax−1.Figure 3
The change of the convergence factor.It can be seen from Figure3 that the original convergence factor a is linearly decreasing and decreases at the same rate during the iterative process, and the improved convergence factor a is a curve based on the exponential law, which is reduced in the initial stage of the iteration. In the initial stage of the iteration, a decreases slowly, so that a keeps a large value for a long time, and A keeps a larger value for a longer period to improve the searching efficiency. At the end of the iteration, a decreases faster, so that a keeps a smaller value for a longer period, and A keeps a smaller value for a longer period to improve the search accuracy. Thus, the global search and local search capabilities of the algorithm are balanced.


## 4.4. Proportional Weight Based on Step of Euclidean Distance

(29)W1=X1X1+X2+X3,(30)W2=X2X1+X2+X3,(31)W3=X3X1+X2+X3,(32)Xt+1=X1×W1+X2×W2+X3×W33,where W1, W2, and W3 represent the learning rates for α, β, and δ, respectively. The proportional weight is used to update the grey wolf position to speed up the convergence speed of the algorithm.The transportation plan update method represented byω is shown in Figure 4. The transportation plan represented by αω in Figure 5 is the result that ω learns from α wolf. ω learns β and δ in the same way, and the most effective scheme in ω, βω, and δω is called ωmax. Compare ωmax with X(t + 1); if ωmax≥Xt+1, update the transportation plan represented by the wolf to the transportation plan; if ωmax≤Xt+1, rerandomly distribute 3PLs between nodes to undertake the proportion of transportation tasks 10 times, and the maximum utility value is ω10. If ω10≥Xt+1, the transportation plan is updated; otherwise, the transportation plan waits for the next iteration to update.Figure 4
The update mode of transportation plan.Figure 5
Optimal route scheme for example I.


## 4.5. Procedures of ACS-IGWO

Step 1: set the maximum number of iterationstmax, A, and C. A and C are randomly generated.Step 2: the population is initialized withNP individuals by the principle of ACSStep 3: calculate the utility value of all grey wolves in the population and sort according to the utility value. Then, select the top three best wolves and record their positions asXα, Xβ, and Xδ.Step 4: update the positions of other grey wolves in the population by equations (29)–(32)Step 5: updateα by equation (28), and update A and C by equations (25) and (26), respectively.Step 6: iftmax is reached, the calculation is stopped, output the current optimal position Xα and its utility value; otherwise, go to step 3.


## 5. Numerical Experiments

A 4PL company receives a logistics order with a total carrying capacity of 100. In an uncertain environment, it is shown that customers are risk aversion to tardiness and overspending. In order to complete transportation tasks successfully, the 4PL integrator provides transportation solutions to meet customer’s requirements for transport time and transport cost. An ACS-IGWO (ant colony system-improved grey wolf optimizer) is applied to solve 4PLRP. To compare the efficiency of the proposed algorithm, ACS, D-ACS (ant colony algorithm for heuristic dynamic pheromone update strategy) [26], and ACS-GWO (ant colony system-grey wolf optimizer) are used to solve the resulting optimization problem.This section presents three experiments with different problem scales of 7-node example, 15-node example, and 30-node example, and they are named example I, example II, and example III, respectively.First, the material data are described. Second, the problem is analyzed from three points of view. The influence of model parameters on the results is analyzed, the customer’s risk attitude is studied, and the transportation modes are comparatively researched. Finally, the algorithm is studied from the point of parameter combination, and the validity and efficiency of the algorithm is verified by example tests between algorithms.


### 5.1. Data Description

For example I, the data of nodes and edges are shown in Tables3 and 4, respectively; there are a total of 33 edges. In addition, when the case scale is fifteen nodes for example II, there are a total of 91 edges, and when the case scale is thirty nodes for example III, there are a total of 197 edges. Assume that the transportation task requires a 3PL supplier with a load capacity of 100 and a credibility requirement of 8.Table 3
Basic data of nodes for example I.
NodeTiCi135234344423534652745Table 4
Basic data of edges for example I.
Node (i)Node (j)Edge (k)TijkSTijkLCijkSCijkLQijkDijk12112.713.411.412.280812219.821.314.715.51201012319.820.4131570913121.623.213.714.3901013216.517.216.918.61201113317.818.518.419.56091349.210.411.413.2801023110.811.312.813.7100142329.710.311.512.5701024112.513.410.511.3901024210.511.29.911.61101225111.712.216.817.51401125210.711.214.915.67083419.710.38.69.5110934211.812.511.812.680934312.713.310.411.7901036112.713.418.519.790936212.913.810.111.780103639.910.58.79.670113641718.86.47.2601545111.812.37.88.970124529.810.113.515.5110745312.813.211.812.6901146113.715.314.51612094628.99.47.99.5110747112.513.211.412.2501447210.411.310.211.3701047311.612.313.815.590857112.513.513.415.1801257212.613.612.713.1701557311.712.114.815.79013671910.59.910.1110106721718.56.47909


### 5.2. Problem Analysis and Discussion

#### 5.2.1. Analysis of Model Parameters

This section will first take the example I as an example to analyze the influence of the model parametersT0,C0,, φ1, and φ2 on the selection of transportation plans. Then, obtain the optimal combination parameters, the optimal route, and its utility. The utility is represented by U on different examples.(1) The Impact of T0 on the Optimal Route. First, set C0=65, ω1=0.4, and ω2=0.6. 1 ⟶ (10.8, 30.2) ⟶ 2 indicates that the route passes through v1 and v2, selects 3PL 1 to complete 80% transportation tasks, and selects 3PL 3 to complete the rest 20% transportation task between the two nodes. And R represents the optimal route. For different value of T0, the result is shown in Table 5.Table 5
The impact ofT0 on the transportation plans.
T0RUCsCLTSTL351 ⟶ (20.5, 30.1, 40.4) ⟶ 3 ⟶ (10.5, 20.1, 30.4) ⟶ 4 ⟶ (20.4, 30.6) ⟶ 7−0.3152.355.155.258.8401 ⟶ (20.3, 30.1, 40.6) ⟶ 3 ⟶ (20.8, 30.2) ⟶ 4 ⟶ (10.4, 20.2, 30.4) ⟶ 7−0.1952.054.354.658.2451 ⟶ (30.5, 40.5) ⟶ 3 ⟶ (11) ⟶ 4 ⟶ (20.4, 30.6) ⟶ 7−0.0753.455.354.858.5501 ⟶ (20.6, 40.4) ⟶ 3 ⟶ (11) ⟶ 4 ⟶ (20.5, 30.5) ⟶ 70.0554.956.953.357.6551 ⟶ (21) ⟶ 3 ⟶ (10.5, 20.5) ⟶ 6 ⟶ (10.3, 20.7) ⟶ 70.1856.858.753.354.0601 ⟶ (20.5, 30.1, 40.4) ⟶ 3 ⟶ (10.5, 20.1, 30.4) ⟶ 4 ⟶ (20.4, 30.6) ⟶ 70.2252.355.155.258.8651 ⟶ (21) ⟶ 3 ⟶ (10.7, 20.1, 30.2) ⟶ 4 ⟶ (20.2, 30.8) ⟶ 70.2855.057.555.658.7701 ⟶ (21) ⟶ 3 ⟶ (10.8, 30.2) ⟶ 4 ⟶ (20.4, 30.6) ⟶ 70.3256.858.553.758.7751 ⟶ (21) ⟶ 3 ⟶ (10.7, 20.1, 30.2) ⟶ 4 ⟶ (20.4, 30.6) ⟶ 70.3655.457.955.658.7801 ⟶ (21) ⟶ 3 ⟶ (20.1, 30.2, 40.5) ⟶ 6 ⟶ (10.4, 20.6) ⟶ 70.4055.857.854.156.8It can be seen from Table5, a series of optimal routes are obtained with variance of T0, and the insight are given as follows. (1) Different values T0 generate different optimal routes. When T0  = 50, the optimal route that satisfy the customer’s requirements for T(R) and C(R) is found for the first time. (2) It can be seen from the grey part of Table 5 that C(R) is extended from [53.4, 55.3] to [54.9, 56.9], and T(R) is reduced from [54.8, 58.5] to [53.3, 57.6]. The optimal route that satisfies the customer’s requirements is obtained. When only considering the tardiness of the route, the customer is willing to choose a route with relatively higher cost, where both T(R) and C(R) meets the requirements. It also shows that the customer is risk aversion and has behavioral characteristics of against loss. (3) If the customer’s requirements for T(R) become looser, the utility of objective function becomes the main factor, and it affects the optimal route.(2) The Impact of C0 on the Optimal Route. Set T0= 65, ω1=0.4, and ω2=0.6. Change the value of C0, and the results are shown in Table 6 and the findings are as follows. (1) Different values C0 generate different optimal routes. When C0  = 40, the optimal route that satisfies the customer’s requirements for T(R) and C(R) is found. (2) It can be seen from the grey part of the table that T(R) is extended from [53.3, 57.6] to [55.0, 58.6], and C(R) is reduced from [54.9,56.9] to [50.8,54.8], and the optimal route meeting the customer’s requirements is found. If the cost of route is only in overspending, the customer is willing to choose a route with relatively longer time, and both T(R) and C(R) satisfy the requirements. This shows that customer is risk aversion and has behavioral characteristics of loss avoidance. (3) With the increase of C0, the customer’s requirements for the optimal route transport cost are looser, and the objective function utility value becomes the main factor determining the optimal route.Table 6
The impact ofC0 on the transportation plans.
C0RUCsCLTSTL401 ⟶ (21) ⟶ 3 ⟶ (10.7, 20.1, 30.2) ⟶ 4 ⟶ (20.2, 30.8) ⟶ 7−0.2355.057.555.658.7451 ⟶ (21) ⟶ 3 ⟶ (20.3, 30.2, 40.5) ⟶ 6 ⟶ (10.4, 20.6) ⟶ 7−0.1255.857.854.156.8501 ⟶ (20.6, 40.4) ⟶ 3 ⟶ (11) ⟶ 4 ⟶ (20.5, 30.5) ⟶ 7−0.0154.956.953.357.6551 ⟶ (20.1, 30.2, 40.7) ⟶ 3 ⟶ (20.7, 30.3) ⟶ 4 ⟶ (10.2, 20.1, 30.7) ⟶ 70.1450.854.855.058.6601 ⟶ (30.5, 40.5) ⟶ 3 ⟶ (11) ⟶ 4 ⟶ (20.4, 30.6) ⟶ 70.2153.455.354.858.5651 ⟶ (21) ⟶ 3 ⟶ (10.7, 20.1, 30.2) ⟶ 4 ⟶ (20.2, 30.8) ⟶ 70.2855.057.555.658.7701 ⟶ (20.6, 40.4) ⟶ 3 ⟶ (11) ⟶ 4 ⟶ (20.5, 30.5) ⟶ 70.3354.956.953.357.6751 ⟶ (21) ⟶ 3 ⟶ (20.3, 30.2, 40.5) ⟶ 6 ⟶ (10.4, 20.6) ⟶ 70.3755.857.854.156.8801 ⟶ (21) ⟶ 3 ⟶ (10.7, 20.1, 30.2) ⟶ 4 ⟶ (20.2, 30.8) ⟶ 70.4255.057.555.658.7851 ⟶ (21) ⟶ 3 ⟶ (10.8, 30.2) ⟶ 4 ⟶ (20.3, 30.7) ⟶ 70.4755.858.353.758.7(3) The Impact of φ1, φ2 on the Optimal Route. Set T0=65 and C0=65, and change the value of φ1 and φ2, the results are shown in Table 7. The best parameters combination for the example I is T0= 65, C0= 65, ω1=0.4, and ω2=0.6. The utility value is 0.28. The best route is 1 ⟶ (20.4, 30.5, 40.1) ⟶ 3 ⟶ (10.5, 20.1, 30.4) ⟶ 4 ⟶ (20.4, 30.6) ⟶ 7.Table 7
The impact ofφ1 and φ2 on the transportation plans.
φ1φ2RUCSCLTSTL0.10.91 ⟶ (30.6, 40.4) ⟶ 3 ⟶ (10.1, 20.3, 30.6) ⟶ 4 ⟶ (20.3, 30.7) ⟶ 70.2848.651.852.857.40.20.81 ⟶ (20.2, 30.2, 40.6) ⟶ 3 ⟶ (20.4, 30.6) ⟶ 4 ⟶ (20.5, 30.5) ⟶ 70.2749.652.853.657.20.30.71 ⟶ (20.5, 30.1, 40.4) ⟶ 3 ⟶ (10.2, 20.5, 30.3) ⟶ 4 ⟶ (20.3, 30.7) ⟶ 70.2750.654.553.457.00.40.61 ⟶ (20.4, 30.5, 40.1) ⟶ 3 ⟶ (10.5, 20.1, 30.4) ⟶ 4 ⟶ (20.4, 30.6) ⟶ 70.2851.954.754.057.60.50.51 ⟶ (21) ⟶ 3 ⟶ (10.8, 30.2) ⟶ 4 ⟶ (20.3, 30.7) ⟶ 70.2855.858.353.758.70.60.41 ⟶ (21) ⟶ 3 ⟶ (20.3, 30.2, 40.5) ⟶ 6 ⟶ (10.4, 20.6) ⟶ 70.2755.857.854.156.80.70.31 ⟶ (21) ⟶ 3 ⟶ (20.3, 30.2, 40.5) ⟶ 6 ⟶ (10.4, 20.6) ⟶ 70.2855.857.854.156.80.80.21 ⟶ (21) ⟶ 3 ⟶ (20.3, 30.2, 40.5) ⟶ 6 ⟶ (10.4, 20.6) ⟶ 70.2855.857.854.156.80.90.11 ⟶ (21) ⟶ 3 ⟶ (10.5, 20.5) ⟶ 6 ⟶ (10.3, 20.7) ⟶ 70.2756.858.753.354.0Through the above analysis, the seven-node instance parameter combination, the optimal utility value, and its optimal path based on ACS-IGWO are obtained. Through the same method, we can get the fifteen and thirty-node case parameter combinations, as shown in Table8. The optimal path of example I is shown in Figure 5:Table 8
Model parameter combination and the best utility of different node cases.
NodeT0C0φ1φ2Utility765650.40.60.2815901200.40.60.30301501800.40.60.30


#### 5.2.2. Customer’s Risk Attitude Analysis

This section analyses the influence of customer’s behavior characteristics on the optimal route selection. In example I, ifλ=1 and υ=1 in equations (7) and (8), the mathematical model is based on RUT that is in line with the situation when the customer is risk-neutral attitude.The ACS-IGWO runs 100 times, and the average utility is obtained, which is organized in Tables9 and 10, respectively. Figure 6 is drawn based on the information of the two tables. RU represents the average utility of the customer based on RUT, and EUM represents the average utility of the customer based on EUT.Table 9
Comparison of RU and EUM (T0 = 65).
C020253035404550556065707580EUM−0.30−0.22−0.17−0.12−0.07−0.020.030.080.130.180.230.280.33RU−0.67−0.56−0.44−0.31−0.19−0.070.050.140.210.280.330.380.43Table 10
Comparison of RU and EUM (T0 = 65).
C0859095100105110115120125130135140145150EUM0.380.430.480.530.580.630.680.730.780.830.880.930.981.03RU0.480.520.560.600.640.670.700.730.750.770.790.800.810.82Figure 6
Comparison of RU and EUM (T0  = 65).The customer has higher requirements for transport cost whenC0 is small, so the two models take a negative value, and the satisfaction of the routes is lower, which means the customer is in a loss state. As C0 increases gradually, the values of the two models are increasing gradually, and customer is more satisfied with the routes, which means the customer transforms from the loss state to the gain state. That is, let the “0” utility be the dividing line; if the customer’s utility is negative, the customer is in a loss state; and if it is positive, the customer is in a gain state. From Figure 6, the insights are as follows:(1)
If the utility is negative, the customer is in a loss state. Compared to EUM, the utility of RU is smaller, which points out that customer is risk aversion. The obtained optimal route will not delay, but it overspends seriously. Customer is risk averse to overspending, so they are less satisfied with the optimal routes.(2)
If the utility is positive andC0≤55, the customer is at the beginning of the gain state. Compared to EUM, RU is smaller, which points out that that customer is risk aversion. In Table 6, when C0  = 55, C(R) has the maximum value 54.8, which is almost equal to C0; in this perspective, the optimal routes will not overspend. When C0<55, in Table 6, it can be found that the optimal routes still overspend and do not meet customer’s expectations.(3)
If55≤C0≤120, the customer is in the continuous gain state. Compared to EUM, RU is larger, which points out that that customer is risk aversion. The optimal routes obtained neither overspend nor delay. Compared with EUM, the RU value is larger, which indicates that the customer is more conservative and risk aversion, and expect for higher utility value.(4)
IfC0≥120, the customer is in the “final” gain state. Compared to EUM, RU is smaller, which points out that that the customer is risk aversion. Although there are no delays and overspend, the rate of increase in customer satisfaction with the transportation plan is gradually decreasing. Therefore, customers are not willing to excessively reduce their cost requirements, resulting in unnecessary cost waste.(5)
Overall, the RU curve of Figure6 shows that the slope of the loss state is higher than the gain state, which indicates that the customer is more sensitive to the loss than the gain.WhenC0 is fixed and C0=65, the trend of the result is the same as in Figure 6, so the analysis will not be repeated here. With the customer’s requirements for transport time being gradually reduced, the customer’s risk attitude is the same as the situation that the customer’s transport cost requirements are gradually reduced. In summary, the model can accurately describe the psychological characteristics of people. It also shows that the consideration and analysis of customer behavior characteristics in 4PLRP have a certain practical value and practical significance.


#### 5.2.3. Comparative Analysis of Transportation Modes

In this section, two transportation modes are considered: single 3PL and multiple 3PLs undertake transport tasks. The impact of different transportation modes in the cases of simple tardiness, simple overspend, no tardiness, and overspend is studied. Among them, simple tardiness, simple overspend, and definitely no tardiness and overspend are represented byF, O, and N, and single 3PL and multiple 3PL are represented by A and M. When T0 = 65, C0 = 65, φ1=0.4, and φ2=0.6, the results are listed in Table 11.Table 11
Transportation plans with different transportation modes.
SituationModeUCSCLTSTLFA−0.2256.858.154.758.0M−0.1952.054.354.658.2OA−0.2557.759.156.058.5M−0.2355.057.555.658.7NA0.2558.962.155.459.0M0.2855.057.555.658.7In Table10, the following conclusions are obtained. In three cases, transportation cost of a single 3PL mode is higher than multiple 3PLs mode, and the time of both is almost equal. A higher transportation cost corresponds to a lower the utility value, that is, the customer’s satisfaction with the optimal route is lower.


### 5.3. Algorithm Analysis and Discussion

In order to test the effectiveness of the proposed algorithm, first, the optimal parameters combination of ACS-IGWO and the optimal results are obtained. Then, the effectiveness of ACS-IGWO is verified by comparing with ACS, D-ACS (ant colony algorithm for heuristic dynamic pheromone update strategy) [26], and ACS-GWO.


#### 5.3.1. Parameters Analysis of ACS-IGWO

In this study, the Taguchi method [27] is used to obtain the optimal parameters combination of the four algorithms.In this section, the Taguchi method is used to obtain the optimal parameters combination of ACS-IGWO. Parameters such asα, β,Q, and ρ are selected to be tested. The application software Minitab 17 is introduced to perform parameter testing on ACS-IGWO for example I. The orthogonal test table is shown in Table 12. The signal-to-noise ratio main effect diagram and the mean main effect diagram are shown in Figures 7 and 8, respectively.Table 12
The orthogonal test table.
Running numberαβQρU122500.60.2800222500.60.2800322500.60.2800423800.70.2800523800.70.2801623800.70.28007251000.80.27998251000.80.27999251000.80.27991032800.80.28001132800.80.28001232800.80.280113331000.60.280114331000.60.279915331000.60.28001635500.70.28011735500.70.27991835500.70.279919521000.70.280020521000.70.280021521000.70.28002253500.80.27992353500.80.28002453500.80.27992555800.60.28012655800.60.28002755800.60.2800Figure 7
Signal-to-noise ratio main effect diagram.Figure 8
Mean main effect diagram.From Figures7 and 8, the average value of the maximum signal-to-noise ratio of the four parameters is −11.0578, −11.0567, −11.0562, and −11.0566, respectively. To make the means of mean of the parameters the largest, α, β, Q, and ρ should take level one. Therefore, α=3, β=5, Q=50, ρ=0.6 is the parameters combination of ACS-IGWO. Other parameters are NG, NP, and q0, which have a fixed value in each algorithm. The parameters combination of the four algorithms is shown in Table 13.Table 13
Parameters combination of the four algorithms.
ExampleAlgorithmαβQρNGNPq0IACS22500.730300.7D-ACS22500.63030—ACS-GWO22500.730300.7ACS-IGWO22500.730300.7IIACS231000.7100500.7D-ACS231000.610050—ACS-GWO231000.7100500.7ACS-IGWO231000.7100500.7IIIACS231000.7150800.7D-ACS231000.615080—ACS-GWO231000.7150800.7ACS-IGWO231000.7150800.7


#### 5.3.2. Comparative Analysis of the Four Algorithms

The four algorithms are run 100 times. To test the performance of the algorithms, the relevant performance parameters of the algorithm are defined. The average maximum utility is best, the average minimum utility is worst, the average utility is mean, the variance of utility isS, and the average running time is time (s). For the three examples, the study is solved by ACS, D-ACS, ACS-GWO, and ACS-IGWO. The results are shown in Table 14.Table 14
Comparison of algorithms of different node cases.
ExampleAlgorithmBestWorstMeanSIACS-IGWO0.28−0.61−0.010.01ACS-GWO0.28−0.62−0.010.03ACS0.28−0.62−0.030.04D-ACS0.28−0.61−0.040.05IIACS-IGWO0.30−0.51−0.020.12ACS-GWO0.30−0.51−0.040.23ACS0.28−0.67−0.020.85D-ACS0.30−0.53−0.050.14IIIACS-IGWO0.30−0.71−0.010.12ACS-GWO0.30−0.73−0.010.32ACS0.25−0.82−0.061.76D-ACS0.30−0.73−0.030.30In Table13, for example I, all algorithms can find the global optimal solution. While for example II and example III, ACS cannot find the global optimal solution and perform worse with the increase of problem scale. We can find that ACS-IGWO performs best among the four algorithms. From the compare of variance value and mean average utility of algorithms in each example, it can be found that ACS-IGWO has more stable convergence interval.In addition, for the three examples, the optimal results are collected and represented by box plot, which is shown in Figures9–11. For example I, in Figure 9, although the four boxes are different in size, the difference of the concentration range is almost zero, which indicates that the four algorithms can find the optimal route accurately. For Figure 10, example II, the optimal results of ACS-IGWO, ACS-GWO, and D-ACS are better and more reliable than ACS, which indicates that ACS-IGWO, ACS-GWO, and D-ACS have strong ability to skip out of local optimum. The solution effects of ACS-GWO and D-ACS are similar. For Figures 10 and 11, example II and example III, the median line of ACS-IGWO is closer to 0.3, and the upper and lower limits are smaller, which shows that ACS-IGWO can find more reliable results. Its feasibility and effectiveness for solving this kind of routing optimization problem are verfied.Figure 9
Box plot of example I.Figure 10
Box plot of example II.Figure 11
Box plot of example III.


#### 5.3.3. Comparative Analysis of Efficiency of the Four Algorithms

The problem size and running time of each algorithm in cases are shown in Table15. In this problem, the transportation path of the goods from the start node to the end point will be planned, and the distribution ratio of 3PL between every two nodes in the path will be obtained. Between the two nodes, the number of 3PL is between 2 and 4, and the proportion of goods allocated by 3PL is between 0 and 1. This ratio needs to be obtained through optimization. It can be found that the scale of the problem is very large, and there are even an infinite variety of results. Therefore, we conservatively estimated the scale of the problem and stipulated that two 3PLs were taken between every two nodes, and the quantity of goods was distributed 1 : 1. For example 1, there are at least 124 transportation possibilities; for example 2, there are at least 761 transportation possibilities; for example 3, there are at least 32,775 transportation possibilities. In Table 15, the rate is the ratio of time to problem size.Table 15
Comparison of algorithms of different node cases.
ExampleProblem sizeAlgorithmTime (s)RateI124ACS-IGWO0.020.000161ACS-GWO0.040.000323ACS0.10.000807D-ACS0.050.000403II761ACS-IGWO0.30.000394ACS-GWO0.50.000657ACS1.50.001971D-ACS0.40.000526III32775ACS-IGWO1.20.000037ACS-GWO5.50.000168ACS22.50.000686D-ACS1.60.000049The running time of ACS-IGWO, ACS-GWO, and D-ACS is shorter than ACS, and ACS-IGWO is the fastest one. For running time in example III, the convergence speed can be arranged as ACS-IGWO > D-ACS > ACS-GWO > ACS.


## 5.1. Data Description

For example I, the data of nodes and edges are shown in Tables3 and 4, respectively; there are a total of 33 edges. In addition, when the case scale is fifteen nodes for example II, there are a total of 91 edges, and when the case scale is thirty nodes for example III, there are a total of 197 edges. Assume that the transportation task requires a 3PL supplier with a load capacity of 100 and a credibility requirement of 8.Table 3
Basic data of nodes for example I.
NodeTiCi135234344423534652745Table 4
Basic data of edges for example I.
Node (i)Node (j)Edge (k)TijkSTijkLCijkSCijkLQijkDijk12112.713.411.412.280812219.821.314.715.51201012319.820.4131570913121.623.213.714.3901013216.517.216.918.61201113317.818.518.419.56091349.210.411.413.2801023110.811.312.813.7100142329.710.311.512.5701024112.513.410.511.3901024210.511.29.911.61101225111.712.216.817.51401125210.711.214.915.67083419.710.38.69.5110934211.812.511.812.680934312.713.310.411.7901036112.713.418.519.790936212.913.810.111.780103639.910.58.79.670113641718.86.47.2601545111.812.37.88.970124529.810.113.515.5110745312.813.211.812.6901146113.715.314.51612094628.99.47.99.5110747112.513.211.412.2501447210.411.310.211.3701047311.612.313.815.590857112.513.513.415.1801257212.613.612.713.1701557311.712.114.815.79013671910.59.910.1110106721718.56.47909


## 5.2. Problem Analysis and Discussion

### 5.2.1. Analysis of Model Parameters

This section will first take the example I as an example to analyze the influence of the model parametersT0,C0,, φ1, and φ2 on the selection of transportation plans. Then, obtain the optimal combination parameters, the optimal route, and its utility. The utility is represented by U on different examples.(1) The Impact of T0 on the Optimal Route. First, set C0=65, ω1=0.4, and ω2=0.6. 1 ⟶ (10.8, 30.2) ⟶ 2 indicates that the route passes through v1 and v2, selects 3PL 1 to complete 80% transportation tasks, and selects 3PL 3 to complete the rest 20% transportation task between the two nodes. And R represents the optimal route. For different value of T0, the result is shown in Table 5.Table 5
The impact ofT0 on the transportation plans.
T0RUCsCLTSTL351 ⟶ (20.5, 30.1, 40.4) ⟶ 3 ⟶ (10.5, 20.1, 30.4) ⟶ 4 ⟶ (20.4, 30.6) ⟶ 7−0.3152.355.155.258.8401 ⟶ (20.3, 30.1, 40.6) ⟶ 3 ⟶ (20.8, 30.2) ⟶ 4 ⟶ (10.4, 20.2, 30.4) ⟶ 7−0.1952.054.354.658.2451 ⟶ (30.5, 40.5) ⟶ 3 ⟶ (11) ⟶ 4 ⟶ (20.4, 30.6) ⟶ 7−0.0753.455.354.858.5501 ⟶ (20.6, 40.4) ⟶ 3 ⟶ (11) ⟶ 4 ⟶ (20.5, 30.5) ⟶ 70.0554.956.953.357.6551 ⟶ (21) ⟶ 3 ⟶ (10.5, 20.5) ⟶ 6 ⟶ (10.3, 20.7) ⟶ 70.1856.858.753.354.0601 ⟶ (20.5, 30.1, 40.4) ⟶ 3 ⟶ (10.5, 20.1, 30.4) ⟶ 4 ⟶ (20.4, 30.6) ⟶ 70.2252.355.155.258.8651 ⟶ (21) ⟶ 3 ⟶ (10.7, 20.1, 30.2) ⟶ 4 ⟶ (20.2, 30.8) ⟶ 70.2855.057.555.658.7701 ⟶ (21) ⟶ 3 ⟶ (10.8, 30.2) ⟶ 4 ⟶ (20.4, 30.6) ⟶ 70.3256.858.553.758.7751 ⟶ (21) ⟶ 3 ⟶ (10.7, 20.1, 30.2) ⟶ 4 ⟶ (20.4, 30.6) ⟶ 70.3655.457.955.658.7801 ⟶ (21) ⟶ 3 ⟶ (20.1, 30.2, 40.5) ⟶ 6 ⟶ (10.4, 20.6) ⟶ 70.4055.857.854.156.8It can be seen from Table5, a series of optimal routes are obtained with variance of T0, and the insight are given as follows. (1) Different values T0 generate different optimal routes. When T0  = 50, the optimal route that satisfy the customer’s requirements for T(R) and C(R) is found for the first time. (2) It can be seen from the grey part of Table 5 that C(R) is extended from [53.4, 55.3] to [54.9, 56.9], and T(R) is reduced from [54.8, 58.5] to [53.3, 57.6]. The optimal route that satisfies the customer’s requirements is obtained. When only considering the tardiness of the route, the customer is willing to choose a route with relatively higher cost, where both T(R) and C(R) meets the requirements. It also shows that the customer is risk aversion and has behavioral characteristics of against loss. (3) If the customer’s requirements for T(R) become looser, the utility of objective function becomes the main factor, and it affects the optimal route.(2) The Impact of C0 on the Optimal Route. Set T0= 65, ω1=0.4, and ω2=0.6. Change the value of C0, and the results are shown in Table 6 and the findings are as follows. (1) Different values C0 generate different optimal routes. When C0  = 40, the optimal route that satisfies the customer’s requirements for T(R) and C(R) is found. (2) It can be seen from the grey part of the table that T(R) is extended from [53.3, 57.6] to [55.0, 58.6], and C(R) is reduced from [54.9,56.9] to [50.8,54.8], and the optimal route meeting the customer’s requirements is found. If the cost of route is only in overspending, the customer is willing to choose a route with relatively longer time, and both T(R) and C(R) satisfy the requirements. This shows that customer is risk aversion and has behavioral characteristics of loss avoidance. (3) With the increase of C0, the customer’s requirements for the optimal route transport cost are looser, and the objective function utility value becomes the main factor determining the optimal route.Table 6
The impact ofC0 on the transportation plans.
C0RUCsCLTSTL401 ⟶ (21) ⟶ 3 ⟶ (10.7, 20.1, 30.2) ⟶ 4 ⟶ (20.2, 30.8) ⟶ 7−0.2355.057.555.658.7451 ⟶ (21) ⟶ 3 ⟶ (20.3, 30.2, 40.5) ⟶ 6 ⟶ (10.4, 20.6) ⟶ 7−0.1255.857.854.156.8501 ⟶ (20.6, 40.4) ⟶ 3 ⟶ (11) ⟶ 4 ⟶ (20.5, 30.5) ⟶ 7−0.0154.956.953.357.6551 ⟶ (20.1, 30.2, 40.7) ⟶ 3 ⟶ (20.7, 30.3) ⟶ 4 ⟶ (10.2, 20.1, 30.7) ⟶ 70.1450.854.855.058.6601 ⟶ (30.5, 40.5) ⟶ 3 ⟶ (11) ⟶ 4 ⟶ (20.4, 30.6) ⟶ 70.2153.455.354.858.5651 ⟶ (21) ⟶ 3 ⟶ (10.7, 20.1, 30.2) ⟶ 4 ⟶ (20.2, 30.8) ⟶ 70.2855.057.555.658.7701 ⟶ (20.6, 40.4) ⟶ 3 ⟶ (11) ⟶ 4 ⟶ (20.5, 30.5) ⟶ 70.3354.956.953.357.6751 ⟶ (21) ⟶ 3 ⟶ (20.3, 30.2, 40.5) ⟶ 6 ⟶ (10.4, 20.6) ⟶ 70.3755.857.854.156.8801 ⟶ (21) ⟶ 3 ⟶ (10.7, 20.1, 30.2) ⟶ 4 ⟶ (20.2, 30.8) ⟶ 70.4255.057.555.658.7851 ⟶ (21) ⟶ 3 ⟶ (10.8, 30.2) ⟶ 4 ⟶ (20.3, 30.7) ⟶ 70.4755.858.353.758.7(3) The Impact of φ1, φ2 on the Optimal Route. Set T0=65 and C0=65, and change the value of φ1 and φ2, the results are shown in Table 7. The best parameters combination for the example I is T0= 65, C0= 65, ω1=0.4, and ω2=0.6. The utility value is 0.28. The best route is 1 ⟶ (20.4, 30.5, 40.1) ⟶ 3 ⟶ (10.5, 20.1, 30.4) ⟶ 4 ⟶ (20.4, 30.6) ⟶ 7.Table 7
The impact ofφ1 and φ2 on the transportation plans.
φ1φ2RUCSCLTSTL0.10.91 ⟶ (30.6, 40.4) ⟶ 3 ⟶ (10.1, 20.3, 30.6) ⟶ 4 ⟶ (20.3, 30.7) ⟶ 70.2848.651.852.857.40.20.81 ⟶ (20.2, 30.2, 40.6) ⟶ 3 ⟶ (20.4, 30.6) ⟶ 4 ⟶ (20.5, 30.5) ⟶ 70.2749.652.853.657.20.30.71 ⟶ (20.5, 30.1, 40.4) ⟶ 3 ⟶ (10.2, 20.5, 30.3) ⟶ 4 ⟶ (20.3, 30.7) ⟶ 70.2750.654.553.457.00.40.61 ⟶ (20.4, 30.5, 40.1) ⟶ 3 ⟶ (10.5, 20.1, 30.4) ⟶ 4 ⟶ (20.4, 30.6) ⟶ 70.2851.954.754.057.60.50.51 ⟶ (21) ⟶ 3 ⟶ (10.8, 30.2) ⟶ 4 ⟶ (20.3, 30.7) ⟶ 70.2855.858.353.758.70.60.41 ⟶ (21) ⟶ 3 ⟶ (20.3, 30.2, 40.5) ⟶ 6 ⟶ (10.4, 20.6) ⟶ 70.2755.857.854.156.80.70.31 ⟶ (21) ⟶ 3 ⟶ (20.3, 30.2, 40.5) ⟶ 6 ⟶ (10.4, 20.6) ⟶ 70.2855.857.854.156.80.80.21 ⟶ (21) ⟶ 3 ⟶ (20.3, 30.2, 40.5) ⟶ 6 ⟶ (10.4, 20.6) ⟶ 70.2855.857.854.156.80.90.11 ⟶ (21) ⟶ 3 ⟶ (10.5, 20.5) ⟶ 6 ⟶ (10.3, 20.7) ⟶ 70.2756.858.753.354.0Through the above analysis, the seven-node instance parameter combination, the optimal utility value, and its optimal path based on ACS-IGWO are obtained. Through the same method, we can get the fifteen and thirty-node case parameter combinations, as shown in Table8. The optimal path of example I is shown in Figure 5:Table 8
Model parameter combination and the best utility of different node cases.
NodeT0C0φ1φ2Utility765650.40.60.2815901200.40.60.30301501800.40.60.30


### 5.2.2. Customer’s Risk Attitude Analysis

This section analyses the influence of customer’s behavior characteristics on the optimal route selection. In example I, ifλ=1 and υ=1 in equations (7) and (8), the mathematical model is based on RUT that is in line with the situation when the customer is risk-neutral attitude.The ACS-IGWO runs 100 times, and the average utility is obtained, which is organized in Tables9 and 10, respectively. Figure 6 is drawn based on the information of the two tables. RU represents the average utility of the customer based on RUT, and EUM represents the average utility of the customer based on EUT.Table 9
Comparison of RU and EUM (T0 = 65).
C020253035404550556065707580EUM−0.30−0.22−0.17−0.12−0.07−0.020.030.080.130.180.230.280.33RU−0.67−0.56−0.44−0.31−0.19−0.070.050.140.210.280.330.380.43Table 10
Comparison of RU and EUM (T0 = 65).
C0859095100105110115120125130135140145150EUM0.380.430.480.530.580.630.680.730.780.830.880.930.981.03RU0.480.520.560.600.640.670.700.730.750.770.790.800.810.82Figure 6
Comparison of RU and EUM (T0  = 65).The customer has higher requirements for transport cost whenC0 is small, so the two models take a negative value, and the satisfaction of the routes is lower, which means the customer is in a loss state. As C0 increases gradually, the values of the two models are increasing gradually, and customer is more satisfied with the routes, which means the customer transforms from the loss state to the gain state. That is, let the “0” utility be the dividing line; if the customer’s utility is negative, the customer is in a loss state; and if it is positive, the customer is in a gain state. From Figure 6, the insights are as follows:(1)
If the utility is negative, the customer is in a loss state. Compared to EUM, the utility of RU is smaller, which points out that customer is risk aversion. The obtained optimal route will not delay, but it overspends seriously. Customer is risk averse to overspending, so they are less satisfied with the optimal routes.(2)
If the utility is positive andC0≤55, the customer is at the beginning of the gain state. Compared to EUM, RU is smaller, which points out that that customer is risk aversion. In Table 6, when C0  = 55, C(R) has the maximum value 54.8, which is almost equal to C0; in this perspective, the optimal routes will not overspend. When C0<55, in Table 6, it can be found that the optimal routes still overspend and do not meet customer’s expectations.(3)
If55≤C0≤120, the customer is in the continuous gain state. Compared to EUM, RU is larger, which points out that that customer is risk aversion. The optimal routes obtained neither overspend nor delay. Compared with EUM, the RU value is larger, which indicates that the customer is more conservative and risk aversion, and expect for higher utility value.(4)
IfC0≥120, the customer is in the “final” gain state. Compared to EUM, RU is smaller, which points out that that the customer is risk aversion. Although there are no delays and overspend, the rate of increase in customer satisfaction with the transportation plan is gradually decreasing. Therefore, customers are not willing to excessively reduce their cost requirements, resulting in unnecessary cost waste.(5)
Overall, the RU curve of Figure6 shows that the slope of the loss state is higher than the gain state, which indicates that the customer is more sensitive to the loss than the gain.WhenC0 is fixed and C0=65, the trend of the result is the same as in Figure 6, so the analysis will not be repeated here. With the customer’s requirements for transport time being gradually reduced, the customer’s risk attitude is the same as the situation that the customer’s transport cost requirements are gradually reduced. In summary, the model can accurately describe the psychological characteristics of people. It also shows that the consideration and analysis of customer behavior characteristics in 4PLRP have a certain practical value and practical significance.


### 5.2.3. Comparative Analysis of Transportation Modes

In this section, two transportation modes are considered: single 3PL and multiple 3PLs undertake transport tasks. The impact of different transportation modes in the cases of simple tardiness, simple overspend, no tardiness, and overspend is studied. Among them, simple tardiness, simple overspend, and definitely no tardiness and overspend are represented byF, O, and N, and single 3PL and multiple 3PL are represented by A and M. When T0 = 65, C0 = 65, φ1=0.4, and φ2=0.6, the results are listed in Table 11.Table 11
Transportation plans with different transportation modes.
SituationModeUCSCLTSTLFA−0.2256.858.154.758.0M−0.1952.054.354.658.2OA−0.2557.759.156.058.5M−0.2355.057.555.658.7NA0.2558.962.155.459.0M0.2855.057.555.658.7In Table10, the following conclusions are obtained. In three cases, transportation cost of a single 3PL mode is higher than multiple 3PLs mode, and the time of both is almost equal. A higher transportation cost corresponds to a lower the utility value, that is, the customer’s satisfaction with the optimal route is lower.


## 5.2.1. Analysis of Model Parameters

This section will first take the example I as an example to analyze the influence of the model parametersT0,C0,, φ1, and φ2 on the selection of transportation plans. Then, obtain the optimal combination parameters, the optimal route, and its utility. The utility is represented by U on different examples.(1) The Impact of T0 on the Optimal Route. First, set C0=65, ω1=0.4, and ω2=0.6. 1 ⟶ (10.8, 30.2) ⟶ 2 indicates that the route passes through v1 and v2, selects 3PL 1 to complete 80% transportation tasks, and selects 3PL 3 to complete the rest 20% transportation task between the two nodes. And R represents the optimal route. For different value of T0, the result is shown in Table 5.Table 5
The impact ofT0 on the transportation plans.
T0RUCsCLTSTL351 ⟶ (20.5, 30.1, 40.4) ⟶ 3 ⟶ (10.5, 20.1, 30.4) ⟶ 4 ⟶ (20.4, 30.6) ⟶ 7−0.3152.355.155.258.8401 ⟶ (20.3, 30.1, 40.6) ⟶ 3 ⟶ (20.8, 30.2) ⟶ 4 ⟶ (10.4, 20.2, 30.4) ⟶ 7−0.1952.054.354.658.2451 ⟶ (30.5, 40.5) ⟶ 3 ⟶ (11) ⟶ 4 ⟶ (20.4, 30.6) ⟶ 7−0.0753.455.354.858.5501 ⟶ (20.6, 40.4) ⟶ 3 ⟶ (11) ⟶ 4 ⟶ (20.5, 30.5) ⟶ 70.0554.956.953.357.6551 ⟶ (21) ⟶ 3 ⟶ (10.5, 20.5) ⟶ 6 ⟶ (10.3, 20.7) ⟶ 70.1856.858.753.354.0601 ⟶ (20.5, 30.1, 40.4) ⟶ 3 ⟶ (10.5, 20.1, 30.4) ⟶ 4 ⟶ (20.4, 30.6) ⟶ 70.2252.355.155.258.8651 ⟶ (21) ⟶ 3 ⟶ (10.7, 20.1, 30.2) ⟶ 4 ⟶ (20.2, 30.8) ⟶ 70.2855.057.555.658.7701 ⟶ (21) ⟶ 3 ⟶ (10.8, 30.2) ⟶ 4 ⟶ (20.4, 30.6) ⟶ 70.3256.858.553.758.7751 ⟶ (21) ⟶ 3 ⟶ (10.7, 20.1, 30.2) ⟶ 4 ⟶ (20.4, 30.6) ⟶ 70.3655.457.955.658.7801 ⟶ (21) ⟶ 3 ⟶ (20.1, 30.2, 40.5) ⟶ 6 ⟶ (10.4, 20.6) ⟶ 70.4055.857.854.156.8It can be seen from Table5, a series of optimal routes are obtained with variance of T0, and the insight are given as follows. (1) Different values T0 generate different optimal routes. When T0  = 50, the optimal route that satisfy the customer’s requirements for T(R) and C(R) is found for the first time. (2) It can be seen from the grey part of Table 5 that C(R) is extended from [53.4, 55.3] to [54.9, 56.9], and T(R) is reduced from [54.8, 58.5] to [53.3, 57.6]. The optimal route that satisfies the customer’s requirements is obtained. When only considering the tardiness of the route, the customer is willing to choose a route with relatively higher cost, where both T(R) and C(R) meets the requirements. It also shows that the customer is risk aversion and has behavioral characteristics of against loss. (3) If the customer’s requirements for T(R) become looser, the utility of objective function becomes the main factor, and it affects the optimal route.(2) The Impact of C0 on the Optimal Route. Set T0= 65, ω1=0.4, and ω2=0.6. Change the value of C0, and the results are shown in Table 6 and the findings are as follows. (1) Different values C0 generate different optimal routes. When C0  = 40, the optimal route that satisfies the customer’s requirements for T(R) and C(R) is found. (2) It can be seen from the grey part of the table that T(R) is extended from [53.3, 57.6] to [55.0, 58.6], and C(R) is reduced from [54.9,56.9] to [50.8,54.8], and the optimal route meeting the customer’s requirements is found. If the cost of route is only in overspending, the customer is willing to choose a route with relatively longer time, and both T(R) and C(R) satisfy the requirements. This shows that customer is risk aversion and has behavioral characteristics of loss avoidance. (3) With the increase of C0, the customer’s requirements for the optimal route transport cost are looser, and the objective function utility value becomes the main factor determining the optimal route.Table 6
The impact ofC0 on the transportation plans.
C0RUCsCLTSTL401 ⟶ (21) ⟶ 3 ⟶ (10.7, 20.1, 30.2) ⟶ 4 ⟶ (20.2, 30.8) ⟶ 7−0.2355.057.555.658.7451 ⟶ (21) ⟶ 3 ⟶ (20.3, 30.2, 40.5) ⟶ 6 ⟶ (10.4, 20.6) ⟶ 7−0.1255.857.854.156.8501 ⟶ (20.6, 40.4) ⟶ 3 ⟶ (11) ⟶ 4 ⟶ (20.5, 30.5) ⟶ 7−0.0154.956.953.357.6551 ⟶ (20.1, 30.2, 40.7) ⟶ 3 ⟶ (20.7, 30.3) ⟶ 4 ⟶ (10.2, 20.1, 30.7) ⟶ 70.1450.854.855.058.6601 ⟶ (30.5, 40.5) ⟶ 3 ⟶ (11) ⟶ 4 ⟶ (20.4, 30.6) ⟶ 70.2153.455.354.858.5651 ⟶ (21) ⟶ 3 ⟶ (10.7, 20.1, 30.2) ⟶ 4 ⟶ (20.2, 30.8) ⟶ 70.2855.057.555.658.7701 ⟶ (20.6, 40.4) ⟶ 3 ⟶ (11) ⟶ 4 ⟶ (20.5, 30.5) ⟶ 70.3354.956.953.357.6751 ⟶ (21) ⟶ 3 ⟶ (20.3, 30.2, 40.5) ⟶ 6 ⟶ (10.4, 20.6) ⟶ 70.3755.857.854.156.8801 ⟶ (21) ⟶ 3 ⟶ (10.7, 20.1, 30.2) ⟶ 4 ⟶ (20.2, 30.8) ⟶ 70.4255.057.555.658.7851 ⟶ (21) ⟶ 3 ⟶ (10.8, 30.2) ⟶ 4 ⟶ (20.3, 30.7) ⟶ 70.4755.858.353.758.7(3) The Impact of φ1, φ2 on the Optimal Route. Set T0=65 and C0=65, and change the value of φ1 and φ2, the results are shown in Table 7. The best parameters combination for the example I is T0= 65, C0= 65, ω1=0.4, and ω2=0.6. The utility value is 0.28. The best route is 1 ⟶ (20.4, 30.5, 40.1) ⟶ 3 ⟶ (10.5, 20.1, 30.4) ⟶ 4 ⟶ (20.4, 30.6) ⟶ 7.Table 7
The impact ofφ1 and φ2 on the transportation plans.
φ1φ2RUCSCLTSTL0.10.91 ⟶ (30.6, 40.4) ⟶ 3 ⟶ (10.1, 20.3, 30.6) ⟶ 4 ⟶ (20.3, 30.7) ⟶ 70.2848.651.852.857.40.20.81 ⟶ (20.2, 30.2, 40.6) ⟶ 3 ⟶ (20.4, 30.6) ⟶ 4 ⟶ (20.5, 30.5) ⟶ 70.2749.652.853.657.20.30.71 ⟶ (20.5, 30.1, 40.4) ⟶ 3 ⟶ (10.2, 20.5, 30.3) ⟶ 4 ⟶ (20.3, 30.7) ⟶ 70.2750.654.553.457.00.40.61 ⟶ (20.4, 30.5, 40.1) ⟶ 3 ⟶ (10.5, 20.1, 30.4) ⟶ 4 ⟶ (20.4, 30.6) ⟶ 70.2851.954.754.057.60.50.51 ⟶ (21) ⟶ 3 ⟶ (10.8, 30.2) ⟶ 4 ⟶ (20.3, 30.7) ⟶ 70.2855.858.353.758.70.60.41 ⟶ (21) ⟶ 3 ⟶ (20.3, 30.2, 40.5) ⟶ 6 ⟶ (10.4, 20.6) ⟶ 70.2755.857.854.156.80.70.31 ⟶ (21) ⟶ 3 ⟶ (20.3, 30.2, 40.5) ⟶ 6 ⟶ (10.4, 20.6) ⟶ 70.2855.857.854.156.80.80.21 ⟶ (21) ⟶ 3 ⟶ (20.3, 30.2, 40.5) ⟶ 6 ⟶ (10.4, 20.6) ⟶ 70.2855.857.854.156.80.90.11 ⟶ (21) ⟶ 3 ⟶ (10.5, 20.5) ⟶ 6 ⟶ (10.3, 20.7) ⟶ 70.2756.858.753.354.0Through the above analysis, the seven-node instance parameter combination, the optimal utility value, and its optimal path based on ACS-IGWO are obtained. Through the same method, we can get the fifteen and thirty-node case parameter combinations, as shown in Table8. The optimal path of example I is shown in Figure 5:Table 8
Model parameter combination and the best utility of different node cases.
NodeT0C0φ1φ2Utility765650.40.60.2815901200.40.60.30301501800.40.60.30


## 5.2.2. Customer’s Risk Attitude Analysis

This section analyses the influence of customer’s behavior characteristics on the optimal route selection. In example I, ifλ=1 and υ=1 in equations (7) and (8), the mathematical model is based on RUT that is in line with the situation when the customer is risk-neutral attitude.The ACS-IGWO runs 100 times, and the average utility is obtained, which is organized in Tables9 and 10, respectively. Figure 6 is drawn based on the information of the two tables. RU represents the average utility of the customer based on RUT, and EUM represents the average utility of the customer based on EUT.Table 9
Comparison of RU and EUM (T0 = 65).
C020253035404550556065707580EUM−0.30−0.22−0.17−0.12−0.07−0.020.030.080.130.180.230.280.33RU−0.67−0.56−0.44−0.31−0.19−0.070.050.140.210.280.330.380.43Table 10
Comparison of RU and EUM (T0 = 65).
C0859095100105110115120125130135140145150EUM0.380.430.480.530.580.630.680.730.780.830.880.930.981.03RU0.480.520.560.600.640.670.700.730.750.770.790.800.810.82Figure 6
Comparison of RU and EUM (T0  = 65).The customer has higher requirements for transport cost whenC0 is small, so the two models take a negative value, and the satisfaction of the routes is lower, which means the customer is in a loss state. As C0 increases gradually, the values of the two models are increasing gradually, and customer is more satisfied with the routes, which means the customer transforms from the loss state to the gain state. That is, let the “0” utility be the dividing line; if the customer’s utility is negative, the customer is in a loss state; and if it is positive, the customer is in a gain state. From Figure 6, the insights are as follows:(1)
If the utility is negative, the customer is in a loss state. Compared to EUM, the utility of RU is smaller, which points out that customer is risk aversion. The obtained optimal route will not delay, but it overspends seriously. Customer is risk averse to overspending, so they are less satisfied with the optimal routes.(2)
If the utility is positive andC0≤55, the customer is at the beginning of the gain state. Compared to EUM, RU is smaller, which points out that that customer is risk aversion. In Table 6, when C0  = 55, C(R) has the maximum value 54.8, which is almost equal to C0; in this perspective, the optimal routes will not overspend. When C0<55, in Table 6, it can be found that the optimal routes still overspend and do not meet customer’s expectations.(3)
If55≤C0≤120, the customer is in the continuous gain state. Compared to EUM, RU is larger, which points out that that customer is risk aversion. The optimal routes obtained neither overspend nor delay. Compared with EUM, the RU value is larger, which indicates that the customer is more conservative and risk aversion, and expect for higher utility value.(4)
IfC0≥120, the customer is in the “final” gain state. Compared to EUM, RU is smaller, which points out that that the customer is risk aversion. Although there are no delays and overspend, the rate of increase in customer satisfaction with the transportation plan is gradually decreasing. Therefore, customers are not willing to excessively reduce their cost requirements, resulting in unnecessary cost waste.(5)
Overall, the RU curve of Figure6 shows that the slope of the loss state is higher than the gain state, which indicates that the customer is more sensitive to the loss than the gain.WhenC0 is fixed and C0=65, the trend of the result is the same as in Figure 6, so the analysis will not be repeated here. With the customer’s requirements for transport time being gradually reduced, the customer’s risk attitude is the same as the situation that the customer’s transport cost requirements are gradually reduced. In summary, the model can accurately describe the psychological characteristics of people. It also shows that the consideration and analysis of customer behavior characteristics in 4PLRP have a certain practical value and practical significance.


## 5.2.3. Comparative Analysis of Transportation Modes

In this section, two transportation modes are considered: single 3PL and multiple 3PLs undertake transport tasks. The impact of different transportation modes in the cases of simple tardiness, simple overspend, no tardiness, and overspend is studied. Among them, simple tardiness, simple overspend, and definitely no tardiness and overspend are represented byF, O, and N, and single 3PL and multiple 3PL are represented by A and M. When T0 = 65, C0 = 65, φ1=0.4, and φ2=0.6, the results are listed in Table 11.Table 11
Transportation plans with different transportation modes.
SituationModeUCSCLTSTLFA−0.2256.858.154.758.0M−0.1952.054.354.658.2OA−0.2557.759.156.058.5M−0.2355.057.555.658.7NA0.2558.962.155.459.0M0.2855.057.555.658.7In Table10, the following conclusions are obtained. In three cases, transportation cost of a single 3PL mode is higher than multiple 3PLs mode, and the time of both is almost equal. A higher transportation cost corresponds to a lower the utility value, that is, the customer’s satisfaction with the optimal route is lower.


## 5.3. Algorithm Analysis and Discussion

In order to test the effectiveness of the proposed algorithm, first, the optimal parameters combination of ACS-IGWO and the optimal results are obtained. Then, the effectiveness of ACS-IGWO is verified by comparing with ACS, D-ACS (ant colony algorithm for heuristic dynamic pheromone update strategy) [26], and ACS-GWO.


### 5.3.1. Parameters Analysis of ACS-IGWO

In this study, the Taguchi method [27] is used to obtain the optimal parameters combination of the four algorithms.In this section, the Taguchi method is used to obtain the optimal parameters combination of ACS-IGWO. Parameters such asα, β,Q, and ρ are selected to be tested. The application software Minitab 17 is introduced to perform parameter testing on ACS-IGWO for example I. The orthogonal test table is shown in Table 12. The signal-to-noise ratio main effect diagram and the mean main effect diagram are shown in Figures 7 and 8, respectively.Table 12
The orthogonal test table.
Running numberαβQρU122500.60.2800222500.60.2800322500.60.2800423800.70.2800523800.70.2801623800.70.28007251000.80.27998251000.80.27999251000.80.27991032800.80.28001132800.80.28001232800.80.280113331000.60.280114331000.60.279915331000.60.28001635500.70.28011735500.70.27991835500.70.279919521000.70.280020521000.70.280021521000.70.28002253500.80.27992353500.80.28002453500.80.27992555800.60.28012655800.60.28002755800.60.2800Figure 7
Signal-to-noise ratio main effect diagram.Figure 8
Mean main effect diagram.From Figures7 and 8, the average value of the maximum signal-to-noise ratio of the four parameters is −11.0578, −11.0567, −11.0562, and −11.0566, respectively. To make the means of mean of the parameters the largest, α, β, Q, and ρ should take level one. Therefore, α=3, β=5, Q=50, ρ=0.6 is the parameters combination of ACS-IGWO. Other parameters are NG, NP, and q0, which have a fixed value in each algorithm. The parameters combination of the four algorithms is shown in Table 13.Table 13
Parameters combination of the four algorithms.
ExampleAlgorithmαβQρNGNPq0IACS22500.730300.7D-ACS22500.63030—ACS-GWO22500.730300.7ACS-IGWO22500.730300.7IIACS231000.7100500.7D-ACS231000.610050—ACS-GWO231000.7100500.7ACS-IGWO231000.7100500.7IIIACS231000.7150800.7D-ACS231000.615080—ACS-GWO231000.7150800.7ACS-IGWO231000.7150800.7


### 5.3.2. Comparative Analysis of the Four Algorithms

The four algorithms are run 100 times. To test the performance of the algorithms, the relevant performance parameters of the algorithm are defined. The average maximum utility is best, the average minimum utility is worst, the average utility is mean, the variance of utility isS, and the average running time is time (s). For the three examples, the study is solved by ACS, D-ACS, ACS-GWO, and ACS-IGWO. The results are shown in Table 14.Table 14
Comparison of algorithms of different node cases.
ExampleAlgorithmBestWorstMeanSIACS-IGWO0.28−0.61−0.010.01ACS-GWO0.28−0.62−0.010.03ACS0.28−0.62−0.030.04D-ACS0.28−0.61−0.040.05IIACS-IGWO0.30−0.51−0.020.12ACS-GWO0.30−0.51−0.040.23ACS0.28−0.67−0.020.85D-ACS0.30−0.53−0.050.14IIIACS-IGWO0.30−0.71−0.010.12ACS-GWO0.30−0.73−0.010.32ACS0.25−0.82−0.061.76D-ACS0.30−0.73−0.030.30In Table13, for example I, all algorithms can find the global optimal solution. While for example II and example III, ACS cannot find the global optimal solution and perform worse with the increase of problem scale. We can find that ACS-IGWO performs best among the four algorithms. From the compare of variance value and mean average utility of algorithms in each example, it can be found that ACS-IGWO has more stable convergence interval.In addition, for the three examples, the optimal results are collected and represented by box plot, which is shown in Figures9–11. For example I, in Figure 9, although the four boxes are different in size, the difference of the concentration range is almost zero, which indicates that the four algorithms can find the optimal route accurately. For Figure 10, example II, the optimal results of ACS-IGWO, ACS-GWO, and D-ACS are better and more reliable than ACS, which indicates that ACS-IGWO, ACS-GWO, and D-ACS have strong ability to skip out of local optimum. The solution effects of ACS-GWO and D-ACS are similar. For Figures 10 and 11, example II and example III, the median line of ACS-IGWO is closer to 0.3, and the upper and lower limits are smaller, which shows that ACS-IGWO can find more reliable results. Its feasibility and effectiveness for solving this kind of routing optimization problem are verfied.Figure 9
Box plot of example I.Figure 10
Box plot of example II.Figure 11
Box plot of example III.


### 5.3.3. Comparative Analysis of Efficiency of the Four Algorithms

The problem size and running time of each algorithm in cases are shown in Table15. In this problem, the transportation path of the goods from the start node to the end point will be planned, and the distribution ratio of 3PL between every two nodes in the path will be obtained. Between the two nodes, the number of 3PL is between 2 and 4, and the proportion of goods allocated by 3PL is between 0 and 1. This ratio needs to be obtained through optimization. It can be found that the scale of the problem is very large, and there are even an infinite variety of results. Therefore, we conservatively estimated the scale of the problem and stipulated that two 3PLs were taken between every two nodes, and the quantity of goods was distributed 1 : 1. For example 1, there are at least 124 transportation possibilities; for example 2, there are at least 761 transportation possibilities; for example 3, there are at least 32,775 transportation possibilities. In Table 15, the rate is the ratio of time to problem size.Table 15
Comparison of algorithms of different node cases.
ExampleProblem sizeAlgorithmTime (s)RateI124ACS-IGWO0.020.000161ACS-GWO0.040.000323ACS0.10.000807D-ACS0.050.000403II761ACS-IGWO0.30.000394ACS-GWO0.50.000657ACS1.50.001971D-ACS0.40.000526III32775ACS-IGWO1.20.000037ACS-GWO5.50.000168ACS22.50.000686D-ACS1.60.000049The running time of ACS-IGWO, ACS-GWO, and D-ACS is shorter than ACS, and ACS-IGWO is the fastest one. For running time in example III, the convergence speed can be arranged as ACS-IGWO > D-ACS > ACS-GWO > ACS.


## 5.3.1. Parameters Analysis of ACS-IGWO

In this study, the Taguchi method [27] is used to obtain the optimal parameters combination of the four algorithms.In this section, the Taguchi method is used to obtain the optimal parameters combination of ACS-IGWO. Parameters such asα, β,Q, and ρ are selected to be tested. The application software Minitab 17 is introduced to perform parameter testing on ACS-IGWO for example I. The orthogonal test table is shown in Table 12. The signal-to-noise ratio main effect diagram and the mean main effect diagram are shown in Figures 7 and 8, respectively.Table 12
The orthogonal test table.
Running numberαβQρU122500.60.2800222500.60.2800322500.60.2800423800.70.2800523800.70.2801623800.70.28007251000.80.27998251000.80.27999251000.80.27991032800.80.28001132800.80.28001232800.80.280113331000.60.280114331000.60.279915331000.60.28001635500.70.28011735500.70.27991835500.70.279919521000.70.280020521000.70.280021521000.70.28002253500.80.27992353500.80.28002453500.80.27992555800.60.28012655800.60.28002755800.60.2800Figure 7
Signal-to-noise ratio main effect diagram.Figure 8
Mean main effect diagram.From Figures7 and 8, the average value of the maximum signal-to-noise ratio of the four parameters is −11.0578, −11.0567, −11.0562, and −11.0566, respectively. To make the means of mean of the parameters the largest, α, β, Q, and ρ should take level one. Therefore, α=3, β=5, Q=50, ρ=0.6 is the parameters combination of ACS-IGWO. Other parameters are NG, NP, and q0, which have a fixed value in each algorithm. The parameters combination of the four algorithms is shown in Table 13.Table 13
Parameters combination of the four algorithms.
ExampleAlgorithmαβQρNGNPq0IACS22500.730300.7D-ACS22500.63030—ACS-GWO22500.730300.7ACS-IGWO22500.730300.7IIACS231000.7100500.7D-ACS231000.610050—ACS-GWO231000.7100500.7ACS-IGWO231000.7100500.7IIIACS231000.7150800.7D-ACS231000.615080—ACS-GWO231000.7150800.7ACS-IGWO231000.7150800.7


## 5.3.2. Comparative Analysis of the Four Algorithms

The four algorithms are run 100 times. To test the performance of the algorithms, the relevant performance parameters of the algorithm are defined. The average maximum utility is best, the average minimum utility is worst, the average utility is mean, the variance of utility isS, and the average running time is time (s). For the three examples, the study is solved by ACS, D-ACS, ACS-GWO, and ACS-IGWO. The results are shown in Table 14.Table 14
Comparison of algorithms of different node cases.
ExampleAlgorithmBestWorstMeanSIACS-IGWO0.28−0.61−0.010.01ACS-GWO0.28−0.62−0.010.03ACS0.28−0.62−0.030.04D-ACS0.28−0.61−0.040.05IIACS-IGWO0.30−0.51−0.020.12ACS-GWO0.30−0.51−0.040.23ACS0.28−0.67−0.020.85D-ACS0.30−0.53−0.050.14IIIACS-IGWO0.30−0.71−0.010.12ACS-GWO0.30−0.73−0.010.32ACS0.25−0.82−0.061.76D-ACS0.30−0.73−0.030.30In Table13, for example I, all algorithms can find the global optimal solution. While for example II and example III, ACS cannot find the global optimal solution and perform worse with the increase of problem scale. We can find that ACS-IGWO performs best among the four algorithms. From the compare of variance value and mean average utility of algorithms in each example, it can be found that ACS-IGWO has more stable convergence interval.In addition, for the three examples, the optimal results are collected and represented by box plot, which is shown in Figures9–11. For example I, in Figure 9, although the four boxes are different in size, the difference of the concentration range is almost zero, which indicates that the four algorithms can find the optimal route accurately. For Figure 10, example II, the optimal results of ACS-IGWO, ACS-GWO, and D-ACS are better and more reliable than ACS, which indicates that ACS-IGWO, ACS-GWO, and D-ACS have strong ability to skip out of local optimum. The solution effects of ACS-GWO and D-ACS are similar. For Figures 10 and 11, example II and example III, the median line of ACS-IGWO is closer to 0.3, and the upper and lower limits are smaller, which shows that ACS-IGWO can find more reliable results. Its feasibility and effectiveness for solving this kind of routing optimization problem are verfied.Figure 9
Box plot of example I.Figure 10
Box plot of example II.Figure 11
Box plot of example III.


## 5.3.3. Comparative Analysis of Efficiency of the Four Algorithms

The problem size and running time of each algorithm in cases are shown in Table15. In this problem, the transportation path of the goods from the start node to the end point will be planned, and the distribution ratio of 3PL between every two nodes in the path will be obtained. Between the two nodes, the number of 3PL is between 2 and 4, and the proportion of goods allocated by 3PL is between 0 and 1. This ratio needs to be obtained through optimization. It can be found that the scale of the problem is very large, and there are even an infinite variety of results. Therefore, we conservatively estimated the scale of the problem and stipulated that two 3PLs were taken between every two nodes, and the quantity of goods was distributed 1 : 1. For example 1, there are at least 124 transportation possibilities; for example 2, there are at least 761 transportation possibilities; for example 3, there are at least 32,775 transportation possibilities. In Table 15, the rate is the ratio of time to problem size.Table 15
Comparison of algorithms of different node cases.
ExampleProblem sizeAlgorithmTime (s)RateI124ACS-IGWO0.020.000161ACS-GWO0.040.000323ACS0.10.000807D-ACS0.050.000403II761ACS-IGWO0.30.000394ACS-GWO0.50.000657ACS1.50.001971D-ACS0.40.000526III32775ACS-IGWO1.20.000037ACS-GWO5.50.000168ACS22.50.000686D-ACS1.60.000049The running time of ACS-IGWO, ACS-GWO, and D-ACS is shorter than ACS, and ACS-IGWO is the fastest one. For running time in example III, the convergence speed can be arranged as ACS-IGWO > D-ACS > ACS-GWO > ACS.


## 6. Conclusion

This study considers the customer’s tardiness and overspending risk aversion, based on the proportional utility theory to establish a multi-3PL common distribution model. ACS-IGWO is designed to solve the final optimization problem. The analysis of the experimental results shows that when multiple 3PLs complete the transportation task between two nodes, the model proposed in this study can accurately describe the customer’s risk attitude. This article compares ACS-IGWO with ACS-GWO, ACS, and D-ACS to verify the effectiveness of ACS-IGWO. Only single point to single point logistics distribution tasks are considered. In extensive research, scholars can consider a distribution network with multiple supply starting points and multiple demand ending points and establish a many-to-many route optimization model to solve and study.



---

*Source: 9864064-2022-03-18.xml*