European Journal of Operational Research 115 (1999) 237±253
Solving the medium newspaper production/distribution problem Michael G. Van Buer a, David L. Woodru b
b,*
, Rick T. Olson
c
a Assimilation Initiatives, Inc., De Kalb, IL 60115, USA Graduate School of Management, University of California at Davis, Davis, CA 95616, USA c University of San Diego, San Diego, CA 92110, USA
Received 1 October 1996; accepted 1 May 1998
Abstract It is becoming increasingly important that the production and distribution of products be carefully coordinated. In this paper we study a problem from the newspaper industry where production and distribution are especially closely coupled since there can be no ®nished goods inventories. We describe the problem, give a mathematical formulation, and develop a solution strategy using heuristic search algorithms. Using data from a particular newspaper and extensive computational experiments, we ®nd that re-using trucks that have completed earlier routes is the most important way to achieve low-cost solutions. We also compare and contrast various heuristic search algorithms. Ó 1999 Elsevier Science B.V. All rights reserved. Keywords: Newspaper; Production; Distribution; Tabu search; Genetic algorithms
1. Introduction As businesses ®nd it increasingly necessary to become more ecient in both their operations and more eective the delivery of their goods and services to their customers, the eective integration of the production and distribution processes becomes imperative [3,5,12,13]. One of the strategies used to become more ecient in the production of goods is to minimize the amount of inventory ± both
* Corresponding author. Present address: TU Braunschweig, Abt. Allg. BWL, Wirtschaftsinfo. und Informationsmanagement, Abt-Jerusalem-Strasse 7, 38106 Braunscheig, Germany. E-mail: d-l.woodru@tu-bs.de
®nished goods inventory and work in process inventory. As these inventories are reduced, the production schedule becomes more closely linked with both the delivery schedule for inputs and the delivery schedule for ®nished goods. It is no longer sucient that the production schedule and the delivery schedule be approximately synchronized. There has been little work done in the area of reconciling the production and distribution schedules, see for example Blumfeld et al. [3]. An extreme example in which production and distribution are inseparable can be found in the newspaper industry. As the aphorism states, ``Today's news is tomorrow's history''. Consequently, having ®nished goods inventory available to de-couple production and distribution is not possible in this
0377-2217/99/$ ± see front matter Ó 1999 Elsevier Science B.V. All rights reserved. PII: S 0 3 7 7 - 2 2 1 7 ( 9 8 ) 0 0 3 0 0 - 2
238
M.G. Van Buer et al. / European Journal of Operational Research 115 (1999) 237±253
industry (Van Buer [20]). Because of this tight coupling of production and distribution in the newspaper industry, we will utilize the newspaper industry in our study of production/distribution problems. It is, in some sense, the limiting example of just-in-time production. 1.1. The newspaper industry The newspaper production and distribution problem is quite large and complex. A typical, large, metropolitan morning newspaper starts production of the newspaper shortly after midnight. The last newspaper must be delivered to home delivery subscribers by approximately 6:00 AM. In this 5±6 h span, close to 1,000,000 newspapers may be produced and delivered. Adding to the complexity of the problem is the fact that there may be more than 20 dierent editions of the newspaper having dierent editorial content and run of paper (ROP) advertisements, and many more dierent editions with dierent advertising inserts. Because the time to shift the inserting machinery from one insert to another is essentially negligible, i.e., set-up times are approximately zero, only the newspaper editions with dierent ROP content and editorial content are considered when determining a production schedule. In addition to the diculties posed by multiple editorial products, there is a natural con¯ict between the production and distribution functions. Taking a myopic view, the production department would tend to create a production schedule that would result in all newspapers being produced in the shortest amount of time. Since the production time, not counting the set-up times, is the same for any schedule, the production department would attempt to minimize the set-up times. Often, this objective would result in a production schedule in which the newspapers for the largest (in circulation) editions are produced ®rst. The largest editions tend to be nearest the city since the population density is higher. Since it is desirable to minimize set-ups, production would next choose to produce the editions that are editorially similar to the previous edition. A general rule of thumb is that editions close to one another geographically
are also close to one another in editorial content, so the newspapers would tend to be produced for the outlying subscribers last. The distribution department would prefer a dierent production sequence. Distribution requires that the newspaper be delivered by the deadline to all subscribers. The time that a newspaper is delivered is equal to: Production Time for the Edition + Production Time for all Previous Editions (with set-ups) + Stem Travel Time to News Carrier + Carrier Delivery Time. Notice that the second term is the only term that changes with dierent production schedules. The last two terms often result in the need to produce newspapers for the regions farthest from the production facility before those for regions closer to the production facility. These editions also tend to be most dierent from the editions closer to the production facility ± both in content and page count. Because of the dierence in page count, set-ups tend to be most time consuming going from the generally smaller distant editions to closer editions. Therefore, the production schedule that distribution would prefer to make certain delivery commitments are met is, in general, quite dierent than the one production would prefer. In addition to a production schedule, vehicle schedules and routes must be designed. Vehicles are loaded with newspapers at the production facility, then driven to either distribution centers or drop-o points where news carriers collect their newspapers. Drop-o points may be located at a street corner, driveway, or a front porch, and consequently, they require no direct capital investment. A single drop-o point may serve more than one carrier. Many large metropolitan newspapers utilize distribution centers where the news carriers come to receive their newspapers. Assembly of the newspapers may be completed at the distribution center by manually inserting any nontimely sections of the paper (such as Sunday supplements) that have been previously shipped to the distribution center. After all papers are fully assembled, they are bundled and distributed to the news carriers who come to the distribution center
M.G. Van Buer et al. / European Journal of Operational Research 115 (1999) 237±253
and collect their bundles, or the bundles are loaded onto small delivery vehicles, vans for instance, and delivered to drop-o points. Each of these vehicle ¯eets ± the ¯eet that serves the distribution centers and the ¯eet that delivers the newspaper bundles to the drop-o points ± must be scheduled and routed. In some situations, the locations and number of distribution centers and drop-o points must also be determined. Carrier routes may also be assigned to the distribution centers. Overall, the problem involves solving for the sequencing of production and vehicle departure times, the routing of vehicles, including the number of vehicle routes, the location of distribution centers and drop-o points, and allocating drop-o points to distribution centers and carrier routes to drop-o points. 2. The medium-sized problem We will consider the NDP for the medium-sized newspaper (m-NDP). This variant of the NDP contains most of the complexity of the NDP for the large newspaper (l-NDP). The salient dierences are in the circulation (fewer than 100,000 copies), number of editions (fewer than 10), and in the complexity of the distribution system. While a large newspaper may have 30±100 distribution centers, a medium newspaper will have at most 1 or 2 distribution centers. Often, a medium sized newspaper's only distribution center is the production facility; however, the main components of the NDP are still present in the m-NDP and experience gained by solving the m-NDP can be used to construct algorithms for the l-NDP. The m-NDP was previously addressed by Hurter and Van Buer [13] and Van Buer [20]. In these works, the authors describe the m-NDP, give a formulation, and develop a solution strategy. In this work, we both develop a new algorithm for this class of problem and extend the solution strategy to include re-using (recycling) trucks. When a truck is re-used, the expense of an additional truck and driver is saved. In order to save a truck, we may be willing to incur additional mileage costs. For instance, to allow a truck to return to the distribution center early enough to be
239
re-used, some routes that may be optimal for the scenario without recycling will have to be completed earlier by shortening them. Additional routes will then be needed to deliver to the remaining drop-o points so there will usually be more routes. Therefore, there will be at the very least additional route stems. Of course, there are other ways to reduce the number of trucks. One way is to utilize larger trucks. However, this strategy may result in no net savings since reducing the number of routes by utilizing larger trucks will result in longer routes. These longer routes may not be feasible because of the delivery time restrictions. Also, in order to utilize larger trucks eectively, some trucks may have to carry more than one newspaper product. Many newspapers require that trucks carry a homogeneous product in order to avoid any chance of delivering the wrong product to a customer. Therefore, recycling may be the only feasible method to reduce the number of trucks required to distribute the newspaper. In the m-NDP, the objective is to deliver the newspapers to all of the subscribers by the speci®ed time as inexpensively as possible. The costs that can be most readily adjusted in the short-term are the costs associated with the distribution of the newspapers, speci®cally, the number of vehicles and drivers required and the utilization of the vehicles (miles driven per day). The m-NPD can be verbally formulated as follows: min Cost of owning trucks (purchase, maintenance) + Truck Operating Costs s.t. Newspapers for a route are ready before the truck leaves, Each drop-o point is served exactly once, Last drop-o point served by deadline, Vehicle capacity limits not exceeded, Each carrier route assigned to a drop-o point.
(1) (2) (3) (4) (5)
Medium-sized newspapers typically utilize a news carrier force that consists of school age children.
240
M.G. Van Buer et al. / European Journal of Operational Research 115 (1999) 237±253
Consequently, the newspapers for a given carrier are most often delivered to the carrier's home. Since the location of the carrier's home is known, this feature eliminates the need to solve a location/ allocation problem for the drop-o points. Another result of utilizing children as carriers is that the routes are relatively small. Newspapers generally ensure that their routes can be completed in less than 2 h. This requirement results in routes that seldom have more than approximately 70 newspapers. In contrast, metropolitan newspapers rely almost exclusively on adult carriers using automobiles. These carrier routes may consist of more than 200 home subscribers. Finally, a typical medium-sized newspaper has only one printing press. Therefore, there is no need for us to allocate the printing jobs among presses. We also assume that the loading rate and the production rate are equal. In this formulation, we assume that there are an unlimited number of loading docks (operationally, unlimited loading docks merely means that there will never be a delay due to loading docks). From experience, we have found that there do not have to be many loading docks to prevent the loading dock from being a bottleneck. This observation can be explained from both an operational and ®nancial point of view. Operationally, since the loading rate is essentially the same as the production rate (we assume that they are equal), it is unlikely that there would ever be a shortage of loading docks (actually, if vehicles could instantaneously change positions, a single loading dock would be sucient). From a ®nancial point of view, loading docks are relatively inexpensive. Therefore, if the loading dock were the bottleneck, it would generally pay to install an additional one. Examining the verbal formulation of the mNDP reveals both a production sequencing problem with a deadline, and a form of the vehicle routing problem with time windows (VRPTW). To best address the m-NDP, these sub-problems must be solved concurrently. Bruno and Downey [4] showed that the question of whether or not a feasible solution exists for a production sequencing problem with set-ups and a deadline is NP-complete [6]. The production sequencing problem contained within the m-NDP has set-ups and a deadline on when production must be completed.
Since there is an overall deadline on the time when papers must be delivered to the customers, the production schedule must have a deadline. However, this deadline is unknown a priori since the routing schedule is not known. Therefore, the sequencing problem is NP-hard. The question of whether or not a feasible solution exists to the VRPTW is also NP-complete (Savelsbergh [18]). Because of the nature of the component problems, we will not attempt to solve the m-NDP optimally. In the appendix we give a mathematical programming formulation for the m-NDP to illustrate both its complexity and to more completely describe the problem. We reformulate for the heuristic search methods that we describe next. 3. Searching and sorting The solution technique must simultaneously address routing and production issues. Problems of this complexity can be addressed by special purpose heuristics or by general purpose heuristic search techniques such as simulated annealing (e.g., Aarts and Korst [1]), local search (e.g., Savage [17]), or tabu search (Glover [8,9]). From a practical standpoint, the advantage of heuristic search is the relative ease with which constraints or objective function terms can be added or modi®ed without modi®cation of the underlying search algorithms. The use of these techniques is aided considerably by the presence of a single vector solution representation. 3.1. Single vector representation Solutions for the m-NDP can be represented by a single vector, r, which gives the order of the drop-o points. We can evaluate the objective function value for a solution vector by assigning drop-o points in sequence to a truck until either the truck is full (the next drop-o point will result in exceeding the capacity), the next drop-o point would violate the time constraint, or the next dropo point in the sequence is for a dierent product. Computation of the time of delivery for drop-o point r(i) is simply the sum of the production times
M.G. Van Buer et al. / European Journal of Operational Research 115 (1999) 237±253
for all drop-o points up to and including i (plus set-up time) plus the total travel time for the truck to get to r(i). We may wish to maintain additional data structures to speed up the computation of the objective function value, but the solution is fully represented by the single vector r. In the mathematical programming formulation given in Appendix A, two types of variables were used to emphasize the separation between routing the trucks (the x variables) and sequencing the trucks (the p variables). The assumption of a homogenous truck ¯eet is critical to the use of a single solution vector. Without it, we would be forced to add a second dimension that indicates the type of truck used. Although it is not always the case, many medium-sized newspapers do have a truck ¯eet that can be reasonably modeled as homogenous. 3.2. Heuristic sort We can also utilize a heuristic sort of the data to aid in the development of good solution. We can sort the data in such a way that our initial solution vector might yield a reasonable solution, then apply heuristic search to try to improve on this initial solution. Computational experience (Van Buer [20], Hurter and Van Buer [13]) and analytical results (Van Buer et al. [19]) suggest that a good strategy to follow is to ship the newspapers to the most remote locations ®rst. Therefore, it is reasonable to produce the newspapers for the most remote regions ®rst. Also, in some sense, demand points that are geographically close together tend to be near each other in the routing sequence. Using this knowledge, we can sort the data in an intelligent manner. We place the product that has its demand area farthest from the production facility ®rst on the list. We could measure the distance of a large geographic region from the production facility by ®nding the centroid of the demand area. We are essentially developing lots. There is a one-to-one correspondence between lots and products. Then, within these lots we sort by the x-coordinate and then the y-coordinate. In this way, demand points that are physically close to each other, and hence may be close to each other
241
in a good routing sequence, are close to each other in the initial solution vector, r. 3.3. Neighborhoods Simulated annealing and tabu search are both based on local search which, in turn, is based on a neighborhood structure that de®nes moves from one solution to another. If a move based on neighborhood N can be made from r to s we say that s 2 N
r. In the m-NDP there is a hierarchical structure which we will exploit in de®ning neighborhoods and designing the solution algorithm. There are four distinct elements in the hierarchy: 1. production cluster 2. production sequence 3. truck cluster 4. truck route The production cluster pertains to the grouping of drop-o point demands into production lots. After these lots are formed, a production sequence (the sequence of the lots) must be determined. The drop-o point clusters must be amenable to routing, must not violate the vehicle capacity, and must not result in routes that take too long. Lastly, the drop-o points in these clusters must be routed. We will de®ne neighborhoods that can in¯uence these hierarchies singly and in combination. Given a permutation vector r that represents the current solution, there are several potentially interesting moves that can be used to transform it into a new solution. The most straightforward moves are referred to as full insertion moves. A drop-o point is moved from its current point in the sequence to a new point. We say that if source drop-o i, is inserted before destination drop-o j, (i < j) to transform r into r0 , then r0jÿ1 ri ; r0k rk1 for k 2 fi; . . . ; j ÿ 2g, and r0k rk for k 2 fi; . . . ; i ÿ 1; j; . . . ; ng. The case where i > j is similar. New neighborhoods are created by restricting the insertion moves to be within lot and within truck. Moves within a lot are the same as full insertion moves, except that for any source the destination must be in the same contiguous set of drop-o points requiring the same edition as the
242
M.G. Van Buer et al. / European Journal of Operational Research 115 (1999) 237±253
insertion source. More formally, for a solution vector r and a given source i, the destination j is valid for a within lot move if and only if same edition is required for ri , rj , and all drop-os rk with k between i and j. To put it another way, moves within a lot do not change the production sequence. Moves within a truck are more restrictive. When a solution is evaluated as described in Section 3.1 each drop-o point is assigned to a truck. To be a valid within-truck move, the move cannot result in the drop-o being assigned to a dierent truck. Moves that require the source to remain in the same truck are essentially TSP insertion moves. These moves cannot reduce the number of trucks needed, but they may reduce the mileage. It may be more useful to make use of moves that operate simultaneously on all locations that have been assigned to a truck. We can elicit this eect by using whole truck moves. This is repeated insertion move with very restricted sources and destinations. All of the drop-o points that have been assigned to one truck are inserted before the ®rst drop-o point for another truck. Note that for a feasible solution, these whole truck moves have little eect on the objective function value unless they move the locations to a point in the sequence that is outside the current production lot, i.e., the number of set-ups may change since a lot could be split, or two lots could form a single lot. The full insertion move results in the largest neighborhood. This move (or neighborhood) can in¯uence each level of the problem hierarchy. That is, full insertion moves will result in new production and truck clusters, and consequently, new production sequences and truck routes. However, for large data sets the resulting neighborhood size may be too large to allow the use of full insertion. The whole truck neighborhood is much smaller. This neighborhood aects all levels of the hierarchy except truck routes. 3.4. Steepest descent A steepest descent algorithm begins with an initial solution, x0 , and selects solutions at iteration k > 0 using the relation
xk argmin f^x x2N
xkÿ1
(a tie breaking rule may be needed). The function f^ may be the objective function, or it may be some other function that assigns lower values to ``better'' solutions. This function is referred to as the evaluation function. The descent ends when there are only higher cost solutions in the neighborhood of the current solution. Such a solution is referred to as a local minimum. Both tabu search and simulated annealing provide means for escaping local minima. A simpler way of escaping local minima is to restart the steepest descent at a random solution.
3.5. Simulated annealing Simulated annealing algorithms (e.g., see Aarts and Korst [1]) escape local minima by using randomness in move selection. At each iteration of a simulated annealing algorithm, a move is selected at random and the change in cost is computed for the move. In order to select moves at random from multiple neighborhoods, we ®rst select the neighborhood at random from among those de®ned and use that neighborhood for all moves at a temperature. When the temperature is reduced, the random selection of a neighborhood is repeated. An alternative is to exclusively use the full insertion neighborhood. Let Df^
r; s f^
s ÿ f^
r. If a move from r to s is considered, it is accepted (the move is made) if s improves on r
Df^ < 0. If s does not improve on r, it is still accepted with probability exp
ÿDf^
r; s=c; where c is the temperature. For high values of c, the acceptance probability is near one for all moves; hence, the algorithm proceeds from solution to solution at random. For very low values of c, only moves that result in improvement (a decrease in the cost function) are likely to be accepted. Simulated annealing algorithms begin with high values of c and slowly (over a large number of iterations) lower the temperature. The exact speci®cation of a sequence of c values is referred to as a cooling schedule. For a more detailed discussion of
M.G. Van Buer et al. / European Journal of Operational Research 115 (1999) 237±253
243
Table 1 Parameters used to control the cooling schedule in simulated annealing IN I T PR O B TE M P FA C T O R SI Z E FA C T O R MI N PE R C E N T
The initial temperature is set so that the probability of accepting moves is approximately IN I T PR O B . Multiplying the current temperature reduces the temperature by TE M P FA C T O R . The number of iterations at each temperature (both during InitializeTemp and cooling) is SI Z E FA C T O R times the neighborhood size. The algorithm terminates when ®ve temperatures in a row result in less than MI N PE R C E N T percent acceptance of moves.
simulated annealing and the eects of the cooling schedule, see Johnson et al. [14]. We will make use of the cooling schedule that they specify. Its parameters are given in Table 1. Using certain cooling schedules, simulated annealing algorithms can be shown to converge in probability to optimal solutions (see, e.g., Hajek [11] or Lundy and Mees [15]). These results are valid only asymptotically. Cooling schedules used in practice often dier signi®cantly from these schedules primarily because convergence in probability of an algorithm is of no practical value when much simpler theory guarantees that the optimal solution will be visited at least once in a ®nite (albeit unknown and potentially very large) number of iterations with probability one. 3.6. Tabu search Tabu search algorithms select moves according to the steepest descent scheme except they make use of a tabu list to force the search away from solutions visited in recent iterations. This tends to help the search to escape from local minima. Moves are rejected if they satisfy conditions given by the tabu list. The list can be thought of as a FIFO list based on certain attributes of the most recent j moves. Solution vectors that are the result of moves having the attributes found on the tabu list are removed from consideration unless they meet an aspiration criteria used to avoid removing very good moves from consideration. (Since attributes of moves are tabu, rather than recent solutions, many trial solutions may be forbidden that have not been recently visited.) For a more detailed introduction to tabu search and a discussion of some applications, see Glover and Laguna [10]. For a more thorough discussion and examples of
possible move attributes that can be used in a particular application, see Malek et al. [16]. For insertion moves, there are a number of possible choices for the tabu restrictions. Suppose location r(i) is inserted before r(j) to create the new solution s. Note that s(j) r(i). The following possibilities are listed in decreasing order of restrictiveness: 1. Forbid moves involving r(i) or r(j). 2. Forbid insertion of r(i). 3. Forbid insertions of r(i) that move i in the opposite direction (e.g. if i < j forbid all insertions of r(i) to be in front of j0 for j0 < j). 4. Forbid insertions of r(i) to r(i + 1). Intuitively, we would expect to use a shorter tabu list for more restrictive options. The algorithm we will be using adjusts its tabu list length dynamically and automatically. This explains why our preliminary experiments reveal no signi®cant difference in performance between the ®rst two choices. The second two proved to be overly permissive and resulted in a moderate degradation in performance. The second choice was the one used by Gendrau et al. [7] in their application of tabu search to the vehicle routing problem. In addition to the tabu list, many tabu search applications make use of long-term memory to force the search into new regions of the search space. This is often referred to as diversi®cation. Intermediate-term memory is sometimes used to intensify the search in the areas of the search graph where good solutions have been found or are likely to be found. We refer to the tabu search algorithm with full insertion neighborhoods as T S . We will also make use of the reactive tabu search (Battiti and Tecchiolli [2]). It is fully general and has a parameterization that can be applied across a broad range of problems and instances. If all moves are tabu, the algorithm makes a random
244
M.G. Van Buer et al. / European Journal of Operational Research 115 (1999) 237±253
move and reduces the size of the tabu list. Longterm memory is accomplished by hashing solutions (Woodru and Zemel [21], Battiti and Tecchiolli [2]). When a solution is repeated, the tabu list length is increased. Furthermore, when CH A O S dierent solutions have been repeated RE P times or more, an attempt is made to radically change the solution without completely disrupting the current solution structure. Good, robust values for the parameters of this algorithm are RE P CH A O S 2 and initial tabu list length solution vector length divided by 4. It should be noted that Battiti and Tecchiolli made use of RE P CH A O S 3 but since our escape mechanism is much less disruptive than theirs, it is intuitive that smaller values of these parameters would be preferred. This need to radically change the solution gives us a nice opportunity to introduce multiple neighborhoods. Battiti and Tecchiolli escape by making a number of random moves, but for this problem we ®nd it useful to invoke new neighborhoods. Every other time a radical change is indicated, we revert to insertions within trucks. When a radical change is indicated while the search is being conducted using insertions within trucks, the search is changed to one of the other neighborhoods. The neighborhood to be switched to is simply the next on a circular list of the remaining neighborhoods. Our scheme implies that when the search is in one of these neighborhoods and a radical change is indicated, the radical change is accomplished by switching the search back to the within trucks neighborhood. We refer to this algorithm as R T S . 4. Recycled trucks In this section we consider trucks that ®nish their routes early enough in order to service a subsequent route. Trucks that service more than one route will be referred to as recycled trucks. It is expected that recycling trucks may result in slightly higher mileage charges since there may be an extra trip to and from the production facility. However, initial mileage charges will be quickly compensated for by the savings resulting from
reduced ¯eet size. Additionally, allowing the recycling of trucks will encourage some of the earlier routes to depart the production facility to be shorter than without recycling and to service less demand ± reducing the time for production and distribution for these routes. Routes most likely will also be encouraged to have their last drop-o point closer to the production facility than when recycling is not allowed. This change will come about because the return stem travel time will be important in determining whether or not a truck will be eligible for recycling. A vehicle is a candidate for recycling if it is able to return to the production facility prior to the start of production for a subsequent route. Using the notation from the m-NDP formulation given in Appendix A, tk Clast drop 0 6 Starting time of some late route the 0th drop-o point is the production facility. We expect the size of the routes to be no greater than in the previous version. Also, if we do not force recycling, the objective function can be no worse than without recycling. Obviously, we choose not to recycle if the objective function increases. In the previous variant, production and distribution were only tied together by the set-up times. In the recycling variant, there is even a tighter linkage. In order for a vehicle to be eligible for recycling it must arrive back at the production facility before production commences for some route. It can then take its place in the loading dock. 5. Computational experiments We compare the solutions obtained with different algorithms and dierent starting points in this section. Besides evaluating how allowing recycled trucks aect the total distribution cost, we compare the relative performance of the heuristic search methods. The data set we are using was initially obtained from a medium-sized Midwestern newspaper, however, the data has been somewhat modi®ed. For instance, we have estimated some vehicle costs, mileage costs, and set-up times. The demands and locations are those given to us by the newspaper. We use the l1 norm to estimate
M.G. Van Buer et al. / European Journal of Operational Research 115 (1999) 237±253
distances. The data used in these experiments is available on request. Table 2 summarizes the algorithms tested. All algorithms were written in C++ using a common library of class structures and functions. The tests were run on a 166 MHz Pentium PC. To reduce the neighborhood sizes when running T S , W H T and R T S , insertion moves were restricted so that only moves of up to one ®fth the sequence length were considered. S A was varied by altering the parameters given in Table 1. In preliminary studies we veri®ed that performance varied with the parameters as predicted by Johnson et al. [14]. In addition to the parameters they recommend (0.4, 0.95, 16, 2), we also considered two other sets of cooling schedule parameters. Four variants of the m-NDP were solved with each algorithm exploring all the combinations of with or without truck recycling, and with or without using the heuristic sort to derive a good initial solution. Where the heuristic sort was not used, ten initial starting solutions were generated and all six algorithms were applied to these solutions. Where the heuristic sort was used, T S , W H T , and R T S were only applied once because they are deterministic algorithms once a starting point and tabu procedure have been de®ned. Simulated annealing was applied to the initial schedule derived from the heuristic sort ten times because of the probabilistic component associated with accepting or rejecting proposed moves. T S , W H T , and R T S were allowed to run for 1 h. Simulated annealing was allowed to run until the termination criteria were met. Depending on the particular cooling schedule parameters, initial solution and random number seed used took from ten minutes to one hour. Tables 3±6 give summaries of the runs. Tables 3 and 4 consider the scenario where recycling
is not allowed. The eect of recycling is seen in Tables 5 and 6. Tables 3 and 5 consider the case where random initial solutions are used. These tables show how the solutions found with each algorithm improved over time. An entry in a column for an S A algorithm does not necessarily mean that the algorithm required the full amount of time. For example, when S A (0.4, 0.95, 8, 2) was applied without recycled trucks to randomly generated starting solutions (Table 3), the algorithm actually converged in approximately 2500 s. The values in the 3600 s column should be interpreted to mean that if an analyst had one hour available, and this algorithm were applied, the average result found would be approximately 451.33 and the best found in our 10 trials was 451.06. Table 7 provides additional information that can be used to see how the algorithms performed. In this table, the typical eort required for each algorithm to ®nd the best solution in a particular run is summarized. For example, when T S was applied to the heuristically derived initial solution and recycling was not permitted, the best solution (450.08) was found after 46.70 s and 449,028 different schedules were considered. No further improvement was found during the remaining 59 min (and 34 million additional function evaluations). On the other hand, S A (0.4, 0.975, 32, 2) required an average of 3500 s and 32,000,000 function evaluations to converge to solutions that averaged 455.19. One thing is very clear, the use or non-use of recycling is much more important than the choice between the better performing search algorithms. Nevertheless, we can make some interesting statements about the search methods. We have conducted numerous additional computation experiments on other data sets and the algorithmic conclusions are qualitatively the same.
Table 2 Summary of algorithms used in computational experiments TS WHT RTS S A (0.4,
0.975, 16, 2) 0.95, 8, 2) S A (0.4, 0.975, 32, 2) S A (0.4,
245
Tabu search with full insertion neighborhoods Tabu search with whole truck moves Reactive tabu search with hierarchical neighborhoods Simulated Annealing w/default Johnson et al. cooling parameters Simulated Annealing w/modi®ed TE M P FA C T O R and SI Z E FA C T O R Simulated Annealing w/modi®ed TE M P FA C T O R and SI Z E FA C T O R
246
M.G. Van Buer et al. / European Journal of Operational Research 115 (1999) 237±253
Table 3 Summary of results when recycling is not allowed and random initial solutions are used Time (sec)
Algorithm
Data
10
60
600
TS
# of trials Minimum objective Average objective Std. dev.
10 537.56 603.29 58.74
10 536.72 580.35 51.42
10 498.05 563.42 45.83
10 495.44 527.94 40.13
WHT
# of trials Minimum objective Average objective Std. dev.
10 500.05 525.32 32.25
10 450.86 460.00 16.47
10 450.44 455.51 12.58
10 450.44 451.40 0.84
RTS
# of trials Minimum objective Average objective Std. dev.
10 500.30 526.13 31.83
10 450.99 452.61 1.35
10 450.99 451.50 0.33
10 450.99 451.48 0.32
3600
S A (0.4,
0.975, 16, 2)
# of trials Minimum objective Average objective Std. dev.
10 473.44 493.28 16.61
10 473.44 481.04 5.16
10 468.41 470.68 1.83
10 450.35 450.93 0.37
S A (0.4,
0.95, 8, 2)
# of trials Minimum objective Average objective Std. dev.
10 474.26 493.45 16.42
10 469.00 484.59 13.21
10 451.52 451.82 0.34
10 451.06 451.33 0.22
S A (0.4,
0.975, 32, 2)
# of trials Minimum objective Average objective Std. dev.
10 471.38 494.73 18.47
10 468.19 481.15 11.31
10 464.96 472.30 4.20
10 454.38 455.37 0.52
· Using the initial starting solution generally made little dierence in the quality of the ®nal solution obtained. It had no eect on the speed with which the S A algorithms found good solutions. It improved the speed with which the tabu searches found good solutions to the problem without recycling, but had mixed results when recycling was allowed. · The S A parameters (0.4, 0.95, 8, 2) were the most eective. In general, the solutions were comparable to the other sets of parameters, but it converged more quickly because it had the smallest value for SI Z E FA C T O R , so it performed fewer (generally nonproductive) iterations at each temperature. · The R T S and S A using the schedule (0.4, 0.975, 16, 2) were able to ®nd some solutions to the problem with recycling that were much better than the other solutions found, but they were
not able to do it consistently. In 93 total replicates of all algorithms and starting points there were three solutions of approximately 132. The other 90 solutions had a cost of at least 169. The low cost solutions required only two trucks (at a cost of $40/truck), while most solutions found required three. A careful examination of the best three-truck solution and the best twotruck solution revealed that there were many dierences in the solutions. Not only were production sequences dierent, but the routes were generally comprised of dierent drop-o locations. Clearly the eort required to ®nd a twotruck solution is more than a simple transformation of a three-truck solution. · Except for these three instances, W H T , R T S and the S A runs gave very consistent results. The average results obtained with T S were somewhat less consistent as their higher standard deviation
M.G. Van Buer et al. / European Journal of Operational Research 115 (1999) 237±253
247
Table 4 Summary of results when recycling is not allowed and the heuristic sort procedure is used to generate the initial starting solution Time (sec)
Algorithm
Data
10
60
600
3600
TS
# of trials Minimum objective Average objective
1 450.11 450.11
1 450.08 450.08
1 450.08 450.08
1 450.08 450.08
WHT
# of trials Minimum objective Average objective
1 451.94 451.94
1 450.85 450.85
1 450.85 450.85
1 450.85 450.85
RTS
# of Trials Minimum objective Average objective
1 452.92 452.92
1 451.17 451.17
1 451.17 451.17
1 451.17 451.17
S A (0.4,
0.975, 16, 2)
# of trials Minimum objective Average objective Std. dev.
10 460.19 473.06 9.82
10 460.19 472.15 9.74
10 460.19 468.54 4.15
10 450.75 451.13 0.31
S A (0.4,
0.95, 8, 2)
# of trials Minimum objective Average objective Std. dev.
10 459.75 476.74 13.05
10 459.75 474.15 11.07
10 451.02 451.69 0.53
10 449.92 451.07 0.72
S A (0.4,
0.975, 32, 2)
# of trials Minimum objective Average objective Std. dev.
10 459.83 471.52 7.01
10 459.83 470.39 6.49
10 459.83 469.91 5.85
10 454.42 455.19 0.54
indicates. Low variance among replicates can be important when the search is being done as part of sensitivity analysis (e.g., trying trucks with diering capacities). It is interesting to compare solutions from different truck use policies. In Table 8 we present data concerning the best solutions that we found when recycling is not allowed, and when it is allowed. These solutions dier in several signi®cant ways. First, as expected, recycling permitted a solution using fewer trucks, but also more routes. In fact, the solution with recycling permitted revealed three additional ``routes'' which consisted of ``deliveries'' to the production plant. This occurs when the ®nal drop-o location for some group of papers is the plant. In this case a truck is not actually needed. In the subsequent discussion these routes have been disregarded. Also as expected, the total route length for the recycle strategy is greatest. This observation can be explained by the fact that there are more routes with the recycle strategy and consequently there is
more stem travel distance. In addition to there being more routes in the recycle strategy, there are more set-ups. Clearly, trade-os are being made between total travel distance, production time, and the number of vehicles. Recall that vehicles and drivers are the most expensive component of the distribution system. Also note that the average route is shorter when vehicles are recycled. This observation is consistent with the previous comments. It is also worthwhile to examine the average return stem length. When recycling vehicles, the vehicles need to return to the production facility before production begins for the next route that the vehicle will service. Consequently, it is bene®cial if the vehicle can complete its previous route relatively close to the production facility. Notice that the average return stem distance is actually shorter when recycling is not permitted than when it is permitted. The return stem distance when trucks are recycled can be further analyzed. When the trucks return from 9 of the 11 routes
248
M.G. Van Buer et al. / European Journal of Operational Research 115 (1999) 237±253
Table 5 Summary of results when recycling is allowed and random initial solutions are used Time (sec)
Algorithm
Data
10
60
600
3600
TS
# of trials Minimum objective Average objective Std. dev.
10 186.19 215.55 23.67
10 176.88 182.79 6.28
10 175.74 179.56 4.87
10 173.78 177.29 2.95
WHT
# of trials Minimum objective Average objective Std. dev.
10 189.33 223.01 23.18
10 170.63 171.22 0.73
10 170.13 170.77 0.65
10 170.10 170.49 0.41
RTS
# of trials Minimum objective Average objective Std. dev.
10 189.33 223.02 23.19
10 170.96 171.53 0.59
10 170.39 171.11 0.40
10 170.31 170.59 0.26
S A (0.4,
0.975, 16, 2)
# of trials Minimum objective Average objective Std. dev.
10 262.16 269.54 4.74
10 181.66 185.07 1.79
10 181.66 184.17 1.72
10 131.56 166.88 12.50
S A (0.4,
0.95, 8, 2)
# of trials Minimum objective Average objective Std. dev.
10 260.65 269.46 4.80
10 183.36 184.47 0.99
10 174.51 181.53 4.67
10 169.98 170.29 0.22
S A (0.4,
0.975, 32, 2)
# of trials Minimum objective Average objective Std. dev.
10 258.82 268.73 5.66
10 197.05 256.93 21.25
10 182.46 183.82 1.05
10 178.57 181.40 2.34
they are returning to pick up more papers to deliver on another route. The average return stem for these trips was 5.34. The last routes performed by the trucks had an average return stem of 13.31 miles indicating that the last deliveries on these routes were relatively far from the production plant. By saving these deliveries until last, the trucks were able to complete their early routes more quickly returning to the plant for more papers that might otherwise have required another truck for delivery. 6. Conclusions and directions for further research In this paper we have described and formulated an important problem that links the production and distribution of newspapers having strict delivery deadlines and the potential for setup times
between editorial products. We have demonstrated that recycling the delivery trucks can have a dramatic impact on costs. We have shown that there are several neighborhood structures that can support heuristic search. The number of potential search algorithms and parameterization is in®nite. The ``best'' heuristic search algorithm probably will never be found. However, we have shown that good solutions can be found. Although the decision concerning recycling trucks clearly has the most bearing on the total distribution cost, we were able to reach a number of insights by comparing and contrasting various algorithms. This paper makes no eort to contribute to the debate concerning the best meta-strategy for heuristic search. The tabu search implementation has internally low variance. No doubt other T S implementation could improve upon it. In fact, a variant of this problem has been
M.G. Van Buer et al. / European Journal of Operational Research 115 (1999) 237±253
249
Table 6 Summary of results when recycling is allowed and the heuristic sort procedure is used to generate the initial starting solution Time (sec)
Algorithm
Data
10
60
600
3600
TS
# of trials Minimum objective Average objective
1 171.14 171.14
1 170.60 170.60
1 170.32 170.32
1 170.32 170.32
WHT
# of trials Minimum objective Average objective
1 172.59 172.59
1 171.44 171.44
1 171.03 171.03
1 170.63 170.63
RTS
# of trials Minimum objective Average objective
1 172.58 172.58
1 171.91 171.91
1 171.20 171.20
1 133.29 133.29
S A (0.4,
0.975, 16, 2)
# of trials Minimum objective Average objective Std. dev.
10 201.01 223.10 22.95
10 181.77 184.23 1.55
10 181.77 183.94 1.26
10 131.40 168.60 13.16
S A (0.4,
0.95, 8, 2)
# of trials Minimum objective Average objective Std. dev.
10 190.53 212.62 25.13
10 182.33 184.89 1.40
10 175.01 183.01 4.26
10 169.89 170.10 0.13
S A (0.4,
0.975, 32, 2)
# of trials Minimum objective Average objective Std. dev.
10 198.19 218.55 23.58
10 198.19 218.55 23.58
10 183.10 183.99 0.52
10 178.92 181.47 2.28
used by Woodru to illustrate new tabu search methodologies (Woodru [22]). Improved search and sort algorithms remain as a topic for future research. The heuristic sort we provided is valuable if trucks are not going to be recycled. Perhaps a similar fast sort can be developed for use when the trucks are going to be recycled. These problems are solved infrequently but the eects are felt daily, so large amounts of computer time can reasonably be employed. On the other hand, improved search heuristic methods that provide faster access to good solutions would be valuable if the search is to be embedded in a decision support system to explore the impact of diering truck costs and capacities. Another area of potential improvement is the use of a TSP algorithm to post-optimize the tours that are generated. A much more fertile area involves recognition that solutions can be improved by removing the requirement that vehicles carry a single product. Removing this requirement will result in a more complex problem since more production
sequence and route sequences exist. However, by allowing heterogeneous vehicle loads, we can more fully explore the trade-o between the number of vehicles and set-ups and routing distances. In practice, the main objection to mixed loads is the additional complexity in loading and unloading. This increased complexity could be captured by assigning an extra cost to a mixed route. However, unlike the other costs required by our model, this one would be dicult to estimate. An alternative would be an expert system that would evaluate trade-os between operating costs and complexity of operating with the schedule. Another important issue that must be explored is the interaction between production and distribution capacity. Our model explicitly considers the cost of distribution capacity but treats production capacity only in the constraints. The lowest cost recycling solutions that we found had the highest number of setups indicating in some sense that there is an excess of production capacity. However, had there been less excess capacity, it is likely that
250
M.G. Van Buer et al. / European Journal of Operational Research 115 (1999) 237±253
Table 7 Summary of eort required for each algorithm to discover best solutions of a run No recycling, random start
No recycling, heuristic start
With recycling, random start
With recycling, heuristic start
TS
# of trials Avg. solution Time to best sol'n Function calls to best sol'n
10 527.94 2014.15 17,681,400
1 450.08 46.70 449,028
10 177.29 1933.86 13,389,996
1 170.32 289.76 2,408,771
WHT
# of trials Avg. solution Time to best sol'n Function calls to best sol'n
10 451.40 687.95 5,179,673
1 450.85 53.92 415,226
10 170.49 613.87 3,911,860
1 170.63 1539.57 9,892,323
RTS
# of trials Avg. solution Time to best sol'n Function calls to best sol'n
10 451.48 774.29 5,861,841
1 451.17 20.01 153,676
10 170.59 1477.84 9,344,043
1 133.29 1470.18 9,306,965
S A (0.4,
0.975, 16, 2)
# of trials Avg. solution Time to best sol'n Function calls to best sol'n
10 450.93 2547.06 22,791,775
10 451.13 2570.54 23,060,869
10 166.88 3238.97 18,409,078
10 168.60 3434.80 17,799,488
S A (0.4,
0.95, 8, 2)
# of trials Avg. solution Time to best sol'n Function calls to best sol'n
10 451.33 659.94 5,877,393
10 451.07 671.79 5,994,745
10 170.29 1049.03 5,718,738
10 170.10 1087.76 5,904,217
S A (0.4,
0.975, 32, 2)
# of trials Avg. solution Time to best sol'n Function calls to best sol'n
10 455.37 3515.41 31,070,514
10 455.19 3504.44 31,116,076
10 181.40 2144.88 10,555,258
10 181.47 2125.32 10,409,900
Table 8 Comparison of solutions with and without recycling Strategy
Optimal cost
# Trucks
# Routes
Total miles
Avg. route time (h)
Avg. return stem (min)
# Setups
w/o Recycle w/ Recycle
449.92 131.04
10 2
10 11
207.98 214.18
0.44 0.38
5.34 6.79
8 14
more vehicles would have been required to distribute the newspaper since fewer set-ups could be incurred. Since production capacity for newspapers is relatively lumpy (i.e., capacity is essentially purchased by acquiring additional presses so that capacity is generally added in rather large increments) there will typically be considerable excess capacity for some time after the purchase of additional capacity. Consequently, while there is considerable excess capacity, the newspaper can utilize a lower cost distribution system utilizing relatively few vehicles and drivers and relatively many set-
ups. As the excess capacity is exhausted, there can be a trade-o between vehicles and set-ups, possibly eventually resulting in a solution that does not utilize any recycling. This capacity interaction could be modeled by considering press capacity to be a ®rst-stage decision and production/distribution sequences to be a second-stage decision enabling the newspaper to consider the cost reduction of the distribution system when exploring the question of production capacity expansion. As businesses increasingly ®nd it necessary to become both more ecient in their operations and
M.G. Van Buer et al. / European Journal of Operational Research 115 (1999) 237±253
more eective in the delivery of their goods and services to their customers, the eective integration of the production and distribution processes becomes more imperative. Our hope is that this paper contributes to research in support of this goal. We have studied an example that constitutes an extreme form of just-in-time production, but the concepts and algorithms apply to every tightly coupled production/distribution system where the production is centralized and planning must be done for distribution that is along pre-determined routes. Appendix A. Mathematical programming formulation for the m-NDP We start by de®ning the notation to be used in the formulation of the m-NDP. First, the parameters are de®ned and then the decision variables. Parameters The demand at drop-o j for edition p Djp R The production rate of the presses in papers per minute The set up time required to shift Spq between editions p and q The time to travel from drop-o point i Cij drp-o point j including unloading time at the drop-o point K The set of available vechicles, all are identical N The set of drop-o points Q The capacity of the vechicle 1 if drop-o j receives edition p Zjp 0 otherwise P The number of editions T The delivery deadline V Converts the net present value of a truck and driver to equal the net present value of operating a truck for one time unit Pevery day. We assume p Zjp 1 for all j. Variables p A permutation vector giving the order of routes (decision)
lk tk xijk ypk
251
The time required to produce the newspapers for route k (implied) The time that the last bundle is delivered on route k (implied) A binary variable, 1 if vehicle k goes directly from i to j (decision) 0 otherwise 1 if route k delivers edition p (implied) 0 otherwise
A decision maker must know which drop-o points are served by which vehicle route and the order in which the vehicles depart the production facility. The p(k) variables are used to determine the order in which the dierent routes are served, i.e. these variables give the production schedule. The argument, k, for p(k) is the position in the sequence that a particular vehicle occupies. For example, p(1) 3 means that vehicle 3 departs the production facility ®rst. In our formulation of the m-NDP, the permutation vector serves to join the production sequencing problem and the vehicle routing problem contained in the m-NDP. For notational convenience, we assign a value of zero to p(0) and to yp0 for all p. Assuming the production characteristics are given, the decisions that must be made in the design and operation of a production-distribution system and the decision variables corresponding to each decision in the model are: 1. Which drop-o points are on which route? Clustering drop-o points, xijk . 2. How are these drop-o points permuted? Routing of the clusters, xijk . 3. In what order are the routes serviced? Sequencing of production and routes, p(k). 4. In what sequence are the production lots produced? Answered by xijk and p(k). 5. What edition is delivered by vehicle k? Clustering and sequencing, ypk . The production±distribution function of a newspaper is often viewed as a cost center. No revenues are directly attributed to it. In fact, if the delivery deadline is not met, subscribers may not extend their subscription. However, if the newspaper is delivered before the deadline, no increase in subscriptions is anticipated, unless an earlier deadline can be promised to all of the potential
252
M.G. Van Buer et al. / European Journal of Operational Research 115 (1999) 237±253
subscribers in the market area. Therefore, the general objective in the newspaper industry is to provide the promised service in the least cost manner. The costs associated with a particular production/distribution schedule are readily determined. There are no signi®cant out-of-pocket costs associated with a production set-up, only opportunity costs. One of the larger costs in distributing the newspaper is the cost of drivers and vehicles, which can be computed with good precision. There is also a fuel cost associated with the length of the truck routings. In order to mix these two costs gracefully in the objective function, we use a factor V that converts the net present value of a truck and driver to equal the net present value of operating a truck for one time unit (assumed proportional to distance) every day. Since the delivery deadline is considered to be a hard deadline, we do not need to consider the costs associated with late delivery of a newspaper. The problem can be formulated as follows: min
N X K K X N X N X X V x0jk Cij xijk j1 k1
s.t. lp
H
1 R
N X N X
P X
i1 j1
p1
P X P X
A:1
k1 i0 j0
Djp xijp
H H 1; . . . ; K;
p1 q1
A:2 tp
H
H X
lp
k
k1
N X N X
Cij xijp
H ;
H 1; . . . ; K;
i0 j1
A:3 tk 6 T ;
K 1; . . . ; K;
A:4
N X N N X N X X Zjp xijk ÿ ypk xijk 0; i0 j0
p 1; . . . ; P ;
i0 j0
k 1; . . . ; K;
i0
A:5
j 0; . . . ; N ;
p1
k 1; . . . ; K; N X K X
xijk 1;
A:6 i 1; . . . ; N ;
A:7
j0 k1 N N X X xijk ÿ xjik 0; i1
j 0; . . . ; N ;
i0
k 1; . . . ; K; XX
xijk 6 jSj
A:8 8S f1; . . . ; N g;
i2s j2s
k 1; . . . ; K; ypk 2 f0; 1;
p 1; . . . ; P ;
xijk 2 f0; 1; i 0; . . . ; N ; k 1; . . . ; K;
A:9 k 1; . . . ; K;
A:10 j 0; . . . ; N ;
A:11
p a vector containing a permutation of the integers from 1 to K:
!
Spq ypp
H ÿ1 yqp
H ;
! N P X X Djp xijk 6 Q;
A:12
The ®rst term of the objective function (A.1) accounts for the cost of acquiring a vehicle and driver. The second term describes the variable cost of operating a vehicle. Constraints (A.2) de®ne the time it takes to produce the newspapers for a given route. The ®rst term is the production time for the route and the second term accounts for any set-up time that may have been incurred. Constraints (A.3) de®ne the time that the last newspaper bundle is delivered on route k. The ®rst term is the departure time of vehicle k, the sum of the lk variables, and the second term is the actual route time. Notice that the route time does not include the return stem since the time that the vehicle returns to the production facility does not aect the feasible delivery of the newspapers. The Cij take into account the service time at drop-o j. We assume that the service time at each drop-o point is constant. Constraints (A.4) are simple bounds on the latest delivery time. This constraint is somewhat analogous to time windows. If the last
M.G. Van Buer et al. / European Journal of Operational Research 115 (1999) 237±253
drop-o point receives its newspapers on time, all preceding drop-os were necessarily serviced on time. Constraints (A.5) simultaneously de®ne the ypk variables and ensure that only one product is delivered on any route. Consider a given k (route) and p (edition). The ®rst sum will be non-zero if any Zjp is non-zero for which a corresponding xijk is non-zero and will be equal to the number of drop-o points on the route, k, if all receive the same edition, p. Considered in isolation, the second sum will always be equal to the number of drop-o points on the route or else zero. Hence, by requiring their equality, the constraints ensure that only one edition will be present on each route and the ypk value for that edition on that route will be one (and otherwise zero). Constraints (A.6) are vehicle capacity constraints. Constraints (A.7) ensure that every drop-o point is on a route. Constraints (A.8) ensure that vehicles do not accumulate at a drop-o point. Constraints (A.9) are the standard subtour elimination constraints. Constraints (A.10), (A.11), and (A.12) are integrality constraints.
References [1] E.H.L. Aarts, J. Korst, Simulated Annealing and Boltzmann Machines, Wiley, New York, 1989. [2] R. Battiti, G. Tecchiolli, The reactive tabu search, ORSA Journal on Computing 6 (1994) 126±140. [3] D. Blumfeld, C. Daganzo, L. Burns, Synchronizing production and transportation schedules, Transportation Research-B 25 (1991) 23±37. [4] P. Downey, J. Bruno, Complexity of task sequencing with deadlines, set-up times and changeover costs, SIAM Journal of Computing 4 (1978) 393±404. [5] M. Daskin, Logistics: an overview of the state of the art and perspectives on future research, Transportation Research-A 19 (1985) 383±398. [6] D. Johnson, M. Garey, Computers and Intractability: A Guide to the Theory of NP-Completeness, Freeman, San Francisco, 1979.
253
[7] M. Gendreau, A. Hertz, G. Laporte, A Tabu Search Heuristic for the Vehicle Routing Problem. [8] F. Glover, Tabu search-part I, ORSA Journal on Computing 1 (3) (1989) 190±206. [9] F. Glover, Tabu search-part II, ORSA Journal on Computing 2 (1) (1990) 4±32. [10] F. Glover, M. Laguna, Tabu Search, Kluwer Academic Publishers, Boston, 1997. [11] B. Hajek, Cooling schedules for optimal annealing, Mathematics of Operations Research 13 (1988) 311± 329. [12] A. Watts, J. Holt, Vehicle routing and scheduling in the newspaper industry, In: B. Golden, A. Assad (Eds.), Vehicle Routing: Methods and Studies, North-Holland, New York, 1988. [13] A.P. Hurter, M. Van Buer, The newspaper production/ distribution: Medium sized newspapers, The Journal of Business Logistics 17 (1) (1996) 85±108. [14] D.S. Johnson, C.R. Aragon, L.A. McGeoch, C. Schevon, Optimization by simulated annealing: An experimental evaluation; Part I, Graph partitioning, Operations Research 37 (1989) 865±892. [15] M. Lundy, A. Mees, Convergence of an annealing algorithm, Mathematical Programming 34 (1986) 111± 124. [16] M. Malek, J. Guruswamy, H. Owens, M. Pandya, Serial and parallel search simulated annealing and tabu search algorithms for the travelling salesman problem, Annals of Operations Research 21 (1989) 59±84. [17] S.L. Savage, Some theoretical implications of local optimization, Mathematical Programming 10 (1976) 354±366. [18] M. Savelsbergh, Local search in routing problems with time window constraints, Annals of Operations Research 4 (1985) 285±306. [19] M. Van Buer, J. Pellissier, A. Hurter, M. Daskin, D. Chyou, Location, production and distribution on M-rays, Working paper, Management Science Department, Loyola University, Chicago, 1996. [20] M. Van Buer, The newspaper distribution problem: The integrated design of a two-stage distribution system, Ph.D. Dissertation, Northwestern University, Evanston, IL, 1992. [21] D.L. Woodru, E. Zemel, Hashing vectors for tabu search, Annals of Operations Research 41 (1993) 123± 137. [22] D.L. Woodru, Proposals for chunking and tabu search, European Journal of Operational Research 106 (1998) 585±598.