4PL routing problem using hybrid beetle swarm optimization

From the view of comprehensive consideration, this considers the uncertainty of transportation time and cost due to seasonality and human factors, and a multi-objective chance constrained programming model with minimum transportation time and cost was established. Due to the complexity of the problem, traditional algorithms are not enough, and more efficient optimization algorithms are still needed. According to the characteristics of the problem, the beetle antennae search (BAS) with one beetle is changed into multiple beetles, and Dijkstra algorithm is embedded, a hybrid beetle swarm optimization algorithm (HBSO) is designed to solve the problem, and the case analysis and algorithm analysis are made for three different examples; the datasets are referenced in literature (Lu et al. in Computer Integrated Manufacturing Systems 5:10–14, 2020). Under the three examples, the influence of the customer on the time and cost is different. The confidence level and weight of cost and time have great influence on the decision results. 4PL provider will not affect the choice of transportation route, but will change the choice of 3PL supplier and transportation mode. HBSO is an effective method to solve this problem. It is proved on efficiency, effectiveness and reliability in experiments by comparing with GA, PSO and BAS.


Introduction
With the rapid development of economy and technology, many companies have outsourced their logistics to special logistics providers known as third party logistics (3PL).3PL supplier integrates various logistics functions and elements and becomes a specialized logistics organization (Berglund et al. 1999).However, most 3PL providers mainly provide transportation and warehousing services.These providers lack cooperation, and their complementary resources have not been effectively utilized.Thus, the integration and management of these 3PL providers to benefit supply chain members and optimize the supply chain presents a difficult and urgent problem, the fourth party logistics (4PL) came into being.4PL provider helps companies to reduce costs by effectively integrating information resources and relies on excellent 3PL suppliers to provide customers with a unique supply chain solution (Gattorna 1998).
With the emergence of 4PL, the study on the fourth party logistics routing problem (4PLRP) has emerged (Cui et al. 2013).Cost, time and loading volume are the main research factors of 4PLRP (Li et al. 2003;Huang et al. 2013;Liu et al. 2014;Demir et al. 2016;Qiu et al. 2018).Among them, cost and time are important factors for different transportation modes.This is because the cost and time are different for different transportation modes.The rapid transportation mode has the characteristics of high cost, and the transportation mode with low cost will prolong the time.
In order to meet customers' needs, 4PL provider reduces transportation time and transportation cost and optimizes decisions by integrating external resources.In previous studies, scholars mostly studied deterministic problems (Yao et al. 2010;Tao et al. 2016).They all choose the transportation route and the 3PL supplier at the same time, considering the transportation time and transportation cost, and establish a mathematical model, using the exact algorithm (such as Dijkstra algorithm) or heuristic algorithm (such as genetic algorithm, ant system) to solve the model, and provide relevant information for the logistics routing optimization problem.
The above research provides a solid foundation for solving 4PLRP.However, the major impact factors are not considered comprehensively, like uncertainty of transportation time and cost due to seasonality and human factors.On other hand, due to the complexity of 4PLRP, which is a NP hard problem, more efficient optimization algorithms a required.In order to enable 4PL provider to meet customer requirements and reduce transportation time and transportation cost, this paper considers the transportation time and 3PL supplier while considering the uncertainty of time and cost due to seasonality, human factors, etc.A multi-objective chance-constrained programming model is established, and a hybrid beetle swarm optimization (HBSO) embedded Dijkstra algorithm is proposed to solve the problem.Study the impact of the customer's emphasis on transportation time and transportation cost on the objective function and determining the transportation plan.Finally, the Taguchi test is used to determine the algorithm parameters of the hybrid algorithm, and the effectiveness of the algorithm is verified by simulation experiments with boxplots.
The remainder of this paper is organized as follows: The relevant literature is reviewed in Sect.2; Description of the fourth party logistics routing problem and mathematical model is given in Sect.3; A hybrid beetle swarm optimization algorithm (HBSO) embedded Dijkstra algorithm is given in Sect.4; Experimental results are presented and discussed in Sect.5; And the analysis of the algorithm is shown in Sect.6; Conclusions and future research appear in the last Section.
2 Literature review

The 4PLRP
In recent years, 4PLRP has become a popular topic for researchers.Many models have been built for 4PLRP with multi-mode of transportation.Zhang and Ling (2002) established a multi-objective 0-1 planning model based on the technical and economic characteristics of different modes of transportation for 4PLRP.Chen et al. (2003) established a directed graph model with multi dimension weights for its optimization problem of how to select the route, the transportation method and the third-party logistics provider.Chen and Su (2010) propose a decisionmaking method based on preemptive fuzzy goal programming for 4PL with preemptive structure.Xiong and Wang (2011) considered the use of multiple modes of transportation as an effective method for logistics companies to reduce transportation time and transportation costs and established a multi-modal transport optimization model with time windows based on graph structure and then, proposed agent selection and path.Cui et al. (2013) and Huang et al. (2013Huang et al. ( , 2016) ) studied the 4PLRP with uncertain environments, such as the duration time.They describe it as fuzzy numbers or uncertain variables, and then, chance-constrained programming model and fuzzy programming model of uncertainty theory are established.Li and Huang (2015) combine the optimization choices of 3PL suppliers for multiple modes of transportation and path selection and establish a fourth-party logistics routing optimization mathematical model for multiple modes of transportation.Huang and Ren (2015) studied 4PLRP with uncertain delivery time under emergency conditions, and an uncertain programming model is proposed based on the uncertainty theory.Tao et al. (2017) established a mixedinteger programming model for the planning problem of 4PLRP with setup cost and edge cost discount policies which are commonly seen in practice.Lu et al. (2020) consider risk preference of customer, and he introduces ratio utility theory (RUT) to build the mathematical model which can describe the 4PLRP with uncertain factors of customer's risk aversion.
In previous studies, scholars mostly studied deterministic problems.They all choose the transportation route and the 3PL supplier at the same time, considering the transportation time and transportation cost, and establish mathematical models.
The main purpose of this paper is to solve the 4PLRP.There are three key issues: (1) transportation route selection; (2) transportation mode selection; (3) 3PL selection.In Chen's work, the objective function of the model is to get 'min F' and find the optimal routes.The route function a ij (k) = (C, T, P, E) consists of many indexes or factors, such as cost, time, ability of 3PL provider, service quality and reputation.There are some similarities between Chen's work and this paper; however, there are many differences.The similarity is that the two works are both selecting 3PL providers under multiple transportation modes.
The differences between Chen's work and this paper are: (1) The factors, transportation time and cost are overall considered in Chen's work, and a single objective model is built.In this paper, transportation time and cost are considered, respectively, and a multi-objective model is built, and the analyses on the experimental results are also different.(2) In Chen's work, the transportation time and cost are taken as determined variables.However, they are uncertain factors in the real-world operations; in Cui's (2013) work, these uncertain factors are studied by fuzzy theory.This paper also considered the uncertain factors and take another unique perspective-stochastic of the factors.In detail, the transit time and transit cost are considered as stochastic factors, and a stochastic programming model is built.

The efficient algorithm
For the 4PLRP, the efficient algorithm is very important for solving the problem.Some exact algorithm and heuristic algorithm are used in the previous works.The Dijkstra algorithm is used by Zhang and Ling (2002) to solve a multi-objective 0-1 planning model with different modes of transportation.Chen et al. (2003) designed a genetic algorithm with adaptive length-variant coding method to solve a 10 nodes 4PLRP.Chen and Su (2010) propose a modified particle swarm optimization to solve the problem of activity assignment of 4PL with preemptive structure.The proposed method does not require weights and prevents results without feasible solutions caused by improper resource setting.Xiong and Wang (2011) design a twolayer optimization algorithm of collaborative optimization of transportation modes, the effectiveness of the model and the optimization algorithm, which is applied to an agent selection and path planning problem.That is the multimode transport optimization model with time windows.Cui et al. (2013) and Huang et al. (2013Huang et al. ( , 2016) ) introduced memetic algorithm (MA), two-step genetic algorithm, improved node-based genetic algorithm (INGA) and improved distance-based genetic algorithm (IDGA) to solve the 4PLRP mode with uncertain variables.Huang et al. (2013) designed a two-step genetic algorithm with the fuzzy simulation to find approximate optimal solutions.Numerical experiments are carried out to investigate the performance of the proposed algorithm.Li and Huang (2015) design a maximum and minimum ant system for a fourth-party logistics routing optimization problem.Tao et al. (2016) proposed a column generation approach combined with graph search heuristic to efficiently solve the 4PLRP.The good performance in terms of the solution quality and computational efficiency of the approach is shown through extensive numerical experiments on various scales of test instances.Fuqiang et al. (2020) design an ant colony system-improved grey wolf optimization (ACS-IGWO) for the 4PLRP with risk preference.The ability of ACS-IGWO is tested in simulation experiments.
In the above works on the solution method of 4PLRP, many researchers use the exact algorithm or heuristic algorithm to solve the model of deterministic problems.But for an uncertain problem of 4PLRP, the exact algorithm like Dijkstra is not enough to finish the work.On the other hand, the traditional heuristic algorithms like genetic algorithm and ant system are very difficult for the further improvement.Therefore, some newly developing algorithms have great potential to become a good solution method.
The HBSO proposed in this paper has many differences with BAS-PSO introduced by Lin and Li (2018).In Lin and Li's work, the PSO is combined with BAS, each beetle is taken as a particle with speed and position, the PSO is used for global search, and the BAS is used for local search.This framework is applied in our HBSO; however, there are many differences between them.Firstly, the Dijkstra is combined in the stage of the population initialization for HBSO, which is also designed according to the characteristics of the problem in this paper.Secondly, the BAS-PSO is only tested by four test functions, but has not been used in any engineering optimization problem.This paper designed HBSO for routing optimization problem in logistics, which has more complexity than the test functions, usually.The work of this paper extends the application field of the BAS-PSO (or HBSO) and provides new solving method for this kind of routing optimization problem.
In summary, researchers have done a lot of work to solve 4PLRP.Some of them considered the uncertainty of transportation time, and others considered the uncertainty of cost.In detail, the uncertainty of cost can be caused by many reasons, like seasonality and human factors.But a comprehensive consideration is not existing, that is what we should do in the next.On other hand, 4PLRP is a complex NP hard problem, and researchers have tried many methods to deal with it, such as GA, PSO, ACS, BAS, GWO and their variants.But more efficient optimization algorithms are still needed.
3 Problem description and modeling

Problem description
A 4PL company undertakes a transportation task, and the transportation network is represented by a multi-graph G(V, T, E), as shown in Figs. 1 and 2 (Cui and Huang 2013).V denotes the set of transportation note, T denotes the set of transportation mode, and E denotes the set of transportation path.There are multiple 3PL suppliers and their transportation modes between the two nodes are shown in Fig. 2, where ''1, 2'' indicates the 3PL supplier number that can be selected between the two nodes and ''(1), ( 2), (3)'' indicates the three transportation modes (highway, railway, and aviation) owned by 3PL suppliers.4PL company needs to choose the appropriate transportation path, 3PL supplier and its transportation mode to transport the goods from the source point V s to the destination point V t .In the process of transportation, V 2 * V n-1 is used as a transit node for converting 3PL supplier and its transportation modes.And each 3PL supplier has multiple transportation modes (Li et al. 2014).Each transportation mode is regarded as an edge between two adjacent nodes in the transportation network.When the transit of transportation mode occurs, a certain time and cost are also required.Based on the above factors, 4PL company should select the optimal transportation route to minimize the transportation cost and time, the optimal transportation route is shown by the thick line in Fig. 1.

Mathematical model
A. Assumption Each 3PL supplier has ability to finish the transportation task, and only one 3PL supplier is selected to undertake the transportation task between two nodes.Only one transportation mode of 3PL supplier can be selected between any two adjacent nodes, and the conversion of transportation mode can only occur at the transit nodes.The transportation cost and transportation time at the transit node are uniformly distributed.

B. Parameters
Model parameters are defined in Table 1.

C. Mathematical model
Due to the traffic environment and human factors, the uncertainty of the transit time and transit cost during the process of transferring from one transportation mode to another transportation mode.Therefore, the transportation time and transportation cost for a route are calculated by Eqs. ( 1) and (2).
This paper aims to find a route from the source to the destination with minimum time and minimum cost, and the solution is also subject to the constraints.A multi-objective chance constrained programming model is established as follows: Equations ( 1) and ( 2) are the objective functions, f 1 represents the minimum value of TðRÞ under confidence level b, b 2 ð0; 1Þ, and f 2 represents the minimum value of CðRÞ under confidence level a, a 2 ð0; 1Þ; Eq. ( 5) is a chance constraint for the transportation cost under Fig. 2 Multi-graph of multiple transportation mode between two adjacent nodes confidence level a; Eq. ( 6) is a chance constrain for transportation time under confidence level b; Eq. ( 7) is the constraint of capacity required by the customer; Eq. ( 8) is the constraint of reputation required by the customer; Eq. ( 9) is used to ensure that the selected nodes and edges is a legal route, which includes transportation path, 3PL supplier and its transportation mode.Equations ( 10), ( 11) and ( 12) represent that x l ijk ðRÞ; y i ðRÞ and z l;u i ðRÞ are 0-1 decision variables, respectively.

Algorithm designed
The beetle antennae search (BAS) algorithm (Jiang and Li 2017;Wang and Li 2022) has been proposed in recent years and belongs to the emerging intelligent optimization algorithm.It has been successfully applied in many fields (Wang et al. 2018;Zhao and Qian 2018;Chen et al. 2018;Sun et al. 2019;Al-qaness et al. 2021; Chen and Zhan 2022; Rajeshwari and Sughasiny 2022) for its high efficiency.However, the performance of the BAS algorithm is not competent to the high-dimensional function problems.This is because the BAS algorithm is used by a beetle to find the optimal solution through continuous iteration, and the initial position of the beetle has a great impact on the results (Wang et al. 2018;Chen et al. 2019;Xu et al. 2019;Abed-Alguni and Alawad 2021;Abed-alguni and Paul 2022).Therefore, the idea of population search from particle swarm optimization (Shi andEberhart 1998, 1999;Yousri et al. 2022;Ewees et al. 2022;Al-qaness et al. 2022) (PSO) is used to improve the BAS algorithm.Multiple beetles are used to replace the single beetle to search the solution of the problem; therefore, a Hybrid Beetle Swarm Optimization (HBSO) algorithm is designed to combine the global search capability of PSO with the local search capability of BAS.It can fully develop the global search ability of the PSO, and at the same time, it can also play the local search ability of the BAS.In this way, the optimization function can be complementary, and the optimization ability of the algorithm can be improved to a high level.
In order to better solve the resulting mathematical model, Dijkstra's algorithm (Yue and Gong 1999) is used to initialize the problem.In HBSO, the solution process is divided into two phases, includes generating an initial population by Dijkstra algorithm and updating the population by HBSO.In the first stage, a simple graph is randomly generated according to the multi-graph of the 4PLRP, and the shortest path is obtained according to the Dijkstra algorithm, and the fitness value is calculated.In the second stage, the population obtained from the previous stage is taken as the initial population of HBSO, and then, the 3PL supplier and its transportation mode are updated; the optimal solution is finally obtained.

Coding and population initialization
According to the characteristics of the 4PLRP, an encoding based on two arrays (Liu et al. 2014) is adopted to represent a beetle in HBSO, where the first array consists of nodes on a route encoded by natural number, and the second array consists of edges (3PL suppliers) and its transportation mode on a route encoded by integer number.Because the quality of the coding directly affects the computational efficiency and the accuracy of the solution.In order to provide a good solution for the HBSO algorithm, the Dijkstra algorithm is used to initialize the beetle population.However, the condition for initialization using the Dijkstra algorithm is that the multi-graph must be a simple graph, that is, only a transportation mode with one 3PL supplier between two adjacent nodes.Therefore, the The number of 3PL supplier (edges) between node i and node j (i, j = 1, 2, …, N) Then, calculate the weight of path according to Eq. ( 13): In Eq. ( 13), W ij represents the weight on the edge between two adjacent nodes, and the calculation of the weight includes the transportation time and transportation cost on the path.

Updating of the solutions
The speed and position updating formula of BSO algorithm is used to determine the 3PL supplier and its transportation mode.The specific operations are as follows: There is a population with NP beetles, and the position of beetle gðg ¼ 1; 2; :::; NPÞ is represented as X gij ¼ ðX gij1ð1Þ ; X gij2ð2Þ ; :::; X gijrðmÞ Þ, and the speed of beetle g is represented as V gij ¼ ðV gij1ð1Þ ; V gij2ð2Þ ; :::;V gijrðmÞ Þ.For 8k 2 f1; 2; :::; rg; l 2 f1; 2; :::; mg, The position of the beetle g is updated by Eq. ( 14): where V b gijkðlÞ is expressed as the speed of beetles, and q b gijkðlÞ represents the increase in beetle position movement; b ¼ 1; 2; :::; NG, NG is the number of iterations;c is a position constants.
Then, the speed formula of beetle g is given: where c 1 and c 2 are positive constants; r 1 and r 2 are random numbers in the range ½0; 1, respectively; P b gijkðlÞ represents the best position of beetle g; P 0 b gijkðlÞ represents best position of the population so far.
The position incremental function of beetle g is calculated as follows: where d b ¼ eta Á d bÀ1 represents the step factor of beetles; f ðXL b gijkðlÞ Þ and f ðXR b gijkðlÞ Þ represent the fitness of scent intensity at the right and left antennae for beetles, respectively; signð Þ is a sign function.
Then, create beetle g's left and right spatial coordinates by Eqs. ( 17) and ( 18), respectively: where d b ¼ d b c 2 represents the distance between two antennas

Fitness function
The objective function is used as the fitness function.The multi-objective model is transformed into a single-objective model by linear weighting method, and the weight coefficients k 1 and k 2 are added to the two objectives, respectively.The objective function is expressed as follows: where The cost and time are adjusted by the weight coefficients k 1 and k 2 , which means that 4PL pays different attention to the cost and time according to customer's degree and provides decision for different needs.

Computational complexity
The computational complexity of HBSO mainly depends on three processes: population initialization, Dijkstra Algorithm and updating of the beetles.Supposing that there are NP beetles, and the computational complexity of population initialization is O(NP).The complexity of the updating mechanism is O(NG 9 NP) ?O(NG 9 NP 9 D), where NG is the maximum number of iterations, and D is the dimension of a specific solution.The complexity of Dijkstra Algorithm is O(|V|^2), where V is the number of vertices in the weighted directed graph.Therefore, the computational complexity of HBSO is O(NP 9 (1 ?NG ?NG 9 D) ?|V|^2).

Procedures of HBSO
The procedures of HBSO are shown as follows: Step 1: Set the parameter value of HBSO, Including NP, NG, c 1 , c 2, d; Step 2: The Dijkstra algorithm is used to obtain the position information of the beetle group, including the node set, 3PL supplier and its transportation mode set; Step 3: According to Eqs. ( 19) and ( 20), the fitness is calculated corresponding to beetle, and the initial individual optimal solution and global optimal solution are obtained; Step 4: Calculate the left and right spatial coordinates of each beetle according to Eqs. ( 17) and ( 18), and obtain the corresponding fitness; Step 5: Update the speed and position increment of the beetle according to Eqs. ( 15) and ( 16), respectively; Step 6: Update the individual optimal solution and the global optimal solution, respectively.
Step 7: If the maximum number of iterations reached, go to step 9, otherwise, go to step 3; Step 8: Stop and output the optimal solution.The flowchart of HBSO is shown in Fig. 3.

Numerical analysis
In this section, numerical analysis is provided to test the effectiveness of the proposed mathematical model and algorithm.First, three examples are given with different problem scale.Then, based on customer's different attention on time and cost, the relationship between the weight of the objective function and the confidence level is analyzed.For the different parameter combinations in the mathematical model, the suitable transportation plan for the customer is obtained.Finally, considering the influence of algorithm parameters on the experimental results, the orthogonal experiment is designed to obtain the optimal parameters combination of the algorithm, and the performance of BAS, GA, PSO and HBSO is compared and analyzed.

Numerical examples
In this paper, three examples are used for numerical analysis, which are example I with 8 nodes, example II with 16 nodes and example III with 32 nodes, respectively.Experimental data are referenced in literature (Lu et al. 2020), and some revisions are made.The settings for the random transportation time and transportation cost of the transit node are shown in Tables 2 and 3, respectively, where the random transit time and transit cost subject to the uniform distribution.Take symbol U(1, 6.5) in Table 2 for example, 3PL 1's transit time is a random number, which is subject to uniform distribution U(1, 6.5).
The numerical experiments are performed by a 3.0 GHz quad-core processor, 8 GB memory and a Win 10 system.The numerical experiment is simulated by MATLAB R2014b software.

Model parameter analysis
The important model parameters are the weight k 1 , k 2 in the objective function and the confidence level a, b in the constraint, respectively.For the weight coefficient, Considering different attention of customers on cost and time, there are 11 combinations of values for k 1 and k 2 .Take example II as an example, the model parameters analysis is divided into two angles, and the confidence level a and b are equal and unequal, respectively.The experimental results are shown in Tables 4 and  5, Figs. 4 and 5, respectively.
1.The confidence level a and b are equal Table 4 and Fig. 4 show the relationship between the objective function and the model parameters when a = b, k 1 þ k 2 ¼ 1.For the same confidence level, the objective function value shows an increasing trend with the increase in the weight k 1 and k 2 .In practice, customers pay more attention to transportation costs than time.Under the same weight combination, as the confidence level a and b increase, the objective function value is also accompanied by an increase.In real life, when 4PL meets customer's requirements under a certain probability, as customer's requirements increase, the transportation time and transportation cost will also increase.
2. The confidence level a and b are unequal Table 5 and Fig. 5 show the values of the objective function and the relationship between them when the model parameters aþb¼ 1 and  2 U(4.5, 10) U(4, 9) U(3, 8) U(1.5, 7.5) U(6.5, 12) U(6, 11.5) Aviation 1 U(7, 12.5) U(8, 13) U(6, 11) U(6.5, 12) U(4, 9) U(5, 10.5) 2 U(8, 13) U(7, 12) U(6.5, 12) U(6, 11.5) U(4, 10.5) U(4, 9.5) can be seen from the graph, under the same confidence level combination, as the value of cost weight k 1 increases and the value of time weight k 2 decreases; the objective function exhibits an increasing trend.Under the combination of confidence levels where a varies from 0 to 1, b varies from 1 to 0, respectively, while k 1 2 ð0; 0:4Þ; k 2 2 ð0:6; 1Þ, as the confidence level a increases and the confidence level b decreases, the objective function shows a decreasing trend; when k 1 2 ð0:6; 1Þ; k 2 2 ð0; 0:4Þ, as the confidence level a increases and the confidence level b decreases, the objective function shows an increasing trend; When k 1 ¼ k 2 ¼ 0:5, no matter how a and b change, due to the randomness of the problem, the objective function maintains a small fluctuation above and below a value.
The above results show that when the confidence level a þ b ¼ 1 and the weighting factor k 1 ¼ k 2 ¼ 0:5, the change of objective function is independent of in the confidence level.When the cost weight k 1 is low and the time weight k 2 is high, the objective function decreases as the confidence level a increases and the confidence level b decreases.Therefore, the objective function is greatly affected by the confidence level b; the cost weight k 1 is higher.When the time weight k 2 is low, the objective function increases as the confidence level a increases and the confidence level b decreases, so the objective function is greatly affected by the confidence level a.

Results and comparative analysis
In this section, take the confidence level a¼b¼ 0.8 as an example, Table 6 shows the experimental results for example I, example II and example III.
As can be seen from Table 6, the customer's preference on transportation time and transportation costs will not  affect the choice of transportation route, but will change the choice of 3PL supplier and its transportation mode.For example, the optimal route for example I, II and III is [1 4 6 8], [1 2 8 14 16] and [1 9 13 18 26 30 32] under different weight parameters k 1 and k 2 , respectively.But the 3PL supplier and its transportation mode are different, and there are three different choices that are [1(2) 1(2) 2(2)], [1(1) 2(2) 2(2)] and [1(2) 1(2) 2(1)] in example I.This is because the 4PL provider can change the weighting factor corresponding to the transportation time and transportation cost when the customer pays different attention to the two factors.So as to make the corresponding decision for providing the most economical transportation route, 3PL supplier and its transportation mode.
6 Algorithm analysis

Parameters setting
The Taguchi method (Fuqiang et al. 2019;Raz et al. 2021) is used to generate the optimal parameters setting of the algorithm.The main purpose of using the Taguchi method is to find the optimal parameters setting by small scale experiments in the most economical and effective way.In this experiment, Minitab software is used to carry out the Taguchi method.The optimal parameters setting is shown in Table 7.
Example II is taken to explain how to use this method to determine the parameters setting.The different parameters setting are tested including NP = 30,40,50;NG = 30,40,50; C 1 = 1.5, 2.0, 2.5; C 2 = 1.5, 2.0, 2.5; c = 0.3, 0.5, 0.7 and d = 1.2, 1.5, 1.8.The parameters are set on three levels based on the results of testing.The optimal parameters setting is generated from the following three steps.
(3) Determination of optimal parameters setting According to the literature (Xiong and Wang 2011), the optimal parameters setting for example II are NP = 50,NG ¼ 50,C 1 ¼ 2,C 2 ¼ 2,c ¼ 0:5 and d ¼ 1:2.In this way, the parameters setting for example I and example III can be obtained, as shown in Table 7.

Performance comparison of algorithms
In order to show the ability of HBSO, BAS, GA and PSO are introduced to solve the resulting optimization problem.The results are shown in Table 8, where the confidence level is a ¼ b ¼ 0:8 and the weight is k 1 ¼ k 2 ¼ 0:5.In order to show the comparison of the algorithms more clearly, the population size and the number of iterations are set as NP = 80, NG = 100.The algorithms were run 100 times, and the best fitness (Best), the worst fitness (Worst), the average fitness (Avg), the variance of fitness (Var) and the average running time (Time) of the algorithm were recorded, See Table 8.
It can be seen from the table, for the running time, since the BAS algorithm uses a beetle to search, BAS is the fastest one among the four algorithms; the HBSO algorithm is faster than the GA and PSO algorithms, and the speed is sorted as BAS [ HBSO [ GA [ PSO.And the best fitness calculated by GA, PSO and HBSO algorithm has little difference.Although the results of the four indicators Best, Worst, Avg and Var are similar for small scale problem like example I.But with the increase in problem scale, for example II and III, the Best, the Avg of HBSO algorithm are superior to the BAS, GA and PSO algorithms, and the four algorithms are comprehensively ranked as HBSO [ PSO [ GA [ BAS.
Boxplots is used to compare the optimization ability of algorithms, due to the poor quality of the BAS's solution, we only consider three algorithms, HBSO, PSO and GA.Based on the experimental results, Figs. 8, 9 and 10 are drawn.Firstly, as can be seen from the figures, for terms of optimization ability, the HBSO finds the best solution among the three algorithms in every example.Secondly, the size of the box reflects the dispersion degree of the samples, and the dispersion degree of HBSO is less than GA and PSO; in other words, the sample data of HBSO is more concentrated on the average.In addition, compared with boxplots of GA and PSO, the boxplot's median and quartile of HBSO are more concentrated, which proves the HBSO's reliability.The reasons for this result are that the HBSO algorithm turns a beetle of the BAS algorithm into multiple beetles, which increases the population size, makes the algorithm easier to find the optimal solution and also performs better by embedding the Dijkstra algorithm to the stage of population initialization.
The best iteration of the four algorithms for example I, II and III is shown in Figs. 11,12,13.It can be seen that HBSO has a relatively fast convergence speed, and it can achieve a solution slightly better than GA.This article also has some shortcomings.For example, this paper only considers the logistics distribution task from single origin to single destination.In future research, a multi-origin to multi-destination distribution network will be considered.The Floyd-Warshall algorithm will be used in future work to calculate the minimum distance between two points.
The lth l (l = 1, 2, …, M) transportation mode of the kth (k = 1, 2, …, K) 3PL supplier between node i and node jC l ijkThe unit cost for the lth transportation mode for the kth 3PL supplier between node i and node jC iThe waiting cost at node i T l ijk The unit time o for the lth transportation mode of the kth 3PL supplier between node i and node j T i The waiting time at node i P l ijk The transportation capacity for the lth transportation mode of the kth 3PL supplier between node i and node j D l ijk The historical reputation for the lth transportation mode of the kth 3PL supplier between node i and node j n l;u i Random transit cost from transportation mode l to transportation mode u(u = 1, 2, …, m) at node i g l;u i Random transit time from transportation mode l to transportation mode u at node i 4PL routing problem using hybrid beetle swarm optimization multi-graph of the 4PLRP is processed into a simple graph.

Fig. 4 Fig. 5
Fig. 4 Influence of model parameters under equal confidence level a and b

Fig. 12
Fig. 12 Best iteration of algorithms for example II

Table 1
Definition of parameters

Table 2
Transit time on transfer node for different mode of 3PL suppliers

Table 7
Parameters setting for the three examples