An improved multi-directional local search algorithm for the multi-objective consistent vehicle routing problem

ABSTRACT This article presents a multi-objective variant of the Consistent Vehicle Routing Problem (MoConVRP). Instead of modeling consistency considerations such as driver consistency and time consistency as constraints as in the majority of the ConVRP literature, they are included as objectives. Furthermore, instead of formulating a single weighted objective that relies on specifying relative priorities among objectives, an approach to approximate the Pareto frontier is developed. Specifically, an improved version of multi-directional local search (MDLS) is developed. The updated algorithm, IMDLS, makes use of large neighborhood search to find solutions that are improved according to at least one objective to add to the set of nondominated solutions at each iteration. The performance of IMDLS is compared with MDLS and five other multi-objective algorithms on a set of ConVRP test instances from the literature. The computational study validates the competitive performance of IMDLS. Furthermore, results of the computational study suggest that pursuing the best compromise solution among all three objectives may increase travel costs by about 5% while improving driver and time consistency by approximately 60% and over 75% on average, when compared with a compromise solution having lowest overall travel distance. Supplementary materials are available for this article. Go to the publishe's online edition of IIE Transactions for datasets, additional tables, detailed proofs, etc.


Introduction
Service consistency has received attention in recent years from a number of service industries whose customers receive repeatable deliveries throughout a planning horizon, such as small package delivery (Groër et al., 2009), vendor-managed inventory systems (Coelho et al., 2012), and home health care. Two primary elements of service consistency are driver consistency and time consistency. Driver consistency refers to using the same driver for a given customer as often as possible, with time consistency refering to making deliveries to customers at roughly the same time each visit. These concepts were first discussed in Groër et al. (2009) with the introduction of a routing problem variant named the Consistent Vehicle Routing Problem (ConVRP). A more recent review of routing problems with consistency considerations is provided in Kovacs, Golden, Hartl, and Parragh (2014).
Reducing transportation costs has historically been the primary objective in managerial decision-making in logistics firms. However, due to the recognized importance of service consistency, it is desirable for some firms to obtain solutions that balance transportation costs and service consistency objectives. These objectives are often conflicting due to variation in customer demands, vehicle capacity, and driver work time limits. Optimizing routing plans with regards to travel cost may result in a customer being visited by multiple drivers at highly varying times over the planning horizon. However, pursuing optimal CONTACT Ashlea Bennett Milburn ashlea@uark.edu Color versions of one or more of the figures in the article can be found online at www.tandfonline.com/uiie.
Supplemental data for this article can be accessed on the publisher's website at www.tandfonline.com/uiie. service consistency with no consideration given to transportation costs will increase travel costs. This dilemma is illustrated in the small example given in Fig. 1, in which three customers require service across a 2day horizon. Node 0 denotes the depot, and the first element of the node label at other points denotes the customer number. Arc labels provide the Euclidean distance between pairs of locations. Note that customers 1 and 3 require service on both days, whereas customer 2 only requires a visit on day 1. Suppose a vehicle has the capacity to visit two customers, and there is no limit on route length. Then on day 1, two drivers will be required. The solution for day 1 that minimizes total distance uses driver 1 to visit customers 1 and 2, in that order, and uses driver 2 to visit customer 3. The resulting tours are illustrated in Fig. 1(a) and the associated arrival times at each customer are provided in the second element of the node labels in parentheses. Now on day 2, if minimizing travel cost is the only objective, the route depicted in Fig. 1(b) can be used, requiring one driver and a distance of 3.5. However, customer 3 is visited by two different drivers across the planning horizon. If, instead, minimizing driver consistency were the only objective, the route depicted in Fig. 1(c) could be used. In this case the distance is higher (4.0 instead of 3.5) but now both customers 1 and 3 are only visited by a single driver across the planning horizon. Note that the day 2 route illustrated in Fig. 1(c) is also better than the route in Fig. 1(b) with respect to time consistency, as the arrival times to both customers 1 and 3 do not change between days 1 and 2. If the route in Fig. 1(b) were used on day 2, customer 3 would be visited at a later time on day 2 than day 1 (time 2.5 compared with time 1.0). Note also that vehicles are not allowed to wait at customer locations to improve time consistency.
The traditional approach used in the literature to study the relationship between travel cost and consistency objectives has been to enforce service consistency by incorporating side constraints in Vehicle Routing Problem (VRP) models while minimizing travel cost using a single-objective (Groër et al., 2009;Tarantilis et al., 2012, Kovacs, Parragh, Hartl, andKovacs, Golden, Hartl, and Parragh, 2015). This approach approximates the increase in travel costs that can be expected when service consistency measures are strictly enforced. However, it may be desirable in some applications to treat service consistency measures as something to strive for rather than to enforce. That is, service consistency measures can be treated as objectives along with travel costs instead of constraints, and compromise (trade-off) solutions can be sought. Modeling service consistency using this new approach requires multi-objective VRP models. Kovacs, Golden, Hartl, and Parragh (2015) make a stride in this direction by using a single weighted objective function to simultaneously consider travel cost minimization and time consistency maximization. Milburn and Spicer (2013) also introduce a multi-objective VRP variant to consider service consistency. They develop an approach to approximate the Pareto frontier for a routing problem variant with travel cost, driver consistency, and workload balance objectives.  analyze the trade-offs between the objectives of minimizing travel cost and improving driver consistency and time consistency in the context of a generalized ConVRP. Customers in their study are associated with AM/PM time windows. However, no paper of which we are aware has studied a multi-objective variant of the traditional ConVRP without time windows where the Pareto frontier is developed explicitly. The primary objective of this article is to fill this gap.
In this article, we use a multi-objective approach to study the relationship between travel cost and service consistency objectives. Specifically, the three objectives considered include minimization of traveling cost, maximization of driver consistency, and maximization of time consistency. There are three primary contributions made by this article. First, to the best of our knowledge, this is the first article to approach the traditional ConVRP with both driver and time consistency from the perspective of developing the Pareto frontier. This enables us to study the explicit trade-offs between travel cost and consistency objectives, and use the results of our computational study to make observations that can facilitate managerial decision-making. Second, we propose an improved variant of the Multi-Directional Local Search (MDLS) framework for general multi-objective optimization problems (Tricoire, 2012). We denote this IMDLS (Improved MDLS). The effectiveness of IMDLS is compared with the original MDLS and three traditional multi-objective algorithms: the Nondominated Sorting Genetic Algorithm II (NSGAII), the Nondominated Neighbor Immune Algorithm (NNIA), and the Strength Pareto Evolutionary Algorithm 2 (SPEA2) (Deb et al., 2002;Zitzler et al., 2002;Gong et al., 2008). Additionally, IMDLS is compared with two more recent multi-objective algorithms: the Nondominated Sorting Genetic Algorithm III (NSGAIII) and the Multi-Objective Evolutionary Algorithm based on Decomposition (MOEA/D) (Qingfu and Hui, 2007;Deb and Jain, 2014). The competitive performance of IMDLS is observed on a variety of benchmark instances taken from the ConVRP literature. Finally, Large Neighborhood Search (LNS) operators are designed to explicitly improve each of the three objectives studied in this article. These are different from the LNS operators in the literature, which are designed to improve travel cost while enforcing service consistency measures as constraints.
The remainder of this article is organized as follows. Section 2 provides a review of the ConVRP literature as it pertains to the application areas of small-package delivery and home health care. The multi-objective problem we study is formally introduced in Section 3. Our proposed solution methodology is detailed in Section 4. Then, the results of our computational study are provided in Section 5 with final remarks given in Section 6.

Related literature
This section provides a brief review of existing research regarding service consistency in routing problems. The papers are classified based on two primary application areas. Groër et al. (2009) are the first to introduce the ConVRP. The problem is presented in the context of the small-package delivery industry. A formulation of ConVRP as a mixed-integer program is provided. The objective is to minimize total transportation costs and there are side constraints to enforce consistency requirements. Specifically, each customer is required to be visited by the same driver at roughly the same time on each day that he/she needs service. The time consistency aspect is modeled by imposing an upper limit on the arrival time differential for each customer. For a single customer, this is measured as the difference between the earliest and latest arrival time to the customer. The authors employ a record-to-record travel heuristic (ConRTR) to quickly obtain near-optimal solutions. Optimal driver consistency is enforced and time consistency is encouraged throughout the search process by utilizing the concept of template routes. Template routes are a set of artificial routes that only include visits to customers that require service on more than one day of the planning horizon. To construct actual routes for a specific day d, customers that do not need service on day d are removed from the template and customers who require service only on day d are added. Two additional papers employ the concept of template routes for the ConVRP. Namely, Tarantilis et al. (2012) and  develop a Template-based Tabu Search heuristic (TTS) and a Template-based Adaptive Large Neighborhood Search heuristic (TALNS) for ConVRP, respectively. Both of these algorithms outperform ConRTR with respect to solution quality, with TALNS offering the best performance among these three approaches. A relaxed variant of the ConVRP is also considered in , in which vehicles are allowed to wait at the depot before beginning their routes. Improved performance with respect to time consistency is observed in this problem variant. Kovacs, Golden, Hartl, and Parragh (2015) generalize the ConVRP in a number of ways. First, instead of requiring that exactly one driver visit each customer throughout the planning horizon, they include a constraint that limits the maximum number of different drivers that visit a customer. Second, the maximum arrival time differential used to model time consistency is penalized in the objective function instead of being included as a constraint. Additionally, customers are each associated with AM/PM time windows and vehicles are allowed to wait at the depot before beginning their routes, as in . The minimization objective considered in this ConVRP generalization is a weighted sum of travel cost and maximal arrival time differential. A large neighborhood search that does not make use of template routes is proposed. Note that the ConVRP is a special case of the new problem these authors propose. Therefore, they use LNS to solve ConVRP and compare results with ConRTR, TTS, and TALNS. The superior performance of the new approach is validated. Based on this research,  recently studied the multi-objective generalized con-VRP, in which improving driver consistency and arrival time consistency and minimizing travel cost are treated as independent objectives. Two exact solution approaches based on the -constraint framework are used to solve small instances with 10-12 customers and a planning horizon of 3 days to optimality. A heuristic multi-directional LNS algorithm is proposed to solve large instances. Trade-off analysis shows that 70% better arrival time consistency can be achieved by increasing travel cost by not more than 3.84%. Smilowitz et al. (2013) designed a modified tabu search heuristic to study the effects of three workforce management principles on routing costs in periodic VRPs. The workforce management principles include driver consistency, customer familiarity, and region familiarity. Three objectives are proposed to model these principles and are combined into a single weighted objective. Experimental results show that trade-off solutions having relatively good performance with respect to the workforce management principles can be obtained at the cost of increasing traveling cost by at most 5.3%. Zhong et al. (2007) develop a two-stage model to improve drivers' familiarity with their service territories in a vehicle dispatching problem for local package deliveries. In the first stage of their model, a set of core areas to serve as service territories for each driver are designed. Then, the optimal sequences of visits to customers in each core area are determined in the second stage. The value of the two-stage model is compared with a single-stage "no-core area" model that simply reoptimizes routes on a daily basis. Computational results show that the two-stage model is able to provide more consistent service than the no-core model. Luo et al. (2015) study a multi-period VRP in which customers require visits within specific time windows over a time horizon and customers can be served by at most a certain number of different vehicles over the planning horizon. The primary objective is to minimize the number of vehicles required, and in the case of ties, solutions are evaluated using a secondary objective of minimizing the total travel distance. A mixedinteger programming model is proposed and the limited visiting quota is enforced using hard constraints. A three-stage heuristic approach is developed to solve the problem. Computational results show that when vehicle capacity is not the primary limiting constraint, consistent customer service can be achieved with slight increases in operational costs. However, when only a small number of customers can be serviced by a vehicle due to its capacity limit, enforcing service consistency leads to more significant increases in operational cost.

Service consistency in the home health care industry
Service consistency is especially important when the customer must be present to receive the delivery or service, as in home health care. Home health care involves using skilled licensed professionals to provide prescribed medical services to homebound patients over a prescribed episode of care. Typically, the duration of an episode of care is 7 to 10 weeks (Dey et al., 2011). In this industry, service consistency is recognized as key to improving customer satisfaction and loyalty (Woodward et al., 2004). Driver (nurse) consistency can increase the familiarity between patients and caregivers and thus reduce the complexity of communication. It also improves the ability of caregivers to make accurate observations, thereby benefiting health outcomes (Woodward et al., 2004). Time consistency allows patients to plan their days more readily, without causing too many disturbances to their daily routines (Woodward et al., 2004). A number of papers in the operations research literature study home health nurse routing and scheduling problems with consistency requirements.
Bennett and Erera (2011) present a rolling horizon myopic planning approach for the single nurse routing and scheduling problem, in which patients are revealed dynamically over a time horizon. In their approach, both driver and time consistency are strictly enforced. Eveborn et al. (2006) consider nurse consistency in the development of LAPS CARE, a decision support system for determining nurse routes in Sweden. The minimization of the total number of different nurses visiting a patient is incorporated into a weighted objective function in their model. The objective function also includes a number of additional elements, such as travel time, travel cost, and patient preferences for particular staff members. They provide a set partitioning formulation of their problem and describe a repeated matching algorithm for its solution. Macdonald et al. (2009) also study a nurse scheduling and routing problem where the objective includes a weighted sum of travel cost and the number of different nurses visiting a patient. Nurse routes with good consistency and routing costs are reported. Nickel et al. (2012) study both short-term and mid-term home health nurse routing problems using an approach that aims to minimize the weighted sum of four objectives: the number of unscheduled tasks, overtime costs, patient-nurse loyalty, and travel distance. Patient-nurse loyalty is analogous to driver consistency in the ConVRP literature, defined as the number of different nurses visiting each patient during the planning horizon. Constraint programming and tabu search methods are developed.
Another paper that uses a multi-objective approach to incorporate consistency considerations in home health nurse routing and scheduling problems is Milburn and Spicer (2013). They consider a problem with nurse consistency, travel cost, and balanced nurse workload in the objective function. They approximate the Pareto frontier for these three objectives using a tabu search-based approach. Table 1 summarizes the reviewed research on routing problems with service consistency. Columns are included for time consistency and driver consistency, and a checkmark indicates whether each type of consistency is considered in each paper. The third column indicates whether the consistency elements are treated as hard or soft constraints. The review of the literature indicates there is currently a gap with respect to simultaneously considering the three objectives of travel cost, driver consistency, and time consistency within an approach aimed at approximating the Pareto frontier in the context of the traditional ConVRP. We aim to fill this gap. Most other multiobjective approaches use models that combine the objectives into a single weighted objective function. Such approaches require choosing a set of weights to represent the relative importance among objectives. Our model is based on the premise that a diverse set of decision makers may have different opinions regarding the relative importance of the three objectives. Instead of using such a set of weights to recommend a single best solution, our approach obtains an effective set of promising compromise solutions.

Problem description
The multi-objective consistent vehicle routing problem (MoConVRP) studied in this article is defined on a com- . . , n} is the set of customers and {0} indicates the depot, and A = {(i, j) | i, j ∈ N 0 , i = j}. A time horizon of D days is considered and each customer i ∈ N has a predetermined non-negative demand q id and service time s id on day d ∈ D. An auxiliary parameter w id is defined such that w id equals 1 if customer i requires service on day d (i.e., q id > 0) and equals 0 otherwise. There are K homogenous vehicles available at the depot, where they start and end their daily operations. Each vehicle is restricted by its physical capacity Q, and a limit on route duration T . The number of vehicles is unlimited; practically, we can set K to n. In this article, the terms driver and vehicle are used interchangeably. The problem is to determine a set of vehicle routes for each day of the planning horizon that are feasible with respect to capacity and route duration constraints. Each route must begin and end at the depot and each customer must be visited by exactly one vehicle on each day that he/she needs service. The objectives are to minimize travel cost and maximize driver and time consistency.
To model this problem, the following decision variables are defined: r a l i : continuous variable describing the arrival time of the latest visit to customer i over the planning horizon. Then, the multi-objective problem is formulated as follows, extended from the mixed-integer program given in Groër et al. (2009): The first objective (1) aims to minimize the total traveling distance of all vehicles on all days of the planning horizon. Objective (2) seeks to minimize the maximum number of different drivers that visit a customer. Objective (3) tries to minimize the maximal arrival time differential over all customers. Constraint sets (4) and (5) require that the depot be departed at time 0 by all vehicles on all days. Note that constraint (5) may be relaxed for the case of flexible vehicle departure time. Constraint set (6) ensures that customers are visited by exactly one vehicle on each day they need service. The vehicle capacity and route duration constraints are given in sets (7) and (8), respectively. Constraint set (9) makes sure that each customer has only one predecessor and successor. Constraint sets (10) and (11) determine the customer arrival times and also serve to eliminate subtours. Constraint sets (12) and (13) compute the maximum number of different drivers that visit a customer over the time horizon. Time consistency is computed in constraint sets (14) and (15). The remaining constraints give the variable types. Note that objectives (2), (3) and constraint sets (12)

Improved MDLS for the MoConVRP
In this section, we first introduce the basic concepts of Paretobased multi-objective optimization. Next, we describe the original MDLS algorithm and comment on its performance. We then present an improved MDLS framework for multi-objective optimization problems, followed by the description of LNS that is used as a subroutine of our proposed framework.

Multi-objective optimization
In the context of multi-objective optimization, solutions are evaluated according to an objective function vector Accordingly, a solution x in the decision space X is said If there exists no solution x ∈ X that dominates x, x is said to be Pareto optimal. The set of Pareto-optimal solutions in X is referred to as the Pareto set and its image in the objective space is called the Pareto frontier. Multi-objective algorithms aim to attain the Pareto frontier if possible. In many cases, heuristic approaches are utilized to approximate the Pareto frontier as closely as possible. In addition, the final nondominated solutions obtained by a heuristic algorithm should be distributed along the Pareto frontier as evenly as possible in order to represent a diverse approximation to the true Pareto frontier.

The original MDLS algorithm
MDLS is proposed by Tricoire (2012) for general multi-objective optimization problems. It is motivated by the concept of Pareto dominance, which implies that a neighbor solution x of x is either (i) dominating x or (ii) non-comparable with x if x is better than x on at least one objective. Therefore, it is sufficient to improve one objective at a time to find desirable neighbor solutions of x. This facilitates the utilization of singleobjective local search algorithms within the overall MDLS framework.
The input to MDLS is an initial set F of nondominated solutions. The set F can be obtained by randomly generating a population of initial solutions, assigning dominance rank using the sorting algorithm proposed in Deb et al. (2002), and deleting all dominated solutions. The solution generation scheme used to generate a population of initial solutions for the MoCon-VRP in this article, will be detailed in a later section. At each iteration, MDLS selects a solution x from F and initiates an empty set G for keeping neighbor solutions of x. Then for each of M objectives, a corresponding local search method LS m (x) is applied on x and the resulting neighbor solution x is saved in G. After that, the nondominated set F is updated by merging solutions in F and G and deleting all dominated solutions. Algorithm 1 in the online supplement summarizes the MDLS framework. Figure 2 illustrates two key steps of MDLS. In Fig. 2(a), the nondominated solution set at iteration t, F t , has four solutions. Solution A 3 is the chosen solution and its two neighbor solutions are obtained by applying local search on each of its two objectives. As the two neighbor solutions both dominate A 3 and are not dominated by any other solutions in F t , both of them enter F t+1 , as shown in Fig. 2 Our preliminary experiments reveal two limitations of MDLS. First, MDLS does not limit the size of F during its search process, which may become very large for complex problems such as MoConVRP. As F has to be updated at each iteration, maintaining a large number of nondominated solutions in F induces computational burden. In addition, MDLS focuses on only one randomly selected solution to improve at each iteration, which may lead to a slow convergence speed and poor diversity of the final nondominated set. Based on these observations, we propose an improved MDLS framework and validate its effectiveness on the MoConVRP studied in this article.

Improved multi-directional local search (IMDLS)
Algorithm 1 shows the improved MDLS framework. In addition to an initial nondominated solution set F, the input to IMDLS includes a size limit on F, denoted F max , which specifies the maximum number of solutions that can be maintained in F throughout the search process. At each iteration, IMDLS begins by exploring neighbor solutions of every solution x ∈ F with respect to each of the M objectives. As with MDLS, a variety of local search methods can be used in IMDLS to find neighbor solutions with respect to each objective. In this article, a LNS heuristic is employed. Every new neighbor solution enters a set G that is later used to update F. Performing local search on all solutions in F instead of only one solution helps IMDLS converge more quickly to the Pareto frontier. The additional computation time introduced by this expanded search is remedied, to some extent, by the time saved in updating the set F at the end of each iteration, as a reduced F is maintained. Next, if the size of the updated F exceeds F max , the crowding distance is computed for solutions in F (Deb et al., 2002). The purpose is to guide the selection of specific solutions in F to delete.

Algorithm 1 IMDLS
Input: A set of nondominated solutions F and its size limit The crowding distance of a solution on the Pareto frontier measures the density of solutions surrounding it. It is defined as the average distance of two neighboring solutions on either side of the solution along each of the objectives. For example, suppose a problem having only two objectives is being considered. To compute the crowing distance of a solution on the Pareto frontier, one must find the distances between the solutions on the Pareto frontier just better than and just worse than it for both of the objectives. Lower values of crowding distance refer to more-crowded solutions. In IMDLS, once the crowding distance has been computed for every solution in F, a function truncate(F ) iteratively removes solutions from F in order of nondecreasing crowding distance until the number of solutions that remain is no greater than F max . With this modification, IMDLS is able to focus attention on less-crowded areas of the nondominated set and hence improves the diversity of F. Figure 3 highlights the two primary differences between IMDLS and MDLS. Comparing Fig. 3(a) with Fig. 2(a), it can be seen that the neighborhoods of all solutions in F t are explored in IMDLS, whereas only a single neighborhood is explored in MDLS. Then, Fig. 3(b) illustrates the process in IMDLS for updating F when it becomes too large. Note that after the local search step, a total of eight nondominated solutions are contained in F t+1 when the size limit F max is set to 7. Therefore, F t+1 must be truncated by one. The crowding distances of all solutions in F t+1 are computed and A 2 is identified as the solution with minimum crowding distance (with a value of 0.9). Therefore, A 2 is removed from set F t+1 .

Initial solution generation
The proposed IMDLS framework requires an initial set of nondominated solutions as input. Although widely used in the literature to generate solutions with good service consistency, template routes, as described in Section 2.2, are not utilized in this article to generate initial solutions. This is because the resulting daily routes derived from template routes always have optimal driver consistency. This limits the possibility to explore trade-offs between the objectives of cost minimization and consistency maximization. Therefore, instead of employing the concept of template routes, our IMDLS framework operates on the entire routing plan. For each day of the planning horizon, the routing plan contains a route for each vehicle in service that day. The set of routes corresponding to a particular day contains visits to all customers requiring service on that day.
An overview of the initial solution generation scheme we use is provided here, with details in Algorithm 2 in the online supplement. The algorithm begins generating an initial solution (i.e., routing plan) s by constructing an empty route for each day of the planning horizon. Any time an empty route is added to the routing plan s, a driver is assigned to it. Then, customers are considered for insertion into s in random order. The first step in inserting a customer c to a routing plan is to identify a "good" driver for the customer. In general, a driver r is a good choice for customer c if they are already operating routes in close proximity to customer c on the days customer c requires service. Selecting such a driver r for customer c should encourage good driver consistency and low travel costs. To this end, a parameter α is used to identify the set of routes that are adjacent to c on each day customer c requires service. A route r is said to be adjacent to customer c if the angle made between the location of customer c and the centroid of all locations in route r, taking the depot as the origin, is no greater than α. Also, all empty routes are considered to be adjacent to all customers. Then, the driver r * appearing most frequently in the set of routes adjacent to c is determined. Next, for each day d on which customer c requires service, the visit to customer c is placed into its cheapest insertion location in route r * . If no feasible insertion locations exist in route r * on day d for the visit to customer c, the visit will be inserted into its cheapest feasible insertion location among all routes. Doing so may require creating a new empty route. Note that although this method of defining adjacency using the angle between two locations with respect to the depot is appropriate for problem instances having Euclidean distances, such as those studied in this article, it may need further investigation for problem instances having distances along real street networks. For example, two locations separated by a natural barrier such as a lake or river may be considered as adjacent using the method employed here, but the road distance between them could be quite large. In this case, the real distance between two customers could be used to decide whether they are adjacent. To decide the closeness of a customer and a route, the average distance of the customer to all the customers on the route could be used.
The objective of this generation scheme is to create an initial set of nondominated solutions that are relatively good. This scheme encourages the creation of solutions with low travel cost and good driver consistency; thus, the initial set of solutions will represent closer approximation to the Pareto frontier. Moreover, randomness is introduced in the scheme to ensure the diversity of generated solutions.

LNS
In this section, we first explain the LNS framework used in IMDLS to find neighbor solutions. Then we describe the removal and reinsertion operators designed for LNS to improve each of the three objectives considered in this article.

... Algorithm framework
Both MDLS and IMDLS are algorithm frameworks designed for general multi-objective optimization problems; they rely on problem-specific local search operators to obtain neighbor solutions. Any local search algorithm can be used for this purpose and we employ LNS in this research for its simplicity and effectiveness (Shaw, 1998). LNS searches for a neighbor of a current solution using two steps: removal and reinsertion. Essentially, removal involves choosing a number of customer visits to remove from the solution based on some specific criteria and reinsertion aims to reinsert removed customer visits back into the solution based on some other criteria. The neighbor solution is the routing plan obtained after one iteration of removal and reinsertion is executed on the current solution.
Algorithm 2 describes the workflow of LNS within the overall algorithm framework of IMDLS. Every time LNS is called with an input solution s, a pair of removal and reinsertion operators, denoted by τ , is selected and applied to s to obtain its neighbor solution s . If s dominates or is noncomparable with s, then s replaces s. The LNS subroutine will return the new neighbor solution to the IMDLS procedure. In this article, operator pairs are randomly selected in each iteration of LNS based on probabilities that depend upon their previous performance. That is, for each pair, we keep a record of the ratio of the number of times the pair has succeeded in improving a solution to the total number of times it has been invoked. This ratio is initialized to 1 for all pairs of operators at algorithm onset. Note that the removal and reinsertion operator performance is updated at every iteration in this research, while in Ropke and Pisinger (2006) and  it is updated only at certain predefined iterations.

... Removal and reinsertion operators
This section describes the removal and reinsertion operators used in LNS. First, there are four types of removal operators and two reinsertion operators. One reinsertion operator aims to achieve good travel cost, and the other to achieve good driver consistency. Any of these removal operators can be paired with the appropriate reinsertion operator if the objective of interest during the current iteration of local search is travel cost or driver consistency. However, these removal and reinsertion pairs are not used to improve the time consistency objective. This is because the removal operators tend to remove many customer visits from a solution at once, as will be described. Customer arrival times are highly interrelated, so it is difficult to reinsert several visits at once if the objective is to improve time consistency. Therefore, a designated pair of removal and reinsertion operators for improving the time consistency objective is described last.

Random removal
We employ two random removal operators to remove randomly chosen customer visits from a solution. This is intended to diversify the search process. Both begin by generating a removal probability p c ∼ U [0, 1] for each customer c. The first random removal operator, denoted Rand 1 , then randomly generates a removal threshold ρ ∼ U [0.2, 0.4]. For every customer c for which p c < ρ, the visits to the customer on every day it needs service are removed from the routing plan s. The second random removal operator, Rand 2 , randomly generates a removal threshold ρ d ∼ U [0.2, 0.4] for each day d of the planning horizon. Then, for every customer and day pair (c, d), the visit to customer c on day d is removed from routing plan s if p c < ρ d . Thus, Rand 1 removes all visits for select customers from s, whereas Rand 2 removes only a subset of visits for select customers.

Adjacent removal
There are two removal operators that aim to deconstruct routes in a particular geographic area. These adjacent removal operators define any two customer b and c as adjacent if the angle made between them, taking the depot as the origin, is no greater than α. The first, Ad j 1 , begins by randomly choosing a seed customer c. Then, the set of customers adjacent to c are identified. Finally, visits to c and all of its adjacent customers are removed from the routing plan s on every day the customers require service. The second removal operator, Ad j 2 , identifies d seed customers, one for each day of the planning horizon. Denote the seed customer on a given day d as c d . The customers adjacent to c d are identified, and then the visits on day d to c d and all of its adjacent customers on day d are removed from the routing plan s. This process is repeated for each day d. For both operators Ad j 1 and Ad j 2 , α is set to 360/K active , where K active is the maximum number of drivers working in a given day in routing plan s.

Worst-driver-consistency removal
The worst-driver-consistency removal operators attempt to remove visits to customers who are seen by the largest number of different drivers over the planning horizon. The first operator, WorstDC 1 , sets a removal threshold K active /2, where K active is as previously defined. All visits to every customer seen by more than K active /2 drivers throughout the planning horizon are removed. The second operator, WorstDC 2 , compares the number of days a particular customer c requires visits (denoted m c ) to the number of different drivers that visit c in routing plan s (denoted n c ). If these two numbers are equal (that is, m c = n c , implying that a different driver visits customer c each time it requires service), then all visits to customer c are removed from routing plan s. Otherwise, if n c < m c , then the the driver d c that visits customer c most frequently is identified. Visits to customer c are removed from s for every day c is visited by a driver different from d c . This process is repeated for every customer c.

Worst-time-consistency removal
The worst-time-consistency removal operators try to remove visits for customers whose arrival times differentials across the planning horizon are high. Both operators set a removal threshold f TC /2, where f TC is the value of the time consistency objective value of routing plan s. Denote the arrival time differential for customer c as a c . Based on the first operator, WorstTC 1 , all visits to every customer c for which a c > f TC /2 are removed from s. Based on the second operator, WorstTC 2 , if a c > f TC /2 for a particular customer c, then the arrival times for all visits to c are placed into a list L c and sorted in either increasing or decreasing order. Next, list L c is partitioned into two lists L 1 c and L 2 c by splitting the list L c at the largest difference between two consecutive arrival times. Finally, for the list with smaller cardinality, the visits to customer c associated with each of the arrival times in the list are removed from routing plan s. This process is repeated for every customer c. This operator is inspired by a similar one in Kovacs, Golden, Hartl, and Parragh (2015).

Reinsertion
Note that after removal has been carried out based on one of the above operators, a partial solution (i.e., partial routing plan) exists. That is, the partial routing plan is missing some of the required customer visits. In what follows, we let σ denote the partial routing plan that remains after removal is carried out on routing plan s.

Reinsertion to reduce travel cost (Reinsert 1 )
This operator focuses on travel cost when reinserting customer visits. For every day d of the planning horizon, a list L d of customer visits on day d that need to be reinserted into σ is maintained. Then, every customer visit i ∈ L d that requires reinsertion is placed into the cheapest feasible insertion location among all routes in σ on day d, updating the partial routing plan after each insertion. If no feasible insertion locations exist for a particular customer visit, a new empty route is added to σ . Each list L d is shuffled to introduce diversity into the order in which customer visits are considered for reinsertion.

Reinsertion to improve driver consistency (Reinsert 2 )
This operator focuses on driver consistency when reinserting customer visits. A list L of customers who need one or more visits reinserted into σ is maintained. For each customer c ∈ L, two possible scenarios exist: customer c requires reinsertion for either all of their visits or a subset of their visits. These scenarios are discussed separately below.

Customer c requires reinsertion for a subset of their visits:
For each driver, a count n d of the number of times they visit customer c in σ is maintained. For a particular day d on which c requires reinsertion, the driver with maximum n d is identified, where the driver has a route on day d and a feasible insertion location for the visit to c exists in that route. If multiple drivers meet these conditions, one is selected arbitrarily. A new empty route is added on day d if no drivers meet these conditions. The visit to customer c is inserted into the selected route on day d in the cheapest feasible insertion location. The driver count (n d ) for the selected driver is updated and the process is repeated for each day that customer c requires reinsertion. The order in which each day is considered is randomized to promote diversity.

Customer c requires reinsertion for all of their visits:
For every day customer c requires service, the cheapest insertion locations within each route for which feasible insertion locations exist for customer c are noted. Each of these locations are recorded using (d, r, j, δ), where d is the day, r is the driver/route, j is the index of the customer that c will immediately precede if inserted, and δ is the insertion cost. Next, for every driver r, a count n r of the number of such locations in which they appear is tabulated. Furthermore, for every driver r, the sum of insertion costs over all such locations in which they appear as the driver is computed and denoted . Note that n r can be at most m, the number of visits customer c requires. Note also that represents the total travel cost associated with using driver r to perform all of the identified visits to customer c. Next, the driver r * with maximal visit count n r * is identified. Visits to customer c are inserted into the previously recorded locations in driver r * 's routes (preceding customer j). If there are multiple such drivers, ties are broken using (the minimum total insertion cost will be selected). Finally, if the maximum visit count n r * is equal to m, then the process of reinserting customer c is complete. Otherwise, for every visit for which customer c still requires reinsertion, the location with minimum δ is selected. If at any point in this process there are no feasible insertion locations on a required day, an empty route is created and a visit for customer c is inserted into it. Note that when a new route is created in any of the reinsertion operators, the driver numbered r + 1 will be assigned to it if there already exist r drivers on the same day. In removal operators, it is possible that all customer visits on a route are removed; in this case, the route simply exists as an empty route and is not deleted from the routing plan. Note also that empty routes are defined to be a neighbor of any customer to be reinserted.

Removal and reinsertion to improve time consistency (RR)
This operator pair identifies the customer c * with the biggest arrival time differential. Then, the visit to customer c * that will be removed is identified as follows. First, the average of the earliest and latest arrival time to customer c is denotedā and the median arrival time to customer c is denoted med(a). Then, if med(a) >ā, the visit with the earliest arrival time is removed; otherwise, the visit with the latest arrival time is removed. Denote the day of this visit d * . The process for reinserting the visit on day d * for customer c * begins by examining each route on day d * . These routes are sorted into a list L in nondecreasing order of the angle between the location of customer c * and the centroid of the route, taking the depot as the origin. Beginning with the first route in L, the operator examines whether there is a feasible insertion location in the route that will improve the time consistency objective. If there are multiple such locations within the route, the one that improves time consistency the most is selected. If there are no feasible locations in the current route, the search proceeds through the list L. If the end of list L is reached, the visit to customer c * on day d * is reinserted back into its original location.
Although the two worst-time-consistency removal operators, WorstTC 1 and WorstTC 2 , are not used explicitly in this removal and reinsertion operator pair, they do aid in improving this objective. This is because they remove the customer with the largest arrival time differential. The removal and reinsertion operators used for each objective are listed in Table 2.

Computational experiments and analysis
In this section, we first introduce the benchmark instances used in this article to validate the performance of our proposed IMDLS for solving MoConVRP. Then, metrics used to compare the performance of the various algorithms are described. Finally, the results of the computational study are presented and a trade-off analysis between the objectives of cost minimization and consistency maximization is provided.

Benchmark instances and experiment setup
A total of 36 instances are taken from the literature to verify the performance of our proposed algorithm. They appeared first in Groër et al. (2009) and  based on the Christofides benchmark instances for VRP (Christofides and Eilon, 1969). All instances are generated in such a way that each customer requires service with a certain probability, namely, 0.5, 0.7, and 0.9, on each day of the planning horizon. These instances are classified into three groups, namely, Group 0.5 , Group 0.7 , and Group 0.9 , based on their service probability.
The performance of our proposed IMDLS is compared against the original MDLS; three classical multi-objective algorithms including NSGAII, NNIA, and SPEA2; and two more recent multi-objective algorithms, namely, NSGAIII and MOEA/D. For MOEA/D, there exist three decomposition approaches that convert a multi-objective optimization problem into a set of scalar optimization subproblems, which results in three algorithms, namely, the MOEA/D with the Weighted Sum approach (MOEA/D-WS), the MOEA/D with the Tchebycheff approach (MOEA/D-TCH), and the MOEA/D with the Penaltybased Boundary Intersection arpproach (MOEA/D-PBI). Neither MDLS nor IMDLS require the recombination of two solutions during execution. However, crossover is an essential element of the comparison algorithms NSGAII, NNIA, SPEA2, NSGAIII, and MOEA/D. The crossover operator used in this article works in a removal and reinsertion fashion as in LNS: given two parent solutions A and B, a child solution is initialized as a copy of parent A. Then all of the customers in the child solution whose driver consistency objectives are worse than that of parent solution B are removed from the child. After removal, Reinsert 2 is employed to reinsert the removed customers back to the partial child solution.
All of the comparison algorithms are implemented following their descriptions in the literature. The same crossover operator and LNS described in the literature are used whenever applicable. Both crossover and LNS operators are applied to parent solutions with probability 1 in the comparison algorithms. The population size in NSGAII is set to 100. The population size and archive size in SPEA2 are both set to 100. The sizes of the dominant population, active population, and clone population in NNIA are set to 100, 20, and 100, respectively. The population size of NSGAIII is determined by the size of the userdefined reference point set. In this article, we use the reference point set generated by Deb and Jain (2014) for three-objective optimization problems. It contains 91 reference points and the NSGAIII population size is set to 91. Similarly, the MOEA/D population size is determined by the number of weight vectors and we use the weight vectors provided by Zhang and Li (2007) for three-objective optimization problems. The MOEA/D population size is therefore set to 351. The number of initial nondominated solutions is set to 100 for both MDLS and IMDLS. The F max in IMDLS is set to 100 in all experiments. A time limit of 30 minutes is set on all the algorithms and 10 replications are run for each instance. All algorithms are implemented in C++ and experiments are conducted on a computer with 2.70 GHz CPU.

Metrics for comparison algorithms
Three metrics are employed to compare the performance of different multi-objective algorithms. First, hypervolume (I H ) is a unary operator that is able to indicate the convergence and diversity of a nondominated approximation set for the Pareto frontier (Coello et al., 2007). It measures the size of the objective space covered by a set of nondominated solutions F. A reference point z is necessary to compute the hypervolume of F. For maximization problems, it is common to set z as the origin (0, 0, 0). For minimization problems, z is usually set to a point with the worst values for each objective. Either way, larger hypervolumes indicate better performance. For our problem with three objectives, each solution x ∈ F in the objective space covers a cuboid defined by its coordinates ( f 1 (x), f 2 (x), f 3 (x)) and the reference point z. The hypervolume is computed as the size of the union of all such cuboids covered by solutions in F. The WFG algorithm descried in While et al. (2012) is used to compute the hypervolume metric. As the three objectives considered in this article do not have the same scale, it is necessary to normalize the objective values before computing the hypervolume metric. To this end, for each instance, we record the best and worst values for each of the three objectives obtained by all algorithms over all 10 replications. Then all objective values are normalized into the range [0, 1] and the reference point is set to (1.1, 1.1, 1.1) for all instances.
Second, coverage (I C ) is a binary operator that compares the convergence of different nondominated solution sets (Zitzler et al., 2000). It measures the extent to which one solution set B is covered by another solution set A by comparing the number of solutions in B that are dominated by solutions in A to the cardinality of B: The  Third, the unary multiplicative epsilon indicator (I ) computes the minimum factor by which each point in the reference set R can be multiplied such that the resulting set is weakly dominated by set A (Zitzler et al., 2003): This indicator is based on the -dominance relation, , which is defined as for a minimization problem with M objectives and assuming that all points are positive in all objectives. A smaller I value indicates a better performance. The reference set R is constructed for each instance by taking the union of all replications of all algorithms and removing dominated solutions. Note that a I value smaller than 1 indicates that A strictly dominates the reference set R.

Parameter tuning
In the proposed initial solution generation scheme, the parameter α is defined to decide whether a route r is considered as the neighbor of a customer c that is to be inserted.  Table 3 that the best performance (the biggest hypervolume value, 0.87) is observed with α = 60; therefore, this value is used in the final experiments.

Algorithm performance comparison
This section presents a comparison of the multi-objective algorithms based on the metrics described in the previous section. Then, the frequencies with which various LNS operator pairs in IMDLS are invoked are presented. Table 4 provides the hypervolumes for the instances in Group 0.5 . The first column indicates the instance number. Each of the next nine columns corresponds to one of the comparison algorithms. Each entry in the table gives the average hypervolume across 10 replications for the instance number and algorithm indicated by the row and column labels. The last two columns give the hypervolumes for the reference set R and the I % H value defined by I H (IMDLS)/I H (R), respectively. The last row provides the averages across all 12 instances for each of the comparison algorithms. The best result for each instance is highlighted in bold.   occurrences are indicated in bold in Table 7. Detailed pairwise comparison results are given in Tables 1, 2, and 3 in the online supplement. Table 8 shows the unary multiplicative epsilon indicator comparisons of IMDLS with other algorithms for instances in Group 0.5 . The first columns gives the instance number and the next nine columns correspond to one of the comparison algorithms. Each entry in the table gives the average I across 10 replications for the instance number and algorithm indicated by the row and column labels. The last row provides the averages across all 12 instances for each of the comparison algorithms. The best value for each instance is shown in bold.
It can be seen from Table 8 that IMDLS outperforms the comparison algorithms for all instances in Group 0.5 except instance 11. These observations are also supported by the comparisons Bold values indicate the best unary multiplicative indicator values among all algorithms. given in Tables 9 and 10, corresponding to instances in Group 0.7 and Group 0.9 , respectively. Note that there exist some large unary epsilon indicator values, sometimes over 300, in Tables 8, 9, and 10. They are due to very small normalized travel cost objective values. Such a scenario occurs any time the final nondominated solution set of an algorithm includes a solution with very good travel distance while another comparison algorithm fails to do so.
It can be seen from the above comparisons that IMDLS outperforms other comparison algorithms with respect to the coverage and multiplicative epsilon indicators. On the other hand, IMDLS is outperformed by other comparison algorithms with respect to hypervolume. These conflicting observations with hypervolume and unary epsilon indicators are due to the fact that these two indicators work on different principles and therefore opposite preference orderings for two approximation sets A ans B may result (Knowles et al., 2006). In the next section, we conduct our trade-off analysis based on the reference set instead of the approximation set produced by IMDLS alone.
Next the frequencies with which each removal/reinsertion operator pairs are invoked in IMDLS are reported. Table 11 provides these results for the instances in Group 0.7 . The first row gives the name of the removal operator and the second row gives the name of the reinsertion operator that is paired with it. Each of the eight removal operators is paired with both of the reinsertion operator designed for Travel Distance (TD) and Driver Consistency (DC), which results in a total of 16 operator pairs. The operator pair designed for TC is not shown in this table because it is invoked with frequency 1.0 in IMDLS by definition. Each element in the table is an average frequency over 10 replications for each instance in Group 0.7 . In each row, the entries in the eight columns corresponding to either TD or DC sum to 1.0. The average application frequency for each operator pair over all instances is shown in the last row of the table.
It can be seen from the last row of Table 11 that with regard to the TD objective, removal operator WorstDC is chosen with frequency 22.72% (18.92% + 3.80%) on average, more than any other operator. The next most frequently invoked removal operator is WorstTC, with frequency 14.26% (5.57% + 8.69%) on average. Moreover, random removal seems to be more effective in reducing travel cost than adjacent removal (Ad j 1 and Ad j 2 ).
With regard to the driver consistency objective, it is clear that WorstDC is the most effective in improving the objective with a frequency of 34.59% on average. Next, WorstTC removal is the next most frequently invoked operator to improve driver consistency. Additionally, random removal is more effective than adjacent removal in improving driver consistency. Table 11 also shows that random removal and adjacent removal are more effective in improving travel cost than driver consistency.

Trade-off analysis
In this section, we aim to analyze the nondominated reference set R obtained by all algorithms in all runs to facilitate managerial decision-making, for there usually exist many solutions in R and it is desirable to identify the best compromise solution to implement in practice. To this end, we employ the level diagram technique proposed by Blasco et al. (2008). First, the objective values of nondominated solutions in R are normalized on a [0,1] scale. Next, the Euclidean norm is computed for each solution s ∈ R, wheref i represents the normalized value of objective i: The original objective values of solutions in R are then plotted against their Euclidean distance in separate graphs corresponding to each of the three objectives. Figure 4 depicts the level diagram of nondominated reference set R for instance 1 in Group 0.9 . The x-axis in each graph provides the original objective values for TD, DC, and TC, respectively. The y-axis provides the Euclidean distances for solutions in R. Taken together, the three points at the same coordinate in each of the three graphs correspond to a particular solution in R. For example, the point denoted by a triangle in black color corresponds to the particular solution in R with a travel distance of 2412.96, driver consistency of two, and time consistency of 107.87. This point corresponds to the solution in R having lowest travel distance. In the graph, a smaller value of Euclidean distance indicates a better compromise solution. Therefore, the best compromise solution in R is denoted by a red square in Fig. 4. It has a travel distance of 2507.75, driver consistency of one, and time consistency 14.05.
Let s BC denote the best compromise solution in the reference set R. It is identified as described above. Furthermore, let s TD , s DC , and s TC denote the solutions with best travel distance, driver consistency, and time consistency value, respectively. This results in 12 values recorded for each test instance. Table 12 provides these 12 values for each test instance. Consider, for example, the best compromise solution s BC for instance 1 of Group 0.5 in Table 12. The travel distance of the best travel distance solution for this instance is 1528.72, with the driver and time consistency values being 3 and 91.00, respectively. For comparison purposes, the travel distance of the best compromise solution (s BC ) is 1642.34, approximately 7.4% higher than the travel distance in the best travel distance solution. However, the increase in travel distance results in better consistency values, with driver and time consistencies of 1 (a 66.67% decrease) and 38.13 (a 58.10% decrease), respectively. From the table, it can be seen that an improvement in time consistency typically comes at a higher cost, in terms of travel distance, than does improvement in driver consistency.    Table 13 provides a summary of the trade-offs between the traditional efficiency objective of travel cost and the consistency considerations. Specifically, it compares the three  when compared with the compromise solution with the best travel distance. That is, an approximate 60% improvement in driver consistency and a 75% improvement in time consistency comes at a cost of a 5% increase in travel distance. Now suppose that instead of seeking the best compromise among all three objectives, the manager is only interested in improving driver consistency. Then, they can expect the travel distance to increase by 5.49%, the driver consistency objective value to decrease by 64.81%, and the time consistency objective value to decrease by 76.30%, when compared with the compromise solution with the best travel distance, on average. This table also supports the conclusion that improvements in time consistency come at a greater increase in travel costs than do improvements in driver consistency.

Removal: R e i n s e r t i o n : T D ( % ) D C ( % ) T D ( % ) D C T D D C T D D C T D D C T D D C T D D C T D D C
In previous sections, it is assumed that all vehicles leave the depot at time 0. This rigidness limits the potential to further improve time consistency. In this article, we also consider the case where vehicles are allowed to wait at the depot to improve time consistency. Note that vehicles are not permitted to wait at customer locations after finishing service. In order to determine the vehicle's departure time, we employ the post-optimization algorithm proposed by Kovacs, Golden, Hartl, and Parragh (2015). This algorithm is called every time before the time consistency objective value is computed. Table 14 shows the four solutions identified from the reference set R generated by all of the mentioned multi-objective algorithms considering flexible vehicle departure time. Table 15 shows the trade-off analysis in the case of a flexible vehicle departure time. Compared with Table 13, in which a 76.47% improvement of time consistency is obtained at the cost of a 5.02% increase in travel cost, the first row in Table 15 indicates that a 73.33% improvement of time consistency could be achieved at the cost of a 4.40% increase in travel cost. This comparison shows that better time consistency could be achieved by allowing a flexible vehicle departure time.

Conclusions
This article presents a MoConVRP. Instead of modeling consistency considerations such as driver consistency and time consistency as constraints as in the majority of the ConVRP literature, they are included as objectives. Furthermore, instead of formulating a single weighted objective that relies on specifying relative priorities among objectives, an approach to approximate the Pareto frontier is developed. Specifically, an improved version of multi-directional local search is developed. The updated algorithm, IMDLS, makes use of LNS to find solutions that are improved based on at least one objective to add to the set of nondominated solutions at each iteration. The performance of IMDLS is compared with MDLS, three classic multi-objective algorithms, and two recent state-of-the-art multi-objective algorithms on a set of ConVRP test instances from the literature. The computational study validates the superior performance of IMDLS. Traditionally, travel distance has been the most common objective in routing problems in the literature. Although customer service considerations, such as driver and time consistency, are becoming increasingly important in certain industries such as small-package delivery and home health care, their associated costs are not well known. There are studies that quantify increases in travel distance when consistency considerations are implemented using hard constraints. However, depending on the application, strict enforcement of consistency considerations may not be required. The relative costs of "good but not perfect" consistency was not previously known. The results of our computational study suggest that pursuing the best compromise solution among all three objectives may increase travel costs by about 5% while improving driver and time consistency by approximately 60% and over 75% on average, when compared with a compromise solution having lowest overall travel distance.
Directions for future work can include studying consistency concerns in the context of other routing problem variants. For example, in the general ConVRP, customers require service on various days throughout a planning horizon, and the specific days each customer requires service are treated as problem input. The days each customer requires service could be modeled using decision variables instead, as in the Periodic VRP. Additionally, the problem studied in this article is static and deterministic. Real problems are often dynamic and/or stochastic in nature. Therefore, dynamic and/or stochastic variants of the ConVRP can be explored.

Notes on contributors
Ashlea Bennett Milburn is an Assistant Professor of Industrial Engineering at the University of Arkansas. She holds master's and doctoral degrees in Industrial and Systems Engineering from Virginia Tech and Georgia Tech, respectively. Her research focuses on applying operations research tools and techniques to large-scale logistics planning problems arising in the healthcare and humanitarian sectors. Much of her work focuses on issues in home health systems and disaster response.
Kunlei Lian received his B.S. and M.S. degrees in Industrial Engineering from Huazhong University of Science and Technology in 2009 and 2012, respectively. He is currently a Ph.D. candidate in the Department of Industrial Engineering at the University of Arkansas. His research interests include developing models and algorithms for various optimization problems arising in transportation, scheduling, logistics, and healthcare sectors. His teaching and research interests center on largescale optimization modeling and algorithms, especially their application in healthcare and energy. He is an award-winning teacher of those topics and co-author of numerous research papers and two comprehensive textbooks: a graduate text Discrete Optimization, published in 1988, and a comprehensive undergraduate/graduate textbook on mathematical programming, Optimization in Operations Research, which received the Institute of Industrial Engineers (IIE) Book of the Year award for 1998 and was recently updated in a second edition. Among his many other honors, he is a Fellow of both IIE and the Institute for Operations Research and the Management Sciences (INFORMS), as well as 2012 winner of the IIE's David F. Baker award for career research achievement.