{'fraction_non_alphanumeric': 0.027098138539848035, 'fraction_numerical': 0.015909705658150343, 'mean_word_length': 4.623340707964601, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A bi-criteria evolutionary algorithm for a constrained multi-depot vehicle routing problem

Most research about the vehicle routing problem (VRP) does not collectively address many of the constraints that real-world transportation companies have regarding route assignments. Consequently, our primary objective is to explore solutions for real-world VRPs with a heterogeneous fleet of vehicles, multi-depot subcontractors (drivers), and pickup/delivery time window and location constraints. We use a nested bi-criteria genetic algorithm (GA) to minimize the total time to complete all jobs with the fewest number of route drivers. Our model will explore the issue of weighting the objectives (total time vs. number of drivers) and provide Pareto front solutions that can be used to make decisions on a case-by-case basis. Three different real-world data sets were used to compare the results of our GA vs. transportation field experts’ job assignments. For the three data sets, all 21 Pareto efficient solutions yielded improved overall job completion times. In 57 % (12/21) of the cases, the Pareto efficient solutions also utilized fewer drivers than the field experts’ job allocation strategies.


Introduction and background of BSL
Amid significant shifts in the socioeconomic landscape, the growing complexity of global supply chains and innovative technological advances, it has become increasingly challenging for transportation companies to maintain their competitive advantage. Transportation-based companies are incessantly looking for ways to cut logistics costs, reduce waste, utilize information technologies, improve operational efficiency and increase overall productivity levels. Transportation companies are ultimately tasked with delivering packages efficiently and exceeding customer expectations.
After extensively researching logistics companies and existing methodologies for vehicle routing problems (VRP), it became apparent that the literature focused on more generalized problems that did not collectively address many of the variants that are common to transportation companies. Specifically, we perused the literature for research about solving VRPs that involved a heterogeneous fleet of vehicles, multi-depot subcontractors (drivers), and pickup/delivery time window and location constraints. Since we were unable to find existing literature that collectively addressed the aforementioned variants in a VRP, the initial phase of our research involved developing an evolutionary algorithm to minimize the total time to complete all jobs with the fewest number of route drivers given pickup/delivery times, locations, and vehicle capacity constraints (Lightner-Laws et al. 2016).
Our results revealed an interesting dichotomy between the objectives of minimizing the overall completion time vs. the number of route drivers needed to complete all jobs. There was an apparent trade-off between improving either the total job completion time or number of drivers. Consequently, in this phase of our research, our primary objective is to minimize this bi-criteria problem by exploring varying weights that prioritize the time vs. the number of drivers needed.
We plan to use a nested bi-criteria genetic algorithm (GA) that computes a Pareto frontier to explore this VRP, with a heterogeneous fleet of 4 vehicle types, multi-depot subcontractors (drivers), and pickup/delivery time window and location constraints. Our model will explore the issues of weighting the objectives and provide Pareto front solutions that can be used to make decisions on a case-by-case basis. Three data sets from BSL, a mid-sized transportation company, will be used in this paper. BSL has customarily used in-house field specialists to assign jobs to each driver based on the pickup/delivery times, location, number of packages and weight of the load. Although this is a common practice, management wants to automate their process and find more efficient methods to allocate jobs for route drivers as a means of improving operations and competitiveness. Ultimately, our model will provide the Pareto front solutions (route assignments) that BSL can use on a case-by-case basis to assign jobs to route drivers. This paper is organized as follows: Sect. 2 provides a review of the literature about VRPs. In Sect. 3, we describe our VRP and the associated constraints. Our GA framework and its specific components are presented in Sect. 4. The experimental results and discussion are given in Sect. 5. Finally, the conclusion and future research directions are presented in Sect. 6.

Literature review
The vehicle routing problem is a combinatorial problem that seeks to find the optimal route that minimizes the total travel distance required to deliver packages for a given set of customers. Typically customer demand is satisfied by a homogeneous fleet of vehicles from a central depot. Dantzig and Ramser (1959) capacitated vehicle routing problem (CVRP) is a VRP in which the homogeneous fleet of vehicles servicing the demand has a limited capacity which cannot be exceeded. The CVRP precipitated variant streams of research about multiple depots (Baldacci and Mingozzi 2009;Lau et al. 2010;Ombuki-Berman and Hanshar 2009;Vidal et al. 2011a), heterogeneous fleets (Baldacci and Mingozzi 2009;Brandao 2011;Choi and Tcha 2007), time windows for making deliveries (VRPTW), and designated pickup/delivery times (Dumas et al. 1991;Ropke and Cordeau 2009;Ropke et al. 2007;Baldacci et al. 2011).
A wide variety of closed form techniques have been used to solve VRPs. However, as more complex constraints are considered, finding an explicit optimal solution becomes computationally expensive and virtually impossible to ascertain; thus researchers have explored a variety of heuristics such as data mining (Chen et al. 2012), evolutionary algorithms (Vidal et al. 2011b;Lau et al. 2010), tabu searches (Cordeau and Maichberger 2011;Brandao 2011), graph theory (Likaj et al. 2013) andsimulated annealing (1993) to solve these complex transportation problems. All of these approaches seek to find near optimal solutions to address routing problem variants.
There has been extensive research on VRPs; however, Pisinger and Ropke (2007) are among the few who have developed a single robust heuristic that can be used to solve five different variants of this transportation problem. Their approach will solve the vehicle routing problem with time windows (VRPTW), capacitated vehicle routing problem (CVRP), multi-depot vehicle routing problem (MDVRP), site-dependent vehicle routing problem (SDVRP) and the open vehicle routing problem (OVRP). While this heuristic has been applied to a variety of VRPs individually, it was not applied to multiple variants simultaneously in a single problem instance (i.e., a CVRP with multi-depots, time window constraints and a heterogeneous fleet) or multi-objective problems. Research has been conducted on general multi-objective transportation route problems using linear programming parametrics (Aneja and Nair 1979;Gal 1975; and adjacent efficient methods (Evans and Steuer 1973;Yu and Zeleny 1974). Researchers have also used heuristic techniques and Pareto front solutions to explore bi-criteria transportation problems (Muller 2010;Prakash et al. 2014;Konak et al. 2006;Abounacer et al. 2012).
Recently, an increasing number of studies have emerged to solve a multi-criteria VRPTW. Tan et al. (2006) implemented a hybrid multi-objective evolutionary algorithm to minimize the total travel distance and the number of vehicles used to meet customer demand. They incorporated Pareto's optimality concepts for determining the best solutions to solve their multi-objective problem. Muller (2010) examined a VRPTW that utilized a homogeneous fleet of vehicles from a central depot. The objective of their research was to minimize overall costs and penalties for not adhering to soft time constraints. Their optimization problem had two competing objectives that made simultaneously optimizing both cost and penalties challenging. The dual objectives formed a Pareto front that the decision maker used to determine which solution best fits their needs. Their multi-objective optimization problem was solved using the ε-constraint method (Coello et al. 2002;Miettinen 1999), graph theory, Solomon's 11 insertion heuristic, ejection chain theory, or-opt and 2-opt procedures (Potvin and Rousseau 1995). Zou et al. (2013) proposed a hybrid swarm optimization model for minimizing the number of vehicles utilized to serve customer demand, the total travel distance and the total waiting times. Their model assumed that there were an unlimited number of service vehicles. Banos et al. (2013) presented a Pareto-based simulated annealing model for solving a VRPTW that minimized the travel distance and the imbalance (travel distance and vehicle loads) of the individual routes.
Most single and multi-objective VRPTW research has been tested using Solomon's benchmark data (1987). This data set consists of 56 different problems with varying fleet sizes, vehicle capacities, and extensive information regarding the location of customer demand and travel time/distances between different location sites. This data set has been used extensively to compare different heuristic methods for solving a VRPTW. The research presented in this paper includes the additional variants of using a heterogeneous fleet and multiple depots, with a multi-objective VRPTW. Although these additional variants are common for a variety of real-world delivery/logistics scenarios, an exhaustive exploration of the literature did not reveal existing research about a VRPTW with all of the variants presented in this paper. As a result, no standardized benchmark datasets are available to test and compare with the results of this research.
Lightner-Laws et al. (2016) developed a heuristic solution for a VRP where multiple variants are presented in a single problem instance. Their approach utilized a GA to solve a VRP with a heterogeneous fleet of vehicles, multi-depot subcontractors, and constraints on pickup/delivery times and locations. The solution aimed to minimize the overall driving time to complete all job assignments while utilizing the fewest number of drivers. While their approach addressed multiple variants of the classic VRP, it only marginally explored the trade-off between overall driving time and the total number of drivers.
This paper further develops the solution proposed by Lightner-Laws et al. (2016) to optimize a multi-objective VRP where multiple variants of the problem exist in a single problem instance. We develop a nested GA, which explores the VRP for a heterogeneous fleet of vehicles with multidepot subcontractors and constraints on pickup/delivery times and locations. We provide Pareto front solutions that allow BSL management to fully explore the trade-off between these two objectives.

Problem formulation
The multiple depot VRP presented in this research is concerned with the execution of a set of jobs where each job represents a package to be delivered. Each job is associated with a set of requirements that specify how the job is to be fulfilled. The following list gives the six requirements specified for each job: PT is the earliest possible time the package can be picked up, while DT is the latest possible time that the package can be delivered. Specialists use their expertise to assign the appropriate vehicle types (based on customer estimates of the weight, size, and number of packages) for each job. Although cars (C), SUVs (S), box trucks (B), and tractor trailers (T) are the vehicle types used in this problem formulation, alternative vehicle types could also be easily incorporated as well. Finally, JW is the sum of the weight of all packages for a given job.
All jobs, J i , (where i = 1. . .N and N is the total number of jobs) to be fulfilled are known in advance. A potential driver list is created for each job, J i , (based on customer input, estimated weight, size, and capacity constraints) from a set of available drivers, D k , (where k = 1. . .M and M is the total number of available drivers). A driver's job completion time is calculated from the time they leave their home and pickup their first job until the time their last job is delivered. The overall total time for all route deliveries is equivalent to the sum of the job completion times for each individual driver. Both mapping software and input from field experts are used to determine the travel time between any home, pickup or delivery location (Lightner-Laws et al. 2016).
BSL wishes to automate the process of allocating jobs to drivers and designating the order in which jobs should be picked up/delivered. The driver assignments should be made in a manner that minimizes total travel time and the total number of drivers required to complete all jobs. Holland (1975) developed a stochastic search technique, a genetic algorithm (GA), which incorporates features of natural selection to solve optimization problems. The main components in the design of a GA include:
A GA starts with a population of initial solutions. Genetic operations, (mutations and/or crossover events) alter initial parent solutions to create new children solutions. The next generation of solutions is selected from the set of parent and children solutions, in such a way that the best solutions have a higher probability of continuing to the next generation. The goal is to produce improved solutions at each generation by ultimately allowing natural selection to filter out the weaker candidates.
Our optimization model features a nested GA. Our main GA seeks to minimize the total distance traveled and the number of drivers required to complete all jobs. This genetic algorithm employs a secondary GA to help compute the fitness of a candidate solution. The secondary GA minimizes the total travel time for an individual driver by finding the optimal ordering of assigned pickups and deliveries. Our primary GA will be referred to as the Outer GA, and the secondary GA will be referred to as the Inner GA.

Candidate solution representation/initial population generation
Candidate solution representation refers to the encoding scheme that is used to represent a candidate solution. Suppose we were given the five jobs and potential driver list for completing each job shown in Table 1.
As mentioned in a previous section, available driver lists for each job are determined by BSL experts based upon the availability of BSL drivers with vehicle types that meet the established requirements. For this example, job J 1 can be serviced by driver D 1 or D 2 ; however, driver D 4 is the only driver available to service job J 5 . Figure 1 depicts our encoding scheme for three possible candidate solutions to this problem. Our encoding scheme includes the pickup and drop off locations for each job since optimal driver routes may include consecutive (nested) pickups and/or drop offs when multiple jobs are assigned to a single driver, in addition to the maximum amount of weight that the driver is carrying at any given time. In the first solution, driver D 2 was assigned job J 1 ; thus the pickup location J 1 PL and drop off location J 1 DL are assigned to driver D 2 . Since driver D 2 only carried job J 1 , the maximum weight that the driver carried at any given time is 300 lbs (the total package weight for job J 1 ). Driver D 3 was assigned jobs J 2 , J 3 , and J 4 , where they picked up and delivered each job consecutively, and the maximum weight that the driver carried at any given time was 850 lbs. Additionally, driver D 4 was assigned job J 5 (600 lbs). In the second solution, driver D 2 was assigned jobs J 1 and J 3 ; driver D 3 was assigned jobs J 2 and J 4 , and driver D 4 was assigned job J 5 . The last candidate solution assigned jobs J 1 , J 2 , J 3 , and J 4 to driver D 2 , and job J 5 to driver D 4 . The maximum weight carried by a driver is recorded at all times.
To create our initial population, we use the aforementioned encoding scheme to randomly assign job pickup and delivery locations to drivers, based upon the potential drivers list for each job. We impose large time penalties in the fitness function for each candidate solution if time conflicts make it impossible for the assigned driver to complete all jobs within the given time constraints, or if the known vehicle capacity was exceeded. These penalties effectively prevent infeasible candidates from being considered when other feasible candidates have been found.

Selection strategy
We adopt the tournament selection strategy to select candidate solutions for breeding. Our "tournament" is conducted by randomly selecting two candidate solutions from the current population, comparing their fitness values, and then declaring the strongest and weakest candidate of the pair (i.e. stronger candidates have high fitness values). Since weaker candidates can sometimes produce strong progenies, we set a selection probability of 0.9 that the strongest candidate solution from each pair will be selected for genetic operation.
Additionally, we employ the elitism strategy, where the best candidate solution from each generation is directly copied into the next population, unaltered. For our model, we select the top three candidate solutions to directly move on to the next generation. We then conduct tournaments, as described above, until our next generation is fully populated.

Genetic operator-mutation
Once a candidate solution is identified via the tournament selection process, our model utilizes mutation as the genetic operator. A candidate solution is mutated as follows:  1. We randomly select three current job assignments. 2. For each job selected, we randomly reassign the job to a different driver in the potential driver list for that job. If a selected job only has one potential driver in its list, no change is made.
For example, suppose we are given the potential drivers list provided in Table 1, and our tournament selection process produced Candidate Solution 1 (from Fig. 1). Our process for mutating Candidate Solution 1 is illustrated in Fig. 2. In step 1 of this figure, we see that jobs J 1 , J 3 , and J 5 are randomly selected to mutate. In step 2, we review the potential driver lists for the three jobs selected, and reassign each job to a different driver from their respective list. Since driver D 4 is the only available driver that can complete job J 5 , this job assignment remained unchanged.

Fitness evaluation
The fitness evaluation function mathematically expresses the value of a candidate solution. Genetic algorithms are primarily concerned with finding the best or "most fit" solutions for the final population. Accordingly, we start our fitness evaluation, by checking the vehicle weight capacities for all proposed driver assignment solutions. We immediately impose a large fitness penalty value to all solutions that violate capacity limits; thus reducing the likelihood of invalid solutions continuing to the next generation (Michalewicz 1995). For all valid solutions we focus on the multi-objective goal of minimizing the travel time and the number of drivers to complete all jobs. Like many multi-criteria problems, it is challenging to optimize both objectives simultaneously. In

Original Candidate Solution 1
Step 1: Randomly select 3 jobs to mutate Supposed jobs J 1 , J 3 , and J 5 Step 2a: Review the potential driver list for each job selected.
From Table 1: • Job J 1 could assigned to drivers D 1 or D 2 • Job J 3 could be assigned to drivers D 2 or D 3 • Job J 5 must be assigned to driver D 4 Step 2b: For each job selected, randomly reassigned the job to a different driver in its potential driver list • Job J 1 was assigned to driver D 2 . We reassign this job to driver D 1 .
• Job J 3 was assigned to driver or D 3 We reassign this job to driver D 2 • Job J 5 must remain with driver D 4 an effort to address potential conflicts, we aim to generate a Pareto optimal set consisting of solutions that are not dominated by either objective. For our problem, we define our first objective as minimizing the travel time and our second objective as minimizing the number of drivers. We use the weighted approach to formulate our multiobjective problem into the following single scalar objective function: where x represents a candidate solution and z i (x) represents the ith normalized objective function.
Our normalized first objective, z 1 (x) is In this equation f 1 (x) is the travel time required for all jobs to be completed using the driver allocation specified by candidate solution x. A separate GA model, which we refer to as our Inner GA (discussed below), is used to compute the minimum time for a driver to complete the assigned jobs. The Inner GA returns the time and optimal order that jobs should be picked up/delivered. Thus, for a single candidate solution, each driver assignment must be processed through our Inner GA, and the sum of the minimum times returned for all drivers will serve as f 1 (x) for the candidate solution. The function f * 1 is the minimal travel time to complete all jobs if we consider time as our only objective. We discuss f * 1 further below (Initial GA Model Parameters).
Conversely, our normalized second objective is computed as the total number of drivers that the candidate solution x allocated to complete all jobs, f 2 (x), divided by the minimal number of drivers, f * 2 . The denominator in this equation represents the minimal number of drivers to complete all jobs if we consider the total number of drivers as our only objective. We discuss f * 2 further below (Initial GA Model Parameters).
One of the challenges in using the weighted sum approach for addressing multi-objective problems is determining the appropriate weights to accurately reflect the priorities of the objectives. Since BSL was unable to definitively prioritize their goals, this issue of how to weight objectives surfaced in our research. Management asserted that determining the best trade-off between minimizing time and the total number of drivers is a judgment call that is made on a case-by-case basis. Accordingly, we choose to vary the weight in our multi-objective fitness function (Eq. 4.1.1) from 0 to 1, using increments of 0.1; the rationale is that varying the weight captures a spectrum of driver allocation solutions by gradually shifting the priority for each objective. We use this approach to determine the Pareto front solution set for our problem.

Initial GA model parameters
Before running our GA using our multi-objective fitness function (described in Eq. 4.1.1) we must determine f * 1 and f * 2 . f * 1 is determined by running our GA using total travel time as our fitness function and completely ignoring the number of drivers required to complete all jobs. The best time produced by our single objective GA, determines our f * 1 value for Eq. 4.1.1. We determine the value of f * 2 by running our GA using the single objective of minimizing the number of drivers needed to complete all jobs, without regard to the requisite travel distance. These parameters are then used in Eq. 4.1.1 to determine driver allocations for our multiobjective fitness function.

Termination condition
The Outer GA continues for a fixed number of generations or until no improvement in solution quality is seen. The final population is then analyzed to determine the solution with the minimum time and then the fewest number of total drivers.

Evolutionary algorithm for inner GA
The job assignments for a single driver are input into the Inner GA. The Inner GA is used to process all driver assignments and return the corresponding minimum time and ordering that jobs should be picked up/delivered. For example, let us consider the following assignment for Driver 2: Originally the driver is scheduled to pickup job J 1 , deliver job J 1 , pickup job J 2 , then deliver job J 2 . However, the minimum time solution may be to pickup job J 1 , pickup job J 2 , deliver job J 2 , then deliver job J 1 . The Inner GA would accept the original driver assignment above, and return the candidate solution below with its respective travel time: In computing the total travel time, we assume that drivers always begin at a known home location, before making their first pickup. The travel time ends after the last job is delivered. For a candidate solution to be considered valid, the pickup for each job must precede its delivery in the assignment ordering. Additionally, a substantial time penalty is added to any potential assignment that cannot (due to travel time conflicts) meet the required pickup and delivery times of all assigned jobs. This penalty significantly reduces the likelihood that these solutions are selected to continue on to later generations.

Candidate solution representation initial population generation
We retain the encoding scheme shown above for encoding candidate solutions for our Inner GA. The initial pool of candidate solutions for the first generation is created by taking the list of pickup and drop off locations for a driver and randomly shuffling their respective positions to form a new initial candidate solution. The generated solution is then checked to ensure that it is valid. Figure 3 shows a sample initial generation that could be generated from an original driver assignment. In this figure the notation D k H denotes the home location of driver D k .

Selection strategy
We employ the same selection strategy used for the Outer GA described above.

Genetic operator-mutation
Once a candidate solution is identified through our tournament selection process, our model utilizes mutation as the genetic operator. To be consistent with common GA notation, we refer to a candidate solution as a chromosome and we refer to a single item within a solution as a gene. A candidate chromosome is mutated as follows: • Step1: Randomly select one drop off or pickup location within the chromosome (i.e., randomly select one gene). • Step2: With a 0.5 probability, exchange the selected gene with its neighbor to the left or right. • Step3: Check the validity of new solution by making sure that following conditions are met: (a) the first position must remain the driver's home location, and (b) drop off location must come after the respective pickup location. If the solution is not valid, repeat steps 1 and 2 above up to three times, in attempt to find a valid mutated solution. After 3 attempts, if a valid solution has not been found, the original candidate solution will be copied unaltered into the next generation. Figure 4 shows a mutation process that yields an invalid candidate solution. For this instance, we must abandon the mutated solution and make up to two additional attempts to mutate the original candidate solution. If the future attempts yield a valid solution, the mutated candidate solution will continue to the next generation. If the future attempts fail to produce a valid mutated solution, the original candidate solution will move on to the next generation, unaltered. Figure  5 illustrates a mutation process that yields a valid candidate solution.

Fitness evaluation
The fitness value for a candidate solution is calculated by adding the travel times to reach each successive location in the solution. For example, consider the following candidate solution: Step 1: Select a gene to mutate. Suppose the gene J 1 DL was selected.
[D 1 H, J 1 PL, J 1 DL, J 2 PL, J 3 PL, J 2 DL, J 3 DL, J 4 PL, J 4 DL] Step 2: With a .5 probability, exchange the selected gene with its neighbor to the right or left. Suppose this probability process selected the left neighbor, J 1 PL. These two genes are exchanged.
This solution is not valid since the ordering has job J1 being dropped off before it was picked up. Therefore we must abandon this solution.

Selected Candidate Solution for Mutation:
[D 1 H, J 1 PL, J 1 DL, J 2 PL, J 3 PL, J 2 DL, J 3 DL, J 4 PL, J 4 DL] Step 1: Select a gene to mutate. Suppose the gene J 3 PL was selected.
[D 1 H, J 1 PL, J 1 DL, J 2 PL, J 3 PL, J 2 DL, J 3 DL, J 4 PL, J 4 DL] Step 2: With a .5 probability, exchange the selected gene with its neighbor to the right or left. Suppose this probability process selected the right neighbor, J 2 DL. These two genes are exchanged.
This solution is valid since the ordering has jobs being picked up before delivery. Thus this mutated candidate solution will move on to the next generation. Its fitness value would be the travel time from driver D 1 's home to job J 1 's pickup location, plus the travel time from job J 1 's pickup location to job J 1 's drop off location, plus the travel time from job J 1 's drop off location to job J 2 's pickup location, and so on.

Termination condition
The Inner GA continues for a fixed number of generations or until no improvement in solution quality is seen. Once the termination condition is met, the Inner GA returns the minimum time solution (and ordering) for the driver's job assignments.

Sample problem using BSL parameters
Initially, we used a small sample problem to test our GA and ensure that all six variants were properly incorporated into the model. We used BSLs Midwest metropolitan area service area for the pickup/delivery locations in the sample problem. To determine the time and distance required to complete jobs within the service area, 141 zones were established based on geographical locations. GPS software was then used to ascertain the distance and time between any two zones within the service area; finally, this information was used to calculate the overall time and distance required to complete all jobs. There were fifty drivers, with unique driver identification numbers (ID#), available to service customers in this sample problem. Specifically, there were 18 car drivers (ID# D1-D18), 16 SUV drivers (ID# D19-D34), 13 Box Truck drivers (ID# D35-D47) and 3 tractor trailers drivers (ID# D48-D50). It is worth noting that SUVs have a unique degree of flexibility because they could be used for jobs, which require cars or SUVs.
The sample input data are displayed in Table 2. Specifically, the table gives the earliest pickup times, the latest delivery times, the pickup/delivery location zones, the vehicle type required, maximum vehicle type weight capacity, and the job weights for 11 different jobs. A perusal of the vehicle types for jobs 2, 4 and 6 raises two questions: (1) why is an SUV required for Job 2 and the job weight is only 180 lbs? and (2) why is a box truck required for jobs 4 and 6, when the weights are only 1000 and 468 lbs, respectively? These jobs highlight the need for transportation field experts to be involved in the process of assigning vehicle types. Although these three jobs weigh considerably less than the maximum vehicle type weight capacity, other factors such as the length, width, surface area, fragility, density, materials handling requirements, etc. must also be considered when choosing the appropriate vehicle type; thus, factors other than simply the weight dictated that a larger vehicle was required to complete jobs 2, 4 and 6. Table 3 displays the four Pareto front solutions and specific coordinates (total travel time, number of drivers) for the sample problem. In this table, the MTPW is the maximum total package weight at any time and the VC is the maximum amount of cargo that the vehicle can carry (or the maximum vehicle type capacity). Let us examine the W = 0.3, Run 2, driver 19 data. For this solution, driver 19 uses an SUV to complete jobs 1, 7 and 8. Specifically, driver 19 would pickup job 8, pickup job 7, drop off job 7, drop off job 8, pickup job 1 and drop off job 1. Jobs 7 and 8 would be in the vehicle at the same time (MTPW = 525 lbs) and job 1 (MTPW = 280 lbs) would be in driver 19's SUV alone. The maximum total package weight at any time would be less than the maximum amount of cargo that an SUV can carry (VC = 2000). Figure 6 shows a graph of the 110 solutions found by our GA, under varying values of w. For each solution, the graph shows the travel time vs. the number of drivers required to complete all jobs. Our Pareto front solutions are depicted as square-shaped data points in the graph. The solutions in Table  3 and Fig. 6 encompass all six variants-pickup/delivery time windows, pickup/delivery locations and vehicle capacity constraints.

Results from GA vs. BSL job allocations
We conducted our experiments using actual BSL data from three different days. In Table 4, there is a complete list of all jobs that needed to be filled on three different days. This table also provides the earliest pickup times, the latest delivery times, the pickup/delivery locations, the type of vehicle needed for the job and the drivers that could potentially do the jobs. There were sixty-eight drivers, with unique driver identification numbers (ID#), available to service customers.
BSL does not currently preserve information about the estimated weight or package dimensions of fulfilled jobs in their records. Thus, for these experiments we assume that BSL vehicle type assignments are made correctly and hence have the available capacity to fulfill all assigned jobs over the course of a day. However, it is worth noting that our model is designed to ensure that vehicle capacity constraints are not violated when package weights are provided (refer to the sample problem in the previous section).
Preliminary experiments for both the inner and outer GAs, which systematically investigate alternative population sizes and maximum number of generations, were executed to determine efficacious GA parameter settings. Upon consideration of total computation time, overall job completion time, and marginal improvement of generated results, the population size and maximum number of generations were both set to 25. Prior to implementing our bi-criteria objective, the GA was run twice-once with the objective of minimizing the total travel time and once to minimize the number of drivers. The results from these two runs were used to determine f * 1 and f * 2 as described in Sect. 4 above. The minimum travel time, f * , 1 was determined to be 1260 for our Day 1 data set, 1673 for Day 2, and 2039 for Day 3. The minimum number of drivers, f * 2 , was determined to be 8 for our Day 1 data set, 6 for Day 2, and 15 for Day 3. Our multi-objective fitness function was investigated for w = [0, 0.1, 0.2, . . . , 1]. For each of the 3 data sets, the GA was run 10 times for each value of w (110 total runs per data set). All experimental results were performed using a Windows 7 Professional 64-bit Operating system with 6 GB of RAM and Inter(R) Xeon(R) CPU W3530 @ 2.8GHz.    For the Day 1 jobs, Fig. 7 shows a graph of the 110 solutions found by our GA, under varying values of w. For each solution, the graph shows the travel time vs. the number of drivers required to complete all jobs. Similarly, Figs. 8 and 9 show the solutions yielded from our GA for Days 2 and 3, respectively. Table 5 displays the specific coordinates (total travel time, number of drivers) for our 7 Pareto front solutions for Day 1, 6 Pareto front solutions for Day 2, and 8 Pareto front solutions for Day 3.
Tables 6, 7 and 8 show the specific job assignments for each Pareto solution, overall completion times and number of drivers for our GA for Days 1, 2 and 3. The final col-umn of these tables displays the job assignments that BSL experts gave to each driver; while the set of jobs are listed for each driver, the order that the jobs were picked up and delivered was not provided. Thus, we ran the job assignments that BSL provided through our program that determines the optimal route (i.e. the order for pickups/deliveries) for a set of jobs. We used these optimal routes to determine the minimum time that it would take the drivers to complete all orders. It is important to note that these times are based on mapping software's estimates of driving times between two locations and may not account for delays due to traffic or weather or alternative route selections made by individual drivers.  If we examine our Day 1 results, there are 7 Pareto optimal solutions with total job completion times ranging from 1243 to 2533 min as compared with 1571 min from the BSL job assignments. This represents an improvement in completion time of up to 21 %. Additionally, while 14 drivers were required for the BSL solution, only 8-13 drivers were needed for the Pareto solutions from the GA. Thus, the GA was able to provide automated solutions that would reduce both travel time and the number of drivers. Specifically, one of our Pareto solutions yielded 1414 min and 9 drivers. This solution reduced the travel time by 10 % and showed a 36 % decrease in the number of drivers when compared to BSL's manual driver allocations.
For Day 2, the Pareto front consists of 6 distinct points. The total job completion times for the Pareto solutions ranged from 1673 to 2479 min for the GA while the BSL solution was 1858 min. This represents time improvements of up to 10 %. Here, 12 drivers were required for the BSL solution in For Day 3, the 8 solutions in the Pareto front had job completion times spanning 2039-3883 min while utilizing 15-26 drivers. BSL's route allocation used 15 drivers who completed the jobs in 3698 min. The GA was able to improve the total completion time by up to 45 %. However, in providing this dramatic improvement in time, it utilized 11 more drivers than the BSL solution. One of our Pareto solutions yielded 2948 min with 18 drivers. With this solution, the travel time was reduced by 20 % with a 20 % increase in the number of drivers.
The results from the GA on Days 1 and 2 revealed that the overall job completion times were markedly better than the BSL solutions while the total number of drivers was equal or better. For the Day 3 results, the GA provides dramatically better completion times but utilizes more drivers in doing so. Overall, our model was able to successfully automate their driver assignment problem while providing options for marked improvement of their combined objectives.

Conclusion and future directions
In our collaborative efforts with BSL, we address their need to automate the route assignment process by developing a nested GA for a heterogeneous fleet of vehicles with multidepot subcontractors, pickup/delivery time windows, varying pickup/delivery locations, job weights and vehicle capac-ity constraints. BSL can input job requests for a given day and the model will yield the Pareto optimal solutions for their problem. Our model specifies the driver assignments, the order that jobs should be picked up/delivered, and the total time and number of drivers needed to complete all jobs. This allows BSL the ability to fully explore the trade-offs between driving time and number of drivers. Once a particular solution representing the best trade-off has been selected, the model will specify the driver assignments that should be designated.
Three different data sets from BSL were used to compare the results of the GA and BSL field experts. Our GA showed that the total job completion times were improved up to 21, 10 and 45 % on Days 1, 2, and 3 respectively. In addition, 12 of 13 Pareto front solutions for Days 1 and 2 reveal a decrease in the total number of drivers needed to complete all jobs. On Day 3, while results from the GA solutions yielded remarkably lower total job completion times, the BSL assignment yielded the fewest number of drivers. The Day 3 results indicate the need for us to continue to work with BSL to further understand how they prioritize both objectives (minimizing time vs. number of drivers required) on a daily basis. Overall, our GA results were promising and provided a substantive alternative to manually allocating job assignments.
In the future, we plan to modify our algorithm so that we can explore how unforeseen changes in the environment or disruptions in transportation can impact route assignments. For instance, if a driver has delivered some jobs and has a breakdown, we will explore which drivers can take over the remaining jobs most efficiently. We also intend to solve this vehicle routing problem using other meta- Table 6 Pareto front solutions from our GA and the BSL job Allocations for the Day 1 jobs W=0.9, Run 1 W=0.9, Run 9 W=0.8, Run 5 W=0.5, Run 2    heuristic algorithms, such as a Tabu Search and Simulated Annealing. We plan to compare these metaheuristics and refine our methods of minimizing both objectives simultaneously.