Multi-task evolutionary optimization of multi-echelon location routing problems via a hierarchical fuzzy graph

Multi-echelon location-routing problems (ME-LRPs) deal with determining the location of facilities and the routes of vehicles on multi-echelon routing tasks. Since the assignment relationship in multi-echelon routing tasks is uncertain and varying, ME-LRPs are very challenging to solve, especially when the number of the echelons increases. In this study, the ME-LRP is formulated as a hierarchical fuzzy graph, in which high-order fuzzy sets are constructed to represent the uncertain assignment relationship as different routing tasks and cross-task operators are used for routing task selection. Then, an evolutionary multi-tasking optimization algorithm is designed to simultaneously solve the multiple routing tasks. To alleviate negative transfer between the different routing tasks, multi-echelon assignment information is considered together with associated routing task selection in multi-tasking evolution optimization. The experimental results on multi-echelon routing benchmark problems demonstrate the competitiveness of the proposed method.


Introduction
Location routing problems (LRPs) [28] originate from transportation distribution systems in city logistics, and many real-world applications can be formulated as LRPs, such as multimodal transportation [12], space science [1], food supply chains [9], and express delivery service [2].An LRP is a complex combinatorial optimization problem with different kinds of resource constraints, aiming to determine a series of routes with a minimal routing cost by handling facility location and vehicle routing at the same time.To solve LRPs, several evolutionary algorithms have been proposed [22,26,27], where facility location and vehicle routing are It has been suggested that the assignment relationships for each echelon should be considered in local search routing strategies [5,23].For instance, the impact of satellite-tocustomer assignment is investigated in two-echelon location routing problems, and the sub-problems on the first and second echelons have different effects on generalized cost [7].After that, a fuzzy correlation based on a two-echelon assignment relationship is designed, and an adaptive neighborhood local search scheme is employed for two-echelon routing planning [21].In addition, a graph-based evolutionary algorithm [36] is presented to deal with the possibilities of the satellite-to-customer assignment on the two-echelon routing problem, which takes advantage of the population-based evolution and refrains from excessive fitness evaluations of unpromising moves in different intermediate facilities.
The graph-based evolutionary search has been demonstrated highly competitive for the routing problems of each echelon in terms of solution quality and search efficiency [33].However, it is non-trivial for these algorithms to model the assignment relationship for each echelon routing subproblem when the number of the echelons in an ME-LRP increases.
Recently, the idea of jointly optimizing a series of routing sub-problems has been proposed [29,30].In addition, the evolutionary multitasking framework has been demonstrated powerful for solving a series of routing problems [10], where specific characteristics of routing problems are incorporated and a range of routing sub-problems are considered as different tasks.Note, however, a simple clustering strategy is used for dealing with assignment relationships in [10].Consequently, the performance of the evolutionary multitasking approach may degrade when dealing with two or more echelon location routing problems.Moreover, the representation of the individual only with distance traits is not conducive to joint optimization of multiple routing subproblems at different echelons [16].Out of the above reasons, the evolutionary multitasking framework will become more beneficial for solving multi-echelon routing sub-problems, if the multi-echelon assignment relationships can be modeled with the idea of joint optimization.
In this work, we propose a multi-task evolutionary optimization approach to solving ME-LRPs with the help of a hierarchical fuzzy graph, which is called MEO-HFG.To model the multi-echelon assignment relationships in ME-LRPs, we adopt a hierarchical fuzzy graph to represent an ME-LRP as a series of sub-problems by considering the similarities and differences between multi-echelon routing sub-problems.Then, the multi-tasking evolutionary optimization framework is used to solve the sub-problems simultaneously, based on the estimation of multi-echelon assignment information with a hierarchical fuzzy graph.
The main contributions of this paper are summarized as follows.First, an ME-LRP is modeled as a series of joint routing sub-problems by employing a hierarchical fuzzy graph, which enables us to estimate the uncertain assignment relationships for multi-echelon routing tasks.Second, a multi-tasking evolutionary algorithm with the hierarchical fuzzy graph representation is proposed, where the associated task selection strategy is combined with multi-echelon assignment information to avoid negative knowledge transfer in different routing tasks.Third, adaptive operators based on the hierarchical fuzzy graph are designed for different search stages to enhance search performance.Finally, we performed extensive experiments to verify the performance of the proposed method, and the results demonstrate that it significantly outperforms the state-of-the-art methods on multi-echelons location routing problems with different routing tasks.
The remainder of this paper is organized as follows."Problem description and background" begins with an introduction to the mathematical model of ME-LRPs.A brief literature review of recent progress in the evolutionary multitask optimization for ME-LRPs is also given."Multi-tasking evolutionary optimization via a hierarchical fuzzy graph" presents the details of the proposed MEO-HFG for ME-LRPs."Experimental settings" describes the experimental results on the benchmark problems.Finally, concluding remarks of this article are drawn in "Conclusion and future work".

Problem description and background
To facilitate the understanding of our proposed method, the definition of ME-LRPs is introduced first in this section.Next, a brief review of multifactorial evolutionary algorithms (MFEAs) for multitasking optimization is presented.

Problem description
An ME-LRP can be formally defined in a weighted and directed graph G = (V , E), where V and E denote the set of vertices and edges, respectively.Set V consists of the subset of the suppliers' nodes V d , the subset of the customers' nodes V c , and the subset of multi-echelon intermediary facilities (satellites) nodes V s = ∪ N −1 e=1 V s e , such as regional centers, distribution centers, and plants.N is the number of routing echelons for the ME-LRP.Set E = N i=1 E i is divided into N parts, representing the edge sets at different echelons.E 1 = {(i, j) : i, j ∈ V s 1 V c } includes the edges connecting the first echelon satellites to the customers, as well as those connecting pairs of customers.
comprises the edges connecting the (k − 1)-satellites to the k-satellites, as well as those connecting pairs of the k − 1 echelon satellites.Multi-echelon refers to a level of multiple facilities of distribution systems, where the freight delivery from the supplies to the customers is often carried out by incorporating one or multiple intermediate facilities (satellites), so the solution of different ME-LRPs can be represented by a series of routes on multi-echelons.Three illustrative examples of feasible solutions to two-echelon, three-echelon, and four-echelon LRPs are given in Fig. 1, where the regional centers, distribution centers, and plants are called 1-satellites, 2-satellites, and 3-satellites set separately.The aim of solving ME-LRPs is to determine a set of vehicle routes together by consolidating the freight through multiple satellites, including location decisions, assignment decisions and routing decisions.
The mathematical model of an ME-LRP is given as follows: subject to where the objective function of the ME-LRP in Eq. ( 1) is to minimize the overall multi-echelon cost, including the fixed cost of opening the facilities (supplies and a set of satellites), the fixed cost for using vehicles on multi-echelon routing, and the transportation cost for multi-echelon routes.The number of echelons for routing sub-problems in the ME-LRP is N .For the k-th routing sub-problem (in the bottom-up order), the set of upper-level nodes is defined as U k , and the set of lower-level nodes is defined as L k .The variable y u (y u ∈ 0, 1) is equal to 1 only if the supplies or the satellites u, u ∈ U k is opened in the solution.F u represents the fixed cost of opening the supply or the satellite u. f k , vc k , and vt k stand for the maximum capacity, the maximum driving distance, and the fixed cost of using a vehicle on the k-th echelon routing sub-problem, respectively.dem_D k (i) represents the actual demand of the i-th lower level node on the k-th echelon routing sub-problem, while cap_D k (i) represents the capacity of the i-th upper level node on the k-th echelon routing sub-problem.q k j denotes the distance cost of each route r j , r j ∈ R k , while t k is the per-unit transportation cost on the k-echelon routing sub-problem.The variable z j (z j ∈ {0, 1}) is equal to 1 only if the route r j is in the solution.The variable Z i j (k)(Z i j (k) ∈ {0, 1}) is equal to 1 only if there is an assignment relationship between the upper level node i and the lower level node j on the k-th echelon routing sub-problem.The variable S i j (k)(S i j (k) ∈ {0, 1}) is equal to 1 only if node i and node j are sequentially adjacent nodes when they are served by the same vehicle on the k-th echelon routing sub-problem.
The vehicle capacity constraint is given in Eq. ( 2), where the total delivery demand of each route should not exceed the maximum capacity of each vehicle for multi-echelon routing sub-problems.The freight balanced constraints are given in Eqs.(3)- (5).The total freight capacity of the upper nodes set U k needs to meet the total delivery demand of the lower nodes set U k for the k-echelon routing sub-problem, which are indicated in Eq. (3).The upper nodes set U k for the k-echelon routing sub-problem is also regarded as the lower nodes set U k+1 for the k + 1-echelon routing sub-problem.Therefore, as shown in (4), the total freight demand of the lower nodes set U k+1 , which should be equal to the total freight capacity of the upper nodes set U k , cannot exceed the total capacity of the upper nodes for k + 1-echelon routing sub-problem.Any change in the assignment relationship on the current echelon routing sub-problem may more or less affect the assignment relationship on the adjacent upper echelon as well as routing sub-problems at higher echelons.In a sense, the diversity of the assignment relationship in multi-echelon nodes is determined by the freight-balanced constraint.The sum of the capacities of the depots in the uppermost routing task should meet the sum of the demands of all customers in the lowermost routing task, which is described in Eq. (5).Equations ( 6) and (7) state that the customer visited constraints, where each customer, which is served by exactly one vehicle, is assigned exactly once by the intermediate facilitates.If a lower-level node has an actual demand in the k(1 < k ≤ N ) echelon routing task, then the corresponding upper-level node must be assigned to it in Eq. ( 8), and a lower-level node can appear at least one route in Eq. ( 9).If a lower-level node does not have an actual demand in the k(1 < k ≤ N ) echelon routing task, any upper-level node cannot be assigned Eq. (10) and the lower-level node cannot appear in any routes in Eq. (12).The flow conservation constraints are defined in Eqs. ( 8)- (12), meaning that the flow from the upper nodes to the lower nodes on the current echelon routing sub-problem depends on whether there is an assignment relationship between the upper and lower nodes on the adjacent lower echelon routing sub-problem.Equation (13) denotes that the distance of each path cannot exceed the maximum distance of each vehicle.
In this work, we consider two-echelon, three-echelon, and four-echelon LRPs.The possibilities of multi-echelon assignment relationships result in the high complexity of solving a number of routing sub-problems on multi-echelon in the presence of constraints.

MFEAs
MFEAs are one class of the most widely used evolutionary multi-tasking optimization algorithms.Recently, some variants of MFEAs [14,19] have been proposed for problem-specific optimization, which can deal with several tasks simultaneously in evolutionary computing by using cross-task similarities and complementary.To facilitate the interaction between multiple tasks in the evolutionary optimization, a deliberately designed objective function needs to be designed to consider all the tasks.Assume that there are N tasks to be solved together, and all tasks T j , j = 1, 2, . . ., N are minimization problems.The objective function for each task is f j : X j −→ R, X j is the search space.The aim of MFEAs is to tackle multiple optimization tasks simultaneously using only one single population, and find a set of solutions x * 1 , x * 2 , . . ., x * N such that Some basic features defined in MFEAs for each individual are explained as follows: 1. Factorial Cost: the factorial cost of individual p i on task T j is denoted as Ψ i j = f i j , which means the factorial cost is the corresponding objective function value of individual p i on task T j when p i is a feasible solution.2. Factorial Rank: The factorial rank r j i of individual p i on task T j can be determined by ascending the factorial cost of the current population on task T j .In particular, if p i is the best individual, it has r j i = 1. 3. Skill Factor: The skill factor τ of individual p i is the index of the task on which the individual p i can perform most effectively.Given the performance ranks of individuals on all N tasks, the skill factor of individual p i can be calculated as τ i = arg min j {r i j }, j ∈ 1, 2, . . ., N .
4. Scalar Fitness: The scalar fitness demonstrates the ability for individual p i to solve a variety of tasks, which can be expressed as φ i = 1/min j r j i , j ∈ 1, 2, . . ., N .
The framework of an MFEA can be divided into four core parts, namely initialization, assortative mating, selective evaluation, and population update.Different from the canonical evolutionary algorithms, a parameter called random mating probability, rmp, is employed to enhance the exploration ability of cross-task evolution for assortative mating.Besides, only a single task is chosen for offspring evaluation for each individual.A multi-tasking evolutionary algorithm conducts an evolutionary search in multiple search spaces (of multiple optimization sub-problems) by leveraging transferable knowledge across different tasks.On the basis of MFEA, some multi-tasking evolutionary algorithms have exhibited better performance on simple routing optimization problems [10] in terms of the search efficiency with the help of knowledge transfer between the individuals for different tasks [15].Recently, some evolutionary many-tasking algorithms [11,34] are proposed to enhance optimization performance via the knowledge transfer between a large number of tasks along the evolutionary search process.However, for more complex routing problems, the multi-tasking and many-tasking evolutionary algorithms need to explore more effective transfer mechanisms to alleviate the negative influences between jointly learned routing tasks.[35,38].In addition, the scalability to the task number [32] needs to be considered by gradually increasing the number of routing tasks performed.

Multi-tasking evolutionary optimization via a hierarchical fuzzy graph
In this section, ME-LRPs are transformed into a set of routing sub-problems with a hierarchical fuzzy graph formulation and solved by a multi-tasking evolutionary algorithm.Considering the flow conservation constraint of ME-LRPs, the possibilities of multi-echelon assignment relationships should be considered in parallel [36].In addition, there exists capacity dependence between each pair of routing tasks on adjacent upper and lower echelons for the freight balanced constraint.In this work, we treat each echelon routing sub-problem as a task and present a hierarchical fuzzy graph formulation of multi-echelon routing sub-problems with a series of high-order fuzzy sets.Next, we solve the multi-echelon routing sub-problems simultaneously using multi-tasking evolutionary optimization, resulting in a multitask evolutionary optimization approach call MEO-HFG.
The overall framework of the proposed MEO-HFG is illustrated in Fig. 2, and the pseudo-code of the proposed MEO-HFG is given in Algorithm 1.The main components of the proposed algorithm include the hierarchical fuzzy graph Obtain the cross-task probability matrix rmp in different routing sub-task feasible regions by Eq. ( 22).8: Generate an offspring population P new across the hierarchical fuzzy graph G h f in Algorithm 2. 9: Evaluate each individuals in offspring population P new for selected optimization routing tasks.10: Update the scalar fitness φ and the skill factor τ of every individual in intermediate population P t .12: Sort the solution with the fittest in P t , and then select N pop better individuals to form current population P cur .13: Sort the solution in P cur by Eq. ( 1), and update the best solution Q r with the individual in P cur with the best global fittest.

14:
Update current the hierarchical fuzzy graph for the next generation with fuzzy statistical method.15: end while formulation, solutions representation and initialization, associated task selection, cross-task operators, evaluation and selection, and the update of the hierarchical fuzzy graph representation.

Hierarchical fuzzy graph formulation of ME-LRPs
To construct the correlation in different routing subproblems, this work formulates the ME-LRP with a hierarchical fuzzy graph representation.Note that the multi-echelon uncertain assignment information plays an important role in solving ME-LRPs.
On the one hand, the variety of multi-echelon assignment relationships increases the number of routes generated in different echelons.If a certain echelon assignment relationship changes, the routes in the adjacent upper routing task and lower routing task may also be changed as well as the routes in the current routing task.On the other hand, for each routing sub-problem (task), there exist upper nodes and lower nodes.Only if there is an assignment relationship between the upper node in the current routing task and those in any lower nodes, meaning that the upper node is available, we need to arrange the reasonable routes under multiple constraints of an ME-LRP.Therefore, the assignment relationship may have a direct or indirect influence on any two routing tasks.Thus, an ME-LRP can be reduced to a multi-echelon uncertain assignment relationship problem if  we ignore the influence of multi-echelon route generated for each routing sub-problems.
To model the multi-echelon uncertain assignment relationships, we use a hierarchical fuzzy graph [25] to formulate the different routing tasks in ME-LRPs.As visualized in Fig. 3, an example of a four-echelon LRP is formulated as a hierarchical fuzzy graph, which consists of one 4-order fuzzy set, two 3-order fuzzy sets, three 2-order fuzzy sets, and four 1-order fuzzy sets.Associated high-order fuzzy sets are employed to express the uncertain assignment in different routing tasks.The weights in the hierarchical fuzzy graph represent the values of the assignment relationship memberships.
The high-order fuzzy sets can be built on a family of generic terms that are not only fuzzy sets but also a family of low-order fuzzy sets.The k-order fuzzy set for multi-echelon uncertain assignment can be defined as follows.
where Γ stands for the k-order fuzzy sets, and n 0 is the number of lower nodes for the routing task_k.Γ l p i 1 i 2 ...i p stands for the (k − p) th order fuzzy set and n i p (1 ≤ p ≤ k − 1) is the number of lower nodes for the routing task_(k− p).μ represents the membership values of the assignment of the

Solutions representation and initialization
During the multi-task evolutionary optimization, the fuzzy assignment of the multi-echelon nodes and a series of routing node orders served in the solution can be illustrated across the hierarchical fuzzy graph representation.With the assignment information of the hierarchical fuzzy graph G HF for each generation, we also define a cross-list Lf = (Lf 1 , Lf 2 , . . ., Lf i , . . ., Lf m ) consisting of a series of routing nodes orders of each routing tasks for an individual.Lf i (1 ≤ i ≤ m) means the sequences of the lower nodes of the routes through the upper node, and the association in different lists Lf i can be represented in the hierarchical fuzzy graph.An illustrative hierarchical fuzzy graph representation for the four-echelon LRP solution is depicted in Fig. 4, where Fig. 4a depicts the hierarchical fuzzy graph information for the individual in the current generation, and Fig. 4b represents cross-list storage of the individual for the four-echelon LRP in Fig. 1 based on the hierarchical fuzzy graph information.The assignment correlation between the upper nodes and lower nodes for each routing task can be obtained by estimating the multi-echelon assignment membership values of a set of high-order fuzzy sets.In particular, when the hierarchical fuzzy graph is updated at each generation, the individuals can be rapidly adjusted in the multi-task evolutionary optimization process.
Based on the hierarchical fuzzy graph representation of the solution, the initial solutions for each individual are generated echelon by echelon from the bottom to the top routing tasks.First, every customer (the lowest echelon nodes) is assigned to the adjacent upper nodes by a roulette wheel selection mechanism based on the corresponding 1-order fuzzy set under the concerned constraints.The routing problems in the lowest echelon are solved by the savings algorithm [24] for each adjacent upper node.With the estimated demand at each adjacent upper node, the adjacent upper echelon routes are constructed in the same way as for the corresponding routing tasks.Similarly, we can obtain the initial solution for the complete routes containing multiple task routes.

Associated task selection
To alleviate possible negative transfer caused by different routing sub-task feasible regions, a task selection strategy is designed by fully utilizing the global potential assignment information based on a hierarchical fuzzy graph.As for the multiple constraints of an ME-LRP, the association between the adjacent echelon tasks is relatively straightforward, while all other associations between tasks need to be passed through other tasks.We define a set of assignment membership matrix A j , j = 1, 2 . . .K for task T j .For each routing task T j , X d j−1 is the demand vectors of the lower n l j nodes, and X c j is the capacity vectors of the upper n u j nodes.Assuming that the demand for freights for multi-echelon tasks is met, we can obtain the systems of linear equations for K tasks as follows: Given the freight balance between each adjacent echelon task, the direction of transfer between task T i and task T j can be determined by analyzing Eq. ( 20).The similarity of feasible areas between task T i and task T j can be represented by Eq. ( 21): Next, we introduce the augmentation matrix to analyze the similarity between the feasible solutions of task T i and task T j , for adapting the cross-task probability, the element of the parameter K × K rmp matrix is defined as Eq. ( 22): where r (.) refers to the rank of the matrix or square matrix.If r (A j A j−1 . . .A i ) < r (B), the system of linear equations in Eq. ( 21) has no solution, the individuals will not consider the task transfer between task T i and task T j .If r (A j A j−1 . . .A i ) = r (B) = n u j , the system of linear equations in Eq. ( 21) has a unique solution, the individuals can directly tackle the transfer issue from a feasible area-related perspective.

Cross-task evolution operators
The hierarchical fuzzy graph for individuals of the current generation, which can indicate the multi-echelon assignment relationship, is applied to adaptively guide the cross-task evolution operators for different routing tasks.The pseudocode of the cross-task operators is presented in Algorithm 2. Based on the hierarchical fuzzy graph representation, this work realizes knowledge sharing across tasks to generate new offspring via the crossover and mutation operator.
Following the idea of a single-point crossover in the genetic evolution process, we randomly choose a crossover echelon position cep, cep ∈ [min(τ p l 1 , τ p l 2 ), max(τ p l 1 , τ p l 2 )], and perform the crossover operator over two selected individual p 1 1 and p l 2 .Figure 5 shows an example of the crossover operator between individuals p l 1 and p l 2 , where individuals p l 1 and p l 2 are based on the hierarchical fuzzy graph and p new l 2 with different crossover echelon position N p (right panel).As shown on the top right, when N p = 3, the upper task associated with the nodes at crossover echelon position is task 3 and the lower task is task 2. Since individuals p l 1 and p l 2 have the same nodes (D 2 and D 4 ) at the crossover echelon position, the new individual p new l 1 generated by the crossover operator retains a part of individual p l 1 for the task 3 and task 4, and parts of individual p l 2 for task 1 and task 2. The new individual p new l 2 generated by the crossover operator retains a partial solution of individual p l 1 for task 1 and task 2, and parts of individual p l 2 for task 3 and task 4. As shown on the bottom right, when N p = 4, the upper task associated with the nodes at crossover echelon position is task 2 and the lower task is task 1.Since individual p l 1 has three nodes (R 1 , R 3 and R 5 ) at the crossover echelon position, and individual p l 2 has two nodes (R 1 and R 5 ) at the crossover echelon position, the crossover operation requires an appropriate adjustment of the solution of the upper task 2, which is associated with the nodes in the crossover echelon position, rather than performing the crossover operation directly as shown on the top right.The new individual p new l 1 generated by the crossover operator retains a partial solution of individual p l 1 for task 3 and task 4, as well as a partial solution of individual p l 2 for task 1.In addition, we close the transfer node R 3 while a partial solution of individual p l 1 for task 2, and then apply the fuzzy neighborhood selection strategy [36] to update the routes for task 2. Similarly, the new individual p new l 2 generated by the crossover operator retains a partial solution of individual p l 1 for task 1, as Moreover, a local search strategy [17] is employed as the mutation operator to generate new offspring.The mutation operator is only applied to the solution in the effective task for the individuals, which is performed for the routing problems within a certain high-order fuzzy subset separately.Therefore, it does not change the demand assigned to upper nodes, there is no impact on other routing tasks.

Evaluation and selection
The pseudocode of fitness evaluation is listed in Algorithm 3. Based on the hierarchical fuzzy graph, an ME-LRP can be represented as a multi-task optimization problem N is the number of the routing tasks in the ME-LRP.As for the evaluation of offspring, only a task on which the individual p(x) performs best is chosen for evaluation during the evolutionary optimization.To facilitate understanding, the corresponding objective function value of individual p(x) for each task T k is provided in Eq. ( 23): where u∈ω F u y u can be calculated by based on the hierarchical fuzzy graph, ω denotes the set of upper nodes, which is open for the k-th routing task.r j ∈Lf k i ( f k +q k j )z j represents the total cost of the routes in the list Lf k i .After the evaluation of offspring individuals, we incorporate the offspring population O into the current population P and select the best n individuals with scalar fitness to update the current population P.Moreover, for the offspring population generated from cross-task operators, the skill factor τ i of the individual is set in imitation of either parent randomly, which can enhance the information sharing between routing tasks.

Update of the hierarchical fuzzy graph
To fully reflect the available hierarchical fuzzy graph information in evolutionary multi-tasking optimization, the multiechelon assignment relationship in the hierarchical fuzzy graph needs to be updated in each generation.We consider the selection of multi-echelon nodes on the hierarchical fuzzy graph G HF as a series of different fuzzy events Z i , 1 ≤ i ≤ N + 1.For each event Z i , the sample space is denoted as Ω i = {y i |1 ≤ i ≤ l i }, where l i is the number of the i-th echelon nodes.To determine membership values of multiechelon assignment relationships for high-order fuzzy sets in the hierarchical fuzzy graph, we use the idea of fuzzy probability statistics, where the membership frequency tends to be a stable value in a large number of experiments, and each individual in the population denotes a random experiment.
The assignment relationship in multi-echelon routing tasks is related to fuzzy membership value in a series of highorder fuzzy sets.As shown in Fig. 3, the potential assignment relationship for the routing task_1, task_2, and task_3 and task_4 corresponds to the construction of 1-order, 2-order, 3-order and 4-order fuzzy sets, respectively.For simplicity, we use the estimated fuzzy membership value in 2-order fuzzy sets as an example to illustrate the update of higherorder fuzzy graphs.We define three fuzzy events Z1 , Z 2 and Z 3 for three echelon nodes of two routing tasks from the bottom to the top.As discussed in the previous section, the assignment relationship between adjacent routing task_1 and routing task_2 are not independent, so we employ the Bayesian formula [13] to estimate the value of the membership values in 2-order fuzzy sets for routing task_1 and task_2 separately.
where P(Z 2 j |Z 3 i ) means the posterior assignment probability where node j on the second echelon is also assigned under the condition that node i on the third echelon is available in the solution.P(Z 1 k |Z 2 j ) denotes the posterior assignment probability, where node k on the first echelon is also assigned under the condition that node j on the second echelon is available in the solution.P(Z 2 j ) = 1/l 2 refers to the post assignment probability where node j on the second echelon is assigned in the solution.P(Z 1 k ) = 1/l 1 refers to the post-assignment probability, where node k (customer node) on the first echelon is assigned in the solution.P(Z 3 i |Z 2 j ) and P(Z 2 j |Z 1 k ), which represent the probability of the lower node selection in the solution on the routing task_2 and routing task_1, respec-tively, can be calculated with a fuzzy statistical method [4] by analyzing the hierarchical fuzzy graph representation of the solution of the current individuals.
refers to the fuzzy conditional probability of the lower node for routing task_2 to be opened for the current generation, where N pop is the number of the current individuals.η 2 i j is equal to 1 when the upper node i and lower node j are linked in the first-row list for individual p k , otherwise, η 2 i j is equal to 0. Similarly, P(Z 2 j |Z 1 k ) can be obtained based on the current individuals for routing task_1.
Once the membership values of all fuzzy sets with different orders are estimated, the hierarchical fuzzy graph is updated while ensuring a set of constraints of an ME-LRP can be met for each generation.It is worth noting that the update of the hierarchical fuzzy graph not only estimates the assignment membership in different higher-order fuzzy sets but also enables the availability of intermediate facilities or suppliers in the ME-LRP.

Experimental results
In this section, to access the performance of the proposed MEO-HFG towards multi-echelon LRPs, comprehensive empirical studies have been conducted on two-echelon LRP instances, three-echelon LRP instances, and fourechelon LRP instances, compared with the state-of-the-art approaches.We also conducted further empirical analysis to evaluate the contribution of the hierarchical fuzzy graph representation to the proposed MEO-HFG.

Experimental settings
In this study, we employ one-echelon LRP, two-echelon LRP, three-echelon LRP, and four-echelon LRP benchmarks to evaluate the proposed MEO-HFG, containing 48 instances with different echelon routing tasks.The number of instances for different echelons of LRPs is 12.Although MEO-HFG is proposed for solving multi-echelon LRPs, one may still be interested in the performance of one or two-echelon LRPs.One-echelon LRP benchmark 1 is a public dataset, and two-echelon LRP, three-echelon LRP, and four-echelon LRP benchmarks are generated based on the data derived from Prins et al. [27], which include the instances with 100 or 200 customers and 10 depots, each instance containing 1000 or 2000 potential assignment relationships.We replace the depots with regional centers for LRP benchmarks with two, three, and four echelons routing tasks.Similar to [9], eight distribution centers, five plants, and three suppliers are generated for LRPs of different echelons, so the num- To evaluate the proposed MEO-HFG, we compare it with GRASP [6], GFEA [36], and EEMTA [10].GRASP is a distance-based assignment heuristic approach with an iterative learning process.GFEA is a fuzzy evolutionary algorithm, where the graph-based assignment strategy is used to guide the fuzzy local search process.EEMTA is a state-of-the-art multi-task optimization algorithm for routing problems, where the assignment relation is obtained by K-means based on distances, and the knowledge transfer for routing generation is based on the best solutions for each routing optimization task from bottom to top.In our MEO-HFG, the population size N pop is set to 100.The maximum number of generations is set to 5000 as the termination condition for MEO-HFG, and GRASP, GFEA, and EEMTA.Other parameter settings are set according to those in the original publications.

Comparison of MEO-HFG with state-of-the-art algorithms
In this section, we present empirical evidence to show the performance of the proposed MEO-HFG and compare it to GRASP, GFEA, and EEMTA on ME-LRP instances with different echelon routing tasks.Table 2  instances with different echelon routing tasks.To identify differences between MEO-HFG and the compared algorithm on all instances over 30 independent experiments, a non-parametric test at 5% significance levels is conducted."w/d/l" indicate the number of instances for "win-draw-lose" of MEO-HFG versus other compared algorithms.When compared with GRASP, MEO-HFG significantly outperforms in almost all instances of ME-LRP.Besides, MEO-HFG significantly outperforms GFEA on 38 instances and outperforms EEMTA on 42 instances.Moreover, Table 3 illustrates more details of the computational results of MEO-HFG compared with GRASP, GFEA, and EEMTA on all ME-LRP instances.

provides the total average cost of statistics performance comparisons on all
To better present the experiment results, we divided the ME-LRP instances into four parts, including one echelon LRPs, two-echelon LRPs, three echelon LRPs, and four echelon LRPs.The columns "Average", "Std", and "T (s)" stand for the average total cost, the standard deviation of the cost, and the average run time.The average cost results with " * " indicate that the corresponding algorithm is significantly better than all other algorithms under comparison.MEO-HFG outperforms GRASP in almost one-echelon LRP instances as well as two-echelon LRPs, three-echelon LRPs, and four-echelon LRPs.MEO-HFG is comparable to GFEA on seven instances of one-echelon LRPs, 10 instances of two-echelon LRPs, 10 instances of three-echelon LRPs, and 11 instances of four-echelon LRPs.Compared to EEMTA, MEO-HFG performs better on 9 instances of one-echelon LRP problems, 10 instances of two-echelon LRP problems, 11 instances of three-echelon LRP problems, and 12 instances of four-echelon LRP problems.These results demonstrate the advantage of MEO-HFG in terms of the average cost for ME-LRPs with different echelon routing tasks.To be specific, as the number of echelons for routing tasks increases, the superiority of our approach becomes even greater.The reason might be that multi-echelon assignment in the hierarchical fuzzy graph may be more effective to improve the solution in different routing tasks.Furthermore, MEO-HFG is also competitive in terms of standard deviations for all instances with different echelon routing tasks in an acceptable computation time.The average standard deviation of all ME-LRP instances is just 118.25 for MEO-HFG, while those obtained by GRASP, GFEA, and EEMTA are 211.38,119.72, and 211.61, respectively.In addition, Fig. 6 gives the boxplots of MEO-HFG, GRASP, GFEA, and EEMTA on the four-echelon LRP instances.Compared with other algorithms under comparison, we can observe that the proposed MEO-HFG can obtain a solution with smaller cost fluctuations for most four-echelon LRP instances.The superior performance of MEO-HFG is expected since knowledge transfer of multiechelon assignment in evolutionary multi-tasking routing optimization makes it possible to avoid unpromising moves  The average results with "*" means that the corresponding algorithm is significantly better than all other compared algorithms 123

Role of hierarchical fuzzy graph
As discussed in "Multi-tasking evolutionary optimization via a hierarchical fuzzy graph", the hierarchical fuzzy graph is employed to estimate the uncertain assignment relationships for multi-echelon routing tasks.In this section, we focus on the performance of multi-echelon assignments with the hierarchical fuzzy graph for multi-echelon location routing problems.For this purpose, MEO-CPLEX as an MEO-HFG variant is employed.The main difference between MEO-HFG and MEO-CPLEX is that the assignment relationship between the upper nodes and lower nodes for each routing task in MEO-CPLEX is certain and obtained by the CPLEX solver [9] instead of using the hierarchical fuzzy graph representation for multi-echelon assignment.
Table 4 provides the results of MEO-HFG and MEO-CPLEX on four-echelon LRP instances.The column "Best" " Average" and "Std" shows the best cost of MEO-HFG and MEO-CPLEX, respectively."+", "−" and "=" indicate that MEO-CPLEX is significantly better than, comparable to, and worse than the results achieved by MEO-HFG.From Table 4, it can be found that MEO-HFG is more capable of achieving the best solutions than MEO-CPLEX in solving instances with four echelon routing tasks, which might be attributed to the guidance of the hierarchical fuzzy graph.For instance, MEO-HFG is more likely to obtain the best solution than MEO-CPLEX on 9 out of 12 four-echelon LRP instances.More importantly, MEO-HFG shows advantages in solution stability compared with MEO-CPLEX.For example, the average standard deviation of four-echelon LRP instances is 126.95 for MEO-HFG, while that obtained by MEO-CPLEX are 253.06.A possible explanation for these observations might be that the formulation of the hierarchical fuzzy graph can facilitate the generation of efficient multi-echelon assignment relationships in different routing tasks.
To further analyze the impact of the hierarchical fuzzy graph for multi-task evolutionary optimization, we evaluate the cost of the solution over the iterations on two typical four-echelon LRP instances.Figure 7 shows the best cost obtained by MEO-HFG and MEO-CPLEX over the generations, where the hierarchical fuzzy graph assignment is employed by MEO-HFG, and a given assignment for each routing task is obtained by MEO-CPLEX in MEO-CPLEX.The experimental results on instances (100-10-8-5-3-6 and 200-10-8-5-3-11) can be divided into three stages.MEO-CPLEX can get solutions with a lower cost than MEO-HFG in the first phase, but MEO-HFG can catch up with MEO-HFG in the second phase.The reason is that the hierarchical fuzzy graph assignment relationship could greatly affect the performance of the population evolution in an iterative learning process.In particular, the results of MEO-HFG in the three stages indicate a better trend than MEO-CPLEX on 200-10-8-5-3-11 in the three phases.In general, with the multi-echelon assignment knowledge of the hierarchical fuzzy graph, MEO-HFG can search for the optimal solution more quickly than MEO-CPLEX.A possible reason might be that MEO-HFG takes advantage of the fuzzy graph assignment knowledge in solving multi-echelon location routing problems based on cross-task information exchange in multitask evolutionary optimization.

Effects of evolutionary multi-tasking for ME-LRPs
In this section, we investigate whether the proposed multitask evolutionary optimization (MEO) scheme is indispensable for MEO-HFG.Two variants called EA-HFG and MEO-R-HFG are considered.As for EA-HFG, we regard ME-LRP as a single task and solve it by a traditional EA instead of using a multi-task evolutionary optimization algorithm in MEO-HFG.MEO-R-HFG is similar to MEO-HFG, but it uses a random task selection strategy.EEMTA, which is also an evolutionary multitasking algorithm, is adopted to compare with MEO-HFG, where knowledge transfer in EEMTA is limited to high-quality solutions of certain echelon routing tasks.For a fair comparison, we use the same initial individuals for the four algorithms under comparison.To further assess the search efficiency of the proposed MEO-HFG, the search convergence traces of EEMTA and two variants of the EA-HFG on representative one-echelon, two-echelon, three-echelon, and four-echelon LRP instances are presented in Fig. 8. Y-axis gives the averaged objective values in log scale, while the X-axis denotes the computational cost incurred in terms of the number of generations.
Owing to the knowledge transfer in different routing optimization tasks based on the representation of hierarchical fuzzy graph, the proposed MEO-HFG is found to converge comparably to or faster than EEMTA, EA-HFG, or MEO-R-HFG on most ME-LRP instances from the results presented in Fig. 8.For instance, on "100-10-3a", MEO-HFG uses only 300 generations to nearly arrive at the solution obtained by EEMTA, EA-HFG, and MEO-R-HFG over 400 generations.On "200-10-8-10", MEO-HFG uses only about 500 generations to achieve the solution found by EEMTA, EA-HFG, and MEO-R-HFG at 600, 650, and 550, respectively.Owing to the increase in the number of the echelon of routing tasks, we note that the advantage of MEO-HFG in convergence speed gradually declines.However, compared to EEMTA and EA-HFG, EA-HFG is still advantageous.For example, on the three-echelon LRP instances "200-10-8-5-11", it takes about 500 generations of MEO-HFG to achieve the solution obtained by EEMTA, EA-HFG over 700 generations.More importantly, when adaptive task selection is adopted for multi-echelon routing tasks, it is observed that in most of the multi-echelon routing instances, MEO-HFG converges much faster than MEO-R-HFG.On a four-echelon LRP instance, 200-10-8-5-3-11, the proposed MEO-HFG takes only around 750 generations to reach the solutions.By contrast, MEO-R-HFG uses about 850 generations.This again confirmed the efficacy of the adaptive task selection in evolutionary multitask optimization compared with MEO-R-HFG.

Conclusion and future work
In this paper, we present an MEO-HFG for multi-echelon location routing problems.MEO-HFG adopts a hierarchical fuzzy graph to represent the multi-echelon uncertain assignment relations in ME-LRPs, thereby avoiding excessive fitness evaluations of unpromising movements in multiple intermediate facilities.In addition, MEO-HFG employs a multi-task evolutionary algorithm for solving LRPs with Having shown the benefit of the hierarchical fuzzy graph based representation in improving the performance of the multi-task optimization methods, our future work will aim to further analyze the impacts between high-order fuzzy sets for solving other complex variants of LRP problems, such as dynamic LRP problems with time windows [31] and multi-objective LRP problems [3,20], as well as employing the hierarchical fuzzy graph representation to some manytasking evolutionary algorithms.In addition, we would also like to investigate an initialization scheme with reinforce-ment learning to improve the multi-echelon assignment in the ME-LRP problem as the hierarchical fuzzy graph.

Fig. 1
Fig. 1 Examples of feasible solutions to a two-echelon LRP, a threeechelon LRP, and a four-echelon LRP.The solution to the three location routing problems can be represented by a series of routes on multi-

Algorithm 1 2 :
Pseudocode of MEO-HFG Input: an ME-LRP instance; the number of individuals N pop Output: The best-so-far solution Q r 1: Formulate the hierarchical fuzzy graph G h f for the ME-LRP instance.Generate an initial population of individuals P init based on the hierarchical fuzzy graph representation.3: Evaluate the individuals with respect to multi-echelon sub-routing tasks in the multitasking optimization environment.4: Compute the skill factor τ of each individual.5: P cur = P init .6: while stopping conditions are not satisfied do 7:

Fig. 2
Fig.2The overall framework of the proposed MEO-HFG Fig.3An example of the hierarchical fuzzy graph representation for four-echelon LRP, which is composed of one 4-order fuzzy set, two 3-order fuzzy sets, three 2-order fuzzy sets, and four 1-order fuzzy sets.Associated high-order fuzzy sets are employed to express the uncertain multi-echelon assignment relationship.The weights in this graph represent the values of assignment relationship membership, respectively

Fig. 4
Fig. 4 Solution representation with a hierarchical fuzzy graph.a An example of the hierarchical fuzzy graph representation for multiechelon assignment relationships.b Cross-list storage of the feasible

Fig. 5
Fig. 5 An example of the crossover operator between the individuals p 11 and p l2 , where N p is the echelon position of the crossover operator.If N p = 3, the selected echelon nodes will be the same, a direct graph structure crossover is performed to generate p new l1 and p new l2 , as shown

Input: n: Population size 1 : 5 :
P: Current population; 2: O: Offspring population Output: S: New Population 3: Calculate the factorial cost of offspring population O by Eq. (23).4: Incorporate the offspring population O into the current population P. Sort individuals with the skill factor in P for each task.6: Update the scalar fitness of all individuals in P. 7: Choose the top n fitness individuals from P based on their scalar.8: Evaluate the individual from P on the original multi-echelon tasks by Eq. (1).9: return P 123

Fig. 6
Fig. 6 Box plots in the instances of four-echelon LRP in comparison with three different algorithms, where 1, 2, 3, and 4 in x-axis stand for MEO-HFG, GRASP, GFEA, and EEMTA, respectively.Y-axis denotes the total cost for multi-routing tasks in the log scale

Fig. 7 Fig. 8
Fig. 7 Average cost of the individuals with different multi-echelon assignment relationships obtained by MEO-HFG and MEO-CPLEX over generations

Algorithm 2
Cross-task evolution operatorsInput: p l1 , p l2 : two random selected individuals; G h f : the hierarchical fuzzy graph Output: p new l1 , p new l2 : two generated individuals 1: Record skill factors of p l1 and p l2 as τ p l 1 and τ p l 2 separately.2:obtain the mating probability rmp in Eq. (22) based on hierarchical fuzzy graph G h f .3:if τ p l 1 == τ p l 2 or rand < rmp(τ p l 1 , τ p l 2 ) then 4:Choose a random integer cep as a crossover echelon position, cep ∈ [min(τ p l 1 , τ p l 2 ), max(τ p l 1 , τ p l 2 )].
l2 representation in Fig. 4a (left panel), and the graph structure crossover is performed to generate two new offspring p new l 1

Table 1
Properties of the ME-LRP instances in multi-echelon routing tasks The properties of the ME-LRP instances in multiechelon routing tasks are given in Table1.It lists instances of benchmarks concerning the number of the lower nodes l k of each echelon routing tasks, per-unit cost transportation costs t k , the fixed costs for using vehicle f k , the fixed cost of upper nodes (supplies or intermediate facilities) F k and the capacities of vehicles U k for different echelon routing task_k.

Table 2
Statistics of performance comparison on ME-LRP instances with different routing tasks

Table 4
Results obtained by MEO-HFG and MEO-CPLEX on four-echelon LRP instances