# 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
 (2020)

---

## Abstract

The fourth-party logistics routing problem (4PLRP) is an important issue in the operation of fourth-party logistics (4PL). In this paper, the study of fourth-party logistics (4PL) path optimization considers that more third-party logistics (3PL) undertake transportation tasks. Under the condition that the 3PL transportation time, transportation cost, node transit time, and transit cost are uncertain, 4PL provides customers with a set of transportation solutions to transport transportation tasks from the initial node to the destination node according to the customer’s risk aversion preference. The transportation scheme not only meets the customer’s time and cost requirements but also meets the carrying capacity and reputation constraints of 3PL. Between the two nodes, one or more 3PLs will undertake the transportation task. The customer’s risk preference will be measured by the ratio utility theory (RUT). An ant colony system-improved grey wolf optimization (ACS-IGWO) is designed to solve the model, and the grey wolf optimization (GWO) is improved by the convergence factor and the proportional weight. Problem analysis is conducted through simulation experiments.


---

## Body

## 1. Introduction

Fourth-party logistics (4PL) provides the customer with satisfactory supply chain solutions by integrating their own and other resources and capabilities. The 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 logistics with a 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. Based on the reliability theory and multigraph, Li et al. [3] proposed a chance-constrained programming model of 4PL, aiming at the 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 the dynamic time planning in the logistics transportation network, Cui et al. [6] proposed a 4PL model considering the transportation, staying, 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 the 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 the EUT, it is difficult to make a reasonable explanation for the phenomena such as the Allais paradox [8]. The 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, 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, Dorigo and Gambardella [17] proposed an 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, and 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 the 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 searchability of the grey wolf, a GWO algorithm RW-GWO based on a random walk was proposed. The algorithm was subjected to a nonparametric test, Wilcoxon test, and performance index analysis, and the improved algorithm has shown 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 the grey wolf, which combines a 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 an increase in transport cost, and the related research of 4PL routing optimization based on risk preference is still relatively rare. In this paper, a 4PLRP model based on RUT is established to indicate the decision problem when the customer is risk aversion. This paper 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.There are many contributions to our work, which are listed as follows:(1)
In the transportation scheme provided to customers, the situation of multiple 3PL cooperative transportations 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 paper establishes the utility function and 4PLRP model from the perspective of customers’ risk appetite for time and cost(3)
Firstly, 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; secondly, the introduction of proportional-based dynamic weighting strategy improves the iteration method of the next generation of individuals, making the algorithm have powerful adaptability and faster convergence speedThe rest of this paper is organized as follows. Section2 is the problem description and parameter definition. Section 3 describes our 4PLRP mathematical model based on the ratio utility theory. In Section 4, the mechanism of population initialization based on ACS is introduced and the improved 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 on the current transportation network, 3PLs, due date, and the total cost is obtained. And one or more 3PLs can be selected between nodes; 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 the 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 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; it also indicates the label of endpoint.rijThe number of edges (3PLs) between nodevi and node vj.eijkThekth edge between node vi and node vj, where i,j∈1,2,…,n, and k is the index number of edges.TiThe time of nodevi, including the time of processing, inventory, loading, and unloading.CiThe cost of nodevi, including the cost of processing, inventory, loading, and unloading.TijkThe time of the 3PL provider in edgeeijk.CijkThe cost of the 3PL provider in edgeeijk.TRThe transport time of the route.CRThe transport cost of the route.TijkSThe minimum time of the 3PL provider in edgeeijk.TijkLThe maximum time of the 3PL provider in edgeeijk.TijSThe shortest transport time between nodevi and node vj.TijLThe longest transport time between nodevi and node vj.TSThe minimum transport time of the route.TLThe maximum transport time of the route.CijkSThe minimum cost of the 3PL provider in edgeeijk.CijkLThe maximum cost of the 3PL provider in edgeeijk.CijSThe maximum transport cost between nodevi and node vj.CijLThe maximum transport cost between nodevi and node vj.CSThe minimum transport cost of the route.CLThe maximum transport cost of the route.T0The due date required by the customer.C0The total cost required by the customer.φ1Weight of the time.φ2Weight of the cost.QijkTransport capacity of the 3PL provider in edgeeijk.DijkReputation of the 3PL provider in edgeeijk.QThe required transport capacity from the task.ΔQiThe amount of task that transits at the node.ptThe delay probability of transportation tasks.pcThe probability of the overspending of transportation tasks.DThe required reputation required by the customer.RThe route from nodevs to node ve.hijkThe proportion of transport tasks undertaken by edgeeijk.xijWhether the 3PL provider is selected between nodevi and node vj to undertake the transportation task, which is a 0-1 variable.xijkWhether the edgeeijk is selected to undertake the transportation task, which is a 0-1 variable.yiWhether 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 certain amount of tasks assigned to the 3PL under the condition that the capacity of each selected 3PL is satisfied, which is shown in picture 2. Assume that three 3PLs are selected between nodesv2 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 paper 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 80%T5 and 80%C5, respectively.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 paper 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, Δ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 the change valueΔG to the final value G0+ΔG indicates consumer utility, named RUT as follows: RU=∑iΔGi/Gi+ΔG, if RU > 0, the consumer has a willingness to receive the option ΔG; otherwise, the 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∑k=1rijhijkQQijkCijkSxijk,(4)CijL=∑i=1n∑j=1n∑k=1rijhijkQQijkCijkLxijkTR, 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=1rijhijkQQijkcijkxijk∑i=1nΔQiQCiyi,CL=∑i=1n∑j=1n∑k=1rijhijkQQijkcijkxijk+∑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<CL1,C0<CS..This paper proposes a decision-making method based on RUT for the 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 the 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,1hij1,…,khijk,…,rijhijrij,vj,…,ve.Equation (9) is the objective function, which means the maximum of the total utility, including the utility of transport time and the utility of transport cost; equation (10) represents the task distribution ratio of the 3PL provider, where 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 whether the 3PL undertakes the transportation task between the two nodes; equation (14) represents whether to choose the 3PL between two nodes; equation (15) represents 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 that 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 certain amount of tasks assigned to the 3PL under the condition that the capacity of each selected 3PL is satisfied, which is shown in picture 2. Assume that three 3PLs are selected between nodesv2 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 paper 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 80%T5 and 80%C5, respectively.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 paper 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, Δ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 the change valueΔG to the final value G0+ΔG indicates consumer utility, named RUT as follows: RU=∑iΔGi/Gi+ΔG, if RU > 0, the consumer has a willingness to receive the option ΔG; otherwise, the 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∑k=1rijhijkQQijkCijkSxijk,(4)CijL=∑i=1n∑j=1n∑k=1rijhijkQQijkCijkLxijkTR, 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=1rijhijkQQijkcijkxijk∑i=1nΔQiQCiyi,CL=∑i=1n∑j=1n∑k=1rijhijkQQijkcijkxijk+∑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<CL1,C0<CS..This paper proposes a decision-making method based on RUT for the 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 the 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,1hij1,…,khijk,…,rijhijrij,vj,…,ve.Equation (9) is the objective function, which means the maximum of the total utility, including the utility of transport time and the utility of transport cost; equation (10) represents the task distribution ratio of the 3PL provider, where 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 whether the 3PL undertakes the transportation task between the two nodes; equation (14) represents whether to choose the 3PL between two nodes; equation (15) represents 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 that 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, 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 the disadvantages of low solution accuracy and slow convergence speed. Therefore, this paper 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 paper adopts the integer number encoding method. Firstly, the multigraph is represented by an adjacency list shown in Figure1. If there is no edge between two nodes, the value of the 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.


#### 4.1.2. Transfer Strategy

The ant determines the transfer direction according to the pseudo-random scale rule of equation (19):(19)j=argmaxj∈JsiτijkNgαηijkβ,q≤q0,Eq.26,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 equation (20):(20)j=argτijkNgαηijkβ∑j∈JsiτijkNgαηijkNgβ,j∈Jsi,0,otherwise,where α and β are 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 be provided in advance. Before and after the node, if the 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 a 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 to 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 in each cycle globally was updated:(21)τijkbsNg=1−ρτijkbsNg+ρΔτijkbsNg.Δτijkbs=Q/TRbs+CRbs is the pheromone increment of the current optimal route. TRbs and CRbs are the transport time and transport cost of the current optimal route. ρ and Q are algorithm parameters. ρ is the pheromone volatilization coefficient, ρ∈0,1, and Q is a certain amount of pheromone concentration.


#### 4.1.5. Local Pheromone Update

After all the ants complete a transfer, a local update of the pheromone was performed according to equation (22). τ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 the 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 prey. The mathematical model of prey behavior is described as follows:(23)D=C⋅Xpt−Xt,Xt+1=Xpt−A⋅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 coefficient vectors:(24)a=2−2ttmax,(25)A=2⋅a⋅r2−a,(26)C=2⋅r,1.r1,r2 are random numbers, and 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 and 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 Xt+1 represents the new position of ω:(27)Dα=C1⋅Xαt−XtX1=Xαt−A1Dα,Dβ=C2⋅Xβt−XtX2=Xβt−A2Dβ,Dδ=C3⋅Xδt−XtX3=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 is searched 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 paper adopts a convergence factor based on the change of exponential law [33], and its modified expression is in 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.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 the αω in Figure 5 is the result that ω learns from α wolf. The ω learns β and δ in the same way, and the most effective scheme in αω, βω, and δω is called ωmax. Comparing ωmax with Xt+1, if ωmax≥Xt+1, the transportation plan represented by the wolf to the transportation plan is update; if ωmax≤Xt+1, 3PLs between nodes are rerandomly distributed 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 the 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 ACS.Step 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: updatea by equation (28), and update A and C by equations (25) and (26), respectively.Step 6: iftmax is reached, the calculation is stopped, and 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 paper adopts the integer number encoding method. Firstly, the multigraph is represented by an adjacency list shown in Figure1. If there is no edge between two nodes, the value of the 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.


### 4.1.2. Transfer Strategy

The ant determines the transfer direction according to the pseudo-random scale rule of equation (19):(19)j=argmaxj∈JsiτijkNgαηijkβ,q≤q0,Eq.26,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 equation (20):(20)j=argτijkNgαηijkβ∑j∈JsiτijkNgαηijkNgβ,j∈Jsi,0,otherwise,where α and β are 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 be provided in advance. Before and after the node, if the 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 a 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 to 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 in each cycle globally was updated:(21)τijkbsNg=1−ρτijkbsNg+ρΔτijkbsNg.Δτijkbs=Q/TRbs+CRbs is the pheromone increment of the current optimal route. TRbs and CRbs are the transport time and transport cost of the current optimal route. ρ and Q are algorithm parameters. ρ is the pheromone volatilization coefficient, ρ∈0,1, and Q is a certain amount of pheromone concentration.


### 4.1.5. Local Pheromone Update

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


## 4.1.1. Coding Method

This paper adopts the integer number encoding method. Firstly, the multigraph is represented by an adjacency list shown in Figure1. If there is no edge between two nodes, the value of the 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.


## 4.1.2. Transfer Strategy

The ant determines the transfer direction according to the pseudo-random scale rule of equation (19):(19)j=argmaxj∈JsiτijkNgαηijkβ,q≤q0,Eq.26,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 equation (20):(20)j=argτijkNgαηijkβ∑j∈JsiτijkNgαηijkNgβ,j∈Jsi,0,otherwise,where α and β are 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 be provided in advance. Before and after the node, if the 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 a 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 to 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 in each cycle globally was updated:(21)τijkbsNg=1−ρτijkbsNg+ρΔτijkbsNg.Δτijkbs=Q/TRbs+CRbs is the pheromone increment of the current optimal route. TRbs and CRbs are the transport time and transport cost of the current optimal route. ρ and Q are algorithm parameters. ρ is the pheromone volatilization coefficient, ρ∈0,1, and Q is a certain amount of pheromone concentration.


## 4.1.5. Local Pheromone Update

After all the ants complete a transfer, a local update of the pheromone was performed according to equation (22). τ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 the 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 prey. The mathematical model of prey behavior is described as follows:(23)D=C⋅Xpt−Xt,Xt+1=Xpt−A⋅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 coefficient vectors:(24)a=2−2ttmax,(25)A=2⋅a⋅r2−a,(26)C=2⋅r,1.r1,r2 are random numbers, and 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 and 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 Xt+1 represents the new position of ω:(27)Dα=C1⋅Xαt−XtX1=Xαt−A1Dα,Dβ=C2⋅Xβt−XtX2=Xβt−A2Dβ,Dδ=C3⋅Xδt−XtX3=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 the 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 prey. The mathematical model of prey behavior is described as follows:(23)D=C⋅Xpt−Xt,Xt+1=Xpt−A⋅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 coefficient vectors:(24)a=2−2ttmax,(25)A=2⋅a⋅r2−a,(26)C=2⋅r,1.r1,r2 are random numbers, and 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 and 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 Xt+1 represents the new position of ω:(27)Dα=C1⋅Xαt−XtX1=Xαt−A1Dα,Dβ=C2⋅Xβt−XtX2=Xβt−A2Dβ,Dδ=C3⋅Xδt−XtX3=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 is searched 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 paper adopts a convergence factor based on the change of exponential law [33], and its modified expression is in 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.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 the αω in Figure 5 is the result that ω learns from α wolf. The ω learns β and δ in the same way, and the most effective scheme in αω, βω, and δω is called ωmax. Comparing ωmax with Xt+1, if ωmax≥Xt+1, the transportation plan represented by the wolf to the transportation plan is update; if ωmax≤Xt+1, 3PLs between nodes are rerandomly distributed 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 the 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 ACS.Step 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: updatea by equation (28), and update A and C by equations (25) and (26), respectively.Step 6: iftmax is reached, the calculation is stopped, and 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 the transport time and transport cost. An ACS-IGWO (ant colony system-improved grey wolf optimizer) is applied to solve the 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.Firstly, the material data is described. Secondly, 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 are verified by example tests between algorithms.


### 5.1. Data Description

For example I, the data of nodes and edges are shown in Tables2 and 3, respectively, and 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 2
Basic data of nodes for example I.
NodeTiCi135234344423534652745Table 3
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 example I as an example to analyze the influence of the model parametersT0,C0,φ1, and φ2 on the selection of transportation plans. Then, the optimal combination parameters, the optimal route, and its utility are obtained. The utility is represented by U on different examples.(1) The Impact ofT0 on the Optimal Route. Firstly, 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% of transportation tasks, and selects 3PL 3 to complete the rest 20% transportation task between the two nodes. And the R represents the optimal route. For different values of T0, the result is shown in Table 4.Table 4
The impact ofT0 on 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 Table 4 that a series of optimal routes are obtained with a variance ofT0, and the insight is given as follows. (1) Different values T0 generate different optimal routes. When T0 = 50, the optimal route that satisfies 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 a relatively higher cost, where both T(R) and C(R) meet the requirements. It also shows that the customer is risk aversion and has behavioral characteristics 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.Table 5
Logistics network adjacency list (7-node).
Node123456710340000200222003000304040000323500000036000000270000000(2) The Impact ofC0 on the Optimal Route. Set T0=65, ω1=0.4, and ω2=0.6. Change the value of C0. 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 the route is only in overspending, the customer is willing to choose a route with a relatively longer time, and both T(R) and C(R) satisfy the requirements. This shows that the 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 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φ1andφ2 on the Optimal Route. Set T0=65, C0=65, and change the value of φ1 and φ2. The results are shown in Table 7. The best parameters’ combination for example I is T0=65, C0=65, ω1=0.4, and ω2=0.6. And 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 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 analyzes 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 a 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 the RUT, and EUM represents the average utility of the customer based on the 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 that the customer is in a loss state. As C0 increases gradually, the values of the two models are increasing gradually, and the customer is more satisfied with the routes, which means that 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 the customer is risk aversion. The obtained optimal route will not delay, but it overspends seriously. The customer is risk averse to overspending, so it is 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 the customer is risk aversion. In Table 6, when C0 = 55, C(R) has a maximum value of 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 the 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 expects a 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 overspending, 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 the 4PLRP have 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 overspending, and definitely no tardiness and overspending are represented byF, O, and N, and single 3PL and multiple 3PLs 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 Table11, the following conclusions are obtained. In three cases, the transportation cost of a single 3PL mode is higher than the multiple 3PLs mode, and the time of both is almost equal. A higher transportation cost corresponds to a lower 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, firstly, the optimal parameters’ combination of ACS-IGWO and the optimal results are obtained. And then, the effectiveness of ACS-IGWO is verified by comparing it 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 paper, 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 the ACS-IGWO. Parameters such asα, β,Q, and ρ are selected to be tested. The application software Minitab19 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 numberQαβρU15022Figure0.60.2801250330.70.2800350550.80.2800450220.60.2801550330.70.2799650550.80.2799750220.70.2801850330.80.2799950550.60.279910100220.80.280011100330.60.280012100450.70.280013100220.70.280014100330.80.279915100550.60.280116100220.80.279917100330.60.279918100550.70.2799Figure 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 the mean of the parameters, the largest, α, β, Q, and ρ should take level one. Therefore, α=3, β=5, Q=50, and ρ=0.6 are the parameters’ combination of ACS-IGWO. Other parameters are NG, NP, and q0, which have 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 isBest, the average minimum utility is Worst, the average utility is Mean, the variance of utility is S, and the average running time is Time (s). For the three examples, the paper 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 Table14, 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. And we can find that ACS-IGWO performs best among the four algorithms. From the comparison of variance value and the mean average utility of algorithms in each example, it can be found that ACS-IGWO has a more stable convergence interval.In addition, for the three examples, the optimal results are collected and represented by the 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 example II in Figure 10, 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 the local optimum. The solution effects of ACS-GWO and D-ACS are similar. For examples II and III in Figures 10 and 11, 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 verified.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 the 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 endpoint 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 number 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 that of ACS, and ACS-IGWO is the fastest one. For the 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 Tables2 and 3, respectively, and 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 2
Basic data of nodes for example I.
NodeTiCi135234344423534652745Table 3
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 example I as an example to analyze the influence of the model parametersT0,C0,φ1, and φ2 on the selection of transportation plans. Then, the optimal combination parameters, the optimal route, and its utility are obtained. The utility is represented by U on different examples.(1) The Impact ofT0 on the Optimal Route. Firstly, 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% of transportation tasks, and selects 3PL 3 to complete the rest 20% transportation task between the two nodes. And the R represents the optimal route. For different values of T0, the result is shown in Table 4.Table 4
The impact ofT0 on 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 Table 4 that a series of optimal routes are obtained with a variance ofT0, and the insight is given as follows. (1) Different values T0 generate different optimal routes. When T0 = 50, the optimal route that satisfies 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 a relatively higher cost, where both T(R) and C(R) meet the requirements. It also shows that the customer is risk aversion and has behavioral characteristics 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.Table 5
Logistics network adjacency list (7-node).
Node123456710340000200222003000304040000323500000036000000270000000(2) The Impact ofC0 on the Optimal Route. Set T0=65, ω1=0.4, and ω2=0.6. Change the value of C0. 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 the route is only in overspending, the customer is willing to choose a route with a relatively longer time, and both T(R) and C(R) satisfy the requirements. This shows that the 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 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φ1andφ2 on the Optimal Route. Set T0=65, C0=65, and change the value of φ1 and φ2. The results are shown in Table 7. The best parameters’ combination for example I is T0=65, C0=65, ω1=0.4, and ω2=0.6. And 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 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 analyzes 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 a 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 the RUT, and EUM represents the average utility of the customer based on the 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 that the customer is in a loss state. As C0 increases gradually, the values of the two models are increasing gradually, and the customer is more satisfied with the routes, which means that 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 the customer is risk aversion. The obtained optimal route will not delay, but it overspends seriously. The customer is risk averse to overspending, so it is 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 the customer is risk aversion. In Table 6, when C0 = 55, C(R) has a maximum value of 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 the 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 expects a 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 overspending, 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 the 4PLRP have 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 overspending, and definitely no tardiness and overspending are represented byF, O, and N, and single 3PL and multiple 3PLs 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 Table11, the following conclusions are obtained. In three cases, the transportation cost of a single 3PL mode is higher than the multiple 3PLs mode, and the time of both is almost equal. A higher transportation cost corresponds to a lower 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 example I as an example to analyze the influence of the model parametersT0,C0,φ1, and φ2 on the selection of transportation plans. Then, the optimal combination parameters, the optimal route, and its utility are obtained. The utility is represented by U on different examples.(1) The Impact ofT0 on the Optimal Route. Firstly, 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% of transportation tasks, and selects 3PL 3 to complete the rest 20% transportation task between the two nodes. And the R represents the optimal route. For different values of T0, the result is shown in Table 4.Table 4
The impact ofT0 on 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 Table 4 that a series of optimal routes are obtained with a variance ofT0, and the insight is given as follows. (1) Different values T0 generate different optimal routes. When T0 = 50, the optimal route that satisfies 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 a relatively higher cost, where both T(R) and C(R) meet the requirements. It also shows that the customer is risk aversion and has behavioral characteristics 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.Table 5
Logistics network adjacency list (7-node).
Node123456710340000200222003000304040000323500000036000000270000000(2) The Impact ofC0 on the Optimal Route. Set T0=65, ω1=0.4, and ω2=0.6. Change the value of C0. 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 the route is only in overspending, the customer is willing to choose a route with a relatively longer time, and both T(R) and C(R) satisfy the requirements. This shows that the 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 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φ1andφ2 on the Optimal Route. Set T0=65, C0=65, and change the value of φ1 and φ2. The results are shown in Table 7. The best parameters’ combination for example I is T0=65, C0=65, ω1=0.4, and ω2=0.6. And 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 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 analyzes 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 a 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 the RUT, and EUM represents the average utility of the customer based on the 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 that the customer is in a loss state. As C0 increases gradually, the values of the two models are increasing gradually, and the customer is more satisfied with the routes, which means that 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 the customer is risk aversion. The obtained optimal route will not delay, but it overspends seriously. The customer is risk averse to overspending, so it is 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 the customer is risk aversion. In Table 6, when C0 = 55, C(R) has a maximum value of 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 the 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 expects a 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 overspending, 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 the 4PLRP have 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 overspending, and definitely no tardiness and overspending are represented byF, O, and N, and single 3PL and multiple 3PLs 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 Table11, the following conclusions are obtained. In three cases, the transportation cost of a single 3PL mode is higher than the multiple 3PLs mode, and the time of both is almost equal. A higher transportation cost corresponds to a lower 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, firstly, the optimal parameters’ combination of ACS-IGWO and the optimal results are obtained. And then, the effectiveness of ACS-IGWO is verified by comparing it 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 paper, 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 the ACS-IGWO. Parameters such asα, β,Q, and ρ are selected to be tested. The application software Minitab19 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 numberQαβρU15022Figure0.60.2801250330.70.2800350550.80.2800450220.60.2801550330.70.2799650550.80.2799750220.70.2801850330.80.2799950550.60.279910100220.80.280011100330.60.280012100450.70.280013100220.70.280014100330.80.279915100550.60.280116100220.80.279917100330.60.279918100550.70.2799Figure 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 the mean of the parameters, the largest, α, β, Q, and ρ should take level one. Therefore, α=3, β=5, Q=50, and ρ=0.6 are the parameters’ combination of ACS-IGWO. Other parameters are NG, NP, and q0, which have 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 isBest, the average minimum utility is Worst, the average utility is Mean, the variance of utility is S, and the average running time is Time (s). For the three examples, the paper 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 Table14, 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. And we can find that ACS-IGWO performs best among the four algorithms. From the comparison of variance value and the mean average utility of algorithms in each example, it can be found that ACS-IGWO has a more stable convergence interval.In addition, for the three examples, the optimal results are collected and represented by the 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 example II in Figure 10, 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 the local optimum. The solution effects of ACS-GWO and D-ACS are similar. For examples II and III in Figures 10 and 11, 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 verified.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 the 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 endpoint 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 number 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 that of ACS, and ACS-IGWO is the fastest one. For the 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 paper, 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 the ACS-IGWO. Parameters such asα, β,Q, and ρ are selected to be tested. The application software Minitab19 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 numberQαβρU15022Figure0.60.2801250330.70.2800350550.80.2800450220.60.2801550330.70.2799650550.80.2799750220.70.2801850330.80.2799950550.60.279910100220.80.280011100330.60.280012100450.70.280013100220.70.280014100330.80.279915100550.60.280116100220.80.279917100330.60.279918100550.70.2799Figure 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 the mean of the parameters, the largest, α, β, Q, and ρ should take level one. Therefore, α=3, β=5, Q=50, and ρ=0.6 are the parameters’ combination of ACS-IGWO. Other parameters are NG, NP, and q0, which have 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 isBest, the average minimum utility is Worst, the average utility is Mean, the variance of utility is S, and the average running time is Time (s). For the three examples, the paper 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 Table14, 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. And we can find that ACS-IGWO performs best among the four algorithms. From the comparison of variance value and the mean average utility of algorithms in each example, it can be found that ACS-IGWO has a more stable convergence interval.In addition, for the three examples, the optimal results are collected and represented by the 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 example II in Figure 10, 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 the local optimum. The solution effects of ACS-GWO and D-ACS are similar. For examples II and III in Figures 10 and 11, 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 verified.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 the 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 endpoint 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 number 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 that of ACS, and ACS-IGWO is the fastest one. For the running time in example III, the convergence speed can be arranged as ACS-IGWO > D-ACS > ACS-GWO > ACS.


## 6. Conclusion

This paper 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 paper 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: 8831746-2020-10-15.xml*