A Novel Multi-Objective Model for the Cold Chain Logistics Considering Multiple Effects

Abstract

This paper focuses on solving a problem of green location-routing with cold chain logistics (GLRPCCL). Considering the sustainable effects of the economy, environment, society, and cargos, we try to establish a multi-objective model to minimize the total cost, the full set of greenhouse gas (GHG) emissions, the average waiting time, and the total quality degradation. Several practical demands were considered: heterogeneous fleet (HF), time windows (TW), simultaneous pickup and delivery (SPD), and a feature of mixed transportation. To search the optimal Pareto front of such a nondeterministic polynomial hard problem, we proposed an optimization framework that combines three multi-objective evolutionary algorithms (MOEAs) and also developed two search mechanisms for a large composite neighborhood described by 16 operators. Extensive analysis was conducted to empirically assess the impacts of several problem parameters (i.e., distribution strategy, fleet composition, and depots’ time windows and costs) on Pareto solutions in terms of the performance indicators. Based on the experimental results, this provides several managerial insights for the sustainale logistics companies.

1. Introduction

With the development of urbanization and the change of client’s life style, the production and consumption of refrigeration-dependent food have changed, which promote the rapid development of the cold chain-based logistics [1]. Cold chain logistic(CCL), as a special type of transportation logistics, is developed to maintain the freshness of temperature-sensitive products by the thermal and refrigerated packaging methods and logistics plans [2,3]. Aiming at keeping supplies fresh and cut waste, CCL consumers more fossil fuels than ordinary logistics to maintain low-temperatures features during transportation, and greenhouse gas(GHG) emissions linearly related to freight fuel consumption account for 5.5% of global GHG emissions [4]. Meanwhile, clients’ satisfaction is another significant indicator not only concerning clients’ feelings about products and services received [5] but also concerning the long-term development of logistics enterprises. Hence, how to simultaneously consider economic, environmental, social, and cargos impacts is very important for the sustainable development of the CCL, which is the motivation of this paper to model a multi-objective model for the CCL aiming at considering the above four effects.

One of the significnt differences between ordinary logistics and CCL is the low-temperature feature that could maintain products’ freshness. The CCL has been successfully applied in frozen food logistics [6], pharmaceutical logistics [7], etc. In supply chains and logistics systems, there are two effective tools: vehicle-routing problem(VRP) and location-routing problem(LRP) [8], which construct two kinds of logistics for the cold chain: VRP-based CCL (VRPCCL) and LRP-based CCL (LRPCCL). The VRPCCL has been extensively researched by scholars, such as [4,9,10]. The VRPCCL didn’t consider a strategic-level problem (i.e., facility allocation problem, FAP), but the LRPCCL merged two decisions: FAP and VRPCCL, and has recently become one of the most investigated versions. Zheng et al, Wang et al and Ghomi et al. [11,12,13] defined a logistic costs-based formulation for the LRPCCL. Only the latter two considered lost sale costs. However, the above three only considered single-objective models. Wang et al. [14] proposed a bi-objective model for LRPCCL to minimize the total cost and distribution time. Leng et al. [3] developed a bi-objective model by defining the total costs and the total waiting time as two objectives concerning multiple practical constraints. The above two handled the carbon emissions and cargo quality decay as penalty functions in the total costs. In Leng et al. [2], a novel bi-objective model for the LRPCCL was proposed, which handles the cargo quality decay as the second objective. Leng et al. [15] developed several novel solution methods (i.e., decomposition-based hyper-heuristic approaches) to solve the bi-objective model proposed by Leng et al. [2].

To solve the limitation and fill gap of existing studies, this study will establish a novel multi-objective model based on sustainable effects of environment, economic, social, and cargo, it is named as GLRPCCL. Development concept of the model show as following: the first objective concerns the total costs consisting of the fixed costs of opening depots and renting vehicles, and the unfixed costs, including drivers’ salaries and the depreciation costs of vehicles (referred to as economic effects). The second is the full set of GHG emissions including CO${}_2$, CH${}_4$, and N${}_2$O emissions (referred to as environmental effects). The third aims to improve the clients’ and drivers’ satisfaction by minimizing the average waiting time of vehicles and clients (referred to as social effects). The last is the total cargos’ quality decay (referred to as cargos’ effects). This paper considers three MOEAs with an effective strategy differing from the above solution methods.

To enhance our model, this study will consider several practical constraints: SPD, hard TW, HF, and a feature of mixed transportation. The mixed transportation is defined by the types of shipments, including ordinary cargo (OC), refrigerated cargos (RC) at 3∼10 ${}^{\circ }$C, and frozen cargos (FC) at −4∼24 ${}^{\circ }$C. We assume that OC and RC could be simultaneously delivered to the clients, but FC must be separately delivered. As the solution methods, this paper proposed an effective framework that combines three well-known MOEAs, that is, non-dominated sorting genetic algorithm-II (NSGA-II) [16], strengthen Pareto evolutionary algorithm2 (SPEA2) [17], indicator-based evolutionary algorithm(IBEA) [18]. We also provided 16 operators, performing first (FI) and best-improvement (BI) search mechanisms, to form a large composite neighborhood for the proposed model. Extensive experiments are conducted to empirically assess the effects of several problem decisions, i.e., distribution strategy, fleet composition, and depots’ time windows and costs, on the performance indicators of Pareto results, which could provide several managerial views for the logistics enterprises.

As far as we know, no multi-objective model with more than two objectives for GLRPCCL has been proposed. The rest of this paper is organized: Section 2 provides the recent literature on green LRP (GLRP), green VRPCCL (GVRPCCL), and GLRPCCL. Section 3 gives out the proposed model for the GLRPCCL; Section 4 designs the MOEAs algorithm for the proposed model, including the solution chromosome, the MOEAs framework, and 16 neighborhood operators. Section 5 provides he detailed computational experiments and result analysis. Finally, the conclusions made of contributions, limitations, and the future study are outlined in Section 6.

2. Literature Review

From the perspectives of sustainable views of economy, environment, society, and cargos, this paper addressed a novel problem by defining multi-objective formulation, namely GLRPCCL, which has several features: location-routing decision, environmental impacts, cold chain, etc. Hence, the following sections provide reviews on GLRP, GVRPCCL, and GLRPCCL.

2.1. Review on GLRP

LRP is an extended version of VRP that takes into account FAP, which involves locating depots and routing a fleet of vehicle to serve a given set of clients, the purpose of which is to minimize the total cost of location and routing [19]. Among the LRP variants in the literature, the most studied are LRP with SPD [8], LRP with HF [20],LRP with TW [21], multi-level LRP [22], and LRP considering multiple constrains [23]. However, the above only considered the economic impacts of logistics.

As the name suggests, the GLRP deals with LRP considering environmental effects. Dukkancj et al. [19] stated that the GLRP is constructed by the classical LRP and pollution-routing problem which is also viewed as green VRP. As far as we know, many scholars have studied GLRP. The first one is by Mohammadj et al. [24], who described a bi-objective stochastic green hub location routing problem with simultaneous pickups and delivers, aiming at minimizing both the logistics cost and environmental impacts.

Table 1. Features of the GLRP without considering the cold chain.

Type	Reference	Emission Model	Solution Method	Features	Method
SO	Koc et al. (2016) [32]	CMEM	Exact & heuristic	HF	Penalty
	Qazvini et al. (2016) [29]	Factor	Exact	SPD, TW	Constraint/Penalty
	Nakhjirkan and Rafiei (2017) [30]	Factor	Exact & heuristic	Inventory	Penalty
	Wang et al. (2017) [33]	Given	–	Inventory	Penalty
	Leng et al. (2018) [34]	CMEM	Heuristic	TW	Penalty
	Benotmane et al. (2019) [27]	Factor	Heuristic	Two echelons	Penalty
	Dukkanci et al. (2019) [19]	CMEM	Exact & heuristic	TW, no vehicles cost	Penalty
	Koc (2019) [23]	CMEM	Heuristic	TW	Penalty
	Zhang et al. (2019) [35]	CMEM	Heuristic	TW, HF, SPD	Penalty
	Zhang et al. (2020) [36]	CMEM	Heuristic	Time-dependent, TW	Penalty
	Pitakaso et al. (2020) [28]	Given	Heuristic	Two echelons	Main objective
	Zhou et al. (2020) [25]	Factor	Exact	Robust and stochastic	Penalty
BO	Mohammadi et al. (2013) [24]	Factor	Heuristic	TW, SPD	Distance as a green object
	Govindan et al. (2014) [37]	Factor	Heuristic	Two echelons, TW	As a green objective
	Validi et al. (2014) [38]	Given	Heuristic	HF	As a green objective
	Tang et al. (2016) [39]	Factor	Heuristic	Inventory	As a green objective
	Tricoire and Parragh (2017) [40]	Factor	Exact & heuristic	HF	As a green objective
	Toro et al. (2017) [41]	Macro	Exact	Net present value	As a green objective
	Wang and Li (2017) [42]	Factor	Heuristic	TW, HF, SPD	Penalty
	Qian et al. (2018) [43]	Factor	Heuristic	Hyperheuristic	As a green objective
	Chen et al. (2018) [44]	Given	Heuristic	Full truckloads	As a green objective
	Faraji and Afshar-Nadjafi (2018) [45]	Given	Heuristic	TW, HF, multiple periods	As a green objective
	Leng et al. (2019a) [46]	CMEM	Heuristic	TW, HF, SPD	Penalty
	Govindan et al. (2020) [31]	Given	–	Reverse logistics, HF	Penalty
TO	Rabbani et al. (2017) [47]	Given	Exact & heuristic	Refueling stations	As a green objective
	Leng et al. (2019b) [48]	CMEM	Heuristic	TW, HF, SPD	Penalty
	Shen et al. (2019) [49]	Factor	Heuristic	Emergency chain, fuzzy	As a green objective
MO	Our study	CMEM	Heuristic	Cold chain, TW, HF, SPD	As a green objective

SO/BO/TO/MO is a model using the single-/bi-/tri-/multi objectives; Factor and Macro models are classified by Demir et al. (2014) [50], CMEM is a micro model.

is an extended one that reviews the papers on GLRP from the used emission models, solution method, feature, and method handling environmental effects, following the methods proposed by Dukkancj et al. [19] and Zhou et al. [25]. Looking at Table 1, following can be detected:

(1)

From the perspectives of the number of objectives defined for the GLRP, there are single- objective, bi-objective, and triple-objective models, but multi-objective models with over three objectives have not been studied. The single and bi-objective model has been widely researched.

(2)

In views of emission model, most of these papers applied factor model to estimate the amount of fuel consumption and carbon emissions, like the fuel consumption rate proposed by Xiao et al. [26], models using fuel efficiency or emission parameter (see “Given” in Table 1). As an efficient micro view, CMEM has been widely used in the GLRP, which is easily applicable and is capable to accurately estimate fuel consumption and carbon emission. Only Benotmane et al. [27] used a macro version to calculate fuel consumption.

(3)

For handling the fuel consumption and carbon emission, several methods have been used as (a) a part of logistics costs, namely penalty function, (b) the main objective like Pitakaso et al. [28], (c)  constraints like Qazvini et al. [29], and (d) an objective in the multi-objective model.

(4)

Most of them developed heuristics to tackle the proposed GLRP model, and a few papers also applied exact approaches. Moreover, there are two papers by Nakhjikan et al. [30] and Govindan et al. [31] that studied the GLRP in terms of theory and applied practical cases to verify their models.

However, the above models only have three objectives, and no paper has investigated the models with over three objectives, such as the total cost, traveled distance/time, GHG emissions, clients’ satisfaction, cargo quality decay.

2.2. Review on GVRPCCL

With the increasing improvement of people’s living standards and the growing demand for high-quality and fresh food, the CCLs industry has developed rapidly [9]. Ma et al. [51] studied the CCL in terms of Industry 4.0 and provided theoretical proof for their models and carried out numerical studies to illustrate the theoretical results. Al et al. [52] defined an integrated mixed integer optimization model that considers inventory allocation problems, VRP, and CCL. Wei et al.  [53] defined a model on the transportation cost for the inventory routing problem combing CCL. However, the environmental effect is not considered by the above papers. As stated by Accorsi et al. [54], as much as 15% the worlds’ total energy already fuels cold-chain infrastructures with 40% of food logistics needing refrigeration, both the growth of food demand and the CCLs’ expansion will therefore hugely increase energy consumption and associated GHG emissions. GVRPCCL concerns economic, environmental, and cargos’ effects.One of the greatest challenges posed by today’s CCLs is to provide high-quality food with minimizing GHG emissions and total logistics costs. In recent years, the GVRPCCL has attracted extensive academic attention, and also provided optimization support for real-world applications, such as fruit-and-vegetable [55,56].

Table 2. Features of the papers on the GVRPCCL.

Type	Reference	Emission Model	Solution Method	Method 1	Method 2	Feature
SO	Ma and Wu (2015) [57]	Given	Heuristic	Penalty	Penalty	Stochastic demands, TW
	Hariga et al. (2017) [58]	Factor	Exact & Heuristic	NA	Penalty	Multi-stage, carbon tax regulation
	Wang et al. (2017) [59]	Factor	Heuristic	Penalty	Penalty	Carbon tax, TW
	Hsiao et al. (2017) [60]	Given	Heuristic	Penalty	Penalty	TW, shelf life
	Hsiao et al. (2018) [55]	Given	Heuristic	Penalty	Penalty	TW, shelf life
	Stellingwerf et al. (2018) [61]	CMEM	Exact	NA	Penalty	Stops, temperature-VRP
	Zhang et al. (2019) [62]	Factor	Heuristic	Penalty	Penalty	Soft TW
	Chen et al. (2019) [9]	CMEM	Heuristic	Penalty	Penalty	HF, TW
	Wei et al. (2019) [53]	Given	Exact & Heuristic	NA	Penalty	Time period
	Qin et al. (2019) [5]	Factor	Heuristic	Penalty	Penalty	TW, Client satisfaction
	Chen et al. (2019) [63]	Factor	Exact & Heuristic	NA	Penalty	TW, multi-compartment
	Li et al. (2019) [4]	Factor	Heuristic	Penalty	Penalty	Soft TW
	Li et al. (2020) [64]	Given	Heuristic	Penalty	Penalty	HF, TW
	Wang and Wen (2020) [14]	CMEM	Heuristic	Penalty	Penalty	Two-echelon, HF, TW
	Wang et al. (2020) [10]	Factor	Heuristic	Penalty	Penalty	Soft TW, risk factor
	Qi and Hu (2020) [65]	Factor	Heuristic	Penalty	Penalty	Emergency CCL, Minimum loss
BO	Ma et al. (2018) [51]	Given	Exact	An objective	Constraint	Shelf life
MO	This study	CMEM	Heuristic	An objective	Green objective	Cold chain, TW, HF, SPD

provides the papers on GVRPCCL in terms of the number of objectives, emission model, solution method, method1 handling cargos’ quality decay, and method2 dealing with GHG emissions.

    Looking at Table 2, the GVRPCCL has been widely researched, and we can obtain the following summaries. (1) Only one paper considered the bi-objective model for the GVRPCCL. Ma et al. [51] defined the model by minimizing transport cost and quality degradation. (2) like the conclusion in Section 2.1, most of the papers applied factor models and special parameters-models (i.e., Given) to estimate the fuel consumption and carbon emissions, and only three papers used CMEM. (3) Most papers applied penalty function to deal with the cargo quality decay, and four papers didn’t consider the cargo quality decay, and only one viewed it as an objective. (4) All papers applied penalty function to handle the fuel consumption/carbon emissions as a part of transport costs, but Ma et al. [51] applied a constraint to restrict the CO${}_2$ emission regulated by the government. (5) Among the papers handling cargo quality decay, two methods have been developed: (a) shelf-life applying piecewise function [55,60]; (b) the variable function of refrigerated goods quality [33].

However, the multi-objective model with three or more objectives for the GVRPCCL is not studied. Meanwhile, LRP, as an extensive version of VRP, should be considered in the CCL, since location and route problems are equally significant for CCLs, it is indispensable to combine the depots’ location selection with route optimization for considering the comprehensive impacts bringing by the entire logistics system [11].

2.3. Review on GVRPCCL

The GLRPCCL is an extensive version of GVRPCCL by considering FAP, the reason could be that: the location of the facility could not only improve the economic, environmental, and social impacts but also could greatly reduce cargos’ quality decay. To our knowledge, only a few studies have considered both LRP and CCL. To our knowledge, only a few studies have considered both LRP and CCL. This lack of studies is seen, for instance, in the field of urban studies as transportation is currently dealing with urban planning and consequently urban renewal and regeneration [66]. If we take into account one of the most innovative way of provide the so-called right to the city [67], i.e., the Barcelona’s Poblenou neighborhood superblocks, it is seen that the most recent studies [68,69,70,71] tackling the former industrial heart of the Catalan capital did not take into account the issues this paper shows slightly referred on the impacts on transportation system introduced by superblock unit in Poblenou.

Shen et al. [49] developed a bi-objective model for the LRPCCL, where the first objective optimizes the total cost concerning transportation, refrigeration, cargo damage, and punish, and the second is to minimize the no-load cost of vehicles. Zheng et al. [11] defined an optimal location-routing model that takes into account the temperature changes among sensitive products in cold chain distribution centers and aims to minimize the sum of fixed costs, variable costs, and transportation costs. Zheng  et al. [72] provided a multi-objective LRP model for the two-echelon cold chain logistics network. The first objective aims to minimize the location cost of depots, and the second is used to minimize the transportation cost. Wang et al. [73] split the GLRPCCL into two parts: the location of depots and the optimization of the transporting route. In the location of depots, the objective is to minimize the total related costs, including installation costs, turnover costs, and transportation costs. The second level is used to optimize the costs comprising vehicle cost, damage cost, and penalty cost for violation of time window. Suraraksa et al. [74] also applied two phases to model their problem. The first phase is to determine the location of fresh cargos distribution and demand distribution. The second phase is to design the delivery route and the number of vehicles rented to minimize the total distance travelled, while transporting the cargos from the wholesaler through the depots to the retailer. However, the above didn’t consider environmental effects. And the latter two separately optimized the FAP and VRPCCL, which may cause suboptimal solutions, like LRP [75]. The following are the latest papers on the GLRPCCL.

Wang et al. [12] developed a single objective model for GLRPCCL, with the lowest total cost as the objective function, which includes carbon emission costs, fixed costs of depots, transportation costs related to refrigerated trucks, refrigeration costs, fine costs, and damage costs.

In Ren et al. [76], a mathematical formulation concerning the minimum distribution cost of fresh food was developed and considered resource/information sharing and soft time window constraints. Their model can meet the requirements of logistics distribution, including high timeliness, low cost, greenness, and resource/information sharing. The objective is the total logistics cost, including the costs of scheduling vehicles, transportation, cargo loss, carbon emissions, and penalty.

Van et al. [77] developed a two-stage distribution location-routing model with the minimum total logistics costs as the objective, considering varying capacity of vehicles according to different delivery stages. Their objective comprises fuel consumption cost of vehicles per kilometer, driver’s cost, damage cost, refrigeration cost (i.e., energy consumption), penalty cost, and transfer cost.

The above three papers defined single-objective models for the GLRPCCL. The following defined bi-objective models for the GLRPCCL.

Wang et al. [10] designed a bi-objective model for the GLRPCCL. In their model, the first objective is used to minimize of total cost comprising of fuel consumption cost, quality damage cost, transport cost, refrigeration cost of refrigerated vehicles, penalty cost of the time window, and opening cost of depots, while the second is the shortest vehicle distribution time.

In the model proposed by Leng et al. [3], the first objective comprises the fixed costs of depots and vehicles, vehicle renting cost, driver salaries, fuel consumption cost, carbon emission costs, and damage costs of cargos that need to be refrigerated or frozen. The second objective consists of the waiting time of clients and vehicles to promote client satisfaction and the efficiency of the CCLs. Leng et al. [15] also developed a bi-objective model for the GLRPCCL, where the first objective is to minimize the total cost of logistics, including the costs of operating depots, renting fleet, fuel consumption, and carbon emissions, and the second objective is to reduce cargo damage, which can improve client satisfaction. The formulation was also solved by decomposition-based hyperheuristic approaches [2].

From the perspectives of the above three papers, several conclusions are drawn: (1) only bi-objective models are defined; (2) the penalty function is used to view the fuel consumption and carbon emissions as a part of the total costs; (3) the constraints SPD, TW, HF, and a feature of the mixed transportation are considered; (4) Only carbon emissions are concerned. However, this paper defines the model from a novel perspective, a multi-objective model is developed for the GLRPCCL by defining four objectives. The first objective is to optimize the economic effects concerning the fixed costs of opening depots and renting fleet and the unfixed costs including drivers’ salaries and the depreciation costs of vehicles. The second is to optimize the environmental effects, which is made of the full set of GHG emissions including CO${}_2$, CH${}_4$, and N${}_2$O emissions. The third is to optimize the social effects characterized by the average waiting time of drivers and clients. The last one is to optimize the cargo’s effect that is affected by the cargos’ quality decay. To the best of our knowledge, the proposed multi-objective formulation for the GLRPCCL has not been researched by others thus far.

3. Formal Problem Description and Mathematical Model

The GLRPCCL is defined on a complete directed graph G=(V, E) where $V.=\{1;;N.;N.+1;;N.+M.\}$ denotes the set of nodesrepresenting both the set of clients ($C=1,\dots ,N$) and the set of candidate depots ($D=\{N+1,...,N+M\}$) to be selected and $E=\left\{\right(i,j):i,j\in V,i\ne j\}\backslash \left\{\right(i,j):i,j\in D\}$ is the set of arcs. For a client $i\in C$, the pickup and delivery demand is $p_i$ and $d_i$, and the time window is [$e_i$, $l_i$]. Each client has one identical type of cargos (OC, RC, and FC) and service time $s{t}_i$ depending on the pickup and delivery demand. For a depot $j\in D$, the fixed costs of operating it, capacity, and closing time window are, respectively, $C{D}_j$, $F{D}_j$, and $TW{D}_j$. A heterogeneous fleet $H=\{h1,h2,h3\}$is used to serve the clients, and each type of vehicle has a fixed cost of renting vehicle $F{V}_h$ and a fixed capacity $C{V}_h$$(h\in H)$. The drivers’ salary for each type of vehicle is denoted as $D{S}_h$. The traveling speed and distance over $(i,j)\in E$ are pre-known and fixed. The goal of the presented model is to address the sustainable development for the CCLs in views of improving the economic, environmental, and social effects, and also reducing the quality degradation of RC and FC.

Before providing our model, several assumptions are given: (1) Each client must be served only once by a depot and vehicle; (2) No overflows on the capacity of a vehicle and depot are allowed; (3) The types of pickup and delivery demand of each client is the same; (4) The information of each arc including speed and distance is known and static; (5) If the vehicle arrives early, it must wait until the time window of each client is open, and must return to the original depot which it departs from before the closing time window; (6) The OC and RC can be stored in the same vehicle, but the FC must be separated.

Moreover, the proposed multi-objective model relies on the following three decision variables: (1)$x_{ijh}=1$ if the vehicle $h\in H$ continuously visited clients i and $j\in C$ and to 0 otherwise; (2) $y_i=1$ if a depot $i\in D$ is open and to 0 otherwise; (3) $z_{ij}=1$ if the client $i\in C$ is served by the depot $j\in D$ and to 0 otherwise.

3.1. Objectives Development

In the proposed multi-objective model for the GLRPCCL, four objectives are defined to concern economic, environmental, and social impacts, as well as to reduce the quality degradation of cargos in terms of RC and FC. The following four aspects are described as following.

(1)

Economic aspect (the first objective)

Economic effects refer to the sustainable development concerning the long-term economic interests of the logistics enterprise, namely the total logistics costs, which includes three parts: the fixed costs of operating the selected depots, the fixed costs of renting fleet, and the variable costs of hiring drivers and vehicle maintenance. Hence, the economic indicator is denoted as

$$TC=\sum _{i\in D}F{D}_i{y}_i+\sum _{i\in D}\sum _{j\in C}\sum _{h\in H}F{V}_h{x}_{ijh+}\sum _{i\in C}\sum _{j\in D}\sum _{h\in H}D{S}_h{AT}_{jh}{x}_{ijh}$$

(1)

where the unit of $D{S}_h$ is Yuan/min; $A{T}_{jh}$ (in seconds) is the moment that the type$h\in H$ of vehicle arrives at the node $j\in V$.

(2)

Environmental aspect(the second objective)

The transportation industry has a major impact on the planet, because of the large amount of fuel used in its daily operations and the environmental consequences and greenhouse effects resulted by consuming fossil fuel and exhausting gas [41]. Moreover, in cold chain logistics, vehicles consume fuel and emit GHG during transportation and refrigeration equipment. Hence, the  environmental impacts in this paper refer to the full set of GHG emissions including carbon dioxide, methane, and nitrous oxide. The fuel consumption during the transportation is

$$TFC(t_1,t_2,s\left(t\right),h)=\lambda (k_h{N}_h{T}_h+P_h^{rc}/\eta )\times (t_2-t_1)+\gamma \alpha G{\int }_{t_1}^{t_2}s\left(t\right)dt+\gamma {\beta }_h{\int }_{t_1}^{t_3}{s\left(t\right)}^{3}dt$$

(2)

where $TFC(t_1,t_2,s\left(t\right),h)$ represents that the total fuel consumption (in liter) of a type of vehicle $h\in H$ traveling at $s\left(t\right)$ from $t_1$ to $t_2$; G(in kg) is the sum of vehicle’s net weight $\left(N{W}_h\right)$ and the weight of load; $P_h^{rc}$is the engine power requirement related to engine wear and the operation of vehicle accessories such as air conditioners and refrigeration compressors used in CCLs; $\lambda =\phi /\kappa \psi$, $\gamma =1/(1000\times n_{tf\eta })$, $\alpha =a+gsin\theta +g{C}_{r}cos\theta$, and ${\beta }_h=0.5C_h^d\rho A_h$; and the detailed parameters in Equation (2) could see Tables S1 and S2 in the Supplementary Materials.

Considering the fourth assumption, in our instances, the speed and distance are no-dynamic. hence, the fuel consumption $TF{C}_h\left(L\right)$ of vehicle type h over an arc $(i,j)\in E$ with a distance $d_{ij}$ (in meter) traveling at constant speed $s_{ij}$ (in m/s) is rewritten as

$${TFC}_{ijh}=\lambda ({(k}_h{N}_h{T}_h+P_h^{rc}\left(index\right)/\eta )\times d_{ij}/s_{ij}+\gamma \alpha G{d}_{ij}+\gamma {\beta }_h{d}_{ij}{s}_{ij}^2)\times x_{ijh}$$

(3)

where $P_h^{rc}$(index) is the extra power need to maintain the freshness of RC and FC, and it equals to 0 if only no RC and FC are loaded. For example, for a type of vehicle $h_1$, if the types of loaded cargos contain RC, $P_h^{rc}$ (equals to 5 kW/s, and if the type of loaded cargos is FC, $P_h^{rc}$ (equals to 8 kW/s. However, Equation (3) is not suitable for the case that the speed of vehicle equal to 0, since if the types of cargos contain RC or FC and the vehicle must wait for opening the time windows and serving the client, then the engine must not shut down. Hence, extra fuel consumption follows Equation (2).

$${EFC}_{ijh}=\lambda (k_h{N}_h{T}_h+P_h^{pc}\left(index\right)/\eta )\times (max(e_j-{AT}_{jh},0)+{st}_j)\times tc\times x_{ijh}$$

(4)

where $t_c$ is the state of the vehicle, in other words, $t_c$=0 if the type of cargo is OC and to 1 otherwise. Hence, the total fuel consumption can be:

$$LFC=\sum _{i\in V}\sum _{j\in V}\sum _{h\in H}({EFC}_{ijh}+{TFC}_{ijh})$$

(5)

According to the “Greenhouse gas emissions accounting method for land transportation enterprises” (National Development and Reform Commission of the People’s Republic of China [78]), GHG emissions are composed of methane, carbon dioxide, and nitrous oxide, and GHG emissions are estimated as follows:

$$E=E_{{CO}_2}+E_{{CH}_4}+E_{N_{2}O}$$

(6)

CO${}_2$ emissions are directly proportional to fuel consumption. In other words, 1 liter of diesel can produce 2.32 kg of carbon dioxide. For the amount of CH${}_4$ and N${}_2$O, we follow the methods used by [4]:

$$E_{{CH}_4}=\sum _{i\in V}\sum _{i\in V}{d}_{ij}\times {EF}_{{CH}_4}\times {GWP}_{{CH}_4}\times {10}^{-9}$$

(7)

$$E_{N_{2}O}=\sum _{i\in V}\sum _{i\in V}{d}_{ij}\times {EF}_{N_{2}O}\times {GWP}_{N_{2}O}\times {10}^{-9}$$

(8)

Hence, the environmental indicator is represented by total emissions:

$$TE=2.32\times LFC+E_{{CH}_4}+E_{N_{2}O}$$

(9)

(3)

Social aspect(the third objective)

The social indicator is extremely significant for the logistics enterprises’ development in long-term cooperation with clients and drivers by earning their loyalty and favor. This paper uses average waiting time (AWT) to represent the social indicator:

$$minAWT=\sum _{i\in V}\sum _{j\in C}\sum _{h\in H}\{(max(e_j-{AT}_{jh},0)+\alpha (max(e_j,{AT}_{jh})-e_j)\}\times x_{ijh}/N$$

(10)

The first part max($e_j-A{T}_{jh}$, 0) represent that the waiting time of vehicle $h\in H$ at client $j\in V$; the second part $max(e_j,A{T}_{jh})-e_j$ is the waiting time of client $j\in V$ before being served; the parameter $\alpha$ is the degree of importance for the waiting time of the client, in this paper, we set $\alpha =2$. Moreover, we set that $l_j$ equals to $e_j+s{t}_j$.

(4)

Quality decay of cargos (the fourth objective)

The purpose of CCL is to improve the quality of RC and FC as much as possible, but the RC (e.g., strawberries, milk, bayberry) and FC (e.g., seafood) can easily deteriorate in the CCL, so this kind of cargos needs to be maintained in a proper low-temperature environment. Perishable goods will gradually decline in quality or lose value over time. When product quality drops to a certain extent, degradation in quality will occur. The quality degradation of cargos consists of two parts: the degradation that accumulates over time during transport and the degradation that occurs when the door is opened during unloading. We used the variable function of refrigerated goods quality proposed by [79] to estimate the quality decay: $D\left(t\right)=D\theta e-\delta t$.

The cargoes damage CQD1 resulted by the refrigerated vehicles in the travel process and the cargoes loss CQD2 during the vehicles serve the clients can be defined as:

$${CQD1}_{ijh}=Q_{ijh}^{*}\ast (1-e^{-{\delta }_{1,{\mathrm{state}}^{({AT}_{jh}-max{e}_i,{AT}_{ih}+max{e}_j-{AT}_{jh},0-{st}_i)}}})$$

(11)

$${CQD2}_{ijh}=(Q_{ijh}^{*}-q_j*r)*(1-e^{{-\delta }_{2,{state}^{ST}j}})$$

(12)

where $Q_{ijh}^{*}$ is the weight of RC/FC cargos loaded by the vehicle $h\in H$ traveling over an arc $(i,j)\in E$, and it equals to 0 if no RC/FC is loaded; ${\delta }_{1,state}$ and ${\delta }_{2,state}$ are the spoilage rates of each type of cargos, and $state\in \{1,2,3\}$; let r equals to 0 if no RC and FC is load and to 1 otherwise. Hence, the total quality degradation of cargos as follows:

$$CQD=\sum _{i\in V}\sum _{j\in V}\sum _{h\in H}{CQD1}_{ijh}{x}_{ijh}+\sum _{i\in V}\sum _{j\in C}\sum _{h\in H}{CQD2}_{ijh}{x}_{ijh}$$

(13)

3.2. Constraint

This section provides the necessary constraints according to the above assumptions, described as follows.

$$\sum _{h\in H}\sum _{i\in V}{x}_{ijh}=1,j\in C$$

(14)

$$\sum _{h\in H}\sum _{i\in V}{x}_{jih}=\sum _{h\in H}\sum _{j\in V}{x}_{ijh},i\in C$$

(15)

Particularly, constraint (14) ensures that each client must be served exactly once. Constraint (15) makes sure that entering and leaving edges to each client are equal.

$$max\{\sum _{i\in C}{d}_i{z}_{ij},\sum _{i\in C}{p}_i{z}_{ij}\}\le {CD}_j{y}_j,j\in D$$

(16)

$$\sum _{h\in H}\sum _{i\in C}{Q}_{jih}=\sum _{i\in C}{d}_i{z}_{ij},j\in D$$

(17)

$$\sum _{h\in H}\sum _{i\in C}{Q}_{ijh}=\sum _{i\in C}{p}_i{z}_{ij},j\in D$$

(18)

$$0\le Q_{ijh}\le C{V}_h{x}_{ijh},i\in V,j\in V,h\in H$$

(19)

$$Q_{ijh}\ge (d_i-p_i)x_{ijh},i\in C,j\in V,h\in H$$

(20)

$$Q_{ijh}\le ({CV}_h-d_j+p_j)x_{ijh},i\in V,j\in C,h\in H$$

(21)

$$\sum _{i\in D}\sum _{j\in C}{Q}_{ijh}=\sum _{i\in C}\sum _{j\in V}{d}_i{x}_{ijh},h\in H$$

(22)

$$\sum _{i\in C}\sum _{j\in D}{Q}_{ijh}=\sum _{i\in C}\sum _{j\in V}{p}_i{x}_{ijh},h\in H$$

(23)

$$\sum _{i\in V}\sum _{h\in H}(Q_{ijh}-d_j)x_{ijh}=\sum _{i\in V}\sum _{h\in H}(Q_{jih}-p_j)x_{jih},j\in C$$

(24)

Constraints (16)–(18) are the limitations on the load of the depots, and constraints (19)–(24) forbid that the load of the vehicle exceeds the vehicle capacity. In detail, constraint (16) is the limitation on the load of the depot. Constraint (17) guarantees that a total load of each depot is equal to the total delivery demand of clients assigned to it before starting to serve clients. Constraint (18) guarantees that a total load of each depot is equal to the total pickup demand of clients assigned to it when the vehicles return to it. Constraints (19)–(21) are bounding constraints for load variables. Constraint (22) imposes that a load of each vehicle must equal the total delivery demand before serving clients. Constraint (23) ensures that a load of each vehicle must equal the total pickup demand after returning to the depot. Constraint (23) is the equilibrium constraint for the load of each vehicle traveled over each arc. $Q_{ijh}$ is the load of a type of vehicle $h\in H$ traveled over an edge $(i,j)\in E$.

$${AT}_{jh}=(max\{{AT}_{ih},e_i\}+{st}_i+d_{ij}/s_{ij})x_{ijh},i\in V,j\in V,h\in H$$

(25)

$${AT}_{ij}{x}_{ijh}\le {AT}_{jh}{x}_{ijh}\le l_j,i\in V,j\in V,h\in H$$

(26)

$${AT}_{jh}\le {TWD}_j,j\in D,h\in H$$

(27)

Constraint (25) is the temporal coherence. Constraint (26) makes sure that the vehicle must reach each client before its closing time. Constraint (27) is the limitation on that each vehicle must return to its departed depot before its closing time window.

$$\sum _{h\in H}{x}_{ijh}\le z_{ij},i\in C,j\in D$$

(28)

$$\sum _{h\in H}{x}_{ijh}\le z_{ji},i\in C,j\in D$$

(29)

$$\sum _{h\in H}{x}_{ijh}+z_{ik}+\sum _{m\in M,m\ne k}{z}_{jm}\le 2,i,j\in C,k\in D$$

(30)

Constraints (28)–(30) forbid illegal routes, i.e., routes that do not start and end at the same depot.

$$\sum _{j\in D}{z}_{ij}=1,i\in C$$

(31)

$$x_{ijh}+\sum _{k\in V}\sum _{p\in H,p\ne h}{x}_{jkp}\le 1,i\in V,j\in C,h\in H$$

(32)

Constraints (30) and (31) ensure that each client is assigned to only one depot and vehicle, respectively.

$$z_{ij}\le y_j,i\in C,j\in D$$

(33)

$$\sum _{i\in C}{z}_{ij}\ge y_j,j\in D$$

(34)

Constraint (33) makes sure that the unselected depot cannot serve the client. Constraint (34) guarantees that the opened depot must serve at least one client.

$${ct}_i\ast {ct}_j\ast x_{ijh}\notin \{3,6\},i\in C,j\in C,h\in H$$

(35)

Constraint (35) limits on the type of cargos loaded by each vehicle. $c{t}_i$ is the type of cargos of client $i\in C$, and $cti\in 1,2,3$ where 1, 2, and 3 represents the OC, RC, and FC, respectively. Hence, the mixed type of cargos loaded by a vehicle must not belong to RC&FC and OC&FC.

4. Solution Method

In the single-objective problems, the purpose is to search one and only one optimal solution. While the multi-objective problem is to obtain a set of Pareto solutions. This section provides the solution method to obtain the Pareto solutions for solving the proposed multi-objective problem. The following describes the solution representation, 16 neighborhood operators, and an effective algorithm framework of MOEAs.

4.1. Solution Representation

One significant decision to develop an algorithm for combinatorial optimization problems is to decide how to represent the solution in an effective way and associate it with the search space [24]. The chromosome used represents a complete solution, i.e., a set of routes. Each routing is stored in the cell array (called as route cell) which is also stored in another cell array (called as depot cell) that has a fixed length equaling the number of depots M. The depot cells contain the information about decision on being selected whether or not and the route cells. If there is an unopened depot, the information is set to ∅ and no routes are assigned. Moreover, we also apply the objective cell which is in line with the objective values of the proposed problem which allowing a fast evaluation of the change of objective values.

Figure 1 is a simple chromosome representation of a solution, which has a set of 15 clients $C\in \{1,2,\dots ,14,15\}$ and a set of 5 possible depots $D\in \{16,17,18,19\}$. The solution in the figure selects {16, 17, 19} depots to serve the clients in the given order by four routes (vehicles).

The applied solution representation, together with the following 16 neighborhood operators, allows obtaining feasible results, which could avoid to using repair methods for restoring feasibility. Hence, the proposed representation and operators can speed up the search over the solution spaces.

4.2. Neighborhood Operator

Aiming at effectively solving the proposed problem, this paper developed a large composite neighborhood formed by 16 operators, which is grouped into three modules: dominated pool (DP), non-dominated pool (NDP), and mutational pool (MP). The purpose of DP is to obtain the children solutions dominating parent solutions, which could speed up convergence in the early stages. The goal of NDP is to achieve the children solutions that can’t be dominated by parent solutions, which could help in obtaining a well-distributed Pareto front in the later period of the algorithm. However, NDP may cause a stagnant state for the search process if the DP is not used. Although MP is usually not sufficient to obtain competitive solutions, it can provide randomization when searching global optimal solution. Note that the DP and NDP are based on local search procedures.

The set of MP consists of 6 mutational operators: “add”, “delete”, “crossover”, “insert”, “swap”, and “decompose”. The first mutational operator was proposed by [80], named “add". It could avoid premature convergence with few depots. The mechanism is to diversify the selected depots by randomly choosing a new one, before reassigning between 1 and 2/3 of the routes to it. The second mutational operator randomly deletes one opened depot and reassigns the routes to other opened depots. The third mutational operator is “crossover” which exchanges one route of an opened depot with another route of the other opened depots until a feasible solution is generated. The “insert" operator is executed by inserting a client into another route. The “swap" operator swaps two clients from different routes. Finally, the “decompose” operator split a solution into sub solutions, which could avoid a long tracing of the route of vehicles with larger load. However, the use of MP is not to search for better results but to provide randomness.

The set of DP and NDP consists of five local search operators: “2-opt”, “swap+”, “insert+”, “fragment-based swap”, and “fragment -based insert”, which is used to search dominated solutions or non-dominated results. The “swap+”/“insert+” operators are similar to the mutational operators “swap”/“insert”, while “fragment-based swap” and “fragment-based insert” are used to dealing with a fragment of two or three client locations rather than one client. The operator “2-opt” removes the two arcs from two different routes before reconnecting the new four fragments created. Figure 2 shows the schematic diagrams of several operators.

Moreover, in this paper, tow search strategies are proposed for DP and NDP operators: first (FI) and best improvement (BI). FI returns with a dominated solution for DP and a non-dominated solution for NDP once the solution is searched. While BI returns with a Pareto solution dominating parent solution for DP and a set of Pareto results that can’t be dominated by the parent solution for NDP until to extra better solutions can be found. However, the decisive issue is to design the stop criteria for providing a fair comparison, described in Section 4.3.

4.3. Optimization Framework

The proposed MOEA algorithm starts by initializing the population (P), calculating fitness (FP), and the archive population (AP) only including Pareto solutions, which is then returned (Step 1). After that, the main loop is executed and stops when the maximum number of iterations $T_{max}$ is reached.

In each iteration of the main loop, tournament selection is used to select individuals to construct the mating pool according to the evaluation function in NSGA-II, SPEA2, and IBEA (Step 4). For each individual in the mating pool, a three-phase evolutionary process is developed to modify it. First, the proposed algorithm randomly selects a mutational operator from MP to perturb the individual (Steps 6 and 7). Second, a local search procedure is randomly selected from DP to search for better solutions that dominated the individual generated in the first phase (Steps 8 and 9). The third phaser randomly select a local search procedure from NDP to achieve non-dominated results (Steps 10 and 11). The reason for this strategy is that DP can speed up the convergence by obtaining better solutions and NDP can achieve many non-dominated solutions, which could search the well-distributed solutions located in the real Pareto front (RPF). Steps 14 and 15 are used to update the population which will be evolved in the next iteration, and the environmental selections of NSGA-II, SPEA2, and IBEA are used. The goal of Steps 16–18 is to update the AP made of $10\times Np$ individuals by the environmental selection used in NSGA-II. The reason returning the AP is that the proposed problem is NP-hard, the population AP contains the best Pareto solutions so far in each independent run. $p_m$ is the mutational probability.

Algorithm 1 MOEA framework for the GLRPCCL

Input:$N_p$, $T_{max}$, $p_m$, and k in IBEA

Output: Archive population ($AP$)

//Initialization

1: T←0; Generate(P); Calculate fitness ($FP$); $AP$←P

//Main loop

2: Repeat

3:       T←$T+1$; $OP$=$\Phi$

4:       Tournament selection (P, $FP$) to construct mating set ($MS$)

5:          while i<$N_p$ do

//Mutation

$CS$←$MS\left(i\right)$

6:          if rand ≤ $p_m$ then

Randomly select an operator ($op$) in $MP$

Apply $op$ to modify $MS\left(i\right)$ and obtain $CS$

7:          end if

//Local search

8 :         Randomly select an operator ($op$) in $DP$

9 :         Apply op to improve $CS$ and obtain $C{S}_2$

10:         Randomly select an operator ($op$) in $NDP$

11:         Apply op to improve $C{S}_2$ and obtain $C{S}_3$

//Construct offspring population ($OP$)

12:          OP←[$OP$, $C{S}_2$, $C{S}_3$]

13:       end while

//Update population

14:         Merge population: $TP$←[P, $OP$]

15:         Update of population: Apply environmental selection of a MOEA

(i.e., NSGA-II, SPEA2, IBEA) to generate new population (NP)→ P

//External archive

16:         Merge population: $TAP$→[$AP$, $OP$]

17:         if the number of non-dominated results in $TAP$ exceeds $10\times Np$ then

Save the best $10\times Np$ individuals of $TAP$ into $AP$

else

Save the nondominated individuals of $TAP$ into $AP$

18:            end if

19: until T← $T_{max}$

20: Return $AP$

The algorithm returns the AP when the main loop stops. Algorithm 1 provides an overview of the pseudo-code of the proposed MOEA framework. However, Algorithm 1 is only suitable for the MOEA implementing FI. To provide fair comparisons between MOEAs using FI and MOEAs using BI, this paper provides the average computing time of Step 5–13 in MOEAs using FI as the termination conditions of the MOEAs using BI.

5. Computational Experiments and Analysis

Implementation aspects, computational experiments, and analysis are presented and discussed in the following sections.

5.1. Implementation and Parameters Setting

The proposed MOEA algorithms are implemented in MATLAB language and run on a PC with an Intel Core i5-6200K CPU at 2.40 GHz and 8 GB memory under the Windows 10 system.

Since the proposed problem is time-consuming and difficult to solve, the population size $N_p$ is set to 100. As to the size of the archive population, we save $10\times N_p$ individuals. For the maximum number of iterations ($T_{max}$) used in MOEAs using FI, it depends directly on the number of clients and depots:

$$T_{max}=\omega \times (N+M)$$

(36)

The multiplier $\omega$ depends on the size of the instance, taking on the value 30 for smaller instances (N = 40), 40 for mid-size instances (N = 50), and 50 for larger instances (N = 60). The larger the instance size, the harder it is to solve. Hence, the performance of the algorithm depends on $\omega$. In the IBEA, we follow the suggested value for the parameter $\kappa$ equalling to 0.05. For the mutational probability $P_m$, we set it by the initial experiments in Section 5.4.

5.2. Performance Indicator

This paper investigates the multi-objective GLRPCCL, hence, how to effectively evaluate the performance of the proposed MOEA framework is a significant issue for the multi-objective GLRPCCL. Li et al. [81] stated that the assessment of an MOEA is mainly from two aspects: the proximity to the RPF and the diversity of the approximate Pareto front (APF). Hence, this paper applied four widely used metrics to estimate the performance of the proposed approaches: hypervolume (HV) [82], inverted generational distance (IGD) [83], coverage (C) [82], and diversity metric (DM) [16].

HV measures the size of the coverage space between the APF and a reference point. The greater the HV value, the better the diversity and distribution. IGD measures the distance between RPF and APF. The lower the IGD value, the better the quality of the APF approximating the entire RPF. IGD is equal to 0, indicating that the obtained Pareto front contains every point of the RPF. The C is to compare the relative coverage of the two solution sets in MOEA. In this paper, we use the sets of RPF and APF to compare the non-dominated relationship between APF and RPF. A smaller value of this metric indicates the APF has better performance. The DM is adopted to evaluate the diversity of a set of non-dominated solutions. A larger DM value indicates a better diversity of APF.

However, in our case, the set of RPF is unknown. Therefore, according to the literature [84], we use all the APF sets obtained by all algorithms to form RPF after removing the repeated and dominated solutions. However, the set RPF is still an approximation of the true Pareto front.

5.3. Test Instances

The instances used in this paper were provided by [2]. In the randomly generated instances, the number of clients is$N\in \{40,50,60\}$, and the number of depots is $M\in \{6,7,8\}$. Clients and depots are randomly distributed in the square 50 km × 50 km. A uniform distribution is applied to generate the pickup and delivery demand for each client in a range [0.1 1.6] tons and the capacity of the depots in a range [10, 15] tons. The opening time window of each client was derived from the instance C101 in [85], then it is divided by 2, and the closing time window of each client relies on the service time (see Equation (37), that is, $l_i=e_i+s{t}_i(i\in C)$.

$${st}_i=1800\times \frac{d_i+p_i}{{\sum }_{i\in N}(d_i+p_i)/N}$$

(37)

For the operation fee of each depot, we set them in a range [500, 100] Yuan, and the closing time windows of all depots are set to 12 h. The speed of each edge is in a range [20, 70] km/h.

5.4. Parameter Tuning for $P_m$

In this section, we apply three well-known MOEAs combined with FI and BI to solve the instance with 40 clients, aiming at analyzing the effect of mutational probability $P_m$ and the performance of six pairs, that is, NSGA-II/FI, NSGA-II/BI, SPEA2/FI, SPEA2/BI, IBEA/FI, and IBEA/BI. A total of 12 independent runs are performed.

To effectively analyze the performance, the comparison analysis has three features: (1) we used the gaps of HV and IGD rather than the raw IGD and HV values since the raw IGD and values can’t reflect the whole performance for all instances and variants (see Equations (38) and (39)); (2) the non-dominated sorting strategies to obtain the best front construed by the gaps of HV and IGD; (3) if the first front organized by the gaps of HV and IGD have multiple points, we used the scoring system of the CHESC cross-domain heuristic search challenge (http://www.asap.cs.nott.ac.uk/external/chesc2011/), in other words, the top eight ones score 10, 8, 6, 5, 4, 3, 2, and 1, respectively. This paper uses the average gaps of IGD and HV, hence, the highest score is 20 for each pair. Aiming at analyzing the effect of mutational probability $P_m$, we set $p_m\in 0,0.2,0.4,0.6,0.8,1$.

Table 3. Average gaps of 12 runs for HV and IGD.

Instance	Algorithm	Average Gaps of 12 Runs (HV)						Average Gaps of 12 Runs (IGD)
		0	0.2	0.4	0.6	0.8	1	0	0.2	0.4	0.6	0.8	1
C40-6-1	NSGAII/FI	13	12.6	12	12.4	12.8	13.1	64.1	56.6	40.1	44.7	48.4	42.7
	NSGAII/BI	14.8	12.7	12.8	13.3	14.1	15.2	72.7	37.3	37.9	33.1	31.2	35.6
	SPEA2/FI	10.7	12.6	13.2	15	15.8	16.5	27.9	42.5	40.1	52	50.6	55.5
	SPEA2/BI	11	13.1	14.2	15.1	15.6	16.8	29.5	40.2	44.1	48.1	46.5	52.8
	IBEA/FI	7.6	7	7	7.1	7.2	7.8	71.5	38	26.1	32	21.3	16.5
	IBEA/BI	7.9	5.9	5.4	6.1	6.5	6.4	72.6	26.8	33.5	23.3	16	16.7
C40-6-2	NSGAII/FI	12.4	11.7	11	12.1	11.8	11.8	147.2	22.1	20	28.2	20	16.5
	NSGAII/BI	12.5	11.9	12.6	12.1	12.4	13.8	161	27.1	27.9	22.7	17	28.7
	SPEA2/FI	11.2	12.3	12.6	13.7	14.2	14.9	150.5	36.3	36.6	49	45.2	46.5
	SPEA2/BI	11.4	12.8	13.7	14.7	14.8	15.6	175.2	40.3	55.9	51.3	58.7	52.9
	IBEA/FI	10.6	8.6	8.3	8.1	8.6	8.6	390.6	200.1	139.9	102.2	39.1	36.1
	IBEA/BI	9.1	6.9	6.4	6.9	7.1	8.1	316.1	260.7	102.5	66.3	46.9	32.9
C40-6-3	NSGAII/FI	13.5	11.2	10.7	11	11.8	12.1	272.8	40.9	25	19.6	25.1	21.1
	NSGAII/BI	13.4	11	10.7	11.5	11.7	12.3	236.5	33.6	29.5	23.4	22.3	24.2
	SPEA2/FI	12.6	11.1	12.1	13.2	15.2	15.2	295.9	34.7	34.7	33.2	47.1	48.7
	SPEA2/BI	12	11.5	12.2	13.4	13.8	15.2	213.1	27.7	28.2	36	35.6	44.5
	IBEA/FI	11.8	8.2	6.8	6.4	6	7.1	400.1	228.7	140.1	94.9	59.8	48.5
	IBEA/BI	10.5	7.3	6.4	5.8	5.6	6.6	332.7	186.6	124.4	81.4	56.9	53.8
C40-6-4	NSGAII/FI	12.9	11.1	12.4	11.8	12.7	13.5	136.5	47.2	39.9	30.7	34.2	36.8
	NSGAII/BI	12.6	12.9	13.1	13.2	14.2	14.2	93.1	42.2	35.3	40.6	43.3	42.3
	SPEA2/FI	12.2	11.9	12.9	14.1	15.2	16.1	114.2	29.5	23.7	23.1	29.6	29.8
	SPEA2/BI	11.6	12.3	13.6	14.4	14.9	15.8	59.7	29.2	29.1	28.4	31.2	34.1
	IBEA/FI	11.3	7.9	9.1	7.9	8.2	8.7	189.4	63.5	84	52.3	65.9	35.6
	IBEA/BI	10.6	7.3	7.2	7.6	6.8	8.6	137.9	99.6	65.3	68.3	60.9	40.2
Average	NSGAII/FI	12.9	11.6	11.5	11.8	12.3	12.6	155.1	41.7	31.2	30.8	31.9	29.3
	NSGAII/BI	13.3	12.1	12.3	12.5	13.1	13.9	140.8	35	32.6	30	28.5	32.7
	SPEA2/FI	11.7	12	12.7	14	15.1	15.7	147.1	35.8	33.8	39.3	43.1	45.1
	SPEA2/BI	11.5	12.4	13.4	14.4	14.8	15.9	119.4	34.4	39.3	40.9	43	46.1
	IBEA/FI	10.3	7.9	7.8	7.36	7.5	8	262.9	132.6	97.5	70.4	46.5	34.2
	IBEA/BI	9.5	6.9	6.3	6.57	6.5	7.41	214.8	143.4	81.4	59.8	45.2	35.9

is the average gap of 12 runs for six MOEAs, and the average gaps of four instances for the HV and IGD values are plotted in Figure 3.

$${Gap}_{HV}\left(i\right)=\frac{{HV}_{RPF}-{HV}_i}{{HV}_{RPF}}\times 100\%$$

(38)

$${Gap}_{IGD}\left(i\right)=\frac{IG{D}_i-{min}_{i\in K}\left\{IG{D}_i\right\}}{{min}_{i\in K}\left\{IG{D}_i\right\}}\times 100\%$$

(39)

Looking at Table 3 and Figure 3, for average gaps of HV, IBEA/BI was the best in terms of the minimum average gaps of HV, regardless of the $P_m$ value, followed by IBEA/FI. When $P_m$ is larger than 0, NSGA-II/FI and NSGA-II/BI are better than SPEA2/FI and SPEA2/BI. In terms of the search strategy, IBEA/BI is better than IBEA/FI, however, NSGA-II/FI was better than NSGA-II/BI. For the SPEA2, the performance of SPEA2/FI using $p_m\in \{0.2,0.4,0.6,1\}$ is better than SPEA2/BI. Moreover, the average gaps of HV obtained by SPEA2 increase as $P_m$ increases, but the others are concave. Hence, from the perspective of average gaps of HV, IBEA/BI using $p_m=0.4$ was the best among six pairs.

For average gaps of IGD values, although the average gaps of HV obtained by IBEA are the best, the average gaps of IGD are the worst when $P_m$ < 0.8, and the IGD gaps of IBEA decrease along with as $P_m$ increases. The average gaps of IGD achieved by IBEA using $P_m$ = 1 are lower than that of SPEA2 but larger than that of NSGA-II. For the performance of SPEA2, the average gaps increase as $P_m$ increases when $P_m$ > 0 (except for SPEA2/FI using 0.2). For the performance of NSGA-II, with the change of $P_m$, the change on average gaps is little. However, the effect of BI and FI is difficult to identify. From the minimum gap, NSGA-II/BI using pm = 0.8 is the best, followed by NSGA-II/FI using $p_m=1$.

Through the above performance analysis, SPEA2 performed the worst. However, it is hard to identify the top four pairs for the following experiments. Hence, we used the non-dominated sorting strategy to find the first front, as shown in Figure 4. The results show that nine points construct the first front. Then we utilized the scoring system to rank nine pairs in the first front, and the scores show that NSGA-II/BI using 0.8 and IBEA/BI using 0.4 obtain 10 scores, and NSGA-II/FI using 1 and IBEA/BI using 0.8 obtain 9 scores, while others achieve 8 scores. Hence, the top four pairs (i.e., NSGA-II/BI using 0.8, IBEA/BI using 0.4, NSGA-II/FI using 1, and IBEA/BI using 0.8) are selected for the following experiments. Moreover, we selected the best Pareto fronts with the maximum HV values and minimum IGD values from 12 runs for each instance to plot the HV values and IGD values as iteration increases, as shown in Figure 5, where the final population is the one produced in the last iteration.

5.5. Contradictory Nature and Necessity of Researching Four Objectives

This section discusses the contradictory nature of four objectives: (1) TC and AWT. The first objective TC contains the fixed costs of vehicle and depots, as well as the drivers’ salary proportional to the total traveled time, and AWT are included in the total traveled time. However, the fixed costs of depots, the fixed costs of vehicles, and the drivers’ salary per minute depending on vehicle type could provide contradictory nature between TC and AWT. (2) TE. The TE contains CO${}_2$, CH${}_4$, and N${}_2$O emissions. And the carbon emission depends on the vehicle type, load, cargo type, speed, and traveled distance, while CH${}_4$ and N${}_2$O emissions are in line with the traveled distance. (3) CQD. The objective CQD depends on the cargo type, the traveled time with $s\ne 0$, and the waiting time with $s=0$, while the spoiled rate is different. However, there is a positive correlation between TC and TE/AWT, since the number of Pareto fronts of the instances with 40 clients is the fewest among six pairs constructed by two of the proposed four objectives, as shown in Figure 6, which also shows that there is certain contradictory nature among the proposed four objectives.

In terms of the necessity of researching four objectives, we provide the maximum and minimum values concerning TC (in Yuan), TE (in kg), AWT (in minutes), and CQD (in kg) in the Pareto front, as shown in

Table 4. Maximum and minimum values in the Pareto front.

Instance	Minimum Value				Maximum Value
	TC	TE	AWT	CQD	TC	TE	AWT	CQD
C40-6-1	3738.48	932.90	23.02	76.87	11,873.98	2376.73	187.43	159.98
C40-6-2	4129.70	953.86	25.42	119.71	13,660.60	2491.31	181.40	269.60
C40-6-3	3761.50	1051.80	22.05	117.86	13,117.77	2585.22	184.89	284.43
C40-6-4	3780.67	931.78	23.07	63.01	13,189.96	2108.46	159.12	169.93
C50-7-1	4912.56	1188.49	18.83	120.80	16,657.81	3292.93	196.43	285.33
C50-7-2	5345.49	1067.50	20.73	114.89	15,978.91	3259.81	192.12	250.48
C50-7-3	5118.51	1259.00	15.46	140.63	16,916.36	3856.75	190.59	336.35
C50-7-4	4580.33	1006.93	15.80	100.03	13,922.51	3258.54	158.12	238.41
C50-7-5	4662.00	1037.09	12.07	153.45	15,261.42	3625.81	205.85	417.56
C60-8-1	5565.93	1257.69	12.59	117.96	17,675.54	4015.26	165.35	301.91
C60-8-2	5206.41	1300.65	16.99	133.74	18,977.81	4809.71	201.07	356.98
C60-8-3	5130.46	1276.07	18.45	150.45	19,731.14	4582.69	198.99	381.13
C60-8-4	5188.25	1292.08	19.71	138.22	18,057.16	4157.19	183.95	386.58
C60-8-5	5434.07	1177.94	11.58	181.35	17,679.90	3684.71	174.90	513.57

. For the TC, TE, AWT, and CQD, the average gaps between the maximum and minimum values reach to 234.71%, 202.34%, 970.30%, and 149.09%, respectively. Hence, it is necessary to study the proposed multi-objective problem for the GLRPCCL. The reason is that if they can’t simultaneously consider economic, environmental, and social effects and the cargos’ quality, the decision-makers may make an imperfect judgment on the designs of the GLRPCCL.

5.6. Effect of Delivery Strategy

This section is conducted to analyze the benefit of mixed transportation (MT) compared to the traditional way named individual transportation (IT) that only one type of cargo can be assigned in a vehicle. The used instances in this section are the cases with 40, 50, and 60 clients.

Table 5. Performance indicators of IT and MT.

Instance	HV		IGD		C		DM
	IT	MT	IT	MT	IT	MT	IT	MT
C40-6-1	$2.1192\times {10}^{11}$	$2.1599\times {10}^{11}$	14.2566	2.8993	0.9468	0.7344	0.8214	0.9227
C40-6-2	$4.5283\times {10}^{11}$	$4.6690\times {10}^{11}$	19.5010	1.0978	0.9625	0.6617	0.7434	0.9675
C40-6-3	$5.3288\times {10}^{11}$	$5.3967\times {10}^{11}$	15.8131	1.7935	0.9680	0.6803	0.7778	0.9589
C40-6-4	$2.2240\times {10}^{11}$	$2.2865\times {10}^{11}$	16.2299	2.5124	0.9567	0.6774	0.8191	0.9332
C50-7-1	$9.8539\times {10}^{11}$	$1.0049\times {10}^{12}$	22.8042	0.7387	0.9852	0.9482	0.7086	0.9839
C50-7-2	$7.3610\times {10}^{11}$	$7.4825\times {10}^{11}$	24.3371	1.9301	0.9557	0.9458	0.7689	0.9522
C50-7-3	$1.4567\times {10}^{12}$	$1.4774\times {10}^{12}$	20.8781	1.7198	0.9701	0.9400	0.7712	0.9524
C50-7-4	$5.4693\times {10}^{11}$	$5.5718\times {10}^{11}$	18.1430	1.6732	0.9301	0.9090	0.7798	0.9529
C50-7-5	$1.8063\times {10}^{12}$	$1.9208\times {10}^{12}$	43.6377	0.4140	0.9961	0.9825	0.7256	0.9918
C60-8-1	$1.2246\times {10}^{12}$	$1.2687\times {10}^{12}$	30.2990	0.5709	0.9885	0.9669	0.7547	0.9844
C60-8-2	$2.6436\times {10}^{12}$	$2.6937\times {10}^{12}$	26.2407	1.7824	0.9760	0.9679	0.7561	0.9634
C60-8-3	$2.6939\times {10}^{12}$	$2.7430\times {10}^{12}$	30.3263	2.3379	0.9683	0.9567	0.7774	0.9511
C60-8-4	$2.0328\times {10}^{12}$	$2.0741\times {10}^{12}$	23.7116	1.3938	0.9715	0.9631	0.7328	0.9716
C60-8-5	$2.1694\times {10}^{12}$	$2.2769\times {10}^{12}$	35.5579	1.8574	0.9770	0.9787	0.7187	0.9608

is the performance indicators of Pareto fronts, that is, raw HV, IGD, C, and DM values.

On average, the proposed strategy could improve the HV value by 2.619%, the IGD value by 92.153%, the C value by 9.215%, and the DM by 26.519%. However, the IT could obtain better C value for the C60-8-5 compared to the proposed MT. For the other performance indicators of other instances, the proposed strategy performs better than IT, showing that mixed transportation can help in the design of the GLRPCCL in terms of environmental, economic, and social effects, as well as maintaining cargos’ quality.

5.7. Effect of Depots’ Time Window and Depot Cost

This section analyses the joint effects of depots’ time windows and depot cost from a novel perspective. Assuming that the depots’ costs depend on the depots’ time windows, that is, the unit of depots’ costs is Yuan/hour. Considering the original instances, the depots’ cost per hour in this section equals the original depots’ cost divided by 12 h. Moreover, we also use the variants with TWD ∈ {9, 15, 18, 21, 24} h. The results of these experiments are compared with the base case.

Table 6. Raw HV values of six variants of TWD.

Instance	V1(9)	V2(12)	V3(15)	V4(18)	V5(21)	V6(24)
C50-7-1	$1.257\times {10}^{12}$	$1.234\times {10}^{12}$	$1.197\times {10}^{12}$	$1.167\times {10}^{12}$	$1.132\times {10}^{12}$	$1.100\times {10}^{12}$
C50-7-2	$7.163\times {10}^{11}$	$6.936\times {10}^{11}$	$6.634\times {10}^{11}$	$6.372\times {10}^{11}$	$6.051\times {10}^{11}$	$5.737\times {10}^{11}$
C50-7-3	$1.550\times {10}^{12}$	$1.519\times {10}^{12}$	$1.465\times {10}^{12}$	$1.416\times {10}^{12}$	$1.376\times {10}^{12}$	$1.327\times {10}^{12}$
C50-7-4	$5.894\times {10}^{11}$	$5.742\times {10}^{11}$	$5.553\times {10}^{11}$	$5.390\times {10}^{11}$	$5.215\times {10}^{11}$	$5.041\times {10}^{11}$
C50-7-5	$1.324\times {10}^{12}$	$1.284\times {10}^{12}$	$1.223\times {10}^{12}$	$1.166\times {10}^{12}$	$1.121\times {10}^{12}$	$1.059\times {10}^{12}$
C60-8-1	$1.102\times {10}^{12}$	$1.075\times {10}^{12}$	$1.046\times {10}^{12}$	$1.010\times {10}^{12}$	$9.778\times {10}^{11}$	$9.480\times {10}^{11}$
C60-8-2	$2.118\times {10}^{12}$	$2.077\times {10}^{12}$	$2.009\times {10}^{12}$	$1.948\times {10}^{12}$	$1.885\times {10}^{12}$	$1.826\times {10}^{12}$
C60-8-3	$1.613\times {10}^{12}$	$1.568\times {10}^{12}$	$1.527\times {10}^{12}$	$1.470\times {10}^{12}$	$1.411\times {10}^{12}$	$1.366\times {10}^{12}$
C60-8-4	$1.832\times {10}^{12}$	$1.790\times {10}^{12}$	$1.740\times {10}^{12}$	$1.695\times {10}^{12}$	$1.646\times {10}^{12}$	$1.597\times {10}^{12}$
C60-8-5	$2.061\times {10}^{12}$	$2.000\times {10}^{12}$	$1.905\times {10}^{12}$	$1.840\times {10}^{12}$	$1.754\times {10}^{12}$	$1.685\times {10}^{12}$

and

Table 7. Raw IGD values of six variants of TWD.

Instance	V1(9)	V2(12)	V3(15)	V4(18)	V5(21)	V6(24)
C50-7-1	68.838	52.983	84.555	153.712	281.527	451.025
C50-7-2	66.469	67.152	102.654	207.666	385.319	640.475
C50-7-3	50.302	39.490	77.792	162.864	300.579	489.882
C50-7-4	48.622	44.106	56.633	119.794	222.746	373.270
C50-7-5	51.867	49.439	107.765	224.534	410.125	598.243
C60-8-1	45.200	50.228	74.601	144.240	279.142	444.744
C60-8-2	47.074	43.352	55.814	119.988	227.824	389.878
C60-8-3	40.975	45.443	67.780	112.527	269.238	440.275
C60-8-4	42.466	47.948	64.009	127.510	228.563	368.378
C60-8-5	67.540	62.766	90.947	193.527	362.582	585.554

detailed the raw values of HV and IGD values for each instance with six versions.

In terms of raw HV values, we find that as the TWD increases, the value of HV gradually decreases, which means that the lower the TWD, the wider the Pareto front is. Moreover, we calculate the gaps of HV values of each variant to the HV of RPF made of six variants, and we plot the average gaps in Figure 7, which shows that the gaps of HV are almost linear with TWD. Hence, in the context of the same operating cost per hour, lower TWD can help in designing the network of the GLRPCCL.

In terms of raw IGD value, for C50-7-1, C50-7-3, C50-7-4, C50-7-5, C60-8-2, and C60-8-5, the raw IGD values of V2(12) are lower than V1(9). The reasons may be that (1) The lower TWD have tighter restrictions on the assignments of clients and vehicles than higher TWD, but the operating costs of lower TWD allow the decisions that a larger set of depots can be opened; (2) The effectiveness of the proposed algorithm is poor. Figure 8 plots the average IGD gaps of each variant to the minimum IGD values. We can find that the average IGD gap is a polynomial or exponential equation of TWD.

From the above analyzes, we could conclude that the Pareto front of the instances using lower TWD is better than those using larger TWD in terms of HV, but there may exist the best TWD values for a special instance to obtain the minimum IGD values. Hence, the logistics enterprises should analyze the joint effects of depots’ cost and closing time windows to improve the performance of the GLRPCCL concerning multiple impacts.

5.8. Effect of Fleet Composition

We now analyze the impact of one of the operations decisions, namely fleet composition. In this part, we generate six variants for each original instance with different fleet composition, including h1, h2, h3, h1&h2, h1&h3, and h2&h3, The first three only consider the homogeneous vehicles, while the latter three are also the versions of HF that only consider two types of vehicles.

Table 8. Raw HV values of seven variants.

Instance	h1	h2	h3	h1&h2	h1&h3	h2&h3	h1&h2&h3
C50-7-1	$9.588\times {10}^{11}$	$9.587\times {10}^{11}$	$7.657\times {10}^{11}$	$9.943\times {10}^{11}$	$9.824\times {10}^{11}$	$9.706\times {10}^{11}$	$1.001\times {10}^{12}$
C50-7-2	$7.401\times {10}^{11}$	$7.192\times {10}^{11}$	$5.875\times {10}^{11}$	$7.485\times {10}^{11}$	$7.482\times {10}^{11}$	$7.265\times {10}^{11}$	$7.549\times {10}^{11}$
C50-7-3	$1.429\times {10}^{12}$	$1.418\times {10}^{12}$	$1.176\times {10}^{12}$	$1.461\times {10}^{12}$	$1.443\times {10}^{12}$	$1.425\times {10}^{12}$	$1.470\times {10}^{12}$
C50-7-4	$5.426\times {10}^{11}$	$5.283\times {10}^{11}$	$4.343\times {10}^{11}$	$5.499\times {10}^{11}$	$5.486\times {10}^{11}$	$5.300\times {10}^{11}$	$5.549\times {10}^{11}$
C50-7-5	$1.479\times {10}^{12}$	$1.827\times {10}^{12}$	$1.655\times {10}^{12}$	$1.825\times {10}^{12}$	$1.792\times {10}^{12}$	$1.886\times {10}^{12}$	$1.867\times {10}^{12}$
C60-8-1	$1.236\times {10}^{12}$	$1.209\times {10}^{12}$	$9.951\times {10}^{11}$	$1.263\times {10}^{12}$	$1.247\times {10}^{12}$	$1.218\times {10}^{12}$	$1.268\times {10}^{12}$
C60-8-2	$2.603\times {10}^{12}$	$2.602\times {10}^{12}$	$2.254\times {10}^{12}$	$2.665\times {10}^{12}$	$2.644\times {10}^{12}$	$2.613\times {10}^{12}$	$2.677\times {10}^{12}$
C60-8-3	$2.666\times {10}^{12}$	$2.647\times {10}^{12}$	$2.254\times {10}^{12}$	$2.712\times {10}^{12}$	$2.689\times {10}^{12}$	$2.650\times {10}^{12}$	$2.725\times {10}^{12}$
C60-8-4	$1.993\times {10}^{12}$	$1.998\times {10}^{12}$	$1.670\times {10}^{12}$	$2.053\times {10}^{12}$	$2.030\times {10}^{12}$	$2.006\times {10}^{12}$	$2.063\times {10}^{12}$
C60-8-5	$1.807\times {10}^{12}$	$2.187\times {10}^{12}$	$1.890\times {10}^{12}$	$2.196\times {10}^{12}$	$2.106\times {10}^{12}$	$2.219\times {10}^{12}$	$2.232\times {10}^{12}$
f/s/t	0/0/0	0/0/0	0/0/0	0/8/2	0/0/8	1/1/0	9/1/0

,

Table 9. Raw IGD values of seven variants.

Instance	h1	h2	h3	h1&h2	h1&h3	h2&h3	h1&h2&h3
C50-7-1	28.3382	28.1789	164.4297	17.1665	21.6414	26.7950	10.4875
C50-7-2	22.2361	31.3098	130.5514	18.8328	21.4808	27.6786	11.4050
C50-7-3	18.2962	25.5281	172.3676	17.9380	21.1502	24.0540	11.2074
C50-7-4	17.4139	23.1922	114.0420	15.3664	15.9547	21.7611	9.7126
C50-7-5	221.2876	32.2058	106.5927	37.6615	40.2009	25.8714	24.1445
C60-8-1	24.4832	30.6912	165.9259	18.7288	22.7786	28.1751	10.8799
C60-8-2	31.9548	28.5600	151.1358	21.5366	24.6615	26.5044	14.2788
C60-8-3	30.3134	31.8175	157.9705	24.1734	27.7285	31.1875	15.4287
C60-8-4	26.9337	25.6547	172.5683	18.2267	23.6548	23.9208	10.7317
C60-8-5	190.7419	29.9549	129.8395	29.5820	47.3324	24.9296	18.9193
f/s/t	0/0/1	0/0/1	0/0/0	0/8/1	0/0/6	0/2/1	10/0/0

and

Table 10. Raw IGD values of seven variants.

Instance	h1	h2	h3	h1&h2	h1&h3	h2&h3	h1&h2&h3
C50-7-1	0.9541	0.9925	0.9990	0.8911	0.9481	0.9853	0.8681
C50-7-2	0.9557	0.9979	0.9987	0.9021	0.8852	0.9918	0.8616
C50-7-3	0.9401	0.9953	0.9998	0.8929	0.8925	0.9861	0.8545
C50-7-4	0.9035	0.9957	0.9997	0.8407	0.7733	0.9903	0.8128
C50-7-5	0.9820	0.9812	0.9936	0.9411	0.9886	0.9365	0.9413
C60-8-1	0.9625	0.9989	0.9994	0.9160	0.9121	0.9948	0.8718
C60-8-2	0.9664	0.9949	0.9999	0.9442	0.9331	0.9894	0.9095
C60-8-3	0.9618	0.9928	0.9994	0.9255	0.8965	0.9858	0.9029
C60-8-4	0.9674	0.9907	0.9999	0.9093	0.9264	0.9752	0.9060
C60-8-5	0.9879	0.9827	0.9971	0.9481	0.9830	0.9609	0.9464
f/s/t	0/0/0	0/0/0	0/0/0	0/4/6	2/4/2	1/0/1	7/2/1

present the raw HV, IGD, and C values for the instances with seven variants, and the last row represents the number of rank places at first, second, and third places of seven variants for 10 instances.

Looking at Table 8, the instances using h1&h2&h3 as fleet could obtain nine maximum HV values and averagely improve the HV values about 7.245%, 3.622%, 23.349%, 0.870%, 2.184%, 2.667% when compared to h1, h2, h3, h1&h2, h1&h3, and h2&h3, respectively. However, the maximum HV value is obtained by h2&h3 for the instance C50-7-5. From a detailed perspective of the HV values obtained by the homogeneous fleet, eight instances prefer the use of h1, followed by h2. However, the instances C50-7-5 and C60-8-5 prefer the use of h2, followed by h3. In terms of the heterogeneous fleet using two types of vehicles, the HV values of the instances using h1&h2 are the maximum among three variants, followed by h1&h3. But the instance C50-7-5 and C60-8-5 uses the h2&h3 to achieve the highest HV values. The reason may be that clients’ pickup and delivery demand and time windows may affect the selection of type of vehicles. Hence, logistics companies should analyze the effects of fleet composition, which is a significant factor in affecting the performance indicators of the GLRPCCL.

Table 9 shows the raw IGD values for the instances using seven variants. We can find that all instances using h1&h2&h3 can obtain the minimum IGD values, which are, respectively, higher than h1, h2, h3, h1&h2, h1&h3, and h2&h3 about 59.40%, 52.67%, 90.61%, 37.75%, 47.79%, and 47.24%. The instance C50-7-5 and C60-8-5 using h1 achieve the maximum IGD value, showing that it prefers the usage of h2 rather than h1. The instances C50-7-1, C50-7-5, C60-8-2, C60-8-4, and C60-8-5 using h2 could obtain lower IGD values than that using h1, and the instances C50-7-5 and C60-8-5 using h2&h3 can achieve the best IGD values when compared to the variants using two types of vehicles. Hence, the above conclusion is examined and effective.

Table 10 presents the raw C values for the instances using seven different fleet compositions. The instance using h1&h2&h3 could obtain the minimum C values except for the C50-7-4, C50-7-5, and C60-8-3. The above similar conclusions are also be drawn for Table 10. The Pareto solutions of the instances using h1&h2&h3 should be made of the other six pairs, but this paper can’t achieve the perfect point. The reasons are that: (1) the proposed multi-objective GLRPCCL is NP-hard and is much more difficult to solve than the LRP and VRP; (2) it is difficult to obtain all individuals located in the Pareto front since the number of individuals of the first front is plentiful and over ten thousand; (3) the efficiency of the proposed algorithm can’t be guaranteed for all instances. However, from the results, the instances using HF could obtain better Pareto fronts when compared to those using a homogeneous fleet.

5.9. Management Insights

This paper investigates the sustainable development for the GLRPCCL from the environmental, economic, social, and cargos’ effects by defining the multi-objective model for the GLRPCCL. Based on the above experiments, several management implications can be drawn as follows.

Section 5.6 carried out the experiments to analyze the effect of delivery strategies, and the results show that the mixed transportation could obtain the Pareto fronts dominated the fronts of the individual transportation in terms of performance indicators. Hence, the mixed transportation is the best choice for cold chain logistics.

Section 5.7 conducted the experiments to analyze the joint effects of depots’ time windows and operating costs, which illustrated that when the operational cost per hour is identical, the lower time windows could obtain better performance indicators for the Pareto solutions. Hence, the logistics companies should consider the joint effects of depots’ time windows and operating costs to obtain a good delivery network before serving the clients.

Section 5.8 investigated the benefits of using a heterogeneous fleet in designing the delivery network for the GLRPCCL. Compared to the homogeneous fleet and heterogeneous fleet made of two types of vehicles, the heterogeneous fleet using three types of vehicles could obtain the best performance indicators of Pareto front. In terms of comparison analysis, the heterogeneous fleet using two types of vehicles also performs well. Moreover, the fleet composition is jointly affected by the time windows, pickup demand, and delivery demand, which results in that there is a special type of vehicle that may be the best choice for an instance. Hence, before assigning the vehicles, the logistics enterprises should determine the appropriate fleet composition to avoid wasting resources.

6. Conclusions

GLRPCCL is an important issue in logistics and is quite crucial to reach the sustainable development concerning economic, environmental, and social effects as well as the cargos’ effect. For this purpose, this paper studied a GLRPCCL with multiple features: the time windows, the heterogeneous fleet, simultaneous pickup and delivery, and mixed transportation. The main contribution of this paper is the investigation of defining the model for the GLRPCCL concerning multiple effects.

A multi-objective MILP formulation was proposed which minimizes four impacts. In the defined model, the first objective is to optimize the total logistics costs including the fixed costs of depots and fleet as well as the costs of drivers’ salaries. The second is defined by the total emissions made of CO${}_2$, CH${}_4$, and N${}_2$O emissions. The third concerns the social impacts by minimizing the average waiting time of vehicles and clients. The last one emphasizes the quality degradation of cargos which is one of the main purposes of cold chain logistics. Since this problem is NP-hard, an effective framework of MOEA was developed. In the framework, two search strategies were utilized, and a large composite neighborhood made of 16 operators was developed to obtain Pareto solutions.

Afterward, the experiments were conducted by benchmark instances. The first experiment was carried out to analyze the sensitivity of mutational probability for obtaining Pareto fronts. The analysis illustrated that different MOEAs favor different parameter settings, and IBEA and NSGA-II are, respectively, the best in achieving maximum HV values and minimum IGD values. Furthermore, we also analyzed the contradictory nature among four objectives and provide the descriptions on the necessity to study the proposed model. The results show that our proposed model could improve the economic, environmental, social effects and could reduce the quality decay of cargos. Finally, we carried out extensive experiments to analyze the effects of the delivery strategy, depots’ costs and closing time windows, and fleet composition on the Pareto solutions of the GLRPCCL. And the conclusions are drawn as management insights.

The paper under consideration has some limitations. First, the algorithms could be improved to obtain good Pareto fronts, which may be affected by the mutational probability. Second, the type of cargos of each client is the same. For future research, we will develop a high-efficiency MOEA framework with free parameters. And we also will improve the model by defining that the cargos’ types of each client are different. Meanwhile, various uncertainties of the problem, such as uncertainties of travel time and clients’ demands, will be considered, and fuel supply at stations may also be concerned.