Vehicle routing problem based on a fuzzy customer clustering approach for logistics network optimization

Abstract

Vehicle routing problem (VRP) is a combinatorial optimization and integer programming problem seeking to service a number of customers with a fleet of vehicles. Customer characteristics are neglected in traditional VRPs in the past due to the heterogeneity and ambiguousness. This study presents a vehicle route optimization model in consideration of customer characteristics with three major components: (1) A hierarchical analysis structure is developed to convert customers’ characteristics into linguistic variables, and fuzzy integration method is used to map the sub-criteria into higher hierarchical criteria based on the trapezoidal fuzzy numbers; (2) A fuzzy clustering algorithm based on Axiomatic Fuzzy Set is proposed to group the customers into multiple clusters; (3) The fuzzy technique for order preference by similarity to ideal solution (TOPSIS) approach is integrated into the dynamic programming approach to optimize vehicle routes in each cluster. A numerical case study in Anshun, China demonstrates the advantages of the proposed method by comparing with the other two prevailing algorithms. In addition, a sensitivity analysis is conducted to capture the impacts of various evaluation criteria weights. The results indicate our approach performs very well to identify similar customer groups and incorporate individual customer’s service priority into VRP.

Keywords

Vehicle routing problem axiomatic fuzzy set fuzzy clustering algorithm dynamic programming approach technique for order preference by similarity to ideal solution (TOPSIS)

1 Introduction

Proposed by Dantzig and Ramser [1], Vehicle routing problem (VRP) has been extensively studied and widely applied to reduce the transportation cost for a large number of logistics companies in reality [2 –4]. The classic VRP is a combinatorial optimization procedure seeking the optimal route to service a number of different customers in the network, where a fleet of vehicles start their trips from the depot and return to the depot after serving all customers with minimized cost manner [5 –7]. Vehicle routing optimization is considered one critical countermeasure to reduce cost and improve service quality for logistics operators [8, 9].

In a large-scale logistics network, a successful implementation of VRP should include both customer clustering and vehicle routing optimization [2 , 10]. However, in the past two decades, the characteristics of customers and products were rarely involved in customer clustering optimization procedure of traditional vehicle routing problem. Recently, the rise of reverse logistics and cold chain logistics study has promoted the depth of the study on clustering customers and products with similar characteristics [11 –13]. Customer clustering can not only improve the logistics system efficiency, but also reduce the operational costs. The customers with similar characteristics can be clustered into several groups where customers share certain common features (i.e. geospatial location, demand etc.), and then, the vehicle routing optimization can be undertaken within each zone.

Each customer’s characteristics can be represented by surrounding transportation condition, geographical environment, demand requirement, goods compatibility etc. However, a certain number of these characteristics cannot be directly measured numerically. To partition delivery regions based on heterogeneous customer characteristics, it is necessary to transform those ambiguous features as inputs for vehicle routing optimization, and then integrate with transportation cost to find the optimal routes within each clustered area. Rather than solving a specific VRP, this study aims to develop a procedure integrating both pre-route customer clustering and en-route VRP with individual customer’s requirement.

The remaining of this paper is organized as follows: The relevant studies on VRP and its variations are initially reviewed, and then, a hierarchical analysis structure for vehicle routing optimization is established to convert customers’ characteristics based on the fuzzy comprehensive evaluation theory. Multiple service distribution regions can be generated by a proposed fuzzy clustering algorithm. Within each clustered region, VRP can be solved using the fuzzy TOPSIS approach and dynamic programming approach. To evaluate the effectiveness of the proposed approach, a case study of vehicle routing problem in Anshun, China is presented, followed by a thorough comparison with different prevailing approaches. In addition, a sensitivity analysis is conducted in capturing the impacts of different evaluation criteria weights. Finally, conclusions and future works are discussed at the end of this paper.

2 Literature review

It is well known that the Vehicle Routing Problem (VRP) is a NP-hard problem, to tackle the VRP in an efficient manner, a large number of study methods have been conducted in the past [14 –16]. Novoa and Storer [17] proposed an approximate dynamic programming approach for the single-vehicle routing problem with stochastic demand from the reoptimization viewpoint. Fuellerer et al. [18] studied the two-dimensional loading vehicle routing problem, which considered the loading rate of each delivery vehicle to satisfy each customer’s demand, and then, an ant colony optimization (ACO) algorithm was presented to combine two different heuristic measures for satisfactory computational results. Following the previous research, Fuellerer et al. [19] further studied the three-dimensional loading capacitated vehicle routing problem, and utilized the ACO algorithm to resolve the problem. Marinakis et al. [20] presented a hybrid particle swarm optimization algorithm to solve the vehicle routing problem and produced very satisfactory results within short computational time. Juan et al. [21] proposed a flexible solution methodology to study the vehicle routing problem with stochastic demands, and the metaheuristics were used to assign different values to the level of safety stocks. The Monte Carlo simulation (MCS) technique was adopted to capture the estimate of solution reliability. Cordeau and Maischberger [22] presented a parallel iterated tabu search heuristic for solving four different vehicle routing problems: classical VRP, periodic VRP, multi-depot VRP, and site-dependent VRP. Vidal et al. [23] provided a survey of heuristics and meta-heuristics for multi-attribute vehicle routing problems, and pointed out that it is an important step for handling the large range of vehicle routing variants. Besides the above various methods, there are a large amount of algorithms to solve VRP including tabu search, local search, metaheuristics, and intelligent algorithms, etc. [24 –26].

By relaxing the constraints imposed on VRP, several expansions of VRP have been proposed and studied by researchers. Ho et al. [27] presented two hybrid genetic algorithms (HGAs) for multi-depot VRP. HGA1 was used to generate the initial solutions, and HGA2 was utilized to optimize the initial solutions. Osvald and Stirn [12] considered the distribution of fresh vegetables into traditional VRP with time windows and time-dependent travel-times, and developed a tabu search approach to solve the problem. Chen et al. [28] presented a nonlinear mathematical model to formulate production scheduling and vehicle routing problems with time windows for perishable food products. Karaoglan et al. [29] proposed an effective branch-and-cut algorithm for solving the location-routing problem with simultaneous pickup and delivery. The algorithm was proven to be very effective for obtaining the satisfied solution. Çetinkaya et al. [30] introduced the two-stage vehicle routing problem with arc time windows, and a mixed integer programming (MIP) formulation and a heuristic approach were proposed to solve this issue. Wang et al. [6] studied the vehicle routing problem of simultaneous deliveries and pickups with split loads and time windows (VRPSDPSLTW) as a mixed integer programming problem. A hybrid heuristic algorithm was proposed to test the algorithm effectiveness by using the modified Solomon datasets. To sum up, the aforementioned VRP variations reflect customers’ actual delivery requirements. For instance, customers may request fresh vegetable or pharmaceutical drugs for online shopping, and this requires a fleet of temperature-controlled vehicles to deliver [31]. Many retailers provide merchandise return service for customers, and customers may choose to return goods that they are not satisfied with. This is also known as reverse logistics. In these cases, delivery and pickup activities for each customer can be conducted simultaneously, and thus vehicle route optimization problem becomes more difficult [7]. Another example is door-to-door shipping. The delivery person needs to get off the delivery vehicle for dropping off the parcel at the door and comes back to the vehicle on foot. The amount of time varies by different types of buildings such as apartments, colleges, arcades, and should be taken into account into VRP [32]. These above examples imply that grouping those customers with similar characteristics or shopping habits should be taken into account into VRP for reducing algorithm complexity.

The fuzzy set theory defines several types of membership functions to transform linguistic attributes into fuzzy numbers, such as trapezoidal fuzzy numbers, triangular fuzzy numbers and interval fuzzy numbers [4 , 34–39]. In this study, trapezoidal fuzzy numbers are used as the general form of fuzzy numbers. Most of the fuzzy-set-based VRP literatures concentrated on modeling either demand or supply uncertainties, very few studies have been undertaken to apply the fuzzy set theory into depicting the vagueness of customers’ characteristics for clustering. As mentioned in the introduction, customer characteristics clustering should be considered during the optimization procedure for logistics distribution networks, because integrating both pre-route customer grouping and en-route goods delivery will improve the effectiveness and efficiency of logistics system [35]. Hu and Sheu [40] proposed a new approach based on fuzzy clustering techniques that can be utilized to cluster customers into multiple groups before executing fleet routing in logistical operations, five major attributes were adopted: safety, transit time, transportation cost, accessibility and service quality. Based on the research by Hu and Sheu [40], Sheu [35] further presented a logistics distribution methodology based on integrated fuzzy-optimization customer grouping for quickly responding to a variety of customer demands, and the multiple attributes of customer demands and static geographic attributes were considered in clustering customers’ orders and planning vehicle routes.

This study aims to construct a vehicle routing optimization model in consideration of heterogeneous customer characteristics, and then develop a fuzzy clustering algorithm and fuzzy dynamic programming approach with the hierarchical analysis for solving this problem. Compared with the previous studies, the main contributions of this paper lie in: (1) Constructing a hierarchical analysis structure, where the rating for each criterion varies according to its level of importance; (2) Designing a fuzzy clustering algorithm to comprehensively evaluate customer characteristics and group customers into different clusters; (3) Utilizing the fuzzy TOPSIS approach to calculate individual customer’s service priority and integrating with traditional VRP; (4) Developing a dynamic programming approach to optimize vehicle routes within each distribution region, and conducting sensitivity analysis to capture the impacts of different evaluation criteria weights.

3 Methodology

The methodology contains four procedures(see Fig. 1): 1) major criteria and sub-criteria are chosen to form a hierarchical analysis structure; 2) linguistics variables used for customer clustering and weight evaluation are then defined and transformed into trapezoidal fuzzy number for further vehicle routing problem optimization; 3) a fuzzy clustering framework is proposed to group customers into multiple service distribution regions with similar attributes; and 4) vehicle routing problem optimization in each region is finally developed based on individual customer’s service priority (e.g. weight evaluation) and travel distance.

3.1 Hierarchical analysis structure for customer clustering and weight evaluation

The hierarchical analysis structure is used for customer clustering and weight evaluation in the urban logistics distribution network. As shown in Fig. 2, six major criteria (P _i) and fifteen sub-criteria (P _ij) are selected to determine to categorize the customers with similar characteristics into several groups and calculate the customer weights. These criteria can be obtained from previous studies [7 , 41] and are considered important through discussions with several decision-makers and members in the local transportation and business management departments.

The decision makers including three urban planners and two senior engineers used linguistic terms to evaluate the each customer based on the above hierarchical analysis structure. The score of a major criterion evaluation can be calculated from these sub-criteria evaluations. Due to space constraints, the detailed explanation of each criterion is not expanded in this paper. Interesting readers can refer to another literature by Wang et al. [7]. Two additional criteria are added: P ₂₄ (indicates whether monopolistic and competitive products can be delivered together) and P ₄₂ (indicates whether customers request split deliveries and pickups, or simultaneously deliveries and pickups services). For the convenience of calculating customer weights, the total 15 sub-criteria can be divided into two categories: cost criteria and benefit criteria. Sub-criteria P ₂₂, P ₂₃ and P ₅₂ fall into cost criteria, and the remainingsub-criteria are grouped as benefit criteria.

3.2 Linguistic variables fuzzification and related definitions

3.2.1 Linguistic variables fuzzification

In this study, natural languages are used to evaluate the characteristics of customers in hierarchical analysis structure. The fuzzy set theory is used to convert natural languages into trapezoidal fuzzy number. Trapezoidal fuzzy number is presented as $\tilde{n} = (n_{1}, n_{2}, n_{3}, n_{4})$ [7 , 42], and the membership function $μ_{\tilde{n}} (x)$ can be then calculated based on the trapezoidal fuzzy number $\tilde{n}$ : $μ_{\tilde{n}} (x) = {\begin{matrix} 0, & x \leq n_{1}, \\ \frac{x - n_{1}}{n_{2} - n_{1}}, & n_{1} \leq x < n_{2}, \\ 1, & n_{2} \leq x < n_{3}, \\ \frac{n_{4} - x}{n_{4} - n_{3}}, & n_{3} \leq x < n_{4}, \\ 0, & x \geq n_{4} \end{matrix}$ (1) where n ₁, n ₂, n ₃, n ₄ are real numbers, x is between n ₂ and n ₃ gives the maximal grade of $μ_{\tilde{n}} (x)$ , and x at n ₁ and n ₄ gives the minimal grade of $μ_{\tilde{n}} (x)$ .

In the fuzzy set theory, the linguistic terms can be converted into fuzzy numbers [38 , 44]. Due to the uncertain nature of the customer clustering and weight evaluation problem in the logistics distribution network, criteria ratings and customer characteristics ratings are considered as the linguistic variables [45, 46]. As shown in Table 1, a scale of AL-AH is used to rate the level of importance for sub-criteria, and AP-AG is used to rate each sub-criterion of each customer. Table 1 provides the linguistic variables and trapezoidal fuzzy numbers for each criterion.

3.2.2 Related definitions

Before conducting the vehicle routing optimization, customer clustering algorithm and weight formulation are constructed, and several related definitions are needed and presented as follows.

Definition 1. Notations for the parameters and model are defined as shown in Table 2 as follows:

In addition, denotes the comprehensive evaluation matrix of customers; denotes the standardization of comprehensive evaluation matrix ; denotes the fuzzy weight matrix, denotes the weighted normalized matrix, A ^* represents positive ideal solution; A ^- represents negative ideal solution; $d_{i^{'}}^{*}$ represents the dimensional Euclidean distance of positive ideal solution; $d_{i^{'}}^{-}$ represents the dimensional Euclidean distance of negative ideal solution; d _i′,j denotes the actual distances between customer i′ and j. F _i′,j denotes the evaluation function of customer i′ and j.

Definition 2. x _i,t = μ _t (x _i) = f _t (x _i) represents the membership degree of the sample x _i associated with attribute f _t; In the clustering process, the attribute f _t can be further split into m _t,1, m _t,2, …, m _{t,s_t}, wherem _{t,s_t} represents the $s_{t}^{th}$ sub-attribute of attribute f _t, thereby, if we set M = {m _1,1, m _1,2, …, m _1,s₁, …, m _t,1, m _t,2, …, m _{t,s_t}, …, m _r,1, m _r,2, …, m _{r,s_r}},which is defined as the fuzzy concept (attribute) set, and the attribute set F can be expressed as F = M, and it can be used to further split the attribute set of X.

Definition 3. According to Axiomatic Fuzzy Set (AFS) theory [7 , 47–52], (M, τ, X) is set as an Axiomatic Fuzzy Set (AFS) structure, X is a sample set, M is a concept (attribute) set, τ is called structure; Let R be the binary relation of the sample set X, and the binary relation contains sub-preference relations and preference relations, and the corresponding concept of sub-preference relation is called simple concept, otherwise, R is called complex concept. In addition, Let m be the simple fuzzy concept [50], and m ∈ τ (x, y), ρ _m : X → R ⁺ = [0, ∞), if ρ _m satisfies the conditions: a) ρ _m (x) =0 ⇔ (x, x) ≠ R _m, x ∈ X; b) (x, y) ∈ R _m ⇒ ρ _m (x) ≥ ρ _m (y), x, y ∈ X; ρ _m is called membership degree of simple fuzzy concept m.

Definition 4. Let η (x _i) be the simple fuzzy concept set, S _t is the number of sub-attributes of f _t, η (x _i) = {m _t,s′|t = 1, 2, …, r ; s′ = 1, 2, …, s _t ; i = 1, 2, …, n}, let ζ _{C
_i} be the fuzzy description of each sample C _i, then, $\begin{matrix} ζ_{C_{i}} = arg min {Π m_{t, s^{'}} | m_{t, s^{'}} \in η (x_{i}), t = 1, 2, \dots, r; \\ s^{'} = 1, 2, \dots, s_{t}; i = 1, 2, \dots, n} \end{matrix}$ (2)

In addition, let ξ _{P
_γ} be the weighted fuzzy description, ξ _{P
_γ} = {S _{P
_γ}, w _{P
_γ}}, where S _{P
_γ} denotes the fuzzy description of cluster P _γ, and $w_{p_{γ}}$ is the weight of attributes belong to cluster P _γ. Let ∧ be the fuzzy characteristic and ∧ can be used to combine the fuzzy description of samples (e.g., ζ _{C
₁} = {m ₁}, ζ _{C
₂} = {m ₂}, then, ζ _{C
₁} ∧ ζ _{C
₂} = {m ₁, m ₂, m ₁ m ₂}).

Definition 5. Let B is fuzzy concept set, and C ⊆ X, B ⊆ M, X ⊆ R, $\bar{m} \subseteq \bar{B}$ , where $\bar{B} (x)$ represents the set of membership degrees which are less than or equal to ρ _m (x) in the fuzzy concept set B, then J _m (x) is the measure of sample x belonging to the simple fuzzy concept m, and μ _B (x) is the membership degree of fuzzy concept B at sample x, they can be expressed as follows: $\bar{B} (x) = {y | y \in ρ_{m} (x), m \subseteq B, \forall x \in C}$ (3) $J_{m} (x) = \frac{\sum_{x \in \bar{m} (x)} ρ_{m} (x)}{\sum_{x \in X} ρ_{m} (x)}$ (4) $μ_{B} (x) = inf_{m \in B} (J_{m} (x)) \in [0, 1], x \in X$ (5)

Definition 6. Let X be the sample set of customers, M be the simple fuzzy concept set of X, B ⊆ M, and n be the sample size, the membership information entropy [53] function E (B) and the membership distribution coefficient function D (B) are defined as follow: $E (B) = - \sum_{x \in X} (μ_{B} (x) ln (μ_{B} (x)))$ (6) $D (B) = - ((\sum_{x \in X} μ_{B} (x) / n) ln (\sum_{x \in X} μ_{B} (x) / n))$ (7)

When E (B) becomes smaller, the membership degree of the sample x within fuzzy concept B becomes closer to the two ends of the interval [0, 1], and the boundaries are clearer. When D (B) becomes smaller, the membership degree of the sample x within fuzzy concept B is more approaching the left end or right end of the interval [0, 1] [52]. Therefore, the evaluation index can be defined as V ^′′ = E (B)/D (B) for comprehensively evaluating the membership information entropy and distribution coefficient. When V ^′′ is smaller, it is more reasonable that the fuzzy concept B describes the sample X.

3.3 Vehicle routing optimization based on customers clustering algorithm

There are three major steps included in this section. The evaluation sub-criteria by m′ decision makers should be properly mapped into the higher hierarchical criteria based on the fuzzy integration method in the first step, and the clustering algorithm is undertaken to group the customers into different clusters in the second step, and then the clustering validity index is designed to determine the reasonable number of clusters and the customers’ clustering units, the last step is to finalize the vehicle routing optimization within each cluster by using the fuzzy TOPSIS and dynamic programming approaches. These steps are presented in detail as follows:

3.3.1 Fuzzy integration method based on sub-criteria

According to Definition 1, $W_{u, i, t_{l}}^{2}$ and $X_{u, t_{l}}^{2}$ are respectively expressed as: $W_{u, i, t_{l}}^{2} = (θ_{u, i, t_{l}}^{2}, h_{u, i, t_{l}}^{2}, g_{u, i, t_{l}}^{2},$ $k_{u, i, t_{l}}^{2})$ and $X_{u, t_{l}}^{2} = (a_{u, t_{l}}^{2}, b_{u, t_{l}}^{2}, c_{u, t_{l}}^{2}, d_{u, t_{l}}^{2})$ . The comprehensive evaluation index for customer i under major criterion t can be calculated as

$\begin{matrix} Y_{t, i}^{1} & = & \frac{1}{m^{'} \times t_{s}} \otimes \sum_{u = 1}^{m^{'}} ((W_{u, i, t_{1}}^{2} \otimes X_{u, t_{1}}^{2}) \\ \oplus & \dots \oplus (W_{u, i, t_{s}}^{2} \otimes X_{u, t_{s}}^{2})) \end{matrix}$ (8) where ⊗ denotes the vector multiplication, and ⊕ denotes the vector addition. $Y_{t, i}^{1}$ denotes the comprehensive trapezoidal fuzzy number for evaluating the customer i under major criterion t, m′ is the number of experts, and t _s is the number of sub-criteria s under major criterion t.

Let Y = (a, b, c, d) be a trapezoidal fuzzy number. The integrated representation of trapezoidal fuzzy number Y is defined as $P (Y) = \frac{1}{8} (a + 3 b + 3 c + d)$ [7 , 55]. Thereby, according to Definition 1, $Y_{t, i}^{1}$ can be expressed as $Y_{t, i}^{1} = (A_{t, i}^{1}, S_{t, i}^{1}, T_{t, i}^{1}, Q_{t, i}^{1})$ , then the integrated membership degree $μ_{t, i}^{1}$ of the customer i under major criterion t is described as $μ_{t, i}^{1} = \frac{1}{8} (A_{t, i}^{1} + 3 S_{t, i}^{1} + 3 T_{t, i}^{1} + Q_{t, i}^{1})$ (9)

3.3.2 Fuzzy clustering algorithm procedure

The integrated membership degree $μ_{t, i}^{1}$ calculated from the Equation (9) is used as the input of clustering approach. The next step is to conduct clustering analysis for grouping the customers into different clusters. Axiomatic fuzzy set (AFS) theory logic has been proven as an effective approach for tackling human perception related clustering problem [47 , 50]. Since the traditional AFS method is not appropriate to include a large number of criteria [51], an improved AFS theory logic algorithm is introduced as follows:

Step 1. According to Definition 2, the attribute f _t can be transformed into five sub-attributes: m _t,1, m _t,2, m _t,3, m _t,4, m _t,5, and their membership degree can be expressed as ρ _{m
_t,1}, ρ _{m
_t,2}, ρ _{m
_t,3}, ρ _{m
_t,4}, ρ _{m
_t,5} respectively, and let $f_{t} (x_{i}) = μ_{t, i}^{1}$ , and then, the membership degrees can be written as $\begin{matrix} ρ_{m_{t, 1}} (x_{i}) = f_{t} (x_{i}) \\ ρ_{m_{t, 2}} (x_{i}) = h_{t, 1} - f_{t} (x_{i}) \end{matrix}$

$\begin{matrix} ρ_{m_{t, 3}} (x_{i}) = | f_{t} (x_{i}) - h_{t, 2} | \\ ρ_{m_{t, 4}} (x_{i}) = h_{t, 3} - | f_{t} (x_{i}) - h_{t, 2} | \\ ρ_{m_{t, 5}} (x_{i}) = | h_{t, 2} - | f_{t} (x_{i}) - h_{t, 3} | | \\ h_{t, 1} = max {f_{t} (x_{1}), f_{t} (x_{2}), \dots, f_{t} (x_{n})} \\ h_{t, 2} = \frac{f_{t} (x_{1}) + f_{t} (x_{2}) + \dots + f_{t} (x_{n})}{n} \\ h_{t, 3} = max_{1 \leq i \leq n} {| f_{t} (x_{i}) - h_{t, 2} | + h_{t, 2}} \\ t = 1, 2, \dots, r \end{matrix}$ (10)

Step 2. Calculate the fuzzy attributes of each sample.

Step 2.1. Calculate j _{m
_t,1} (x _i), j _{m
_t,2} (x _i), j _{m
_t,3} (x _i), j _{m
_t,4} (x _i), j _{m
_t,5} (x _i) by following step 1 and the Equation (4).

Step 2.2. Let f _{m
_t,1} (x _i) = J _{m
_t,1} (x _i), f _{m
_t,2} (x _i) = J _{m
_t,2} (x _i), f _{m
_t,3} (x _i) = J _{m
_t,3} (x _i), f _{m
_t,4} (x _i) = J _{m
_t,4} (x _i), f _{m
_t,5} (x _i) = J _{m
_t,5} (x _i), then, fuzzy concept of the maximum of membership fuzzy values can be expressed as: $\begin{matrix} η_{t, s^{'}} (x_{i}) = {m_{t, s^{'}} | f_{m_{t, s^{'}}} (x_{i}) = max {f_{m_{t, 1}} (x_{i}), f_{m_{t, 2}} (x_{i}), \\ f_{m_{t, 3}} (x_{i}), f_{m_{t, 4}} (x_{i}), f_{m_{t, 5}} (x_{i})}, t = 1, 2, \dots, r, \\ s^{'} = 1, 2, \dots, 5, i = 1, 2, \dots, n} \end{matrix}$ (11)

Step 2.3. Calculate the ratio of membership information entropy and distribution coefficient function. In addition, let η (x _i) be the simple fuzzy concept set, and the evaluation index and η (x _i) can be calculated as follows: $V_{η_{t, s^{'}} (x_{i})} = E (η_{t, s^{'}} (x_{i})) / D (η_{t, s^{'}} (x_{i}))$ (12) $\begin{matrix} η (x_{i}) & = & {η_{t, s^{'}} (x_{i}) | t = 1, 2, \dots, r; \\ s^{'} & = & 1, 2, 3, 4, 5; i = 1, 2, \dots, n} \end{matrix}$ (13)

Step 2.4. Find the smallest value μ _λ (x _i) corresponding to attribute λ and the second smallest value μ _β (x _i) corresponding to attribute β. μ _λ (x _i) and μ _β (x _i) can be defined as $μ_{λ} (x_{i}) = inf_{σ \in η (x_{i})} (J_{σ} (x_{i})) \in [0, 1]$ (14) $\begin{matrix} μ_{β} (x_{i}) & = & inf_{σ \in η (x_{i})^{'}} (J_{σ} (x_{i})) \in [0, 1], \\ η (x_{i})^{'} & = & {η (x_{i}) - {λ}} \end{matrix}$ (15) Step 2.5. Select the evaluation index V _{η_t,s′(x_i)} corresponding to attribute α as well as the evaluation index V _{η_t,s′(x_i)′} corresponding to attribute β, and then compare V _{η_t,s′(x_i)} and V _{η_t,s′(x_i)′}. If V _{η_t,s′(x_i)′} ≥ V _{η_t,s′(x_i)}, the attribute λ can be eliminated. Then, the remaining sample attributes become η (x _i) ′ = {η (x _i) - {λ}}, and return to step 2.4. This recursion continues until V _{η_t,s′(x_i)′} < V _{η_t,s′(x_i)}. Consequently, the final remaining attributes η (x _i) ′ for each sample x _i is returned, and can be used to describe the sample x _i. Finally, the fuzzy attributes of each sample are generated from step 2.1–step 2.5.

Step 3. Calculate the final clusters based on the fuzzy attributes of each sample.

Step 3.1. Given that the fuzzy description of each sample C _i is ζ _{C
_i}, and we construct the fuzzy sample relation q _ij = min {μ _{ζ_{C
_i}∧ζ_{C
_j}} (C _i) , μ _{ζ_{C
_i}∧ζ_{C
_j}} (C _j)}, Liu [47] has demonstrated that the integer θ exists for $(Q_{\land}^{θ})^{2} = (Q_{\land}^{θ})$ . We can deduct the equivalence relation based on the fuzzy sample relation matrix $Q = Q_{\land}^{θ} = (q_{ij}^{θ})_{n \times n}$ .

Step 3.2. Define the diagonal element of fuzzy sample relation matrix as $q_{ii}^{θ}$ . It can be verified as $q_{ij}^{θ} = q_{ji}^{θ} = q_{ii}^{θ}$ . The diagonal elements are expressed as α _e (e = 1, 2, …, g), and these values gradually increase according to the sequence of e. Find $q_{ij}^{θ} = q_{ji}^{θ} = q_{ii}^{θ}$ in the fuzzy sample relation matrix; the corresponding samples can be grouped into one or multiple clusters, and the remaining samples where $q_{ij}^{θ} > α_{e}$ can be grouped as another cluster. This procedure will recursively continue until e = g, then, we can obtain the initial clusters $Φ_{1}^{'}, Φ_{2}^{'}, \dots, Φ_{c^{'}}^{'}$ .

Step 3.3. Calculate the weights of fuzzy description of $Φ_{r^{'}}^{'}$ : $ξ_{Φ_{r^{'}}^{'}} = {S_{Φ_{r^{'}}^{'}}, w_{Φ_{r^{'}}^{'}}}$ , $ξ_{x} = {m_{t, s^{'}} | t \in {1, 2, \dots, r}, s^{'} = 1, 2, 3, 4, 5, x \in Φ_{r^{'}}^{'}}, ξ_{x}$ denotes the fuzzy attributes set of $Φ_{r^{'}}^{'}, S_{Φ_{r^{'}}^{'}} = {m | m \in ξ_{x},$ $x \in Φ_{r'}^{'}}$ , $w_{Φ_{r^{'}}^{'}} = {w_{m} | m \in S_{Φ_{r'}^{'}}}, w_{m} = \frac{|Φ_{r^{''}}^{m}|}{\sum_{x \in Φ_{r^{''}}} |ξ_{x}|}$ , $| Φ_{r^{''}}^{m} | = | {x | x \in Φ_{r^{'}}^{'}, m \in ξ_{x}} |$ , where | · | is defined as the number of elements, and $W_{m} w_{Φ_{'}^{'}}$ is the weight of $m \in S_{Φ_{r^{'}}^{'}}$ .

Step 3.4 Calculate the weighted membership degree $μ_{(c, w, r^{'})} (x) = \sum_{m \in S_{P_{r^{'}}^{'}}} w_{m} J_{m} (x)$ , x ∈ X, where J _m (x) is the measure of initial sample, and choose the maximum membership degree as . Then, we can obtain $x \in Φ_{r^{'}}^{'}$ and the final clusters are generated as Φ ₁, Φ ₂, …, Φ _c.

Step 4. Clustering validity index based on the final clusters.

According to different thresholds α ∈ [0, 1], we can select the optimal result from the clustering validity index [56]. According to definition 1, V = {v _r|r = 1, 2, 3, …, c} is the set of sample attribute center in each cluster, $v_{γ} = {v_{γ}^{m} = \sum_{x \in P_{γ}}$ ρ _m (x)/ n _{P
_γ}|m ∈ S _{P
_γ} } , 1 ≤ γ ≤ c, where ρ _m is the membership degree of simple fuzzy concept m, and n _{P
_γ} is the number of samples in cluster γ. Therefore, the clustering validity index I _α can be expressed as follows: $I_{α} = \frac{c \times (c - 1) \sum_{γ = 1}^{c} \sum_{x \in P_{γ}} \sum_{m \in ξ_{P_{γ}}} {∥ ρ_{m} (x) - v_{γ}^{m} ∥}^{2}}{2 \times α \times (\sum_{γ = 1}^{c} \sum_{i = 1, i \neq γ}^{c} \sum_{m \in ξ_{P_{γ}}} {∥ v_{γ}^{m} - v_{k}^{m} ∥}^{2})}$ (16) where $\frac{2 \sum_{γ = 1}^{c} \sum_{i = 1, i \neq γ}^{c} \sum_{m \in ξ_{P_{γ}}} {∥ v_{γ}^{m} - v_{k}^{m} ∥}^{2}}{c (c - 1)}$ denotes the dispersion degree between clusters based on their attributes, which can be used for merging different clusters, and $\sum_{γ = 1}^{c} \sum_{x \in P_{γ}} \sum_{m \in ξ_{P_{γ}}} {∥ ρ_{m} (x) - v_{γ}^{m} ∥}^{2}$ describes the closeness of attributes from different samples within each cluster. The clustering validity index I _α becomes smaller when the closeness within each cluster becomes smaller and the dispersion degree between clusters becomes greater. The cluster becomes the most distinct from other clusters when I _α is the smallest value.

3.3.3 Vehicle routing optimization within each cluster

The Fuzzy TOPSIS approach [36 , 57] is used to calculate the weights of customers within each cluster, and the approach is combined with dynamic programming approach to calculate the evaluation function value for vehicle routing optimization. Based on the previous definitions, the procedures using Fuzzy TOPSIS and dynamic programming approaches can be described as follows:

Step 1. Calculate the aggregated fuzzy ratings for the criteria and customers.

Let $\tilde{w_{t^{'}}}$ be the aggregated fuzzy weights for sub-criteria within each cluster, and $\tilde{w_{t^{'}}}$ can be expressed as $\tilde{w_{t^{'}}} = (w_{t_{1}^{'}}, w_{t_{2}^{'}}, w_{t_{3}^{'}}, w_{t_{4}^{'}})$ . Each element of $\tilde{w_{t^{'}}}$ can be calculated as:

$\begin{matrix} w_{t_{1}^{'}} & = & min_{u} {w_{u, t_{1}^{'}}}, w_{t_{2}^{'}} = \frac{1}{m^{'}} \sum_{u = 1}^{m^{'}} w_{u, t_{2}^{'}}, \\ w_{t_{3}^{'}} & = & \frac{1}{m^{'}} \sum_{u = 1}^{m^{'}} w_{u, t_{3}^{'}}, w_{t_{4}^{'}} = max_{u} {w_{u, t_{4}^{'}}} \end{matrix}$ (17)

Let ${\tilde{x}}_{t^{'}, i^{'}}$ be the aggregated fuzzy rating of customer within each cluster, and ${\tilde{x}}_{t^{'}, i^{'}}$ can be expressed as ${\tilde{x}}_{t^{'}, i^{'}} = (a_{i^{'}, t^{'}}, b_{i^{'}, t^{'}}, c_{i^{'}, t^{'}}, d_{i^{'}, t^{'}})$ . Each element of ${\tilde{x}}_{t^{'}, i^{'}}$ can be calculated as:

$\begin{matrix} a_{i^{'}, t^{'}} & = & min_{u} {a_{u, i^{'}, t^{'}}}, b_{i^{'}, t^{'}} = \frac{1}{m^{'}} \sum_{u = 1}^{m^{'}} b_{u, i^{'}, t^{'}}, \\ c_{i^{'}, t^{'}} & = & \frac{1}{m^{'}} \sum_{u = 1}^{m^{'}} c_{u, i^{'}, t^{'}}, d_{i^{'}, t^{'}} = max_{u} {d_{u, i^{'}, t^{'}}} \end{matrix}$ (18)

Step 2. Calculate the comprehensive evaluation matrix and the fuzzy weight matrix $\tilde{W}$ .

Based on step1, the comprehensive evaluation matrix and the fuzzy weight matrix can be defined as follows: $\tilde{R} = [\begin{array}{c} {\tilde{x}}_{1, 1} & {\tilde{x}}_{1, 2} & \dots & {\tilde{x}}_{1, n^{'}} \\ {\tilde{x}}_{2, 1} & {\tilde{x}}_{2, 2} & \dots & {\tilde{x}}_{2, n^{'}} \\ \dots & \dots & \dots & \dots \\ {\tilde{x}}_{k, 1} & {\tilde{x}}_{k, 2} & \dots & {\tilde{x}}_{k, n^{'}} \end{array}], \tilde{W} = [{\tilde{w}}_{1}, {\tilde{w}}_{2}, \dots, {\tilde{w}}_{k}]$ (19)

Step 3. Normalize the comprehensive evaluation matrix .

Let be the normalized comprehensive evaluation matrix, and can be shown as , t′ = 1, 2, …, k, i′ = 1, 2, …, n′. In addition, the normalized cost criteria and normalized benefit criteria can be calculated as follows:

For the normalized cost criteria: ${\tilde{r}}_{t^{'}, i^{'}} = (\frac{a_{i^{'}, t^{'}}}{c_{t^{'}}^{*}}, \frac{b_{i^{'}, t^{'}}}{c_{t^{'}}^{*}}, \frac{c_{i^{'}, t^{'}}}{c_{t^{'}}^{*}}, \frac{d_{i^{'}, t^{'}}}{c_{t^{'}}^{*}}), c_{t^{'}}^{*} = max_{i^{'}} (c_{i^{'}, t^{'}})$ (20)

For the normalized benefit criteria: ${\tilde{r}}_{t^{'}, i^{'}} = (\frac{a_{i^{'}}^{-}}{d_{i^{'}, t^{'}}}, \frac{a_{i^{'}}^{-}}{c_{i^{'}, t^{'}}}, \frac{a_{i^{'}}^{-}}{b_{i^{'}, t^{'}}}, \frac{a_{i^{'}}^{-}}{a_{i^{'}, t^{'}}}), a_{t^{'}}^{-} = min_{i^{'}} (a_{i^{'}, t^{'}})$ (21)

Step 4. Calculate the weighted normalized matrix by multiplying normalized comprehensive evaluation matrix with the fuzzy weight matrix , and can be expressed as $\tilde{V} = {[{\tilde{v}}_{t^{'}, i^{'}}]}_{k \times n^{'}}$ , among which $\begin{matrix} {\tilde{v}}_{t^{'}, i^{'}} = {\tilde{r}}_{t^{'}, i^{'}} (.) \tilde{w_{t^{'}}} = (v_{t^{'}, i_{1}^{'}}, v_{t^{'}, i_{2}^{'}}, v_{t^{'}, i_{3}^{'}}, v_{t^{'}, i_{4}^{'}}), \\ t^{'} = 1, 2, \dots, k, i = 1, 2, \dots, n^{'} \end{matrix}$ (22)

Step 5. Calculate the positive ideal and the negative ideal solutions, and the positive ideal solution A ^* and the negative ideal solution A ^- can be calculated as follows: $\begin{matrix} A^{*} = ({\tilde{v}}_{1}^{*}, {\tilde{v}}_{2}^{*}, \dots, {\tilde{v}}_{k}^{*}), where, {\tilde{v}}_{t^{'}}^{*} = max_{i^{'}} (v_{t^{'}, i_{4}^{'}}), \\ t^{'} = 1, 2, \dots, k, i^{'} = 1, 2, \dots, n^{'} \end{matrix}$ (23) $\begin{matrix} A^{-} = ({\tilde{v}}_{1}, {\tilde{v}}_{2}, \dots, {\tilde{v}}_{k}), where, {\tilde{v}}_{t^{'}} = min_{i^{'}} (v_{t^{'}, i_{1}^{'}}), \\ t^{'} = 1, 2, \dots, k, i^{'} = 1, 2, \dots, n^{'} \end{matrix}$ (24)

Step 6. Calculate the dimensional Euclidean distance of positive ideal solution and negative ideal solution.

The dimensional Euclidean distance of positive ideal solution can be computed as $d_{i^{'}}^{*} = \sum_{t^{'} = 1}^{k} d_{v} ({\tilde{v}}_{t^{'}, i^{'}}, {\tilde{v}}_{t^{'}}^{*}), i^{'} = 1, 2, \dots, n^{'}$ (25)

The dimensional Euclidean distance of negative ideal solution can be computed as $d_{i^{'}}^{-} = \sum_{t^{'} = 1}^{k} d_{v} ({\tilde{v}}_{t^{'}, i^{'}}, {\tilde{v}}_{t^{'}}), i^{'} = 1, 2, \dots, n^{'}$ (26)

Step 7. Calculate the evaluation value of customers’ priority.

Define V _i′ as the evaluation value of individual customer’s priority as shown in Equation (27), the higher V _i′ is, the higher a customer’s delivery service priority is. $V_{i^{'}} = \frac{d_{i^{'}}^{-}}{d_{i^{'}}^{-} + d_{i^{'}}^{*}}, i^{'} = 1, 2, \dots, n^{'}$ (27)

Step 8. From the definition 1, d _i′,j denotes the actual transportation distances between customer i′ and j. The evaluation function formula is used to calculate the comprehensive evaluation value, and is defined as distribution cost. This cost integrates both transportation cost and individual customer’s service priority, and is used as the weight inputs for VRP:

$\begin{matrix} F_{i_{, j}^{'}} & = & \frac{d_{i^{'}, j}}{V_{i^{'}}} = d_{i^{'}, j} V_{i^{'}}^{'}, i^{'} = 1, 2, \dots, n^{'}, \\ j & = & 0, 1, \dots, n \end{matrix}$ (28)

Step 9. The dynamic programming approach [17] is adopted to optimize vehicle routes based on the comprehensive evaluation value F _i′,j within each cluster. The route with the minimum total comprehensive evaluation value (i.e. total distribution cost) can be chosen as the optimal route.

4 Implementation and comparisons

4.1 Data source

A case study in Anshun, China is used to evaluate the effectiveness of the proposed approach. A distribution center owned by a local logistics company is located in Anshun city, and it aims to satisfy local customers’ demands in an effective and efficient fashion. To study vehicle routing problem based on customers’ attributes, several typical customers are chosen in the entire distribution network. As shown in Fig. 3, 21 representative customers (expressed as C1, C2, … , C21) have been chosen for vehicle routing optimization. The circleunder the distribution center indicates the CentralBusiness District (CBD) and the gray lines are the major arterials of Anshun City on the map.

Five decision makers E = {E ₁, E ₂, E ₃, E ₄, E ₅} were invited to evaluate the criteria and customer characteristics in Anshun City. Three are from the department of market management, one is from the local third party logistics company, and the other one is from the bureau of commerce in Anshun City. All five have many years management experience, and very familiar with the characteristics of these customers, local transportation and economic conditions.

Due to the space limit, the language variables used for evaluating the characteristics of customers and the ratings for 15 sub-criteria are not listed here, but according to Equations (8) and (9), we can integrate the sub-criteria into the major criteria and obtain the comprehensive evaluation values. The evaluation values are presented in Table 3. In addition, the actual Distances from customers to Distribution Center (DDC) are shown in Table 4 (the distance unit is mile).

4.2 Result analysis

The aggregated evaluation matrix is used as the inputs in the fuzzy clustering algorithm procedure. Following Steps 1 through 3 in the clustering algorithm procedure, the diagonal elements α _e in steps 3.2 can be generated, and we can obtain the final clustering results based on the initial clustering results and membership degree selection in the Step 3 of the Section 3.3.2. Finally, we can calculate the clustering validity index to determine the final clustering results and use threshold α to represent the value of α _e to find the optimal number of clusters. The results are demonstrated as follows:

When threshold α = 0.4165 I _α = 3.52, I _α = 4.22, there are two clusters: $\begin{matrix} Φ_{1} & = & {C 13, C 14, C 15, C 19, C 20} \\ Φ_{2} & = & {C 1, C 2, C 3, C 4, C 5, C 6, C 7, C 8, C 9, \\ C 10, C 11, C 12, C 16, C 17, C 18, C 21} \end{matrix}$

When threshold α = 0.4683, I _α = 3.71, there are three clusters: $\begin{matrix} Φ_{1} & = & {C 13, C 14, C 15, C 19, C 20} \\ Φ_{2} & = & {C 4, C 5, C 8, C 9, C 21} \\ Φ_{3} & = & {C 1, C 2, C 3, C 6, C 7, C 10, C 11, C 12, \\ C 16, C 17, C 18} \end{matrix}$

When threshold α = 0.5216, I _α = 2.65, there are four clusters: $\begin{matrix} Φ_{1} & = & {C 13, C 14, C 15, C 19, C 20} \\ Φ_{2} & = & {C 4, C 5, C 8, C 9, C 21} \\ Φ_{3} & = & {C 1, C 2, C 3, C 6, C 7} \\ Φ_{4} & = & {C 10, C 11, C 12, C 16, C 17, C 18} \end{matrix}$

When threshold α = 0.5609, I _α = 3.24, there are five clusters: $\begin{matrix} Φ_{1} & = & {C 13, C 14, C 15, C 19, C 20} \\ Φ_{2} & = & {C 4, C 5, C 8, C 9, C 21} \\ Φ_{3} & = & {C 1, C 2, C 3, C 6, C 7} \\ Φ_{4} & = & {C 12, C 18} \\ Φ_{5} & = & {C 10, C 11, C 16, C 17} \end{matrix}$

When threshold α = 0.6163, I _α = 3.52, there are six clusters: $\begin{matrix} Φ_{1} & = & {C 13, C 14, C 15, C 19, C 20} \\ Φ_{2} & = & {C 4, C 5, C 8, C 9, C 21} \\ Φ_{3} & = & {C 1, C 2, C 3, C 6, C 7} \\ Φ_{4} & = & {C 12, C 18} \\ Φ_{5} & = & {C 10, C 11} \\ Φ_{6} & = & {C 16, C 17} \end{matrix}$

Comparing all of the above final clustering results, we can find when α = 0.5216, I _α = 2.65 is the smallest value, thus, the option with four clusters is the best choice in this situation. In addition, the clustering validity index value decreases and then increases during the clustering optimization procedure. Therefore, the lowest value can be chosen as the optimal clustering result.

Then, we use TOPSIS approach to calculate the weight evaluation value of customers in each cluster. The initial weight evaluation value and comprehensive evaluation value can be also generated using Equations (17)–(28) shown in Table 5.

The dynamic programming approach [17] is used to optimize the vehicle routes, and the comprehensive evaluation values are continuously computed during the iterative computation procedure. The optimal routes can be achieved based on the minimized total comprehensive evaluation values. The distribution center is defined as DC. Thus, the optimization results in each cluster can be shown as follows.

Cluster Φ ₁: The optimization route is DC⟶ C13⟶ C19⟶ C20⟶ C15⟶ C14⟶DC, and the total distribution cost is 55.651;

Cluster Φ ₂: The optimization route is DC⟶ C8⟶ C21⟶ C9⟶ C5⟶ C4⟶DC, and the total distribution cost is 63.844;

Cluster Φ ₃: The optimization route is DC⟶ C7⟶ C3⟶ C2⟶ C1⟶ C6⟶DC, and the total distribution cost is 54.426;

Cluster Φ ₄: The optimization route is DC⟶ C12⟶ C11⟶ C10⟶ C17⟶ C16⟶ C18⟶ DC, and the total distribution cost is 66.243.

The comprehensive evaluation value can be calculated in each cluster for vehicle routing optimization. This is different from traditional VRP studies where minimizing total travel distance or travel time is often considered but little attention is paid to adjust delivery services for customers’ diverse requirements. In reality, the customers’ diverse requirements should be taken into account and included in distribution cost. Therefore, the comprehensive distribution cost value is of practical significance.

4.3 Algorithm comparisons

To further demonstrate the superiority of the proposed approach, two recent vehicle routing optimization approaches were implemented and compared with the proposed approach. Awasthi and Chauhan [36] presented a hybrid approach integrating Affinity Diagram, AHP and fuzzy TOPSIS for sustainable city logistics planning, and the proposed hybrid approach can be applied to rank the logistics customers for generating the vehicle routings. Novoa and Storer [17] developed an approximate dynamic programming approach for vehicle routing problem. For comparison purposes, we implemented Awasthi and Chauhan’s approach [36] as well as Novoa and Storer’s approach [17] in the same context within these four clusters, the results are demonstrated as follows.

Through Awasthi and Chauhan’s approach [36], we can obtain the optimization route in each cluster. The optimization results are shown as follows:

Cluster Φ ₁: The optimization route is DC⟶ C13⟶ C19⟶ C20⟶ C14⟶ C15⟶DC, and the total distribution cost is 61.923;

Cluster Φ ₂: The optimization route is DC⟶ C8⟶ C4⟶ C5⟶ C9⟶ C21⟶DC, and the total distribution cost is 70.747;

Cluster Φ ₃: The optimization route is DC⟶ C7⟶ C6⟶ C3⟶ C2⟶ C1⟶DC, and the total distribution cost is 60.235;

Cluster Φ ₄: The optimization route is DC⟶ C18⟶ C12⟶ C11⟶ C10⟶ C17⟶ C16⟶DC, and the total distribution cost is 71.245.

Through Novoa and Storer’s approach [17], we can obtain the optimization route in each cluster. The optimization results are shown as follows:

Cluster Φ ₁: The optimization route is DC⟶ C13⟶ C14⟶ C15⟶ C20⟶ C19⟶DC, and the total distribution cost is 57.944;

Cluster Φ ₂: The optimization route is DC⟶ C4⟶ C5⟶ C9⟶ C21⟶ C8⟶DC, and the total distribution cost is 66.129;

Cluster Φ ₃: The optimization route is DC⟶ C7⟶ C6⟶ C1⟶ C2⟶ C3⟶DC, and the total distribution cost is 57.208;

Cluster Φ ₄: The optimization route is DC⟶ C18⟶ C17⟶ C16⟶ C10⟶ C11⟶ C12⟶DC, and the total distribution cost is 68.913.

As shown in Fig. 4, the proposed approach outperforms the other two approaches in terms of the total distribution cost. The following findings can be drawn:

The cost function of the proposed method includes both the physical distance between customers and each customer’s service priority, and the other two approaches only consider either transportation cost [17] or customer service priority [36]. For example, when both the travel distance and service priority are integrated in Cluster Φ ₁, customer C15 is visited prior to customer C14 using the proposed approach. However, if the distance and customer priority are considered separately, customer C14 is visited prior to customer C15. This increases the total distribution cost.

The customer characteristics ratings are also considered in the clustering process. These ratings capture the homogeneity among different customers and heterogeneity between individual customer’s attributes, so the vehicle routing optimization can be executed in these groups with similar customer characteristics, and it is effective to reduce the computational complexity.

4.4 Sensitivity analysis

To further investigate the impacts of different criteria weights on the vehicle routing optimization procedure, a sensitivity analysis was also conducted. A total of 28 experiments were performed in each cluster. The sub-criteria P ₂₂, P ₂₃ and P ₅₂ are cost criteria, while the remaining sub-criteria are benefit criteria. In each cluster, the weights of all criteria were equally set to AL, VL, B.VL&L, L, FL, M, FH, H, B.H&VH, VH, AH from experiment 1 to experiment 11; The weight for one criterion is set highest as AH and the remaining criteria are set lowest as AL from experiment 12 to experiment 26; Experiment 27 set the cost criteria (P ₂₂, P ₂₃, P ₅₂) highest as VH, with the remaining criteria weights lowest as AL; In experiment 28, the cost criteria weights are set lowest as AL, with the remaining criteria highest as AH. The sensitivity analysis results can be observed in Fig. 5(a-d).

In experiment 13, 16 and 28, the customer visit sequence for cluster Φ ₁ is changed from C ₁₅⟶C ₁₄ to C ₁₄⟶C ₁₅ due to F ₁₅ > F ₁₄. In addition, in experiment16 and 28, the delivery order for customer C ₁₉ and C ₁₃ is reversed due to F ₁₃ > F ₁₉; Similarly, in experiment 2, 13, 16 and 27, the service order for customer C ₄ and customer C ₅ is changed within cluster Φ ₂; In experiment 13, 16 and 28, the customer visit sequence for C ₈ and C ₁ in cluster Φ ₃ is reordered due to F ₁ > F ₈; For cluster Φ ₄, in experiment 12, 13, 16 and 28, the delivery order for customer C ₁₇ and C ₁₆ is reversed due to F ₁₇ > F ₁₆. Through these experiments, we can conclude that the customer visit sequence is relatively insensitive to benefit criteria weights; however they are sensitive to cost criteria weights. From Fig. 5a-d, the cost criteria weights in clusters Φ ₂ and Φ ₃ are more sensitive than clusters Φ ₁ and Φ ₄, which implies that decision-makers need to assess the cost criteria weights more seriously than the benefit criteria for vehicle routing optimization.

To evaluate the effectiveness of customer clustering results, a comparison with actual commercial area distribution has been made as shown in Fig. 6. The cluster Φ ₁ corresponds to an apparel trading area; the cluster Φ ₂ belongs to a restaurant area; the cluster Φ ₃ falls into an electronic business district; the cluster Φ ₄ is located in a bulk material district [7]. Based on the mapping relationship and characteristics between customer clusters and realistic business regions, an interesting phenomenon can be found that customers living in electronic business district and restaurant area are more sensitive to the cost criteria than those in apparel trading area and bulk material district.

5 Conclusions

Integrating pre-route customer clustering with traditional VRP is considered as an effective countermeasure to enhance system efficiency, low operational cost and increase customer satisfaction. This is especially critical for those industries with high delivery requirements, such as food cold-chain logistics and special drug distribution center. This paper presents a novel approach to optimize the vehicle routing problem based on customer fuzzy clustering with similar characteristics under a hierarchical analysis structure. The hierarchical analysis structure can describe each customer’s characteristics by major and minor criteria, followed by a fuzzy clustering algorithm based on an improved AFS approach for grouping customers. Furthermore, a fuzzy TOPSIS approach is integrated into the dynamic programming approach and is used to seek the minimum total comprehensive evaluation value for optimizing the vehicle routing within each cluster.

A numerical case study to optimize vehicle routes in Anshun City, China was undertaken to demonstrate the effectiveness of the proposed method. Two recent VRP algorithms were tested: a hybrid approach used to rank the logistics customers and generate the vehicle routes by Awasthi and Chauhan [36] and an approximate dynamic programming approach to minimize the total travel cost by Novoa and Storer [17]. Compare with these two algorithms, the proposed approach incorporates both individual-level customer service priority and travel distance, and thereby can lower the total distribution cost up to 10% within a certain cluster. A sensitivity analysis is conducted in capturing the impacts of different evaluation criteria weights. The results demonstrate the proposal approach performs very well as a cost-effective decision-support tool to provide efficient goods delivery service.

This study can be further improved to enhance its usability and depth. For example, it should be desired to investigate the performance of applying the proposed method into a large-scale logistics network. Under this circumstance, conducting pre-route customer clustering will save a considerable amount of computational resources compared with directly solving VRP. In addition, it also warrants in-depth research to incorporate time-dependent traffic condition into model formulation. This will provide more insights for logistics operators to implement the proposed approach into real scenarios.

Footnotes

Acknowledgments

The authors thank the helpful comments from the two anonymous reviewers are gratefully acknowledged. This research is supported by the National Natural Science Foundation of China (Project No. 71402011, 51408019, 71471024, 71301180, 51138003), the National Social Science Foundation ofChongqing of China (No. 2013YBJJ035), the Natural Science Foundation of Chongqing of China (No. cstc2015jcyjA30012), the Scientific and Technological Research Program of Chongqing Municipal Education Commission (No. KJ1400307), China Postdoctoral Science Foundation (No. 2014M560711). In addition, Special thanks are given to Jianxin Fan from Chongqing Jiaotong University, China, Yunteng Lao and Matt Dunlap from University of Washington, USA for his valuable suggestions and revisions.

References

10.

11.

12.

13.

14.

15.

16.

17.

18.

19.

20.

21.

22.

23.

24.

25.

26.

27.

28.

29.

30.

31.

32.

33.

34.

35.

36.

37.

38.

39.

40.

41.

42.

43.

44.

45.

46.

47.

48.

49.

50.

51.

52.

53.

54.

55.

56.

57.