Abstract
In recent years, low-carbon supply chain network design has been the focus of studies as the development of low-carbon economy. The location-routing problem with full truckloads (LRPFT) is investigated in this paper, which extends the existing studies on the LRP to full truckloads problem within the regional many-to-many raw material supply network. A mathematical model with dual objectives of minimizing total cost and environmental effects simultaneously is developed to determine the number and locations of facilities and optimize the flows among different kinds of nodes and routes of trucks as well. A novel multiobjective hybrid approach named NSGA-II-TS is proposed by combining a known multiobjective algorithm, NSGA-II, and a known heuristics, Tabu Search (TS). A chromosome presentation based on natural number and modified partially mapping crossover operator for the LRPFT are designed. Finally, the computational effectiveness of the hybrid approach is validated by the numerical results and a practical case study is applied to demonstrate the tradeoff between total cost and CO2 emission in the LRPFT.
1. Introduction
Facility location is one of the important problems in supply chain network design (SCND), while the cost and carbon emissions associated with the driving route of vehicles contribute a lot to the total values of a logistics system. Therefore, it is a vital necessity to combine the facility location problem (FLP) and the vehicle routing problem (VRP) in SCND.
Location-routing problem (LRP), a combination of the FLP and the VRP, has been discussed in many literatures [1–3], and the corresponding model and solving approach can be used as a principal tool for making strategy and tactics decisions. Subsequently, various variants of the LRP emerge [4], such as multi-echelon LRP [5], LRP with pickup and delivery [6], LRP with time windows [7], transportation LRP [8], and dynamic LRP [9]. However, similar to the VRP, those studies almost focus on distribution system, particularly, non-fully loaded shipments.
In practice, full truckload shipment is more common in raw material supply logistics system like timber or mineral. Although they may be transported to plants from harvesting/mining sites directly (direct transportation mode) or to logistics centers (LCs) firstly and then to plants (transit transportation mode) [10], full truckload shipment is the main form of transportation. Besides, because of the dependence on raw material, plants usually receive goods from both direct supply sites and LCs simultaneously.
The incorporation of carbon emissions into SCND is a relatively recent topic addressed in research works and most related research just took one-way trips into consideration [11–13]. However, one-way trips are ideal situation, since undesirable empty mileages can be hardly eliminated totally in trucks routes [14]. Usually, shipments are served individually and a truck is sent back empty to its next shipment after unloading at the destination of its previous trip [15].
Compared with most LRP research that employs vehicles to serve the distribution system, we study a variant of LRP where trucks are utilized to serve the regional raw material supply logistics system that has two kinds of transportation modes: direct and transit. We call this variant of LRP as location-routing problem with full truckloads (LRPFT).
The LRPFT arises from bulk raw material transportation in supply logistics system. There are several supply points (suppliers), demand points (plants), and a few potential LCs in the network. The LRPFT differs from classic LRP mainly in two aspects. Firstly, the structure of supply chain network (SCN) is different (see Figure 1). Raw material may be transported to plants from supply sites directly or transited by the LC and a plant may be served by different LCs simultaneously. This kind of SCN usually appears when demand points are resource-dependent businesses (i.e., biomass processing plants). Demand points need a large amount of raw material. To acquire more, they not only receive raw material with preprocessing (i.e., Chips) from LCs but also receive raw material without any preprocessing (i.e., fuel wood) from supply sites.
(a) Classic SCN
(b) Bulk raw material SCN
Secondly, the truck routing problem incorporated in the LRPFT is a variant of classic VRP. It belongs to full truckload routing problem. Each route starts from an opened LC and ends at the same LC, while each node can be visited by different trucks many times, which is forbidden in classic VRP. Figure 2 illustrates a solution example of six routes for the truck routing problem addressed in this study. The six routes are the following: route 1, L1-H1-P1-H2-P2-L1; route 2, L1-H1-L1-H1-L1; route 3, L1-P3-H4-P4-L1; route 4, L3-P4-H3-L3; route 5, L3-P4-L3; and route 6, L3-P4-H5-P1-H5-P1-L3. Bold arcs indicate loaded trucks and others unloaded. A truck probably runs several times between two nodes (e.g., L1 and H1 in route 2 and H5 and P1 in route 6).
The LRPFT considered in this study is a deterministic decision problem. Demand points can receive raw material from suppliers and/or LCs. The raw material that comes from suppliers directly will receive preprocessing at the demand point while others are processed at LCs. At the beginning, all trucks are assigned to the opened LCs where they must return to after each route. Finally, due to some limitations such as drivers’ working hours, a route must not exceed a given distance-span. The LRPFT is to determine which LCs to open, how to assign flows, and the scheduling and routing of trucks. Unlike most LRPs, the LRPFT is based on arc routing.
In this study, we investigate the LRPFT with two objectives: total cost minimization and total carbon emissions minimization. A comprehensive mathematical model is developed firstly. Since the LRPFT is a more complex problem than LRP that is a well-known NP-hard problem, a solving method for the LRPFT is designed as well. The remainder of this paper is organized as follows: Section 2 begins with a comprehensive review of related literatures. Section 3 introduces the mathematical formulations for the problem. Pareto optimality concept is provided and the algorithm is proposed in Section 4. Numerical experiments to validate the mathematical model and evaluate the algorithm are presented in Section 5. Finally, conclusions are drawn in Section 6.
2. Literature Review
The classical LRP combines two basic planning tasks in logistics and is defined as an approach to modelling and solving locational problems [1]. For detailed information about LRP, see [1, 2, 4].
The LRP has many real-life applications so that more than 10 kinds of its variants can be summarized [4]. Compared with the large number of LRP papers that address node routing, those considering arc routing are rare. Although the arc routing problem can be transformed into an equivalent node routing problem, this type of transformation increases the problem size and it needs additional calculation for a VRP [16], and this is even more true for a LRP [17]. Levy and Bodin [18] first worked on location-arc routing problem (LARP). They tackled a practical problem arising in the scheduling of postal carriers in the United States postal service. Some researches defined the LARP as an extension of arc routing problem and developed different heuristic algorithms for it [17, 19, 20]. Muyldermans [21] presented a variant of the LARP: the dead mileage problem, in which clients were allowed to be serviced multiple times, and the objective was to minimize dead-mileage.
Multi-echelon LRPs have only very recently attracted the interest of researchers [22], and the two-echelon LRP (2E-LRP) is close to the problem studied in this work to some extent. In the 2E-LRP, three disjoint sets of nodes representing suppliers (first-echelon facilities), LCs (second-level facilities), and customers are considered [5]. Decisions are the location of a subset of suppliers and LCs and vehicle routes for both the 1st-echelon and 2nd-echelon vehicles. Each 1st-echelon vehicle starts from the supplier to LCs while each 2nd-echelon vehicle is routed from the LC to customers [23]. Obviously, they are two mutually independent routing problems. Even though Nguyen et al. [24] took into account directive delivery customers in their study on 2E-LRP. The routing problems are separated completely for different parts. The LRPFT differs from the 2E-LRP mainly on three aspects. Firstly, the routing problem in the LRPFT is a comprehensive one. Trucks may visit suppliers, LCs, and customers in one route. Secondly, the LRPFT has a more complicated decision structure. The allocation of flows between nodes can be divided into two echelons in the 2E-LRP while the LRPFT integrates them. Finally, the routing problems are usually the non-fully loaded ones in the 2E-LRP.
The full truckload routing problem derived from the LRPFT has been considered in literatures, although it is rare. For example, Liu et al. [25] proposed a mathematical programming model for the full truckloads multidepot vehicle routing problem and designed a two-phase greedy algorithm to solve it. Derigs et al. [26] applied various variants of multilevel neighborhood search strategies for solving two specific classes of full truckload routing problems arising in timber transportation. However, the full truckload routing problem has a wide application. The truck container transport is a typical example of the full truckload routing problem since a container can represent a full truckload [27].
As concerns on low-carbon supply chain management increase, the incorporation of the minimization of carbon emissions into the LRP becomes an important topic addressed in researches in recent years. Govindan et al. [23] studied a multiobjective LRP which arises from the application in perishable food SCND. Two minimization objectives are pursued: total cost including fixed facility opening cost and variable routing cost and total CO2 emission of the system. Toro et al. [28] formulated a biobjective problem, named green capacitated LRP (G-CLRP), considering the minimization of operational costs and the minimization of total emissions associated with the fuel consumption.
3. Mathematical Formulation
Let be a directed graph, where is the vertex set and is the arc set. , , and are sets of supply sites, potential LCs, and plants, respectively. Each supplier has a specific supply volume; each LC has a fixed capacity and the demand of each plant must be satisfied. The location of a subset of LCs and the truck routes among the three kinds of nodes are to be determined. The LC conducts preprocessing for plants, which causes weight losses of raw material. Each plant has its own loss rate. Assumptions, parameters, and variables are defined below.
3.1. Assumptions
The underlying assumptions of the LRPFT model include the following:(i)The LRPFT is a deterministic decision problem.(ii)A plant can be served by multiple LCs and suppliers.(iii)A truck can visit suppliers, LCs, and plants multiple times in a single route.(iv)A route must not exceed a given distance-span.(v)All trucks are assigned to the opened LCs where they must return to after each route.(vi)Homogeneous vehicle type is considered.
3.2. Notations
(1) Sets and Parameters : the set of vehicles : the supply of supplier , kg : the demand of plant , kg : unit transportation cost for a vehicle travelling from node to node , $/trip : unit CO2 emission for a loaded vehicle travelling from node to node , kg CO2/trip : unit CO2 emission for an unloaded vehicle travelling from node to node , kg CO2/trip : unit handling cost in LC , $/kg : unit cost for preprocessing at node , $/kg : capacity of the vehicle in use, kg : fixed cost of a vehicle for a single use, $/route : fixed cost for LC , $ : max capacity of LC , kg : fixed CO2 emission of LC , kg CO2 : unit CO2 emission for preprocessing at node , kg CO2/kg : processing loss rate of raw materials for plant : maximum distance-span for a route, km : a sufficient large positive number
(2) Decision Variables : the total amount of load from node to node , kg : times of loaded vehicle moving along arc , trips : times of unloaded vehicle moving along arc , trips
3.3. Formulations
s.t.
We explicitly consider two objective functions. OBJ1 minimizes the total cost of the SCN. The first term is the transportation cost, and the second term is the fixed cost of using vehicles. The third term is the fixed cost of opening LCs, and the fourth term is the total handling cost. The fifth term is the total cost of preprocessing at both LCs and plants. OBJ2 measures the total CO2 emissions associated with logistics activities throughout the network. The first term is the total emissions related to transportation activities, which are arc-dependent. In this study, CO2 emissions on the same arc are distinguished between heavy and empty vehicle. The second term is the fixed emissions associated with opening LCs, and the third term is the variable emissions associated with preprocessing at LCs and plants. The conflict between OBJ1 and OBJ2 exists mainly because of the following two aspects. (1) opening more LCs reduces the transportation needed in the network, and the fixed cost of a LC is much higher than the transportation cost, while it is the opposite case in carbon emissions [13]. (2) Because the transportation cost is not proportional to the transportation emissions even on the same arc, routes decisions under minimum cost and minimum emissions are quite different [29]. Transportation emissions depend on the fuel consumption of vehicles totally and the fuel cost is only a fraction of the transportation cost.
Constraint (2) ensures that the amount of raw material flowing out of supplier should not exceed its supply capacity, while constraint (3) restricts the demand satisfaction of each plant. Constraint (4) states that a LC is disabled if it is not open and the total raw material flowing into an opened LC should not exceed its capacity. Constraint (5) guarantees that the trips of heavy vehicles on an arc match the volumes of raw materials to be transported. Constraints (7), (8), and (9) guarantee the continuity of each route. Equation (10) indicates that no empty vehicle travels to plants and (11) indicates that no heavy vehicle travels to suppliers. Inequality (12) prohibits subtours. In this case, isolated subtour elimination constraints are different from famous subtour elimination constraints in classic VRP. For example, route 2 in Figure 2 is legal in full truck-load VRP, but it is forbidden in classic VRP [25]. Constraint (13) implies the maximum distance-span of each route and constraint (14) guarantees that vehicles travelling between nodes must be in use. Variables’ types are defined in constraints (15)–(18).
The model was validated by a simple example in Section 5.1 and it is observed that, even for small instances, the solver consumes a lot of time. Therefore, a novel hybrid heuristic algorithm is designed to solve the problem.
4. Solution Approach
The Pareto nondominated solution set depicts the tradeoff between objectives in multiobjective optimization and the concept of nondominance can be found in Moslemi and Zandieh [30]. NSGA-II, developed by Deb et al. [31], is one of the representative evolutionary algorithms for multiobjective optimization problems to obtain Pareto solution set and has been modified and applied to different fields successfully [32, 33]. However, several operators in the algorithm need to be designed elaborately when applied to a specific problem, such as encoding and decoding of the chromosome, crossover operator, and mutation operator. The TS is a metaheuristic technique characterized as a local search method. For complex combinatorial optimization problems, combining certain aspects of different algorithms to form a hybrid algorithm is a promising approach due to increases in computational power. Accordingly, the combination of TS as an algorithm enhances the ability of local search of the search process, wherein TS search is restarted with the most preferred solutions. NSGA-II is a popular metaheuristic technique used for its capability of exploring a wide solution space, which paves the ground to avoid falling into the local optimal. In single objective optimization, the hybrid algorithm based on genetic algorithm and TS has shown its excellence [34, 35]. However, the application of such hybridization in multiobjective optimization is rare. In this study, we adopt NSGA-II to the proposed model and develop a hybrid NSGA-II named NSGA-II-TS, in which the TS is incorporated to improve the quality of nondominated solution set.
4.1. Main Loop of NSGA-II-TS
NSGA-II is fundamentally a genetic algorithm with special characteristics in the selection process, which is based on a partial order. The partial order is determined by two metrics: nondomination rank and crowding distance.
The aim of nondominated sorting is to assign a nondominated rank to every solution. Within a population, the solution with the higher nondominated rank has higher priority to be selected and evolved to next generation.
Crowding distance is an index of diversity, which measures the proximity of a solution to its nearest neighbors. A solution with smaller crowding distance is more crowded and the solution located in a less crowded region is with higher preference in the selection process. The introduction of crowding distance helps the Pareto front formed by evenly distributed nondominated solutions, which is preferred in multiobjective optimization.
In evolution operator, the parents and children are combined together in a set firstly; then the solutions with higher nondominated ranks are evolved to the next generation, and if some solutions have the same ranks, the solution with bigger crowding distance has higher priority.
Based on the basic framework of the NSGA-II, we incorporate the TS procedure. In each generation, some chromosomes in current population are selected by certain criterions (see Section 4.5.1), to which the TS procedure are applied. Since the TS procedure would improve the diversity of solutions, the mutation operator is omitted in this study for efficiency. Figure 3 shows the flow diagram of our algorithm and details will be discussed below. It is noteworthy that solutions resulted from the TS procedure are added into the current generation rather than replacing original solutions, because they may be nondominated solutions to original solutions.
4.2. Encoding and Decoding
4.2.1. Encoding
As for the proposed model, a feasible solution contains three parts: which LCs to open, how material flows distribute among the three kinds of nodes, and what the routes of trucks in use. In the study, we use nodes in the sequence of trucks’ visit to form a chromosome, because truck routes imply opened LCs and material allocations can be deduced from truck routes as well.
Furthermore, for truck loading efficiency, the number that a plant needs to be visited by trucks is calculated as . Similarly, the number that a supplier may be visited is . However, the visit number of LCs varies according to the transportation mode (Figure 4). Obviously, the transit mode results in more visit of LCs, which is just twice as many as the total of other nodes in a route. Since the origin and destination of a route must be the same LC, the destination LC is default in the chromosome to control its length. We set the total number that all LCs are visited as in our chromosomes. Consequently, the length of the chromosome is .
(a) Transit transport
(b) Direct transport
Taking a SCN with 4 suppliers (denoted by 1~4), 2 potential LCs (denoted by 5~6), and 4 plants (denoted by 7~10) as example, it may have a chromosome as shown in Figure 5 when a truckload supply/demand for each supplier/plant is assumed.
Although the proposed encoding method would result in redundant genes for LCs, it is intuitive, intelligible, and easy to implement local search and evolutional operators.
4.2.2. Decoding
The decoding procedure begins at the first gene in a chromosome. A new route is built firstly. Then the genes are checked one by one. If it is a redundant gene, go to the next directly. Otherwise, whether it can be added into the current route is checked. If yes, insert it into current route. If no, finish current route with current LC and start a new route. Repeat the above steps until all genes are checked.
For the resulted route, the truck travelling on each arc can be determined to be loaded or empty according to its from-to nodes. A truck from LCs to plants is a heavy truck, and that from LCs to suppliers is an empty truck. If an arc’s initial node is a supplier, trucks travelling on it are loaded and if it is a plant, the trucks are empty.
In terms of the chromosome shown in Figure 5, the routes obtained by decoding may be 5-1-5, 5-2-10-5, 5-3-7-5, 5-4-8-5, and 5-9-5. It is noted that a LC (node 5) is inserted between node 7 and node 4 because of the maximum route length constraint.
Because constraints on flow conservation in LCs and demand satisfaction of plants are ignored in the decoding procedure, an extra objective is introduced to manage the infeasibility of solutions, which is defined in the following.where and represent the volume of material out of and into LC , respectively, and represents the unsatisfied volume of plant .
4.3. Initialization
To generate an initial chromosome, three sets, , , , and three records introduced in Section 4.2.1 need to be kept: visit times of every demand point, ; visit times of every supply points, ; and visit times of LCs, . Then genes are generated one by one following the rules below and the pseudocode of randomly generating an initial chromosome is shown in Algorithm 1.(1)The first gene must be a LC, since vehicles originate at LCs.(2)If current gene is a LC, the next gene may be a LC, a supplier, or a plant.(3)If current gene is a supplier, the next gene may be a LC or a plant.(4)If current gene is a plant, the next gene may be a LC or a supplier.
| ||||||||||||||||||||||||||||||||||
4.4. Crossover Operator
Because of the special structure of our chromosome, it is hard to implement any conventional crossover operator directly. We modify the partially mapping crossover (PMX) [36] to adapt to the chromosome in this study. Given two parents (P1, P2) for crossover, the PMX has 4 steps, as follows.
Step 1. Randomly generate two split points for parent chromosomes, respectively. Two segments are formed: sublist 1 for P1 and sublist 2 for P2, respectively, as shown in Figure 6.
Step 2. Exchange sublist 1 and sublist 2 to format two prototypes of offspring, O1 and O2.
Step 3. Mark the genes in O1 apart from sublist 2 that are the same with genes in sublist 2, and mark the genes in O2 apart from sublist 1 that are the same with genes in sublist 1.
Step 4. Substitute the marked genes with the missing genes in O1 in the order they are encountered in sublist 1 and substitute the marked genes with the missing genes in O2 in the order they are encountered in sublist 2.
The PMX is designed for the sequential coding with no duplication code. The application of PMX on our chromosomes generates invalid chromosomes inevitably. To avoid this, we propose a modified PMX, named MPMX. There are mainly two differences between PMX and the MPMX. The first is the formation of the exchange segments. In the MPMX, the sublist 1 is generated randomly firstly, and then the sublist 2 is generated on the premise of the following exchange constraints.
As depicted in Figure 6, let and denote the immediate preceding and the succeeding gene of crossover segment in the parent chromosome, respectively; let and represent the beginning and the ending gene, respectively. Crossover operator can be applied with sublist 1 and sublist 2 if the following conditions are satisfied.(1)If is a supplier (plant), is not a supplier (plant).(2)If is the beginning gene of P1 (P2), is a LC.
The repairing procedure after the segments exchange is the second difference between the PMX and the MPMX. The infeasibility in O1 and O2 is that the frequency of a node in a chromosome does not match what is stipulated in the encoding phase. Therefore, the repairing procedure is applied to maintain the numbers of each supplier, each plant, and total LCs in an offspring prototype, which ensures its feasibility.
In terms of O1, genes in sublist 1 are the missing ones and genes in sublist 2 are the redundant ones. Missing genes and redundant genes are paired off and the pairs are checked one by one. For each pair, if the missing gene is exactly the same gene as the redundant, nothing needs to be done. Otherwise, whether the missing gene and the redundant gene represent the same kind of nodes (suppliers or LCs or plants) is checked. If yes, a gene that is the same to the redundant gene in the offspring prototype is randomly replaced by the missing one. If no, a randomly generated LC is used as media. The redundant gene is replaced by the LC and then an appropriate LC in the offspring prototype is replaced by the missing gene under the prerequisite of guaranteeing feasibility.
4.5. TS Procedure
4.5.1. Solutions Selection
For efficiency, it is necessary to determine which solutions the TS procedure is applied with. In this study, there are two criteria. Firstly, the solutions are of the highest nondominated level in current population. Besides, the following two parts are considered: solutions with lowest value of OBJ1 or OBJ2 and solutions with higher crowding distance. The former aims to enhance the search ability in extremum spaces of the algorithm, while the later focuses on the distribution of solutions. Because the solution space is discrete, a well-distributed solution set is important to depict the Pareto front precisely. Finally, a TS probability is introduced to determine the specific number of the solutions.
4.5.2. Solution Evaluation
The evaluation function of solutions in the TS is calculated using the following.where is the objective value of and is the weighted of objective . The configuration of guides the search direction. In this paper, we set for solutions with low value of OBJ1, for solutions with low value of OBJ2, and a randomly generated configuration of is considered for solutions with high crowding distance.
4.5.3. Neighborhoods Generation
A neighborhood is generated as follows: randomly choosing a gene as major exchange gene firstly, then checking its type (LC/supplier/plant), and finally randomly choosing another gene with the same type and swapping those two genes.
5. Numerical Experiments
To validate the proposed model and algorithm, a small instance is first constructed. Then, we test the NSGA-II-TS on 36 instances of different characteristics of network and compare its results with those of the NSGA-II. The algorithms are implemented in Matlab R2014a and run on PCs with 2.4 GHz and 4 GB RAM memory.
5.1. A Simple Example
The simple example contains 4 suppliers, 2 potential LCs, and 4 plants. The visit times of each supplier and plant is one. The model is solved by Lingo 11.0 using the normalized normal constraint (NNC) method [37], an exact algorithm, and the NSGA-II-TS, respectively. Figure 7 illustrates the solutions.
On the simple instance, solutions of the NSGA-II-TS are almost the same with exact solutions. Three out of five are exactly the same including two extreme points. The Inverted Generational Distance (IGD) defined in the following equation is usually used to measure the performance of a multiobjective optimal algorithm [38].where is a set of points uniformly distributed on the Pareto front and is the approximation obtained by the algorithm; is the minimum Euclidean distance between and the points in .
Given that the points in set can appropriately represent the Pareto front, IGD can measure both the diversity and convergence of . The smaller the value of IGD is, the better the solution set toward the Pareto optimal front is. When all obtained solutions lie exactly on chosen solutions, IGD takes a value of zero. In this example, the IGD metric of the NSGA-II-TS is 0.003966 that is small enough to prove the validity of the proposed model and algorithm.
5.2. Test Instances Generation
Lack of publicly available test instances for the LRPFT prevents us from comparing our results to similar results from the literature. To examine the proposed algorithm (NSGA-II-TS) further, we compare it with the NSGA-II on a range of randomly generated instances. Those instances can be classified into three scales: large, medium, and small, according to the nodes number in the network and their average visits. Accordingly, each problem is characterized by suppliers LCs lants visits . The distribution of the three types of nodes is generated randomly and the Euclidean distance is used. Following Govindan et al. [23], some coefficient matrices are introduced to calculate basic input parameters, such as , , , and . The ranges of those parameters are summarized in Table 1. The travelling cost is calculated as , and the emissions . The fixed cost of potential LC is calculated as , where and are the benchmarks of fixed cost and capacity of potential LCs ( = 35000, = 12 truckloads) and the fixed emissions of LCs .
5.3. Evaluation Metrics
Convergence and diversity are the two most important aspects in evaluating multiobjective algorithms. Convergence measures the degree to which the nondominated solution set approaches the exact Pareto front, while diversity measures the distribution uniformity of the solutions in the Pareto optimal set. This paper uses the following two performance assessment metrics: mean ideal distance (MID) for convergence and spread extent (Δ) for diversity.
MID[23]. MID calculates the mean distance between the nondominated set and ideal point . The algorithm with a lower value of MID has better convergence. The MID is computed as follows:where is the number of nondominated solutions.
[31]. measures the extent of spread achieved among the obtained solutions, which is calculated as follows:where and are the Euclidean distance between the extreme solutions and the boundary solutions of the obtained nondominated set, respectively, and is the Euclidean distance between consecutive solutions; is the average of all , assuming that there are solutions in the first front. Apparently, an ideal distribution is the one that makes all distance equal to and make . At this point, the metric Δ would be zero. For any other distribution, the value of the metric would be greater than zero. For two distributions having identical values of and , the metric Δ takes a higher value with worse distributions of solutions within the extreme solutions.
5.4. Computational Results
To determine the different performances of the NSGA-II and the NSGA-II-TS, the values of the metrics of solution sets obtained by those two algorithms are calculated as reported in Table 2. The same maximum solving time of algorithms is set and solutions with OBJ3 > 0 are excluded.
The experimental results show that the NSGA-II-TS obtains better value of both metrics than the NSGA-II does on 34 out of 36 test instances. The NSGA-II-TS gets higher value of MID or Δ in the remaining two test instances. However, it is arbitrary if any of the two metrics is used alone to measure the performance of an algorithm. We present the nondominated solutions obtained by the two algorithms on those two instances in Figure 8.
(a) The Pareto solutions of instance 4 × 2 × 4 × 5
(b) The Pareto solutions of instance 10 × 5 × 10 × 2
In the case (Figure 8(a)), the Pareto optimal front obtained by the NSGA-II-TS consists of three segments and the neighboring sections are separated by a distinct distance, which results in large values of and subsequently a pretty large Δ. Meanwhile, the Pareto optimal frontier obtained by the NSGA-II has only one segment; thus the solutions are more concentrated, which results in a smaller Δ. However, the Pareto optimal front obtained by the NSGA-II-TS gives a more complete view of the Pareto front, which provides decision makers with more diverse solutions.
Instance (Figure 8(b)) is another situation. NSGA-II achieves more and better solutions in the middle of the Pareto optimal frontier, so the value of MID is relatively small. However, the NSGA-II-TS exhibits stronger ability in searching extreme points, especially the points in the top left corner.
In general, the computational results prove that the proposed algorithm NSGA-II-TS can solve the proposed model effectively. Besides, in all test instances, the NSGA-II-TS outperforms the NSGA-II. The incorporation improves the search ability for whole solution space efficiently. Meanwhile, the distribution of solutions is improved as well.
5.5. A Practical Case Study
A practical case is motivated by the regional timber supply in Jiangle County, Fujian province, China (Figure 9). The supply network contains 7 suppliers, 4 potential LCs, and 5 demand points. In the planning horizon, the average visit of nodes is 8.
The shortest paths between each pair of nodes were computed, and the corresponding costs were determined through field investigation and the emissions were computed with road conditions (i.e., gradient, rolling resistance, and distance) and truck type in use based on the method proposed in [39].
All other data related to cost calculation is determined based on the combination of logistics practices and the information obtained from the forest service department in Jiangle County.
Table 3 provides the results of the NSGA-II and the NSGA-II-TS. The MID values are 0.90736 and 0.88217, respectively, and the Δ values are 1.03448 and 0.80734, respectively. As expected, the NSGA-II-TS outperforms the NSGA in both criterions.
Furthermore, we can see tradeoffs between the two objectives here. The solutions can be divided into groups according to opened LCs. There are several nondominated solutions in each group, which shows the tradeoff between the objectives under a certain configuration of LCs. As the cost of a route is not proportional to its emissions, an economic truck route is not necessarily the carbon emissions minimum. Emissions are mainly caused by the fuel consumption, while the cost contains fuel cost, driver salary, vehicle depreciation, and even road toll. For example, a short route with bad road condition causes higher emissions but lower cost. On the contrary, because of the toll of high-grade road in China, a route of low emissions may have high cost. Besides, as Table 3 illustrates, among groups, solutions with more opened LCs are of lower level of emissions. As also indicated in [13], more depots are opened to reduce the kilometers travelled and the fixed cost of a LC is sufficiently higher than the transportation cost.
6. Conclusion
In this paper, the biobjective location-routing problem with full truckloads is addressed. This problem is an extension of the existing studies on the LRP in distribution system to supply logistics system. The problem involves two transportation modes, direct and transit, and integrates trucks’ routes for different parts, which is more realistic for perfecting raw material supply logistics system. The mathematical model with dual aims of minimizing total cost and environmental effects is proposed to investigate the tradeoff between the objectives.
Due to the complexity of the problem, a novel hybrid metaheuristic algorithm named NSGA-II-TS that combines the NSGA-II and the TS is developed to obtain Pareto solutions, in which a new representation and a new crossover operator are designed. To validate the proposed model, a simple example is solved by the NSGA-II-TS and the exact algorithm, NNC. After that, 36 instances with different scales are generated randomly to test the performance of the NSGA-II-TS. Finally, a practical case motived in timber supply logistics system in Jiangle County, Fujian province, China, is applied as well. The computational results validate the effectiveness and efficiency of the NSGA-II-TS on the LRPFT. Besides, the numerical solution results on practical case study give the general essential qualitative relationships between total cost and carbon emissions, which provide managerial insights for decision makers.
Further research could focus on an extension of the LRPFT by including, for example, multiple vehicle types or time windows, to enhance its applicability to real-life scenarios. On the other hand, efficient multiobjective algorithms for the LRPFT can also extend the research.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This research was financially supported by the Natural Science Foundation of Fujian Province under Grant no. 2017J01788 and the University Development Foundation of Fujian Agriculture and Forest University (no. 612014018).