UFDC Home  myUFDC Home  Help 



Full Text  
PAGE 1 VERY LARGESCALE NEIGHBORHOOD SEARCH HEURISTICS FOR COMBINATORIAL OPTIMIZATION PROBLEMS By KRISHNA CHANDRA JHA A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLOR IDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2004 PAGE 2 This work is dedicated to my father PAGE 3 iii ACKNOWLEDGMENTS First of all, I would like to thank Dr. Ravindra K. Ahuja, who has been a great supervisor and mentor throughout my four y ears at the University of Florida. His inspiring enthusiasm and energy have been contagious, and his advice and constant support have been extremely helpful. I would also like to express my gratitude to Dr. Joseph P. Geunes, Dr. Elif Akcali, and Dr. Anand Rangrajan for serving as my supervisory committee members and giving me thoughtful advice. My special acknowledgement goes to my collaborators and colleagues at University of Florida and at other places, especially, Dr. Jian Liu, Dr. Claudio Cunha, Arvind Kumar, Guvenc Sahin, Michelle Hann a, Onur Seref and Kamalesh Somani. My special thanks go to Ameetkumar Motwani fo r the help in editing and proofreading of this dissertation. My utmost appreciation goes to my father, the source of all my strength, and to my family for their encouragement and eternal suppo rt. I would like to thank my wife Nili for her patience, kindness and continuous suppor t during my four years of study at University of Florida. Finally, my thanks go to my daughter Archita for the hardships faced by her during my study and to my son An ish, whose arrival pushed me to complete this dissertation at the earliest. PAGE 4 iv TABLE OF CONTENTS Page ACKNOWLEDGMENTS.................................................................................................iii LIST OF TABLES............................................................................................................vii LIST OF FIGURES.........................................................................................................viii ABSTRACT....................................................................................................................... ..x 1 INTRODUCTION........................................................................................................1 1.1 Dissertation Summary............................................................................................8 1.1.1 Quadratic Assignment Problem....................................................................9 1.1.2 WeaponTarget Assignment Problem........................................................10 1.1.3 University Course Timetabling Problem....................................................11 1.1.4 BlocktoTrain Assignment Problem.........................................................12 1.1.5 Train Design Scheduling Problem..............................................................13 1.2 Contribution Summary.........................................................................................14 1.2.1 Quadratic Assignment Problem..................................................................14 1.2.2 WeaponTarget Assignment Problem........................................................14 1.2.3 University Course Timetabling Problem....................................................15 1.2.4 BlocktoTrain Assignment Problem.........................................................15 1.2.5 Train Schedule Design problem.................................................................16 2 QUADRATIC ASSIGNMENT PROBLEM..............................................................18 2.1 Introduction...........................................................................................................18 2.2 Improvement Graph..............................................................................................25 2.3 Identifying Profitable MultiExchanges...............................................................30 2.4 Specific Implementations.....................................................................................34 2.5 The Neighborhood Search Algorithm..................................................................36 2.6 Accelerating The Search Algorithm.....................................................................38 2.7 Expressions for the Asymmetric QAP..................................................................39 2.8 Computational Testing..........................................................................................40 2.8.1 Accuracy of the Solution............................................................................42 2.8.2 Effect of Neighborhood Size......................................................................44 2.8.3 Effect of the Speedup Technique...............................................................46 2.9 Conclusions...........................................................................................................47 PAGE 5 v 3 WEAPONTARGET ASSIGNMENT PROBLEM...................................................49 3.1 Introduction..........................................................................................................49 3.2 Lowerbounding Schemes....................................................................................52 3.2.1 A Lower Bounding Scheme using an Integer Generalized Network Flow Formulation.......................................................................................52 3.2.2 A Minimum Cost Flow Based Lower Bounding Scheme..........................57 3.2.3 Maximum Marginal Return Based Lower Bounding Method...................61 3.3 Branch and Bound Algorithm...............................................................................64 3.4 A Very LargeScale Nei ghborhood Search Algorithm........................................65 3.4.1 A Minimum Cost Flow formulati on based Construction Heuristic...........65 3.4.2 The VLSN Neighborhood Structure...........................................................67 3.5 Computational Results..........................................................................................70 3.5.1 Comparison of the Lower Bounding Schemes...........................................70 3.5.2 Comparison of Branch and Bound Algorithms..........................................71 3.5.3 Performance of the VLSN Search Algorithm............................................72 3.6 Conclusions...........................................................................................................73 4 UNIVERSITY COURSE TIMETABLING PROBLEM...........................................75 4.1 Introduction...........................................................................................................75 4.2 Problem Definition...............................................................................................77 4.3 Literature Review.................................................................................................79 4.4 Mathematical Formulation....................................................................................80 4.5 Random Instance Generation................................................................................83 4.6 Building Initial Solution.......................................................................................85 4.6.1 Greedy Heuristic 1......................................................................................85 4.6.2 Greedy Heuristic 2......................................................................................86 4.7 Very Large Scale Neighborhood Search (VLSN)................................................87 4.8 VLSN with Tabu Search.......................................................................................90 4.9 Implementation.....................................................................................................91 4.10 Conclusion and Further Scope............................................................................93 5 BLOCKTOTRAIN AS SIGNMENT PROBLEM....................................................94 5.1 Introduction...........................................................................................................94 5.2 Mathematical Programming Formulations...........................................................99 5.2.1 The SpaceTime Network...........................................................................99 5.2.3 Pathbased IP Formulation.......................................................................106 5.3 Path Enumeration Algorithms............................................................................109 5.3.1 Path Enumeration Restrict ed by Number of Trains.................................110 5.3.2 Path Enumeration Algorithm Rest ricted by Number of Trains and Number of Labels..............................................................................113 PAGE 6 vi 5.4 Solution Approaches...........................................................................................115 5.4.1 Lagrangian Relaxation Based Heuristic...................................................116 5.4.2 The Greedy Heuristic Algorithms............................................................119 5.5 Computational Results........................................................................................121 5.6 Summary and Conclusions.................................................................................124 6 TRAIN SCHEDULE DESIGN PROBLEM.............................................................126 6.1 Introduction.........................................................................................................126 6.2 Problem Description...........................................................................................129 6.2.1 Input..........................................................................................................130 6.2.2 Objective Function...................................................................................130 6.2.3 Decision Variables....................................................................................132 6.2.4 Constraints................................................................................................132 6.3 Problem Decomposition.....................................................................................134 6.3.1 Phase I......................................................................................................134 6.3.2 Phase II.....................................................................................................135 6.3.3 Phase III....................................................................................................136 6.4 Phase I: Designing Train Segments....................................................................138 6.4.1 Constructing an Initia l Feasible Solution.................................................144 6.4.2 Finding an Improved Neighbor................................................................145 6.5 Phase II: Forming Trains....................................................................................146 6.6 Computational Results........................................................................................152 6.7 Conclusion and Further Scope............................................................................153 7 FURTHER RESEARCH AND CONCLUDING REMARKS.................................154 7.1 Contribution Summary.......................................................................................154 7.2 Future Scope.......................................................................................................156 APPENDIX: COMPUTATIONAL EXPE RIMENT ON QAP INSTANCES................158 LIST OF REFERENCES.................................................................................................163 BIOGRAPHICAL SKETCH...........................................................................................172 PAGE 7 vii LIST OF TABLES Table page 2.1 Number of iterations with di fferent implemented cycle length...............................46 2.2 Effect of accelerated path enumeration scheme.......................................................47 3.1 Comparison of four lower bounding schemes..........................................................71 3.2 Comparison of branch and bound algorithms..........................................................72 3.3 Results of the construction heuris tic and the VLSN search algorithm.....................73 4.1 Performance of VLSN and Tabu search..................................................................92 5.1 Performance analysis of different solution approaches..........................................122 6.1 Illustrating the cost of fl ow in the blocking network.............................................141 A.1 Experimental Results for Symmetric Instances.....................................................161 A.2 Experimental Results for Asymmetric Instances...................................................162 PAGE 8 viii LIST OF FIGURES Figure page 2.1 Facilities Locations: (a ) Initial assignment and (b) Assignment after the cyclic exchange 3763...................................................21 2.2 The generic search procedure for identifying profitable multiexchanges...............32 2.3 The neighborhood search algorithm for the QAP....................................................37 2.4 Effect of cycle length on time taken and solution quality for 100 runs on problem sko100af....................................................................................................45 2.5 Effect of number of paths in each stage on time taken and solution quality for 100 runs on problem sko100af..........................................................................46 3.1 Formulating the WTA problem as Integer generalized network flow problem.......54 3.2 Approximating a convex function by a lo wer envelope of linear segments............56 3.3 A network flow formulation of the WTA problem..................................................58 3.4 Combinatorial lower bounding algorithm................................................................63 3.5 The VLSN search algorithm for the WTA problem................................................69 4.1 Algorithm for assignment of c ourses using Greedy Heuristic 2..............................87 5.1 A part of the spacetime network...........................................................................102 5.2 Path enumeration algorithm restricted by the number of trains.............................111 5.3 Example of path enumeration algorithm................................................................111 5.4 Path enumeration algorithm restrict ed by number of trains and number of labels..................................................................................................................114 5.5 Lagrangian relaxation ba sed lower bounding scheme...........................................119 5.6 Greedy Heuristic for the BTA problem..................................................................120 5.7 Analysis of performance with respect to maximum number of labels at a node...123 PAGE 9 ix 6.1 An example of a train segment network.................................................................139 6.2 The VLSN search algorithm for the TSD problem................................................144 6.3 Graph for minimum cost flow formulation............................................................149 PAGE 10 x Abstract of Dissertation Pres ented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy VERY LARGESCALE NEIGHBORHOOD SEARCH HEURISTICS FOR COMBINATORIAL OPTIMIZATION PROBLEMS By Krishna Chandra Jha May 2004 Chair: Ravindra K. Ahuja Major Department: Industrial and Systems Engineering Combinatorial optimization plays an impor tant role in decisionmaking since optimal decisions often depend on a nontrivial combination of various factors. Most combinatorial optimization problems are NP hard and sharp bounds on the optimal value are typically hard to derive. This means that partial enumeration based exact algorithms have a slow convergence rate, and they can solve only small instances to optimality. But reallife combinatorial optimization problems are typically of big size and since exact approaches are inadequate, heuristics are co mmonly used in practice. Over the last 15 years much of the research effort has been concentrated on the development of local search heuristics. In local search methods, an intensive exploration of the solution space is performed by moving at each step from the current solution to another promising solution in its neighborhood. As a rule of thum b, there is a better probability of finding a promising solution, if the neighborhood size is larger. In this dissertation, we have developed algorithms, which perform intens ive searches in a very large size PAGE 11 xi neighborhood. We have used the concepts of network flow algorithms to make our algorithms computationally efficient. In Chapter 2, we propose very larges cale neighborhood (VLSN) search algorithms for the quadratic assignment problems (QAP), and we propose a nov el search procedure to enumerate good neighbors heuris tically. Our search procedur e relies on the concept of an improvement graph, which allows us to evaluate neighbors much faster than the existing methods. We present extensive com putational results of our algorithms on a standard benchmark QAP instances. In Chapter 3, we suggest linear programming, integer programming, and network flow ba sed lower bounding; using these methods, we obtain several branch and bound algorithms fo r the weapontarget assignment problem. We also propose a network flow based cons truction heuristic and a very largescale neighborhood (VLSN) search al gorithm. In Chapter 4, we have developed VLSN search heuristics in conjunction with tabu search to solve a uni versity coursetimetabling problem. We have developed a generic m odel for the problem and have performed computational experiments on randomly generated test instances. In Chapter 5, we have given integer pr ogramming formulations of the blocktotrain assignment problem; propose several exac t and heuristic algorith ms to solve this problem, and present computati onal results on the data provided by a major US railroad. In chapter 6, we have proposed a decompositi on based approach to solve another railroad scheduling problem, called the train schedule design problem. We have developed VLSN search heuristics to solve the subproblem in first phase and have formulated the problem in second phase as a minimum cost flow problem. PAGE 12 1 CHAPTER 1 INTRODUCTION Optimization has been expanding in all direc tions at an astonishing rate during the last few decades. New algorithmic and theo retical techniques have been developed, the diffusion into other disciplines has proceeded at a rapid pace, and our knowledge of all aspects of the field has grown even more profound. At the same time, one of the most striking trends in optimization is the constantly increasing emphasis on the interdisciplinary nature of the field. Optimization today is a basic research tool in all areas of engineering, medicine and the sc iences. The decision making tools based on optimization procedures are successfully appl ied in a wide range of practical problems arising in virtually any s phere of human activities, including biomedicine, energy management, aerospace research, telecommunications and finance. The problems of finding the best and th e worst have always been of great interest. For example, given n sites, what is the fastest way to visit all of them consecutively? In manufacturi ng, how should one cut plates of a material so that the waste is minimized? Some of the first optim ization problems were solved in ancient Greece and are regarded among the most significan t discoveries of that time. In the first century A.D., Alexandrian mathematician Heron solved the problem of finding the shortest path between two points by way of the mirror. This result, also known as Herons theorem of the light ray, can be viewed as the origin of the theory of geometrical optics. The problem of finding extreme valu es gained a special importance in the PAGE 13 2 seventeenth century when it served as one of motivations in the invention of differential calculus. The soonafterdeveloped calculus of variations and the th eory of stationary values lie in the foundation of the modern mathematical theory of optimization. Many wellstructured problems encountered in business (and other domains) can be defined as optimization problems P : maximize (minimize) f ( x ) s.t. x S where S represents the set of admissible solutions (satisfying all hard constraints) of an underlying search space X The function where R represents the set of real numbers, is often called goal function or objec tive function. A large subclass of these problems is the class of combinatorial optimization problems. An optimization problem is called combinatorial, if X and S are combinatorial or discrete in some sense. Although ther e does not seem to be the formal common understanding of what is 'combinato rial', we consider here that X and/or S are 'combinatorial' if they are discrete sets of fi nite elements or countably infinite elements (Ibaraki [1987]). Combinatorial optimization plays an impor tant role in decisionmaking since optimal decisions often depend on a nontrivial combination of various factors. Most combinatorial optimization problems are NPhard, and sharp bounds on the optimal value are typically hard to derive. This means that partial enumeration based exact algorithms have a slow convergence rate, and they can solve only small instances to optimality. But reallife combinatorial optimization problems are typically of big size, and since exact approaches are inadequate, heuristics are co mmonly used in practice. Etymologically, the PAGE 14 3 word heuristic comes from the Greek heuriskein (to find). There has been a steady evolution over the past forty years in th e development of heuristics, which produce solutions of reasonably good quality in a reason able amount of time. The first proposed heuristics tried to systemize decisionmaking processes done by hand. With the help of computers, which can test a huge amount of combinations in a short amount of time, solutions could easily be generated which have turned out to be of much better quality than what an expert in the field could produ ce by hand. In these early heuristics, much of the emphasis was put on quickly obtaining a fe asible solution and possibly applying to it a postoptimization procedure. The literature devoted to heuristic al gorithms often distinguishes between two broad classes: constructive algorithms a nd improvement algorithms. A constructive algorithm builds a solution from scratch by as signing values to one or more decision variables at a time. An improve ment algorithm generally starts with a feasible solution and iteratively tries to obtain a better solution. Neig hborhood search algorithms (alternatively called local search algorithms) are a wide class of improvement algorithms where at each iteration an improving solu tion is found by search ing the neighborhood of the current solution. A procedure of this t ype usually terminates when the first local optimum is obtained. Randomization and restar ting approaches used to overcome poor quality local solutions are often ineffective. More general strategies known as metaheu ristics usually combine some heuristic approaches and direct them towards solutions of better quality th an those found by local search heuristics. Over the last 15 years much of the research effort has concentrated on the development of metaheuristics, using mainly two principles: local search and PAGE 15 4 population search. The term metaheuristic was coined by Glover[1986] and has come to be widely applied in the literature, both in the titles of comparative studies and in the titles of collected research pape rs. A metaheuristic refers to a master strategy that guides and modifies other heuristics to produce solutions beyond those that are normally generated in a quest for local optimality. The heuristics guided by su ch a strategy may be a high level procedure or ma y embody nothing more than a description of available moves for transforming one solution to anothe r, together with an associated evaluation rule. The evolution of metaheuristics during th e past decade has taken an explosive upturn. Metaheuristics in their modern forms are based on a variety of interpretations of what constitutes intelligent search. These inte rpretations lead to design choices that in turn can be used for classification purpos es. However, a rigorous classification of different metaheuristics is a difficult and ri sky enterprise, because the leading advocates of alternative methods often differ among them selves about the essential nature of the methods they espouse. Metaheuristics are ofte n viewed as composed of processes that are intelligent, but in some instances the intelligence belongs more to the underlying design than to the particular characte r (or behavior) of the method itself. In local search methods, an intensive exploration of th e solution space is performed by moving at each step from the current so lution to another promising solution in its neighborhood. Simulated annealing (Kirkpatrick et al. [1983]), Tabu search (Glover [1986]) and variable neighborhood search (Miladenovic and Hanson [1997]) are the most famous local search methods. Population sear ch consists of maintaining a pool of good solutions and combining them in order to produce hopefully better solutions. Classical PAGE 16 5 examples are genetic algorithms (Holland [1975]) and adaptive memory procedures (Rochat and Taillard [1995]). Metaheuristics are genera l combinatorial optimization techniques, which are not dedicated to the solu tion of a particular pr oblem, but are rather designed with the aim of being flexible e nough to handle as many different combinatorial problems as possible. These general tech niques have rapidly demonstrated their usefulness and efficiency in solving hard problems. Success storie s are reported in many papers. While metaheuristics can handle in theory any combinatorial optimization problem, it is often the case th at an important effort must be put on finding the right way to adapt the general ingredient s of these methods to the pa rticular considered problem. We think that in order to be successful in the adaptation of a metaheuristic to a combinatorial optimization problem, it is necessa ry to follow some basic principles. We give in the next sect ions some guidelines, which may he lp in producing su ch successful adaptations of local search and population search methods for the solution of difficult combinatorial optimization problems. This dissertation concentrates on neighborhood search algorithms where the size of the neighborhood is very large with respect to the size of the input data. For large problem instances, it is impractical to s earch these neighborhoods explicitly, and one must either search a small portion of th e neighborhood or else develop efficient algorithms for searching the neighborhood implicitly. A critical issue in the design of a neighbor hood search approach is the choice of the neighborhood structure, that is, the manner in which the neighborhood is defined. This choice largely determines whether the nei ghborhood search will develop solutions that are highly accurate or whether they will deve lop solutions with very poor local optima. PAGE 17 6 As a rule of thumb, the larger the neighborhood the better is the qua lity of the locally optimal solutions, and the greater is the accuracy of the final solution that is obtained. At the same time, the larger the neighborhood, th e longer it takes to search the neighborhood at each iteration. Since one generally pe rforms many runs of a neighborhood search algorithm with different starti ng points, longer execution times per iteration lead to fewer runs per unit time. For this reason a larg er neighborhood does not ne cessarily produce a more effective heuristic unless one can sear ch the larger neighbor hood in an efficient manner. Some very successful and widely used methods in operations research can be viewed as very largescale neighborhood search techniques. For example, if the simplex algorithm for solving linear programs is view ed as a neighborhood search algorithm, then column generation is a very largesc ale neighborhood search method. Also, the augmentation techniques used for solving many network Bows problems can be categorized as very largescale neighborhood search methods The negative cost cycle canceling algorithm for solving the min cost bow problem and the augmenting path algorithm for solving matching problems are two such examples. The combinatorial optimization problems, we studied in this di ssertation, are very largescale optimization problems and are im possible to solve to optimality using the stateoftheart algorithmic approaches. In this dissertation, we focus on developing heuristic algorithms for this problem that do not guarantee optimality of the solution but find a nearly optimal solution within a reas onable computational time. The literature devoted to heuristic algorithms often di stinguishes between two broad classes: constructive algorithms and neighborhood search algorithms (alternatively called local PAGE 18 7 search algorithms ). A constructive algorithm builds a solution of an optimization problem from scratch by assigning values to one or more decision variables at a time. In contrast, a neighborhood search algorith m starts with a feasible solution x of the problem, defines a set of neighboring solutions ( x ) of the solution x called the neighborhood of x evaluates neighboring solutions and if it finds a solution y ( x ) that is better than x then it replaces x with y It next defines the neighborhood with respect to the new solution x and repeats this pro cedure until the solution x is at least as good as any other solution in ( x ). The solution x is then declared as a local optimal solution and the algorithm terminates. The use of neighborhood search algorithms in solving optimization problems has a long history (Croes [1958], Lin [1965], Lin and Kernighan [1973], Kirkpatrick et al. [1983], Reeves [1993], Osman and Laporte [1996], and Aarts and Lenstra [1997]). However, in the past fifteen years, the fi elds of operations research, mathematical programming and computer science have all witnessed a strong renewed interest in the development and analysis of neighborhood search based approaches. This interest can be attributed in large part to the intuitive appe al of these algorithms, flexibility, and the ease of implementations of these algorithms to solve many complex realworld problems. Neighborhood search algorithms are now widely regarded as an important tool to solve difficult optimization problems effectively, and in particular, combinatorial optimization problems. A critical issue in the design of a neighbor hood search algorithm is the choice of the neighborhood structure ; that is, the manner in which the neighborhood is defined. In general, the larger the neighborhood, the better is the quality of loca lly optimal solutions PAGE 19 8 and greater is the accuracy of the final so lution. At the same time, the larger the neighborhood, the longer it takes to enumerate the neighborhood and evaluate solutions in it. Hence, a larger neighborhood does not necessarily produce a more effective heuristic until one can efficientl y search the larger neighborhood. Recently, numerous efforts are made in developing very largescale neighborhood (VLSN) search algorithms wher e the size of the neighborhood is very large. In fact, it is so large that enumerating all neighbors and evaluating them is prohibitively expensive. Using concepts from network flow theory (A huja et al. [1993]), a VLSN search algorithm implicitly enumerates the neighborhood to identify an improved neighbor. A VLSN search algorithm can consider neighborhoods with trillions of neighbors while the time to identify improved neighbor is fairly small, often less than a fraction of a second (Ahuja et al. [2002a). Some of the inte resting applications of VLSN algorithms are for airline fleet scheduling problems (Ahuja et al. [ 2001b, 2003b, 2003c]), the locomotive scheduling problem (Ahuja et al. [2002b]), the capacita ted minimum spanning tree problem (Ahuja et al. [2001a, 2003a]), the capac itated facility location probl em (Ahuja et al. [2002c]). 1.1 Dissertation Summary As the title suggests, in this dissertation efforts have been made to solve a few hard combinatorial optimization problems by deve loping VLSN search heuristics. The problems studied are the quadratic assignmen t problem, the weapontarget assignment problem, the university coursetimetabling problem, the blocktotrain assignment problem, and the train design pr oblem. We give below a brie f summary of our efforts to solve these problems. PAGE 20 9 1.1.1 Quadratic Assignment Problem The Quadratic Assignment Proble m (QAP) consists of assigning n facilities to n locations so as to minimize the total weighted cost of interactions between facilities. In this dissertation, the QAP arises in many di verse settings, is known to be NPhard, and can be solved to optimality only for fairly small size instances (typically, n 30). Neighborhood search algorithms are the most po pular heuristic algorit hms to solve larger size instances of the QAP. The most exte nsively used neighbor hood structure for the QAP is the 2exchange neighborhood. This neighborhood is obtained by swapping the locations of two facilities and thus has size O( n2). Previous efforts to explore larger size neighborhoods (such as 3exchange or 4exchange neighborhoods) were not very successful, as it took too long to evaluate the larger set of neighbors. In this dissertation, we propose very largescale neighborhood (VLSN) search al gorithms where the size of the neighborhood is very larg e and we propose a novel search procedure to enumerate good neighbors heuristically. Our search pr ocedure relies on the concept of an improvement graph, which allows us to evalua te neighbors much faster than the existing methods. We present extensive computational re sults of our algorithms on standard benchmark instances available in QAPLIB (http://www.opt.math.tugraz.ac.at/qaplib). There are two types of quadrat ic assignment problems: (i) fl ow and distance matrices are symmetric (for example, distance from location A to B is same as that from B to A), and (ii) flow and distance matrices are asymmetric (for example, distance from location A to location B is different than that from B to A). In this dissertati on, we have tested our algorithms on both types of problems. The computational results show that our PAGE 21 10 algorithms are very robust and produce very good quality solutions for all kinds of problem instances. 1.1.2 WeaponTarget Assignment Problem The WeaponTarget Assignment (WTA) probl em is a fundamental problem arising in defenserelated applications of operations research. This problem consists of optimally assigning n weapons to m targets so that the total expect ed survival value of the targets after all the engagements is minimum. Th e WTA problem can be formulated as a nonlinear integer programming problem a nd is known to be NPcomplete. The researchers have studied this problem for more than thr ee decades; however, there do not exist any exact methods for the WTA problem, which can solve even small size problems (for example, with 20 wea pons and 20 targets). Though several heuristic methods have been proposed to solve the WTA problem, due to the absence of exact methods, no estimates are available on the quality of solu tions produced by such heuristics. In this dissertation, we suggest linear programm ing, integer programming, and network flow based lower bounding methods using whic h we obtain several branch and bound algorithms for the WTA problem. We have used piecewise approximation to obtain linear objective function for the problem. We have developed algorithms to determine the bounds on number of weapons assigned to a targ et. These bounds play a critical role in formulating the problem as a minimum cost flow problem. We also propose a network flow based construction heuristic and a VLSN search algorithm. The combination of network flow based construction heuristic and VLSN search algorithms is very efficient and produces ex ceptionally good quality solutions for the problem instances even of si ze 200 weapons and 200 targets. We present computational results of our algorithms, which indicate that we can solve moderately large size PAGE 22 11 instances (up to 80 weapons and 80 targets) of the WTA problem optimally and obtain almost optimal solutions of fairly large instances (up to 200 weapons and 200 targets) within a few seconds. 1.1.3 University Course Timetabling Problem The university C ourse Timetabling Problem (CTP) is to develop a weekly course timetable containing the time and the meeti ng place for each course over a given time horizon (generally, a semester or trimester). The university coursetimetabling problem is much different than the schooltimetabling pr oblem. In the school timetabling problem, students are generally required to take a fixed set of courses, whereas students have lot of freedom in selecting courses in a university Also, teachers are mainly engaged in teaching courses at a school, whereas they ha ve dual responsibilities of teaching and research in a university. The objectives of th e university coursetimetabling problem are to maximize the preferences of teachers for the periods at which they want to teach, and to minimize the time clash among the courses a student wants to take. Generally, teachers give their preferences at the beginning of the planning period, whereas the students interest data are collected from the past enrollment. The course timetable problem is a highly constrained optimization problem, in wh ich a constraints set widely varies from one case to another. Previous efforts to solv e the problem are mainly case specific and are limited to small instances. In this dissert ation, we propose a generic model of the problem, which can be applied to a wide variety of reallife instances. In our model, the periods are grouped into slots and we assign the courses to the slots to optimize the objective function. We have formulated the coursetimetabling problem as an integerprogramming problem. We have developed a co nstruction heuristic and have developed a VLSN search algorithm to solve the problem We have used a popular metaheuristic, PAGE 23 12 tabu search, in conjunction with VLSN search algorithms. This approach permits us to search for the better solution in wide solution space. We report the performance of our algorithm when compared to an optimal solu tion on randomly generated test instances. 1.1.4 BlocktoTrain Assignment Problem Railroads classify shipments into blocks to reduce intermediate handlings of cars. Once a shipment is assigned to a block, it is not reclassified at the intermediate yards on the route from the origin of the bloc k to the destination of the block. The blocktotrain assignment (BTA) problem consists of assigning th ese blocks to trains so that the shipments reach their destinations such that the transportation cost and transit time are minimum. The transportation cost of a block consists of distance traveled by the block and cost of swapping whenever it changes the train on the route. The transit time for a block is the time taken in traveling from its origin to its destination. The solution of the BTA problem must honor the cap acity constraints on the trains. This is a large size hard combinatorial optimization problem. In the past very few efforts have been made to solve the BTA problem. In this dissertation, we formulate this problem as an integer programming problem on a spacetime network. We have developed label enumerating algorithms to enumerate possible routes for a block in the spacetime network. Our label correcting algorithm is a VLSN s earch technique, which very efficiently produces a set of best routes out of enormous possible routes for a block. We also propose several exact and heuristic algorithms to solve this problem We have tested our algorithms on a real life BTA problem of a major US railroad company. PAGE 24 13 1.1.5 Train Design Scheduling Problem The train design scheduling (TDS) problem is one of the most difficult and complex scheduling problems in the railroad. Due to the ever changing nature of the transportation demand, railroads are required to come up with a completely new train schedule or will have to change the existing train schedule. The train schedule design problem is to design a zerobased (do not c onsider existing solution) train schedule, in which decisions are made with respect to trains to be made, their origins and destinations, their routes and frequency, the time schedule, and the blocktotrain assignment. The objectives of the TDS problem are to minimize the overall transportation cost and transit time of shipments. There are many logical and practical constraints, which a solution of the TDS problem is required to honor. For a r eallife instance of any major US railroad, the TDS problem consists of millions of integr al decision variables and constraints. With the current age computational capability, it is not possible to solv e a problem of this magnitude using any of the stateofart optim ization techniques. In the past, few efforts were made to solve this problem heuristically, however, they were implemented on a much simplified version of the problem or to small instances. In this dissertation, we have developed a decomposition based approach to solve the problem. Our approach consists of solving the problem sequentially in thr ee phases. We have developed algorithms to solve the subproblems in Phase I and Phase II. The problem in Phase I is a network design problem and we have developed a VLSN search algorithm to solve this problem, which sequentially solve the problem in this phase. We have formulated the problem in Phase II as a minimum cost flow problem on a suitably constructed graph, in which we define the cost of arcs in minimum cost flow formulation in accordance with the objectives of the problem in this phase. The computational expe riment on a real life PAGE 25 14 instance indicates that the problems in Phase I and II can be solved very efficiently using a normal computational facility. 1.2 Contribution Summary We summarize below the contributio ns made in this dissertation. 1.2.1 Quadratic Assignment Problem We have developed a VLSN search al gorithm for the qua dratic assignment problem, a nonpartitioning optimization pr oblem. The application of VLSN search algorithm is different than that of partitioning problems, to which this algorithm was applied earlier. We have an approximated an improvement gr aph, in which the cost of the cycle is a very good approximation of the cost of the multiexchange, and allows us to enumerate good neighbors quick ly. Typically, evaluating a k exchange neighbor for the QAP takes O( nk ) time; but using the improvement graph we can do it in O( k ) average time per neighbor. We have developed a generic search procedure to enumerate neighbors using improvement graphs. We have also deve loped several implementations of the generic search procedure, which enumer ate the neighborhoods exactly as well as heuristically. We present a detailed computational inve stigation of local improvement algorithms based on our neighborhood sear ch structures. In our e xperimentation, we have considered 98 symmetric instances and 34 as ymmetric test instances available in the QAPLIB (QAP Library). 1.2.2 WeaponTarget Assignment Problem We formulate the WTA problem as an inte ger linear programming problem, that is, as a generalized integer network flow problem on an appropriately defined network. The linear programming relaxati on of this formulation gives a lower bound on the optimal solution of the WTA problem. We propose a minimum cost flow formula tion that yields a different lower bound on the optimal solution of the WTA problem This lower bound is, in general, not as tight as the bound obtained by the lin ear programming form ulation described above but it can be obtained in much less computational time. PAGE 26 15 We propose a third lower bounding scheme, which is based on simple combinatorial arguments and uses a greedy approach to obtain a lower bound. We develop branch and bound algorithms to solve the WTA problem employing each of the three bounds described above. We propose a very largescale neighborhood (VLSN) search algorithm to solve the WTA problem. The VLSN search algorithm is based on formulating the WTA problem as a partition problem The VLSN search starts with a feasible solution of the WTA problem and performs a sequence of cyclic and path exchanges to improve the solution. We describe a heur istic method that obtains an excellent feasible solution of the WT A problem by solving a sequence of minimum cost flow problems, and then uses a VLSN search algorithm to iteratively improve this solution. We perform extensive computational inve stigations of our algorithms and report these results. Our algorithms solve moderately large size instances (up to 80 weapons and 80 targets) of the WTA probl em optimally and obtain almost optimal solutions of fairly large instances (up to 200 weapons and 200 targets) within a few seconds. 1.2.3 University Course Timetabling Problem We have developed a generic model of the university course timetabling problem. We expect that this model can be appl ied to many reallife course timetabling problems. We formulate the course timetabling probl em as an integer programming problem. We have been able to solve moderate size problems to optimality using CPLEX optimizer. The course timetabling problem is a partitio ning problem with side constraints. We develop an improvement graph based VL SN search algorithm to solve this problem. Due to the constraints in course timetabling problem, solution space is very sparse. We developed an algorithm in conjunction with tabu search, which permits the algorithm to search for the good quality solution in wide solution space. We generated the test instances randomly and performed computational experiment of the algorithms. Our heuristic algorith ms perform extremely well on the test instances and ware able to produce optimal solution in almost all the cases. 1.2.4 BlocktoTrain Assignment Problem We develop a spacetime network, which al lows us to formulate the BTA problem as a multicommodity flow problem with additional side constraints. PAGE 27 16 We develop two integer programming form ulations of the BTA problem. The first formulation formulates the BTA probl em as a nodearc version of multicommodity flow problem in the spacetime network. The second formulation conceives the BTA problem as a pathflo w version of the multicommodity flow problem in the spacetime network. The integer programming formulations can also be used to solve the BTA problem approximately. Using these formulations we suggest a Lagrangian relaxation algorithm that can solve the BTA problem to nearoptimality within a few seconds. We also propose a greedy construction al gorithm that proceeds by assigning blocks to trains one by one until all blocks are assigning. By choosing different rules for selecting blocks and assigning them to trains, we can obtain different implementations of this generic approach. We perform extensive com putational investigations of these algorithms. Our integer programming formulations can solv e the BTA problem to optimality within a few minutes of computer time using CPLEX 8.1. The Lagrangian relaxation algorithm can solve the BTA problem to nearoptimality within a few seconds. The computational performance of greedy c onstruction algorithms is also very attractive. 1.2.5 Train Schedule Design problem We have developed a decomposition based ap proach to solve reallife train design problem. The sub problem in phase I is a network design problem and we have developed a VLSN search heuristic to solve the problem in this phase. We have formulated the problem in Phas e II as a minimum cost flow problem on a suitably constructed graph, which can be solved very effici ently using network flow optimization technique. We perform the computational experime nt on the algorithms using reallife problem instance. Our results show that our approach can produce substantial saving in terms of numb er of trains required. To summarize, our contributions in this dissertation are in de veloping very large scale neighborhood search heuristics, which use efficient network flow algorithms to obtain good quality solutions of different combinatorial optimization problems. Our strength lies in proper formul ations of the problems and st ateofart implementation of PAGE 28 17 algorithms. We have focused on developing algorithms such that near optimal solution can be obtained very quickly using a normal com putational facility. We have also tried to develop robust algorithms, which can be a pplied across different reallife scenarios. PAGE 29 18 CHAPTER 2 QUADRATIC ASSIGNMENT PROBLEM 2.1 Introduction The Quadratic Assignment Problem (QAP) is a classical combinatorial optimization problem and is widely regarded as one of the most difficult problems in this class. Given a set N = {1, 2, ..., n}, and n x n matrices F = { fij}, D = { dij}, and B = { bij}, the QAP is to find a permutation of the set N which minimizes: z() = 11 nn ij fijd(i)(j) + 1 n i bi(i) 2.1 The QAP arises as a natural problem in facility layout. In this context, the set N represents a set of n facilities (numbered 1 through n ) that need to be assigned to locations (numbered 1 through n ). The matrix F = { fij} represents the flow between different facilities, and the matrix D = { dij} represents the distance between locations. For example, if the facilities are depa rtments in a campus, then the flow fij could be the average number of people walking daily from department i to department j The decision variable (i) 1 i n represents the location assigned to facility i Since there are n facilities and n locations and a facility can be assi gned to exactly one location, there is a onetoone correspondence between feasible solutions of QAP and permutations Observe that Equation 2.1 consists of two terms. The first term is the sum of n2 flow costs between n facilities (the term fijd(i)(j) represents the cost of flow from facility i to facility j ). The second term considers the cost of erecting facilitie s, which may be PAGE 30 19 locationdependent. The matrix B = { bij} represents the cost of creating facility i at location j Hence, the QAP is to find an assignment of facilities to locations so as to minimize the total cost of flow between th e facilities and the cost of erecting the facilities. The matrices F and D are typically symmetric matrices but are not required to be so. In our algorithms, we allow asymmetr ic instances and thus do not assume that fij = fji or dij = dji. However, for the sake of simplicity, we will assume in Sections 2.2 through 2.6 that we are working with symmetric instances. In Section 2.7, we will study asymmetric instances. In addition to the facility layout, the QAP arises in many other applications, such as the allocation of plants to candidate loca tions, the backboard wiring problem, design of control panels and typewriter ke yboards, turbine balancing, or dering of interrelated data on a magnetic tape, and others. The detail s and references for these and additional applications can be found in Malucelli [1993], Pardalos, Rendl and Wolkowicz [1994], Burkard et al. [1998], and Cela [1998]. Given the wide range of applications and the difficulty of solving the problem, the QAP has been investigated extensively by the research community. The QAP is known to be NPhard, and a variety of exact and heuristic algorithms have b een proposed. Exact algorithms for solv ing QAP include approaches based on (i) dynamic programming (C hristofides and Benavent [1989]); (ii) cutting planes (Bazar aa and Sherali [1980]); and (iii ) branch and bound (Lawler [1963], Pardalos and Crouse [1989]). Among these, the branch and bound algorithms are the most successful, but they are generally unable to solve problems of size larger than 25 facilites. PAGE 31 20 Since the applications of the QAP have often given rise to problems of size far greater than 25, there is a need for good heuristics for QAP that can solve larger size problems. A wide variety of heuristic appr oaches have been developed for the QAP. These can be classified into the following categories: (i) construction methods (Buffa, Armour and Vollmann [1964], MullerMerbach [1970]); (ii) limited enumeration methods (West [1983], Burkard an d Bonniger [1983]); (iii) GRASP (greedy randomized adaptive search procedure) (Li, Pa rdalos, and Resende [1994]); (iv) simulated annealing methods (Wilhelm and Ward [1987]); (v) tabu search methods (SkorinKapov [1990], Taillard [1991]); (vi) genetic algorithms (Fleurent and Ferland [1994], Tate and Smith [1985], Ahuja, Orlin, and Tewari [ 1998], Drezner [2003]); and (vii) ant systems (Maniezzo, Colorni, and Dorigo [1994]). Th e tabu search method of Taillard [1991], the GRASP method of Li, Pardal os, and Resende [1994], a nd the genetic algorithm by Drezner [2003] are the most accurate heuristics among these methods. As observed in the survey paper of Burk ard et al. [1998], the current neighborhood search metaheuristic (tabu search and simu lated annealing) algorithms for the QAP use the twoexchange neighborhood structure in th e search. A permutation is called a 2exchange neighbor of the permutation if it can be obtained from by switching the values of two entries in the permutation It is easy to see that the number of twoexchange neighbors of a permutation is O( n2). There has been very limited effort in the past to explore larger nei ghborhood structures for the QAP as the time needed to identify an improved neighbor becomes too high. In this chapter, we inves tigate the neighborhood structure based on multiexchanges which is a natural generalization of the 2exchanges. A multiexchange is specified by a cyclic sequence C = i1 i2 ik i1 of facilities PAGE 32 21 such that ip iq for p q This multiexchange implies that facility i1 is assigned to the location ( i2), facility i2 to ( i3), and so on, and finally facility ik is assigned to ( i1). The location of all other facilities is not changed. We denote by C, the permutation obtained by applying the multiexchange C to the permutation In other words, C( i ( i ) for i N \{ i1, , ik}, 2.2 C( ip) = ( ip +1) for p = 1 , , k 1 and C( ik) = ( i1). We define the length of a multiexchange as the number of facilities involved in the corresponding cyclic sequence. For example, the cyclic sequence C = i1 i2 ik i1 has length k We also refer to a multiexchange of length k as a kexchange Figure 2.1 illustrates an example of a 3exchange We note that a kexchange can be generated by k different cyclic sequences. For example, the 3exchange shown in Figure 2.1 can be generated by any of the sequences 3763, 7637 and 6376 Figure 2.1 Facilities Locations: (a) Initi al assignment and (b ) Assignment after the cyclic exchange 3763. Given a positive integer 2 K n the Kexchange neighborhood structure consists of all the neighbors of a permutation obtaine d by using multiexchanges of length at most 1 3 2 3 4 5 6 7 8 2 8 7 5 6 4 1 1 6 2 3 4 5 6 7 8 2 8 3 5 7 4 1 (a) (b) 1 3 2 3 4 5 6 7 8 2 8 7 5 6 4 1 1 6 2 3 4 5 6 7 8 2 8 3 5 7 4 1 (a) (b) PAGE 33 22 K We note that the twoexchange nei ghborhood structure is contained in the K exchange neighborhood structure. The num ber of neighbors in the K exchange neighborhood structure is ((1)! n K K ). This number is very large even for moderate values of K For example, if n = 100 and K = 10 then the Kexchange neighborhood may contain as many as 6x1018 neighbors. This neighborhood structure falls under the category of very largescale neighborhood (VLSN) structures where th e size of the neighborhood is too large to be searched explicitly and we use implic it enumeration methods to identify improved neighbors. Algorithms based on very largescal e neighborhood structures have been successfully used in the context of several combinatorial optimization problems (see Ahuja et al. [2002a], and Deineko and Woegi nger [1999] for surveys on this type). One of the tools used in performing search over very largescale neighborhood structures is the concept of the improvement graph In this technique, we a ssociate a graph, called the improvement graph G() with each feasible solution of the combinatorial optimization problem. The improvement graph G() is constructed such th at there is a onetoone correspondence between every neighbor of to some directed cycle (possibly satisfying certain constraints) in the improvement graph G() We also define arc costs in the improvement graph so that the difference in the objective function value of a neighboring solution and the solution is equal to the cost of the c onstrained cycle corresponding to the neighbor. This transforms the problem of finding an improved neighbor into the problem of finding a negative cost constraine d cycle in the improvement graph (assuming that the combinatorial optimization problem is a minimization problem). The concept of PAGE 34 23 improvement graph was first proposed by Th ompson and Orlin [1989] for a partitioning problem, where a set of elements is partitioned into several subsets of elements so as to minimize the sum of the objective functions of the subsets. This tec hnique has been used to develop several VLSN search algorithms for specific partitioning problems such as the vehicle routing problem (Thompson and Psar aftis [1993]) and the capacitated minimum spanning tree problem (Ahuja, Orlin, Sh arma [2003a, 2001a]). The concept of improvement graph was also used by Talluri [1996] and Ahuja et al [2001b] to search very largescale neighborhoods as in flee t assignment problems arising in airline scheduling. Ergun [2001] also proposed seve ral improvement graphs for the vehicle routing problem and mach ine scheduling problems. In this chapter, we study the use of the improvement graph for the multiexchange neighborhood structure for th e QAP. However, our curre nt application of the improvement graph is different than previous applications. In previ ous applications, the improvement graph satisfied the property that the cost of the multiexchange was equal to the cost of the corresponding (constrained) cycle in the improvement graph. This property is not ensured for the improvement graph for the QAP. Rather, the cost of the cycle is a very good approximation of the cost of the multiexchange, and allows us to enumerate good neighbors quickly. The improvement graph also allows us to evaluate the cost of a neighbor faster than usi ng a normal method. Typically, evaluating a k exchange neighbor for the QAP takes O( nk ) time; but using the improveme nt graph we can do it in O( k ) average time per neighbor. We developed a generic search pro cedure to enumerate neighbors using improvement graphs. We also developed seve ral implementations of the generic search PAGE 35 24 procedure, which enumerate the neighborhoods exactly as well as heuristically. We present a detailed computational investiga tion of local improvement algorithms based on our neighborhood search structures. Our inve stigations yield the following conclusions: (i) locally optimal solutions obtained using multiexchange neighborhood search algorithms are superior to those obtained using 2ex change neighborhood search algorithms; (ii) generally, in creasing the size of the nei ghborhood structure improves the quality of local optimal solutions but after a certain point there are diminishing returns; and (iii) enumerating a restricted subset of neighbors is much faster than enumerating entire neighborhood and can develop improvements that are almost as good. This chapter is organized as follows. In Section 2.2, we describe the improvement graph data structure for the QAP. We present a generic heuristic search procedure for the K exchange neighborhood structure for the QAP in Section 2.3. In Section 2.4, we describe several specific implementations of the generic search procedure. In Section 2.5, we describe the neighborhood search algorithm based on the generic search procedure. Section 2.6 describes an acceleration technique we use to speed up the performance of the algorithm. For clarity of concepts, in a ll theses sections we analyze and discuss the algorithms for symmetric cases only, which can be easily generali zed for asymmetric cases. To keep the expressions simple, we assu me that diagonal elements are zero either in flow matrix or cost matrix, and hence do not account completely the interaction among diagonal elements. In Section 2.7, we pres ent the corresponding e xpressions for the general case (both symmetric and asymmetric instances). We provi de and analyze the computational results from our implementa tions in Section 2.8. Section 2.9 summarizes our contributions. PAGE 36 25 2.2 Improvement Graph One of the main contributions of th is chapter is the development of the improvement graph to enumerate multiexchanges for the QAP. In this section, we describe how to construct the improvement graph, and how it may help us in evaluating multiexchanges quickly. This section as well as the following sections requires some network notations, such as cycles and paths. We will use the graph notation given in the book by Ahuja, Magnanti, and Orlin [1993] and refer the reader to this book for the same. Given a permutation and a k exchange C we denote the cost of the cyclic exchange by Cost (, C ). This cost term represents the difference between the objective function values of C and that is, Cost (, C ) = z (C) z () = 2 ()() ()() 1CCn ijij ij iCjfdd 2.3 Clearly, the cost of the k exchange C can be computed in time O( kn ). We will show that using improvement graphs, the cost of C can be computed in O( k2) time. This time can be further reduced to an average of O( k ) time. Since we choose k to be much smaller than n the improvement graph allows us to eval uate multiexchange substantially faster than standard methods. In fact, it also leads to dramatic improvements in the running time to identify traditional 2exchanges. We associate an improvement graph G () = ( N A ) with which is a directed graph comprising of the node set N and arc set A The node set N contains a node i for every facility i and the arc set A contains an arc ( i j ) for every ordered pair of nodes i and j in N Each multiexchange C = i1 i2  ik i1 defines a (directed) cycle i1 i2  PAGE 37 26 ik i1 in G () and, similarly, eac h (directed) cycle i1 i2  ik i1 in G () defines a multiexchange i1 i2 ik i1 with respect to Thus, there is onetoone correspondence between multiex changes with respect to and cycles in G (). We will, henceforth, use C to denote both a multiexchange and a cycle in G (), and its type will be apparent from the context. An arc ( i j ) A signifies that the facility i moves from its current location to the current location of facility j In view of this interpretation, a cycle C = i1 i2 ik i1, signifies the following changes: facility i1 moves from its current location to the location of facility i2, facility i2 moves from its current locati on to the location of facility i3, and so on. Finally, facility ik moves from its current location to the location of facility i1. We now associate a cost ijc with each arc ( i j ) A Ideally, we would like to define arc costs so that the cost of the multiexchange C with respect to the permutation is equal to the cost of cycle C in G (). However, such a possib ility would imply that P = NP because the multiexchange neighborhood struct ure includes all feasible solutions for an instance of the QAP. We will, instead, define arc costs so that the cost of the multiexchange is close to the cost of the corresponding cycle. We define ijc as follows: it is the change in the cost of the solution when facility i moves from its current location to the location of facility j and all other facilities do not m ove. Observe that this change indicates that after the change there is no facility at location ( i ) and the location ( j ) has two facilities. Thus, to determine the cost of the change, we need to take the difference between the costs of interactions between facility i and other facilities before and after the change. Let denote the solution after the change. Then, ( l ) = ( l ) for l i and ( i ) PAGE 38 27 = ( j ). Note that is not a permutation because ( i ) = ( j ). We define ijc= z () z (). Then, ijc = z () z () = 2 ()()()() 1 n iljlil lfdd, 2.4 which captures the change in the cost of interaction from facility i to other facilities. The manner in which we define arc costs in the improvement graph does not ensure that the cost of the cycle C in G (), given by (,)C ij ijc will equal Cost (, C ). The discrepancy in these two cost terms aris es because when defining the arc cost ijc we assume that the facility i moves from its current locati on to the location of facility j but all other facilities do not move But in the multiexchange C several facilities move and we do not correctly account for the cost of flow between facilities in C We, however, correctly account for the cost of flow between any two f acilities if one of the two facilities is not in C We show next that the cost term Cost (, C ) can be computed by adding a corrective term to (,)C ij ijc Cost ( C) = z (C) z () = 2 ()() ()() 1CCn ijij ij iCjfdd = 2 ()() ()()Cijij ij iCjCfdd + 2 ()() ()()CCijij ij iCjCfdd = 2 *(,)C ij ijc 2 ()() ()()Cijij ij iCjCfdd + 2 ()() ()()CCijij ij iCjCfdd 2.5 The equation 2.5 shows that we can determine the cost of a multiexchange C by first determining the cost of the cycle C in G (), which is (,)C ij ijc and then correcting PAGE 39 28 it using the second term given in equation 2.5. This corrective term can be computed in time O( k2) if C is a k exchange. Let us now remark on the usefulness of the improvement graph. First, it allows us to determine the approximate cost of a multiexchange C quickly. The cost term (,)C ij ijc is a reasonable estimate of the cost of the multiexchange C To see this, observe from equation 2.5 that the corrective term Cost (, C (,)C ij ijc contains O( k2) interactions between facilities. However, n facilities have O( n2) interactions between them. If we choose k to be a relatively small fraction of n then the corrective term (on the average) will be substantially smaller than the total cost and the cost of the cycle C in G () will be a good estimate of th e cost of the multiexchange C For example, if n = 100 and k = 5 then there are 9,900 interactions between facilities and only 20 of them are counted incorrectly. If we use k = 10 then about 100 of them are counted incorrectly which is only 1% of the total interactions between facilities. Thus, the improvement graph allows us to enumerate extremely larg e set of neighbors quickly using approximate costs, and the approximation in costs is quite small. The improvement graph also allows us to determine the correct cost of a multiexchange faster than it normally takes to co mpute its cost. Normally, to compute the cost of a multiexchange takes O( kn ) time as we would need to update the cost interactions between k facilities (that move) w ith other facilities. However, using equation 2.5 we can compute the cost of a multiexchange in O( k2) time. For example, if n = 100 and k = 10 then we can compute the cost of a multiexchange about 10 times faster which can make substantial difference in an algorithms performance. PAGE 40 29 The benefits we derive from the use of improvement graph come at a cost; we need to construct the improvement graph and ca lculate arc costs. It follows from equation 2.4 that we can construct the impr ovement graph from scratch in O( n3) time. But we need to compute the improvement graph from scratc h just once. In all subsequent steps, we only update the improvement graph as we perform multiexchanges. We show in the next lemma that updating the impr ovement graph following a k exchange takes only O( kn2) time. We also show in Section 2.8 that our neighborhood search algorithms need to use k (4 and 5 ) only as on the benchmark tests higher values do not add extra benefit. Hence, it takes O( n2) time to update the improvement graph, wh ich is quite efficient in practice. Thus, the time needed to construct and update the improvement graph is relatively small, and is well justified by the savings we obtain in enumerating and evaluating multiexchanges. Lemma 1: Given the improvement graph G() and a kexchange C with respect to the improvement graph G(C) can be constructed in O(kn2) time. Proof: The improvement graphs G () and G (C) have the same set of nodes and arcs. They differ only in arc costs. Each arc ( i j ) G (C) is one of the following two types: (i) either i C or j C and (ii) i C and j C. There are 2 k ( n k ) = O( nk ) arcs of type (i), and O( n2) arcs of type (ii). Using equation 2.5, we can determine the cost of a type (i) arc in O( n ) time, thus giving a total time of ( n2k ) to compute the cost of all type (i) arcs. We show next that we can dete rmine the cost of a type (ii) arc in O( k ) time, which also yields a total time of O( n2k ) to compute the costs of all type (ii) arcs. Cijc= 2 ()()()() 1CCn il j lil lfdd PAGE 41 30 = 2 ()()()() iljlil lCfdd + 2 ()()()()CCil j lil lCfdd = 2 ()()()() 1 n iljlil lfdd 2 ()()()() iljlil lCfdd + 2 ()()()()CCil j lil lCfdd = 2 ijc 2 *()()()() iljlil lCfdd+ 2 ()()()()CCil j lil lCfdd. 2.6 Since we already know ijc and C is a k exchange, we can eval uate equation 2.6 in O( k ) time, which establishes the lemma. 2.3 Identifying Profitable MultiExchanges Our algorithm for the QAP is a neighbor hood search algorithm and proceeds by performing profitable multiexchanges. To keep the number of multiexchanges enumerated manageable, we first enumerate 2exchanges, followed by 3exchanges, and so on, until we reach a specified value of k denoted by K which is the largest size of the multiexchanges we wish to perform. This enumeration scheme is motivated by the consideration that we look for larger size multiexchanges when smaller size multiexchanges cannot be found. In this section, we describe a generic search procedure for enumerating and identifying multiexchanges using improvement graphs. Our method for enumerating multiexc hanges with respect to a solution proceeds by enumerating directed paths of increasing lengths in the improvement graph G (), where the length of a path is the numbe r of nodes in the path. Observe that each path P = i1 i2 ik in the improvement graph has a corresponding cycle in the improvement graph i1 i2 ik i1 obtained by joining the last node of the path with PAGE 42 31 the first node in the path; this cycle also defines a multiexchange with respect to Let C ( P ) denote the multiexchange defined by the path P Our method for enumerating cycles of in creasing lengths performs the following three steps repeatedly for increasing values of k starting with k = 2. Let Sk denote a set of some paths of length k in G (). We start with S1 = { 1 2 , n }, which is the set of n paths of length 1, each consisting of a singleton node. Path extension: We consider each path P Sk1 one by one and extend it by adding one node to it. To extend a path P = i1 i2 ik 1, we add the arc ( ik 1, ik) for each ik N \{ i1, i2, , ik}, and obtain several paths of length k Let E ( P ) denote the set of all paths obtained by extending the path P Further, let k =1PSE(P)k Cycle evaluation: Each path P k yields a corresponding multiexchange C ( P ). We evaluate each of these multiexchanges and determine whether any of them is a profitable multiexchange. If yes, we return the best multiexchange and stop; otherwise we proceed further. Path pruning: In this step, we prune several paths in the set k, which are less likely to lead to profitable multiexchanges. We call a procedure, PathSelect ( k), that takes as an input the set of paths k enumerated in the previous step and selects a subset Sk of it. This subset of paths will be extended further in the next iteration for the next higher value of k We describe in Section 2.4 several ways to implement the PathSelect procedure. Path pruning is cr itical to keep the number of paths enumerated manageable. The following algorithmic description su mmarizes the steps of our heuristic search procedure, which we call the Kexchange search procedure. PAGE 43 32 procedure Kexchange search; begin k 1; let S1 N be the set of paths of length 1; C* and W* 0; while Sk is nonempty and k < K and W* 0 do begin k k+1; k 1PSE(P)k ; let Pmin Sk be the path such that Cost(,C(Pmin)) = min{Cost(,C(P)): P k}; if W* > Cost(,C(Pmin)) then W* Cost(,C(Pmin) and C* C(Pmin); Sk PathSelect(k); end; return C*; end. Figure 2.2 The generic search procedure for identifying profitable multiexchanges. Observe that in this procedure, the value of K is a parameter a nd can be specified by the user. Increasing the value of K may in general improve th e quality of local optimal solutions obtained, but our computational inve stigations show that there are diminishing returns after K = 4; hence K = 4 is a good value to be used in the search procedure. For another implementation (Implementation 4) of PathSelect as discussed in Section 2.4, we keep the value of K = 5. Also observe that the algorithm terminates in two ways: C* is empty or C* is nonempty. If C* is empty, then it implies that the algorithm has failed to find a profitable multiexchange and the current solution is locally optimal. If C* is nonempty, then it implies that the algo rithm found a profitable multiexchange C*. PAGE 44 33 We now analyze the complexity of the algorithm. Let p denote the maximum number of paths in any Sk. The while loop executes at most K times. In each execution of the while loop, it takes O( pn ) time to compute the set k and it may contain as many as pn paths. Since computing the cost of k exchange for each P k takes O( k2) time, we require O( k2pn ) time to find a profitable k exchange, if any. We shall show in Section 2.4 that the subroutine PathSelect takes O( pn log( pn )) time. Since for most situations considered by us log( pn ) < k2, the running time of the algorithm is O( k2pn ). It is easy to see that if we ignore the time taken by the procedure PathSelect then the bottleneck operation in the generic sear ch procedure is to evaluate the cost Cost (, C ( P )) of each path P k. Since C ( P ) is a k exchange with respect to the solution using equation 2.6 we can determine its cost in O( k2) time. We will next show that we can determine the cost of k exchange C ( P ) in O( k ) time. The generic search procedure pr oceeds by enumerating paths in G (). Each path P = i1 i2 ik in G () defines a path exchange with respect to the solution in an obvious manner, which is the same as the k exchange C = i1 i2 ik i1 except that we do not perform the last move of shifting facility ik from its current location to the location of facility i1. Alternatively, P( il) = ( il +1) for all l = 1, 2, , k 1, and P( i ) = ( i ) for all i N \ { i1, i2, , ik 1}. We denote the cost of the path exchange P with respect to the solution by Cost (, P ). Hence, Cost (, P ) = z (P) z () = 2 PP()() ()() P1 n ijij ij ijfdd 2.7 PAGE 45 34 Observe that P and C ( P ) differ only in the location of the facility ik. This observation allows us to compute th e cost of the cyclic exchange C ( P ) from the cost of the path exchange P in O( k ) time using the following expression: Cost (, C ( P )) Cost (, P ) = 2 1kiic + 2 CP 1 1()()()() ()()()() Ckk kijijij ijij jfdddd 2.8 Now suppose that we extend the path P to P = i1 i2 ik ik +1 by adding the node ik +1. Then, we can determine the cost of the path P from the cost of the path P in O( k ) time using the following expression: Cost (, P) Cost (, P ) = 2 1 kkiic + 2 PP 1 1()()()() ()()()() Pkkk kkijijij ijij jfdddd 2.9 In our enhanced version, we maintain the cost of each path P enumerated by the algorithm. Given the cost of path P we can determine the cost of the cycle C ( P ) in O( k ) time. Further, when we extend any path P then the cost of the extended path too can be computed in O( k ) time. Thus, the running time of the generic search procedure is O( K 2K k k P), plus the time taken by the subroutine PathSelect 2.4 Specific Implementations In Section 2.3, we presented a generic sear ch algorithm to identify a profitable multiexchange. We can derive several specific implementations of the generic version by implementing the procedure PathSelect ( k) differently. The procedure PathSelect ( k) PAGE 46 35 accepts as an input a set of paths k and returns a subset Sk of these paths. We describe next several ways in which PathSelect can be implemented. Implementation 1 (all paths): In this version, we define PathSelect ( k) to be k itself; that is, we select all the paths to be taken to the next stage. This version guarantees that we will always find a profitable multiexch ange if it exists. However, the number of paths enumerated by the algorithm increase exponentially with k and it takes too long to find profitable k exchanges for k 6 even for n = 25. Implementation 2 (negative paths): In this version, the subroutine PathSelect ( k) returns only those paths which have negative cost; that is, PathSelect ( k) = { P k: Cost (, P ) < 0} where is the current solution. This vers ion is motivated by the intuition that if there is a profitable multiexchange C = i1 i2 ik i1, then there should exist a node in this sequence, say node il, so that each of the paths il il +1, il il +1 il +2, , il il +1 il +2 il + k has a negative cost. Though results of this type are valid for many combinatorial optimization problems, it is not true for the QAP. However, it is a reasonable heuristic to eliminate paths that are less likely to yield profitable multiexchanges. Implementation 3 (best n2 paths): In this version, we sort all the paths in k in the nondecreasing order of path costs, and select the first n2 paths, where is a specified constant. For example, if = 2, then we select the best 2 n2 paths. This version is motivated by the intuition that the paths with lower cost are more likely to yield profitable multiexcha nges. The choice of allows us to strike a right tradeoff between the running time and the solution quality. Higher values of will increase the chances of finding profitable multiexchanges but also incr ease the time needed to find a profitable PAGE 47 36 multiexchange. Our computational results presented in Section 2.8 indicate that = 1 is a good choice considering both the running time and solution quality. We have used max heap data structure to keep n2 paths in a stage. Hence if there are pn possible paths (as discussed in Section 2.3), it takes pn log( pn ) time to store n2 best paths in a heap. Implementation 4 (best n paths): In this implementation, we select the best path in k starting at node i for each 1 i n Therefore, the set Sk contains at most one path starting at each node in N. No te that in Implementation 3, it is possible that many low cost paths contain the same set of arcs making the search less diverse. Allowing each node to be the starting point of a different path can add some diversity to the heuristic search process. We implemented Implementation 3 and 4 of PathSelect on the test problems as our experiment on the test problems indicates that it is difficult to use Implementation 1 even for moderate size problems and in all case s Implementation 3 gives better result than Implementation 2. We present in Section 2.8 the computational results of these implementations. 2.5 The Neighborhood Search Algorithm In this section, we describe our neighbor hood search algorithm (Figure 2.3) for the QAP. Our algorithm starts with a rando m permutation (obtained by generating pseudorandom numbers between 1 and n and rejecting the numb ers already generated) and successively improves it by performing profitable multiexchanges obtained by using the K exchange search procedure, until the proc edure fails to produce a profitable multiexchange. PAGE 48 37 algorithm QAPneighborhoodsearch; begin generate an initial random permutation ; construct the improvement graph G(); while Kexchange search returns a nonempty multiexchange C do begin replace the permutation by the permutation C; update the improvement graph; end; return the permutation ; end; Figure 2.3 The neighborhood search algorithm for the QAP. Let us perform the running time analysis of the algorithm. The initial construction of the improvement graph takes O( n3) time. The time needed by the procedure K exchange search is O( K2p ), where p is the maximum number of paths maintained by the procedure during any iteration (see S ection 2.3). For Implementation 3 of PathSelect p n2 and this procedure requires O( K2n2) time per iteration (that is, per improvement). For the Implementation 4 of the PathSelect p n2, and the procedure again takes O( K2n2) time. Updating the improvement graph takes O( n2K ) time (see Section 2.2). Each execution of the QAPneighborhoodsearch algorithm yields a locally optimal solution of the QAP with resp ect to the neighborhood defined by the Kexchange search procedure. The solution obtained de pends upon the initial random permutation and the version of the PathSelect procedure we use. We refer to one execution of the algorithm as one run Our computational investigations re vealed that if we apply only one run of the algorithm, then the solution method is not very robust. The QAP in general has an extremely large number of locally optim al solutions even if the size of the neighborhood is very large. Each run produces a locally optimal solution, which is a PAGE 49 38 random sample in the solution space of lo cally optimal solutions To obtain a robust locally optimal solution, we need to perfor m several runs of the algorithm and use the best locally optimal solu tion found in these runs. 2.6 Accelerating the Search Algorithm In this section, we describe a method to speedup the performance of the generic search algorithm and also its specific implementations. The speedup uses the fact that several paths give the same multiexchange. For example, all the paths i1 i2 i3 i4, i2 i3 i4 i1, i3 i4 i1 i2, and i4 i1 i2 i3 imply the same multiexchange i1 i2 i3 i4 i1 when we connect the last node of these path s to the first node of the path. In general, a k exchange can be represented by k different paths. Since our generic search algorithm enumerates k exchanges by enumerating paths, we may obtain the same k exchange several times during the search process th rough different paths. To avoid repeated enumeration of the multiexchanges, our sear ch algorithm maintains certain kinds of paths, called valid paths defined as follows: Valid Paths : A path i1 i2  ik is a valid path if i1 ij for every 2 j k. Our generic search algorithm enumerates only valid paths. The following lemma shows that we do not miss any multiexc hanges by maintaini ng valid paths only. Lemma 2 Any multiexchange can be enumerated by maintaining only valid paths. Proof : Consider a multiexchange j1 j2  jk. Let jl = min { jh: 1 h k }. Now define i1 = jl, i2 = jl +1, , ik = jl + k, where all subscript ma thematics is modulo ( k +1). It follows from the definition of jl that each of the paths i1, i1 i2, i1 i2 i3, i1 i2  ik is a valid path. He nce starting at node i1 we can gradually build i1 i2  ik by PAGE 50 39 maintaining only valid paths, and joining node ik to node i1 gives us the desired multiexchange. We can easily modify the generic sear ch algorithm so that it only enumerates valid paths. In this modified algorit hm, when we consider adding the arc ( ik, ik +1) to the path i1 i2  ik, we compare i1 with ik +1. If i1 ik +1, we add the arc; otherwise we do not add it. It can be noted that above lemma holds if we enumerate all paths. However, as we keep only n2 paths in each stage, there may be the cases when we might miss a profitable multiexchange. Our experiment s hows that losses in missed improvements are well compensated by the gain in time. The co mputational results pr esented in Section 2.8 show that enumerating only valid paths decr eases the running time of the generic search algorithm substantially. 2.7 Expressions for the Asymmetric QAP In the previous sections, we gave expr essions for calculating various cost terms for symmetric instances of QAP. In this sec tion, we give expressions for the asymmetric QAO. We state the expressions without proof since the logic is similar to those for the symmetric case. For the asymmetric case, we will replace the expressions 2.3 2.9 by the following expressions 2.10a 2.10g respectively. Cost(, C) = z(C) z() = ()()()() ()()()() 1CCCCn ijijjiji ijji iCjfddfdd 2.10a ijc = z() z() =()()()()()()()() 1 n iljlilliljli lfddfdd. 2.10b PAGE 51 40 Cost(, C) = (,)C ij ijc ()()()() ()()()()CCijijjiji ijji iCjCfddfdd + ()() ()()CCijij ij iCjCfdd + ()() ()()CC j iji ji iCjCfdd 2.10c Cijc= ijc ()()()()()()()() iljlilliljli lCfddfdd + ()()()()()()()()CCCCilli jlilljli lCfddfdd 2.10d Cost(, P) = z(P) z() = PPPP()()()() ()()()() P1 n ijijjiji ijji ijfddfdd 2.10e Cost(, C(P)) Cost(, P) =1kiic+ CP 1 1()()()() ()()()() Ckk kijijij ijij jfdddd + CP 1 1()()()() ()()()() Ck k kjijiji jiji jfdddd 2.10f Cost(, P) Cost(, P) =1kkiic+ PP 1 1()()()() ()()()() Pkk kijijij ijij jfdddd + PP 1 1()()()() ()()()() Pk k kjijiji jiji jfdddd 2.10g 2.8 Computational Testing In this section, we describe computat ional results of th e neighborhood search algorithms developed by us. We implemented al l of our algorithms in C and ran them on IBM SP machine (model RS6000) with a proce ssor speed of 333 MHz. We tested the algorithms on 132 benchmark instances availa ble at the QAPLIB, the library of QAP instances maintained by the Institute of Mathematics, Gr az University of Technology (http://www.opt.math.tugraz.ac.at/qaplib/). Our computational results include analyzing the PAGE 52 41 CPU times taken by our algorithms, quality of the solutions obtained by them as well as understanding the behavior of the VLSN search algorithms. Neighborhood search algorithms need some feasible solution as the starting solution. We generated random permutations of n numbers and used them as starting solutions. Further we implemented a multistart version of the neighborhood search algorithm, where we apply the neighborhood search algorithm multiple times with different starting solutions, called different runs, and select the be st solution found in these runs that the number of runs de pend on the size of the problem instance. In Section 2.4, we proposed four implem entations of the generic VLSN search algorithm for the QAP. The first implementatio n maintains all the paths enumerated in the search process. We found that the nu mber of paths grows very quickly with k and the algorithm runs very slowly even when we go up to k exchanges with k = 6. For example to solve a QAP with n = 42 (problem sko42), each run of this implementation takes about 8 seconds for k = 4 whereas Implementation 3 takes only 0.025 seconds per run. Additional preliminary tests yielded that th is implementation is not as competitive as other implementations and we decided not to perform a thorough testing of the algorithm. In the second implementation of the VLSN search algorithm, we maintain only those paths, which have negative costs. Fo r many combinatorial optimization problems, maintaining only negative cost paths is sufficient to enumerate negative cost cycles (improved neighbors), but this is not true for the QAP due to the nonlin earity in the cost structure. Our computational testing revealed that maintaining only negative cost paths is not a good heuristic to enumer ate negative cost cycles. Thus, we did not perform a thorough testing of this implementation. PAGE 53 42 Our preliminary testing revealed that Im plementation 3 and 4 exhibited the best overall behavior and deserved a thor ough testing. The following details of implementation 3 are worth mentioning. Recall from Section 2.4 that we keep only n2 best paths in k. We used the Max Heap data struct ure (Cormen et al. [2001]) to store these paths. We found that = 1 gives fairly good results and hence used this value. In addition, we used only those paths whose path cost is not more than 0.5% of the best objective function value of the QAP found so fa r. We found that using higher cost paths rarely leads to negative cost cycles. Finally, when we examined paths in k to enumerate cycles of length k and find several negative cost cycl es, we use the least cost negative cycle to obtain the next solution. As far as Implementation 4 is concerned, we implemented it in the straightforward fashion but before enumerating paths, we eliminate all negative cycles of length 2 by performing 2exchanges. 2.8.1 Accuracy of the Solution We applied Implementation 3 and 4 to 132 benchmark instances in QAPLIB, of these 98 were instances of symmetric QAP and the remaining were for the asymmetric case. We applied multiple runs of each im plementation and ran them for a specified amount of time. For the symmetric instances, we ran our algorithm for 1 hour for n 40 and for 2 hours for n > 40. The running times for the as ymmetric instan ces were 1.5 hours for n 40 and for 3 hours for n > 40. Tables 1 & 2 in Appendix 1, respectively, give the results of these algorithms for symmetric and asymmetric instances and compare our solutions with the solutions obtained by the 2exchange algorithm (2OPT) and the best known solutions (BKS). The columns titled BestGap, AvgGap, nRuns, %Best, respectively, give the percent deviation of the best solution found in all runs with respect PAGE 54 43 to the best known solution, av erage deviation over solutions found in all runs, number of runs, and the percentage of the solutions found were best known solutions. User can derive the following conclusions from these tables. Implementation 3 exhibited the best overa ll performance. It obtained the best known solutions in 74 out of 98 symmetric instances and in 24 out of 34 asymmetric instances. Its av erage error was the lowest and it found the best known solutions with the maximum frequency. Implementation 3 is found to exhibit superi or performance compared to 2OPT in terms of gap of best solution found by al gorithm with the best known solution. For 25 symmetric instances implementation 3 obt ained better solutions than 2OPT, and for only 2 symmetric instances 2OPT obtained better solutions than Implementation 3. Similarly, for 10 asym metric instances implementation 3 obtained better solutions than 2OPT, a nd for only 1 asymmetric instance 2OPT obtained better solutions than Implementation 3. Implementation 3 was also found to be bette r than 2OPT in terms of average gap and frequency of finding best known so lution. The average of AvgGap of Implementation 3 was 7.6%, whereas th is number for 2OPT was 11.05% for symmetric instances and these numbers were 6.42% and 7.49% respectively for asymmetric instances. Finally, wher eas Implementation 3 found best known solution with an average frequency of 17.13% in symmetric case, this number for 2OPT was 11.85% in symmetric case. Fo r asymmetric case, Implementation 3 found best known solution with an aver age frequency of 1.97%, whereas this number for 2OPT is 0.46%. Implementation 4 also exhibited superior pe rformance with respect to 2OPT, but its overall performance was worse than Impl ementation 3. Implementation 4 runs very fast and it terminates in a fraction of second for most problem sizes, but the solutions obtained using this method are not as robust as those obtained using Implementation 3. Above results seem to suggest that very largescale neighborhood is overall more effective than the traditiona l 2exchange neighborhood. When both the algorithms are run for the same time, the 2OPT performs many mo re runs but still the best solution found is, on the average, not as good as found by VLSN se arch in fewer number of runs. Hence the extra time taken by VLSN search algorithm is more than justified by the better quality of the solutions obtained. PAGE 55 44 We will now describe some computational investigations we performed to understand the behavior of our implementations. 2.8.2 Effect of Neighborhood Size In our approach, the size of the nei ghborhood critically depends upon (i) the maximum cycle length, and (ii) the number of paths, of a given length, maintained. The larger the cycle length and th e number of paths maintaine d, larger is the neighborhood, more is the running time, and better is the quality of the solution obtained (in general). Hence it is worthwhile to examine the effect of these two paramete rs on the running time and the solution quality. In our first experiment, we considered six problems of the same size; sko100a, sko100b, sko100c, sko100d, sko100e, sko100f, and applied 100 runs of Implementation 3 with cycle lengths varying from 2 to 7 and noted the average running time taken by the algorithm (per run) and the average gap (per run). We kept the number of paths maintained by the algorithm as fixed at n2. Figure 2.4 plots these two values as a function of cycle length. It can be seen that the average gap decreases significantly with the increase in cycle length until cycle length is 4, and after that th e average gap does not vary much. We also observe that the running time of the algorithm in creases linearly with the increase in the cycle length. We think that the cycle length of 4 strikes a right balance between the solution accuracy and solution time and hence we used this value in the computational results presented earlier. PAGE 56 45 1.6 1.7 1.8 1.9 2 2.1 2.2 234567 Cycle LengthAvg. Gap (%)0 200 400 600 800 1000Time (Sec) A vera g e Ga p ( % ) Time ( Sec ) Figure 2.4 Effect of cycle length on ti me taken and solution quality for 100 runs on problem sko100af Our second experiment was similar to th e first experiment but we varied the number of paths maintained by the algorithm while keeping the cycle length fixed at 4. Figure 2.5 gives a plot of the average gap a nd average time per run when we performed 100 runs of Implementation 3 on the six problems sko100af. We observe that the solution accuracy gradually improves as the number of paths increase as well as the running time of algorithm incr eases linearly with the number of paths maintained. We believe that maintaining n2 paths is a good compromise between solution quality and solution time and we used this value in our experiments. In another experiment, we counted the number of improvement iterations with cycle length 2, 3 and 4. Recall th at our algorithm pe rforms a 3exchange when it fails to find 2exchange, and performs a 4exchange wh en it fails to find a 3exchange. Table 2.1 gives these values for 10 benchmark in stances on which we apply 100 runs of Implementation 3. We observe that there ar e many more iterations with 2exchanges compared to 3exchanges, and many more 3exchanges compared to 4exchanges. PAGE 57 46 1.73 1.74 1.75 1.76 1.77 1.78 1.79 1.8 1.81 0.5n2n21.5n22n22.5n23n2 # of pathsAvg. Gap (%)0 100 200 300 400 500 600 700 800Time(Sec) Avg Gap(%) Time(Sec) Figure 2.5 Effect of number of paths in each stage on time taken and solution quality for 100 runs on problem sko100af Table 2.1 Number of iterations with different implemented cycle length # of iterations of Cycle Length Problem 2 3 4 chr22a 1614 151 49 kra30a 2283 123 30 kra30b 2306 125 32 nug30 2580 104 42 ste36a 3537 142 44 tho40 3839 124 35 wil50 5298 81 38 sko42 4246 150 64 sko100a 13232 270 70 tai100a 7267 274 75 2.8.3 Effect of the Speedup Technique The reader may recall from Section 2.6 that we used a speedup technique to reduce redundant enumeration of cycles. In this t echnique, we maintain only those valid paths i1,i2, ., ik for which ik > i1. Lemma 2 showed that we would not miss any negative cycles even if we only maintain valid paths. This proof relied on the assumption that we maintain all valid paths. Since our algorithm maintains only n2 paths, we might miss PAGE 58 47 some negative cycles and the speedup tec hnique may deteriorate the quality of the solutions obtained. We performed an experi ment to assess the effect of the speedup technique on the solution quality and solution time. Table 2.2 gives these values for 10 benchmark instances. We applied 100 runs on each benchmark instance and obtained the average values. We observe that speedup t echnique not only decreases the running time substantially but also worsens the solution quality. We believe th at overall it is advantageous to use the speedup technique sin ce the saved time can be used to perform more runs of the algorithm and improve the overall performance of the algorithm. Table 2.2 Effect of accelerated path enumeration scheme Using Speedup Technique Without Speedup Technique Problem Average Gap Time in SecondsAverage GapTime in Seconds Chr22a 10.00 1 9.13 5 Kra30a 6.44 5 6.27 12 Kra30b 4.26 5 4.05 12 N ug30 3.19 5 2.92 12 Ste36a 9.34 10 8.37 27 Tho40 3.87 12 3.76 28 Wil50 1.54 24 1.41 70 Sko42 2.68 17 2.57 40 Sko100a 1.87 283 1.78 962 tai100a 2.88 280 2.48 1191 2.9 Conclusions In this chapter, we develop a very la rgescale neighborhood structure for the QAP. We show that using the concept of impr ovement graph, we can easily and quickly enumerate multiexchange neighbors of a gi ven solution. We develop a generic search procedure to enumerate and evaluate neighbors and propose several specific implementations of the generic procedur e. We perform extensive computational PAGE 59 48 investigations of our implementations and they suggest that multiexchange neighborhoods add value over the comm only used 2exchange neighborhoods. Our implementations of multiexchange neighborhood search algorithms are local improvement methods. We wanted the focus of our research effort more on neighborhood structure and less on specific implementati ons. Further possibilities for improvement could possibly be obtained us ing ideas from tabu search (G lover and Laguna [1997]). We leave it as a topic of future research. Neighborhood search algorithms have also been used in genetic algorithms to improve the quality of the individuals in the population (Ahuja, Orlin, and Tiwari [2000], and Drezn er [2003]). Our neighborhood structure may be useful in these genetic algorithms too. PAGE 60 49 CHAPTER 3 WEAPONTARGET ASSIGNMENT PROBLEM 3.1 Introduction The WeaponTarget Assignment (WTA) problem is a fundamental problem arising in defenserelated applications of operations research. The problem consists of optimally assigning weapons to the enemytargets so that the total expected surv ival value of the targets after all the engagements is mini mized. There are two versions of the WTA problem: static and dynamic In the static version, all the inputs to the problem are fixed; that is, all targets are known, all weapons are known, and all weapons engage targets in a single stage. The dynamic version of the pr oblem is a multistage problem where some weapons are engaged at the targ ets at a stage, the outcome of this engagement is assessed and strategy for the next stag e is decided. In this chapter, we study the static WTA problem; however, our algorithms can be used as important subroutines to solve the dynamic WTA problem. We now give a mathematical formulation of the WTA problem. Let there be n targets, numbered 1, 2, , n and m weapon types, numbered 1, 2, , m Let Vj denote the value of the target j and Wi denote the number of weapons of type i available to be assigned to targets. Let pij denote the probability of destroying target j by a single weapon of type i Hence qij= 1 pij denotes the probability of survival of target j if a single weapon of type i is assigned to it. Observe that if we assign xij number of weapons of type i to target j then the survival probability of target j is given by ijx ijq A target may be PAGE 61 50 assigned weapons of different types. The WT A problem is to determine the number of weapons xij of type i to be assigned to target j to minimize the total expected survival value of all targets. This problem can be formulated as the following nonlinear integerprogramming problem: Minimize 1 x ij m n j j ij i1Vq 3.1a Subject to 1 n iji j x W, for all i = 1, 2, , m 3.1b ij x 0 and integer, for all i = 1, 2, , m, and for all j = 1, 2, , n. 3.1c In the above formulation, we minimize the expected survival va lue of the targets while ensuring that the total number of weapons used is no more than those available. This formulation presents a simplified vers ion of the WTA problem. In more practical versions, we may consider adding additional co nstraints, such as (i) lower and/or upper bounds on the number of weapons of type i assigned to a target j; (ii) lower and/or upper bounds on the total number of w eapons assigned to target j; or (iii) a lower bound on the survival value of the target j. The algorithms proposed in th is chapter can be easily modified to handle these additional constraints. Research on the WTA problem dates back to the 1950s and 1960s where the modeling issues for the WTA problem was i nvestigated (Manne [1958], Braford [1961], Day [1966]). Lloyd and Witsenhausen [1986] established the NPcompleteness of the WTA problem. Exact algorithms have been pr oposed to solve the WTA problem for the following special cases: (i) when all the wea pons are identical (DenBroder et al. [1958] and Katter [1986]) or (ii) when the targets ca n receive at most one weapon (Chang et al. PAGE 62 51 [1987] and Orlin [1987]). Some of the heuris tics proposed to solv e the WTA problem are based on nonlinear network flow (Castanon et al. [1987]), neural networks (Wacholder [1989]), and genetic algorithms (Grant et al [1993]). Green et al. [1997] applied a goal programmingbased approach to the WTA problem. Metler and Preston [1990] have studied a suite of algo rithms for solving the WTA problem efficiently, which is critical for realtime applications of the WTA problem. Maltin [ 1970], Eckler and Burr [1972] and Murphey [1999] provide comprehensiv e reviews of the literature on the WTA problem. Research to date on the WTA probl em either solves the WTA problem for special cases or develops heuristics for the WTA problem. Moreover, since no exact algorithm is available to solve the weapon target assignment problems, it is not known how accurate are the solutions obtaine d by these heuristic algorithms. In this chapter, we propose several exact and heuristic algorithms to solve the WTA problem. Our branch and bound algorithms are the first implicit enumeration algorithms that can solve moderate size instances of the WTA problem optimally. We also propose heuristic algorithms, which generate almost optimal solutions within a few seconds. This chapter makes the following contributions: We formulate the WTA problem as an inte ger linear programming problem, that is, as a generalized integer network flow problem on an appropriately defined network. The linear programming relaxati on of this formulation gives a lower bound on the optimal solution of the WTA pr oblem. We describe this formulation in Section 3.2.1. We propose a minimum cost flow formula tion that yields a different lower bound on the optimal solution of the WTA problem This lower bound is, in general, not as tight as the bound obtained by the lin ear programming form ulation described above but it can be obtained in much le ss computational time. We describe this formulation in Section 3.2.2. We propose a third lower bounding scheme in Section 3.2.3, wh ich is based on simple combinatorial arguments and uses a greedy approach to obtain a lower bound. PAGE 63 52 We develop branch and bound algorithms to solve the WTA problem employing each of the three bounds described above These algorithms are described in Section 3.3. We propose a very largescale neighborhood (VLSN) search algorithm to solve the WTA problem. The VLSN search algorithm is based on formulating the WTA problem as a partition problem The VLSN search starts wi th a feasible solution of the WTA problem and performs a sequence of cyclic and path exchanges to improve the solution. We desc ribe in Section 3.4 a heuristic method that obtains an excellent feasible solution of the WTA problem by solving a sequence of minimum cost flow problems and then uses a VLSN search algorithm to iteratively improve this solution. We perform extensive computational inve stigations of our algorithms and report these results in Section 3.5. Our algorithms solve moderately large size instances (up to 80 weapons and 80 targets) of the WTA problem optimally and obtain almost optimal solutions of fairly la rge instances (up to 200 weapons and 200 targets) within a few seconds. 3.2 Lowerbounding Schemes In this section, we describe four lower bounding schemes for the WTA problem, using linear programming, integer programmi ng, minimum cost flow problem and a combinatorial method. These four approaches produce lower bounds with different values and have different running times. 3.2.1 A Lower Bounding Scheme using an Integer Generalized Network Flow Formulation In this section, we formulate the WT A problem as an integerprogramming problem with a convex objective function value. This formulation is based on a result reported by Manne [1958] who attributed it to Dantzig (personal communications). In formulation 3.1, let sj = 1 x ij m i ijq. Taking logarithms on bot h sides, we obtain, log(sj) = 1()m ijij i x logqor log(sj) = 1()m ijij i x logq Let yj = log(sj) and dij = log(qij). Observe that since 0 qij 1, we have dij 0. Then yj = 1 m ijij idx Also observe PAGE 64 53 that 1 x ij m i ijq= jy2. By introducing the terms ijdand yj in formulation 3.1, we get the following formulation: Minimize 1jn y j jV2 3.2a Subject to n j ijx1 Wi for all i = 1, 2, ..., m, 3.2b m i ij ijx d1 = yj for all j = 1, 2, ..., n, 3.2c ij x 0 and integer for all i = 1,..., m and for all j = 1,..., n, 3.2d yj 0 for all j = 1, 2, ..., n. 3.2e Observe that formulation 3.2 is an inte gerprogramming problem with separable convex objective function. This integer program can also be viewed as an integer generalized network flow problem with conve x flow costs. Generalized network flow problems are flow problems wher e flow entering an arc may be different than the flow leaving the arc (see, for example, Ahuja, Ma gnanti, and Orlin [1993]). In a generalized network flow problem, each arc (i,j) has an associated multiplier ij and the flow xij becomes ijxij as it travels from node i to node j. The formulation 3.2 is a generalized network flow problem on the network show n in the Figure 3.1. We give next some explanations of this formulation. PAGE 65 54 Figure 3.1 Formulating the WTA proble m as Integer generalized network flow problem The network contains m weapon nodes, one node corresponding to each weapon type. The supply at node i is equal to the number of weapons available, Wi for the weapon type i. The network contains n target nodes, one node corresponding to each target. Further, there is one sink node t whose demand equals the sum of all supplies. The supplies/demands of target nodes are zero. We now describe the arcs in the network. The network contains an arc conn ecting each weapon node to each target node. The flows on these arcs are given by xij, representing the number of weapons of type i assigned to the target j. The multipliers for these arcs are dijs. Since there is no cost coefficient for xijs in the objective function, the co st of flow on these arcs is zero. The network contains an arc from each of the target nodes to the sink node t. The flow on arc (j, t) is given by yj and the cost of flow on this arc is Vj.jy2 Finally, there is a loop arc (t, t) incident on node t with multiplier . An appropria te flow on this arc is sent so as to satisfy the mass balance constraints at node t. In formulation 3.2, the cost of the flow in the network equals the objective function 3.2a; the mass balance constraints of weapon node s are equivalent to the constraint 3.2b; : 1 2 m Weapons 1 2 n Targets t xij i : j : : yj W1 W2 W3 Wm dij 1/2 PAGE 66 55 and mass balance constraints of target nodes are equivalent to th e constraint 3.2c. It follows that an optimal solution of the above generalized network flow problem will be an optimal solution of the WTA problem. The generalized network flow formulation 3.2 is substantially more difficult than the standard generalized netw ork flow problem (see Ahuja et al. [1993]) since the flow values xijs are required to be integer numbers (ins tead of real numbers ) and the costs of flows on some arcs is a convex function (instead of a linear function). We will approximate each convex function by a piecewise linear convex function and relax the integer flows by realvalued flows so that the optimal solution of the modified formulation gives a lower bound on the optimal solution of the generalized formulation 3.2. We consider the cost function jy jV2 at values yj that are integer multiples of a parameter p > 0, and draw tangents of jy jV2 at these values. Let Fj(p, yj) denote the upper envelope of these tangents. It is easy to see that the function Fj(p, yj) approximates jy jV2from below and for every value of yj provides a lower bound on .jy jV2 Figure 3.2 shows an illustration of this approximation. PAGE 67 56 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1012345678910jy 2 jy Figure 3.2 Approximating a convex function by a lower envelope of linear segments. Thus, in formulation 3.2 if we replace the objective function 3.2a by the following objective function: 1(,)n jj jFpy, 3.2a we obtain a lower bound on the optimal objectiv e function of 3.2a. Using this modified formulation, we can derive lower bounds in two ways: LP based lower bounding scheme : Observe that the preceding formulation is still an integer programming problem because are flows xijs are required to be integer valued. By relaxing the integrality of the xijs, we obtain a mathema tical programming problem with linear constraints and pi ecewise linear convex objective functions. It is wellknown (see, Murty [1976]) that linear programs w ith piecewise linear conve x functions can be PAGE 68 57 transformed to linear programs by introducing a variable for every linear segment. We can solve this linear programming problem to obtain a lower bound for the WTA problem. Our computational results indicate the lower bounds generated by this scheme are not very tight. MIP based lower bounding scheme : In this scheme, we do not relax the integrality of the xijs, which keeps the formulation to be an integer programming formulation. We, however, tr ansform the piecewise linear c onvex functions to linear cost functions by introducing a variable for every li near segment. We then use cutting plane methods to obtain a lower bound on the optimal objective function value. We have used the builtin routines in the software CPLEX 8.0 to gene rate Gomory and mixed integer rounding cuts to generate fairly ti ght lower bounds for the WTA problem. We summarize the discussion in this section as follows: Theorem 1. Both the LP and MIP based lo wer bounding schemes give a lower bound on the optimal objective func tion value for the WTA problem. 3.2.2 A Minimum Cost Flow Base d Lower Bounding Scheme The objective function of the WTA problem can also be interpreted as maximizing the expected damage to the targets. In this section, we deve lop an upper bound on the expected damage to the targets. Subtract ing this upper bound on the expected damage from the total value of the targets (that is, n j jV1) will give us a lower bound on the minimum survival value. We will formulate the problem of maximizing the damage to targets as a maximum cost flow pr oblem. We show the underlying network G for the maximum cost flow formulation in Figure 3.3. PAGE 69 58 Figure 3.3 A network flow form ulation of the WTA problem This network has three laye rs of nodes. The first la yer contains a supply node i for every weapon type i with supply equal to Wi. We denote these supply nodes by the set N1. The second layer of n odes, denoted by the set N2, contains nodes corresponding to targets, but each target j is represented by several nodes j1, j2, , jk, where k is the maximum number of weapons that can be assigned to target j. A node jp represents the pth weapon striking the target j. For example, th e node labeled 31 represents the event of the first weapon being assigned to target 3, the node labeled 32 represents the event of the second weapon being assigned to target 3, and so on. All nodes in th e second layer have 1 2 3 11 Weapons Targets 1213 212223313233 t W3 W2 W1 m i i=1W arc capacity = 1 arc cost = c(i, jk) arc capacity = 1 arc cost = c(jk, t) PAGE 70 59 zero supplies/demands. Finally, the th ird layer contains a singleton node t with demand equal to 1m i iW. We now describe the arcs in this ne twork. The network contains an arc (i, jk) for each node iN1 and each node jk N2; this arc represents the assignment of a weapon of type i to target j as the kth weapon. This arc has a unit capacity. This network also contains an arc (jk, t) with unit capacity for each node jk N2. We call a flow x in this network a contiguous flow if it satisfies the property that if x(i,jk) = 1, then x(i,jl) = 1 for all l=1,2, ,k1. In other words, the contiguous flow implies that a weapon i is assigned to target j as the kth weapon provided that (k1) weapons have already been as signed to it. The following resu lt directly follows from the manner we have constructed the network G: Observation 1 There is onetoone correspondence between feasible solutions of the WTA problem and contiguous flows in G. While there is a onetoone correspondence between feasible so lutions, it is not a cost preserving correspondence if we require costs to be linea r. We instead provide linear costs that will overestimate the true nonlinear costs. We define our approximate costs next. The arc (i, jk) represents the assignmen t of a weapon of type i to target j as the kth weapon. If k = 1, then the cost of this arc is the damage caused to the target: c(i, j1) = Vj(1qij) 3.3 which is the difference between the survival value of the target before strike (Vj) and the survival value of the target after strike (Vjqij). Next consider the cost c(i, j2) of the arc (i, j2) which denotes the change in the survival value of target j when weapon i is PAGE 71 60 assigned to it as the second weapon. To determ ine this, we need to know the survival value of target j before weapon i is assigned to it. But this cost depends upon which weapon was assigned to it as the first weapon. The first weapon striking target j can be of any weapon type 1,2, ,m and we do not know its type a priori. Therefore, we cannot determine the cost of the arc (i, j2). However, we can determine an upper bound on the cost of the arc (i, j2). We will next derive the expression for the cost of the arc (i, jk), which as a special case includes (i, j2). Suppose that the first (k1) weapons assigned to target j are of weapon types i1, i2, ..., ik 1, and suppose that the type of the kth assigned weapon is of type i. Then, the survival value of target j after the first (k1) weapons is 121...kjijijijVqqqand the survival value of the target j after k weapons is 121...kjijijijijVqqqq. Hence, the cost of the arc (i, jk) is the difference between the two terms, which is c(i, jk) = 121...(1)kjijijijijVqqqq 3.4 Let qj max = max{ qij: 1, 2, ..., m }.Then, we can obtain an upper bound on c ( i, jk) by replacing each qij by qj max. Hence, if we set c ( i ,jk) = max1()(1)k jjijVqq, 3.5 we get an upper bound on the to tal destruction on assigning weapons to targets. It directly follows from equation 3.5 that c ( i, j1) > c ( i, j2) > . .> c ( i, jk1) > c ( i, jk), 3.6 which implies that the optimal maximum cost flow in the network G will be a contiguous flow. It should be noted here that si nce this is a maxi mization problem, we solve it by first multiplying all arc costs by 1 and then using any minimum cost flow algorithm. Let z* represent the upper bound on destruction cause d to targets after all the PAGE 72 61 assignments obtained by solving this maxi mum cost flow problem. Then, the lower bound on the objective function of formulation 3.1 is 1 n j jV z*. We can summarize the preceding discussion as follows: Theorem 2. If z* is the optimal objective fu nction value of the maximum cost flow problem in the network G, then1 n j jV z* is a lower bound for the weapon target assignment problem. 3.2.3 Maximum Marginal Return Based Lower Bounding Method In this section, we describe a different relaxation that provi des a valid lower bound for the WTA problem. This approach is base d on underestimation of the survival of a target when hit by a weapon as we assume that each target is hit by the best of the weapons. Let qj min be the survival probability for target j when hit by the weapon with the smallest survival probability, i.e., qj min = min{ qij: i = 1, 2, , m }. Replacing the term qij in formulation 3.1 by qj min, we can formulate the WTA problem as follows: Minimize min 1 1()ijm n x jj j iVq 3.7a Subject to 1 n ij j x Wi, for all i = 1, 2, , m, 3.7b ij x 0 and integer for all i = 1, 2, , m and for all j = 1, 2, , n 3.7c 1Let ,m j ij i x x and if we let min()()jx jjjjgxVq then we can rewrite (3.7) as: Minimize 1()n jj jgx 3.8a Subject to PAGE 73 62 1 n ij j x Wi, for all i = 1, 2, , m, 3.8b 1 n ijj i x x for all i = 1, 2, , m, 3.8c ij x 0 and integer for all i = 1, 2, , m and for all j = 1, 2, , n 3.8d It is also possible to eliminate the variables xij entirely. If we let 1 m i iWW then we can rewrite formulation 3.8 as an equivalent integer program 3.9: Minimize 1()n jj jgx 3.9a Subject to 1 n j j x W 3.9b j x 0 and integer for all j = 1, 2, , n 3.9c It is straightforward to transform a soluti on for 3.9 into one for 3.8 since all weapon types are identical in formulation 3.8. Observe that the formulations 3.1 and 3.7 have the same constr aints; hence, they have the same set of solutions. However, in the formulation 3.7, we have replaced each qij by qj min= min{ qij: i =1, 2, , m }, where qj min < qij Noting that the optimal solution value for problems 3.7, 3.8 and 3.9 are all id entical, we get the following result: Theorem 3. An optimal solution value for 3.9 is a lower bound for the WTA problem. Integer program 3.9 is a special case of the knapsack problem in which the separable costs are monotone decreasing and co ncave. As such it can be solved using the greedy algorithm. In the following, assigned( j) is the number of weapons assigned to PAGE 74 63 target j, and value(i, j) is the incremental co st of assigning the next weapon to target j. algorithm combinatoriallowerbounding; begin for j := 1 to n do begin assigned(j) := 0; value(j) := gj(assigned(j)+1) gj(assigned(j)); end for i = 1 to m do begin find j corresponding to the minimum value(j); assigned(j) := assigned(j) + 1; value(j) := gj(assigned(j)+1) gj(assigned(j)); end; end. Figure 3.4 Combinatorial lower bounding algorithm This lower bounding scheme is in fact a variant of a popular algorithm to solve the WTA problem, which is known as the maximum marginal return algorithm In this algorithm, we always assign a weapon with maximum improvement in the objective function value. This algorithm is a heuristic algorithm to solve the WTA problem but is known to give an optimal solution if all weapons are identical. We now analyze the running time of our lower bounding algorith m. If we store value(j) in a Fibonacci heap. We point that after the initiali zation of the heap with the initial values, we perform W findmin opera tions and W decreasekey steps, for a total running time of O( W ) time. In our implementation, we used binary heaps, which run in O( W log W ) time but are comparably fast for th e problem sizes that we considered. PAGE 75 64 3.3 Branch and Bound Algorithm We developed and implemented four branch and bound algorithms based on the four lower bounding schemes described in th e previous section. A branch and bound algorithm is characterized by the branching, lower bounding and sear ch strategies. We now describe these strategies for our approaches. Branching strategy: To keep the memory requirement low, the only information we store at any node is which variable we branch on at that node; and the lower and upper bounds at the node. To r ecover the partial so lution associated with a node of the branch and bound tree, we trace back to the root of the tree. The branching strategy we have used in our implementation is based on the maximummarginal return. For each node of the branch and bound tree, we find th e weapontarget combination which gives best improvement and set the corresponding vari able as the one to be branched on next. Ties are broken arbitrarily. Lower bounding strategy: We used the three lower bounding strategies described in Section 3.2. We provide a comparative analysis of these bounding schemes in Section 3.5. Search strategy: We implemented both the breadthf irst and depthfirst search strategies. We found that fo r smaller size problems (i.e., up to 10 weapons and 10 targets), breadthfirst strategy gave overall better results; but for larg er problems, depthfirst search had a superior performance. We re port the results for the depthfirst search in section 3.5. PAGE 76 65 3.4 A Very LargeScale Neighborhood Search Algorithm In the previous two sections, we desc ribed branch and bound algorithms for the WTA problem. These algorithms are the first exact algorithms that can solve moderate size instances of the WTA problem in reasonable time. Nevertheless, there is still a need for heuristic algorithms, which can solve larg escale instances of the WTA problems. In this section, we describe a neighborhood s earch algorithm for the WTA problem, which has exhibited excellent computational results. This algorithm is an application of very largescale neighborhood (VLS N) search to the WTA problem. A VLSN search algorithm is a neighborhood search algorithm where the size of the neighborhood is very large and we use some implicit enumeration algorithm to identify an improved neighbor. We refer the reader to the paper by Ahuja et al. [2002a] for an overv iew of VLSN search algorithms. A neighborhood search algorithm starts with a feasible solution of the optimization problem and successively improves it by replacing it by an improved neighbor until it obtains a locally optimal solution. The quality of the locally optimal solution depends both upon the quality of the st arting feasible solution and the structure of the neighborhood, that is how we define the neighborhood of a given solution. We next describe the method we used to construct the starting feasible solution followed by our neighborhood structure. 3.4.1 A Minimum Cost Flow formulatio n based Construction Heuristic We developed a construction heuristic, wh ich solves a sequence of minimum cost flow problems to obtain an excellent soluti on of the WTA problem. This heuristic uses the minimum cost flow formulation shown in Figure 3.3, which we used to determine a lower bound on the optimal solution of the WTA problem. Recall that in this formulation, PAGE 77 66 we define the arc costs ( i j1), ( i j2), , ( i jk), which, respectively, denote the cost of assigning the fi rst, second and kth weapon of type i to target j Also recall that only the cost of the arc ( i j1) was computed correctly, and for the other arcs, we used a lower bound on the cost. We call the arcs whos e costs are computed correctly as exactcost arcs and the rest of the arcs as approximatecost arcs This heuristic works as follows. We fi rst solve the minimum cost flow problem with respect to the arc costs as defined earli er. In the optimal solu tion of this problem, exactcost arcs as well as approximatecost arcs may carry positive flow. We next fix the part of the weapon target assi gnment corresponding to the flow on the exactcost arcs and remove those arcs from the network. In othe r words, we construct a partial solution for weapontarget assignment by as signing weapons only for exact cost arcs. After fixing this partial assignment, we again compute th e cost of each arc. Some of the earlier approximatecost arcs will now become exactcost arcs. For example, if we set the flow on arc ( i j1) equal to 1, we know that that weapon i is the first wea pon striking target j and hence we need to update the costs of the arcs ( l jk) for all l = 1, 2, , m and for all k 2. Also observe that the arcs ( l j2) for all l = 1, 2, , m now become exact cost arcs. We next solve another minimum cost flow problem and again fix the flow on the exactcost arcs. We recompute arc costs, make so me additional arcs exactcost, and solve another minimum flow problem. We repeat th is process until all weapons are assigned to the targets. We tried another modification in the mini mum cost flow formulation, which gave better computational results. The formulati on we described determines the costs of approximatecost arcs assuming that the wo rst weapons (with th e largest survival PAGE 78 67 probabilities) are assigned to targets. However, we observe d that in any nearoptimal solution, the best weapons are assigned to th e targets. Keeping this observation in mind, we determine the costs of valid arcs assumi ng that the best weapons (with the smallest survival probabilities) are assigned to targets. Hence, the cost of the arc (i, jk) which is c ( i ,jk) = 121...(1)kjijijijijVqqqq is approximated by c(i, jk) = Vj [qmin(j)]k1(1qij) Our experimental investigation shows that this formulation generates better solutions compared to the previous formulation. We present computationa l results of this formulation in Section 3.5. 3.4.2 The VLSN Neighborhood Structure The WTA problem can be conceived of as a partition problem defined as follows. Let S = { a1, a2, a3, . an} be a set of n elements. The partition problem is to partition the set S into the subsets S1, S2, S3, , SK such that the cost of the partition is minimum, where the cost of the partition is the sum of the cost of each part. The WTA problem is a special case of the partition problem where the set of all weapons is partitioned into n subsets S1, S2, , Sn, and subset j is assigned to target j 1 j n Thompson and Orlin [1989] and Thompson and Psaraftis [1993] proposed a VLSN search approach for partitioning problems which proceeds by performing cyclic exchanges Ahuja et al. [2001a, 2003a] proposed further refinements of this approach and applied it to the capacitated minimum spanning tree problem. We will present a brief overview of this approach when applied to the WTA problem. Let S = ( S1, S2, , Sn) denote a feasible solution of the WTA problem where the subset Sj, 1 j n denotes the set of wea pons assigned to target j Our neighborhood search algorithm defines neighbors of the solution S as those solutions that can be PAGE 79 68 obtained from S by performing multiexchanges A cyclic multiexchange is defined by a sequence of weapons i1i2i3 iri1 where the weapons i1, i2, i3, , ir belong to different subsets Sjs. Let t(i1) t(i2) t(i3) , t(ir) respectively, denote the targets to which weapons i1, i2, i3, , ir, are assigned. The cyclic multiexchange i1i2i3 iri1 represents that weapon i1 is reassigned from target t ( i1) to target t ( i2), weapon i2 is reassigned from target t ( i2) to target t ( i3), and so on, and finally weapon tr is reassigned from target t ( ir) to target t ( i1). We can similarly define a path multiexchange by a sequence of weapons i1i2i3 ir which differs from the cyclic multiexchange in the sense that the last weapon ir is not reassigned and re mains assigned to target t ( ir). The number of neighbors in the multiex change neighborhood is too large to be enumerated explicitly. However, using the concept of improvement graph a profitable multiexchange can be identified using network algorithms. The improvement graph G (S) for a given feasible solution S of the WTA problem contains a node r corresponding to each weapon r and contains an arc ( r l ) between every pair of nodes r and l with t ( r ) t ( l ). The arc ( r l ) signifies the fact that weapon r is reassigned to target (say j) to which weapon l is currently assigned and weapon l is unassigned from its current target; the cost of this arc, crl, is set equal to the change in the survival value of the target. Let j V denote the survival value of the target j in the current solution. Then, the cost of the arc ( r, l ) is crl = (/)1jrjljVqq We say that a directed cycle W = i1i2i3 iki1 in G (S) is subsetdisjoint if each of the weapons i1, i2, i3, , ik is assigned to a different target. Thompson and Orlin [1989] show ed the following result: Lemma 1. There is a onetoone correspondence between multiexchanges with respect to S and directed subsetdis joint cycles in G(S) and both have the same cost. PAGE 80 69 This lemma allows us to solve the WTA problem using the following neighborhood search algorithm: algorithm WTAVLSN search; begin obtain a feasible solution S of the WTA problem; construct the improvement graph G(S); while G(S) contains a negative cost subsetdisjoint cycle do begin obtain a negative cost subsetdisjoint cycle W in G(S); perform the multiexchange corresponding to W; update S and G(S); end; end; Figure 3.5 The VLSN search algorithm for the WTA problem. We now give some details of the VLSN search algorithm. We obtain the starting feasible solution S by using the minimum cost flow base d heuristic described in Section 3.4.1. The improvement graph G (S) contains W nodes and O( W2) arcs and the cost of all arcs can be computed in O( W2) time. We use a dynamic pr ogramming based algorithm (as described by Ahuja et al. [2003]) to obtain subsetdisjoint cycles This algorithm first looks for profitable twoexchanges involving two targets only; if no profitable twoexchange is found, it looks for profitable th reeexchanges involving th ree targets; and so on. The algorithm either finds a profitable multiexchange or terminates when it is unable to find a multiexchange involving k targets (we set k = 8). In the former case, we improve the current solution, and in the latte r case we declare the cu rrent solution to be locally optimal and stop. The running time of the dynamic programming algorithm is O( W2 2k) per iteration, and is typical ly much faster since most cyclic exchanges found by the algorithm are swaps. PAGE 81 70 3.5 Computational Results We implemented each of the algorithm described in the previous section and extensively tested them. We tested our al gorithms on randomly generated instances as data for the reallife instances is classifi ed. We generated the data in the following manner. We generated the target survival values Vjs as uniformly distributed random numbers in the range 25100. We generated the kill probabilities for weapons engaging with the targets as uniforml y distributed random numbers in the range 0.600.90. We performed all our tests on a 2.8 GHz Pentiu m 4 processor computer with 1 GB RAM PC. In this section, we present the results of these investigations. 3.5.1 Comparison of the Lower Bounding Schemes In our first investigation, we compared the tightness of the lower bounds generated by the lower bounding algorithms developed by us. We tested the three lower bounding schemes described in Section 3.2: (i) LP based lower bounding scheme; (ii) MIP based lower bounding scheme; (iii) the minimum co st based lower bounding scheme; and (iv) the maximum marginal return based lower bounding scheme. We tested an additional lower bounding scheme which is a va riant of the LP based scheme. Table 3.1 gives the computational results of these four lower bounding schemes. For each of these schemes, the first column gives the % gap from the optimal objective function value and the second column give s the time taken to obtain bound. The following observations can be derived from this table: (i) the MIP lower bounding scheme gives the tightest lower bounds but al so takes the maximum computational time; (ii) the minimum cost flow based bounding scheme gives fairly tight lower bounds when the number of weapons is less than or equal to the number of targets; and (iii) the PAGE 82 71 maximum marginal return algorithm takes th e least amount of time to obtain lower bounds. Table 3.1 Comparison of four lower bounding schemes. LP Scheme MIP Scheme Min Cost Flow Scheme Max Marginal Return Scheme # of Weapons # of Targets % Gap Time (in secs) % Gap Time (in secs) % Gap Time (in secs) % Gap Time (in secs) 5 5 8.03 0.015 0.21 0.016 1.66 <0.001 10.61 <0.001 10 10 3.63 0.015 0.12 0.031 0.00 <0.001 11.01 <0.001 10 20 19.70 0.015 0.04 0.062 0.00 <0.001 1.45 <0.001 20 10 11.88 0.015 0.53 0.156 21.32 <0.001 19.00 <0.001 20 20 7.28 0.031 0.25 0.109 1.32 <0.001 6.40 <0.001 20 40 23.35 0.046 0.04 0.296 0.00 <0.001 1.57 <0.001 40 10 14.79 0.015 2.12 0.609 42.41 <0.001 46.89 <0.001 40 20 7.06 0.031 0.45 0.359 25.52 0.015 13.53 <0.001 40 40 6.83 0.078 0.11 0.703 1.63 0.015 3.05 <0.001 40 80 21.69 0.14 0.03 1.812 0.00 0.046 0.88 <0.001 3.5.2 Comparison of Branch and Bound Algorithms We developed branch and bound algorith ms using the preceding lower bounding schemes. Table 3.2 gives the results of thes e algorithms. The bran ch and bound algorithm using the LP based lower bounding scheme di d not perform well at all and we do not present its results. We replaced this algor ithm by another algorithm, which we call the hybrid algorithm The hybrid algorithm computes lo wer bounds using both the minimum cost flow based and the maximum margin al return based lower bounding schemes and uses the better of these two bounds. We fi nd that the branch and bound algorithm using the MIP based lower bounding gives the most consistent results and able to solve the largest size problems (containing 80 weapons and 80 targets). We also find that the hybrid algorithm also gives excellent result s for those instances where the number of weapons is less than or equal to the number of targets. PAGE 83 72 Table 3.2 Comparison of branch and bound algorithms. MIP Based B&B Algorithm Min Cost Flow Based B&B Algorithm Maximum Marginal Return Based B&B Algorithm Hybrid Algorithm # of Weapon s # of Targets Nodes Visited Time (in secs) Nodes Visited Time (in secs) Nodes Visited Time (in secs) Nodes Visited Time (in secs) 5 5 15 0.14 11 <0.001 23 <0.001 11 <0.001 10 10 29 0.56 1 <0.001 181 <0.001 1 <0.001 10 20 23 0.83 1 <0.001 83 <0.001 1 0.015 20 10 101 7.27 2,8611 1.34 20,251 2.52 20 20 109 6.56 2,383 4.39 15936 0.94 1,705 2.50 20 40 105 16.58 1 <0.001 111,603 10.14 1 0.015 40 10 1,285 327.27 40 20 205 35.19 ~108 13,651.9 ~107 25,868.9 40 40 211 50.96 ~106 10,583.62~106 943.03 38,3275 1,891.83 40 80 385 235.41 1 0.031 1 0.031 80 40 117,227 43,079.5580 80 44905 58,477.3180 160 1055 3,670.49 1 0.062 1 0.062 3.5.3 Performance of the VLSN Search Algorithm We now present computational results of the minimum cost flow based construction heuristic and the VLSN search algorithm. Table 3.3 gives the objective function values of the solutions obtained by the construction heuris tic and the improved values when VLSN search algorithm is app lied to these solutions. We observe that the construction heuristic obtained optimal solu tions for over 50% of the instances and for the remaining instances the VLSN search algorithm converted them into optimal or almost optimal solutions. The computationa l times taken by these algorithms are also very small and even fairly large in stances are solved within 3 seconds. PAGE 84 73 Table 3.3 Results of the construc tion heuristic and the VLSN search algorithm. Construction Heuristic VLSN Algorithm # of Weapons # of Targets Optimality Gap Time (in seconds) Optimality Gap Time (in seconds) 10 5 0% <0.001 0% <0.001 10 10 0% <0.001 0% <0.001 10 20 0% <0.001 0% <0.001 20 10 0% <0.001 0% <0.001 20 20 0% <0.001 0% <0.001 20 40 0% 0.015 0% 0.015 20 80 0% 0.015 0% 0.031 40 10 1.79% 0.015 0% 0.031 40 20 0.33% 0.015 0% 0.015 40 40 0% 0.015 0% 0.015 40 80 0% 0.031 0% 0.078 40 120 0% 0.062 0% 0.109 80 20 2.33% 0.109 0% 0.156 80 40 0.10% 0.062 0% 0.109 80 80 0.0003% 0.093 0.0003% 0.156 80 160 0% 0.172 0% 0.219 80 320 0% 0.390 0% 0.625 100 50 0.79% 0.120 0.0015% 0.437 100 100 0.001% 0.187 0.0009% 0.250 100 200 0% 0.375 0% 0.609 200 100 0.01% 0.656 0.0059% 0.828 200 200 0.001% 0.921 0.0008% 1.109 200 400 0% 1.953 0% 2.516 3.6 Conclusions In this chapter, we consider the w eapon target assignment problem, which is considered to be one of th e classical operations research problems that has been extensively studied in the literature but still has remained unsolved. Indeed, this problem is considered to be the holy grail of de fenserelated operations research. Though weapon target assignment problem is a nonlinear integerprogra mming problem, we use its special structure to develop LP, MIP, ne twork flow, and combinatorial lower bounding PAGE 85 74 schemes. Using these lower bounding schemes in branch and bound algorithms gives us effective exact algorithms to solve the WT A problem. Our VLSN s earch algorithm also gives highly impressive results and gives either optimal or almost optimal solutions for all instances it is applied to. To summarize, we can now state that the WTA problem is a wellsolved problem and its largescale inst ances can also be solved in realtime. PAGE 86 75 CHAPTER 4 UNIVERSITY COURSE TIMETABLING PROBLEM 4.1 Introduction Timetabling is an important problem en countered by almost everybody either in terms of using it or making it. One of the major categories of tim etabling which is of importance in an academia is university or school timetabling. In practice, the school or university timetable is prepar ed manually or with the help of simulation software or customized optimization software. Researcher s in OR field have proposed number of algorithms for timetabling problem to mi nimize the conflict and to maximize the preferences in timetable. A university or schooltimetabling problem, in broad sense, can be defined as the problem of finding a sequence of lectures or examinations be tween teachers and students in a prefixed period of time satisfying a set of constraints of various types. We can categorize timetabling problem in three pa rts based on use and context, namely, Course timetabling (in university): Scheduling of courses, i.e., fixing the time and place for lectures, with the objective of optimizing two important goals, (i) minimizing students interest c onflict and (ii) maximizing teachers preference. Teachers preference is important as he/she has got dua l responsibility of teaching and research and students interest is important as there are lot of freedom to stude nts in selecting the course. Course timetabling is generally done for a week, which is repeated through out PAGE 87 76 semester or trimester. The problem may also include assignment of courses to classrooms so that distance covered by students in between the classes can be minimized. Examination timetabling: Scheduling of examinations over a prespecified period of time with the objective of eliminating the conflict in students schedule (i.e. a student cannot be scheduled to take two exams at the same time) over a set of constraints type. Here constraints are mainly in terms of minimum gap between two exams of a student, different requirements for compulsory and elective courses and rooms availability. School timetabling: Scheduling of courses, i.e., fixing the place and time for different courses (lectures) w ith the objective of minimizing t eacher clash, i.e., a teacher cannot be asked to take two lectures at the same time. Although, it looks similar to university course timetabling, the basic differe nce lies in terms of teachers assignment and clash of students interest. The number of lectures taken by a teacher in a school is normally much more than that in university whereas all students in a class or group follow same schedule and assigned to same room as against there are much more electives and students have to roam in campus to attend different classes. One can easily notice from above that objectives for timetabling problem are nothing but set of constraints, e.g., we can make the constrai nts that a teacher should be assigned to teach in most preferred periods and there should not be any conflict in students interest. Usually, timetabling problems are defined with constraints only, which are further divided into two categories soft and hard. Hard constraints are those, which must be honored by a timetable whereas soft constraints may be violated, which are subsequently made part of the objective with a cost, associated with every violation of soft constraints. For example, in our case we are considering teachers preference and PAGE 88 77 students interest conflict as soft constrai nts and make it part of objective functions. An ideal timetable should have zero objective f unction value, i.e., all the constraints are strictly followed. Here in this chapter, we are restricti ng our focus to course timetabling in a university and more specifically to the fixing of time for th e courses for a week. We are not considering here the assignment of c ourses to the room, which can be taken separately as a subproblem with the obj ective of minimizing the distance traveled between classrooms by a student. Our approa ch is based on application of a metaheuristic Very Large Scale Neighborhood S earch, which has shown very encouraging results for number of combinatorial optimization problem. We have structured this chapter in sectio ns with section 4.2 dealing with problem definition, section 4.3 deal ing with literature survey section 4.4 dealing with mathematical formulation, section 4.5 deali ng with random instance generation approach, section 4.6 dealing with heuristics for initial solution, section 4.7 de aling with very large scale neighborhood search, sectio n 4.8 dealing with VLSN a nd tabu search, section 4.9 dealing with implementation and finally s ection 4.10 dealing with conclusions drawn from the above chapter. 4.2 Problem Definition As against the available literatures, which are generally specialized to the case specific algorithms for course timetabling prob lem, in this chapter, our attempt is to develop an algorithm applicable to genera l cases. The terminology for the course timetabling as used in this chapter are follows: PAGE 89 78 Period: An interval of time in which a teacher gives lecture to a group of students in a single sitting. Different periods can ha ve different duration and are characterized by starting time. For this chapter we have numbe red the periods in increasing order of the days and time for whole week. We have also assumed that number of periods is same everyday, i.e., if there are 5 periods in a day then there will be 25 periods in week. Slot: Slot is the set of periods, to which a cour se can be assigned, i.e., if a course is assigned to slot s and s consist of periods p1 p2 and p3 then this course will be taught in periods p1 p2 and p3 every week. We assume that we know allpossible slots i.e., all possible combinations of periods apriory. Course: A course is the material, which is taught. We have not categorized the courses here; however, generally this belongs to either com pulsory or elective category. For a course we know the teacher and how ma ny times in a week the course is to be taught (i.e., credit hours ). We also know apriori number of students goi ng to take the course, which is generally derived from past data. Permissible slots: Set of slots to which a course can be assigned. The number of periods in a permissible slot must be equal to credit hours for cour se it belongs to. We assume that the permissible slots are given apriori for every course. Now, we can define the coursetimetabli ng problem for this chapter as: Minimize the sum of the penalty of assigning a teacher (i .e., a courses the teacher is taking) to an undesired slot and penalty for assigning two courses, in both of which a student is interested, to common period; subject to following constraints. Each course must be assigned to a sl ot in the set of permissible slots. A teacher cannot teach more th an one course in a period. The number of courses to be taught in a period must not exceed the available number of rooms in that period. PAGE 90 79 All preassignments of course to a slot and unavailability of teacher/room in a slot (period) must be honored. There may be additional constraints like co mpactness of timetable (i.e., minimum number of periods required), consideration of multiple sections of a course etc., which are not considered here in this chapter. 4.3 Literature Review For the last three decades researchers have been working on this problem and have produced number of algorithms for the course timetabling problem. The survey of the automated timetabling problem by Schaerf [1999] has very systematically categorized the timetabling problems and has given the math ematical formulation for each one along with a general solution appro ach. Schaerf has discussed both exact approach as well as heuristics for the problem. Metaheuristics such as genetic algorithms, tabu search, simulated annealing etc. are successfully applied for timetabling problem. The reduction of timetabling problem to graph coloring probl em and application of heuristics for graph coloring problem has been proved very effici ent. Another survey by Carter and Laporte [1998] gives excellent reference to diffe rent timetabling problems and algorithms available for the same. The paper by Werra [1997] has represen ted the problem as a graph coloring problem and has used generic approach to solve the problem. The approach by Stallaert [1997] has used Integer Programming Formula tion to solve the tim etabling problem for the Anderson School of Management at UC LA. Author has successfully reduced the number of variables by exploiting the featur es of the problem. Thompson and Downsland [1996] have used a variant of simulated a nnealing to solve examination timetabling. PAGE 91 80 Their approach dynamically changes the neighborhood structure and also includes infeasible solution in neighborhood, however with penalty. The paper by Dimopoulou and Miliotis [2001] have solved the course timetabling and examinationtimetabling problem for a specific case and have form ulated the problem as integer program. The work in this chapter is mainly influenced by the application of tabu search for large scale timetabling problem by Hertz [1991] and tabu sear ch approach by Valdes et al. [2002]. Hertz [1991] has also considered the problem of grouping of students in different sections. Authors approach include s three phases, (a) getting initial feasible solution, (b) applying tabu sear ch on timetabling problem and (c) applying tabu search on grouping problem. This paper al so explores infeasible solutions in neighborhood structure, however, keeps a se parate account of it. Valdes et al. [ 2002] has solved both course timetabling as well as room assignment for University of Valencia, Spain with the application of tabu search. Distin ction of authors approach lie s in considering not only of pairwise exchanges in neighborhood search but that of multiexchanges as well. Authors also discuss briefly the approach for assigning courses to rooms. In this chapter we are applying very large scale neighborhood (VLSN) technique, which has shown very good results for capaci tated telecommunicatio n network design (Ahuja et al. [2001a]) as well as on other combinatorial optimization problems. VLSN uses very efficiently the concept of im provement graph to find an improving cycle consists of multiple moves. 4.4 Mathematical Formulation The coursetimetabling problem as defined above can be formulated as quadratic integer program as given below. PAGE 92 81 Parameters C : set of courses CTt: set of courses for teacher t T : set of teachers S : set of slots D: set of periods di: set of slots which contains period i R : set of room types CRr: courses requiring room of type r mrs: number of rooms of type r available in slot s Si: set of permissible slots for course i PSt: set of preferred slots for teacher t ij: cost of assigning teacher (course) i to slot j this may be made zero for the preferred slot for a teacher qij: number of students taking both the courses i and j ij: penalty if courses i and j are assigned to same slots, this may be made very high for the cases when courses i and j are compulsory for a class, i.e., courses i and j cannot be taught in same slot ij: number of periods common in slots i and j ij: 1 if ij > 0, 0 otherwise F : set of preassigned courses pi: slot to which course i is preassigned CSp: set of courses taken by student p PAGE 93 82 P : set of students Decision Variable: xis: 1, if course i is assigned to slot s Si 0, otherwise Objective: min is C iS s isxi + 1212 12 ij s sijijisjs iCjisSsS jCqxx 4.1a subject to: Each course i must be assigned to a slot s iS s isx= 1, i C 4.1b A teacher cannot be assigned to mo re than one course at a time '' 'tjisssjs jCTsS ji x x 1, i CTt, t T, s Si, 4.1c Room availability in each slot ri pis iCRsS sd x mrs, p D, r R 4.1d Preassignment of courses to slots iipx = 1, i F 4.1e Integrality constraint 1 0isx s Si, i C 4.1f Although literature suggests many appro aches to solve timetabling as integer program, we assume here that when considered in general way, it may not be possible to reduce the number of variables, which is ve ry crucial for the successful implementation PAGE 94 83 of this approach. In our opinion nonlinear terms in objective function and binary variables make the above formulation difficult. We implemented the integerprogramming equi valent with linear objective function value for above formulati on by introducing variable yj for each quadratic term in objective function value and a dding following constraints. 1is x + 2js x ky 1, for each quadratic term in objective function value for which 12 s sijijq > 0 4.2a 0 ky 1 4.2b The experimental results on the test inst ances show that convergence to optimal solution for bigger size inst ances are very difficult. 4.5 Random Instance Generation In literature, applications of algorithm s are specific to the problem set up of a university or school and it is very difficult to obtain problem instances to compare the results. Therefore, we have developed here one approach for random instance generation. For this chapter, we have made followi ng assumptions in problem defined above and before generating instance. The durations of all the peri ods are same; hence a course can be assigned to a slot having number of periods in it equal to cr edit hours. This can be generalized by making the condition, as total duration of th e periods in a slot should be equal to credit hours for the course. All rooms are having same capacity and hen ce a constraint for room availability is limited to total number of rooms available in a period rather than number of rooms available of different capacity. The constraint for preassignment of cour ses and unavailability of room is not considered. PAGE 95 84 In our approach, we take three inputs for instance generation, namely, total number of periods, total number of slots and total number of courses. The number of periods in a week should be multiple of 5 as we have assu med that there are 5 working days in a week and number of periods per day is same. The f easibility of the instan ce highly depends on the number of courses, which user has to sp ecify judiciously. With these inputs, the algorithm for instance generation works as follows. 1. The number of periods in each slot takes value 2, 3 or 4 randomly. The options are fixed as we have assumed that credit hour s for a course can be 2, 3 or 4 only. 2. Assignment of periods to slot s has been done as follows: If number of periods is 2 then first pe riod can be anyone but last period of the week, and then, second period can be either in continuation (on same day) with first period with 50 % probability, otherwise must be from other day (days after the day of first period). If number of periods in a slot is 3 th en first period can be anyone from the set of periods but those of last day. Second period can be in continuation with first period with 30 % probability. Similarly, if first and second periods are not in continuation than with 30 % probability second and third periods can be in continuation, otherwise be long to different days. If number of periods are 4, then it consists of periods from two days. First two periods are continuous on one day and last two periods are continuous on another day. 3. Each course is assigned cred it hours 2, 3 or 4 randomly. 4. Teacher is assigned to the course ra ndomly. The number of courses taught by a teacher is kept random in the range 1 to 2. 5. The number of permissible slots for each cour se is kept in the range of 25% to 75% of the total number of slot s having number of periods e qual to credit hour. The set of permissible slot for a course is gene rated randomly. It can be noted that two courses may have same set of permissible slots. 6. We keep number of rooms equal to 1. 5 (total credit hour s/total number of periods). In ideal case the number of rooms should be equal to total credit hours/total number of periods; however, we keep a slack of half of required rooms to take care of diffe rent constraints. 7. The penalty for scheduling two courses in sa me period when a stude nt is interested in both the courses is randomly generated in range (5, 20). PAGE 96 85 8. The percentage of pairs of courses, whic h may have common students, is fixed at 10 % (it can be varied by user). We genera te such pair of courses randomly. The number of common students in these pair of courses is generated randomly in the range of (0, 40%) of minimum number of students in a cour se between the two courses. 9. The penalty for assigning a teacher to a slot is kept in the range (10, 50). In this algorithm number of parameters can be varied ba sed on the specific application. We expect that this procedure of instance generation represent a wide variety of course timetabling setup. 4.6 Building Initial Solution Getting a feasible solution in case of course timetabling is not a trivial task as in case of many combinatorial optimization probl ems. We have developed two construction heuristics to obtain initial solution. For course timetabling, construc tion heuristic has to meet two requirements namely, it should be able to produce feasible solution, and solution should be of good quality. In both the approaches courses are assigned to slots one by one based on greedy criteria to impr ove the quality of the solution. One case decision of assigning course to a slot is solely based on minimum increase in objective function value, whereas in another approach we consider the gap in increase in objective function value if assign ed to best available slot or s econd best availabl e slot. In every procedure we check, at the time of every assign ment, the availability of room in a period and conflict of teacher in terms that a teach er cannot teach two courses at same time. 4.6.1 Greedy Heuristic 1 The procedure for obtaining initial solution, in the case when we consider only the increase in objective function value at each stage can be described as follows. PAGE 97 86 Sort the courses in nondecreasing order of number of permissible slots. This ordering gives priority to courses that have less flexibility in assigning to the slots. Take courses in above order and calculat e the incremental cost in assigning this course to its permissible slots, without violating room availability and teachers conflict constraint. If found such slot(s), assign the course to slot with least incr emental cost. Repeat Step 2 for all the courses. If a course (c1) could not get assigned in Step 2 and Step 3, find the slot (s1) for this course in its set of perm issible slots by shifting course (c2) already assigned to s1 to some other feasible slot (s2). We explore all such possibilities and assign c1 to s1 and c2 to s2 with least incremental cost. The incremental cost here is change in cost by assigning c1 to s1 (considering c2 has been removed from s1) plus change in cost by assigning c2 to s2 minus change in cost by removing c2 from s1. This procedure does not always guarantee th e generation of feasib le solution. While applying this heuristic in conj ecture with VLSN, we need different initial solutions in different runs for an inst ance. We obtain this by assi gning a course to one among p (a parameter) best slots randomly instead of that to the best slot. 4.6.2 Greedy Heuristic 2 As applied in greedy heuristi c 1, the selection of a cour se for next assignment only based on number of available sl ots may lead to a solution, which is not good in terms of quality. We propose here another approach, which has shown very good result (Kiaer et al. [1992]). Let the set S contains the courses, wh ich are yet to be assigned to a slot. At start it contains all the courses. Let set Ai contains the available slots to which course i can be assigned at any stage. It implies, in periods of these slots, there is no conflict of teacher and room is available to t each the course. Initially, set Ai contains all the permissible slots for course i. PAGE 98 87 Algorithm: 1. Create vector cost and fdeg for each course in S, where cost(i, j) represents the change in objective function value if course i S is assigned to the slot j Ai and fdeg(i) is the cardinality of Ai. Initially, cost(i, j) is set to the cost if teacher of course i is assigned to slot j and fdeg(j) is set to the number of permissible slots for course i. 2. For each course in S, if fdeg(i) = 0, mark the course as unassigned and remove it from set S else if fdeg(i) = 1, assign the course i to only slot available. Remove i from set S and update the cost vector for all other courses in S. 3. If fdeg(i) > 1 for each i S, calculate del(i), the difference in two lowest c(i, j) for each course and select the course k for next assignment for which del(i) is maximum over all i S, use lower fdeg(i) in case of tie. Assign the slot with minimum cost(k, j) to the course k. Remove the course k from S and update the vectors cost and fdeg. 4. Repeat Step 2 and 3 until S = Figure 4.1 Algorithm for assignment of courses using Greedy Heuristic 2 Above algorithm in figure 4.1 assumes that if we will not assign a course with maximum del(i) in current stage, then it may result in assigning this course to the slot having high increase in objectiv e function value. For all the courses, which could not be assigned to any slot, we perfor m backtracking and try to allocate slot to these courses. To randomize the initial solution, we will select randomly one among best 2 or 3 courses, as got in Step 3. 4.7 Very Large Scale Neighborhood Search (VLSN) We use the partitioning feature of the timetabling problem, which is expected to give a very efficient VLSN approach. In cour se timetabling problem we can think of slots as different sets and courses assigned to slots as elements Then the problem becomes the partitioning of all the courses into available slots. Our proposed heuristic for the coursetimetabling problem is based on the fo llowing structure of neighborhood search. PAGE 99 88 Cyclic exchange: We call a cycle c(1) c(2) c(3) . c(k) c(1) consisting of courses c(1),c(2),c(3), ., c(k) such that course c(1) is assigned to slot of course c(2), course c(2) is assigned to slot of course c(3) and so on and finally c(k) to slot of c(1), as Cyclic Exchange. If the change in objective function value is negative due to such exchange, we call it Negative Cycle. Please also note that this is a subsetdisjoint cycle as elements in cycle are assigned to different slots. However, periods in th ese slots may not be disjoint. Path exchange: We call a set of courses c(1),c(2),c(3), ., c(k) such that course c(1) can be assigned to slot of course c(2), course c(2) can be assigned to slot of course c(3) and so on, as path exchange. It can be noted that in path exchange, cardinality of slot of c(k) increases by one and that of c(1) decreases by one. Neighborhood function N (x): Set of all feasible solutions that can be reached from x (current solution) by making cyclic exchange with maximum cycle length k. Improvement graph G (N, A): This is a graph consists of nodes N, a set of all courses; arcs A, a set of all arcs (i, j) such that, (i) course i and course j are currently assigned to different slots, (ii) assignment of course i to the slot of course j should not violate any constraint (considering course j has been removed from its current slot) and, (iii) the slot of course j is in the set of permissible slots for course i. One can think of arc (i, j) as representing that the course i is assigned to curr ent slot of course j and course j will leave its current slot. We also add dummy nodes and dummy arcs in the Improvement Graph, to consider the possibility of path exchange. The details of this are omitted here, however, one can refer to Ahuja et al.[2001a] for details. Cost (i, j): This represents the cost of arc (i, j) in improvement graph. This can be visualized as the change in objective function if course i is assigned to slot of course j and PAGE 100 89 course j leaves its current slot. It can be noteed that cost (i, j) depends only on the courses assigned to periods in the current slots of courses i and j. Above formulation of Improvement Graph and Cyclic exchange is based on the following theorem by Thompson and Orlin [1989]. There is one to one corresponding betw een cyclic exchange with respect to partition S (of courses in slots) and subs et disjoint directed cycles in the improvement graph G, and both have the same cost. With the above structure, we propose fo llowing heuristic for course timetabling problem. Start with an initial solution generated with randomized greedy method. Construct the Improvement Graph with resp ect to above initial solution. If the initial solution is infeasible, i.e., some courses are yet to be assigned to a slot, exclude the unassigned courses from Improvement Graph. It implies that we will look for negative cycle consisting of course s, which are already assigned to slots. Find subset disjoint negative cycle in above improvement graph of maximum length k (a parameter). We apply dynamic programming approach to find such cycle, which enumerate the paths in stages of different length, i.e ., it starts with the paths of length 2 in improvement graph, once enumerated all, goes to paths of length 3 and so on. At every stage we check the cycle cost by making cycle by adding arc from last node to first node in path. While calculating path cost we require adjust ment in it to take care of moves of nodes in the path. Suppose there is a path i>j>k, where courses i and j have common students. Also suppose that there is a common pe riod in current slots of courses j and k, which is not there for the current slots of courses i and j. Therefore, students conflict component of cost (i, j) must be zero as there is no common period in current slots of i and j. However, when i will move to the slot of j and course j will move to that of k, then this component will also get value, which is required to be added in path cost at the time of path enumeration. Similarly, conflict in teachers assignment is also required to be checked. If found any negative cost cycle, implem ent the cycle and update the improvement graph with respect to change in timetable due to implementation of negative cycle. If there are unassigned courses (i.e., curren t solution is infeasible), check whether an unassigned course can be assigned to a permissible slot. If so assign the course to a permissible slot resulting in least in crease in objective function value. Repeat Step 3. PAGE 101 90 Else return the current solution as local optimal solution and stop. 4.8 VLSN with Tabu Search Our intuition says that due to hard cons traints in course timetabling problem, the feasible region is highly scattered and sparse Therefore, even with very largescale neighborhood, the local search can converg e to a solution far from global optimal solution. To overcome this we also apply tabu search (Gover et al. [1997]), a memory based popular metaheuristic, which can come out a local optimal point. We have used shortterm memory to avoid cycling and use longterm memory to diversify the search. Our implementation of tabu search for course timetabling problem can be described as follows. Short term memory: In tabu search to minimize the cycling of search, shortterm memory is kept for the moves made in recent past. For a cycle implemented, we prohibit the algorithm to consider same move for certain number of iterations, called tabu tenure. For example, if course i is assigned to slot s in current iteration, then for a prespecified number of ite rations, we will not assign course i to any other slot. We call a move tabu if it results in assigning course i to a slot other than s for prespecified number of iterations. We call a cyclic exchange tabu, if any of the moves in cycle is tabu. We assign thre e different tabu tenures for different tabu moves, which is based on th e desirability of the move. 1. Long tabu tenure for an assignment, wh ich results into most preferred/ second most preferred slot for the teach er. We keep it double of short tabu tenure. 2. Medium tabu tenure if not the abov e case and move do not induce new conflict in students preference. We k eep it 1.5 times the short tabu tenure. 3. Short tabu tenure if assignment results into new conflict in students preference. We keep it 0.1 times the number of courses. Aspiration Criteria: In tabu search, if we find a cyclic exchange which produce lucrative solution, we accept it irrespectiv e of tabu status For course timetabling, we have kept 2 such aspiration criteria as indicated below. 1. If objective function value after implementing the multiexchange is better than best objective function value achieved so far, accept the multiexchange irrespective of its tabu status. PAGE 102 91 2. If a multiexchange results into reduced number of students clash than that achieved so far, accept the multiexcha nge irrespective of its tabu status. Long Term Memory: In tabu search, along with searching solution on a prespecified path, we may need a diversificatio n strategy, so that s earch can be done in wider solution space. To achieve, we keep longterm memory, which keeps the count of number of iterations, in which a cour se is assigned to a slot. This type of longterm memory is called residence frequency. In our search procedure, if there is no improvement in best solution found in a run for certain number of iterations, we restart the tabu search with new initial solution. In algorithm for getting initial solution, we penalize the assignments, wh ich are having high residence frequency. Stopping criteria: We stop the tabu search after 10 number of courses iterations. 4.9 Implementation We have implemented the procedures fo r instance generation, initial solution generation and VLSN search al gorithm. In VLSN, we have implemented two variants of neighborhood structure. In one we do not make arc (i, j) if moving of course i to the slot of course j and taking out course j from its current slot creates either teachers conflict or room unavailability. Here we assume that all ot her courses will remain to be assigned to their existing slot. In another approach, we do not check the teachers conflict and room availability while creating an arc; rather, we check these constraints while enumerating cycles. This approach gives larger neighborhood th an compared to first variant, as even if an arc is infeasible, the resulting cycl e may be feasible, as the courses causing infeasibility might have been moved to so me other slot. The ne ighborhood size is also influenced by the maximum cycl e length and number of path s enumerated in a stage. Experimental results show that maximum cy cle length beyond 5 does not benefit much to the solution quality, whereas it increases running time consider ably. Similarly, a path of cost more than 100 or 10 % of current objec tive function value does not result into a negative cycle. Therefore, for our implemen tation we kept maximum cycle length limited PAGE 103 92 to 5 and enumerate paths with maximum cost min (100, 10 % of current objective function value). Table 4.1 shows the solution obtained and time taken for implementation of integer programming formulation and VLSN algorithm. We appl y 100 runs of VLSN local search on each test instance starti ng with initial solution obtained by greedy heuristic 2. We apply 5 runs of VLSN tabu search on each test instance starting with the best initial solution obtained in 20 runs of greedy heuristic 2. We keep stopping criteria of maximum (8 number of courses) iterati ons in each run and apply diversification strategy if there is no improvement in local be st in (3 number of courses) iterations. Table 4.1 Performance of VLSN and Tabu search Instance Statistics Optimal Solution VLSN Local (100 runs) VLSN Tabu (5 runs) Nperiods NSlots nCourses Obj. T(Sec) Best % Dev Avg % Dev T(Sec)Best % Dev Avg % Dev T(Sec) 20 20 10 299 0.188 0 2.34 0.375 0 0 0.094 20 30 20 447 0.406 11.63 45.41 0.625 0.45 10.45 0.453 20 30 15 298 0.109 1.34 1.34 0.438 1.34 1.34 0.281 20 20 30 2011 0.547 0 12.53 0.718 0 2.39 0.562 20 30 30 984 0.875 0 20.22 0.718 0 5.69 0.891 30 30 15 300 0.109 0 2.00 0.469 0 0.67 0.297 30 30 25 806 0.313 1.36 12.41 1.157 1.86 8.31 1.047 30 30 20 378 0.234 0 6.35 0.656 0.79 0.79 0.547 40 35 20 557 0.156 0 16.88 0.782 0 4.13 0.641 40 40 30 622 0.578 0 12.54 1.843 0 5.47 1.500 50 60 40 598 1.25 1.67 18.90 4.532 0.17 2.69 4.610 50 50 50 822 1.141 0 5.16 3.782 0.97 2.92 6.652 50 100 100 2036 2299 21.4541991 2047 43.828 50 100 150 4229 5105 30.6094164 4485 87.890 50 100 200 8482 10282 60.6108349 9105 157.73 In table 4.1, the columns nPer, nSlots, and nCour represent number of periods, number of slots and number of courses respec tively. We report the best objective function PAGE 104 93 value in column Obj. Fn. and average of final objective fu nction values in each run in column Avg. Obj. The total times taken by different approaches are shown in column T(Sec). From table 4.1 we can deduce following: Implementation of integer programming form ulation gives optimal solution in very short time for small size test instances. However, it fails to converge to optimal solution for larger test instances. VLSN local search has been able to achieve optimal solution in most of the cases; however, the average solution quality by this approach is not very robust. VLSN tabu search gives optimal solutio n in almost all the cases. The average solution quality is also much more superior to that in local search. However, this approach also takes more time as compared to local search. 4.10 Conclusion and Further Scope In this report we have outlined the general framework of course timetabling problem and has proposed two heuristics for ge tting initial solution and a heuristics for local search and tabu search using VLSN We observed that integer programming formulation can be used to obtain optimal solu tion for smaller test instances; however, as size of the instance increases one should use VLSN with tabu search, as it gives very robust solution in manageable time. Further wo rk is required in setting parameters for random instance generator to make it better representative of the real life problem instances. PAGE 105 94 CHAPTER 5 BLOCKTOTRAIN ASSIGNMENT PROBLEM 5.1 Introduction Transportation of goods by railroads is an integral part of the US economy. Railroads carry millions of shipments annually from their origins to their respective destinations. Planning and executing the trans portation of these shipments involves many complex decisions at different levels of opera tion. In the last few years, a growing body of advances concerning several aspects in rail freight, such as train scheduling, blocking plan, blocktotrain assignment, locomotiv e scheduling, etc., has appeared in the operations research literature (see, for exam ple, Cordeau et al. [1998], Newman et al. [2002], and Kraft [1998]). In this paper, we study a keyplanning problem, called the blocktotrain assignment (BTA) problem, which arises afte r shipments are classified into blocks and a train schedu le has been defined. We first gi ve the reader a brief overview of the railroad blocking probl em, the trainscheduling problem, and the blocktotrain assignment problem. The blocking problem: The first and foremost planning problem to be solved in railroad is the blocking probl em, which involves the grouping of shipments into blocks. In rail transportation, a shipment may pass through many classification ya rds on its route from its origin to its destination. To preven t shipments from being reclassified at every yard they pass through, se veral shipments may be gr ouped together to form a block. After a shipment is placed in a block, it is not recla ssified until it reaches the destination of that PAGE 106 95 block. The blocking problem is to identify this classification plan for all shipments in this network, called a blocking plan, so that the total shipment cost is minimal. The total shipment cost is the sum of two cost terms: the cost of classification of shipments and the car miles traveled by as they are carried fr om their origins to their destinations. Some recent references on blocking problem are: Ne wton et al. [1998], Bar nhart et al. [1998] and Liu [2003]. The train scheduling problem: Given a blocking plan, developing a train schedule (also called train timetable) is perhaps the ne xt most important opera tional planning task faced by a railroad. The trainscheduling pr oblem is to identify train routes, their frequencies (how often to run in a week), and their timetable s so as to minimize the cost of carrying cars from their orig ins to their destinations. Th is problem also includes the synchronization of realtime tr ain movement on the links of the physical railway network. The train schedule once developed is repeated every period (mostly weekly) of a given time horizon (may be a quarter year). Some re cent references on this problem are due to Farvolden and Powell [1994], Campbell [ 1996], Kraft [1998] a nd Barnnlund et al. [1998]. The blocktotrain assignment problem: Once the blocking plan and the train schedule have been developed, the next step is to determ ine which trains should carry which blocks. This problem is known as the blocktotrain assignment problem and this paper is concerned with the modeling and developing algor ithms to solve this problem. The goal in this problem is to have the block transported at minimum cost and with minimum delay from due dates while meeting the operational constraints. Based on the service level of a block, railroads may im pose penalty on the deviation from due date. PAGE 107 96 Ideally, in an optimal solution a block should take th e shortest train rout e from its origin to its destination. However, there are capacity and operational constraints, which must be satisfied by a solution. These re strictions make the blocktot rain problem difficult to be solved to optimally. Most railroads solve the BTA problem manually. We are aware of only few efforts in literature to solve the BTA problem exclusively. No zick et al. [1997] deals with the finite horizon, discrete time problem of minimizing the total vari able cost of moving cars given a fixed train sche dule and satisfying due dates. They have developed a procedure that involves iter atively solving a linear progr amming relaxation and rounding some of the resulting fractiona l values until a feasible inte gral solution is found. Another paper, Kwon et al. [1998] deals with the al gorithm to improve a gi ven blocking plan and blocktotrain assignment. They formulate the problem as a linear multicommodity flow problem and column generation technique is used as a solution approach. There are few papers, which consider the BTA problem indirectly in their problem formulation. Early work in this area goes to Thomat [ 1971]. It develops a cancellation procedure that gradually replaces direct sh ipments by a series of intermediate train connections to minimize operation and delay co sts. Crainic et al. [1984] has proposed a nonlinear, mixed integer, multicommodity flow model that deals with the interaction between blocking, makeup, and train and traffic routing deci sions. Haghani [1989] proposed a formulation and heuristic approa ch for a combinatorial train routing and makeup and empty car dist ribution problems. Keaton [ 1989] proposed a model and a heuristic method based on Lagrangian rela xation for the combined problem of car blocking and train routing and makeup. Th e subsequent paper by Keaton [1992a] also PAGE 108 97 considers the constraints for blocking a nd maximum transient time for each origindestination pair. An instance of the BTA problem is define d by blocks that must be delivered; the train schedule and the physical rail system will be used to deliver the shipments. The physical rail system consists of undirected links (one or more railro ad tracks) connecting stations (rail yards) and the trains carrying the shipments on this rail system. Based on the tractive power of the locomotives pulling a trai n and the links on whic h it runs, there are limitations on the carrying capacity the trains. Hence a solution of the BTA problem must honor these constraints. These constraints may be in the fo rm of maximum number of cars in the train, maximum length of the trai n, and maximum weight of the train and may differ on different section of the trains route. It can be noted that these constraints are interlinked, i.e., if number of cars in a trai n will be high then the length of the train will also be high and so the weight of the train. To make the BTA probl em somewhat simpler, railroads make some of these constraints as hard constraints while keeping other constraints as soft constraint s (that is, we allow some violation but there a penalty for each violation). In our discussions in this pape r, we have included one type of constraints in our formulation while assuming that others will be soft constraints and hence can be included in the objective functi on. Apart from these, ther e are many other operational constraints, which will be discussed in detail in the appropriate sections. Two of the inputs for the BTA problem are the train schedu le and the blocking plan. Generally a train schedule is for a weekly period and different tr ains have different running frequency in the schedule. For exampl e, one train runs ev eryday of the week while another runs only on Monday, Wednesday and Thursday. Similarly, not all the PAGE 109 98 blocks are made everyday and they may have different frequency. Therefore, ideally the blocktotrain assignment problem should be formulated for a peri od equal to maximum of the period of train sche dule and blocking plan Obviously, size of the BTA problem depends on this period and the longer the peri od the larger will be the problem. To keep the size of the problem manageable, we assu me that all the trains run every day and similarly all the blocks are made every day. This assumption allows us to define the BTA problem on the trains and blocks, which is repe ated every day. It can be noted that this assumption does not restrict the travel time of a train and the total time taken by a block to reach its destination to one day. To so lve the problem with longer periods, we can either generalize the approaches developed in this paper or can use heuristics to get the solution of the problem with longer period in the postprocessing phase. In the rest of the paper, we refer the daily BTA pr oblem simply as the BTA problem. In this paper, we give several inte gerprogramming formul ations of the BTA problem and develop several algorithms to so lve this problem. This paper makes the following contributions. We develop a spacetime network, which al lows us to formulate the BTA problem as a multicommodity flow problem with additional side constraints. We develop two integerprogramming formul ations of the BTA problem. The first formulation formulates the BTA probl em as a nodearc version of multicommodity flow problem in the spacetime network. The second formulation conceives the BTA problem as a pathflow version of the multicommodity flow problem in the spacetime network. The integer programming formulations can also be used to solve the BTA problem approximately. Using these formulations we suggest a Lagrangian relaxation algorithm that can solve the BTA problem to nearoptimality within a few seconds. We also propose a greedy construction al gorithm that proceeds by assigning blocks to trains one by one until all blocks are assigning. By choosing different rules for selecting blocks and assigning them to trains, we can obtain different implementations of this generic approach. PAGE 110 99 We perform extensive computational inve stigations of these algorithms. Our integer programming formulations can solve the BTA problem to optimality within a few minutes of computer time using CP LEX 8.1. The Lagrangian relaxation algorithm can solve the BTA problem to n earoptimality within a few seconds. The computational performances of greedy c onstruction algorithms are also very attractive. To summarize, this paper reports the first se rious effort to utilize the stateoftheart mathematical programming techniques to solve the BTA problem. We demonstrate that this problem can be solved to optimality or nearoptimality within reasonable computational times. 5.2 Mathematical Programming Formulations In this section, we describe the spaceti me network so that the BTA problem can be formulated as a flow problem in this network. This spacetime network is similar to the spacetime network developed by Ahuja et al. [2002b] to model the flow of locomotives in a train network while modeling the locomotive scheduling problem. 5.2.1 The SpaceTime Network We define the BTA problem on a spacetime network, in which every node is distinctly identified by place, time and tr ain. We denote the spacetime network as G = (N, A), where N denotes the node set and A denotes the arc set. We construct the spacetime network as follows. A train runs from its origin to its destination at prespecified time schedule and stops at many stations on the route to perform different activities, like, maintenance checkup, refueling, dropping or adding shipment s, etc. However, a block can be attached to or detached from a train at a station, onl y if the station is origin of the train, or destination of the train, or a stop where ship ments can be added or dropped. We call all PAGE 111 100 such stations as validstops for the train and we call the itinerary of a train between two consecutive validstops as a trip. Hence, the route of a train consists of several trips. In practice, a train route may contain as many as 20 trips. The trip l of a train is represented by a trainarc (l l ) in spacetime network. The tail node l of the arc denotes the event for the departure of train at the or igin of the trip and is called the traindeparture node. The head node l denotes the arrival event of the train at the destination of the trip and is called the trainarrival node. To connect two consecutive trip s of a train, we create trainconnecting arcs from the trainarrival node of a train at a station to the traindeparture node of the same train at the same station. The time difference between these two nodes represents the stay time of the train at the station. To allow a block to be added or dropped from a train on its route, we introduce ground nodes and connecting arcs. For each trainarrival node, we create a corresponding groundarrival node with same station, time and train attributes as that of trainarrival node. Similarly, for each traindeparture node, we create a grounddeparture node with the same station, time and train attr ibutes as that of tr aindeparture node. We connect each trainarrival node to the asso ciated groundarrival node by a directed arc called arrivalconnection arc. We connect each grounddeparture node to the associated traindeparture node by a directed arc called departureconnection arc. To facilitate the switching of a block from one train to another train at a station, we sort all the ground nodes at each station in the chronological orde r of their time attri bute, and connect each ground node to the next ground node in th is order by direct ed arcs called ground arcs (we assume here without any loss of generality that ground nodes at each station have distinct PAGE 112 101 time attributes). Further, to allow a block to take a train in next period, we connect the latest ground node at a station with the earliest ground node at the same station. To summarize, the spacetime network G = (N, A) has three types of nodes, namely trainarrival nodes, traindeparture nodes and ground nodes (comprising of two subtypes: groundarrival nodes a nd grounddeparture nodes). The ar c set A consists of three types of arcs, namely train arcs, connection arcs and ground arcs. The set of connection arcs further consists of three types of arcs : arrivalconnection arcs departureconnection arcs and trainconnection arcs. We also asso ciate a cost with each arc, which depends on the block flowing on the arc and is discussed in detail in Se ction 5.2.2. In the spacetime network, there is one incoming arc and at mo st two outgoing arcs at each trainarrival node, at most two incoming arcs and one outgo ing arc at each traindeparture node, and at most two incoming arcs and two outgoing ar cs at ground nodes. It implies that number of arcs in the spacetime ne twork cannot be more than do uble of the number of nodes, which makes the network G highly sparse. A part of the spacetime network for BTA problem is shown in Figure 5.1. With respect to the spacetime network G, we define the BTA problem as follows. We say a grounddeparture node as the originnode of a block, if it is the earliest grounddeparture node after the release time of the block at the origin station of the block. At the destination station of a block, we sa y the set of groundarrival nodes as the destinationnode set of the block. With these notations, we define the BTA problem as to find a path for each block from its originnode to any of the node in destinationnode set in spacetime network G, such that the overall cost is minimum and all the operational requirements are satisfied. We claim that a ny solution to the above problem will also be PAGE 113 102 the solution of the BTA pr oblem. To prove this, let p is such a path in the spacetime network, which is assigned to a block in th e solution. As only tr ainarcs are having different station attributes at its tail node and head node, the path p must consists of trainarcs along with other arcs. Fu rther, two trainarcs in p must be connected by one or more connection arcs and/or ground arcs, which imp lies that a block will take a train at the origin station, will go to another station and fr om there it will either take the same train or will be switched to some other train, and so on, until it reaches the destination station. Therefore, path p in spacetime network represents a feasible train rout e for the block, which proves our claim. Figure 5.1 A part of the spacetime network Train arcs: Connection arcs: Ground arcs: Train Arcs Train 1 Train 2 TrainArrival Nodes Train 2 Train 3 Train 4 TrainDeparture Nodes Train Arcs Ground Nodes Time PAGE 114 103 We have developed two integerprogram ming (IP) formulations of the BTA problem in the spacetime network G. In one formulation, decision variables are whether a block will flow on an arc or not; whereas in the other formulation the decision variables are whether a block will take a path in spacetime network or not. In Sec tions 5.2.2 and 5.2.3, we discuss these two formulations. 5.2.2 ArcBased IP Formulation In Section 5.2.1, we have defined the BTA problem as to find a path in the spacetime network, from originnode of each bloc k to one of the nodes in its destinationnode set. This can be formulated as an IP in wh ich a decision has to be made whether an arc will be in the path for a block. We use the following notation in our arcbased IP formulation. Parameters: G: spacetime network N: set of nodes in G and we denote a node by i A: set of arcs in G and we denote an arc by a B: set of blocks and we denote a block by b vb: number of cars in block b B Ob: originnode of block b in G Ai +: set of outgoing arcs at node i N Ai : set of incoming arcs at node i N Db: destinationnode set of block b in G Ca: capacity of arc a A a b: incidence taking value 1 if block b B can flow on arc a A PAGE 115 104 ca b: unit cost of flowing block b B on arc a A Decision variables: 1 if block bB flows on arc an aA 0, otherwiseb ax Formulation (5.1) min B bA a b a b a bx c v 5.1a subject to 1 b OA a b a b ax B b 5.1b 0 i iA a b a b a A a b a b ax x b bD i O i N i B b , 5.1c 1 b iD i A a b a b ax B b 5.1d a b a B b bC x v A a 5.1e 1 0 b ax B b A a 5.1f In formulation 5.1, the constraints 5.1b, 5.1c and 5.1d ensure that there exists a path for each block from its originnode to a node in its de stinationnode se t. Constraints 5.1e ensure that the total number of cars flow ing on an arc is no more than its capacity. If an arc is a trainarc, we set its capacity equa l to the capacity of associated train on the segment represented by station at tribute of the tail node and th e head node of the arc. For all other arcs, we set its capacity equal to infi nity. It can be noted th at in the absence of constraints 5.1e, the BTA problem reduces to finding a shortest path for each block. PAGE 116 105 Now, we will describe the objective functi on 5.1a in formulation 5.1. We have associated a cost with each arc and block comb ination. This cost represents the cost of flow of a car in a block on an arc. In the case of trainarcs this cost depends on the carmiles, transittime and service level of th e block. In the case of ground arcs, the cost ca b depends on the waiting time (duration of arc a) of block b at the associated station and the service level of block b. We do not associate any cost with any of the arrivalconnection arcs and trainconnection arcs. Howe ver, at all the stations other than the origin station of a block, if the deci sion variable corresponding to the departureconnection arc takes value 1, it signifies that the block is switching th e train and a cost is incurred in detaching the block from a train an d in attaching it to other train. Therefore, we associate a cost equal to the train switc hing cost with all such departureconnection arcs. In formulation 5.1, we have also associated an incidence scalarb a with each arcblock combination. This is re quired to prevent a block bein g assigned to an arc, which violates any of the operational requireme nts. Generally, to meet the government regulations or safety requirements (like, wi dth, height or weight limitations, hazardous material limitation, etc.) a bloc k may not be associated to a tr ain or to any train traveling on a particular link. We set the incidence b a = 0, for all such arcs and blocks combination and set b a = 1 to all other combinations. Railroads may also set the incidence b a = 0, if they do not want to assign a bl ock to a trip for hi storical reasons or obligation to customers. We will now briefly discuss the size of th e problem 5.1 for a typical US railroad application. The number of st ations in the physical networ k of the railroad is around PAGE 117 106 6,000 and there are around 350 trains running da ily. Altogether, these trains are having around 2,000 valid stops, which results in to a spacetime network with around 7,000 nodes and 10,000 arcs. The railroad makes around 1,200 blocks every day. It implies that there will be around 10,000 1,200 = 12 million binary decision variables and around 1,200 7,000 + 10,000 8.5 million constraints. It is almost impossible to solve a problem of this magnitude with commer cially available optimization software. In addition to the difficulty of this pr oblem due to its magnitude, there are many operational requirements, which are difficult to be included in formulation 5.1. For example, a block is not allowed to travel on mo re than a certain number of trains from its origin to its destination. It is difficult to include this constraint in the formulation. Similarly, instead of preventing a block from being associated with a trip, railroad may want that a block should not take certain routes, or they may want that a block should take a route one amongst the limited number of routes. It is difficult to include this constraint in the formulation (1). We next describe another IP formulation, which overcomes many of these difficulties. 5.2.3 Pathbased IP Formulation This formulation requires that for each bl ock we have already enumerated a set of feasible train paths and decision to be made is which path to use for each block. We call a path a validpath for a block if it is from the origin node of the block to a node in the destinationnode set of the block in the sp acetime network and meets all the operational requirements, except the capacity constraint on the train. In other words, a block cannot be assigned to a validp ath only if it violates the capacity constraint on tr ainarcs in the spacetime network. In the rest of the paper, we will interchangeably use validpaths and PAGE 118 107 paths. In addition to the notati ons used before, this formulation uses some more notations. We describe these new no tation next followed by the pathbased formulation. Parameters: Pb: set of feasible paths for block b B in G P: set of all the paths, i.e., P = b bBP P a: {path p: p P and arc a A is on p} cp: cost of assigning path p P to the corresponding block Decision Variables: 1 if path isassignedtothecurrentblock 0, otherwiseppP x Formulation (5.2) min B bP p p pbx c 5.2a subject to, 1 bP p px B b 5.2b a P p b P p p bC x vb a : A a 5.2c 1 0 px B b P pb 5.2d In the above formulation, constraints 5.2b ensure that exactly one path is assigned to each block and constraints 5.2c ensure that the train capacity is not violated. In this formulation, the number of decision variab les depends on the possible number of paths for each block. In Section 5.3, we discuss the path enumeration algorithms and for the problem under consideration the number of decision variables can go up to 500,000. As we have constraints for each block and each arc in the spacetime network, the number of PAGE 119 108 constraints can be at most 12,000. It can be obs erved that the number of constraints in the above formulation is much less th an that in formulation 5.1. An important advantage of the formulation 5.2 over formulation 5. 1 is in defining the objective func tion. One of the important goals fo r the railroads is to deliver the shipments before due date to the customer. In the pathbased formulation, we know that at what time a block is reaching its destinati on and hence we can define the path penalty proportional to the deviation from the due date. Hence, the cost associated with each path (cp) takes into account the carmiles, deviation from due date, service level of block and switching cost of blocks. This cannot be done in the arcbased fo rmulation of the BTA problem. Another important advantage of this formulation is that we can eliminate those paths, which do not meet operational constraint s by simply not listing them in the set of valid paths. We cannot eliminate such possibilities in formulation 5.1. However, in order for this formulation to be of reasonable size, we should not allow the set of valid paths to be too large. The total number of train paths for a block in the spacetime network can be exponentially larg e, but we need to identify a small set of valid paths that the block is likely to use. Using some preprocessing, we can easily identify the set of valid paths for each block. This is consistent with the current practice at railroads. Railroads do not allow blocks to be sent along arbitr ary paths in the train network, but there are many restrictions th at these paths must follow which keeps the number of possible paths manageable. In the rest of this paper, we will use the pathbased formulation of the blocktot rain assignment problem. PAGE 120 109 5.3 Path Enumeration Algorithms Our pathbased formulation for the BTA problem requires that we be able to identify the set of valid paths for each block. In this section, we describe an algorithm to determine these sets of valid paths. In our algorithms, we first identify all the distinct originnodes (of blocks) and then identify all th e feasible paths from each of the distinct originnodes to all other nodes in network. Therefore, instea d of identifying valid paths for each block separately, we identify valid paths for all the blocks originating at a specific originnode in one go. To restrict th e number of valid paths for a block, we use the following criteria used by the railroads: (i ) limit the number of arcs on the valid path to a number less than a specified parameter (say, example, 3); and (ii) limit the length of the valid path to a specific percentage over th e shortest path in the train network. These restrictions are motivated by the curre nt practices followed by US railroads. Our algorithms for path enumeration are the variation of the algorithms to find shortest paths from a source node in a gra ph (see Ahuja et al. [1993] for detail). The algorithms to find shortest paths from a node, we calculate the cost in reaching other nodes from the source node through all possible paths and pick up the path for which this cost is minimum. In these algorithms, we k eep a label with each node, which denotes the cost of some path from the source node to a node and algorithm updates it with lower value whenever a better path is found. In our algorithms, we keep more than one label at a node and every time we find a new path to a node from the source node, we add a new label. A label is uniquely identified by the node with which it is associated and its predecessor label. The attributes of each label and their sign ificance are discussed later. PAGE 121 110 5.3.1 Path Enumeration Restricted by Number of Trains Generally, railroads put an upper bound on the number of trains (say t) a block can take from its origin to its destination. In our algorithm, we enumerate all the paths and ignore all those paths, which take more than t number of trains. Our algorithm works as follows. We start from a di stinct originnode (say, the source node) and visit the head nodes of all the arcs emanating from this nod e, then visit the head nodes of all arcs emanating from the nodes just visited, and so on, until the number of trains visited is more than t. In the above procedure, every time we visit a node, we create a new label at the node; a label signifie s that there exists a unique path from the source node to that node. With each label, we keep a list of attributes, such as, the node associated with this label, the distance from the source node to the node, transit time, number of trains visited, predecessor label (the label using which the cu rrent label is created), etc. Since a node can be visited multiple times, a node can have multiple labels. To initialize the algorithm, we create a label at the sour ce node and set its predecessor to null signifying that no path exists from the source node to itself. Figure 5.2 describes this path enumeration algorithm. Figure 5.3 gives an example of the path enumeration algorithm restricted by the number of trains. Suppose that node 1 is the or iginnode and we need to enumerate all the paths from node 1. We create a label l11 with null predecessor at node 1 and with respect to this label, we examine arcs (1, 2) and (1, 3) and create labels l21 and l31 respectively. Next, we take label l21 and with respect to this label create labels l31 and l41 at nodes 3 and 4, respectively. Finally, labels l42 and l43 are created with predecessor labels at node 3. All the labels and their predecessor paths are also shown in Figure 5.3. It can be seen that by PAGE 122 111 tracing the predecessor path of each label, we get the unique pa th from node 1 to the node associated with label and the algorithm has been able to enumerate all the paths. algorithm path enumerationI; begin for each originnode (say s) create a label with null predecessor at s; make a set S := {s}; while S do remove a label l from S; for each of the outgoing arcs (say a) from the node associated with l if number of trains taken to reach head node of arc a by the path represented by l and a is not more than t then create a new label at the head node of a and add the new label to S; end. Figure 5.2 Path enumeration algorithm restricted by the number of trains. Figure 5.3 Example of path enumeration algorithm. Theorem 1: Algorithm in Figure 5.2 enumerates all the paths from the originnode of a block to the nodes in the destinati onnode set, on which the number of trains taken is not more than t. 1 2 3 4 l11 l21 l31 l32 l42 l41 l43 null PAGE 123 112 Proof: We prove the above theorem by contra diction. Suppose that there exists a path p (p1p2p3 pk) from origin node (say p1) to a node (say pk) in the destinationnode set, for which the number of trains is not mo re than t but our algorithm fails to find this path. As there exists an arc (p1, p2) and we created a label at p1, we must also have created a label at p2. Similarly, as we must have created a label at p2, the algorithm examines the arc (p2, p3) and must have created a label at p3, and so on, unless the number of trains taken in reaching a node (say pi in path p) is more than t. Now, if the number of trains taken to reach pi is more than t then we cannot reach pk using path p and by taking trains not more than t, which contradicts our claim of existence of the path p. Our experiments on the reallife data pr ovided by a major US railroad indicate that the above algorithm is able to enumerate all the feasible paths for all the blocks in few seconds of computational time on Pentium IV computer. However, we found that for some blocks the numbers of valid paths cr eated by the algorithm described earlier were too large. For such blocks, we needed anot her criterion to reduce the number of valid paths further. We found that though severa l valid paths rode on no more than the specified number of trains, their paths are quite circuitous in terms of the distance of the path and/or the time taken by those paths. Railroads typically do not allow blocks to follow paths, which take more than a specific percentage over the shortest path (either distance or time). We next de scribe in Section 5.3.2 a me thod that keeps track of the distance of each valid path (or the time taken by it) and eliminates the paths, which do not meet the specified criteria. PAGE 124 113 5.3.2 Path Enumeration Algorithm Restricted by Number of Trains and Number of Labels In the path enumeration algorithm discu ssed in Section 5.3.1, we make a simple modification and will keep only prefixed number of labels (say k) at a node. Hence, in this algorithm, whenever we create a new label at a node already having k labels, we replace one of the existing labels with new label, if it is of better quality. Now if the label to be replaced is already the predecessor of ot her labels, then the re placement of this label with the new label will make the entire succe ssors of old label obs olete and hence should be removed. However, computationally it is very difficult to trace and remove all the successors of a label. Therefore, we need a mechanism in which a label is not required to be replaced, if it is the predecessor of other labels. In other words, before making a label predecessor of other labels, we ensure that the label is perm anent and will not be replaced by any other label at any point of time. In the BTA problem all the cost figures are nonnegative, which motivated us to develop an algorithm analogous to Dijkstras algorithm. In D ijkstras algorithm, we keep a distance label with each node, which denotes the cost of some path from source node to the node, and we iteratively make the distance label of a node perman ent if this is the smallest distance label among a ll the distance labels, which are not permanent. In our algorithm, we also make a label permanent and enumerate paths from this label, if this is the best label among all the labels from which paths are requir ed to be further enumerated. This path enumeration algorithm is given in Figure 5.4. PAGE 125 114 algorithm path enumerationII; begin for each originnode (say s) create a label with null predecessor at s; make a set S := {s}; while S do take out the best label from S; for each of the outgoing arcs (say a) from the node associated with l if number of trains taken to reach head node of arc a by the path represented by l and a is not more than t then create a new label at the head node of a and add the new label to S; end. Figure 5.4 Path enumeration algorithm restri cted by number of tr ains and number of labels. Theorem 2: In the above algorithm, a perman ent label will never be replaced by any new label. Proof: We will prove this theorem by contradiction. Let us assume that we can keep at most k labels and there are k labels at a node i. Let us claim that there is k+1th label (say l1) at node i, which is of better quality than a permanent label (say l2) at node i. Now there exist two scenarios. In scenario 1, both the labels l1 and l2 have already been there at the time l2 was made permanent and as we decided to make l2 permanent, l1 can not be better than l2, which contradicts our claim. In scenario 2, label l1 was not there when l2 was made permanent. It implies that predecessor of l1 must be made permanent after l2 was made permanent and hence quality of predecessor of l1 cannot be better than that of l2. Also as all cost figures are nonnegative, quality of l2 cannot be better than that of any of its predecessor. It implies that l1 cannot be of better quality than l2, which again contradicts our claim. As our spacetime network is very spar se, the size of the labels to be made permanent remains manageable if the maximum nu mber of labels to be kept at a node is PAGE 126 115 not too high. However, this algorithm may e xhibit inferior performance than that in Section 5.3.1, if the upper bound on the number of labels at a node is too high or is unrestricted. This is because of the mainte nance of expensive data structure in the implementation of this algorithm. In the algorithms discussed in Sections 5. 3.1 and 5.3.2, we have been able to put some of the operational requirements while en umerating the paths. However, there are many other constraints (like hei ght, width and weight restriction, etc.), which must be honored by a valid path. We can easily meet th ese requirements by simply ignoring all those paths, which violate any of these re quirements. Therefore, in the pathbased formulation there is no need of incidence v ector and we can eliminate a large number of decision variables in the preprocessing stage. In following section, we describe few approaches to solve the pathbased formulation of the BTA problem. 5.4 Solution Approaches We have seen that the pathbased form ulation is more manageable than the arcbased IP formulation. Therefore, all of our approaches to solve the BTA problem are based on solving the pathbased IP formulati on. We have developed three approaches to solve the BTA problem: (i) solving to optima lity using CPLEX, (ii) solving heuristically using Lagrangian relaxation, and (ii) so lving heuristically using greedy algorithm. To solve the BTA problem optimally, we used the commercially available CPLEX optimizer. It can be recal led that there are two types of constraints in pathbased formulation: (i) each block mu st be assigned to a path, and (ii) a trai n cannot carry more cars than its capacity (called capacity constraint). Although, capacity constraints make a formulation hard, CPLEX has been able to so lve a reallife BTA problem to optimality in PAGE 127 116 a very reasonable amount of time. This ma y be because of abunda nt train capacity on a large number of paths. Therefore, we cannot ge neralize this behavior of optimizer to all instances of the problem. Also, the bran ch and bound algorithm of CPLEX may not be scaleable to large problem size and may not be able to solve the problem of higher magnitude to optimality. In Section 5.4.1, we describe a heuristic to solve the BTA problem based on Lagrangian relaxation an d in Section 5.4.2; we discuss the neighborhood search heuristic. 5.4.1 Lagrangian Relaxation Based Heuristic It should be noted that in the absen ce of capacity constraints, the BTA problem reduces to assigning the least cost path to each block, which is trivially solvable. In other words, capacity constraints make the problem hard to solve. One common way of dealing with such constraints of an IP problem is Lagrangian relaxation (F isher [1981]). In this approach, capacity constraints are relaxed a nd a penalty is assigned for each unit capacity violation on each train arc. The relaxed problem is solved iteratively with different value of penalty. If in iteration, no arc capacity is violated, the solution is optimal. On the other hand, if capacity constraints are violated on some arcs, we try to get feasible solution using heuristic approach and repeat the proce ss with new penalty. Our heuristic to obtain feasible solution is similar to the greedy al gorithm described in Section 5.4.2. In our heuristic, we first identify all the trains ar cs on which capacity is violated and do not assign any path to the blocks flowing on these arcs. This step frees the capacity on few arcs and reduces the BTA problem to assigni ng path to few bloc ks. We apply greedy algorithm as in Section 5.4.2. to assign path s to these blocks. In the remainder of the section, we only describe Lagrangian relaxa tion approach as greedy algorithm is already discussed in Section 5.4.2. PAGE 128 117 In the Lagrangian relaxation approach, th e capacity constraints are relaxed and we assign penalty ua 0 for each unit capacity violation on every train arc a and add the constraint to the objective func tion of pathbased IP formulation. This may be interpreted as putting penalty ua for each unit of violation of trai n capacity in a trip. The relaxed problem is given in formulation 5.3. Formulation (5.3): min B bA a a P p b P p p b a P p p pC x v u x cb a b) (: 5.3a subject to 1 bP p px B b 5.3b 1 0 px B b P pb 5.3c Let us call the above problem as Pu. It can be observed that for a fixed penalty vector u, problem Pu separates into independent subproblems for each block b and is given in formulation 5.4. Formulation (5.4): min p a p b a P p p px v u x cb) ( 5.4a subject to, 1 bP p px 5.4b 1 0 px ,bP p 5.4c For the above problem, we only need to assign the best path with respect to objective function (4a) to a block, which is very easy and can be solved trivially. Let (u) is the optimal solution to Pu for u 0, then it is standard that (u) provides a lower bound PAGE 129 118 to the BTA problem. Now the problem of fi nding the best possible lower bound is the so called Lagrangian dual problem and is given in formulation 5.5. Formulation (5.5) max (u) 5.5a subject to u 0. 5.5b The approach to solve the problem 5.5 works as follows. We perform multiple iterations with different value of u. At iteration k, using the multipliers, uk, we first solve the relaxed problem (Pu k). Thereafter, we update the multi pliers to get closer to the solution of 5.5 and repeat the procedure. In an iteration, we update the multipliers using the subgradient optimization technique. The basic idea of updating multiplier is simple. If in an iteration, capacity is violated on some arc (say a), then we set ua k+1 > ua k making the flow on this arc una ttractive and else if ua k > 0 and there is no capacity violation then we set ua k+1 < ua k making the arc more attractive. We use the following formula to update multiplier ua for an arc a in iteration k. ua k+1 = ua k + k ) (:a P p b P p p bC x vb a where k is the step length in iteration k and is calculated using following formula: k = 2 :  ) (a P p b P p p b k kC x v UBb a where k is a scalar UB is an upper bound on the optimal solution of original problem. y = 2 1 2) (iy. PAGE 130 119 Above algorithm to obtain lower bound of the BTA problem is summarized in Figure 5.5. Our computational experiment on the reallife data indicate that very good lower bound can be obtained for the BTA problem in few seconds of computational time. algorithm lower bounding scheme; begin set k := 1 and initialize u1 by setting ua 1 := 0 for all a A; set stopping criteria; until stopping criteria is met do assign the best path to a block using objective function (4a) find the new multiplier uk+1; set k := k+1; end. Figure 5.5 Lagrangian relaxation based lower bounding scheme. 5.4.2 The Greedy Heuristic Algorithms As noted earlier, the optimal soluti on obtained by the in teger programming software CPLEX may not be scaleable to th e size of the problem or the computational time taken in solving the problem may not be acceptable to the railroads. Further, railroads are more interested in using an algorithm which is time efficient, is flexible in adjusting parameters and requirements and produces good quality solution (not necessarily optimal solution). Lagrangian rela xation based heuristic is one approach to get such solution. In this section, we ha ve developed another approach called greedy algorithm, which is able to produce nearoptimal solution very efficiently. Greedy heuristics are the si mplest and the most popular algorithmic approach to solve combinatorial optimization problems. In this approach, values are assigned to decision variables based on certa in greedy criteria. There ar e many variations of this algorithm in use. For BTA problem, we have developed an approach, which is able to produce reasonably good quality solution in frac tion of second. In addition, our algorithm PAGE 131 120 provides flexibility to users to specify the greedy criteria at run time, a feature highly desired by the railroads. In our greedy algorithm, we first order the blocks, based on certain priority and then assign the best feasible path to blocks in that order. At each stage of assignment of path, we ensure that train capacity is not violated on any trainarc in the spacetime network. Figure 5.6 describes our greedy algorithm. algorithm greedy heuristic; begin order the blocks in descending order of some priority; for each block b in above order order the paths in ascending order of path cost; for each path p in above order if none of the train capacity is violated on path p by assigning block b then assign path p to block b; update the available train capacity of train arc on path p; exit for; end. Figure 5.6 Greedy Heuristic for the BTA problem. In the above algorithm, the user has the fl exibility to specify the criteria on which blocks can be prioritized. Generally, there are many factors, which govern the criteria selection. Some of these criteria are listed below. Give high priority to the blocks having less number of possible paths, so that, for each block there will be sufficien t options while selecting a path. Give high priority to the bloc ks having more number of cars, so that, best path can be assigned to them and cost can be minimized. Give high priority to the blocks cont aining shipments with high service level (service level signifies the importance of the shipment), so that, paths with least transit time can be assigned to high priority blocks. Give high priority to the blocks for whic h shortest distance from their origins to their respective destinations is maximum. This will give the possibility of assigning best paths to the blocks traveling long distance. As discussed earlier, there are other capacity constraints on trains, which are soft constraints (such as, the maximum length of a tr ain, the maximum weight of a train, etc.). PAGE 132 121 However, railroads may want to make these constraints hard on all the trainarcs or on some of them. The greedy algorithm gives the flexibility to users to impose these restrictions strictly. This can be simply incorporated by adding extra checks at the time we check the feasibility of a path in the gr eedy algorithm. In Section 5.5, we report the computational result for the case when we pr ioritize the blocks by the number of possible paths in the spacetime network. 5.5 Computational Results In this section, we presen t our computational results. We were supplied data for reallife problem by a major US railroad. All of our algorithms were implemented in C++ and we used ILOG CPLEX 8.1 to solve the BTA problem to optimality. The algorithms were run and tested on a Xeon PC with a speed of 2.4 GHz. Railroad network of the company, which provided us the data, consists of over 6,000 stations and over 6,000 li nks connecting the stations. To transport the shipments, they run around 350 trains daily. The total number of stops in the network, where a block can be attached or detached from a train is around 2,000. On the average, the railroad makes over 1,200 blocks every day. It can be recalled from Section 5.4 th at to solve the BTA problem, we have developed three approaches: (i) solving to optimality, (ii) solving heuristically using Lagrangian relaxation, and (i ii) solving heuristically usi ng greedy method. In our first investigation, we measured the quality of solution obtained by implementing these approaches. For each of them, Table 1 gives the percentage deviation from the optimal solution and the computational time taken. Th is analysis corresponds to the case when we put a restriction of ma ximum 50 labels at each node in path enumeration algorithm. It can PAGE 133 122 be seen that the Lagrangian re laxation approach seems to be the best approach, as it gives very good quality solution in few seconds of computational time. The time efficiency of Lagrangian relaxation approa ch can be attributed to th e fast convergence of the algorithm. For the instance under considerati on, it converges to be st lower bound in 143 iterations. It is worthy to note here that the lower bound of the optimal solution is only 0.03 % away from the optimal solution. Table 5.1 Performance analysis of different solution approaches. Approach % deviation from optimal sol. Computational Time Optimal Approach 0% 97 sec Greedy Approach 1.43 % 1 sec Lagrangian Relaxation Approach 0.78 % 6 sec In Section 5.3, we have discussed two algorithms for path enumeration in spacetime network. We also indicated that both al gorithms have their own merit, depending on the restriction imposed on the number of path s to be enumerated. We analyze here the performance of these two algorithms. It can be recalled that to restrict the number of paths enumerated, we put the restriction on the maximum numbe r of labels at a node in spacetime network. Figure 5.7(a) and Figure 5.7(b) give the total number of paths enumerated and time taken in path enumerati on respectively for differe nt level of number of labels at a node. It can be seen that when we put restric tion on the number of labels, the time taken in path enumeration proportiona tely increases with the number of paths enumerated. However, as noted earlier, if the restriction on maximum number of labels at a node is too high, it is a bett er option to remove this rest riction and enumerate all the paths. For example, the time taken in path enumeration is around 30 seconds when we PAGE 134 123 keep at most 70 labels at each node, whereas this time is only 22 seconds in unrestricted case. We have also analyzed the quality of op timal solution where we put the restriction on maximum number of labels at a node. Figure 5.7(c) gives the percentage deviation of the optimal solution in each case from the op timal solution of the problem, when there is no restriction on number of paths. In Figur e 5.7(d), we also report the time taken by CPLEX in each case. It can be seen that alt hough the quality of the solution is the best when we do not put any restriction on number of labels at a node; how ever, it also takes more computational time. In a reallife implementation, us ers can tradeoff between the optimal objective function value and computationa l time. Further, even if the quality of optimal solution is slightly better in the unrestricted case, railroads may put restriction on the number of paths as they do not want to assi gn an inferior path in terms of transit time to any of the block. This may be required for the customer satisf action or goodwill of the organization. # of Paths Vs # of Labels0 100000 200000 300000 400000 5000001 0 2 0 30 40 5 0 6 0 70 u nrestricted# LabelsTotal Number of Paths Time in path enumeration Vs # of Labels0 5 10 15 20 25 301 0 20 30 4 0 5 0 60 70 unres t rict e d# LabelsTime in Seconds (a) Total number of paths enumerated (b) Time taken in path enumeration Figure 5.7 Analysis of performance w ith respect to maximum number of labels at a node. PAGE 135 124 % deviation from best sol Vs # of Labels0.00 0.40 0.80 1.20 1.60 2.001 0 2 0 30 40 50 60 70 unrestric t ed# LabelsPercentage deviation Optimization time Vs # of Labels0 50 100 150 200 250 30010 20 30 40 50 60 70 unrestr i cted# LabelsTime in Seconds (c) Percentage deviation from best solution (d) time taken in solving to optimality Figure 5.7Continued Analysis of pe rformance with respect to maximum number of labels at a node. 5.6 Summary and Conclusions In this chapter, we consider the bloc ktotrain assignment problem in rail transportation. We defined the problem on a spacetime network and gave two integer programming formulations. Of these two formulations, the pathbased integer programming formulation is more flexible and contains much less variables and constraints. One basic assumpti on in the pathbased formula tion is that a block can be assigned only to a predetermined set of paths; however, the user can specify the number of such paths. This formulation is manag eable enough to be solv ed to optimality in reasonable time using a commercial IP solv er. We also developed several heuristic algorithms to solve the BTA problem, including a Lagrangian relaxation based heuristic and several greedy construction heuristics. A ll of these heuristics proceed by enumerating feasible shortest paths and augmenting fl ows along those paths. We performed an extensive computational testi ng of this algorith m on the data provided by a major US railroad and report the result s of this testing. We found that the intege r programming formulation of the BTA problem can be solv ed to optimality within two minutes using CPLEX 8.1, the Lagrangian heuristic algor ithm can obtain a solution within 1% of the PAGE 136 125 optimal solution in about six seconds, and greedy construction heur istics can obtain solutions within 2% in just one second. Hence, we find that the daily version of the BTA problem, we considered in this paper, is quite efficiently solvable. We plan to pursue research to solve the weekly BTA problem, which is substantially larger, the daily BTA problem. PAGE 137 126 CHAPTER 6 TRAIN SCHEDULE DESIGN PROBLEM 6.1 Introduction As described in Chapter 5, railroads are full of large size and complex optimization problems. The arena of optimization problems in railroad application consists of the blocking problem, the train scheduling probl em, the locomotive assignment problem, the crew scheduling problem, the blocktotrain assi gnment problem, etc. In this chapter, we study one of the most important railroad op timization problems, called, train scheduling problem. This problem can be defined in th ree contexts: (i) zero based train scheduling problem in which decisions are made at the planning level about the trains to be made, their frequency, timings, routes, etc., (ii) incremental train scheduling problem, in which decisions are made about the change in existin g train schedule to im prove efficiency in transportation for given demand of shipments on planning level, and (iii) real time train scheduling problem, in which decisions are ta ken based on the positions of trains at a given time. In this chapter, we study zero based train schedul ing problem and we call this problem as train schedule design problem. In the train schedule design problem, the inputs are the blocking plan (refer chapter 5 for details) and the physical rail network. The decisions to be made are with respect to formation of train and deciding their schedule. The constraints are on the capacities of a train and number of trains that can stop at a station. Apart from these, there are many practical constraints which mu st be honored by a solution of the problem. Considering a PAGE 138 127 real life instance as in chapter 5, there will be millions of decision variables and constraints. Therefore, almost all the efforts in the past ar e made in developing heuristics to solve the train schedule design problem to nearoptimality. We next describe the literatures available to solve th e train schedule design problem. The train schedule design problem consists of two parts (a) developing the train schedule, (b) routing the shipments. The problem of routing the shipments is a multicommodity flow problem with side constraint s. Generally a blocking plan is developed first and then the routings of the blocks ar e developed. Chapter 5 of this dissertation addresses the problem of block to train assi gnment and we have given the related past research to solve this problem in same chapte r. In this chapter, we briefly present the literatures discussing the deve lopment of train schedule. Thomet [1971] considers the tradeoffs of replacing direct trains with a series of connections ignoring classifi cation strategies. Morlok a nd Peterson [1970] develop timetables for a railroad without defining trai n connections and frequencies. Assad [1982] is concerned with optimal stop schedules of trains, but gives no consideration to routing of shipments. Van Dyke [1986] and Bodi n et al. [1980] are concerned with the classification of shipments given a set train sc hedule, but do not consider the effect of classification strategy on the schedule. Assad [1980] was the first to consider the train itinerary and makeup problems in the same study. He suggests Benders decomposition as a possible approach to solving his model; however, he does not solve it. More recent studies have attempted to model and solve the joint train sc heduling, shipment routi ng problem: Keaton [1989, 1991, 1992a, 1992b], Crainic et al. [1984], Crainic a nd Rosseau [1986], and Haghani [1989]. In PAGE 139 128 each study, the problem is broken into two co mponents: a trainscheduling problem and a shipment routing problem. An iterative procedure is used to solve each of the two parts of the problem in succession, using the solution from the other part to guide the next iteration. Each author relies on the separa bility of the train scheduling problem and shipment routing problem. None of the mode ls are meant to produce a detailed schedule for the railroad; they produce only train fre quencies for a represen tative days schedule which provides a rough framewor k for creating a schedule. Keaton [1992a, 1992b] uses Lagrangian re laxation to simplify the problem. He incorporates the train capacity, travel time and demand flow constraints into the objective with a Lagrangian multiplier. Relaxing the constraints allows him to decompose the problem into separable train and shipment routing problems. By relaxing the train capacity constraint in his model, the ship ment routing problem may be viewed as a collection of shortest path problems. Using a dual adjustment approach, the author arrives at an infeasible lower bound to the problem, then uses a simp le heuristic to arrive at a feasible solution. Haghani[1989] incorporates empty equipment flows into his model, which are ignored by the other studies. He schedules them ag ainst variable daily demands but uses a fixed daily trai n frequency. The other two models schedule against a representative days demands. He applies his model to a smaller ne twork of eight nodes and ten twoway arcs. His heuristic decomposition technique depends on the special structure of his formulation and becomes unt enable for larger problems. Crainic and Rosseau [1986] present a genera l model of multimode freight transportation, and Crainic et al. [1984] examine the spec ific case of freight rail. They solve the problem through decomposition and column generation techniques A modified shortest path algorithm in PAGE 140 129 which the capacity constraints are integrated in to the objective with a penalty term is used to find the best route of th e shipments given a schedule. The information from the demand flow is used to generate a new tr ain schedule. Gorman [1998] applies metaheuristics: tabu search and ge netic algorithm to solve an integrated train scheduling problem. The strength of his approach lies in considering weekly train schedule and in modeling the problem with broader scope. Our work in this dissertation, consider the problem instance in totality, as against the past research on solving a subset of problems. We have developed a decomposition approach to solve the problem. However, our approach does not separate the train scheduling and shipment routing. We have developed heuristic algorithms to solve each of the decomposed problems. The strength of our approach is in decomposing the problem in sequential phases and developing e fficient algorithms to solve the problems in each stage. Rest of the chapter is organized as follows In section 6.2, we describe the train schedule design problem in more details. In section 6.3, we present the decomposition strategy. Sections 6.4 and 6.5 give the soluti on approaches to solve the problems in Phase I and Phase II respectively. In section 6.6 we present the computati onal experience with a reallife instance and in section 6.7, conclusions and scope for further research are given. Also, in rest of the paper we invariably us e train schedule design pr oblem as train design problem. 6.2 Problem Description As discussed in Section 6.1, the train de sign problem is closely related to other railroad scheduling problems. This problem takes input from the solutions of other PAGE 141 130 scheduling problems and produces a solution, which optimizes the objectives and honors many constraints, which are essential for th e feasible execution of the schedule. The constraints set of this problem includes many historical and contractual requirements, which are very difficult to be represented mathem atically. In this section we have tried to include as many constraints as possible. We describe below the inputs, constraints, objective functions and decision variab les of the train design problem. 6.2.1 Input The train design problem is to design and to schedule trains to transport given shipments using existing physical infrastructure like, physical networ k, locomotives, etc. As described in section 6.1, generally blocking plans are deve loped first, which consolidate the shipments into blocks and th e demand for the train design problem can be specified in terms of transpor tation of blocks. Therefore, th e essential inputs for the train design problem are the existi ng blocking plan and the exis ting physical network, which includes the stations and links in the network. Apart from these, we also need input on the capacities available in network, costs of tr ansportation and use of infrastructure, and matrices to measure the goodness of the solu tion. Also, very often matrices are in different units of measurement, we need so me conversion factors to weigh them against each other. 6.2.2 Objective Function In the train design problem, there are main ly two types of objectives: one which minimizes the transportation cost and transit time of shipme nts, and another which helps in smooth operation of the train schedule. We present below these objectives and also discuss their importance in the context of overall objective PAGE 142 131 A big part of railroad expenditure goes in transporting shipments from their origins to their respective destinations. This cost consists of expenditure in operating the trains, maintaining the trains, managing gr ound operations, fuel cost and cost of capital utilized. The transportation cost of shipments greatly affects the ground level economic performance of the railroad and hence minimizing the transportation cost is the most important objective of the railroad. Railroads are facing stiff competition from the trucking industry and the greatest advantage of the transportati on by the truck over that by rail is less transit time. From the supply chain perspective and cust omer satisfaction, transit time is an important factor. Apart from this, costs proportional to the transit time are also incurred in terms of the capital of the customer stuck up dur ing transportation. Therefore, one of the important objectiv es of the train de sign problem is to minimize the transit time of shipments. Whenever, a train is engaged in the trans portation of shipments, railroads have to incur some cost of infrastructure a nd crews. Therefore, to have economic advantage, railroads want to utilize minimum number of trains and this is one of the important objective of the train design problem. As discussed earlier, the so lution of the train design pr oblem should not warrant for additional physical infrastructure. It implies that solution should be such that it does not cause congestion and infeas ibility at the stations or on the links. To achieve this, one of the objectives of the train design problem is to minimize the number of stops of a train as well as to minimize the number of trains stopping at a station. For the smooth operation of the trains, it is highly desired by the railroads that there should be consistency in the operation of trains. For example, if a block is made everyday of the week, then the railroad will like to assign these blocks to same train sequence everyday. Therefore, one of the objectives of the trai n design problem is to minimize the inconsistency in blocktotrain assignment. It can be noted that some of these objectives may be rephrased and can be considered as the constraints. However, in cluding them in objective function makes the mathematical formulation easier to solve. It ca n also be noted that some of objectives are conflicting in nature. For example, to minimize the transit time we can create direct trains from the origin to the destination of each bloc k. However, this will adversely affect the objectives of minimum number of trains and minimum trans portation cost. Therefore, a tradeoff is needed to balance these objectiv es and input is requi red to prioritize the objectives. PAGE 143 132 6.2.3 Decision Variables The train design problem is a complex di screte optimization problem in which many decisions are required to be made. We present below the decision variables for this problem. As we are considering zerobased train design problem here, decisions are required to be made with respect to the trains to be made for the trans portation of shipment. Then for each train, decisions are required to be made with resp ect to the origin and the destination of the train, physical rout e of the train, stops on the route of the train, frequency of the train and its time schedule. As the trains are designed to carry the shipments from their origins to their respective destinations, decisions are to be made with respect to train rou for each shipment. However, as the shipments are already consolidated into blocks, the decision variable for the trai n design problem is to find train route for each block. This problem is also called the blocktot rain assignment problem and is studied in previous chapter, when we determine the train route for blocks using given train schedule. 6.2.4 Constraints The constraints set for the tr ain design problem consists of logical requirements for the feasible execution of the train schedul e and practical require ments for the smooth implementation of the solution. We next describe these two types of constraints separately. Logical requirements Each train must have a valid route and valid schedule, i.e., two consecutive stops of the train must be physically connected and timings at the stops must be logical and feasible. Each block must be assigned to a valid se quence of trains. For example, if a block changes train at a station, then this sta tion must be valid stops for both the trains. Practical Requirements The load on a train should not exceed th e capacity restriction. The capacity restrictions on a train are given in terms of maximum number of cars in the train, maximum weight of the trai n, and maximum length of th e train. Depending on the case, railroads may need all or some of these restrictions im posed while assigning PAGE 144 133 blocks to the trains. Also, the capacity of a train may vary on different parts of its route. If prespecified, a block should not flow on particular links. Sometime, this restriction is imposed to prevent the fl ow of special kind of shipments (like hazardous material) on certain route. If prespecified, a block mu st flow on a particular rout e. This restriction may be imposed to meet the specific customer requ irement or to meet legal requirements. Apart from above, there may be many additional requirements which must be honored. For example, certain trains must be made, a train must follow given timings, a block must be swapped at a particular yard, etc. In this chapter, we have tried to develop algorithms, in which these requireme nts can be easily incorporated. We now discuss the size of a typical US railroad train design problem in terms of number of decision variables and number of constraints. Let us suppose that there are 6000 stations and 12000 links in the physical network and around 1300 blocks are made everyday. Now let us assume that we have a upper limit on number of trains to be formed (say 500). Then for each of these trains, there are 6000 possible origin and 6000 possible destination. Also, for each train we will have to select a proper set of links to form the route of trains. Once the train route is decide d, we have 24 X 4 = 96 options for each train as starting time at origin. We assume here th at once the starting time of a train is known, we can determine its timing on the route. Need less to say that altogether the number of decision variables will be in millions. It can be recalled that another important decision variable for the train design problem is to fi nd the blocktotrain assignment. It can be recalled from Chapter 5, that once the tr ain schedule is known, number of decision variables for the blocktotrain assignme nt problem was around 12 millions. Now when we are solving an integrated problem, thes e decision variables will be much more. The reallife train schedule of a ra ilroad indicates that the numbe r of stops on the route of a PAGE 145 134 train varies from 2 to 200. If we treat th e block to train assignment problem as multicommodity flow problem, then we need mass constraint for each of the block at each station (around 1300 X 6000). Apart from this, we also need the logical constraint for each of the train at each stop on its route. As we do not know the train route in advance, the number of flow constraints for trains will be around 500 X 6000. Besides these, we will have constraints for cap acity requirements and logical timings on the route. Considering altogether, the number of constrai nts will be in millions. With normally used computational facility, it is extremely difficult to solve any mathematical formulation of the train design problem to optimality. Therefore, in this chapter, we are only concentrating on development of heuristics. 6.3 Problem Decomposition As described in previous section, it is extremely difficult to solve the train design problem in totality using normally used comput ational resources. Therefore, it is prudent to develop some algorithms, which can produce nearly optimal solution. However, as there are millions of decision variables involve d, it is again very difficult to solve the train design problem in totality to obtain ne arly optimal solution. To keep the problem size manageable, we have developed an appr oach, which decomposes the train design problem in three stages. In our decompositi on technique, the problem in Phase I is a network design problem, that in Phase II is a network flow problem and problem in Phase III is a scheduling problem. We next brie fly describe the problem in each phase. 6.3.1 Phase I It can be recalled that one of the input to the train design problem is the blocking plan and decisions are to be made with respec t to assigning blocks to train sequences so PAGE 146 135 as to minimize the transportation cost. In a train sequence, a block is picked by a train at a station and then dropped at another station, where it is picked up by another train and so on until the block reaches its de stination. Now, if we call the train path between two consecutive stop as the train segment, then in a train sequence, a block rides on one or more train segments of a train. If we treat all the train segments of same type irrespective of the trains to which they are associated then finding a train sequence of a block is equivalent to finding the sequence of train se gments. In view of train segment concept, we can define the train design problem as to designing the train segments, routing the blocks on train segments and then combining th e segments to create the trains. Now, in the train design problem, there is time associ ated with each train and hence there should be time component associated with each segmen t as well. However, in this formulation, as we do not know the train of which a tr ain segment is part of, deciding the time component of a segment is impossible. To simplify this, we do not decide the time component of a train segment in this phase Once we know of which train a segment is part of, the decisions regardi ng its timing are taken. Thus th e problem in this phase is defined as to design the train segments and to route the blocks on a network consisting of train segments. The objectives in this phase ar e to minimize the transportation cost and to minimize the number of train segmen ts terminating at a station. 6.3.2 Phase II From Phase I, we get the train segments to be made and the route for each block. At this point we may combine train segments to create the trains. When combined, the end points of train segments denote the stops of a train. If we construct a network consisting of train segments then creating trains are equi valent of creating paths in the network such that all arcs are covered and each arc can be the member of only one path. Now, flow of PAGE 147 136 blocks can be represented as the flow on path s. Whenever, a block flows from one path to another, there is change in trai n and any transition in trains results into increase in transit time and transportation cost. Therefore, one of the objectives in this phase is to construct the trains in such a way that block swappi ng (traintransition) can be minimized. In a train we also want to keep the consistency in train load, i.e., a train should carry same load on different segments as same set of locomotives are assigned throughout its itinerary and there will be underutilization of power it pulls less load on part of its itinerary. Again as in Phase I, we do not determine the time component of trains in this phase. 6.3.3 Phase III From Phase I and Phase II, we know the tr ains, their physical route and blocktotrain assignment. In Phase III, we take rest of the decisions for the train design problem. This includes the time schedule of trains and their frequency. The obj ectives in this phase are to minimize the transit time of shipment s and maximize the consistency in blocktotrain assignment problem. The constraints in this phase are to honor the capacity restrictions and many othe r practical restrictions. It can be noticed that by using decom position technique, we can have better formulation of the problem in each stage and as described in remainder of the chapter, we can develop time efficient algorithm for the probl ems in different phases. It can be also noted that performance of an algorithm in a phase greatly affect s the performance in subsequent phases. For example, a good design of train segments is very important to construct good trains and a good set of trains is very important for scheduling. To improve the performance of an algorithm in a phase, we relax some of the constraints and impose them in subsequent phases. For exam ple, we relax the capacity constraints in PAGE 148 137 Phase I and Phase II and impose them in Phas e III. Now if there are capacity violations from Phase I and Phase II, we can address them by splitting the trains or by adjusting the frequency of the train. We now briefly desc ribe the simplifications in each phase to achieve efficiency. It can be recalled from Section 6.1, that generally a train schedule is repeated every week on the planning horizon. It implies that the solution of the train design problem should be weekly train schedule. However, if we consider we ekly time horizon in Phase I and Phase II, then it may result into different train route for a bloc k on different days of week, which is undesirable for the railroads as they want to maintain consistency in blocktotrain assignment. Therefore, we c onsider only the daily ve rsion of the problem in Phase I and Phase II, which produces the tr ains with the assumpti on that they will run everyday of week. It implies that before appl ying Phase III, we assume that frequency of each train is seven, which is adjusted by the algorithms developed in that phase. Also, in Phase I and Phase II, we do not impose the capa city restrictions on trains, as the loads on train are not the correct repres entation of train load on weekly time horizon. Also, we can always split up a train if there is too much lo ad on it and conversely, combine two trains in one, if loads on trains are too low. Relaxi ng the capacity constraints on trains in Phase I and Phase II allowed us to develop efficient algorithm. Now, we will briefly discuss the limita tions in developing algorithms for each phase. In Phase I and II, i nputs required are the blocking plan, physical network, limitations on number of trains and number of train segments at a station. These inputs are easy to obtain from the railroads and they do not have to make ex tra effort to generate the data. However, in Phase III, we need mu ch more data regarding the timings of the PAGE 149 138 trains, travel time, crew policy, maintenance requirements, capacity restrictions on each train segments, etc. Needless to say, high leve ls of details are required for these data and deep coordination with railroads is required to develop the algorithm in Phase III. There are many more practical requirements are there, which must be incorporated in the final solution and again mutual work is required to be done with railroads to precisely and quantitatively define these additional constrai nts. Therefore, in this chapter, we only address the problems in Phase I and Phase II and keep the development of algorithms in Phase III as the future work. In remainder of this chapter, we will only discuss the train design problem in scope of problems defined in Phase I and II. In section 4, we describe the solution approach to solve the decom posed problems in Phase I and Phase II. 6.4 Phase I: Designing Train Segments As described in previous section, the problem in Phase I is defined as to create the train segments and to route the blocks on the networ k created. The objectives in this phase are to minimize the transportation cost and to mi nimize the train segments emanating from a station. The constraints for this problem ar e the maximum number of train segments that can be created at a node and valid route for each block. We call the problem in this phase as Train Segment Design (TSD) Problem. The TSD problem is analogous to the blocking problem described in previous chapter. In bl ocking problem, the deci sion variables are to create blockarcs and then to route the ship ment over the network consisting of stations and blockarcs so as to minimize distance tr aveled and number of intermediate handlings of blocks. In the blocking problem, the constr aints are on the number of blockarcs that can be created at a station. Liu [2003] describes the blocking problem and solution PAGE 150 139 approaches to solve the problem. Our approach to solve Phase I problem is very similar to that of Liu [2003]. Mathematically, the TSD problem is designing a network, called a train segment network, and routing blocks ove r this network so as to optimize the total shipment cost. Figure 6.1 gives a sample train segment netw ork, where there are three types of nodes: origins (where blocks originate), yards (where blocks are swapped), and destinations (where blocks terminate). We show here a si mplified network as in practice yards can be origins as well as destinations, and nodes can send as well as receive shipments. Each arc in the network represents a segment. At the tail of the arc, the train can pick up the blocks and at the head of the arc, it can drop the bl ocks. This signifies that no activity by the train will be performed at the intermed iate stations on a train segment. Figure 6.1 An example of a train segment network. O O r r i i g g i i n n s s D D e e s s t t i i n n a a t t i i o o n n s s Y Y a a r r d d s s 3 1 2 4 5 7 8 9 10 11 12 6 S1 S2 S3 S4 T1 T2T3 T4 PAGE 151 140 Blocks travel over the train segment ne twork. We show four such blocks S1, S2, S3, and S4 which originate at the nodes 1, 2, 3, and 4, respectively, and are destined for the nodes 9, 10, 11, and 12, respectively. A block goes from its origin to some yard, then goes from one yard to another yard one or more number of times, and then goes to its destination. A block is picked up at its or igin and subsequently swapped each time it changes train. For example, block S2 can follow two paths 2510 or 25810, or 26810 from its origin to its destin ation. For the path 2510, it is picked up at node 2 and swapped once at node 5, if the train on segment 510 is different than that on 25. For the path 25810, it is picked up at node 2 and swapped twice at the nodes 5 and 8, if trains on different segments are different. We refer to a path of a block over the train segment network as a train segment path. The cost of routing a block is the sum of the cost of flow over the train segment network and the number of swapping it goes th rough. The cost of flow on a segment is linearly proportional to the length of the arc (i, j) in miles, where the length of a arc (i, j) can be considered to be the shortest rout e in the physical railr oad network from node i to node j. If we assume that a block is swapped every time it changes train segment, each car is picked up at its origin and is swapped at each intermediate node of the train path it visits. It can be note d that number of swappi ng thus calculated is an upper bound on the actual number of swapping as if two segments ar e part of same train then there will be no swapping. Since the cost of picking up of each block at its origin is independent of the train segment path, we do not consider th e picking up cost bu t consider only the swapping costs in our model. Each arc (i, j) in the train segment network has an associated cost cij denoting the cost of flow per ca r; this cost might depend upon the PAGE 152 141 distance between the nodes i to node j in the physical railroad. Further, each node i in the train segment network has an associated swapping cost di. In terms of these notations, the costs of the train segment paths 259, 2589, or 2689 ar e as given below. Table 6.1 Illustrating the cost of flow in the blocking network. Path Cost of path per unit flow Cost of swapping per unit flow 2510 c25 + c5,10 d5 25810 c25 + c58 + c8,10 d5 + d8 26810 c25 + c68 + c8,10 d6 + d8 We determine the transportation cost of each block by adding two cost terms: (i) the cost of each segment path multiplied by th e number of cars in and the average cost per mile per shipment, and (ii) the cost of swapping of a shipment multiplied by the number of cars in it and the average handling cost per swapping. The TSD problem is to design the train segment network and route all the blocks over the train segment network so that the total cost of transportation is minimum. The TSD problem can be formulated as an Integer Programming problem, with two kinds of binary decision variables: N: a set of nodes denoting the stations wh ere blocks originate, terminate or are swapped. We will use the index i to denote a station in the node set. A: a set of potential train segment arcs in NN, i.e., (i, j) A if a segment can be built from node i to node j. G = (N, A): the blocking network. +(i): the set of arcs in A emanating from i. PAGE 153 142 (i): the set of arcs in A entering node i. K: a set of blocks. We will use index k to denote a block. o(k): the origin of block k K. d(k): the destination of block k K. vk: the number of cars in block k K. mij: the cost of flow per car on (i, j) A. In general, this cost is proportional to the length of the block. bi : the maximum number of train segments that can be made at node i N. hi: the cost of swapping a car at node i N. cij: the cost of creating a train segment (i, j) A. The TSD problem requires that (i) each block k K consisting of vk cars must be sent along a single path in the train segmen t network; and (ii) the number of train segment arcs emanating from node i be at most bi. The objective f unction in the TSD problem is to minimize the weighted sum of flow cost and swapping cost. The TSD problem has two sets of binary decision variables: yij and k ij x The variable yij takes value 1 if we decide to build the arc (i, j) A, and is 0 otherwise. The decision variable k ij x is 1 if block k flows on the arc (i, j) A, and is 0 otherwise. The PAGE 154 143 TSD problem can be formulated as the following mathematical programming problem: (,)(,)()minkk ijijij kKijAiNkKiji imxhx 6.1a subject to (,)(,)() 0()(), ()iik kk ijij ijij kvifiok x xifiokordk vifidk for all k K, 6.1b ,k ijijij kK x uy for all (i, j) A, 6.1c (,)(),iji ijiyb for all i N, 6.1d yij = 0 or 1 for each (i, j) A, and k ij x = 0 or vk for all (i, j) A and all k K. 6.1e In the above formulation, the constraint 6.1b in conjunction with the constraint 6.1f ensures that a block flows on a single path in the train segment network, and the constraint 6.1c ensures that block can flow on an arc only if it is The constraint 6.1d restricts the number of bl ocks created at a node. Even the TSD problem is very largescal e problem and it is not possible to find an optimal or a nearoptimal optimal solution for this problem using currently available commerciallevel software and computing power. Hence, we focused our effort on developing heuristic algorithms for this problem that do not guarantee optimality of the solution but find a nearly optimal solution within a reasonable computational time. Neighborhood search algorithms ar e perhaps the most effectiv e heuristic algorithms to solve difficult combinatorial optimizati on problems and, hence, we develop a neighborhood search algorithm to solve the TSD problem. Figure 6.2 describes our VLSN s earch algorithm for TSD problem. PAGE 155 144 algorithm VLSNTSD; begin construct an initial train segm ent network and send all blocks along shortest paths in the train segment network; while the current solution is not locally optimal do for each node i N do {one pass} reoptimize the segments emanating from node i; send all blocks along shortest paths in the updated train segment network; end; {one pass} end; Figure 6.2 The VLSN search algorithm for the TSD problem. This algorithm is a neighborhood search al gorithm. It starts with a feasible solution of the TSD problem and iteratively improves the current solution by replacing it by its neighbor until the cu rrent cannot be improved. The neighbor of a solution is defined as all the solutions we can obtain by changing the segments at one of the nodes in the train segment network. This neighborhood search algorithm is a VLSN search algorithm because the number of neighbors is very large. First, there are n ways to select the node for which we change the segments emanating from it. Secondly, the number of ways we can select k segments out of a node with p arcs emanating from it is pCk which grows exponentially with p and k. This algorithm has two important subroutines: to construct the initial feasible solution; and to reoptimize the train segmen ts emanating from a node. We now describe these steps in greater detail. 6.4.1 Constructing an Initial Feasible Solution We use a fairly simple construction heur istic to create the in itial solution of the train segment network. A solution of the trai n segment network must ensure that for each block its origin node is connected to its de stination node through a tr ain segment path. To ensure this, we first construct a directed cycle passing through all the yards exactly once PAGE 156 145 and returning to the first yard (that is, a Hamiltonian cycle). Then, we create a train segment from each origin node (that is, wher e some block starts) to the nearest yard. Next, we connect each destination node (that is, where some block terminates) to the nearest yard, which still has some spare ar c emanating capacity. Since each origin is connected to some yard, and each destination is connected to some yard, and each yard connected to every other yard (through the Ha miltonian cycle), this construction process will ensure that each origin node is connected to each de stination node. After a train segment network has been constructed, we rout e blocks along their sh ortest paths in the train segment network. We may point that our algorithm to constr uct the initial solution is a fairly simple algorithm and does not do any optimization. We have found in our computational testing that the final solution obtained by our algor ithm is fairly insensitive to the starting solution and improving the initial solution does not have much impact on the final solution. We therefore did not spend much effort in tryi ng more sophisticated methods for constructing the initial solution. 6.4.2 Finding an Improved Neighbor As described earlier, we define the nei ghborhood of a solution as all the solutions that can be obtained by changing the segments made at one node only and rerouting all shipments along their updated shortest path s. Our algorithm performs passes over the nodes (all origins and yards where segments are built) one by one by selecting a node (say, node i) and reoptimizes the segments made at that node assuming that the train segments at other nodes do not change. To reoptimize the segments made at node i, we have developed an algorithm, called maximum saving algorithm. The maximum saving algorithm performs the following steps: PAGE 157 146 It removes the train segmen ts already made at node i and reroutes the blocks passing through this node. It then builds the segments at node i one by one using the following maximum savings rule. For every potential train segmen t arc that can be built at node i, it computes the savings in the total cost if th at train segment arc is built and all blocks are rerouted to take advant age of this new train segment arc, selects the train segment arc with the maximum savings, build s it and reroutes all the blocks for which costs would decrease by using the newly built train segment arc. We repeat this process until all the required trai n segments at this node are built. We point out that it is necessary to perform several passes over the nodes since when we reoptimize the segments made at node i, we assume some segments at node i+1, i+2, , n. But later when train segments at nodes i+1, i+2, , n are reoptimized (and possibly changed), the segments made at node i may need to be changed in which case segments built at node i need to be reoptimized agai n. Thus, the algorithm performs passes over the nodes 1, 2, , n, in order, in the train segm ent network and reoptimizes the train segments built at those nodes. The algorithm terminates when the solution converges, that is, in an entire pass the objective function value does not change significantly. We observed in our computational results that the solution converges within 10 passes over the nodes. 6.5 Phase II: Forming Trains Form the phase I, we get a train segment ne twork with the routings of blocks. In the train segment network, trains run on the segm ents with stops at the end points of train segment arc. For the test instance in our cas e, there are around 1100 segments created in phase I. As the train design problem is to de sign trains and their schedule, one straight forward solution for the train problem is to create train for each of the train segment. However, obviously this solution will not be viable economically and practically. This PAGE 158 147 option also means that there will be as many block swaps for a block as many train segments are there on its route. To make th e solution of the train design problem more efficient, train segments are combined to cr eate the trains. The benefits of combining trains segments are in terms of savings lo comotives and crew requirements, and reduced number of block swaps. We call the problem of combining the train segments to form the trains as Train Formation Problem (TFP). As discussed in previous section, we achieve an upper bound on the number of block swaps in the solution of the train desi gn problem. In this phase our objective is to combine segments in such a way that there will be minimum number of blocks swaps. For example, in the train segment network give n in Figure 6.1, let a block takes the train segment path 26810. Now if we create train for each segment separately, there will be two block swaps at stations 6 and 8. Now if we combine the segments (2, 6), (6, 8) and (8, 10) in a train, then there will be no block swap as this block will be carried by only one train from its origin to its destination. Now as a segment may carry more than one block and they may be having different tr ain segment paths, the number of blocks depends on the way we combine the segments to form trains. For example, let there are three blocks (say k1, k2 and k3), which are routed on paths 159, 159, and 15810 respectively. Now if we decide to create the trains 159 (say T1) and 5810 (say T2), then there will be one block swap, as block k3 will be transferred from train T1 to train T2 at station 5. On the other hand, if we decide to create the trains 15810 (say T3) and 59 (say T4), then there will be two block swaps as the blocks k1 and k2 will change train at station 5. Therefore, the number of block swaps depends on the way train segments are PAGE 159 148 combined to form trains. One of objectives of the train formation problem is to minimize the number of block swaps. Generally, a set of locomotives is assigned to a train, which carries it from its origin to its destination. To maximize the locomotiv e utilization, another important objective of the train design problem is to maintain unifo rmity of load in a train. For example, suppose we form a train consisting of segments 15 and 59. Now suppose that there are 30 cars on the segment 15 and 60 cars on the segment 59 then the train must be assigned with locomotives which can pull 60 cars. However, this will result into underutilization of trac tion power on segment 15. For the railroads, locomotives are very vital resource and proper utilization of them is very important. The train formation problem can be defined as path covering problem in which we need to create minimum number of paths in the train segmen t network such that all the arcs are covered at minimum cost. The cost of creating a path consists of fixed charge for a train and interrelation among train segments in the path. The interrelation among train segments implies the common blocks flowi ng on these segments and car load on these segments. If there are more common blocks on the segments in a path, there is less number of block swaps as blocks are riding sa me train for a larger part of its itinerary. Similarly more uniformity of car load implies better utilization of locomotives attached to the train. In the train segment network, a valid train can also be formed by creating a directed cycle, which implies th at a train starts from a statio n and terminates at the same stop after visiting other stations in the network. In this section we will denote such cycles as paths, which starts and terminates at same station. PAGE 160 149 We formulate the train formation problem as a minimum cost flow problem on a suitably constructed graph (say G). In this graph, we capture th e block swapping cost and penalty on the nonuniformity in load on a train as the arc cost. We create the network G as follows (refer Figure 6.3). Figure 6.3 Graph for minimum cost flow formulation. We create an arc for each of the train se gment in Figure 6.1 and create nodes at the tail and head of these arcs. We ca ll the set of nodes thus created as N and set of arcs thus created as train arcs T. Let there is an arc a in G corresponding to an arc a in train segment network. Let Aa + and Aa are the sets of arcs emanating from the head node of the arc a and entering into th e tail node of arc a in the train segment network. Now, for each of the arc a in G we create arcs from th e head of this arc to the tail node of all the arcs which represents the arcs in Aa + in the train segment netw ork. We denote the set of all such arcs as connecting arcs C Let there are m arcs in train segment network, then as A' B' C' D' E' F' A" B" C" D" E" F" s t PAGE 161 150 there is a connecting arc for each pair of in coming and outgoing arcs at a node, there can be O(m2) connecting arcs in G Also, the number of connecti ng arcs increas es with the density of the train segment network. However, in practice, the train segment network is very sparse and hence the number of connecting arcs is linear to the number of train arcs. Next, we create a dummy source node and a dummy sink node in the network G We create arc from source node to the tail of each of the train arc and similarly, create arc from the head of each of the train arc to the dummy sink node. We also create an arc from the dummy source node to the dummy sink node. Thus there will be m + 1 arcs from the source node and there will be m + 1 arcs to the sink nodes If we add the number of different type of arcs, it will be of the O(m2). As we create two nodes in G for every train arc, there will be 2m + 2 nodes in the network G Figure 6.3 represents the graph G We now set the capacity a nd cost of the arcs in G and set the supply/demand at the nodes to formulate the train formation problem as minimum cost flow problem. We set the upper bound of 1 and lower bound of 0 on the arcs from source node to the tail nodes of train arcs. We also set 0 cost of flow on these arcs. The upper and lower bounds on the flow on train arcs are set to 1. This implies a compulsory unit flow on train arcs. The cost of flow on train arcs are also set to 0. We do not set any capacity limit on the flow on connecting arcs and in the following paragraph, we discuss the cost of flow on these arcs. Similar to the arcs from the source node, we set the upper bound and lower bound on flow to 1 and 0 respectively and set the cost of flow to zero. We do not put and capacity restriction on flow on arc between the source node and the sink node and set the cost of flow to 0. We put a supply of m units at the source node and put a demand of m units at the sink node. PAGE 162 151 Now, we describe the cost of flow on the connecting arcs in the network G Any flow on these arcs denotes that the train arcs at the tail node and at the head node are associated to the same train. Therefore, for every possibility of combining two train segments at a node, there is a connecting arc in the network G Let us assume that a1 and a2 are the respective train arcs at the tail node and head node of a connecting arc. Any block which is flowing on arc a1 and not flowing on a2 and not terminating at the head node of the arc a1 in the train segment network must go for block swap as it needs to take some another train for onward itinerary. As we started with an upper bound on number of block swaps from Phase I, we calculate the savings in block swaps if segments a1 and a2 are combined into same train and associat e with the correspondi ng connecting arc in G It can be noted that this cost is independent of the other tr ain segments in a train as a block swap affects only at th e connecting terminal of tw o segments. Similarly, we calculate the penalty in terms of difference in number of cars in segments a1 and a2 and add to the cost of flow on corresponding connecting arc. In case of practical implementation, railroad may also specify their likings for combining certain pair of segments. This may be required to improve the customer satisfaction. Theorem: There exists a onetoone correspon dence between the solution of the train formation problem and the solution of minimum cost flow problem defined on network G Proof: Any solution of the minimum cost flow problem can be represented by the flow along the paths from the source node to the sink nodes and by the flow along the cycles in the network G We claim that each of these pa ths and cycles represent a valid train and all the train segments in the trai n segment network are covered by a train in G PAGE 163 152 As there must be positive flow on each train arc, a train arc must have been covered by some path or cycle and as the flow is of un it quantity, the train arcs must be covered by single path or cycle. Now in any path or cy cle in the solution, there must be connecting arc between two consecutive train arcs, which im plies that there must be common station, which connects two segments. Thus every path and cycle with positive flow represents a valid train for the train forma tion problem. In case of paths, the origins and destinations of the trains are given by the tail node of the first train arc and head node of the last train arc respectively. In case of cyclic flow, a ny node on the cycle can be made as the origin and destination of the train. We have incorporated the objectives of th e train formation problem into the flow cost on connecting arcs in network G the optimal solution of the minimum cost flow formulation gives the optimal solution of th e train formation probl em. The minimum cost flow problem can be solved very efficiently to optimality i.e. we can get the optimal solution of the train formation very efficiently. 6.6 Computational Results We implemented the Phase I and Phase II of our algorithm on a real life instance of a major US railroad. The railroad is engaged in carrying around 50,000 shipments per day. They build around 1,300 blocks per da y to carry these shipments. The physical network of railroad consists of around 6,000 stations a nd 12,000 links. We implemented the algorithms on Visual Studio .Net plat form using C++ programming language. The testing of the algorithms is done on a PC with 1 GB RAM and 3 GHz of processor speed. In Phase I, we need the input on the number of train segmen ts that can be created at a station. We get this input from the ex isting train schedule of the railroad. We PAGE 164 153 considered only those stops of a train where block swapping is permissible. To minimize the number of train segments, we put a pena lty on the creation of an arc in the train segment network. Our algorithm produces a solution with 1,116 train segments, while honoring the constraints on the number of train segments. The computational time required to produce the solution is within 5 seconds. When compared with the existing train schedule, there is a major saving in th e number of train segments, which is around 2,000. We take the solution of Phase I as the i nput to Phase II. We created a graph to formulate the problem in this phase as mini mum cost flow problem. The graph in this phase consists of 10,456 arcs and 2,234 nodes. It can be noted that the number of arcs is much less than the upper bound of O(m2 = 1.2 million). We used CPLEX optimizer to solve the minimum cost flow problem. The co mputational time take n to obtain optimal solution is around 5 seconds and produces a so lution with 325 trains. Out of these, 66 trains starts and ends at same station. In the actual train schedule, there are 339 train schedules to carry same shipments. Therefor e, our solution gives a saving of 14 trains, which may result in substantial economical saving for the railroad. 6.7 Conclusion and Further Scope In this chapter we have tried a decom position approach to solve the train design problem. Our solution approach consists of three phases. We have developed algorithms to solve the problems in Phase I and Phase II. Further work with the coordination of a railroad is required to be done to develop algorithms for Phase III. PAGE 165 154 CHAPTER 7 FURTHER RESEARCH AND CONCLUDING REMARKS In this dissertation, we have tried to solve few hard combinatorial optimization problems using very largescale neighborhood (VLS N) search technique. Our effort is to use efficient network flow techniques in deve loping heuristics to solve these problems. In this chapter, we will first summarize the c ontribution and then will present the scope of future research. 7.1 Contribution Summary Prior to our application of VLSN search heuristic to the quadratic assignment problem, this technique was only applied to partitioning problems. In chapter 2, we present the application of VL SN search technique to the quadratic assignment problem, which is a nonpartitioning optimization probl em. We have developed an improvement graph based approach, in which we approximate the cost of multiexchange with the cost of a cycle in the improvement graph. We have also developed an efficient mechanism to update the improvement graph, whenever, we make some change in the current solution of the problem. We performed extensive computational experiment on the problems available in QAPLIB and our algorithm performs well all across the test instances. In chapter 3, we studied a classical and important defense related optimization problem, called the weapontarg et assignment problem. This problem is to assign given set of weapons to given targets such that lo ss to the enemy can be maximized. This is a PAGE 166 155 nonlinear optimization problem and previous efforts are limited to solve only small size problems. In this dissertation, we suggest linear programming, inte ger programming, and network flow based lower bounding methods us ing which we obtain several branch and bound algorithms for the WTA problem. We ha ve developed a network flow based heuristic and a VLSN search heuristic. The co mbination of network flow based heuristic and VLSN search produces exceptionally good quality solution in very manageable computational time. We have been bale to solve the problems of the size, which were historically intractable. In chapter 4 of the dissertation we have developed algorithm s to solve another combinatorial optimization problem called th e university coursetim etabling problem. In this problem we fix the meeting time of c ourses in a university with the objectives of maximizing the teachers preferences and mini mizing the conflict in students interests. Most of the literatures on this problem are case specific, where algorithms are developed to solve a given test instance. We have trie d to develop a generic model for the coursetimetabling problem. We have developed algo rithms in conjunction with tabu search, which performs very well on the randomly generated test instances. Chapter 5 and Chapter 6 of this disserta tion deal with two reallife railroad scheduling problems. In chapter 5, we have developed algorithms to solve the blocktotrain assignment problem. This problem is to find train path for a block such that transportation cost and transit time are mini mum, and there is no capacity violation. We formulate this problem as an integer program on a spacetime network. We have developed very largescale neighborhood search based heuristics to enumerate the best train paths for each block. We have devel oped numerous approaches to solve this PAGE 167 156 problem, when possible train pa ths are given for each block. In chapter 6, we study a broader railroadscheduling problem called the train schedule design problem. This problem is to design trains a nd their schedules such that ove rall transportation cost can be minimized. This is very large scale and co mplex optimization problem. Enormous efforts have been made in the past to solve this problem. We have developed a decomposition based solution approach, in which the tr ain schedule design problem is solved sequentially in phases. We have developed a VLSN search heuristic to solve the problem in Phase I and have formulated the problem in Phase II as minimum co st flow problem on a suitably constructed graph. Our algorith ms produce very good quality solution when tested on a reallife instance. 7.2 Future Scope Our algorithm to solve the qua dratic assignment problem is an enhanced version of local search heuristic. We expect that wh en applied in conjunction with other metaheuristic, the VLSN search algorithm may perfo rm better. Further work is required to be done to develop algorithms in th is regard. In chapter 3, we have considered the static version of the weapontarget assignment probl em. However, in reallife application the problem is of dynamic nature. Our algorithms in chapter 3 can be easily modified to solve the dynamic version of the problem. We expect further work to modify our algorithms to accommodate dynamic nature of the problem. With respect to th e coursetimetabling problem, further work can be done in developi ng better algorithms to generate random test instances and to incorpor ate any additional constraint in the VLSN search algorithms. In chapter 5, we have consid ered daily version of the bloc ktotrain assignment problem. In this case, we assume that a block is made everyday and a train runs each day of the PAGE 168 157 week. However, in general train schedules are on weekly basis and a train may not run everyday. Further work is required to devel op proper postprocessing routines to extend our approach to the seven day blocktotrain a ssignment problem. In chapter 6, the future work is required to be done to develop the algorithm in coordination with railroads to solve the problem in Phase III. We also expect a thorough te sting of algorithms developed by experimenting with the reall ife data of more than one railroad. PAGE 169 158 APPENDIX: COMPUTATIONAL EX PERIMENT ON QAP INSTANCES PAGE 170 159 Table A.1 Experimental Resu lts for Symmetric Instances 2OPT Implementation 3 (Best n2 paths) Implementation 4 (Best n paths) SN Name n BKS BestGap AvgGapnRuns %BestBestGap AvgGap nRuns %BestBestGapAvgGapnRuns %Best 1 Chr12a 12 9552 0 45.769,436,2200.89031.212,824,7412.34028.983,807,8901.45 2 Chr12b 12 9742 0 52.788,402,6586.42 0 45.373,002,7108.26 0 44.384,572,6628.64 3 Chr12c 12 11156 0 34.929,823,1500.32 0 24.952,957,3610.76 0 23.823,907,7860.52 4 Chr15a 15 9896 0 50.104,704,6520.07 0 34.061,420,4160.38 0 32.432,021,0110.41 5 Chr15b 15 7990 0 60.114,468,5910.17 0 47.091,537,8170.25 0 46.642,266,8580.37 6 Chr15c 15 9504 0 61.614,920,5120.04 0 43.861,395,4820.09 0 41.792,015,7130.19 7 Chr18a 18 11098 0 66.832,730,1900.01 0 50.01779,6530.17 0 49.761,285,0440.06 8 Chr18b 18 1534 0 14.913,061,7780.57 0 8.53749,0022.84 0 14.342,699,9020.92 9 Chr20a 20 2192 0 47.812,102,1130.00 0 33.83531,1980.02 0 43.111,537,4700.00 10 Chr20b 20 2298 0 41.482,248,1470.00 0 27.59545,3280.00 0 33.181,304,3950.00 11 Chr20c 20 14142 0 83.481,638,1260.08 0 63.75537,3040.20 0 68.171,004,5250.23 12 Chr22a 22 6156 0 13.751,409,3410.00 0 10.03321,2330.02 0 12.731,094,9690.00 13 Chr22b 22 6194 0.581 13.981,524,7840.00 0 9.55306,9670.00 0 11.86935,4200.00 14 Chr25a 25 3796 0 57.88921,9460.00 0 44.08261,4260.00 0 52.94700,7770.00 15 Els19 19 17212548 0 25.631,703,9531.91 0 14.47298,48823.78 0 24.861,501,6741.91 16 Esc16a 16 68 0 3.546,915,00635.00 0 0.361,187,39394.63 0 3.546,849,54335.00 17 Esc16b 16 292 0 0.019,060,25498.44 0 0.001,063,243100.00 0 0.018,988,69198.44 18 Esc16c 16 160 0 1.155,753,64953.58 0 0.251,140,44184.19 0 1.155,702,56853.58 19 Esc16d 16 16 0 9.187,089,06643.13 0 0.721,165,76894.29 0 9.187,026,22443.13 20 Esc16e 16 28 0 10.128,572,39120.17 0 3.251,165,22360.37 0 10.128,502,35020.17 21 Esc16f 16 0 0 0.0035,611,664100.00 0 0.004,918,942100.00 0 0.0035,420,396100.00 22 Esc16g 16 26 0 7.867,348,48844.39 0 0.231,124,15297.03 0 7.867,287,15244.39 23 Esc16h 16 996 0 0.007,922,830100.00 0 0.001,051,718100.00 0 0.007,857,457100.00 24 Esc16i 16 14 0 1.098,107,02596.53 0 0.001,712,511100.00 0 1.098,041,28896.53 25 Esc16j 16 8 0 15.499,128,13054.33 0 1.28940,74894.90 0 15.479,037,78654.34 26 Esc32a 32 130 0 23.16591,4840.00 0 13.48105,1970.13 0 23.16589,3190.00 27 Esc32b 32 168 0 28.83573,4800.49 0 15.60119,8283.89 0 28.83571,8740.49 28 Esc32c 32 642 0 0.40806,97082.10 0 0.01116,68099.87 0 0.40804,16182.10 29 Esc32d 32 200 0 6.42789,9864.49 0 3.57115,08221.64 0 6.42787,7134.49 30 Esc32e 32 2 0 0.001,862,801100.00 0 0.00118,660100.00 0 0.001,858,335100.00 31 Esc32f 32 2 0 0.001,861,132100.00 0 0.00118,654100.00 0 0.001,858,965100.00 32 Esc32g 32 6 0 0.641,868,07198.08 0 0.00116,371100.00 0 0.641,865,82098.08 Contd ... PAGE 171 160Table A.1 (contd.): Experimental Results for Symmetric Instances 2OPT Implementation 3 (Best n2 paths) Implementation 4 (Best n paths) SN Name n BKS BestGap AvgGapnRuns %BestBestGap AvgGapnRuns %BestBestGapAvgGapnRuns %Best 33 Esc32h 32 438 0 3.93732,0211.55 0 2.0082,9429.08 0 3.93728,7411.55 34 Esc64a 64 116 0 1.48302,10248.8800.1421,57195.6901.48301,98648.88 35 Esc128 128 64 0 13.1038,9948.2705.611,85033.89013.0939,1018.26 36 Had12 12 1652 0 1.318,903,6174.9700.721,571,86318.0101.206,737,8576.69 37 Had14 14 2724 0 0.874,681,89013.6800.48979,81532.9200.823,890,03114.49 38 Had16 16 3720 0 0.903,046,48313.7500.48674,61026.1900.872,658,71213.81 39 Had18 18 5358 0 1.102,224,9492.6600.78407,7885.3201.072,016,2322.88 40 Had20 20 6922 0 1.191,562,4912.2700.84289,1896.7601.171,435,7512.29 41 Kra30a 30 88900 0 7.70508,2860.0205.9982,8460.2007.70507,3630.02 42 Kra30b 30 91420 0 5.76505,3390.0004.1383,3490.0405.76504,5160.00 43 Nug12 12 578 0 5.2910,480,6161.3503.482,387,4426.5004.856,507,0512.36 44 Nug14 14 1014 0 5.026,045,9030.4403.421,096,6240.9504.703,920,2300.86 45 Nug15 15 1150 0 4.464,738,2361.4602.791,005,8094.9404.023,224,1361.97 46 Nug16a 16 1610 0 4.863,835,4990.2803.43698,4541.6104.562,729,9220.53 47 Nug16b 16 1240 0 5.403,854,4263.1203.64756,4496.9804.932,615,3184.06 48 Nug17 17 1732 0 4.243,103,2730.0802.73553,2631.0804.012,283,7450.14 49 Nug18 18 1930 0 4.442,633,3710.1703.09467,1700.8504.151,876,4540.25 50 Nug20 20 2570 0 4.191,856,5940.2203.04342,8350.6703.971,386,3500.26 51 Nug21 21 2438 0 4.581,481,6780.1103.14283,3270.2904.351,188,5650.14 52 Nug22 22 3596 0 3.721,195,6120.5002.71244,2871.3703.591,015,2300.53 53 Nug24 24 3488 0 4.67967,8960.1003.31181,8180.7004.45787,5860.26 54 Nug25 25 3744 0 3.90842,3540.0502.64160,1320.4203.80726,5710.06 55 Nug27 27 5234 0 4.33626,4720.0403.27115,6440.4004.16527,5710.11 56 Nug28 28 5166 0 4.51579,2780.0103.30100,5790.1904.32481,5220.02 57 Nug30 30 6124 0 4.19453,5520.0103.0686,1790.0304.10398,4080.01 58 Rou12 12 235528 0 5.6410,263,9300.6904.252,348,1741.8304.654,692,4851.09 59 Rou15 15 354210 0 6.905,230,1450.3705.261,081,2790.9306.062,631,3060.57 60 Rou20 20 725522 0 4.862,114,5400.0103.34344,8480.0204.361,256,4720.01 61 Scr12 12 31410 0 7.039,411,2215.0504.672,337,37511.8204.824,371,64216.64 62 Scr15 15 51140 0 10.244,273,1221.7607.661,070,2116.8007.652,114,3905.51 63 Scr20 20 110030 0 9.861,709,0330.0606.22378,0050.2106.52896,5900.18 64 Sko42 42 15812 0.101 3.53294,0680.0002.7353,9250.010.1013.48270,2620.00 65 Sko49 49 23386 0.162 3.08176,4420.000.1542.4433,0020.000.1623.03165,9440.00 66 Sko56 56 34458 0.430 2.96108,9710.000.1632.3920,6230.000.2612.92103,0330.00 67 Sko64 64 48498 0.421 2.7068,9010.000.2932.2012,1620.000.4002.6866,5370.00 68 Sko72 72 66256 0.420 2.6545,8170.000.5462.217,8320.000.4202.6444,5700.00 69 Sko81 81 90998 0.510 2.2731,0870.000.4401.915,2730.000.5102.2630,5210.00 Contd ... PAGE 172 161Table A.1 (contd.): Experimental Results for Symmetric Instances 2OPT Implementation 3 (Best n2 paths) Implementation 4 (Best n paths) SN Name n BKS BestGap AvgGapnRuns %BestBestGap AvgGapnRuns %BestBestGapAvgGapnRuns %Best 70 Sko90 90 115534 0.535 2.2321,4730.000.4601.89 3,4250.000.5352.2221,1260.00 71 Sko100a 100 152002 0.617 2.0814,9550.000.5291.74 2,1800.000.6172.0714,7500.00 72 Sko100b 100 153890 0.516 2.0215,1100.000.3771.71 2,2220.000.5162.0214,9080.00 73 Sko100c 100 147862 0.580 2.3014,8430.000.5111.96 2,1820.000.5802.2914,6330.00 74 Sko100d 100 149576 0.646 2.0715,1700.000.5131.74 2,1790.000.6462.0715,0380.00 75 Sko100e 100 149150 0.581 2.3114,8200.000.4571.94 2,1870.000.5812.3014,6150.00 76 Sko100f 100 149036 0.780 2.0315,2990.000.5731.71 2,2150.000.6662.0215,0830.00 77 Ste36a 36 9526 0.252 12.02238,2390.0009.07 44,9390.010.25212.02239,8270.00 78 Ste36b 36 15852 0 21.46213,1960.01015.95 47,3380.18021.46214,1630.01 79 Ste36c 36 8239.11 0.132 9.46226,5670.0007.24 45,4280.000.1329.46227,5120.00 80 Tai12a 12 224416 0 8.9310,073,8662.1206.88 1,715,4034.7907.864,651,0373.09 81 Tai15a 15 388214 0 4.815,387,0380.1203.41 986,2140.4004.072,661,7770.15 82 Tai17a 17 491812 0 5.693,620,2560.0504.22 627,4780.2104.891,846,4510.17 83 Tai20a 20 703482 0 6.032,214,3710.0004.37 340,5570.0205.111,147,6830.01 84 Tai25a 25 1167256 0 5.571,109,3100.0004.08 173,7440.0004.76609,1590.00 85 Tai30a 30 1818146 0.456 5.06615,7030.000.5323.84 78,6380.0004.47371,3440.00 86 Tai35a 35 2422002 1.083 5.10387,1540.000.6873.72 46,2360.001.0754.47237,9490.00 87 Tai40a 40 3139370 1.358 5.05254,3040.000.9063.68 28,8210.001.3574.56169,1850.00 88 Tai50a 50 4941410 1.897 4.97253,4800.001.4373.71 26,5590.001.8634.46170,0790.00 89 Tai60a 60 7208572 2.177 4.72138,6430.001.6583.54 13,2000.001.9024.2091,4410.00 90 Tai64c 64 1855928 0 0.44416,5186.0900.42 27,9856.1900.44418,3916.09 91 Tai80a 80 13557864 2.037 3.7656,2680.001.5582.78 4,4840.001.7133.2236,1220.00 92 Tai100a 100 21125314 1.969 3.4227,5340.001.6082.49 1,7890.001.7073.0319,1840.00 93 Tai256c 256 44759294 0.175 0.433,5260.000.1750.41 1680.000.1750.433,3920.00 94 Tho30 30 149936 0 4.85432,9140.0103.72 86,3140.0304.60336,2870.01 95 Tho40 40 240516 0.341 4.67174,6380.000.0913.70 35,0560.000.0434.46142,4000.00 96 Tho150 150 8133484 0.877 2.423,5780.000.8142.12 4690.000.9002.392,5890.00 97 Wil50 50 48816 0.213 1.6618,3310.000.0861.37 30,1010.000.1231.65149,6640.00 98 Wil100 100 273038 0.363 1.1015,0050.000.3470.95 2,0310.000.3631.1014,9210.00 PAGE 173 162 Table A.2 Experimental Results for Asymmetric Instances 2OPT Implementation 3 (Best n2 paths) Implementation 4 (Best n paths) SN Name n BKS BestGapAvgGapnRuns %BestBestGapAvgGapnRuns %BestBestGapAvgGapnRuns %Best 1 bur26a 26 5426670 00.32114,6720.3200.2485,6341.1100.32114,6060.32 2 bur26b 26 3817852 00.41122,3080.2900.3087,2930.5800.41122,2470.29 3 bur26c 26 5426795 00.41112,6250.3100.2682,9121.0500.41112,5560.31 4 bur26d 26 3821225 00.47119,8090.2800.2987,7600.5400.47119,7430.28 5 bur26e 26 5386879 00.37110,9530.3900.2083,7872.6500.37110,8900.39 6 bur26f 26 3782044 00.46119,9570.6600.2487,3913.0700.46119,8910.66 7 bur26g 26 10117172 00.36109,9200.4700.2583,3472.4400.36109,8290.47 8 bur26h 26 7098658 00.46116,7440.9500.3593,5404.5000.46116,6770.95 9 lipa20a 20 3683 02.84386,8560.3702.52238,0811.3702.80285,7140.50 10 lipa20b 20 27076 015.22364,5844.27012.66277,40712.94014.80233,9704.96 11 lipa30a 30 13178 02.02104,9870.0301.8354,0750.2402.02104,9390.03 12 lipa30b 30 151426 016.79102,5981.88014.9656,8567.34015.9461,5504.03 13 lipa40a 40 31538 0.9071.5243,0180.0001.3618,0820.010.9071.5136,5660.00 14 lipa40b 40 476581 018.6441,1861. 26016.9019,0345.80018.2329,4782.03 15 lipa50a 50 62093 0.8921.3142,9780.000.8491.1715,3720.000.8921.3037,0430.00 16 lipa50b 50 1210244 018.6741,6330. 46017.5216,6442.33018.4030,9640.70 17 lipa60a 60 107218 0.8441.1124,2480.000.7870.997,1660.000.8441.1021,3390.00 18 lipa60b 60 2520135 020.0824,0460.09019.227,9380.42019.7116,7250.22 19 lipa70a 70 169755 0.7750.9615,0160.000.7250.864,0560.000.7730.9513,1870.00 20 lipa70b 70 4603200 020.8014,6440.05019.944,3580.53020.379,7670.17 21 lipa80a 80 253195 0.6780.8510,0260.000.6280.752,3350.000.6780.838,5100.00 22 lipa80b 80 7763962 021.699,7010.03020.922,4510.08021.236,2140.06 23 lipa90a 90 360630 0.6360.786,8360.000.5810.691,4120.000.6360.776,2310.00 24 tai12b 12 39464925 011.251,619,2583.0308.941,682,70412.55010.35953,5113.56 25 tai20b 20 122455319 017.45266,2890.60014.22207,1186.30015.13183,1130.88 26 tai25b 25 344355646 016.24121,8620.03012.1083,6190.59015.5389,5290.08 27 tai30b 30 637117113 013.7265,7000.0009.0638,0470.13012.6528,9450.00 28 tai35b 35 283315445 08.7240,0950.0006.4723,2130.0308.1231,9320.01 29 tai40b 40 637250948 010.5924,6960.0008.2913,5270.24010.1819,7420.01 30 tai50b 50 458821517 0.0737.2423,7340.000.0585.7311,3150.000.0737.1020,2530.00 31 tai60b 60 608215054 0.0087.9912,7130.000.0386.164,8630.000.0897.7510,8530.00 32 tai80b 80 818415043 1.5946.114,9880.000.8445.271,5970.001.5946.014,5630.00 33 tai100b 100 1185996137 0.8755.312,2980.000.6534.435630.000.7625.202,1160.00 34 tai150b 150 498896643 1.4733.496240.001.4673.10920.001.4733.445840.00 PAGE 174 163 LIST OF REFERENCES Aarts, E. H. L., and J. Korst. 1989. Simulated Annealing and Boltzman Machines: A Stochastic Approach to Combinatoria l Optimization and Neural Computing. Wiley, Chichester. Aarts, E., and J.K. Lenstra. 1997. Local Search in Combinatorial Optimization. John Wiley, New York. Ahuja, R.K., O. Ergun, J.B. Orlin, and A.P. Punnen. 2002a. A survey of very largescale neighborhood search techniques. Discrete Applied Mathematics 123, 75102. Ahuja, R.K., J. Liu, J.B. Orlin, D. Sharma and L.A. Shughart. 2002b. Solving reallife locomotive scheduling problems. Submitted to Transportation Science. Ahuja, R.K., J. Liu, J. Goodstein, A. M ukherjee, J.B. Orlin, and D. Sharma. 2003b. Solving multicriteria combined th roughfleet assignment models. In the Operations Research in Space and Air, Edited by T. A. Ciriani, G. Fasano, S. Gliozzi, and R. Tadei, Kluwer Academic Publishers, Boston pp. 233256. Ahuja, R.K., J. Liu, J. Goodstein, A. Mukherjee, and J.B. Orlin. 2003c. A neighborhood search algorithm for the combined through and fleet assignment model with time windows. Submitted to Networks. Ahuja, R.K., T.L. Magnanti, J.B. Orlin. 1993. Network Flows: Theory, Algorithms, and Applications. Prentice Hall, Upper Saddle River, New Jersey. Ahuja, R.K., J.B. Orlin, S. Pallottino, M.P Scaparra, and M.G. Scutella. 2002c. A multiexchange heuristic for the single source capacitated facility location. Submitted to Management Science. Ahuja, R. K., J. B. Orlin, D. Sharma 2001a. Multiexchange neighborhood search algorithms for the capacitated minimum spanning tree problem. Mathematical Programming 91, 7197. PAGE 175 164 Ahuja, R. K., J. B. Orlin, D. Sharma. 2001b. A very largesc ale neighborhood search algorithm for the combined throughfl eet assignment model. Submitted to INFORMS Journal on Computing. Ahuja, R.K., J.B. Orlin, and D. Shar ma. 2003a. A composite very largescale neighborhood structure for the capacita ted minimum spanning tree problem. Operations Research Letters 31, 185194. Ahuja, R. K., J. B. Orlin, and A. Tiwa ri. 2000. A greedy genetic algorithm for the quadratic assignment problem. Computers and Operations Research 27, 917934. Armour, G. C., and E. S. Buffa. 1963. Heur istic algorithm and si mulation approach to relative location of facilities. Management Science 9, 294309. Assad, A.A. 1980. Models for rail transportation. Transportation Research A14, 205220. Assad, A.A. 1982. A class of trainscheduling problems. Transportation Science 16. 281310. Barnhart, C., H. Jin, and P.H. Vance. 2000. Railroad blocking: A network design application. Operations Research 48, 603614. Bazara, M. S., and M. D. Sherali. 1980. Bende rs' partitioning scheme applied to a new formulation of the quadratic assignment problem. Naval Research Logistics Quarterly 27, 2941. Bean, J. C. 1994. Genetic algorithms and ra ndom keys for sequencing and optimization. ORSA Journal on Computing 6, 154160. Bodin, L.D., B.L. Golden, A. Schuster, and W. Romeg. 1980. A m odel for the blocking of trains. Transportation Research 14B, 115120. Braford, J.C. 1961. Determination of optimal assignment of a weapon system to several targets. AEREITM9, Vought Ae ronautics, Dallas, Texas. Brannlund U., P.O. Lindberg, A. Nou, and J.E. Nilson. 1998. railway timetabling using Lagrangian relaxation Transportation Science 32, 358369. Buffa, E. S., G. C. Armour, and T. E. Vollmann. 1964. Allocating facilities with CRAFT. Harvard Business Review 42, 136158. PAGE 176 165 Burkard, R. E. 1991. Location with spatia l interactions: The quadratic assignment problem. Discrete Location Theory, Edited by P. B. Mirchand ani and R. L. Francis, John Wiley. Burkard, R. E., and T. Bonniger. 1983. A heuristic for quadra tic boolean programs with applications to quadratic assignment problems. European Journal of Operations Research 13, 374386. Burkard, R. E., E. Cela, P. M. Pardalos, and L. S. Pitsoulis. 1998. The quadratic assignment problem. In Handbook of Combinatorial Optimization 3, Edited by D. Z. Zhu and P. M. Pardalos, Kluwer, Boston, pp. 241337. Burkard, R. E., S. E. Karisch, and F. Rendl. 1997. QAPLIB A quadratic assignment program library. Journal of Global Optimization 10, 391403. Campbell K.C. 1996. Booking and revenue mana gement for rail intermodal services. PhD Dissertation. Department of Systems Engi neering, University of Pennsylvania, Philadelphia, PA. Carter M. W., and Laporte G. 1998. Recent Developmen ts in Practical Course Timetabling. PATAT97, LNCS 1408, 319 Castanon D.A. 1987. Advanced weapontarget assignment algorith m quarterly report. TR337, ALPHA TECH In c., Burlington, MA. Cela, E. 1998. The Quadratic Assignment Problem Theory and Algorithms, Kluwer Academic Publishers, Boston. Chang S.C., R.M. James, and J.J. Shaw. 1987. Assignment algorithm for kinetic energy weapons in boost defense. Proceedings of the IEEE 26th Conference on Decision and Control, Los Angeles, CA. Christofides, N., and E. Benavent. 1989. An exact algorithm for the quadratic assignment problem. Operations Research 37, 760768. Cordeau J.F., P. Toth, and D. Vigo. 1998. A survey of optimization models for train routing and scheduling. Transportation Science 32, 380404. Crainic T.G., J.A. Ferland, and J.M. Rousseau 1984. A tactical planning model for rail freight transportation. Transportation Science 18, 165184. PAGE 177 166 Crainic, T.G., and J.M. Rousseau. 1986. Multicommodity, miultimode freight transportation: A general modeling and algorithmic framework for the service network design problem. Transportation Research 208. 225242. Croes, G.A. 1958. A method for solving traveling salesman problems. Operations Research 6, 791812. Davis. L. 1991. Handbook of Genetic Algorithms. Van Nostrand, New York. Day, R.H. 1966. Allocating weapons to target complexes by means of nonlinear programming. Operations Research 14, 9921013. Deineko V. G., and G.J. Woeginger. 2000. A st udy of exponential neighborhoods for the travelling salesman problem and for the quadratic assignment problem. Mathematical Programming 87(3), 519542. DenBroader G. G., Jr., R. E. Ellison, and L. Emerling. 1958. On optimum target assignments. Operations Research 7, 322326. Dimopoulou M., Miliotis P. 2001. Theory and methodology implementation of a university course and examination timetabling system. European Journal of Operational Research 130, 202513. Drezner, Z. 2003. A new genetic algorithm fo r the quadratic assignment problem. INFORMS Journal on Computing 15(3), 320330. Eckler, A.R., and S.A. Burr. 1972. Mathematical models of target coverage and missile allocation. Military Operations Re search Society, Alexandria, VA. Farvolden J.M., and W.B. Powell. 1994. Subgr adient methods for the service network design problem. Transportation Science 28, 256272. Fisher M.L. 1981. The Lagrangian relaxatio n method for solving integer programming problems. Management Science 27, 118. Fleurent, C., and J. A. Ferland. 1994. Ge netic hybrids for the quadratic assignment problem. DIMACS Series in Discrete Mathematics and Theoretical Computer Science, American Mathematical Society 16, 173187. Glover, F. 1986. Future paths for integer prog ramming and links to ar tificial intelligence. Computers and Operations Research 13, 533549. PAGE 178 167 Glover, F., and M. Laguna. 1997. Tabu Search. Kluwer Academic Publishers, Norwell, MA. Gorman, M.F. 1998. The freight rail road operating plan problem. Annals of Operations Research 78, 5169. Grant, K.E. 1993. Optimal resource al location using genetic algorithms. 1993. Naval Review, Naval Research Laboratory, Washington, DC, pp. 174175. Green D.J., J.T. Moore, and J.J. Borsi. 1997 An integer solution heur istic for the arsenal exchange model (AEM). Military Operations Research Society, 3(2), 515. Haghani A.E. 1989. Formulation and solution of a combined train routing and makeup, and empty car distribution model. Transportation Research 23B, 433452. Hertz A. 1991. Tabu Search for large scale timetabling problems. European Journal of Operational Research 54, 3947. Holland, J.H. 1975. Adaptation in Natural and Artificial Systems. University of Michigan Press, Ann Arbor, MI. Ibaraki, T. 1987. Enumerative Approaches to Combinatorial Optimization. Baltzer, Basel. Katter, J.D. 1986. A solution of the multiweapon, multitarget assignment problem. Working paper 26957, MITRE. Keaton M.H. 1989. Designing optimal railroad ope rating plans: Lagrangian relaxation and heuristic approaches. Transportation Research 23B, 415431. Keaton, M.H. 1991. Servicecost tradeoffs for carload freight traffic in the U.S. rail industry. Transportation Research 25A, 363374. Keaton, M.H. 1992a. Designing optimal railro ad operating plans: A dual adjustment method for implementing Lagrangian relaxation. Transportation Science 26, 262279. Keaton, M.H. 1992b. The impact of train timetables on average car time in rail classification yards. Journal of the Transportation Research Forum 32, 345354. Kiaer L., and Yellen J. 1992. Weighted gr aphs and university course timetabling. Computers and Operations Research 19(1), 5967. PAGE 179 168 Kirkpatrick, S. C.D. Gellatt Jr., and M.P. Vecchi. 1983. Optimization by simulated annealing. Science 220, 671680. Kraft E.R. 1998. A reservationbased railway network operations management system. PhD dissertation. Department of Systems E ngineering, University of Pennsylvania, Philadelphia, PA. Kwon O.K., C.D. Martland, and J.M. Sussman 1998. Routing and scheduling temporal and heterogeneous freight car traffic on rail networks. Transportation Research 34E, 101115. Lawler, E. L. 1963. The quadr atic assignment problem. Management Science 9, 586599. Li, Y., P. M. Pardalos, and M. G. C. Resende. 1994. A greedy randomized adaptive search procedure for the quadratic assignment problem. In Quadratic Assignment and Related Problems, Edited by P. M. Pardalos and H. Wolkowicz, DIMACS Series in Discrete Mathematics and Theoretical Computer Science, American Mathematical Society, Providence, RI, pp. 237261. Lin, S. 1965. Computer solutions of the traveling salesman problem. Bell System Technical Journal 44, 22452269. Lin, S., and B.W. Kernighan. 1973. An eff ective heuristic algorit hm for the travelingsalesman problem. Operations Research 21, 498516. Liu J. 2003. Solving reallife tran sportation scheduling problems. PhD dissertation. Department of Industrial and Systems E ngineering, University of Florida, Gainesville, FL. Lloyd, S.P., and H.S. Witsenhausen. 1986. Weapons allocation is NPcomplete. In Proceedings of the 1986 Summe r Conference on Simulation, Reno, NV. Maltin, S.M. 1970. A review of the litera ture on the missile allocation problem. Operations Research 18, 334373. Malucelli, F. 1993. Quadratic assignment problems: solution met hods and applications. PhD Doctoral Dissertation, De partimento di Informatica, Universita di Pisa, Italy. Maniezzo, V., A. Colorni, and M. Dorigo. 1994. The ant system applied to the quadratic assignment problem. Tech. Rep. IRIDIA/ 9428, Universit Li bre de Bruxelles, Belgium. PAGE 180 169 Manne, A.S. 1958. A targetassignment problem. Operations Research 6, 346351. Metler, W.A., and F.L. Preston. 1990. A suite of weapon assignment algorithms for a SDI midcourse battle manager. NRL Memorandum Report 671, Naval Research Laboratory, Washington, DC. Miladenovic, N., and P. Hansen. 1997. Variable neighborhood search. Computers and Operations Research 34, 10971100. Morlok, E.K., and R.B. Petersen. 1970. Railroad Freight Train Scheduling, a Mathematical Programming Formulation. Transportation Center and Technological Institute, No rthwestern University, Evanston, IL. MullerMerbach, H. 1970. Optimale Reihenfolgen. Springer Verlag, Berlin, 158171. Murphey, R.A. 1999. Targetbased w eapon target assignment problems. Nonlinear Assignment Problems: Algor ithms and Applications, Edited by P. M. Pardalos and L. S. Pitsoulis, Kluwer Academic Publishers, Boston, pp. 3953. Murty, K.G. 1976. Linear and Combinatorial Optimization. John Wiley, New York. Newman A.M., L. Nozick, and C.A. Yano. 2002. Optimization in the rail industry. Handbook of Applied Optimization, Edited by P. M. Pardalos and M. Resende, Oxford University Press, New York. Newton, H.N., C. Barnhart, and P.M. Vance. 1998. Constructing rail road blocking plans to minimize handling costs. Transportation Science 32, 330345. Nozick L.K., and E.K. Morlok. 1997. A mode l for mediumterm operations planning in an intermodal railtruck service. Transportation Research 31A, 91107. Nugent, C. E., T. E. Vollman, and J. Ru ml. 1968. An experimental comparison of techniques for the assignment of facilities to locations. Operations Research 16, 150173. Orlin, D. 1987. Optimal weapons allo cation against layered defenses. Naval Research Logistics 34, 605616. Osman, I.H., and G. Laporte. 1996. Metaheuristics: A bibliography. Annals of Operations Research 63, 513623. PAGE 181 170 Pardalos, P. M., and J. Crouse. 1989. A para llel algorithm for the quadratic assignment problem. In Proceedings of the Supe rcomputing 1989 Conference, ACM Press, pp. 351360. Pardalos, P.M., F. Rendl, and H. Wolkowicz. 1994. The quadratic assignment problem. In Quadratic Assignment and Related Problems, Edited by P. M. Pardalos and H. Wolkowicz, DIMACS Series, American Math ematical Society, Providence, RI, pp. 142. Reeves, C.R. 1993. Modern Heuristic Techniques fo r Combinatorial Optimization. Blackwell Scientific Publications, Oxford. Resende, M. G. C., P. M. Pardalos, and Y. Li 1994. Fortran subroutines for approximate solution of dense quadratic assignment problems using GRASP. ACM Transactions on Mathematical Software 22, 104118. Rochat, Y., and E.D. Taillard. 1995. Probabilis tic diversification and intensification in local search for vehicle routing. Journal of Heuristics 1, 147167. Schaerf A. 1999. A survey of Automated Timetabling. Artificial Intelligence Review 13, 87127 SkorinKapov, J. 1990. Tabu search applie d to the quadratic assignment problem. ORSA Journal on Computing 2, 3345. Stallaert J. 1997. Automated timetabling improves course scheduling at UCLA. Interfaces 27, 6781. Taillard, E. 1991. Robust tabu search for the quadratic assignment problem. Parallel Computing 17, 443455. Talluri, K.T. 1996. Swapping applicatio ns in a daily fleet assignment. Transportation Science 31, 237248. Tate, D. E. and A. E. Smith. 1985. A genetic approach to the quadratic assignment problem. Computers and Operations Research 22, 7383. Thomet, M.A. 1971. A user oriented freight railroad operating policy. IEEE Trans. Systems, Man Cybernet 1, 349356. Thompson, P.M., and J.B. Orlin. 1989. The theo ry of cyclic transfer s. MIT Operations Research Center Report 200289, Cambridge, MA. PAGE 182 171 Thompson, P. M. and H.N. Psaraftis. 1993. Cy clic transfer algorith ms for multivehicle routing and scheduling problems. Operations Research 41, 935946. Thompson J. M. and Dowsland K. A. 1996. Va riants of simulated annealing for the examinationtimetabling problem. Annals of Operations Research 63, 105128. Valdes R.A., Crespo E., Tamarit J.M. 2002. Design and implementation of a course scheduling system using tabu search. European Journal of Operational Research 137, 512523. Van Dyke, C.D. 1986. The automated blocking model: a practical approach to freight railroad blocking plan development. Proceedings of the 27th Annual Meeting of the Transportation Research Forum 27, 116121. Wacholder, E. 1989. A neural networkbased optimization algorithm for the static weapontarget assignment problem. ORSA Journal on Computing 4, 232246. Werra D. D. 1997. The combinatorics of timetabling, European Journal of Operational Research 96, 504513. West, D. H. 1983. Algorithm 608: approximat e solution of the quadratic assignment problem. ACM Transactions on Mathematical Software 9, 461466. Wilhelm, M. R., and T. L. Ward. 1987. So lving quadratic assignment problems by simulated annealing. IEEE Transactions 19, 107119. PAGE 183 172 BIOGRAPHICAL SKETCH Krishna Jha is a native of India and got his early educat ion at Darbhanga, India. He obtained his bachelors degree in mechanical en gineering from B. I. T. Sindri, India, and his masters degree in industrial and management engineering from I. I. T. Kanpur, India. After graduation, he worked in the largest private integrated steel plant in India, Tata Steel, for four years. Since 2000, he has been pursuing his PhD degree in the department of Industrial and Systems Engineerin g at University of Florida. 