Research Article  Open Access
Cheng Chen, Rongzu Qiu, Xisheng Hu, "The LocationRouting Problem with Full Truckloads in LowCarbon Supply Chain Network Designing", Mathematical Problems in Engineering, vol. 2018, Article ID 6315631, 13 pages, 2018. https://doi.org/10.1155/2018/6315631
The LocationRouting Problem with Full Truckloads in LowCarbon Supply Chain Network Designing
Abstract
In recent years, lowcarbon supply chain network design has been the focus of studies as the development of lowcarbon economy. The locationrouting 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 manytomany 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 NSGAIITS is proposed by combining a known multiobjective algorithm, NSGAII, 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 CO_{2} 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.
Locationrouting 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 multiechelon 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, nonfully 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 oneway trips into consideration [11–13]. However, oneway 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 locationrouting 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 resourcedependent 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, L1H1P1H2P2L1; route 2, L1H1L1H1L1; route 3, L1P3H4P4L1; route 4, L3P4H3L3; route 5, L3P4L3; and route 6, L3P4H5P1H5P1L3. 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 distancespan. 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 wellknown NPhard 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 reallife 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 locationarc 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 deadmileage.
Multiechelon LRPs have only very recently attracted the interest of researchers [22], and the twoechelon LRP (2ELRP) is close to the problem studied in this work to some extent. In the 2ELRP, three disjoint sets of nodes representing suppliers (firstechelon facilities), LCs (secondlevel facilities), and customers are considered [5]. Decisions are the location of a subset of suppliers and LCs and vehicle routes for both the 1stechelon and 2ndechelon vehicles. Each 1stechelon vehicle starts from the supplier to LCs while each 2ndechelon 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 2ELRP. The routing problems are separated completely for different parts. The LRPFT differs from the 2ELRP 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 2ELRP while the LRPFT integrates them. Finally, the routing problems are usually the nonfully loaded ones in the 2ELRP.
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 twophase 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 lowcarbon 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 CO_{2} emission of the system. Toro et al. [28] formulated a biobjective problem, named green capacitated LRP (GCLRP), 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 distancespan.(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 CO_{2} emission for a loaded vehicle travelling from node to node , kg CO_{2}/trip : unit CO_{2} emission for an unloaded vehicle travelling from node to node , kg CO_{2}/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 CO_{2} emission of LC , kg CO_{2} : unit CO_{2} emission for preprocessing at node , kg CO_{2}/kg : processing loss rate of raw materials for plant : maximum distancespan 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 CO_{2} emissions associated with logistics activities throughout the network. The first term is the total emissions related to transportation activities, which are arcdependent. In this study, CO_{2} 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 truckload VRP, but it is forbidden in classic VRP [25]. Constraint (13) implies the maximum distancespan 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]. NSGAII, 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. NSGAII 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 NSGAII to the proposed model and develop a hybrid NSGAII named NSGAIITS, in which the TS is incorporated to improve the quality of nondominated solution set.
4.1. Main Loop of NSGAIITS
NSGAII 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 NSGAII, 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 fromto 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 515, 52105, 5375, 5485, and 595. 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 welldistributed 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 NSGAIITS on 36 instances of different characteristics of network and compare its results with those of the NSGAII. 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 NSGAIITS, respectively. Figure 7 illustrates the solutions.
On the simple instance, solutions of the NSGAIITS 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 NSGAIITS 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 (NSGAIITS) further, we compare it with the NSGAII 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 NSGAII and the NSGAIITS, 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 NSGAIITS obtains better value of both metrics than the NSGAII does on 34 out of 36 test instances. The NSGAIITS 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 NSGAIITS 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 NSGAII has only one segment; thus the solutions are more concentrated, which results in a smaller Δ. However, the Pareto optimal front obtained by the NSGAIITS gives a more complete view of the Pareto front, which provides decision makers with more diverse solutions.
Instance (Figure 8(b)) is another situation. NSGAII achieves more and better solutions in the middle of the Pareto optimal frontier, so the value of MID is relatively small. However, the NSGAIITS 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 NSGAIITS can solve the proposed model effectively. Besides, in all test instances, the NSGAIITS outperforms the NSGAII. 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 NSGAII and the NSGAIITS. The MID values are 0.90736 and 0.88217, respectively, and the Δ values are 1.03448 and 0.80734, respectively. As expected, the NSGAIITS 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 highgrade 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 locationrouting 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 NSGAIITS that combines the NSGAII 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 NSGAIITS and the exact algorithm, NNC. After that, 36 instances with different scales are generated randomly to test the performance of the NSGAIITS. 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 NSGAIITS 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 reallife 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).
References
 G. Nagy and S. d. Salhi, “Locationrouting: issues, models and methods,” European Journal of Operational Research, vol. 177, no. 2, pp. 649–672, 2007. View at: Publisher Site  Google Scholar  MathSciNet
 C. Prodhon and C. Prins, “A survey of recent research on locationrouting problems,” European Journal of Operational Research, vol. 238, no. 1, pp. 1–17, 2014. View at: Publisher Site  Google Scholar  MathSciNet
 G. Laporte, Y. Nobert, and D. Arpin, “An exact algorithm for solving a capacitated locationrouting problem,” Annals of Operations Research, vol. 6, no. 9, pp. 291–310, 1986. View at: Publisher Site  Google Scholar
 M. Drexl and M. Schneider, “A survey of variants and extensions of the locationrouting problem,” European Journal of Operational Research, vol. 241, no. 2, pp. 283–308, 2015. View at: Publisher Site  Google Scholar  MathSciNet
 C. Contardo, V. Hemmelmayr, and T. G. Crainic, “Lower and upper bounds for the twoechelon capacitated locationrouting problem,” Computers & Operations Research, vol. 39, no. 12, pp. 3185–3199, 2012. View at: Publisher Site  Google Scholar  MathSciNet
 I. Karaoglan, F. Altiparmak, I. Kara, and B. Dengiz, “The locationrouting problem with simultaneous pickup and delivery: Formulations and a heuristic approach,” Omega , vol. 40, no. 4, pp. 465–477, 2012. View at: Publisher Site  Google Scholar
 M. H. F. Zarandi, A. Hemmati, S. Davari, and I. B. Turksen, “Capacitated locationrouting problem with time windows under uncertainty,” KnowledgeBased Systems, vol. 37, pp. 480–489, 2013. View at: Publisher Site  Google Scholar
 I. A. MartínezSalazar, J. Molina, F. ÁngelBello, T. Gómez, and R. Caballero, “Solving a biobjective transportation location routing problem by metaheuristic algorithms,” European Journal of Operational Research, vol. 234, no. 1, pp. 25–36, 2014. View at: Publisher Site  Google Scholar  MathSciNet
 S. Gao, Y. Wang, J. Cheng, Y. Inazumi, and Z. Tang, “Ant colony optimization with clustering for solving the dynamic location routing problem,” Applied Mathematics and Computation, vol. 285, pp. 149–173, 2016. View at: Publisher Site  Google Scholar  MathSciNet
 J. J. Troncoso and R. A. Garrido, “Forestry production and logistics planning: an analysis using mixedinteger programming,” Forest Policy and Economics, vol. 7, no. 4, pp. 625–633, 2005. View at: Publisher Site  Google Scholar
 A. Diabat and M. AlSalem, “An integrated supply chain problem with environmental considerations,” International Journal of Production Economics, vol. 164, pp. 330–338, 2015. View at: Publisher Site  Google Scholar
 F. Wang, X. Lai, and N. Shi, “A multiobjective optimization for green supply chain network design,” Decision Support Systems, vol. 51, no. 2, pp. 262–269, 2011. View at: Publisher Site  Google Scholar
 S. Elhedhli and R. Merrick, “Green supply chain network design to reduce carbon emissions,” Transportation Research Part D: Transport and Environment, vol. 17, no. 5, pp. 370–379, 2012. View at: Publisher Site  Google Scholar
 K. Haridass, J. Valenzuela, A. D. Yucekaya, and T. McDonald, “Scheduling a log transport system using simulated annealing,” Information Sciences, vol. 264, pp. 302–316, 2014. View at: Publisher Site  Google Scholar  MathSciNet
 N. El Hachemi, M. Gendreau, and L.M. Rousseau, “A heuristic to solve the synchronized logtruck scheduling problem,” Computers & Operations Research, vol. 40, no. 3, pp. 666–673, 2013. View at: Publisher Site  Google Scholar
 H. Li, T. Lv, and Y. Li, “The tractor and semitrailer routing problem with manytomany demand considering carbon dioxide emissions,” Transportation Research Part D: Transport and Environment, vol. 34, pp. 68–82, 2015. View at: Publisher Site  Google Scholar
 R. B. Lopes, F. Plastria, C. Ferreira, and B. S. Santos, “Locationarc routing problem: heuristic approaches and test instances,” Computers & Operations Research, vol. 43, pp. 309–317, 2014. View at: Publisher Site  Google Scholar  MathSciNet
 L. Levy and L. Bodin, “The Arc Oriented Location Routing Problem,” INFOR: Information Systems and Operational Research, vol. 27, no. 1, pp. 74–94, 2016. View at: Publisher Site  Google Scholar
 G. Ghiani and G. Laporte, “LocationArc Routing Problems,” OPSEARCH, vol. 38, no. 2, pp. 151–159, 2001. View at: Publisher Site  Google Scholar
 S. H. Hashemi Doulabi and A. Seifi, “Lower and upper bounds for locationarc routing problems with vehicle capacity constraints,” European Journal of Operational Research, vol. 224, no. 1, pp. 189–208, 2013. View at: Publisher Site  Google Scholar  MathSciNet
 L. Muyldermans, “Routing, districting and location for arc traversal problems,” Quarterly Journal of the Belgian, French and Italian Operations Research Societies, vol. 1, no. 2, pp. 169–172, 2003. View at: Publisher Site  Google Scholar
 D. Ambrosino and M. G. Scutellà, “Distribution network design: new problems and related models,” European Journal of Operational Research, vol. 165, no. 3, pp. 610–624, 2005. View at: Publisher Site  Google Scholar  MathSciNet
 K. Govindan, A. Jafarian, R. Khodaverdi, and K. Devika, “Twoechelon multiplevehicle locationrouting problem with time windows for optimization of sustainable supply chain network of perishable food,” International Journal of Production Economics, vol. 152, pp. 9–28, 2014. View at: Publisher Site  Google Scholar
 V.P. Nguyen, C. Prins, and C. Prodhon, “A multistart iterated local search with tabu list and path relinking for the twoechelon locationrouting problem,” Engineering Applications of Artificial Intelligence, vol. 25, no. 1, pp. 56–71, 2012. View at: Publisher Site  Google Scholar
 R. Liu, Z. Jiang, R. Y. K. Fung, F. Chen, and X. Liu, “Twophase heuristic algorithms for full truckloads multidepot capacitated vehicle routing problem in carrier collaboration,” Computers & Operations Research, vol. 37, no. 5, pp. 950–959, 2010. View at: Publisher Site  Google Scholar
 U. Derigs, M. Pullmann, U. Vogel, M. Oberscheider, M. Gronalt, and P. Hirsch, “Multilevel neighborhood search for solving full truckload routing problems arising in timber transportation,” Electronic Notes in Discrete Mathematics, vol. 39, pp. 281–288, 2012. View at: Publisher Site  Google Scholar
 J. Nossack and E. Pesch, “A truck scheduling problem arising in intermodal container transportation,” European Journal of Operational Research, vol. 230, no. 3, pp. 666–680, 2013. View at: Publisher Site  Google Scholar  MathSciNet
 E. M. Toro, J. F. Franco, M. G. Echeverri, and F. G. Guimarães, “A multiobjective model for the green capacitated locationrouting problem considering environmental impact,” Computers & Industrial Engineering, vol. 110, pp. 114–125, 2017. View at: Publisher Site  Google Scholar
 T. Bektaş and G. Laporte, “The pollutionrouting problem,” Transportation Research Part B: Methodological, vol. 45, no. 8, pp. 1232–1250, 2011. View at: Publisher Site  Google Scholar
 H. Moslemi and M. Zandieh, “Comparisons of some improving strategies on MOPSO for multiobjective (r, Q) inventory system,” Expert Systems with Applications, vol. 38, no. 10, pp. 12051–12057, 2011. View at: Publisher Site  Google Scholar
 K. Deb, A. Pratap, S. Agarwal, and T. Meyarivan, “A fast and elitist multiobjective genetic algorithm: NSGAII,” IEEE Transactions on Evolutionary Computation, vol. 6, no. 2, pp. 182–197, 2002. View at: Publisher Site  Google Scholar
 Z. Moravej, F. Adelnia, and F. Abbasi, “Optimal coordination of directional overcurrent relays using NSGAII,” Electric Power Systems Research, vol. 119, pp. 228–236, 2015. View at: Publisher Site  Google Scholar
 K. li, L. Lin, and S. Song, “Hybrid NSGAII algorithm on robust multiobjective inventory management problem,” Transactions of the Institute of Measurement and Control, vol. 37, no. 7, pp. 909–918, 2015. View at: Publisher Site  Google Scholar
 R. J. Labios, S. Kim, H. Song, and B. Lee, “Hybrid Multistarting GATabu Search Method for the Placement of BtB Converters for Korean Metropolitan Ring Grid,” Mathematical Problems in Engineering, vol. 2016, Article ID 1546753, 2016. View at: Publisher Site  Google Scholar
 W. Sukkerd and T. Wuttipornpun, “Hybrid genetic algorithm and tabu search for finite capacity material requirement planning system in flexible flow shop with assembly operations,” Computers & Industrial Engineering, vol. 97, pp. 157–169, 2016. View at: Publisher Site  Google Scholar
 P. W. Jones, “Rectifiable sets and the traveling salesman problem,” Inventiones Mathematicae, vol. 102, no. 1, pp. 1–15, 1990. View at: Publisher Site  Google Scholar  MathSciNet
 A. Messac, A. IsmailYahaya, and C. A. Mattson, “The normalized normal constraint method for generating the Pareto frontier,” Structural and Multidisciplinary Optimization, vol. 25, no. 2, pp. 86–98, 2003. View at: Publisher Site  Google Scholar  MathSciNet
 D. Kundu, K. Suresh, S. Ghosh, S. Das, B. K. Panigrahi, and S. Das, “Multiobjective optimization with artificial weed colonies,” Information Sciences, vol. 181, no. 12, pp. 2441–2454, 2011. View at: Publisher Site  Google Scholar  MathSciNet
 M. Barth and K. Boriboonsomsin, “Energy and emissions impacts of a freewaybased dynamic ecodriving system,” Transportation Research Part D: Transport and Environment, vol. 14, no. 6, pp. 400–410, 2009. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2018 Cheng Chen et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.