Customer-Oriented Vehicle Routing Problem with Environment Consideration: Two-Phase Optimization Approach and Heuristic Solution

Abstract

With the fastest consumer demand growth, the increasing customer’s demands trend to multivarieties and small-batch and the customer requires an efficient distribution planning. How to plan the vehicle route to meet customer satisfaction of mass distribution as well as reduce the fuel consumption and emission has become a hot topic. This paper proposes a two-phase optimization method to handle the vehicle routing problem, considering the customer demands and time windows coupled with multivehicles. The first phase of the optimization method provides a fuzzy hierarchical clustering method for customer grouping. The second phase formulates the optimization en-group vehicle routing problem model and a genetic algorithm to account for vehicle routing optimization within each group so that fuel consumption and emissions are minimized. Finally, we provide some numerical examples. Results show that the two-phase optimization method and the designed algorithm are efficient.

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.

1. Introduction

The vehicle routing problem (VRP) is a well-known NP-hard problem which has introduced by Dantzig and Ramser [1]. Generally, the VRP is to generate sets of vehicle routes to reduce travel cost. It is worth noting that the VRP has important practical significance of social benefits and enterprisers benefits. Especially, with the development of the economic and electrical business, the logistic delivery acts as a key role. From the survey data, we can find out that the carbon emission of logistic accounts for 80% of total road transport. With the continuous development of commodity economy, it becomes important to process customer orders and deliver them via suitable freight vehicles. Furthermore, customers need to be quickly responded for good delivery, which also become a hot topic.

Previous researches on VRP mostly take into account the economic factors, such as distance, cost, and time windows. And varieties of heuristic algorithm are designed to solve these problems [2, 3]. In recent years, how to reduce the energy consumption and carbon emissions has become a hot topic [4, 5]. Especially, with the rapid development of online shopping and express delivery, there is a sharply incasement of tail gas emissions from fuel logistics vehicles. Therefore, reducing the energy consumption and carbon emission has great significance for the ecological system. For logistics transportation enterprisers, the fuel consumption cost is one of the main operating costs, which is directly related to the carbon platoon, which will affect the operators’ profits and social benefits. Enterprises can reduce the cost of fuel consumption and greenhouse gas emission to achieve the goal of energy conservation and emission reduction. Therefore, how to optimize the routes of distribution vehicles and how to reduce energy consumption and carbon emission in the premise of meeting customers demand are the focus of research for logistics transportation enterprisers.

However, with the continuous development of the commodity economy, customers have put forward higher requirements of service quality of logistics enterprises. So how to improve the quality of distribution service has become the focus of enterprises. In particular, customer demands tend to be diversified, which proposes higher requirements for delivery service and will lead to increased costs. Providing appropriate service for each customer can improve the satisfaction, gain customers recognition, and loyalty to the enterprisers as well as improving the market competitiveness of enterprise. Deliver the required quantity of goods by reasonable arrangement of distribution vehicle routes. That is the last ring of the logistics service and the last part of the customer orders, which will directly affect the customers' satisfaction. Inappropriate distribution services result in customers' terrible experience and maybe lead to loss of customers and even worse increase the hidden costs of the enterprises. So it is equally important to reduce transportation costs as well as improve the service quality. To achieve a win-win situation for the vehicle distribution path designing, an appropriate clustering method is proposed to divide the customers into reasonable groups. That can effectively improve the efficiency of the people on both sides of the enterprisers and customers.

Therefore, this paper proposes a two-phase optimization method to design the vehicle delivery paths with comprehensive consideration of customer demand, energy conservation, and emission reduction. For the first phase, the fuzzy system clustering method (FSCM) is applied to reasonably group the customers according to the qualitative and quantitative indicators. The second phase is en-group path optimization; an en-group vehicle routing optimization model is established with the target of reducing energy consumption and emission reduction. To solve the problem, a genetic algorithm (GA) is designed. Based on the proposed method, we can improve the quality of distribution service, as well as reducing energy consumption and carbon emission in the distribution process.

Our paper makes a number of contributions to the literature.

$(1)$ We consider the customer requirements for vehicle routing problem. On this basis, a fuzzy hierarchical cluster method is designed, which can process quantitative and qualitative data. In order to guarantee the delivery quality, the customers can be clustered into different groups based on the order information. Then the problem can be converted into an en-group path optimization problem based on clustering results, so that can guarantee the delivery quality.

$(2)$ Energy saving and emission reduction are also the targets of our study. To characterize the problem, we formulate the model under both deterministic and dynamic traffic conditions, separately. As we know, the energy consumption is affected by the different vehicle speed. So considering the variety of the vehicle speed is more meaning for the study of reduced energy consumption. On this basis, we design a GA to solve the problem and generate an approximate optimal vehicle routing which can gain better solution results.

$(3)$ A two-phase model is proposed to meet the customer demands as well as reduce the fuel consumption. Contrasting to the common model, the two-phase model can improve the solving efficiency and more suitable to deal with customer requirements. The first phase model is a customer-oriented mathematics model which can process customer demands, and at this phase we group the customers into appropriate group, which can improve customers' satisfaction. The second phase model is an en-group VRP that aims to minimize the comprehensive costs including fuel and carbon costs and time penalty costs.

The remainder of this paper consists of five sections. Section 2 discusses previous related work and highlights the contributions of this paper. Section 3 provides the description of the proposed model. In Section 4, we provide the fuzzy hierarchical cluster method and en-group GA algorithm to solve the problem. We describe the design of experiments and numerical results in Section 5. Finally, we conclude our work and provide the propositions for further study in Section 6.

2. Literature Review

The VRP was proposed by Dantzig and Ramser and has been used to solve the optimal distribution route plan of trucks. As the question is raised, it has received wide attention from experts and scholars, especially in the field of operations research. Our work is closing to the multivehicles and green VRP, which considers fuel consumption and time windows. According to the research of different problems, various types of VRP are generated, and corresponding mathematical models are established [6, 7]. Brandao studied the Multitypes VRP and used a tabu search algorithm (TS to solve the problem, analyzed the influence of the vehicle type on the solution effect). To gain optimal results, types of algorithms were designed to solve the problems [8–10]. Specifically, Yazgi designed a greedy stochastic adaptive memory search algorithm to solve mixed fixed vehicle routing problems. Demir et al. [11] designed a large neighborhood search algorithm (VLSN) to solve the problem. In this literature, we review the existing research from three perspectives.

In order to reduce the influence of the environment, in recent years, some researchers turn on to study the vehicle routing problems with consideration of fuel consumption (VRPFC). Bektas et al. [12] studied the carbon emissions from VRP issues in several articles, summed up the relevant calculation model for vehicle energy consumption, and gave the effective application scope of the model. Palmer [13] established an integrated vehicle routing model and proposed a measurement model to calculate the $C{O}_{\mathrm{2}}$ emission. Miguel [3] aimed to reduce carbon emissions and fuel consumption and defined these problems as emissions vehicle routing problem (EVRP). Demir et al. [14] proposed a biobjective model to formulate the pollution routing problem (PRP), with the objective of minimizing the fuel consumption as well as travel time. Kwon et al. [15] studied the problem of vehicle path optimization with consideration of fuel consumption and carbon emissions based on time-varying velocity. The authors used the time-insertion algorithm to solve the single path problem and designed a column generation algorithm based on tabu search algorithm to solve the multipath problem. Huang [4] studied the vehicle routing problem of energy conservation and emission reduction and proposed the concept of green vehicle routing problem. Qian et al. [16] established a heterogeneous vehicle routing model by considering the carbon emissions during the vehicle operation; a tabu search algorithm was designed to solve the problem. Most of the VRP were formulated as mixed interlinear models, and heuristic algorithms are designed to solve the problem. The existing studies use the emission model which include MEET and National Atmospheric Emission Inventory [5, 17].

The vehicle routing problem with time windows was proposed by Avelsbergh [18] and widely studied by the researchers. Most of the studies focus on design better algorithms to solve the problem, such as GA, TS, etc. [8, 19, 20]. Specifically, Calvete [21] established a multiobjective optimization model for the vehicle routing problem (VRPSTW), which considered soft time windows and designed a goal planning method to solve the problem. Miguel [3] proposed an iterative route construction and improvement algorithm (IRCI) based on path generation algorithm to solve the VRPTW problem. Cihan et al. [2] proposed a time-dependent dual-target VRP. In the established model, it was assumed that the path selection will be influenced by the time windows.

While many contributions have been made to the VRP with time windows and the green vehicle problem, to our knowledge, fewer studies have examined the customer demands. However, processing the customer orders and quick response to varieties of customer demands for delivering the goods via suitable freight vehicles appears to be a key element for efficient logistics service. In previous literature [22], it has been pointed out that customer-based order clustering as well as the vehicle routing strategies should be coordinated to solve the vehicle problem. With the development of the VRP, a lot of studies focus on mathematically formulating the corresponding problems. In our research, we first aim to reasonably group the customers based on the differences of customer demand. On this basis, we target to optimize the vehicles’ travel cost including fixed cost, energy consumption, carbon emission, and punishment for time windows. For simplicity, the problem in this study is called the customer-oriented vehicle routing problem considering energy conservation and emission reduction (C-VRPESER).

In this research, we point out seven factors which will affect the customer satisfaction based on the studies [23], such as incorporating location analysis, service quality, service time, external similarity factor, etc. Furthermore, we establish the fuzzy similar matrix to cluster the customers with similar product demand. On the basis of the clustering results, we formulate the en-group optimization model and design a GA to solve the model. This study aims to make the following theoretical contribution to VRP which considers customer demand as well as energy conservation and emission reduction. Then the following two-phase optimization model is proposed to solve the C-VRPESER.

3. Problem Description and Model Formulation

In this paper, we consider a transportation network, that is represented by a directed complete graph $M=(N,A)$. The vertex set $N$ is composed by a depot 0 and a set of customers $N_{\mathrm{0}}=N\setminus \mathrm{0}=\{\mathrm{1,2},\mathrm{3},\dots n\}$. The number of homogeneous vehicles is a deterministic exogenous parameter and the set of vehicles is represented by $K=\{\mathrm{1,2},\mathrm{3},\dots m\}$. The capacity of each vehicle is equal to $Q_e$. Each vehicle $k$ has the earliest departure time $e_k$ from the depot after unloaded and the latest arrival time $l_k$ at the depot after finishing the services of customers. In practically, we can get the customer $i$ demand for goods from the previous order, denoted by ${\phi }_i^{\mathrm{1}}$, and we also can get the other demand for different attributes ${\phi }_i^{\mathrm{2}}$, ${\phi }_i^{\mathrm{3}}$, ${\phi }_i^{\mathrm{4}}$, ${\phi }_i^{\mathrm{5}}$, ${\phi }_i^{\mathrm{6}}$, and ${\phi }_i^{\mathrm{7}}$. And to simplify the problem, we divide the seven attributes into two types, quantitative and qualitative. For the qualitative demand, we use triangle fuzzy number (${\alpha }_{\mathrm{1},t}^p$, ${\alpha }_{\mathrm{2},t}^p$, and ${\alpha }_{\mathrm{3},t}^p$) to present five demand levels (“very high", “high", “medium", “low", and “very low"), in which $t$ presents the customer number and $p$. On the other hand, mathematical method is used to eliminate the dimensional and evaluate the similarity. On this basis, we cluster the customers into different groups, as shown in Figure 1. The customers can be clustered into three groups, which can be presented by different shapes, triangle, quadrilateral, and pentagon. And en-route goods delivery paths are designed according to the customers' location and time windows, which decide the delivery sequence. And then we connect each customer with a line with arrowheads. And in each group, the vehicle should start from the depot (the vehicle remained at the depot) and return to the depot after serving all the en-group customers.

[figure omitted; refer to PDF]

The optimization goal of the problem is to minimize the total cost on the basis of customers' group that realizes the en-group distribution vehicle routing optimization. By reasonably scheduling different types of vehicles for each customer group, the utilization rate of vehicles can be improved. And it also can meet the restriction of vehicle load and time limit. The optimal objective includes the vehicle fixed costs, time penalty cost, fuel cost, and carbon emissions. On this way, we can realize minimizing the total cost considering energy conservation and emissions reduction and reducing the cost of comprehensive optimization.

In the next sections, this paper focuses on formulating a two-phase model to solve the problem. The first model is generated to cluster the customers into groups and the second phase model aims to give out an optimal route en-group which can reduce the energy consumption and emission.

3.1. First Phase Model

The customers demand attributes consist of several indicators; part of them is quantitative indicators, such as time and location. Others are qualitative indexes, such as external and internal characteristics of commodity, customer satisfaction, etc. Through considering the related factors of distribution enterprisers service quality and customer satisfaction, we choose seven indicators as decision variables, $(1)$ physical properties, $(2)$ geographical position, $(3)$ service quality, $(4)$ product value, $(5)$ type of goods, $(6)$ security of goods, and $(7)$ time windows. In order to quantify the decision variables, we define the seven decision variables in Table 1. And the above indicators can represent most of the customer’s demand attributes.

Table 1

Definition of customer’s demand attributes.

Variables	Variable definitions	Variable explanation
$${\phi }_i^{\mathrm{1}}$$	Physical properties	Number of delivery of goods, weight, volume, this article take the quality as reference, the vehicle’s rated load determines the number of delivery vehicles.
$${\phi }_i^{\mathrm{2}}$$	Geographical position	Customers’ location, customers with similar geography have larger possibility be divided into same group,so that realize the joint distribution.
$${\phi }_i^{\mathrm{3}}$$	Service quality	Including reaction time, service attitude and service ability. Customer demands for service quality are different,and customers with similar service can be divided into same groups.
$${\phi }_i^{\mathrm{4}}$$	Product value	Value of products, generally refers to the market value. The high value cargo can be separated from others, this article define the value of the market price of five levels.
$${\phi }_i^{\mathrm{5}}$$	Type of goods	Types of goods and external similarity with other goods. If the external similarity is higher, that it convenient to load and unload the cargo which will improve the delivery efficient.
$${\phi }_i^{\mathrm{6}}$$	Security of goods	This property is mainly denote the customer security requirements(e.g., fragile). Especially for the dangerous cargo,the safety is very important.
$${\phi }_i^{\mathrm{7}}$$	Pressed for time	Present the customer request delivery time limit of time, the latest deadline, we can arrange the customer according to their time windows.

However the decision variables should eliminate dimension before using. For modeling convenience, this paper uses the fuzzy clustering method to cluster the customers into appropriate groups, which have the similar demand attributes. Taking into considering the above seven decision variables, the qualitative indicators include ${\phi }_i^{\mathrm{3}}$, ${\phi }_i^{\mathrm{4}}$, ${\phi }_i^{\mathrm{5}}$, and ${\phi }_i^{\mathrm{6}}$ and the other three variables are quantitative indicators.

In order to calculate the qualitative data, the linguistic terms “very high", “high", “medium", “low", and “very low" are introduced to represent the qualitative decision variables. Specifically, language “very high" can be represented by triangular fuzzy number (0.75, 1, 1), which indicates that the customer’s demand of the service level is much higher than other customers. Similarly, the language “high", “medium", “low", and “very low" can be represented by (0.5, 0.75, 1), (0.25, 0.5, 0.75), (0, 0.25, 0.5), and (0,0,0.25). Based on the language description, we translate the linguistic variables into triangular fuzzy numbers, so that we can calculate the comprehensive similarity between the customers. The linguistic variables represent the correlation between this indicator and customer satisfaction. In this way, each customer’s demand indicators can be denoted by triangular fuzzy number (${\alpha }_{i,\mathrm{1}}^p$, ${\alpha }_{i,\mathrm{2}}^p$, ${\alpha }_{i,\mathrm{3}}^p$), which can be denoted as follows:$$\begin{array}{}\text{(1)}& {\phi }_i^p=\left[{\alpha }_{i,\mathrm{1}}^p,{\alpha }_{i,\mathrm{2}}^p,{\alpha }_{i,\mathrm{3}}^p\right]\end{array}$$

We can get the demand attributers from the orders, which will be pregiven by the customers. However, we should transfer the linguistic variables and the quantitative data as standardized data, so that we can use them to calculate the similarity. And the procedure is shown in below.

(1) The Processing of Qualitative Decision Variables. In the first step we use the standard deviation transformation to process the qualitative data, which can guarantee the data availability. And for the second step, we introduce the range conversion to standardize the quantitative variables. The constraint (1) can be firstly transformed as follows:$$\begin{array}{}\text{(2)}& {\tilde{\phi }}_i^p=\left[{\tilde{\alpha }}_{i,\mathrm{1}}^p,{\tilde{\alpha }}_{i,\mathrm{2}}^p,{\tilde{\alpha }}_{i,\text{,}\mathrm{3}}^p\right],\hspace{1em}i\in N_{\mathrm{0}}\end{array}$$

where$$\begin{array}{}\text{(3)}& {\tilde{\phi }}_i^p=\frac{\left|{\alpha }_{i,t}^p-\overline{{\alpha }_{i,t}^p}\right|}{S_t^p},\hspace{1em}i\in N_{\mathrm{0}}\text{(4)}& \overline{{\alpha }_{i,t}^p}=\frac{{\sum }_{i=\mathrm{1}}^N{\alpha }_{i,t}^p}{N},\hspace{1em}t=\left\{\mathrm{1,2},\mathrm{3}\right\},\hspace{0.17em}\hspace{0.17em}i\in N_{\mathrm{0}}\end{array}$$

However the indicator $\overline{{\alpha }_{i,t}^p}$ is the mean value, $S_t^p$ is the standard deviation, and $\left|N\right|$ is the total customer number. The value range of $p=\{\mathrm{3,4},\mathrm{5,6}\}$.$$\begin{array}{}\text{(5)}& S_t^p={\left[\frac{{\sum }_{i=\mathrm{1}}^N{\left({\alpha }_{i,t}^p-\overline{{\alpha }_{i,t}^p}\right)}^{\mathrm{2}}}{\left|N\right|-\mathrm{1}}\right]}^{\mathrm{1}/\mathrm{2}},\hspace{1em}t=\left\{\mathrm{1,2},\mathrm{3}\right\},\hspace{0.17em}\hspace{0.17em}i\in N_{\mathrm{0}}\end{array}$$

Then we can calculate the similarity of qualitative decision-making variables with Hemingway’s distance method, as follows:$$\begin{array}{}\text{(6)}& X_{i,j}^p=\mathrm{1}-\frac{\left({\sum }_{t=\mathrm{1}}^{\mathrm{3}}\left|{\tilde{\alpha }}_{i,t}^p-{\tilde{\alpha }}_{j,t}^p\right|\right)}{\mathrm{3}},\hspace{1em}i,j\in N_{\mathrm{0}}\end{array}$$

$X_{i,j}^p$ represents the qualitative correlation between the different customers, and the range of the value is during $[\mathrm{0,1}]$. $X_{i,j}^p$ indicates relationship of attributors ${\tilde{\phi }}_i^p(k)$ and ${\tilde{\phi }}_j^p(k)$.

(2) The Processing of Quantitative Decision Variables. Quantitative decision variables are determined values. However, the dimension of the index is different. In order to eliminate the dimension, the quantitative data are needed to be normalized and converted into standardized data. Then, the normalized variables are treated as a kind of special fuzzy number, which can be used to calculate the similarity of the corresponding quantitative decision variables between two customers.$$\begin{array}{}\text{(7)}& {\tilde{\phi }}_i^p=\frac{{\phi }_i^p-{\min }_{\mathrm{1}\le n\le \left|N\right|}\left\{{\phi }_n^p\right\}}{{\max }_{\mathrm{1}\le n\le \left|N\right|}\left\{{\phi }_n^p\right\}-{\min }_{\mathrm{1}\le n\le \left|N\right|}\left\{{\phi }_n^p\right\}},\hspace{1em}\mathrm{f}\mathrm{o}\mathrm{r}\hspace{0.17em}\hspace{0.17em}p\in \left\{\mathrm{2,7}\right\},\hspace{0.17em}\hspace{0.17em}i\in N_{\mathrm{0}}.\end{array}$$

${\phi }_i^p$ and ${\phi }_n^p$ are the actual value of the quantitative decision variables corresponding to customers $i$ and customers $n$. $N$ is the number of customers. On the basis of the standardization of the quantitative decision, the similarity can be calculated by $$\begin{array}{}\text{(8)}& X_{i,j}^p=\mathrm{1}-\frac{\left|{\tilde{\phi }}_i^p-{\tilde{\phi }}_j^p\right|}{{\max }_{\mathrm{1}\le n\le \left|N\right|}\left\{{\tilde{\phi }}_n^p\right\}-{\min }_{\mathrm{1}\le n\le \left|N\right|}\left\{{\tilde{\phi }}_n^p\right\}},\hspace{1em}\mathrm{f}\mathrm{o}\mathrm{r}\hspace{0.17em}\hspace{0.17em}p\in \left\{\mathrm{2,7}\right\},\hspace{0.17em}\hspace{0.17em}i\in N_{\mathrm{0}}.\end{array}$$

Additionally, ${\phi }_i^{\mathrm{1}}$ represents the weight of the goods; it can determine the number of customer groups. And it is used as the termination condition of the customer clustering algorithm in this article. So we do nothing for this variable.

However, the comprehensive similarity for customers is not simply adding up the similarity for each decision variable. The variables present different influence for the delivery quality. In this paper, we consider the importance of each decision variable and use a mathematical method to calculate the weight coefficient for each parameter. A fuzzy mathematics method is introduced to determine the weight of each decision variable.

(3) Weight Calculation of Decision Variables. After obtaining the order information, the managers (experts) of the distribution center are required to submit evaluation reports about the delivery service index. The evaluation value of the impact of each decision variable for the customers is also given out by the managers. The experts' evaluation language is generally described as the influence of decision variables on customer satisfaction. We divide the areas of influence into five evaluations, specifically, largest, large, medium, low, and lowest. To calculate the numerical value, the language evaluation is converted to triangular fuzzy numbers and shown in Table 2.

Table 2

Membership function of evaluation linguistic variables.

Language terminology	Triangular fuzzy number	Judgment scale
largest	(0.75,1,1)	1
large	(0.5,0.75,1)	0.75
medium	(0.25,0.5,0.75)	0.5
small	(0,0.25,0.5)	0.25
smallest	(0,0,0.25)	0

Based on the importance of the decision variables for customer satisfaction, the weight of the corresponding decision variables can be calculated. For convenience, we denote some key parameters as follows.

$W_{u,p}(u=\mathrm{1,2},\dots ,r,\hspace{0.17em}\hspace{0.17em}p=\mathrm{1,2},\dots ,m)$ represents the evaluation of decision variables $p$ for decision maker $u$, where $r$ is the total number of decision makers and $m$ is the number of decision variables. And it can be denoted by triangular fuzzy number $W_{u,p}=\{a_{u,p},b_{u,p},c_{u,p}\}$. ${\omega }^p$ is the weight value for decision $p$. $Z_{i,p}=\{A_{i,p},B_{i,p},C_{i,p}\}(t=\mathrm{1,2},\dots ,\left|N\right|,\hspace{0.17em}\hspace{0.17em}p=\mathrm{1,2},\dots ,m)$, where $\left|N\right|$ is the total number of the covering customers. $Z_{i,p}$ represents the comprehensive evaluation index of decision variables $p$ by all decision makers for customer $i$.$$\begin{array}{}\text{(9)}& A_{i,p}=\frac{\mathrm{1}}{r}\otimes {\displaystyle \sum }_{u=\mathrm{1}}^r{a}_{u,p}\text{(10)}& B_{i,p}=\frac{\mathrm{1}}{r}\otimes {\displaystyle \sum }_{u=\mathrm{1}}^r{b}_{u,p}\text{(11)}& C_{i,p}=\frac{\mathrm{1}}{r}\otimes {\displaystyle \sum }_{u=\mathrm{1}}^r{c}_{u,p}\end{array}$$

The comprehensive evaluation index of decision variables $p$ for customer $i$ can be calculated as$$\begin{array}{}\text{(12)}& Z_{i,p}=\frac{\mathrm{1}}{r}\otimes {\displaystyle \sum }_{u=\mathrm{1}}^r{W}_{u,p}\end{array}$$

To simplify, we define the comprehensive membership of decision variable $p$ for customer $i$ as follows:$$\begin{array}{}\text{(13)}& p_i=\frac{\mathrm{1}}{\mathrm{5}}\left(A_{i,p}+\mathrm{2}{B}_{i,p}+\mathrm{2}{C}_{i,p}\right)\end{array}$$

The comprehensive evaluation value of a decision variable $p$ denotes the average value of all customers, such as$$\begin{array}{}\text{(14)}& {\tilde{\omega }}^p=\frac{\mathrm{1}}{N}{\displaystyle \sum }_{i=\mathrm{1}}^N{p}_i,\hspace{1em}\mathrm{f}\mathrm{o}\mathrm{r}\hspace{0.17em}\hspace{0.17em}{\omega }^p>\mathrm{0}\end{array}$$Considering that ${\omega }^p$ should satisfy the constraint ${\sum }_{p=\mathrm{1}}^m{\omega }^p=\mathrm{1}$, so we calculate ${\omega }^p$ as follows:$$\begin{array}{}\text{(15)}& {\omega }_p=\frac{{\tilde{\omega }}^p}{{\sum }_{p=\mathrm{1}}^m{\tilde{\omega }}^p},\hspace{1em}p=\mathrm{1,2},\dots ,m\end{array}$$

Then, the similarity between customer $i$ and customer $j$ can be denoted as follows:$$\begin{array}{}\text{(16)}& S_{ij}={\displaystyle \sum }_{p=\mathrm{1}}^m{\omega }^p{X}_{i,j}^p,\hspace{1em}i,j\in N_{\mathrm{0}}\end{array}$$

Finally, we get the comprehensive similarity between customers which can be represented by a similarity matrix just as follows:$$\begin{array}{}\text{(17)}& F={\left(S_{ij}\right)}_{N\times N}=\left(\begin{array}{cccc}\mathrm{1}& S_{\mathrm{12}}& \cdots & S_{\mathrm{1}N}\\ S_{\mathrm{21}}& \cdots & S_{\mathrm{2}N-\mathrm{1}}& S_{\mathrm{2}N}\\ ⋮& ⋮& \ddots & ⋮\\ S_{N\mathrm{1}}& ⋮& S_{NN-\mathrm{1}}& \mathrm{1}\end{array}\right)\end{array}$$

3.2. Second Phase Model

For the second phase model, we should assign the appropriate route for the customer en-groups. In this part we consider the energy conservation and emission reduction and target to reduce the comprehensive costs. There, we use the comprehensive emission measurement model to calculate the energy consumption and carbon emission [24]. At first, we regard the travel speed as fixed variables and establish a conventional model. On this basis, we consider the speed as time variation variables and reform the model, which is consistent with the actual situation. To solve the problem, we design a genetic algorithm.

To provide a precise statement of this problem, we define the parameters and indices shown in Table 3.

Table 3

Sets, indices, and parameters in the C-VRPESER.

Notations	Detailed Definition
$N$	Where $N=\{\mathrm{0,1},\mathrm{2},\cdots ,n\}$, 0 represents the depot, $N_{\mathrm{0}}=N\setminus \mathrm{0}=\{\mathrm{1,2},\cdots ,n\}$
	is the set of the customers which should be delivered.
$$A$$	$A=\{(i,j):i,j\in N\}$ represents the arcs.
$T_R$	Set of of customers in the same customer group $R$.
$l$	Type of the vehicle.
$l^k$	$k$ th vehicle with $l$ type.
$Q_e^l$	Maximum capacity of each type vehicles
$$i\in N_{\mathrm{0}}$$	Number of customer.
$q_i$	Demand of customers $i$.
$s_i$	Time of vehicle arrive at customer $i$.
$\left[a_i,b_i\right]$	Desired preferred time window for customer $i$ to be serviced.
$t_i$	Service time of customer $i$ need.
$d_{ij}$	Distance between the customer $i$ and $j$.
$v_{ij}$	Speed of the vehicle travel on the link $(i,j)$.
$q_{ij}^{l^k}$	Load of the vehicle $l^k$ travel on the link $(i,j)$.
$$Q_c^l$$	Light weight of the $l$ type vehicle.
$$c_o$$	Per unit cost of fuel.
$$c_e$$	Per unit cost of carbon emissions (carbon tax).
${\delta }_c$	Fuel emission factor.
$c_l$	Fixed costs of the vehicle used.
$e_k$	Departure time of the vehicle k from the depot.
$l_k$	Latest arrival time at the depot after services the customers.
$c_{\mathrm{1}}$	Penalty coefficient for earliness arrival.
$c_{\mathrm{2}}$	Penalty coefficient for delay delivery.
$S_i$	Service time of customer $i$.

For the en-group route optimization, we consider two decision variables as follows:$$\begin{array}{}\text{(18)}& x_{ij}^{l^k}=\left\{\begin{array}{cc}\mathrm{1},\hfill & \mathrm{i}\mathrm{f}\hspace{0.17em}\hspace{0.17em}\mathrm{t}\mathrm{h}\mathrm{e}\hspace{0.17em}\hspace{0.17em}\mathrm{v}\mathrm{e}\mathrm{h}\mathrm{i}\mathrm{c}\mathrm{l}\mathrm{e}\hspace{0.17em}\hspace{0.17em}{l}^k\hspace{0.17em}\hspace{0.17em}\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{v}\mathrm{e}\mathrm{l}\hspace{0.17em}\hspace{0.17em}\mathrm{o}\mathrm{n}\hspace{0.17em}\hspace{0.17em}\mathrm{t}\mathrm{h}\mathrm{e}\hspace{0.17em}\hspace{0.17em}\mathrm{l}\mathrm{i}\mathrm{n}\mathrm{k}\hspace{0.17em}\hspace{0.17em}\left(i,j\right)\hfill \\ \mathrm{0},\hfill & \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\hfill \end{array}\right.\text{(19)}& z_i^{l^k}=\left\{\begin{array}{cc}\mathrm{1},\hfill & \mathrm{i}\mathrm{f}\hspace{0.17em}\hspace{0.17em}\mathrm{t}\mathrm{h}\mathrm{e}\hspace{0.17em}\hspace{0.17em}\mathrm{c}\mathrm{u}\mathrm{s}\mathrm{t}\mathrm{o}\mathrm{m}\mathrm{e}\mathrm{r}\hspace{0.17em}\hspace{0.17em}i\hspace{0.17em}\hspace{0.17em}\mathrm{n}\mathrm{e}\mathrm{e}\mathrm{d}\hspace{0.17em}\hspace{0.17em}\mathrm{t}\mathrm{o}\hspace{0.17em}\hspace{0.17em}\mathrm{b}\mathrm{e}\hspace{0.17em}\hspace{0.17em}\mathrm{s}\mathrm{e}\mathrm{r}\mathrm{v}\mathrm{i}\mathrm{c}\mathrm{e}\mathrm{d}\hspace{0.17em}\hspace{0.17em}\mathrm{b}\mathrm{y}\hspace{0.17em}\hspace{0.17em}\mathrm{t}\mathrm{h}\mathrm{e}\hspace{0.17em}\hspace{0.17em}k\mathrm{t}\mathrm{h}\hspace{0.17em}\hspace{0.17em}\mathrm{v}\mathrm{e}\mathrm{h}\mathrm{i}\mathrm{c}\mathrm{l}\mathrm{e}\hspace{0.17em}\hspace{0.17em}\mathrm{f}\mathrm{o}\mathrm{r}\hspace{0.17em}\hspace{0.17em}\mathrm{t}\mathrm{y}\mathrm{p}\mathrm{e}\hspace{0.17em}\hspace{0.17em}l\hfill \\ \mathrm{0},\hfill & \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\hfill \end{array}\right.\end{array}$$

The model can be formulated as follows:$$\begin{array}{}\text{(20)}& \min u_R=c_{\mathrm{1}}{\displaystyle \sum }_{i\in T_R}\max \left[\left(a_i-s_i\right),\mathrm{0}\right]+c_{\mathrm{2}}{\displaystyle \sum }_{i\in T_R}\mathrm{m}\mathrm{a}\mathrm{x}\left[\left(s_i-b_i\right),\mathrm{0}\right]+{\displaystyle \sum }_{\left(i,j\right)\in T_R}{c}_o{O}_{F,ij}\left(Q_c^l+q_{ij}^{l^k},v_{ij},d_{ij}\right)+{\displaystyle \sum }_{\left(i,j\right)\in T_R}{c}_e{\delta }_c{O}_{F,ij}\left(Q_c^l+q_{ij}^{l^k},v_{ij},d_{ij}\right)+c_l\end{array}$$

s.t.$$\begin{array}{}\text{(21)}& {\displaystyle \sum }_{l^k\in L^K}{\displaystyle \sum }_{j\in T^R}{x}_{\mathrm{0}j}^{l^k}=\mathrm{1}\text{(22)}& {\displaystyle \sum }_{l^k\in L^K}{\displaystyle \sum }_{j\in T^R}{x}_{j\mathrm{0}}^{l^k}=\mathrm{1}\text{(23)}& {\displaystyle \sum }_{l^k\in L^K}{\displaystyle \sum }_{j\in T^R}{x}_{ij}^{l^k}=\mathrm{1},\hspace{1em}\forall i\in T_R\text{(24)}& {\displaystyle \sum }_{l^k\in L^K}{\displaystyle \sum }_{j\in T^R}{x}_{jh}^{l^k}-{\displaystyle \sum }_{l^k\in L^K}{\displaystyle \sum }_{j\in T^R}{x}_{hj}^{l^k}=\mathrm{0},\hspace{1em}\forall h\in T_R\text{(25)}& {\displaystyle \sum }_{i\in T_R}{q}_i{\displaystyle \sum }_{l^k\in L^K}{\displaystyle \sum }_{j\in T^R}{x}_{ij}^{l^k}\le Q_e\text{(26)}& a_i{\displaystyle \sum }_{l^k\in L^K}{\displaystyle \sum }_{j\in T^R}{x}_{ij}^{l^k}\le s_i^{l^k}\le b_i{\displaystyle \sum }_{l^k\in L^K}{\displaystyle \sum }_{j\in T^R}{x}_{ij}^{l^k},\hspace{1em}\forall i\in T_R\text{(27)}& x_{ij}^{l^k}\left(s_i^{l^k}+t_i+\frac{d_{ij}}{v_{ij}}\right)\le s_i^{l^k}\text{(28)}& x_{ij}^{l^k}\in \left\{\mathrm{0,1}\right\}\text{(29)}& y_i^{l^k}\ge \mathrm{0}\end{array}$$

The objective function (20) is to minimize the total routing cost, including fixed cost, energy consumption cost, carbon emission cost, and punishment for time windows. Constraints (21) and (22) guarantee that all service vehicles start from the depot and finally return to the depot for once. Constraint (23) ensures each vehicle travels on the arc $(v_i,v_j)$ no more than once. Constraint (24) is the flow balance and shows that the number of the vehicles which arrive at the customer is the same as those which depart. Constraint (25) ensures the loading customer demands no more than the vehicle capacity. Inequality (26) is the time window constraint. Constraint (27) represents the time limit for vehicle travel horizon. Constraint (28) is the binary variable constraint. Constraint (29) represents the nonnegative of time.

3.3. Model Reformulation

The above model is generated based on the assumption that the drive speed is invariable. However in reality the speed is influenced by a lot of elements. To make the model more practical, we assume that the speed is a ladder type change variable. This paper supposes the vehicles traveling speed as shown in Figure 2. On this basis, we build the time-varying velocity distribution vehicle routing optimization model considering the velocity change.

[figure omitted; refer to PDF]

The speed is assumed as a step function of time as Figure 2. And during different time intervals, the vehicles' travel speed is different. To model convince, the running time is divided into $I$ interval according different extents speed, which can be represented by $T^{\mathrm{1}}$, $T^{\mathrm{2}}$,..., $T^I$. The $m$th time interval is indicated as $[t_m,t_{m+\mathrm{1}}]$; in each time interval velocity is constant.

The time of vehicles running on the path $(i,j)$ is a function of the time $l{t}_i$ that is the vehicles' departing time from the customer $i$. The speed of vehicles running on the path often crosses more than one speed range. It is assumed that the vehicle will cross $p+\mathrm{1}$ time intervals. It is clear that the vehicle has $p$ different speed. Assume the vehicle derived the customer $i$ in the $m$ time interval with the speed $v_{ij}^m(l{t}_i)$, then the vehicle speed turns $v_{ij}^{m+\mathrm{1}}(l{t}_i)$ in the next time interval, and we denote the vehicle speed as $v_{ij}^{m+p}(l{t}_i)$ when the vehicle derives the link $(i,j)$. And the distance for the link $(i,j)$ is denoted by $d_{ij}^{m+p}(l{t}_i)$. And the total time $t{t}_i$ for the vehicle drived across the link $(i,j)$ can be represented by the time function $l{t}_i$ as follows:$$\begin{array}{}\text{(30)}& t{t}_{ij}\left(l{t}_i\right)=\frac{d_{ij}^m}{v_{ij}^m\left(l{t}_i\right)}+\frac{d_{ij}^{m+\mathrm{1}}}{v_{ij}^{m+\mathrm{1}}\left(l{t}_i\right)}+\cdots +\frac{d_{ij}^{m+p}}{v_{ij}^{m+p}\left(l{t}_i\right)}\end{array}$$

Considering the variety of the vehicle traveling speed, we can represent the fuel consumption function as follows:$$\begin{array}{}\text{(31)}& O_{F,ij}\left(Q_c^l+q_{ij}^{l^k},{\overline{v}}_{ij},d_{ij}\right)={\displaystyle \sum }_{q=\mathrm{1}}^I{O}_{F,ij}\left(Q_c^l+q_{ij},v_{ij}^{m+q},d_{ij}^{m+q}\right)\end{array}$$

If we assumed that the $k$th vehicle of the $l$ kind begins the task at time $y_{l^k}$, the objective function can be denoted as follows:$$\begin{array}{}\text{(32)}& \min u_R=c_{\mathrm{1}}{\displaystyle \sum }_{i\in T_R}\max \left[\left(a_i-s_i\right),\mathrm{0}\right]+c_{\mathrm{2}}{\displaystyle \sum }_{i\in T_R}\max \left[\left(s_i-b_i\right),\mathrm{0}\right]+{\displaystyle \sum }_{\left(i,j\right)\in T_R}{c}_o{O}_{F,ij}\left(Q_c^l+q_{ij}^{l^k},{\overline{v}}_{ij},d_{ij}\right)+{\displaystyle \sum }_{\left(i,j\right)\in T_R}{c}_e{\delta }_c{O}_{F,ij}\left(Q_c^l+q_{ij}^{l^k},{\overline{v}}_{ij},d_{ij}\right)+c_l\end{array}$$

The value of ${\overline{v}}_{ij}$ depends on the time interval. The time constraint is represented as follows:$$\begin{array}{}\text{(33)}& x_{ij}^{l^k}\left(s_i^{l^k}+t_i+t{t}_{ij}\left(s_i^{l^k}+t_i\right)\right)\le s_j^{l^k}\end{array}$$

The vehicle service time cannot be earlier than the start time of the distribution center, not later than the time of the delivery center.$$\begin{array}{}\text{(34)}& t_{\mathrm{1}}\le y_{l^k}\le t_I\end{array}$$

4. Algorithm

To present the customers-oriented vehicle routing problem with environment consideration, a two-phase model is formulated. We design two algorithms to solve the problem. First, we design a fuzzy system clustering algorithm to realize the customers grouping. That algorithm is mainly used to consider the customers demand attributes, including quantitative and qualitative attributes. The second phase model is an en-group VRP that need to design a heuristic algorithm to gain the en-group optimal delivery routes. Here, we compare the performance of several heuristic algorithms, as shown in Table 4. In view of the robustness and global optimization of genetic algorithm, we design a genetic algorithm to solve the problem.

Table 4

Heuristic algorithm performance comparison.

Algorithm	Parameters	Features
Tabu search algorithm (TS)	Tabu size Neighborhood structure Number of solutions	Local optimum Higher solving efficiency
Ant colony algorithm (ACO)	Initial pheromone Information increase Information accumulation Disappearance factor Inspiration factor	Good parallelism Global optimization features Time-consuming
Simulated annealing algorithm (SA)	Initial temperature Temperature return function Mode of state generation Sampling principle	Simple implementation Flexible application The algorithm converges slowly
Genetic algorithm (GA)	Population number Selection operator Crossover operator Mutation operator	Strong global optimization ability Parallelism and strong robustness

4.1. Fuzzy System Clustering Algorithm

On the basis of the similarity matrix, we design a fuzzy system clustering algorithm to cluster the customer into different groups, which are often used to process the data with quantitative and qualitative data [22]. And the specific algorithm procedure is summarized as in Algorithm 1.

Algorithm 1: Fuzzy system clustering algorithm.

1 Step 1 : Data initialization, input the fuzzy equivalent matrix $F(k)$. Let iterations $\pi$=1.

2 Step 2 : Find out the largest elements in the matrix $F_s(k{)}^T$, and chose the customer with the smaller number as the target

customer and denote it with $s$.

3 Step 3 : Mark the selected target customer and delete the corresponding row $F_s(k{)}^T$. Specifically, once a customer is chosen as

the target customer, its column is marked by $F_s(k)$, mark the corresponding elements in row $F_s(k{)}^T$ as 0.

4 Step 4 : Find out the largest element of matrix $F_s(k{)}^T$, denoted the similarity value as ${\omega }_{rs}$.

5 $(1)$ If ${\omega }_{rs}\le \lambda$(pre-given threshold value), then go to Step 5. Otherwise, calculate the total weight of the customers $Q$.

6 $(2)$ Make judgement, if $Q\ge Q_e^l$, go to Step 5. Otherwise, put customers $s$ and $r$ in same group, then delete the corresponding

row $F_r(k{)}^T$ of matrix $F(k)$, return to Step 4.

7 Step 5: Set the customer clustering termination conditions; if the matrix is empty, then clustering process stops; else let

$\pi$ = $\pi$ + 1 and got to Step 2.

To obtain the cluster solutions, we apply the procedure in Figure 3.

[figure omitted; refer to PDF]

4.2. En-Group Genetic Algorithm

The proposed model is a mixed-integer linear programming model. When the scale of the problem is expended, it is difficult to solve the problem with accurate algorithm. However, the GA is suitable to solve multivariables complex problem with multiple parameters. It is widely applied in solving VRP due to its extensive generality high efficiency and strong robustness. Therefore we deign GA to solve the problem. The key steps are designed as follows.

Encoding. Natural number cod approach is used to characterize chromosomes. Chromosome genes' encoding is based on the results of the customer groups with discontinuous customer numbers. There, 0 expresses the distribution center, and customers are expressed by numbers 1, 2, 3, $\cdots$, $N$. Different groups are corresponding different chromosomal genes. The beginning and the end of the genes in the chromosome cod sequence must be "0", on behalf of the vehicle which departs from the distribution center and finally returns to the distribution center after completing the delivery tasks.

Initial Solution. Construct initial population chromosome cod using MATLAB programming software under the premise of customer grouping. Specific we randomly give the initial solution as $pop(i,\mathrm{2}:n-\mathrm{1})=KHz(\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}(numel(KHz)))$,$pop(i,\mathrm{1})=\mathrm{0}$, and $pop(i,n)=\mathrm{0}$, such as $KHz(\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}(numel(KHz)))$=$[\mathrm{0,1},\mathrm{4,5},\mathrm{9,12,18,0}]$. And generate the initial population as a matrix with the dimension $popsize\ast (n+\mathrm{2})$.

Fitness Function. In this paper, the objective function value is used to measure the fitness value as follows:$$\begin{array}{}\text{(35)}& f_k=\mathrm{1}-\frac{\left({\mu }_k-\min \mu \right)}{\left(\max \mu -\min \mu \right)}\end{array}$$

Selection Operator. It is individual selection according to roulette wheel; if the fitness of an individual is denoted by $f_k$, then the probability of being selected can be represented by the following:$$\begin{array}{}\text{(36)}& p_k=\frac{f_k}{{\sum }_{k=\mathrm{1}}^K{f}_k}\end{array}$$

Crossover Operator. In order to improve the convergence speed of the algorithm, we chose adaptive crossover probability as follows:$$\begin{array}{}\text{(37)}& P_c=\left\{\begin{array}{cc}p-k\times \frac{f_{max}-f}{f_{max}-\overline{f}}\hfill & \left(f\ge \overline{f}\right)\hfill \\ p\hfill & \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\hfill \end{array}\right.\end{array}$$

The values of the parameters $p$ are between the values $\mathrm{0.8}~\mathrm{1}$, and the value of the parameter $k$ is between the value of 0.2$~$0.4. In this paper, a partial matching crossover strategy is selected, which means that the genes of two crossing points are interchanged by “two points crossing", but the first “0" does not participate in the intersection, as described in Figure 4.

[figure omitted; refer to PDF]

Mutation Operator. This paper uses the basic bit variation of random exchange, which means random exchange of the genes encoding of chromosomes to produce new chromosomes.

Termination. In this paper, the termination criterion is considered as only if the maximum number of iterations is reached, end the algorithm and output the optimal solution.

Based on the aforementioned analysis, the procedure of this heuristic algorithm is summarized By algorithm represented as follows.

Step 1.

Obtain the customer group results and determine the customer number corresponding to the genetic code.

Step 2.

Set the number of initial populations $n$, adaptive crossover probability $P_c$, mutation probability $P_m$, and maximum number of iterations $N$.

Step 3.

Randomly generate the initial solutions $X$, and let iter $m$=1.

Step 4.

By selecting, crossing, and mutating, generate new strategies $X_m$, calculate the objective, and let $m=m+\mathrm{1}$.

Step 5.

If the algorithm terminates the condition (e.g. reaches the set maximum iteration number), the algorithm is terminated, otherwise returned to Step 4.

To obtain the optimal solutions, we applied the procedure illustrated in Figure 5.

[figure omitted; refer to PDF]

5. Case Study

In this section, in order to verify the effectiveness of the proposed model and algorithms, we conducted experiments with real data. The proposed cluster algorithm and GA algorithm are coded in the MATLAB 2012a. All numerical tests are carried out on a Windows Server 2007 with Intel CPU 64 G RAM to improve the solving efficiency.

5.1. Experiment Description and Parameter Settings

The experimenting setting provides the experimental results and analysis which are shown below. As we know, the distribution center is equipped with several Dong Feng vehicles with the rated load for 2 tons and 3 tons and the parameters are shown in Table 5. There are 24 customers that should be serviced whose demands can be known in advance, and the customers are supermarkets and the markets in the distribution area. The demand information and location of the customer are shown in Table 6. To calculate the cost, we assume that the driving speed of the vehicle is 40km/h for the fixed speed, and the distance between each customer and distribution center is calculated according to their coordinates. The fuel price is 5.55 Yuan/L, and the cost of carbon emission is 51.52 Yuan/T. The fixed cost of using a vehicle is estimated at 100 RMB per visit.

Table 5

The parameters of the covered vehicles.

parameters	2T vehicle	3T vehicle
Gross vehicle weight $w_z$	2T	2.5T
Vehicle emissions $V_l$	3.2	3.8
Type of fuel $f$	diesel	diesel
Load capacity $w_e$	2T	3T
The front surface area of a vehicle $S_c$	2.4$m^{\mathrm{2}}$	3$m^{\mathrm{2}}$

Table 6

The demand information of the covered customers.

Number	Demand	Distance	External similarity	Product value	Time window	Quality of Service
1	0.5	(-10.3, 9.7)	high	medium	$[$ 8:00, 9:00$]$	high
2	0.2	(-4.4,-6.6)	low	low	$[$ 10:00, 12:00$]$	lowest
3	0.6	(-8.5, 6.3)	medium	high	$[$ 12:00, 12:30$]$	high
4	0.4	(2.7, 9.3)	low	medium	$[$ 10:00, 11:00$]$	medium
5	0.3	(-5.3, -4.5)	lowest	high	$[$ 14:00, 14:30$]$	low
6	0.1	(-0.6, 5.2)	low	low	$[$ 8:00, 9:30$]$	medium
7	0.5	(18.6, -13.2)	high	highest	$[$ 11:00, 13:00$]$	highest
8	0.3	(-11.9, 3.7)	highest	medium	$[$ 15:00, 16:30$]$	medium
9	0.1	(3.3, 7.8)	medium	low	$[$ 14:00, 15:30$]$	high
10	0.1	(-1.8, -11.1)	medium	medium	$[$ 13:00, 15:00$]$	medium
11	0.2	(8.6, 2.5)	lowest	lowest	$[$ 14:00, 14:30$]$	low
12	0.3	(-6.9, 1.6)	medium	high	$[$ 14:00, 15:00$]$	high
13	0.4	(7.8, -7)	medium	highest	$[$ 9:00, 11:00$]$	medium
14	0.2	(16.8, -18.4)	low	high	$[$ 15:00, 16:00$]$	high
15	0.1	(12.7, -1.8)	high	medium	$[$ 11:00,11.30$]$	highest
16	0.6	(14.3, -21.3)	medium	medium	$[$ 8:00, 11:00$]$	medium
17	0.2	(8.2, -15.8)	low	lowest	$[$ 14:00, 17:00$]$	lowest
18	0.3	(11.5, 10.4)	high	medium	$[$ 13:00, 14:30$]$	high
19	0.1	(1.1, -8.5)	high	highest	$[$ 10:00, 15:00$]$	high
20	0.4	(20.7, -22.7)	medium	low	$[$ 11:00, 12.30$]$	high
21	0.3	(6.4, -0.9)	low	medium	$[$ 8:00, 9:30$]$	low
22	0.2	(-7.3, 11.6)	highest	high	$[$ 13:00, 14:30$]$	medium
23	0.4	(-12.2, -2.4)	medium	low	$[$ 10:00, 12:00$]$	medium
24	0.5	(10.4, 21.7)	high	high	$[$ 12:00, 13:00$]$	highest

5.2. Customers Group

For simplicity, we choose five parameters as the main variables to cluster the customers into different groups, in which include geographical position, pressed for time, physical properties, product value, and service qualities. Furthermore, we use the customer demand, service time windows, and the latest delivery time as the input data to decide the service authorities and generate the delivery sequence. The corresponding fuzzy evaluation values of the five parameters are shown in Table 7.

Table 7

The qualitative data of triangle fuzzy evaluation.

Customer number	Product external similarity	Product value	Service quality
1	(0.5, 0.75, 1)	(0.25, 0.5, 0.75)	(0.5,0.75,1)
2	(0, 0.25, 0.5)	(0, 0.25, 0.5)	(0, 0, 0.25)
3	(0.25, 0.5, 0.75)	(0.5, 0.75, 1)	(0.5, 0.75, 1)
4	(0, 0.25, 0.5)	(0.25, 0.5, 0.75)	(0.25, 0.5, 0.75)
5	(0, 0, 0.25)	(0.5, 0.75, 1)	(0, 0.25, 0.5)
6	(0, 0.25, 0.5)	(0, 0.25, 0.5)	(0.25, 0.5, 0.75)
7	(0.5, 0.75, 1)	(0.75, 1, 1)	(0.75, 1, 1)
8	(0.75,1,1)	(0.25,0.5,0.75)	(0,0.25,0.5)
9	(0.25, 0.5, 0.75)	(0, 0.25, 0.5)	(0.5, 0.75, 1)
10	(0.25,0.5,0.75)	(0.25,0.5,0.75)	(0.25,0.5,0.75)
11	(0, 0, 0.25)	(0, 0, 0.25)	(0, 0.25, 0.5)
12	(0.25, 0.5, 0.75)	(0.5, 0.75, 1)	(0.5,0.75,1)
13	(0.25, 0.5, 0.75)	(0.75, 1, 1)	(0.25, 0.5, 0.75)
14	(0.5, 0.75, 1)	(0.5, 0.75, 1)	(0.5, 0.75, 1)
15	(0.5, 0.75, 1)	(0.25, 0.5, 0.75)	(0.75, 1, 1)
16	(0.25,0.5,0.75)	(0.25,0.5,0.75)	(0.25,0.5,0.75)
17	(0, 0.25, 0.5)	(0, 0, 0.25)	(0, 0, 0.25)
18	(0.5,0.75,1)	(0.25, 0.5, 0.75)	(0.5,0.75,1)
19	(0.5, 0.75, 1)	(0.75, 1, 1)	(0.5,0.75,1)
20	(0.25, 0.5, 0.75)	(0, 0.25, 0.5)	(0.75, 1, 1)
21	(0, 0.25, 0.5)	(0.25, 0.5, 0.75)	(0, 0.25, 0.5)
22	(0.75, 1, 1)	(0.5, 0.75, 1)	(0, 0, 0.25)
23	(0.25, 0.5, 0.75)	(0, 0.25, 0.5)	(0.25,0.5,0.75)
24	(0.5, 0.75, 1)	(0.5, 0.75, 1)	(0.5, 0.75, 1)

In this paper, different types of vehicles are used for distribution services, which will lead to the fact that group results are different, and thus different vehicle scheduling schemes are available, just as shown in Tables 8 and 9. Take consideration of the two tables, we can point out that the vehicle capacity influences the cluster results. The overall utilization rates of the two programs are 91.25%. However, the two programs produce different vehicle fixed costs and different resources. In contrast to the two programs, taking into account the fixed costs of labor costs, we chose the first program as the actual cluster result, in which intragroup path can be obtained. And we use Figure 6 as the demonstration.

Table 8

The clustering results of vehicle restrictions for 2 tons.

Cluster condition	Group No.	Customer No.	Cargo volume (tons)
$Q_e$  $\lambda$=2t	Group 1	6,21,2,4,17,14,8,22	1.9
	Group 2	7,24,15;,1,19,5	2.0
	Group 3	9,12,10,3,13,23	1.9
	Group 4	11,18,16,20	1.5

Table 9

The clustering results of vehicle restrictions for 3 tons.

Cluster condition	Group No.	Customer No.	Cargo volume (tons)
$Q_e$  $\lambda$=3t	Group 1	6,21,2,4,17,14,8,22,5,11	2.4
	Group 2	7,24,15,1,19,18	2.0
	Group 3	9,12,10,3,13,23,16,20	1.9

[figure omitted; refer to PDF]

5.3. Delivery Vehicle Path Plan

Combining the customer grouping results in Section 5.1, we use the designed genetic algorithm to solve the intragroup delivery vehicle routing model. Considering that the grouped customers only need one delivery vehicle to complete the delivery service, the en-group path solution was converted to the TSP problem. Set the initial population to 100, the mutation probability is 0.09, and the maximum number of iterations was set to 200. Finally, we get the approximate optimal delivery vehicle path plan with a short time consumption. The convergence of the algorithm is shown in Figure 7.

[figure omitted; refer to PDF]

Considering the energy cost, carbon emission cost, time penalty cost, and fixed cost in the transportation process as objective functions, the distribution vehicle route within the group is optimized. As a heuristic algorithm, the genetic algorithm is generally not an optimal solution but an approximately optimal solution. Therefore, it needs to execute the program multiple times to obtain a better approximately optimal solution. Finally we get the optimal solution shown in Table 10. To show the specific vehicle routing, we give out the delivery path of group 1, as shown in Figure 8.

Table 10

The optimization paths results with customer grouping.

Group No.	The optimal path en-group	Distance(Km)	Cost
Group 1	$\mathrm{0}\to \mathrm{21}\to \mathrm{6}\to \mathrm{4}\to \mathrm{22}\to \mathrm{8}\to \mathrm{5}\to \mathrm{2}\to \mathrm{17}\to \mathrm{14}\to \mathrm{11}\to \mathrm{0}$	109.22	403.3
Group 2	$\mathrm{0}\to \mathrm{1}\to \mathrm{19}\to \mathrm{7}\to \mathrm{15}\to \mathrm{24}\to \mathrm{18}\to \mathrm{0}$	104.56	329.96
Group 3	$\mathrm{0}\to \mathrm{13}\to \mathrm{16}\to \mathrm{20}\to \mathrm{10}\to \mathrm{23}\to \mathrm{3}\to \mathrm{12}\to \mathrm{9}\to \mathrm{0}$	106.44	359.11

[figure omitted; refer to PDF]

5.4. Group Optimization Effect Analysis

The proposed models consider the customer demand and group the customers on the basis of considering the diversity of customer needs. However, the traditional two-phase vehicle path optimization model is generally based on the customer’s geographical location and time window and does not fully consider the customer’s other demand attributes. Compared with other methods, the customer-oriented distribution vehicle routing model proposed in this paper has two advantages. For one hand, it can improve the efficiency of problem solving. For the other hand, customers with similar demands are divided into same groups, so that the distribution process can provide more targeted services. In order to prove the optimization effect of the grouping method, this paper solves the case by designing the genetic algorithm without grouping and compares the solution results with the calculation results under the grouping situation. In the case of without grouping, the optimal path of delivery vehicles is shown in Table 11.

Table 11

The optimization paths results with no grouping.

Group No.	The optimal path en-group	Distribution quality(tons)
Path 1	$\mathrm{0}\to \mathrm{3}\to \mathrm{1}\to \mathrm{22}\to \mathrm{6}\to \mathrm{9}\to \mathrm{4}\to \mathrm{18}\to \mathrm{11}\to \mathrm{21}\to \mathrm{0}$	2.7
Path 2	$\mathrm{0}\to \mathrm{13}\to \mathrm{19}\to \mathrm{16}\to \mathrm{14}\to \mathrm{7}\to \mathrm{17}\to \mathrm{10}\to \mathrm{2}\to \mathrm{5}\to \mathrm{20}\to \mathrm{0}$	3
path 3	$\mathrm{0}\to \mathrm{12}\to \mathrm{8}\to \mathrm{23}\to \mathrm{24}\to \mathrm{15}\to \mathrm{0}$	1.6

(1) Improve Solution Efficiency. Without grouping, we should calculate the number of vehicles used in the task. In this paper, $m=[\sum {\phi }_i^{\mathrm{1}}/Q_e]+\mathrm{1}$ is a generally chosen and the number of vehicles is given out 3 of 3t Dongfeng vehicle. The genetic algorithm coding method adopts the natural number coding method and finally determine the chromosome of length 28 as $(\mathrm{0,2},\mathrm{4,13},\dots \mathrm{0,7}\dots ,\mathrm{0,16,23},\dots \mathrm{0})$. The customers demand between two zeros on the chromosome cannot exceed the vehicle load. We can point out that the chromosome coding length is longer than en-group, which makes the cross-variation complexity, and easy to fall into the local optimum. However, the position of 0 needs to be adjusted after each iteration to ensure vehicle loading constraints. The complexity of the algorithm is high and the solution is difficult. Through the solution, it is found that the convergence effect of the algorithm is less than the en-group algorithm, and it is easy to fall into the local optimum. The calculation results are quite different; the solution time without grouping is about 537s; however the total solution time after grouping is about 300s. It easy to point out that the computational efficiency is improved by 79%.

(2) Improve the Quality of Delivery Services. In order to measure the optimization effect of customer grouping on distribution service. We calculate the similarity standard deviation of customer demand attribute values within each group, which are regarded as the comprehensive satisfaction evaluation value of this group. The smaller the standard deviation, the greater the similarity of customer demand, which helps to arrange targeted delivery services and improve service quality. This paper takes the parameter service quality as an example to calculate the standard deviation of the evaluation value. Calculate the similarity standard deviation for two conditions and compare the results. The calculation result is shown in Figure 9.

[figure omitted; refer to PDF]

According to Figure 9, it can be seen that the customer has a high similarity to the quality of service requirements on each of the distribution vehicle routes based on the customer grouping, and the variance is small. On the contrary, the solution without grouping has lower similarity in service quality requirements on each path.

5.5. Analysis of Optimization Effect of Energy Saving and Emission Reduction

The energy consumption and carbon emissions are affected by many factors; the main influencing factors are path length, vehicle load, and vehicle speed. Existing research mainly aims to reduce transportation costs by shortening the length of vehicle path. The optimization goal set in this paper is no longer simple considering the shortest path, but comprehensively considering reducing energy consumption and carbon emission costs and achieving comprehensive optimization goals based on energy conservation and emission reduction. In order to prove the validity of the model, the classic shortest path model is selected, and the optimization effects of the two types of models are analyzed through case comparison.

The optimal path of each group is solved by the shortest path model and the minimum energy consumption and carbon emission model, respectively. The change of objective function causes the change of fitness function. On this basis, we obtain the optimal path for each evaluation mechanism. Furthermore, the corresponding path length, energy consumption, and carbon emission are calculated, respectively. To simplify, we denote the minimum energy consumption and carbon emission path model as ECCM, and the shortest path model is remarked as SPM. Calculate the two models separately and get their corresponding distribution vehicle path plan, as shown in Tables 12 and 13.

Table 12

The data results of the ECCM.

Group No.	The optimal path en-group	Distance (Km)	Energy consumption (L)	Carbon emission (Kg)
Group 1	$\mathrm{0}\to \mathrm{21}\to \mathrm{11}\to \mathrm{4}\to \mathrm{6}\to \mathrm{22}\to \mathrm{8}\to \mathrm{5}\to \mathrm{2}\to \mathrm{17}\to \mathrm{14}\to \mathrm{0}$	107.04	14.64	39.36
Group 2	$\mathrm{0}\to \mathrm{1}\to \mathrm{24}\to \mathrm{18}\to \mathrm{13}\to \mathrm{7}\to \mathrm{19}\to \mathrm{0}$	101.21	12.53	34.21
Group 3	$\mathrm{0}\to \mathrm{13}\to \mathrm{20}\to \mathrm{16}\to \mathrm{10}\to \mathrm{23}\to \mathrm{3}\to \mathrm{12}\to \mathrm{9}\to \mathrm{0}$	105.76	14.47	39.5

Table 13

The data results of the SPM.

Group No.	The optimal path en-group	Distance(Km)	Energy consumption(L)	Carbon emission(Kg)
Group 1	$\mathrm{0}\to \mathrm{21}\to \mathrm{11}\to \mathrm{14}\to \mathrm{17}\to \mathrm{2}\to \mathrm{5}\to \mathrm{8}\to \mathrm{22}\to \mathrm{6}\to \mathrm{6}\to \mathrm{0}$	103.26	16.33	44.56
Group 2	$\mathrm{0}\to \mathrm{1}\to \mathrm{24}\to \mathrm{18}\to \mathrm{13}\to \mathrm{7}\to \mathrm{19}\to \mathrm{0}$	101.21	12.53	34.21
Group 3	$\mathrm{0}\to \mathrm{13}\to \mathrm{20}\to \mathrm{16}\to \mathrm{10}\to \mathrm{23}\to \mathrm{12}\to \mathrm{3}\to \mathrm{9}\to \mathrm{0}$	101.94	15.49	42.30

It can be seen that, in the distribution vehicle path scheme calculated by the two models, only the optimal paths of the second group are the same; the other two groups are different. This shows that the shortest path of the delivery vehicle is not necessarily the optimal route of the delivery vehicle corresponding to the minimum energy consumption and carbon emissions. Taking the first group as an example, the demand of customer 4 is much larger than other customers. In order to reduce the overall cost, the vehicle will serve the customer 4 as early as possible if the path length does not change much. Finally, we get the total energy consumption of the vehicle distribution calculated by the optimization model ECCM which is 41.64L, and the carbon emission is 113.67kg. Compared with the optimization model SPM, the energy consumption and carbon emissions are reduced by 6.5%.

5.6. Sensitivity Analysis of Velocity Variations

The impact of speed changes on vehicle energy consumption and carbon emissions was studied by analyzing the impact of peak hours on road vehicle speeds. According to the investigation, the distribution area generally has traffic congestion from 8:00 to 9:00 am, 12:00 to 14:00, and 17:00 to 19:00; during these periods the speed of vehicle is lower. Assume that the speed of the vehicle varies in steps over one day, as shown in Table 14.

Table 14

The relational table of velocity and time.

Period	$[$ 0:00,6:00$]$	$[$ 6:00,9:00$]$	$[$ 9:00,12:00$]$	$[$ 12:00,14:00$]$	$[$ 14:00,17:00$]$	$[$ 17:00,19:00$]$	$[$ 19:00,22:00$]$
Speed(Km/h)	50	30	40	35	40	30	40

The distribution activity starts at 8:00 am; taking the optimization of the delivery vehicle path of customer group 1 as an example, the optimal path is obtained, as shown in Figure 10. And we also give out the contrast of the three groups' optimal path with and without speed variation, just as Table 15. Affected by the speed change, it may cause the variation of the delivery order. For example, the delivery path of the customer group 3 is changed. This is because the vehicle will try to choose the distribution vehicle path with the lowest total cost. If the vehicle travels according to the original route, the cost is 373.9 Yuan, which is higher than the new delivery vehicle path. In this example, the total cost of distribution considering time-varying speed is 1160.36 Yuan, which is 6.3% higher than the cost of 1091.71 Yuan by fixed speed, and the energy consumption is increased by 5.4%. This shows that the speed change has a certain impact on the distribution cost and fuel consumption, and it needs to be considered when designing the distribution line.

Table 15

The data results of the SM.

Group No.	The optimal path en-group	Whether the path direction changes	Total cost(Y)	Carbon emission(Kg)
Group 1	$\mathrm{0}\to \mathrm{21}\to \mathrm{6}\to \mathrm{4}\to \mathrm{22}\to \mathrm{8}\to \mathrm{5}\to \mathrm{2}\to \mathrm{17}\to \mathrm{14}\to \mathrm{11}\to \mathrm{0}$	NO	445.84	15.89
Group 2	$\mathrm{0}\to \mathrm{1}\to \mathrm{19}\to \mathrm{7}\to \mathrm{15}\to \mathrm{24}\to \mathrm{18}\to \mathrm{0}$	NO	343.52	13.04
Group 3	$\mathrm{0}\to \mathrm{13}\to \mathrm{16}\to \mathrm{20}\to \mathrm{10}\to \mathrm{23}\to \mathrm{3}\to \mathrm{12}\to \mathrm{9}\to \mathrm{0}$	Change	371.00	14.95

[figure omitted; refer to PDF]

6. Conclusion and Further Study

The fastest consumer demand growth requires an efficient distribution planning. Nowadays, special attention is given to routes in which fuel consumption and emissions can be reduced. The manuscript proposes a two-phase approach to handle vehicle routing problem with time windows and multivehicles coupled with environmental concerns and customer demand. First phase is concerned to establish costumer clusters and second phase is devoted to formulate the en-group mathematical formulation of the problem and develop a genetic algorithm to account for vehicle routing optimization within each group so that fuel consumption and emissions are minimized. The comprehensive emission measurement model from Barth and Scora (2004) is used to calculate the energy consumption and carbon emission. The algorithm accounts for similarity between customers which is represented in a fuzzy similarity matrix; on this basis, the customers can be cluster in appropriate groups. Considering the energy, carbon emission, time penalty, and fixed costs in the transport process incorporated in the objective function, the distribution vehicle route en-group is optimized. Numerical results are presented for real data from Shandong province. Sensitivity analysis of velocity variations is also presented, which verified that the proposed two-phase approach could improve the service quality and reduce the fuel consumption and carbon emission.

Further studies can focus on the following aspects. $(1)$ Electric logistics vehicles are not considered in this paper, which will be meaningful for the green vehicle routing problem. The electric vehicle is a critical topic for further study. $(2)$ This study designs a GA algorithm to generate the near-optimal solution; in the future, designing more effective algorithms for larger scale case will be a direction.

Conflicts of Interest

The authors declare that they have no conflicts of interest regarding the publication of this paper.

Acknowledgments

The work described in this paper was jointly supported by the National Natural Science Foundation of China (no. 71571011), the Fundamental Research Funds for Central Universities (no. 2018YJS191) and the State Key Laboratory of Rail Traffic Control and Safety, Beijing Jiaotong University (no. RCS2017ZJ001)

References

[1] G. B. Dantzig, J. H. Ramser, "The truck dispatching problem," Management Science, vol. 6 no. 1, pp. 80-91, DOI: 10.1287/mnsc.6.1.80, 1959.

[2] C. Çetinkaya, I. Karaoglan, H. Gökçen, "Two-stage vehicle routing problem with arc time windows: a mixed integer programming formulation and a heuristic approach," European Journal of Operational Research, vol. 230 no. 3, pp. 539-550, DOI: 10.1016/j.ejor.2013.05.001, 2013.

[3] F. Miguel, "Emissions minimization vehicle routing problem," The 89th Annual Meeting of the Transportation Research Board, pp. 224-231, 2010.

[4] Y. Huang, L. Zhao, T. Van Woensel, J.-P. Gross, "Time-dependent vehicle routing problem with path flexibility," Transportation Research Part B: Methodological, vol. 95, pp. 169-195, DOI: 10.1016/j.trb.2016.10.013, 2017.

[5] J. F. Ehmke, A. M. Campbell, B. W. Thomas, "Vehicle routing to minimize time-dependent emissions in urban areas," European Journal of Operational Research, vol. 251 no. 2, pp. 478-494, DOI: 10.1016/j.ejor.2015.11.034, 2016.

[6] B. Golden, A. Assad, L. Levy, F. Gheysens, "The fleet size and mix vehicle routing problem," Computers & Operations Research, vol. 11 no. 1, pp. 49-66, DOI: 10.1016/0305-0548(84)90007-8, 1984.

[7] J. Brandão, "A tabu search algorithm for the heterogeneous fixed fleet vehicle routing problem," Computers & Operations Research, vol. 38 no. 1, pp. 140-151, DOI: 10.1016/j.cor.2010.04.008, 2011.

[8] B. Ornbuki, M. Nakamura, M. Osamu, "A hybrid search based on genetic algorithm and tabu search for vehicle routing," Proceedings of the 6th International Conference on Artificial Intelligence and Soft Computing, vol. 7, pp. 176-181, .

[9] J. Cordeau, M. Gendreau, G. Laporte, "A tabu search heuristic for periodic and multi-depot vehicle routing problems," Networks, vol. 30 no. 2, pp. 105-119, DOI: 10.1002/(SICI)1097-0037(199709)30:2<105::AID-NET5>3.0.CO;2-G, 1997.

[10] G. Y. Tütüncü, "An interactive GRAMPS algorithm for the heterogeneous fixed fleet vehicle routing problem with and without backhauls," European Journal of Operational Research, vol. 201 no. 2, pp. 593-600, DOI: 10.1016/j.ejor.2009.03.044, 2010.

[11] E. Demir, T. Bektaş, G. Laporte, "An adaptive large neighborhood search heuristic for the pollution-routing problem," European Journal of Operational Research, vol. 223 no. 2, pp. 346-359, DOI: 10.1016/j.ejor.2012.06.044, 2012.

[12] T. Bektaş, G. Laporte, "The pollution-routing problem," Transportation Research Part B: Methodological, vol. 45 no. 8, pp. 1232-1250, DOI: 10.1016/j.trb.2011.02.004, 2011.

[13] A. Palmer, The development of an integrated routing and carbon dioxide emissions model for goods vehicles, 2007.

[14] E. Demir, T. Bektaş, G. Laporte, "The bi-objective pollution-routing problem," European Journal of Operational Research, vol. 232 no. 3, pp. 464-478, DOI: 10.1016/j.ejor.2013.08.002, 2014.

[15] Y.-J. Kwon, Y.-J. Choi, D.-H. Lee, "Heterogeneous fixed fleet vehicle routing considering carbon emission," Transportation Research Part D: Transport and Environment, vol. 23 no. 8, pp. 81-89, DOI: 10.1016/j.trd.2013.04.001, 2013.

[16] J. Qian, Fuel emission optimization in vehicle routing problems with time-varying speeds, 2011.

[17] J. Qian, R. Eglese, "Fuel emissions optimization in vehicle routing problems with time-varying speeds," European Journal of Operational Research, vol. 248 no. 3, pp. 840-848, DOI: 10.1016/j.ejor.2015.09.009, 2016.

[18] M. W. P. Savelsbergh, "Local search in routing problems with time windows," Annals of Operations Research, vol. 4 no. 1–4, pp. 285-305, DOI: 10.1007/BF02022044, 1985.

[19] O. Bräysy, M. Gendreau, "Vehicle routing problem with time windows, part I: route construction and local search algorithms," Transportation Science, vol. 39 no. 1, pp. 104-118, DOI: 10.1287/trsc.1030.0056, 2005.

[20] O. B. Gendreau, "Vehicle routing problem with time windows, Part II: metaheuristics," Transportation Science, vol. 39 no. 1, pp. 119-139, DOI: 10.1287/trsc.1030.0057, 2005.

[21] H. I. Calvete, C. Galé, M.-J. Oliveros, B. Sánchez-Valverde, "A goal programming approach to vehicle routing problems with soft time windows," European Journal of Operational Research, vol. 177 no. 3, pp. 1720-1733, DOI: 10.1016/j.ejor.2005.10.010, 2007.

[22] J. B. Sheu, "A hybrid fuzzy-optimization approach to customer grouping-based logistics distribution operatio," Applied Mathematical Modelling, vol. 31 no. 2, pp. 1048-1066, 2007.

[23] T. Hu, J. B. Sheu, "A fuzzy-based customer classification method for demand-responsive logistical distribution operations," Fuzzy Sets and Systems, vol. 139 no. 2, pp. 431-450, DOI: 10.1016/S0165-0114(02)00516-X, 2003.

[24] M. Barth, G. Scora, T. Younglove, "Modal emissions model for heavy-duty diesel vehicles," Transportation Research Record, vol. no. 1880, pp. 10-20, 2004.