Abstract
This paper presents a saving-based heuristic for the vehicle routing problem with time windows and stochastic travel times (VRPTWSTT). One of the basic ideas of the heuristic is to advance the latest service start time of each customer by a certain period of time. In this way, the reserved time can be used to cope with unexpected travel time delay when necessary. Another important idea is to transform the VRPTWSTT to a set of vehicle routing problems with time windows (VRPTW), each of which is defined by a given percentage used to calculate the reserved time for customers. Based on the above two key ideas, a three-stage heuristic that includes the “problem transformation” stage, the “solution construction” stage, and the “solution improvement” stage is developed. After the problem transformation in the first stage, the work of the next two stages is to first construct an initial solution for each transformed VRPTW by improving the idea of the classical Clarke-Wright heuristic and then further improve the solution. Finally, a number of numerical experiments are conducted to evaluate the efficiency of the described methodology under different uncertainty levels.
1. Introduction
Many real-life logistics and distribution problems can be attributed to a kind of vehicle routing problems with time windows which has been widely studied. Various kinds of algorithms have been developed to deal with the problem [1], which narrows the gap between theory and practical application. In the extant theory approaches, it is usually assumed that the travel time between a pair of points is a definite constant. However, in real traffic conditions of city areas, there are a large number of uncertainties resulting in that the travel time of a pair of points will change with the different traffic conditions and thus is not always the same. Ignoring the travel time fluctuations when developing route plans for pickup and/or delivery vehicles can result in inefficient and suboptimal solutions [2]. The extant researches that incorporate variant travel times are very limited in the assumption of travel time distribution and computational efficiency. Therefore, this study is undertaken to better handle the stochasticity of travel times.
The VRPSTT is arguably one of the most challenging and practical variants of the VRP [3], and only a few studies in extant literature have addressed the VRPSTT. Laporte et al. proposed a variant of VRPSTT chance-constrained programming model and compensation model [4], where a penalty function was introduced to prorate the delay time. Park and Song subsequently constructed three new heuristic algorithms based on an extension of VRP algorithms [5]. Kenyon and Morton developed two stochastic programming models [6]. The first model minimizes the expected completion time, whereas the second model maximizes the probability of completing the project by a prespecified deadline. More recently, Lecluyse et al. introduced the variability in traffic flow into the model [7], which was used to evaluate the routes based on the uncertainty involved. At the same time, the objective function was extended by standard deviation of the travel times. Connors and Sumalee [8] and Chen and Zhou [9] studied the stochasticity of travel times from the view of travelers’ equilibrium. Zhang et al. studied the stochastic travel time vehicle routing problem with simultaneous pickups and deliveries and developed a new scatter search approach for the problem by incorporating a new chance-constrained programming method [3]. Wang and Lin presented a flexible solution methodology for the capacitated vehicle routing problem with stochastic travel times [10]. Although these studies have advanced the research on VRPSTT, practical algorithms and solution approaches are still needed because the stochasticity of travel times poses a very challenging combinatorial optimization problem.
To solve the VRPTWSTT, we present a three-stage saving-based heuristic, in which the issue of solving VRPTWSTT is transformed into a new issue that consists of solving multiple VRPTW. Also, a number of experiments are conducted to test the efficiency of the heuristic. The remainder of this paper is arranged as follows. First, the VRPTWSTT is described and analyzed in Section 2. Second, the main solution framework and detailed solution procedures are developed in Section 3. After the experimental settings and the results are presented in Section 4, conclusions are provided in the end.
2. Problem Description
This paper focuses on the VRPTWSTT which can be described in the following.
Let be a graph, where is the set of customer nodes and depot node and is the set of edges. Each customer node () is associated with a positive demand , a service time , and a time window [] with , where and denote the earliest service start time and the latest service start time that customer has preferred. Each edge () is associated with a positive stochastic travel time whose distribution function is known and can follow any known statistical distribution, either theoretical or empirical, whose expected value exists. Define as a set of identical vehicles available to deliver goods from the depot to the customers. Each vehicle has the same capacity (). Some additional constraints associated with the problem are as follows:(1)Each vehicle can perform only one route in a scheduling period (e.g., a day).(2)Each customer node can be visited only once by exactly one vehicle.(3)All vehicles begin and end their routes at the depot node .(4)No vehicle can be loaded exceeding its maximum capacity .(5)Let be a vehicle’s arrival time at customer if the vehicle services customer and let be the service beginning time of customer . Then the two rules defined in formula (1) must be observed: the vehicle has to wait until the time if it arrives earlier than ; otherwise, the vehicle can immediately begin its service when it arrives at customer . If , a penalty that is charged to the logistics service provider will be incurred and customer is called a “poorly serviced customer.” The penalty can be considered as compensation for customer as a result of delayed service:It should be noted that the problem addressed is essentially different from the vehicle routing problem with soft time windows (VRPSTW), in which customer time window limitations may be violated even in an original schedule so that the number of vehicles required and/or the total travel distance can be significantly reduced. In the VRPTWSTT, customer time window limitations must be satisfied in original routing plans given expected values for all edges but may be violated during the plan execution process, because there are many possibilities that the estimated travel time in the original plan may be delayed in an actual travel process.(6)Once a routing plan is made, vehicles must run according to the plan no matter what disruption events occur until they return to depot in order to meet delivery commitments for customers.
The route of vehicle , which is defined as (), where and represent the depot and the other nodes represent customers, is a cycle starting and ending at the depot and serving a given subset of sequential customers. Let if the edge () is included in the route of vehicle () or 0 otherwise. The route of vehicle is called “deficient route” if there is at least one poorly serviced customer among all its serviced customers.
The VRPTWSTT studied in this paper is a single-objective problem in which the total cost, including the travel cost of employed vehicles and the penalty cost of all poorly serviced customers, should be minimized. The objective function for VRPTWSTT can be defined as formula (2), where and are used to transform the travel time and the vehicle’s lag time of poorly serviced customers to travel cost and penalty cost, respectively. Generally, , implying a penalty of violating the customer time window limitation, and hence formula (2) can be naturally simplified to formula (3), where indicating a penalty coefficient when the customer time window limitation is violated. The value of can be much larger than 1, for example, 10, 20:
Formula (3) can also be understood as the total cost when . Thus the total cost contains two components, the travel cost defined as and the penalty cost defined as . Let if customer is a poorly serviced customer or 0 otherwise. Then the reliability of a solution can be defined as , indicating what percentage of the total number of customers the non-poorly serviced customers account for, and it is an important index in analyzing a solution.
Note that in formula (3) is not a constant but a stochastic value whose statistical distribution and expected value ) are known. Thus the value depends not only on the routing solution, indicated by and , but also on the actual travel times, indicated by . Different travel times have different values even for the same routing solution. For example, any accidental travel time delay may result in large value while the unobstructed traffic situation can lead to small value. In order to exclude the impact of accidental factors, the objective of the problem can be described as minimizing the expected value, which is formulated and transformed in
3. The Three-Stage Saving-Based Heuristic
The VRPTWSTT is NP-hard because it can be reduced to the general VRP, which is a well-known NP-hard problem [11]. Furthermore, according to the problem’s definition, how to balance the expected travel cost and the expected penalty cost is an important factor in constructing a solution. In general, a solution with less expected travel cost usually has more expected penalty cost. This is mainly because a route in a solution with less expected travel cost is often fully engaged in delivery tasks and hence there is often not much flexible time left to cope with unexpected transportation congestion. Once the unexpected event occurs, the penalty cost will probably be incurred. Therefore, an algorithm that is devoted to minimizing the expected total cost should be able to integrate the stochasticity of travel time and generate a solution that can well balance the expected travel cost and the expected penalty cost and hence minimize the total cost. As a result of the high complexity level of the VRPTWSTT, solution techniques capable of producing high-quality solutions in a limited time, that is, heuristics, are of prime importance.
3.1. Main Framework
Our heuristic for the VRPTWSTT consists of three stages: the “problem transformation” stage, the “solution construction” stage, and the “solution improvement” stage. The objective of the first stage is to transform the VRPTWSTT into multiple VRPTW in which the latest service time of each customer is advanced by a period of time in order to cope with unexpected travel time delays. As a result of the problem transformation, penalty costs of solutions to the transformed VRPTW can be greatly reduced. But on the other hand, if too much time is reserved, travel costs will also increase rapidly. Therefore, another problem of how much time should be reserved for each customer under different routing schemes in order to minimize the expected total costs is derived from the first stage. To solve the problem, we develop the next two stages to first construct an initial solution for each transformed VRPTW and then to further improve the solution. Finally, for all transformed problems, we can get their improved solutions, the best of which indicates the best balance between the travel cost and the penalty cost and will be reported to decision-makers.
The proposed three-stage heuristic (TSH) is described in Algorithm 1, where procedures of the three stages are described in the following subsections, respectively.
| ||||||||||||||||||||
3.2. The Problem Transformation Stage
In the problem transformation stage, the VRPTWSTT is transformed into a deterministic problem which is different from VRPTWSTT in the following two points.(1)The stochastic travel time between any pair of nodes is replaced with its expected value according to its statistical distribution.(2)The latest service start time of each customer is supposed to be advanced by a period of time, and thus vehicles should arrive at customers no later than their advanced times.
The reason for the second point comes from the fact that decision-makers of logistics service providers are always reluctant to bear the consequences of penalties but would rather schedule vehicles to arrive at customers earlier than their latest service start time in the original plans so that there are enough reserved times left for vehicles to cope with any unexpected travel time delay during practical travel process, such as road congestion. Otherwise, once the travel time delays are encountered, the penalty cost incurred will be far greater than the cost of earlier arrivals.
Furthermore, it can be easily concluded from common sense that the earlier the arrivals at customers, the smaller the penalty cost but the greater the travel cost. This conclusion can be also validated from our experiments in Section 4. Therefore, how long should be advanced for each customer so that the total cost can be minimized is an important problem.
To decide the advanced time for each customer, we use the variable (), the loop variable in the main framework of Algorithm 1, to represent the advanced percentage of time and calculate the new latest service start time for each customer, for example, customer , by the formula , where avgLeni represents the average route length from all the other nodes to customer . Figure 1 is an illustration of advancing to for customer .
Note that represents the latest time that a vehicle should arrive at customer without incurring any penalty cost. However, in our transformation version of the problem, will be advanced to so that vehicles should arrive at customer before (not ) in the original plan and the time that is equal to ) is reserved to cope with any unexpected travel time delay during practical travel process of vehicles. In the formula , is a percentage of the average route length, , from all the other nodes to customer and its value may be 0.1, 0.2, and so forth. The rationale of the formula comes from the observation that the longer the average route length is, the more likely it is to have accidents during the travel process of vehicles and thus the greater the expected penalty cost is. Meanwhile, our heuristic is essentially a simulation-based method, in which various percentages are tried and the one that generates the solution with the least total cost will be reported to decision-makers.
It should also be noted that we do not guarantee must be larger than . There are possibilities that . Once , the vehicle that serves customer also needs to arrive at customer before in the original plan. And in the practical travel process, if the vehicle arrives at customer earlier than , the vehicle also needs to wait until the time .
3.3. The Solution Construction Stage
In this subsection, a new procedure is presented to construct an initial solution to the transformed problem. The procedure is based on the classical Clarke and Wright Savings (CWS) algorithm [10]. Algorithm 2 describes the procedure, which changes the saving selection strategy of the classical CWS to a long-term perspective.
| ||||||||||||||||||||||||||
The concept of saving in this paper expresses the saving of the travel cost and penalty cost obtained by joining two routes into one route.
There is a deficiency in the saving selection work of the CWS, which prevents the CWS from generating high-quality solutions. This is because the saving selection strategy of the CWS is a greedy one, which suggests selecting customers pairs in the descending order of saving values and hence may not accurately identify the best customers pair in some cases. To deal with this problem, a new saving value selection strategy that can guide computers to high-quality solutions is developed.
Different from the descending order of saving values in CWS, our evaluation strategy is illustrated in Figure 2, where a forestlike data structure is demonstrated. In the data structure, each point represents a saving and may have its child saving points, in which each customers pair can be connected without violating any problem constraints if the pair of customers of its father point has been connected. For example, point 1 is the father point of point 1.1 and point 1.1 is its child point. The relationship of the two points indicates that the two customers of point 1.1 are feasible to be connected if the two customers of point 1 have been connected. As illustrated in Figure 2, points in the same generation belong to the same layer, for example, points 1, 2, and 3 in layer 1.
From the illustration of Figure 2, the evaluation values of point 1 can be calculated based on the evaluation values of its children whose evaluation values can be further calculated based on their children. Therefore, the evaluation values of the first layer’s points may be calculated layer by layer from the bottom layer to the top layer, which can be called a “roll-up evaluation strategy.”
However, the deeper the layer that the calculation begins from is, the more complex the calculation procedure is. Even for a small-scale problem with only hundreds of savings, the number of points in the forest is extremely large and hence it is very time-consuming and impossible to calculate the evaluation values of the top-layer points from the bottom layer. Furthermore, we do not exactly know the depth of the forest. Therefore, calculating the evaluation values of the top-layer points from an intermediate layer, for example, layer 2 or 3, not the bottom layer, becomes a feasible choice.
The evaluation values of a top-layer point and an intermediate-layer point can be calculated by formula (5), where , , and represent the evaluation value, the saving value, and the children number of point , respectively. And in formula (5) represents the weight that the saving value accounts for the evaluation value. The evaluation value of a point in the layer from which the calculation begins is equal to its saving value:
The roll-up evaluation strategy uses a long-term perspective in that it can evaluate a pair of customers not only on its saving value but also on the saving values of its feasible subsequent connections if the pair of customers is connected, while, differently, the CWS only suggests selecting customers pairs in the descending order of their saving values, which may be regarded as a greedy strategy.
The procedure which adopts the recursion idea to calculate the evaluation values of a top-layer point is detailed in Algorithm 3.
| ||||||||||||||||||||||||||||||||||||
3.4. The Solution Improvement Stage
The solution obtained from the solution construction stage can be improved because we have found from a large number of observations that the service start times of some customers in the obtained solution are easily delayed due to occasional travel time delays while others are relatively difficult. The reason that produces the result of the imbalance can be explained in the following. If is small but an edge () that has a large length is selected into the route of a solution, the reserved time for customer is smaller than expected and customers like are prone to be delayed. On the contrary, if is large but an edge () that has a small length is selected into the route of a solution, the reserved time for customer is larger than expected and customers like are not prone to be delayed. This imbalance can result in large penalty costs that are mainly incurred by the customers prone to be delayed. From a number of experiments, we find that the smaller the number of the customers that are prone to be delayed, the lower the penalty cost.
To reduce the number of the customers that are prone to be delayed, we develop a method to balance the possibilities of delaying service start time among multiple customers so that the penalty cost and the associated total cost of a solution could be reduced simultaneously. The key idea behind the method is to iterate the process of deleting from a route a customer node, which is prone to be delayed, inserting it into a relatively idle route, and evaluating the new solution by formula (4) until the total cost cannot be reduced. The method is described in Algorithm 4.
| ||||||||||||||||||||||||||||||
In Algorithm 4, we apply a procedure to iteratively delete a customer with the minimum reserved time and insert the customer into the other routes or the other locations of the same route in order to obtain a solution with less total cost. In each iteration, this procedure can be repeated. Figure 3 gives an example of one iteration, where there are two routes and each route has three customers. Customer 6 in Figure 3 has the minimum reserved time 0 and should be deleted for reinsertion. The reserved time 0 indicates that the service of customer 6 is very easy to be delayed once travel times are larger than expected during the process of delivery. There are 6 new insertion locations for customer 6 to be inserted, among which location 4 is the best one because the new solution obtained from inserting customer 4 into location 4 has the minimum total cost. Therefore, customer 6 should be moved from their route to another and the new solution can be illustrated in Figure 3(b).
4. Computational Results
In this section, we present the computational results of the TSH, which has been coded in Java and run on a laptop with an Intel® Core™ 2 Duo CPU T7100 at 1.8 GHz and 2 GB RAM.
4.1. Test Instances and Parameter Setting
Because there are no standard benchmarks for SVRP, different sets of instances have been presented in the papers of the SVRP literature [3, 4, 6], which makes it difficult to compare the different methodologies that have been presented for SVRP. Therefore, we generalized the well-known Solomon instances and obtained 49 VRPTWSTT instances by replacing the fixed travel time of each pair of customers in VRPTW with a random travel time for VRPTWSTT. These generalized instances are solved by the method presented in this paper and compared with the best-known results of VRPTW. The details of Solomon instances and their best-known solutions can be found at the website http://neo.lcc.uma.es/vrp and the method used to randomly generate travel times is given in the following paragraphs.
Since the travel times in VRPTWSTT can follow any distribution and our heuristic does not have any limitation as long as its mean value is known, we choose the log-normal distribution to model the travel time for each edge as a result of its nonnegative values. The two parameters of the distribution, the location parameter and the scale parameter , are formulated by the following two expressions, respectively:
For each instance, we changed the travel time of an edge () (, ) in VRPTW instances to a random value by setting to the Euclidean distance between the two nodes according to their coordinates. The original VRPTW instances can be particularly defined with . But for VRPTWSTT instances, we have the inequation Var[] > 0. A larger variance of an edge indicates a higher uncertainty when traveling the edge. Our heuristic is tested under three different variance levels, the relatively low variances, the medium variances, and the relatively high variances. And these variance levels are further illustrated in Table 1, in which the relationship between the variance and the mean value for each variance level is given.
Note that there are 56 Solomon VRPTW instances, among which only 49 ones are selected. Some instances are not selected because penalty costs of solutions to these instances are almost close to 0. We test the penalty costs of all the best-known VRPTW solutions under the medium variance level and filter out 7 instances (C109, C201, C202, C203, C204, C205, and C208) by using the condition , which indicates time windows of customers in these instances are wide enough to tolerate travel time delays and hence the time window limitations of these instances are greatly weakened.
In order to decide the parameter values of our heuristic, preliminary experiments were performed on some instances with different parameter configurations to determine which configuration would be effective in solving the VRPTWSTT. The preliminary experiments showed a satisfactory performance when the parameter values were set to , , , , , and , where depth represents the number of the layer from which the calculation of the evaluation of savings begins.
4.2. Detailed Results for One Instance
In this subsection, we present the numerical results for only one of the 49 instances in order to better illustrate our heuristic. Take the generalized version of the classical RC101 instance under the medium variance level () for an example, and we run our heuristic 11 () times in which different runs have different and values. But in one run we set in order to make the heuristic generate a solution for only one percentage of the reserved time. Table 2 shows the alternative solutions (for the generalized instance with medium variances) that we obtained for different values of and . For each of the obtained solutions, the following data is provided: travel costs, penalty costs, total costs, and reliabilities.
In order to better demonstrate the character of alternative solutions to the generalized RC101 instance, a graphical representation is given in Figure 4 to illustrate the changing trends of travel costs, penalty costs, total costs, and reliabilities for all values of . Definitions and solving procedures of these four indexes are given in Sections 2 and 3, respectively. From Figure 4, the following five points can be observed and discussed.(1)The index of travel costs shows a continuously growing trend. This is mainly because, with the growth of the value, time windows of customers are continuously tightened in order to reserve more time for vehicles to cope with unexpected travel time delays, but generally there may be some customers whose time window constraints are violated after it has been tightened. Therefore, a vehicle cannot complete the tasks that have been assigned when the time windows are not tightened, and hence some good routes with more tasks have to be modified.(2)On the contrary to the trend of travel costs, the index of penalty costs shows a continuously decreasing trend. This also conforms to the common sense that the more the time reserved is, the less likely it is to violate time window limitations.(3)As the weighted sum of travel costs and penalty costs, the trend of total costs is not particularly obvious. The trend of the total costs can be divided into three stages: when value is less than 0.2, the total costs demonstrate a decreasing trend, which indicates the penalty costs play a greater role in the contribution to total costs; when value is between 0.25 and 0.45, the total costs fluctuate; when value is greater than 0.45, the total costs grow as a result of the rapid growing trend of travel costs and a relatively stable reduction in penalty costs in this range of values. Note that when is 0.45, neither the travel cost, 2007.21, nor the penalty cost, 48.12, is the minimum, but their weighted sum, the total cost, reaches the minimum value of 2488.41.(4)The two indices of both penalty costs and reliabilities can illustrate the level of a solution’s probability of suffering deficient routes. As is expected, the index of reliabilities shows a stable growing trend, which is just the opposite to the overall trend of penalty costs. This is because the index of reliabilities denotes the percentage of the total number of poorly serviced customers while the penalty cost is the sum of penalty costs of these customers. As the value grows, the reserved time and the reliabilities increase simultaneously. When the reliability reaches its maximum value (0.99), there are two solutions of different values (0.45 and 0.5). Relative to the penalty cost index, a solution’s reliability provides decision-makers with another index to evaluate the reliability of a solution. If a decision-maker cares about the cost of penalties, he can consider the index of penalty cost; and otherwise, if a decision-maker wants to minimize the number of poorly serviced customers, the index of reliabilities may come in handy.(5)As for the running time of our heuristic, the process of generating the results for one alternative solution to the generalized RC101 instance can be completed in just fourteen seconds (average 14359 milliseconds on our laptop whose configuration is relatively low). Therefore, it seems completely feasible to implement the process for most real cases.
4.3. Results for 49 Instances
The average results obtained for all instance sets (c1, c2, r1, r2, rc1, and rc2) under three variance levels are given in Table 3. The meaning of each column is as follows:(1)BKS expected total cost: the expected total cost of best-known solution (BKS) obtained from the website http://neo.lcc.uma.es/vrp/, calculated according to formula (4).(2)TSH expected total cost: the expected total cost of the solution generated by the TSH formulated by (4).(3)Gap: the gaps in expected total costs between BKS and TSH, calculated as , where BKS and TSH stand for the values of their expected total costs.
The results of each generalized instance under any of the three uncertainty levels have similar statistical characteristics to those of the generalized RC101 instance discussed in Section 4.2.
From Table 3, the following four points can be observed and discussed.(1)In a low-variability level defined by , results show that the expected total costs of TSH tend to provide a better average reliability level of 0.9929 when applied to the VRPTWSTT. As is discussed before, higher reliability levels imply fewer numbers of poorly serviced customers and thus less expected penalty costs. The average gap is 2.07% between expected total costs of BKS and TSH. The expected total costs of BKS are the costs when travel times are deterministic instead of stochastic. From the table, the average value is 0.292; that is, averagely, the latest service start time of each customer can be advanced by 29.2% of its average routing time from all the other nodes to it to attend unexpected congestion.(2)As for the case of the medium-variability level with , the average reliability level of the BKS is about 0.9923, which is lower than that of the low-variability level, indicating a higher degree of uncertainty. And the average gap in expected total costs between BKS and TSH is about 63.04%, which is much larger than 2.07% of the low-variance level. Important information about all positive gaps is also given, implying an absolute advantage for TSH over BKS under this level.(3)The average reliability is only 0.920 for BKS under the high-variability level and a big average gap of 204.16% exists in expected total costs between BKS and TSH. The average reliability index is 0.9920 for TSH. Also, the average value is 0.387, which is larger than the previous levels.(4)From the table, we can also observe that all reliabilities of TSH are higher than those of BKS. A histogram is given in Figure 5 to graphically demonstrate and compare the average reliabilities of BKS and TSH solutions under three uncertainty levels.
5. Conclusions
In this paper, we present a three-stage heuristic for the vehicle routing problem with time window and stochastic travel time (VRPTWSTT). The heuristic, which does not require complex fine-tuning processes, solves VRPTWSTT by transforming it into multiple VRPTW, solving the VRPTW, and selecting the best one from the solutions of the VRPTW. The main contributions of this paper include the following.(1)The heuristic reduces a complex VRPTWSTT, where no efficient heuristics have been developed yet, to a set of VRPTW, which can be solved by many excellent heuristics, by advancing the latest service start time of each customer ahead of its schedule. The time reserved for a customer is used for vehicles to cope with any unexpected travel time delay during their real delivery process. This idea can be applied to any variant of vehicle routing problems with stochastic travel times after it is appropriately adjusted.(2)The heuristic is designed as a three-stage procedure to enhance the probability of generating high-quality solutions. The three stages are devoted to transforming the problem, constructing its initial solution, and improving the solution, respectively, so that a good balance between the travel cost and the penalty cost of a solution can be achieved and the total cost of a solution can be minimized. A number of conducted experiments show that the results produced by the heuristic for the generalized instances are better than the best-known solutions of the corresponding VRPTW in which the travel time between any pair of nodes is regarded as its expected value.(3)The heuristic is flexible in that it can apply to a VRPTWSTT with any distribution of travel times with a known mean and the travel times of different edges can have different distributions and even possible dependence. Moreover, our heuristic can provide decision-makers with various solutions to practical problems, each of which considers a different period of time by which the latest service start time of each customer is advanced, indicating different levels of risk or reliability.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work is partially supported by the grants from the National Natural Science Foundation of China (nos. 71271037, 71201055, 71201014, and 71272093) and the Fundamental Research Funds for the Central Universities (no. DUT15QY38). Those supports are gratefully acknowledged.