Optimizing Earthwork Block Removal in Road Construction by Valentin Raphael Koch B.Sc. Hons., The University of British Columbia, 2008 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in The College of Graduate Studies (Mathematics) THE UNIVERSITY OF BRITISH COLUMBIA (Okanagan) April 2010 c Valentin Raphael Koch 2010 Abstract In road construction, earthwork operations account for about 25% of the construction costs. Existing linear programming models for earthwork logistics optimization are designed to minimize the hauling costs and to balance the earth across the construction site. However, these models do not consider the removal of physical blocks that may influence the earthwork process. In this thesis, we extend the linear programming model of Mayer and Stark [MS81] with the addition of a block removal schedule. The resulting model is a mixed-integer linear program. We analyze the model size and the schedule search space in order to make conclusion about the use of the model. Based on structural observations, we introduce a set of algorithms that significantly reduce the solving time of the model. Finally, we conduct numerical experiments to compare our solutions with the solutions of a traditional earthwork process that makes use of linear programming. From our numerical results, we conclude that an optimal removal schedule produces solutions that are 4.1% cheaper on average than a traditional method, with savings that can go as high as 19%. We conclude our discussion with possible extensions to the model, that can help an engineer to design roads that are more economical and ecological with respect to the earthwork operations. ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi Abstract List of Tables Glossary of notation . . . . . . . . . . . . . . . . . . . . . . . . . . vii Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . viii Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 Introduction . . . . . . . . . . . . . . . 1.1 Road design . . . . . . . . . . . . . 1.1.1 Road design and construction 1.1.2 Road design optimization . . 1.1.3 Earthwork optimization . . . 1.2 Basic linear programming model . . 1.2.1 Notation . . . . . . . . . . . 1.2.2 Model parameters . . . . . . 1.2.3 Decision variables . . . . . . 1.2.4 Objective and constraints . . 2 Mixed integer programming model 2.1 Notation . . . . . . . . . . . . . . 2.2 Model parameters . . . . . . . . . 2.3 Decision variables . . . . . . . . . 2.4 Objective and constraints . . . . . 2.4.1 Objective function . . . . . 2.4.2 Volume constraints . . . . 2.4.3 Block removal indicator . . . . . . . . . . ix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1 1 2 4 7 7 7 9 10 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 14 14 15 16 17 17 17 . . . . . . . . . . . . . . . . iii Table of Contents 2.4.4 2.4.5 2.4.6 2.4.7 Move feasibility constraint Stockpile constraints . . . Block removal enforcement Observations . . . . . . . . . . . . . . . . 18 18 20 20 3 Model analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Space complexity . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Schedule search . . . . . . . . . . . . . . . . . . . . . . . . . 22 22 26 4 Algorithms . . . . . . . . . . . . . . . . . 4.1 Mass schedule as starting point . . . . 4.2 Maximal allowed moving distance dmax 4.3 Warmstarting with increasing dmax . 4.4 Warmstarting with mipgap check . . 5 Numerical results . . . . . . . . . . . . 5.1 Experimental setup . . . . . . . . . 5.1.1 Problems . . . . . . . . . . . 5.1.2 Benchmark . . . . . . . . . . 5.2 Results . . . . . . . . . . . . . . . . 5.2.1 Analysis of the dmax -heuristic 5.2.2 Performance . . . . . . . . . 5.2.3 Economical analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 36 37 39 40 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 42 42 42 46 46 48 54 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 57 57 59 59 59 60 61 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 . . . . . . . . . . . . . . . . . . . . . . . . . . 6 Conclusion . . . . . . . . . . . . . . . 6.1 Accomplishments . . . . . . . . . 6.2 Research outcomes . . . . . . . . . 6.3 Future research . . . . . . . . . . 6.3.1 Warmstart schemes and εLP 6.3.2 Vertical alignment . . . . . 6.3.3 Schedule enumeration . . . 6.3.4 Lagrangian relaxation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Appendices A Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 B Figures iv List of Tables 1.1 Major components of construction costs . . . . . . . . . . . . 3 2.1 Stockpile example data . . . . . . . . . . . . . . . . . . . . . . 18 3.1 Schedule example data . . . . . . . . . . . . . . . . . . . . . . 27 5.1 5.2 Test configurations for dmax . . . . . . . . . . . . . . . . . . . Economical test and performance test configurations . . . . . 45 45 A.1 Problem configurations . . . . . . . . . . . . . . . . . . . . . . A.2 Statistics for objective function error with dmax . . . . . . . . A.3 Statistics for solving time with dmax . . . . . . . . . . . . . . A.4 Histogram statistics for 1% tolerance . . . . . . . . . . . . . . A.5 Histogram statistics for 5% tolerance . . . . . . . . . . . . . . A.6 CPU time in seconds for ε = 0.01, problems 1-22 . . . . . . . A.7 CPU time in seconds for ε = 0.01, problems 23-44 . . . . . . A.8 CPU time in seconds for ε = 0.05, problems 1-22 . . . . . . . A.9 CPU time in seconds for ε = 0.05, problems 23-44 . . . . . . A.10 Objective value in $ for direct solve, ε = 0.01 . . . . . . . . . A.11 Objective value in $ for direct solve with x◦ , ε = 0.01 . . . . . A.12 Objective value in $ for warmstart (d0max = 25, α = 10), ε = 0.01 A.13 Objective value in $ for warmstart with ǫ (d0max = 25, α = 10), ε = 0.01 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.14 Objective value in $ for warmstart with ǫ (d0max = 5, α = 5), ε = 0.01 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.15 Objective value in $ for direct solve, ε = 0.05 . . . . . . . . . A.16 Objective value in $ for direct solve with x◦ , ε = 0.05 . . . . . A.17 Objective value in $ for warmstart (d0max = 25, α = 10), ε = 0.05 A.18 Objective value in $ for warmstart with ǫ (d0max = 25, α = 10), ε = 0.05 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.19 Objective value in $ for warmstart with ǫ (d0max = 5, α = 5), ε = 0.05 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 67 67 68 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 v List of Figures 1.1 1.2 Example of vertical alignment with mass diagram . . . . . . . Example of vertical alignment with sections . . . . . . . . . . 5 8 2.1 Earthwork scheme with removal schedule (black squares indicate blocks, gray squares indicate removed blocks). . . . . . 19 3.1 3.2 3.3 3.4 Potential schedules . . . . . . . Earthwork scheme with removal Schedule tree . . . . . . . . . . (6 × 4)-lattice in Z2 . . . . . . . . . . . . schedule. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 28 29 30 5.1 5.2 5.3 5.4 5.5 5.6 Objective and solving time with increasing dmax . . . . Algorithms performance profiles with log2 scale . . . . Algorithms performance profiles with small τ . . . . . Algorithms performance profiles with large τ . . . . . Savings with warmstart with ǫ (d0max = 25, α = 10) . . Stratification for warmstart with ǫ (d0max = 25, α = 10) . . . . . . . . . . . . . . . . . . . . . . . . 47 50 52 53 55 56 B.1 B.2 B.3 B.4 B.5 B.6 B.7 B.8 Savings with direct solve . . . . . . . . . . . . . . . . Savings with direct solve and x◦ . . . . . . . . . . . Savings with warmstart (d0max = 25, α = 10) . . . . . Savings with warmstart and ǫ (d0max = 5, α = 5) . . . Stratification for direct solve . . . . . . . . . . . . . Stratification for direct solve with x◦ . . . . . . . . . Stratification for warmstart (d0max = 25, α = 10) . . . Stratification for warmstart with ǫ (d0max = 5, α = 5) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 84 84 84 85 85 86 86 . . . . . . . . vi Glossary of notation a∈A R R+ φ:A→B [a, b] ∞ f′ |·| |A| {x : P (x)} Z+ A∪B dist(i, j) R++ min f A×B ∧ ∨ A\B ∅ Θ (·) O (·) Z++ A⊆B n k Rn XT Rm×n a←b a is an element of set A. 4 Real numbers. 4 Nonnegative real numbers. 5 The mapping A into B by φ. 5 Interval {x ∈ R : a ≤ x ≤ b}. 5 Infinity. 5 Derivative of function f . 6 Integral. 5 Absolute value. 7 Cardinality of set A. 7 The set of all x such that P (x) is satisfied. 8 Nonnegative integers. 8 The union of sets A and B. 8 Distance between section i and j. 8 Positive real numbers. 8 The minimum of function f . 10 Summation. 10 Cartesian product {(x, y) : x ∈ A, y ∈ B}. 12 Logical connective representing and. 14 Logical connective representing or . 14 Set difference {x : x ∈ A and x ∈ B}. 16 The empty set. 16 Asymptotic tight bound. 22 Asymptotic upper bound. 22 Positive integers. 22 A is a subset of B. 22 Binomial coefficient (n choose k). 30 Real n-vectors (n × 1 matrices). 35 Transpose of matrix X. 35 Real m × n matrices. 38 Assign value b to a. 37 vii Acknowledgements Many people were involved in the creation of this thesis. First, I would like to express my sincere gratitude to my supervisors, Dr. Yves Lucet and Dr. Warren Hare, for their help and support throughout this research. I greatly enjoyed our lively discussions and conversations. It was a great privilege to have supervisors that helped me to overcome so many challenges, whether they were in mathematics, computer science, or simply personal. I would also like to thank the members of my graduate committee, Dr. Heinz Bauschke and Dr. Shawn Wang, for their earnest support and guidance in my graduate studies, and for providing me insight into the beautiful subject of convex analysis. It is a pleasure to thank my colleagues, especially Jim Nastos and Mason Macklem, for their patience to listen and the many ideas and input they contributed. A very special thanks goes out to our industrial partner, for the professional assistance and valuable help on engineering questions. I owe my deepest gratitude to my parents. Without their prayers and support, my academic career would have been impossible. Finally, I would like to thank my wife Delphine, and my children Jordan and Tabitha, for their love and patience during these intense years. This work was partially supported by the Ministry of Advanced Education of the Province of British Columbia through a Pacific Century Graduate Scholarship, and by the Mathematics of Information Technology and Complex Systems network through a MITACS Accelerate internship. viii Dedication To Delphine, Jordan, and Tabitha. ix Chapter 1 Introduction La tendance la plus profonde de toute activit´e humaine est la marche vers l’´equilibre.1 Jean Piaget (1896–1980) 1.1 1.1.1 Road design Road design and construction Transportation of people and goods over land is, and has always been, an important part of civilization. Construction of the first roadways dates back as early as 5000 BCE [Lay99]. The invention of the wheel created the need for better roads and more sophisticated forms of road construction. Complex road design, including the alignment, earthwork, and pavement was already performed during the Roman Empire. The first scientific approach to road design was by Tr´esaguet in France at about 1764 [Lay99]; Tr´esaguet also introduced the first system of continuous road maintenance. The first comprehensive guide on road design and construction was published by McAdam in 1816 [McA24], introducing the concept of the elevation of roads through embankment. This resulted in much more stable roads, but also increased the cost of construction. Roadways are an essential part of the economy in most countries. For example, in Canada in 2008, trucking accounted for 54% of the value of trade with the United States and commercial transportation services accounted for 4.1% of Canada’s value-added GDP [Stu09]. At the same time, all levels of Canadian government combined spent $20.9 billion, or 1.6% of Canada’s GDP, to build and maintain its 1 million kilometer road network. In view of these numbers, it is a natural consequence that governments of countries that depend heavily on road networks have a vast interest in keeping a 1 The most profound tendency of all Human activity is progression toward equilibrium. 1 1.1. Road design high quality road infrastructure, while at the same time trying to lower the expenditures on construction and maintenance. From an economic perspective, it is therefore important to try to minimize the cost and maximize the quality and safety of the road infrastructure. The design of a road largely depends on safety requirements, which in turn depend on the traffic volume and vehicle speed. In 2007, in Canada, about 140,000 road collisions occurred, with about 16,000 seriously injured persons and about 3,000 fatalities. More than 100 million tonnes of dangerous goods were transported on Canadian roads during the same year [Stu09]. While the economical aspects in road construction are important, safety requirements are essential. Depending on the location, the construction of a road can have a great impact on the environment. Maki [MKV02] shows that poorly planned road construction, especially in sensitive areas like the Amazon region, often harms the local economy instead of helping it. According to Maki, it is important to optimize the land use through environmental and social impact analysis. Appropriate classification and zoning of land should therefore be taken into account, and environmental impacts, like deforestation or earthwork, should be minimized when a road is built. It follows from the above discussion that the road infrastructure needs to be optimized with respect to economic and social costs, safety requirements, and environmental impact. 1.1.2 Road design optimization Multiple factors determine the total transportation cost of a road. Jha et. al. [JSJ06] list five main types of costs that affect the total transportation cost: − planning, design, and administrative costs; − construction costs; − operation and maintenance costs; − user costs; and − social and environmental costs. Planning and administrative costs are usually not considered in road design optimization. In [CGF89], Chew et. al. divide construction costs in six major categories, shown in Table 1.1 with an approximate associated 2 1.1. Road design Table 1.1: Major components of construction costs Pavement Earthwork (half of this is haul cost) Bridges Drainage Miscellaneous items Land 30% 25% 20% 10% 10% 5% percentage of the total construction cost. Operation and maintenance costs include roadway surface, roadside features, bridges, tunnels, snow and ice control, and traffic control devices and amount at about 5% of the construction costs, over a period of 30 years. The user costs, which include vehicle operating costs, travel time, and traffic accident costs, range from 300% to 1000% of the construction costs over a period of 30 years [JSJ06]. Finally, social and environmental costs can be the most critical factor in highway design. It is, however, difficult to quantify these costs and they are usually considered before the optimization of the design. The design of the road can roughly be divided into geometric design and pavement design. The geometric design determines the alignment of the road and affects the construction costs, user costs, and social and environmental costs. The pavement design affects operation and maintenance, as well as user costs. Optimization of the geometric design is traditionally done in two stages: the horizontal alignment design and the vertical design. The horizontal alignment is the road trajectory from a satellite’s eye view. The primary objects that influence a horizontal alignment are forbidden areas like lakes, ocean shores, steep mountains, nature reserves, and private property. Cost parameters in horizontal alignments might include land costs and road distance. Early attempts at horizontal alignment optimization were based on dynamic programing [Tri87]. More recent approaches use genetic algorithms [JK06] or local neighborhood heuristics [LTL09]. Once the horizontal alignment is done, one can draw a line from the starting point to the end point of the road and look at the vertical ground profile along this line. The vertical alignment optimization tries to fit a road to the vertical ground profile as close as possible, while respecting grade constraints and other road specifications. The closer the road is to the ground profile, the less earthwork needs to be done in order to cut or fill sections of the road. Usually, vertical alignment optimization includes earthwork allocation and earth-moving optimization. In civil engineering, 3 1.1. Road design the vertical alignment optimization is commonly done through a manual iteration of computations that estimate the earthwork volumes for a given alignment. The alignment is then changed and the volumes are recomputed until a satisfying solution is found. Traditional mathematical models that automate the optimization of the vertical alignment use the method of linear programming (LP) [Eas88b, Mor96, Mor09, KL10]. Recently, the road alignment optimization research community has focused on the simultaneous optimization of horizontal and vertical alignment. In [CGF89], Chew et. al. modeled the road as a three-dimensional cubic spline. Jong and Schonfeld [JS03] use genetic algorithms to first solve for the three-dimensional alignment as a three-dimensional piecewise linear spline, and then fit circular curves to the horizontal and vertical projections of the spline. Parametric curves were studied by Makanae [Mak99] to model the road in three dimensions. Cheng and Lee [CL06] use a heuristic bi-level approach to iteratively optimize horizontal and vertical alignments. All of the above optimization methods depend on the cost evaluation of the earthwork. In this thesis, we focus on the aspect of earthwork optimization for a given alignment. The optimal cost of earthwork is essential for the assessment of a given road alignment. 1.1.3 Earthwork optimization Earthwork depends significantly on the terrain and the geometric design of the road. A road in a mountainous area requires more cutting and filling of earth than a road on the plain. An optimal vertical grade line minimizes the earthwork and satisfies all the design constraints for the road. To determine the optimality with respect to the earthwork, one needs to be able to compute the earthwork volumes for a given alignment. In order to compute the earthwork volumes for a given grade line, crosssections of the road profile are taken at regular intervals. The point along the horizontal axis where the cross-section is taken is called a station. Stations are usually arranged at an equal distance from one another along the horizontal axis of the grade line. The distance Lij is the actual length between station i and j along the horizontal center line of the road. The area Ai of a cross-section at station i is determined and the volume V ∈ R between two stations i and j is computed as V = Lij (Ai +Aj )/2. A traditional approach uses the mass diagram to depict the net accumulation of cut or fill between any two stations along the horizontal axis of the grade line. In [GH08], Garber and Hoel describe the mass diagram as follows. The ordinate of the mass diagram is the net accumulation 4 1.1. Road design road profile, g(¯ x) ground profile, f (¯ x) balance point balance point x ¯ x ¯ cumulative mass, h(¯ x) = 0 f (u) − g(u) du Figure 1.1: Example of vertical alignment with mass diagram in cubic meter from an arbitrary starting point. Thus, the difference in ordinates between any two stations represents the net accumulation of cut or fill between these stations. If the first station of the roadway is considered to be the starting point, then the net accumulation at this station is zero. For this thesis, we will use the first station as the starting point. We can therefore define the mass diagram as follows. Definition 1.1. Consider a vertical alignment. Let x ¯ be a station. If f (¯ x) represents the height of the ground profile at x ¯, and g(¯ x) represents the height of the road profile at x ¯, where f : R+ → R and g : R+ → R are piecewise continuous functions on [0, +∞[, then the mass diagram is the graph of the Riemann integral x ¯ h(¯ x) = 0 (f (u) − g(u)) du, where h(¯ x) represents the cumulative earthwork mass at x ¯. In Figure 1.1, we show an example of a vertical alignment with the corresponding mass diagram beneath. From Definition 1.1 and Figure 1.1, we observe the following characteristics. 5 1.1. Road design − If h′ (¯ x) < 0, then x ¯ is the position of a fill. If h′ (¯ x) > 0, then x ¯ is the position of a cut. − Any extremum of h(¯ x) represents a switch from cut to fill, or from fill to cut. − The net accumulation between a station x ¯i and a station x ¯j , where x ¯i < x ¯j is h(¯ xj ) − h(¯ xi ). − For any two stations x ¯i and x ¯j , if h(¯ xi ) = h(¯ xj ), then the net accumulation between the two stations is zero. These points are called balance points, since they show points of balance between cut and fill volumes. − A consequence of the last point is that if h(0) = 0 and x ¯end is the last station of the road, then an optimal alignment with respect to the balance of earth quantities has h(¯ xend ) = 0. With the above observations, we can derive a rule of thumb to determine the hauling directions. Remark 1.2. Let h(0) = 0. If h(¯ x) > 0 and x ¯i < x ¯ < x ¯j , where x ¯i and x ¯j are balance points, then earth needs to be hauled ahead. The distance between x ¯i and x ¯j is the maximum push or haul distance. If h(¯ x) < 0 and x ¯i < x ¯<x ¯j , then earth needs to be hauled back. The values at the extrema represent the volumes of material that must be hauled along the road. The mass diagram and the above rule of thumb have several shortcomings. − The mass diagram does not take into account haul costs that are not proportional to the haul distance. − The mass diagram does not take into account different soil characteristics. − The mass diagram does not take into account borrow pits or waste pits that can be used to borrow additional quantities of soil or dispose of excessive material at off-the-roadway sites. − The mass diagram does not take into account blocks like trees, rocks, and rivers that are present in the haul route. To address some of the above shortcomings, Mayer and Stark [MS81] presented a LP model that not only balances the earth quantities, but also 6 1.2. Basic linear programming model minimizes the hauling cost and incorporates soil material distribution, borrow pits, and waste pits. Based on [MS81], Easa developed a mixed integer linear program (MILP) that uses non-constant unit costs [Eas87], and a quadratic program (QP) that uses linear unit costs [Eas88a]. Jayawardane and Harris [JH90] extended the model in [MS81] by incorporating project duration and dynamic selection of borrow and waste pits. In this thesis, we extend a basic version of the model by Mayer and Stark, in order to overcome the problem of blocks in the haul route. 1.2 Basic linear programming model In order to better explain our model of earth-moving, it is helpful to begin by considering the LP model in Mayer and Stark [MS81]. An extensive description of the model is given in [MS81] and [JH90]. We give a shorter outline in this section with some changes to the notation and refer the reader to the original paper for further details. 1.2.1 Notation For the remainder of this thesis, the cardinality of a set S is denoted |S|. The same notation |·| is used to indicate the absolute value and the meaning of the notation should be clear from the context. 1.2.2 Model parameters In this section we provide information on the basic LP model parameters. Let a section be the region between two fixed stations along the horizontal axis. We divide the vertical alignment into ms sections of equal length. Consider Definition 1.1. If for a section i, x ¯i represents the station at the beginning of the section, and x ¯i+1 represents the station at the end of the section, then x ¯i+1 h(¯ x) = x ¯i (f (u) − g(u)) du is the required change in volume between the ground profile and the road profile in the section. The required change in volume for a section is positive if the section represents a cut, negative if the section represents a fill. As in [MS81, Eas87, JH90], we introduce mb borrow pits and mw waste pits in order to obtain a feasible solution in the case of an unbalanced vertical alignment. A borrow pit or a waste pit has a volume that represents the capacity. The following input parameters are therefore given: 7 1.2. Basic linear programming model section road profile ground profile 1 2 3 4 5 6 Figure 1.2: Example of vertical alignment with sections − ms ∈ Z+ is the number of sections in the road; − mb ∈ Z+ is the number of borrow pits; − mw ∈ Z+ is the number of waste pits. To simplify notation, we introduce the following index sets: − S = {1, . . . , ms } is the set of indices that is used to index sections; − B = {ms + 1, . . . , ms + mb } is the set of indices that is used to index borrow pits; − W = {ms + mb + 1, . . . , ms + mb + mw } is the set of indices that is used to index waste pits; − M = S ∪ B ∪ W. Finally, we need the required change in volume and the moving costs, which are defined by the following parameters: − Vi ∈ R is the required change in volume of a section i ∈ S in cubic units, or the capacity of a borrow pit i ∈ B or waste pit i ∈ W; − cij ∈ R++ ∪ {+∞} is the cost of moving one cubic unit of material from i to j, where i, j ∈ M, i = j. In Figure 1.2, we show the alignment from Figure 1.1 that is divided into ms = 6 sections of equal length. We denote the distance between section i and j as dist(i, j), and evaluate this using dist(i, j) = |¯ xi − x ¯j |, (1.1) where x ¯i is the midpoint of section i and x ¯j is the midpoint of section j. 8 1.2. Basic linear programming model Remark 1.3. The distance dist(i, j) can be defined arbitrarily. A more sophisticated definition than Equation (1.1) would incorporate the distance between the sections on the horizontal alignment as well as the difference of altitude. We leave this as a future research topic. Note that there are many different ways to define the cost cij of moving earth. Examples of non-constant unit costs are given in [Eas87, Eas88a]. We use the following definition. Definition 1.4. The cost cij ∈ R++ ∪ {+∞} to move one cubic unit of material from section i to section j, where i, j ∈ M and i = j, is defined as cexc + dist(i, j)chau + cemb , if i, j ∈ S, c + dist(i, j)c bor hau + cemb , if i ∈ B, j ∈ S, (1.2) cij = cexc + dist(i, j)chau + cwas , if i ∈ S, j ∈ W, +∞, otherwise, where cexc > 0 is the unit cost of excavation, chau > 0 is the cost to haul one cubic unit of material over a unit distance, cemb > 0 is the unit cost for embankment, cbor > 0 is the unit cost to borrow material from a borrow pit, and cwas > 0 is the unit cost to dump material in a waste pit. The following lemma follows clearly from Definition 1.4. Lemma 1.5. Consider any three sections i, j, k ∈ S. If i = k, i = j, and j = k, then the strict triangle inequality cij < cik + ckj holds whenever cij < +∞. The hauling distance from or to borrow pits or waste pits needs to account for dead-haul distances (the distance of the road from the pit-joining section to the pit). Excavation costs on a section are usually smaller than borrowing costs from a pit. This is usually balanced by the lower embankment costs of borrowed material. Also material that is dumped into a waste pit does not need a lot of embankment work and the embankment cost to dump material into a pit is much smaller than the embankment costs on a section. 1.2.3 Decision variables The decision variable xij ∈ R represents the number of cubic units of earth that will be hauled from section i to section j, where i, j ∈ M and i = j. 9 1.2. Basic linear programming model 1.2.4 Objective and constraints The minimization of the earth-moving costs is modeled by LP (1.3) below min cij xij (1.3a) i,j∈M i=j s.t. j∈M j=i j∈M j=i j∈M j=i (xij − xji ) = Vi for all i ∈ S, (1.3b) xij ≤ Vi , for all i ∈ B, (1.3c) xji ≤ Vi , for all i ∈ W, (1.3d) for all i, j ∈ M, i = j. (1.3e) xij ≥ 0, Objective function The objective function in (1.3a) minimizes the total hauling cost. Illogical moves (e.g., moving earth into a borrow pit) are penalized in Equation (1.2) of the cost Definition 1.4 by +∞. In practice, to reduce memory size, decision variables that represent illogical moves are not implemented. Constraints For every road section i that has either a demand or a supply, we want the amount of incoming and outgoing earth movements to balance the demand or supply of the section. Constraints (1.3b) balance the demand and supply between sections and are called balance constraints. The model has |S| balance constraints. We also need to account for the capacity of borrow and waste pits. We enforce the capacity constraints with inequalities (1.3c) and (1.3d). Since one constraint is needed for each borrow and waste pit, we have |B| + |W| capacity constraints. Finally, we can only move earth quantities that are greater or equal to zero. This is ensured by Equation (1.3e). 10 1.2. Basic linear programming model Observations Remark 1.6. Note that shrink and swell factors that are mentioned in [MS81] are not included in this thesis since the model works with banked units and the mentioned factors are only taken into account to compute haul and compaction costs. Proposition 1.7. Let Vs = i∈S Vi , Vb = i∈B Vi , and Vw = i∈W Vi . If Vs ≥ 0 and Vw − Vs > 0, then the LP (1.3) has always a solution. If Vs < 0 and Vb + Vs > 0, then the LP (1.3) has always a solution. Proof. The conditions Vw − Vs > 0 or Vb + Vs > 0 ensure that a feasible point exists, as the waste or borrow pit provide sufficient slack. Since for all i, j ∈ M, cij ∈ R++ ∪ {+∞} and xij ≥ 0, we know that min i,j∈M cij xij ≥ 0. i=j Hence, the LP (1.3) is bounded below. Since there exists a feasible point, and the LP is bounded below, there exists a solution to the LP (1.3). Proposition 1.8. Let x∗ be a feasible solution to the LP (1.3). If x∗ik > 0 and x∗kj > 0, where i = k and k = j, then x∗ is not an optimal solution. ˜ij = x∗ij + x∗ik . ˜ik = 0, x ˜kj = x∗kj − x∗ik , and x Proof. If x∗ik ≤ x∗kj , consider x Then Vi = x∗ij + x∗ik + s∈S s=i,j,k x∗is − x∗si = x ˜ij + x ˜ik + s∈S s=i,j,k s∈S s=i x∗is − x∗si , s∈S s=i (1.4) Vk = x∗kj − x∗ik + s∈S s=j,k x∗ks − s∈S s=i,k ˜kj − x ˜ik + x∗sk = x s∈S s=j,k x∗ks − x∗sk , s∈S s=i,k (1.5) Vj = −x∗ij − x∗kj + s∈S s=j x∗js − s∈S s=i,j,k x∗sj = −˜ xij − x ˜kj + s∈S s=j x∗js − x∗sj , s∈S s=i,j,k (1.6) and the volume constraints are satisfied. But from Lemma 1.5 we know that cik + ckj > cij , and thus the total cost is cheaper, because ˜ik +ckj x ˜kj +cij x ˜ij . cik x∗ik +ckj x∗kj +cij x∗ij > (cij −ckj )x∗ik +ckj x∗kj +cij x∗ij = cik x 11 1.2. Basic linear programming model ˜kj = 0, and x ˜ij = x∗ij + x∗kj . Then, ˜ik = x∗ik − x∗kj , x If x∗ik > x∗kj , consider x the same arguments show that equations (1.4), (1.5), and (1.6) hold, and thus the volume constraints are satisfied, while the total cost is cheaper. Proposition 1.9. The LP (1.3) contains |M|2 −|M| decision variables and |M|2 constraints. Proof. The number of moves from any section i ∈ M to any other section j ∈ M is equivalent to the number of pairs of section indecies, (i, j), i, j ∈ M, that are formed by the Cartesian product M × M. Now |M × M| = |M|2 . Since we do not count moves where i = j, we have |M|2 − |M| decision variables. The constraint in (1.3b) appears for every section in S. Also, constraints (1.3c), and (1.3d) appear for each borrow pit in B and each waste pit in W, respectively. Constraint (1.3e) appears for each of the |M|2 − |M| variables. It follows that we have |S| + |B| + |W| + |M|2 − |M| = |M|2 constraints. Remark 1.10. In practice, since we do not implement illogical moves, the number of decision variables and constraints are less than the number given in Proposition 1.9. If we count only logical moves, we have |S|2 −|S| potential moves between sections, |S||B| potential moves between sections and borrow pits, and |S||W| potential moves between sections and waste pits. Thus in the actual implementation of the model, we have |S|(|M| − 1) decision variables. Accordingly, we have |S| + |B| + |W| + |S|(|M| − 1) = |M|(1 + |S|) − |S| constraints in the actual implementation. In the next chapter we will use these results in our discussion of the extension of the basic LP model (1.3). 12 Chapter 2 Mixed integer programming model He has barred my way with blocks of stone; he has made my paths crooked. Lamentations 3:9 (Holy Bible, NIV) The model presented in Section 1.2 is a major improvement over the typical mass diagram due to the minimization of hauling distances and the incorporation of borrow and waste pits. However, in practice, a move from one section to another might not be possible due to natural blocks between the two sections. In this thesis, a block is defined as follows. Definition 2.1. A block is an obstacle like a river, a group of trees, a rock-bed, or some other topographical feature that needs to be removed or otherwise dealt with before earth can be moved over the location of the block. In a traditional earthwork construction process, clearing of blocks is usually treated as a separate task before excavation and embankment occurs [Trb90]. In such a case, the use of a mass diagram or a LP model, as the one described in Section 1.2, is perfectly valid, since after the obstacle removal, all the sections are accessible. This, however, does not account for earthwork that may occur during the obstacle removal (i.e. filling of a ditch or blasting of rock beds). During obstacle removal, stockpiles may be created and material may be borrowed elsewhere. This makes it difficult to assess the overall costs of obstacle removal, excavation, embankment, and haul in a traditional setting. Moreover, the order of removal of blocks might affect the final cost of the earthwork. We therefore extend the model in [MS81] with the introduction of blocks. 13 2.1. Notation 2.1 Notation For the remainder of this thesis, the symbol ∧ is a logical connective representing and, and the symbol ∨ is a logical connective representing or. Hence A ∧ B means “A and B”, where A ∨ B means “A or B” 2.2 Model parameters We extend the basic LP formulation in Section 1.2 to develop a model that incorporates nb blocks. Definition 2.2. Let I = {1, . . . , nb }. Each block is assigned to a road section. We let ς: I → S be the mapping that maps the block index to the assigned section index. Also, we add the access roads that lead to the construction site of the road. Each access road is connected to one of the ms sections. Access roads are typically attached to sections at the beginning and end of the road, but may also occur at any point along the road and in between blocks. The number and position of access roads influences the optimal schedule for the removal of the blocks. Definition 2.3. Let R = {1, . . . , nr } be the set of all access roads. Each access road is assigned to a road section. We let ̺: R → S be the mapping that maps the access road index to the assigned section index. Since each access road allows for transporting earth out to an external waste pit or importing earth from an external borrow site, we attach an external borrow and waste pit to every access road with infinite capacity each. These pits are accounted for in set B and W. We introduce discrete time steps to remove the blocks. The removal of a block requires one time step and the block is removed at the end of the time step. Multiple blocks may be removed within the same time step. In a worst-case scenario, where all the blocks need to be removed in sequential order, nb time steps are needed to remove the blocks. Once all the blocks are removed, we may need an additional time step to finish all the earthwork. We therefore introduce nb + 1 time steps t with t ∈ T = {0, . . . , nb }. 14 2.3. Decision variables The cost coefficients from Equation (1.2) remain the same in our model. However, we need a cost coefficient for a move for every time step. Hence, for each t ∈ T we set cijt = cij . Excavation and embankment costs on external pits should be higher than for sections. Remark 2.4. If more than one access road is present, then any movement of earth between two sections could also go over the access roads. One possible way to allow such moves is to represent each access road as a section. Each section that has an access road attached could be cloned and the clone be added to the set S. In this thesis, however, we do not account for such moves and we leave it as a further research topic. In Section 1.2.2, we stated that the temporary stocking of earth between sections is avoided due to Lemma 1.5. With the introduction of blocks, however, this is not guaranteed anymore (as we demonstrate in Example 2.7). We therefore introduce for each section i a stockpile tolerance parameter ǫstk i ∈ R+ . The significance and usage of this parameter is shown in Section 2.4.5. 2.3 Decision variables Similar to the cost coefficients, we keep the same variables representing the quantities of earth that are moved between sections. However, since we now have the flexibility of moving earth between any two sections at any time step, we multiply the variables in Section 1.2.3 by the number of time steps. The variable xijt ∈ R represents the quantity of earth that is moved from section i to section j at time step t. A block k examined at time step t is represented with a binary variable. A zero value for a block shows that no or only partial work was done to remove the block during that time step, and a one indicates that the block was cleared at the end of time step t. Hence for block k with k ∈ I, t ∈ T , we define the variable that indicates the removal of block k at time t via ykt ∈ {0, 1}, such that if block k is still present at time t, then ykt = 0. 15 2.4. Objective and constraints 2.4 Objective and constraints The minimization of the earth-moving costs, including the removal of blocks, can therefore be stated as MILP (2.1) below cijt xijt min (2.1a) t∈T i,j∈M i=j s.t. (xijt − xjit ) = Vi for all i ∈ S, (2.1b) xijt ≤ Vi , for all i ∈ B, (2.1c) xjit ≤ Vi , for all i ∈ W, (2.1d) t∈T j∈M j=i t∈T j∈M j=i t∈T j∈M j=i u (xς(k)jt − xjς(k)t ) ≥ yku Vς(k) for all k ∈ I, u ∈ T , Vς(k) ≥ 0, t=0 j∈M j=ς(k) (2.1e) u (xς(k)jt − xjς(k)t ) ≤ yku Vς(k) for all k ∈ I, u ∈ T , Vς(k) < 0, t=0 j∈M j=ς(k) (2.1f) xijt ≤ M yk,t−1 , for all t ∈ T \ {0}, i, j ∈ M, k ∈ I, (i ≤ ς(k) ≤ j) ∨ (j ≤ ς(k) ≤ i) (2.1g) xijt ≤ M yk,t−1 + M yk+1,t−1 , for all t ∈ T \ {0}, i, j ∈ M, k ∈ I \ {nb }, r ∈ R, ς(k) ≤ i, j ≤ ς(k + 1) ∧ {r : ς(k) ≤ ̺(r) ≤ ς(k + 1)} = ∅ (2.1h) u (xijt − xjit ) ≤ Vi + ǫstk i , t=0 j∈M j=i for all i ∈ S, u ∈ T , Vi ≥ 0 (2.1i) 16 2.4. Objective and constraints u (xijt − xjit ) ≥ Vi − ǫstk i , t=0 j∈M j=i for all i ∈ S, u ∈ T , Vi < 0 (2.1j) u t=0 k∈I yki ≥ u, for all u ∈ T , (2.1k) xijt ≥ 0, for all i, j ∈ M, i = j, t ∈ T , (2.1l) ykt ∈ {0, 1}, for all k ∈ I, t ∈ T , (2.1m) where M ∈ R is a large constant that exceeds the quantity of any possible earth move. 2.4.1 Objective function As in the LP model, the MILP objective function (2.1a) minimizes the total hauling costs, but this time over all the time steps. Note that the cost coefficients for the block removal variables are set to zero. 2.4.2 Volume constraints The balance constraints (2.1b) and the capacity constraints (2.1c) and (2.1d) have the same form as in the LP model, with an additional summation over all the time steps. 2.4.3 Block removal indicator The decision variable ykt , that indicates the removal of a block k at time step t, is binary. If we want to find an optimal removal schedule, we need to set some of these removal indicators to one. Once we decided to remove block k at time step t, we let ykt = 1. The removal of block k in a section ς(k) may involve earthwork. In order to effectively enforce the removal of the block at time step t, we need to satisfy any potential demand or supply of earth in this section at time t. For a supply, this is done by Equation (2.1e), for a demand, by Equation (2.1f). Therefore, the block removal indicator constraints relate the binary variables with the continuous variables. Remark 2.5. Note that once a block k is removed at time step t, there is no need to enforce the monotonicity ykt ≤ yk,t+1 in order to represent the status of a block in the later time steps. If we want to move earth over the 17 2.4. Objective and constraints block at some later time step u > t, we can just set the variables yk,u−1 to 1 and the corresponding removal indicator constraint will already be satisfied. In fact, numerical tests showed that enforcing the monotonicity slows down the solving process. 2.4.4 Move feasibility constraint Consider a move from section i to j. If there is a block k, such that i ≤ ς(k) ≤ j, then we need to remove the block k before the move can happen. Therefore, if yk,t−1 = 0, then the move is not possible. However, if yk,t−1 = 1, then any (reasonable) quantity can be moved. In order to indicate any reasonable quantity, we use a very large constant M as upper bound for such a move. Equation (2.1g) models this relationship. Consider a move from section i to j. If there are blocks k and k + 1, such that ς(k) ≤ i, j ≤ ς(k + 1), and there is no access road between ς(k) and ς(k + 1), then we need to remove either block k or block k + 1 first, before the move can happen. This constraint is modeled in (2.1h). The requirement that there is no access road in between ς(k) and ς(k+1) is shown in the condition of (2.1h) as {r : ς(k) ≤ ̺(r) ≤ ς(k + 1)} = ∅. 2.4.5 Stockpile constraints Remark 2.6. Proposition 1.8 from the LP model does not apply to the block removal MILP model in (2.1). An explanation follows from Example 2.7 below. Example 2.7. Consider the example with sections and volumes given in Table 2.1. The road section 3 represents a block with demand of one cubic unit of earth. Road section 4 has a supply of one, but is not accessible due to the block at section 3. The volume for road sections 1 and 2 is zero. There is one access road at the beginning of the road that joins section 1. For simplicity, let us consider cexc = chau = cemb = cbor = cwas = 1. Since Table 2.1: Road sections and volumes (bold font represents a block, ≈ an access road). Section i Volume Vi ≈1 0 2 0 3 -1 4 1 there is nb = 1 block, we have 2 time steps. It is easy to see that the optimal solution to this problem is to first remove the block in section 3 by borrowing 18 2.4. Objective and constraints material from the adjacent section 2, followed by filling the newly created demand in section 2 with the cut of section 4. Hence the solution is x230 = 1 x421 = 1 y10 = 1, which is illustrated in Figure 2.1. The squares represent the road sections. The time steps are shown as the space between two layers of squares. A black colored square indicates the presence of a block, a gray colored square shows a block that was removed in a previous step. The arrows in between the layers represent earth movements. With the given cost model, any other 1 2 3 4 1 1 2 3 4 1 1 2 3 4 Figure 2.1: Earthwork scheme with removal schedule (black squares indicate blocks, gray squares indicate removed blocks). solution would be more expensive. If we disregard the time steps and add the quantities for each move over time, that is xij = xij0 + xij1 , for all i, j ∈ {1, 2, 3, 4}, then we have x42 = 1 x23 = 1 and the solution moves one cubic unit of earth from section 4 to 2, and from section 2 to 3. Hence, in this example, Proposition 1.8 does not apply to the MILP in (2.1). 19 2.4. Objective and constraints From Remark 2.6, we conclude that a constraint is required to limit the amount of earth that can be stocked or borrowed from a road section. This can be enforced by a similar constraint like the balance constraint (2.1b), where the equality is replaced by a bound that represents the volume of a section plus or minus the stockpile tolerance. However, since stockpiling can occur at any time step, we need the constraint to hold at every time step. The resulting equations are (2.1i) for an upper bound, and (2.1j) for a lower bound. 2.4.6 Block removal enforcement In order to provide the MILP solver with as much information as possible, we add the following constraint. At each time step t = 1, 2, . . . , nb , we want to remove at least one block. This is guaranteed with constraint (2.1k). This constraint also ensures that no more than nb + 1 time steps are needed. 2.4.7 Observations Proposition 2.8. Let Vs = i∈S Vi , Vb = i∈B Vi , and Vw = i∈W Vi . If Vs ≥ 0 and Vw − Vs > 0, then the MILP (2.1) has always a solution. If Vs < 0 and Vb + Vs > 0, then the MILP (2.1) has always a solution. Proof. The conditions Vw − Vs > 0 or Vb + Vs > 0 ensure that constraints (2.1b) to (2.1f), and (2.1i) to (2.1j) can be satisfied, since the waste or borrow pit provide sufficient slack. Since for all i, j ∈ M, t ∈ T , cijt ∈ R++ ∪ {+∞} and xijt ≥ 0, we know that min t∈T i,j∈M cijt xijt ≥ 0. i=j In Section 4.1, we show how to set the binary variables ykt in order to obtain a feasible point. Furthermore, in Section 3.2, we show that there are a finite number of ways to select the binary variables, and therefore a finite number of feasible points to the MILP (2.1). Altogether, since the MILP (2.1) has a finite number of feasible points that are bounded below, and there exists at least one such point, there exists a solution to the MILP (2.1). Proposition 2.9. The MILP (2.1) contains (|I| + 1)(|M|2 − |M|) continuous, and (|I| + 1)|I| binary decision variables. 20 2.4. Objective and constraints Proof. In order to get the number of continuous variables, we multiply the |M|2 − |M| variables from Proposition 1.9 with |I| + 1 time steps. We have therefore (|I| + 1)(|M|2 − |M|) continuous variables. Since there are |I| blocks and |I| + 1 time steps, we have (|I| + 1)|I| binary variables. Altogether, we have (|I| + 1)((|M|2 − |M|) + |I|) variables. Remark 2.10. As in Remark 1.10, the number of decision variables in the practical implementation of the model is smaller than the number given in Proposition 2.9. If we multiply the number of continuous variables from Remark 1.10 with |I| + 1 time steps, we get |S|(|M| − 1)(|I| + 1) continuous variables. The (|I|+1)|I| binary variables remain the same as in Proposition 2.9. In the next section we examine some of the properties of Model (2.1). 21 Chapter 3 Model analysis Aut inveniam viam aut faciam.2 Hannibal (247–183 BCE) 3.1 Space complexity In this section, we discuss the space complexity of Model (2.1). In particular, we look at the growth rate of the model size, with respect to the input parameters, that is required to store the model in the memory of a computer. In this context, we employ the usual asymptotic notation for space complexity. The following definition is a reminder of the most common asymptotic bounds that we are using. Definition 3.1 (Cormen et al. [CLRS09, p. 43]). For a given function g(n), we denote by Θ (g(n)) the set of functions Θ (g(n)) = {f (n) : there exist c1 > 0, c2 > 0, and n0 ∈ Z++ such that 0 ≤ c1 g(n) ≤ f (n) ≤ c2 g(n) for all n ≥ n0 }, and we say that g(n) is an asymptotic tight bound for f (n). We denote by O (g(n)) the set of functions O (g(n)) = {f (n) : there exist c > 0, and n0 ∈ Z++ such that 0 ≤ f (n) ≤ cg(n) for all n ≥ n0 }, and we say that g(n) is an asymptotic upper bound for f (n) and Θ (g(n)) ⊆ O (g(n)). Lemma 3.2. The growth rate of the total number of variables for the model in (2.1) is of order Θ |I||M|2 + |I|2 . 2 I shall either find a way or make one. 22 3.1. Space complexity Proof. From Proposition 2.9 we know that there are (|I| + 1)(|M|2 − |M|) + (|I| + 1)|I| = |I||M|2 + |I|2 + |I| − |I||M| − |M|2 − |M| variables. Let |I| = 1, |M| = 2. Clearly, 2(|I||M|2 + |I|2 ) ≥ |I||M|2 + |I|2 + |I| − |I||M| − |M|2 − |M|, and 1 (|I||M|2 + |I|2 ) ≤ |I||M|2 + |I|2 + |I| − |I||M| − |M|2 − |M|. 2 We therefore conclude that the growth rate for the size of the total variables is Θ |I||M|2 + |I|2 . Lemma 3.3. The growth rate of the number of balance constraints (2.1b), borrow pit capacity constraints (2.1c), and waste pit capacity constraints (2.1d) is of order Θ (|M|). Proof. Each section needs a balance constraint (2.1b). Hence there are |S| balance constraints. Each borrow pit and each waste pit needs a capacity constraint. Thus there are |B| + |W| capacity constraints. Altogether, there are |S| + |B| + |W| = |M| balance and capacity constraints. The growth rate of the balance constraints is therefore of order Θ (|M|). Lemma 3.4. The growth rate of the number of removal indicator constraints (2.1e) and (2.1f) is Θ |I|2 . Proof. We have one block removal indicator constraint for each indicator variable. From Proposition 2.9 we know that the number of indicator variables is |I|(|I| + 1). Hence the number of indicator variable constraints is |I|(|I| + 1) which is in Θ |I|2 . Lemma 3.5. The growth rate of the number of feasible move constraints (2.1g) and (2.1h) is of order O |M|2 |I|2 23 3.1. Space complexity Proof. We consider the (|M|2 − |M|) potential moves between all the sections over |I| + 1 time steps. Each move needs at most a total of |I| feasible move constraints. If the move occurs over one or more blocks, we use constraint (2.1g) at most |I| times. If the move occurs in between two blocks, we use the constraint (2.1h) at most |I| times. If the move occurs over n blocks and between two blocks at the same time, we use constraint (2.1g) n times, and constraint (2.1h) at most |I| − n times. Altogether, we use a total of at most |I| constraints. We have therefore (|M|2 − |M|)(|I| + 1)|I| feasible move constraints, which is of order O |M|2 |I|2 . Lemma 3.6. The growth rate of the number of stockpile constraints (2.1i) and (2.1j) is of order O (|M||I|). Proof. At each time step, we need two stockpile constraints for each section. Hence we have 2|T ||S| = 2(|I| + 1)|S| ≤ 2|I||M| + 2|M| stockpile constraints which is in O (|M||I|). Lemma 3.7. The growth rate of the number of removal enforcement constraints (2.1k) is of order Θ (|I|). Proof. At each time step, we need a constraint that sets the least number of blocks that must be removed at that time step. We therefore have |T | = |I| + 1 constraints which is in Θ (|I|). Lemma 3.8. The growth rate of the number of non-negativity constraints (2.1l) is of order Θ |M|2 |I| . Proof. We need a non-negativity constraint (2.1l) for each of the continuous variables. From Proposition 2.9, we know that there are (|I|+1)(|M|2 −|M|) continuous variables, which is in Θ |M|2 |I| . Lemma 3.9. The growth rate of the number of binary constraints (2.1m) is of order Θ |I|2 . Proof. We need a binary constraint (2.1m) for each of the binary variables. From Proposition 2.9, we know that there are (|I| + 1)|I| binary variables, which is in Θ |I|2 . 24 3.1. Space complexity Proposition 3.10. The growth rate of the memory size for the MILP described in (2.1) is of order O |M|7 . If the number of blocks |I| is fixed, then the memory size for the MILP (2.1) is of order O |M|4 . Proof. We add the asymptotic bounds for all the constraints from Lemma 3.3 to Lemma 3.9. The growth rate of the sum of all the constraints is therefore Θ (|M|) + Θ |I|2 + O |M|2 |I|2 + O (|M||I|) + Θ |M|2 |I| + Θ |I|2 , (3.1) which is of order O |M|2 |I|2 + |I|2 . From Lemma 3.2, we know that the number of variables is of order Θ |I||M|2 + |I|2 . Hence, the size of the problem is bounded by O |M|2 |I|2 + |I|2 × Θ |I||M|2 + |I|2 which is of order O |M|2 |I|4 + |M|4 |I|3 . Since |B| ≥ 1, and |W| ≥ 1, we know that |I| < |M|, and hence |I| ∈ O (|M|). The asymptotic upper bound for the memory size is therefore O |M |7 . If we fix the number of blocks |I|, the dimension of the problem is bounded by O |M|2 × Θ |M|2 , which is of order O |M|4 . Remark 3.11. The asymptotic bound of O |M |7 in Proposition 3.10 is a worst-case bound, where the number of blocks |I| could grow linearly with |M|, the sum of all sections, borrow pits, and waste pits. We mentioned in Chapter 2 that |I| is typically between 5 and 20. Therefore if |M| is large, then |I| is usually small compared to |M|. In that case, we can fix |I| and, as shown in Proposition 3.10, the memory size of the MILP (2.1) grows with order O |M|4 . Regardless of the fact that we fix the number of blocks in practice, the memory growth rate in O |M|4 remains a concern for roads with hundreds of sections. Most MILP solvers take advantage of the sparsity of (2.1) and the bound may not be attained. However, in our numerical tests, we still observed a growth rate that quickly exhausted the memory of a standard desktop computer. From the given bound and from practical experience using a 64 bit desktop computer with 8 GB of RAM and the Gurobi Optimizer 2.0 MILP solver, we recommend to keep the number of sections |S| ≤ 200. 25 3.2. Schedule search 3.2 Schedule search In Section 3.1, we showed that the memory size grows in polynomial time with the input. However, the difficult part in our model is the scheduling for the block removal. In this section we analyze the scheduling problem more closely. In particular, we look at the schedule search space and the set of feasible schedules in the schedule search space. For all time steps t ∈ T , we list the variables ykt as a binary string in the order of block position. Since |T | = |I| + 1, we have |I| + 1 binary strings of length |I|. At time t = 0, all the blocks are in place, thus yk0 = 0, and the binary string contains zeros only. At a time step t > 0, the corresponding binary string contains t or more ones that represent the removal of blocks. A binary string that contains ones only, represents a complete removal. Definition 3.12. A schedule is a set of |I| + 1 binary strings where each string belongs to a distinct time step. A road with five blocks, for example, has an initial binary string of 00000. If one block is removed, say the first one to the left, then we have a new string which is 10000. We continue the removal until the binary string is 11111. Note that complete removal of all blocks may occur before time step t = |I| + 1. In that case, the digits of the remaining |I| + 1 − t strings can be considered all 1 for all remaining time steps (the constraints affected by the binary variables for the remaining time steps are already satisfied). Definition 3.13. The schedule search space is the number of potential schedules in the MILP (2.1) as defined by the binary selection. This includes infeasible schedules. Lemma 3.14. The schedule search space is asymptotically bounded by 2 O 2|I| 2|I| . Proof. Since there are |I| + 1 binary strings of length |I|, we have 2|I|(|I|+1) potential schedules. It follows that the schedule search space is in 2 O 2|I| 2|I| . 26 3.2. Schedule search Table 3.1: Road sections and volumes (bold font represents a block, ≈ an access road). Section i Volume Vi 1 2 2 -2 3 -1 4 -1 ≈5 5 6 -1 7 -2 00000 00100 00110 00111 01111 11111 00000 10000 10001 11001 11011 11111 (a) Feasible (b) Infeasible Figure 3.1: Potential schedules Not all schedules in the schedule search space are feasible. Blocks can only be removed if they have an access road attached or if an adjacent block that is accessible is removed first. We are interested in the number of feasible schedules that lead to a complete removal. Definition 3.15. The feasible schedules are the schedules that only remove a block if the block is accessible at the appropriate time. Hence, if k ∈ I is a block at time step t ∈ T , then ykt = 1 implies that either there exists an access road r ∈ R, such that ς(k − 1) ≤ ̺(r) ≤ ς(k) ∨ ς(k) ≤ ̺(r) ≤ ς(k + 1), or we have yk−1,u = 1 ∨ yk+1,u = 1, for some u ∈ T , u < t. Example 3.16. Consider the example with the sections given in Table 3.1. Sections 1, 3, 4, 6, and 7 are blocks. There is a single access road that joins section 5. Figure 3.1 shows two potential schedules. The schedule in Figure 3.1a is feasible, with a possible earthwork scheme shown in Figure 3.2. However, the schedule in Figure 3.1b is not feasible, since the blocks at either end of the road are not accessible. There are other blocks between the access road and the end of the road. Blocks can be removed sequentially or simultaneously. From an access road with blocks on either side, there are three possible ways to remove 27 3.2. Schedule search 1 2 3 4 5 6 7 1 1 2 3 4 5 6 7 1 1 2 3 4 5 6 7 2 1 2 3 4 5 6 7 1 1 2 3 4 5 6 7 2 1 2 3 4 5 6 7 Figure 3.2: Earthwork scheme with removal schedule. 28 3.2. Schedule search 0000 0100 0110 0010 1100 1110 0110 1110 1111 0111 0110 0111 0011 1110 1111 1110 1111 0111 1111 1111 1110 1111 0111 1111 0111 1111 1111 1111 1111 1111 1111 Figure 3.3: Tree structure of all feasible block removal schedules for 4 blocks with an access road joining between the two middle block adjacent blocks. We can remove the block to the left, to the right, or both at once. If the access road is attached to the block, the block needs to be removed first. Using edges to represent left, right, or simultaneous removal, we can draw a tree that contains all the feasible schedules for a given configuration of access points. An example for 4 blocks with an access road attached to a section between the two center blocks is shown in Figure 3.3. We start by approximating the number of all feasible schedules for a single access road from below. We therefore restrict the set of feasible schedules by analyzing schedules that allow left and right removal only, using one single access road. Let the access road be at a section between block r and block k = |I| − r. Consider the (k × r)-lattice of integral points in Z2 with coordinate (0,0) at the lower left corner and coordinate (k,r) at the upper right corner. We associate a left block removal to a step upward and a right block removal to a step to the right. We call the vertical steps (1,0)-steps and the horizontal steps (0,1)-steps to represent the increase of the y− and x−axis respectively. A lattice path consisting of (0,1)- and (1,0)-steps that starts at (0,0) and ends at (k,r) is in bijection with a schedule starting at the root node with a zero binary string and ending at a leaf with a binary string of ones. Figure 3.4a depicts such a path. The next Lemma provides a formula to count all the possible lattice paths starting at (0,0) and ending at (k,r). We provide the proof as outlined in [Aig07]. Lemma 3.17 (Aigner [Aig07, p. 17]). Consider a (k × r)-lattice of integral points in Z2 . The number of paths consisting of (1,0) and (0,1) steps from 29 3.2. Schedule search (6,4) (6,4) k r (0,0) (0,0) (a) (1,0),(0,1)-lattice path (b) Delannoy path Figure 3.4: (6 × 4)-lattice in Z2 , representing 10 block configuration with a single access road between block 6 and 4 (0,0) to (k,r) is L(k, r) = k+r . r (3.2) Proof. Assign the symbol L(eft) to a (1,0)-step and R(ight) to a (0,1)-step. Every lattice path is in bijection with a (k + r)-word over {L, R} containing possible combinations of (k + r) exactly r L’s. This corresponds to k+r r words containing r number of L symbols. Remark 3.18. For the block removal schedules, since k + r = |I|, for |I| blocks and an access road between block k and r, we have |I| r L(k, r) = (3.3) feasible schedules for the restricted set of feasible schedules for a single access road that contain left and right removals only. Remark 3.19. Equation (3.3) is the formula for the binomial coefficients [Aig07, p. 16]. Consider the first five rows of Pascal’s triangle 1 1 1 1 1 1 2 3 4 1 3 6 1 4 1. 30 3.2. Schedule search Let the rows represent the number of blocks |I|, and the columns represent block r, where a single access road is between block k and r. The binomial coefficients in the triangle represent the number of feasible left and right removal schedules. We notice that the central binomial coeffcients along the central axis of the triangle grow the fastest. They correspond to an access road between k and r, where k = r. If k = r, with the formula given in [Wei02, p. 369], then L(r, r) = 2r r = |I| . r Using Stirling’s approximation [Wei02, p. 2868] for the factorial function, we have |I| 4r , for all r ≥ 1. ≤√ r 3r + 1 Since |I| = 2r, we conclude that the number of feasible schedules for a single access road and left and right removal only is in O 2|I| . We now extend our restricted set of feasible schedules by adding another removal step to the schedules with one access. So far we considered leftremoval and right-removal steps only. But if the access road is in between two blocks, we can also remove the left and right block at the same time, as shown in Figure 3.3. Remark 3.19 is therefore a lower bound on the worstcase. Instead of having only (0,1)- and (1,0)-steps in our lattice path, we now add a diagonal (1,1)-step. This results in Delannoy paths as shown in Figure 3.4b. The proof for the next Lemma is already given in [SAV+ 03]. We prove it in a slightly different way. Lemma 3.20. Consider a (k × r)-lattice of integral points in Z2 . The number of Delannoy paths consisting of (1,0),(0,1) and (1,1) steps from (0,0) to (k,r) is k D(k, r) = i=0 k i k+r−i . k (3.4) Proof. Assign the symbol L(eft) to a (1,0)-step, R(ight) to a (0,1)-step, and B(oth) to a (1,1)-step. Every Delannoy path corresponds therefore to a string in the {L, R, B} alphabet. Consider a Delannoy path with i diagonals. Hence the corresponding string contains r−i L’s and k−i R’s. Now project3 the (1,1)-steps onto the horizontal x-axis. The projected {L, R} string has 3 The projection of the diagonals onto a lattice is based on an idea by Jim Nastos. 31 3.2. Schedule search now r − i L’s and k − i + i = k R’s. Using the formula for (1,0),(0,1)-lattice paths in Lemma 3.17, Equation (3.2), we have L(k, r − i) = k+r−i r−i = k+r−i k (3.5) (1,0),(0,1)-paths. However, there are ki possible ways we can choose the diagonals that we project. Also, we have i ∈ {0, . . . , k}. Hence we multiply Equation (3.5) by ki and sum i from 0 to k in order to get the formula in Lemma 3.20. Proposition 3.21. The number of feasible schedules for a road with a single access road on either end of the road is in O (1). Proof. We use Equation (3.4) with k = |I|, r = 0 and obtain D(|I|, 0) = 1. We are now able to compute the exact number of feasible schedules for a given road with one access road. However, a lot of construction sites have more than one access road. We are particularly interested in the very common case of two access roads, one at the start and one at the end of the road. Lemma 3.22. The number of feasible schedules for a road with two access points, one at the start and one at the end of the road, is |I| S(|I|) = j=1 D(|I| − j, j), (3.6) where D is defined in Equation (3.4). Proof. Consider a binary string with |I| zeros that represents the time step at time t = 0. If we remove a block at the start of the road, we set the leftmost digit from 0 to 1. If we remove a block at the end, we set the right-most digit to 1. Likewise, if we simultaneously remove on both sides, we set the outer-most digits to 1. Since there is no access point anywhere else, this procedure continues until all digits are set to 1. The last digit or the last two digits that are set to 1 represent the meeting point of the construction crews, once all the blocks are removed. For any removal schedule with a meeting point j ∈ {1, . . . , |I|}, we can invert the schedule in the sense that zeros are set to one, and ones to zero. A inverted schedule represents a schedule that has the meeting point as a single access road. For a fixed 32 3.2. Schedule search meeting point, there is a bijection between all the schedules with two access roads that meet at the meeting point and the schedules with a single access road attached to the meeting point. In order to get all feasible schedules, we consider all |I| possible meeting points. Hence we sum j over {1, . . . , |I|} and compute the inverse schedules with (3.4), where r = j and k = |I| − j. This leads to the formula in Lemma 3.22. As presented in Banderier and Schwer [BS05], the central Delannoy numbers along the diagonal from (0,0) to any (r,r) grow the fastest. To get a worst-case scenario for a configuration with |I| blocks and one access road, we let the access road join the center of the binary string, thus k = r = |I|/2. In the same paper [BS05], Banderier and Schwer also give an approximation to compute an upper bound on the central Delannoy numbers. It is given as D(k, k) = 5.82842709k (0.57268163k−1/2 − 0.06724283k−3/2 + 0.00625063k−5/2 + O k−7/2 ). (3.7) Proposition 3.23. The number of feasible schedules for a road with a single access road between the center blocks of the road is in O 3|I| . Proof. Consider Equation (3.4). If the single access road is between the center blocks, we have r = k. The dominating term in the sum of Lemma 3.22 is the number for the central Delannoy paths D(k, k) when k = |I|/2. We therefore use Equation (3.7) as an approximation. Since √ 5.83|I|/2 (|I|/2)−1/2 2.42|I| 2|I|−1/2 , we conclude that the number of feasible schedules is in O 3|I| . Proposition 3.24. The number of feasible schedules for a road with two access roads, one on each end of the road, is in O 3|I| . Proof. Use Equation (3.6) with |I| = 0, 1, 2, . . .. We obtain the sequence {0, 1, 2, 5, 12, 29, 70, 169, . . .}. These are the Pell-Lucas numbers [Wei02, p. 2185], given by the closed form √ √ (1 + 2)n − (1 − 2)n √ . P (n) = 2 2 √ √ √ Since P (n) ≤ (1 + 2)n 2/4 ≈ 2.42n 2/4, we conclude that the number of feasible schedules is in O 3|I| . 33 3.2. Schedule search Remark 3.25. Consider Proposition 3.21, 3.23, and 3.24. If there is one access road only and the number of blocks |I| is small, then the number of feasible schedules is clearly tractable. In this case, it is reasonable to enumerate all the feasible schedules and to solve a LP for each of the feasible schedules. The optimal removal schedule is the schedule for the LP solution with the lowest objective value. The enumeration of the feasible schedules can be done with the use of a Depth-First search algorithm [CLRS09]. The Depth-First algorithm runs linearly in the size of the tree. The solving time is therefore exactly proportional to the tree size, or the number of paths. If there is one access road in the middle of the road, then there is an exponential number of feasible schedules. Whenever |I| is large, the enumeration of the feasible schedules is not tractable and one should solve a single MILP instance of the form in (2.1) instead. If there are two access roads, one on either side of the road, then there is an exponential number of feasible schedules. Again, if |I| is large, one should solve the MILP in (2.1). 34 Chapter 4 Algorithms An algorithm must be seen to be believed. Donald Knuth (1938–) Model (2.1) is sufficient to solve the earth-moving problem. In this section we propose algorithmic techniques that improve the speed of convergence. In order to simplify the notation, we formulate (2.1) in standard form. Let m be the number of constraints in (2.1). Let nc denote the number of continuous variables and n = nc + nb denote the total number of variables in (2.1). We rewrite the inequality constraints and group the left-hand sides in a matrix A ∈ {0, ±1}m×n , and the right-hand side in a vector b ∈ Rm , such that we can write (2.1) in the canonical form min cT x s.t. A(x, y)T ≥ b x≥0 y ∈ {0, 1}nb , (4.1) c is the cost vector, and x ≥ 0 means component where x ∈ Rnc , c ∈ Rn++ wise xi ≥ 0, for all i ∈ {1, . . . , nc }. We can solve (4.1) directly as a MILP, using a MILP solver software [Fou05]. Although this is not an algorithm in the common sense, we call it here the direct solve algorithm for future comparison. Algorithm 1: Direct solve Input: x,y,c,A,b (from (2.1)) Output: Global minimum z ∗ = cT x∗ for (4.1) begin Build and solve MILP (4.1) for x∗ ; end 35 4.1. Mass schedule as starting point 4.1 Mass schedule as starting point Definition 4.1. The LP obtained by relaxing the constraints y ∈ {0, 1}nb with y ∈ [0, 1]nb , is called the LP relaxation of the MILP (4.1). Fact 4.2 (Dantzig [Dan98, p. 516]). The LP relaxation is a lower bound on the MILP. Most modern MILP solvers use a branch-and-cut algorithm [Fou05]. Branch-and-cut is a combination of the branch-and-bound method with the cutting-plane algorithm of Gomory [Gom58]. A classical branch-and-bound algorithm solves first the LP relaxation to the MILP. The algorithm then partitions the LP relaxation into subproblems which have one of the previously relaxed integer variables fixed to closest integer lower- and upperbounds. The continuation of this process results in a branch-and-bound tree with nodes that represent subproblems. The bounding part of the algorithm looks for the best complete solution z m in the tree. At each node k, where the solution z k to the subproblem is greater than z m , the algorithm can discard any further subproblem created by this node, since all following solutions will be greater than z m . For a more detailed description of the branch-and-bound algorithm, we refer to [PS98]. A common technique to accelerate the solution process of a MILP is to provide an initial upper bound z ◦ on the solution value z m . This allows for early bounds on the branch-and-bound tree. An even better technique is to provide the solver with a feasible starting point (x◦ , y ◦ ), where z ◦ = cT x◦ . The number of nodes that need to be expanded can be reduced significantly with a starting point (x◦ , y ◦ ) that produces a sufficiently low upper bound z ◦ . In this section, we describe a way to find a reasonable starting point. In a traditional mass diagram procedure, work is performed in a unidirectional way from the start of the road to the end (see Remark 1.2). For a road with access attached to the starting section, a feasible block removal schedule is obtained by removing the blocks sequentially from the start of the road to the end, one block at each time step. If there is no access at the starting section of the road, we select the section with an attached access that is closest to the starting section and first remove sequentially all the blocks on the left of the access, and then all the blocks on the right. Once we have the schedule, we can fix the binary variables in (4.1) and solve the problem as a LP. This will give us feasible values for the continuous variables, and together with the fixed binary variables, a feasible solution to (4.1). 36 4.2. Maximal allowed moving distance dmax Procedure FindFeasiblePoint(x,y,c,A,b) Input: x,y,c,A,b (from (2.1)) Output: Feasible point (x∗ , y ∗ ) for (4.1) begin a ← FindFirstAccess(S); for k ← a to 1 do ∗ yk,a−k+1 ← 1; end for k ← a + 1 to |I| do ∗ ← 1; ykk end for k ← 1 to |I| do for t ← 1 to |I| do ∗ ← 1; ∗ = 1 then ykt if yk,t−1 end end ∗ as LP; Solve (4.1) with fixed ykk end We can now extend Algorithm 1 with procedure FindFeasiblePoint. We call this the direct with x◦ algorithm. Algorithm 2: Direct with x◦ Input: x,y,c,A,b (from (2.1)) Output: Solution z ∗ = cT x∗ for (4.1) begin Build MILP (4.1); (x◦ , y ◦ ) ←FindFeasiblePoint(x,y,c,A,b); Solve MILP (4.1) with (x◦ , y ◦ ) for (x∗ , y ∗ ); end 4.2 Maximal allowed moving distance dmax To reduce the model size, and hence reduce the solving time of the MILP, we introduce a heuristic that eliminates a significant number of variables. In order to do this, we introduce the parameter dmax . Definition 4.3. The Maximal Allowed Moving Distance, dmax ∈ S, is the largest number of sections over which earth is allowed to be moved. There37 4.2. Maximal allowed moving distance dmax fore, for all i, j ∈ M, t ∈ T , xijt = 0, if |i − j| > dmax . (4.2) The motivation behind dmax is that from practical experience, it makes little sense to move earth from a supply at one end of the road to a demand at the other end of the road. Intuition tells us that if there are other supplies along the road in between the two sections, then the closest supply to the demanding section should be chosen. Hence if we disallow moves that go over too long a distance, then we restrict the search space to shorter, and thus probably cheaper moves. We enforce dmax with a set of constraints in the form of Equation (4.2). If there are p such constraints, we can group them in a matrix Rdmax ∈ {0, 1}p×n and reformulate the MILP (4.1) in a dmax -extended form min cT x s.t. A(x, y)T ≥ b (4.3) Rdmax x = 0 x≥0 y ∈ {0, 1}nb . Remark 4.4. It is clear from (4.2) that if dmax is small, then Rdmax becomes large (although very sparse). So instead of reducing the model size, we increase it. However, with each constraint in Rdmax , we can remove a large number of feasible move constraints in A. This would be done before solving the model; potentially while building the model. Therefore, the model size is actually decreased. In fact, not only is the model size reduced, but also the problem becomes easier to solve due to the sparsity of Rdmax . Lemma 4.5. For any optimal solution value z ∗ = cT x∗ to the MILP (4.1) and any optimal solution value zd∗max = cT x∗dmax to the extended form (4.3), the inequality z ∗ ≤ zd∗max is always true. Proof. Let C = {(x, y) : A(x, y)T ≥ b, x ≥ 0, y ∈ {0, 1}nb } and Cdmax = {(x, y) : A(x, y)T ≥ b, Rdmax x = 0, x ≥ 0, y ∈ {0, 1}nb }. Clearly, Cdmax ⊆ C. It follows that z ∗ ≤ zd∗max . From Lemma 4.5 we can see that the optimality of the solution with respect to (4.1) is not guaranteed when using dmax , unless dmax ≥ |S|. For this reason, we call (4.3) the dmax -heuristic to the original problem (4.1). Note that the demands and supplies for a particular vertical alignment are 38 4.3. Warmstarting with increasing dmax approximations. As such, a solution to (4.3) that is not optimal with respect to the original problem, may still be a good enough solution for an earthmoving problem. We have therefore an approximation algorithm that we call direct with x◦ and dmax . Algorithm 3: Direct with x◦ and dmax Input: x,y,c,A,b (from (2.1)) Output: Solution z ∗ = cT x∗ for (4.3) begin Build Rd0max using (4.2); (x◦ , y ◦ ) ←FindFeasiblePoint(x,y,c,A,b,Rd0max ); Solve (4.3) for (x0 , y 0 ) with starting point (x◦ , y ◦ ); end Note that in this algorithm, we modified FindFeasiblePoint in order to take Rd0max into account. The change is implicit and we do not show the adapted procedure here. 4.3 Warmstarting with increasing dmax In Section 4.2, we mention that the dmax -heuristic does not guarantee an optimal solution for (4.1). However, the solution from (4.3) can be used as a starting point to solve (4.1) using a warmstart procedure. Further improvement can be made by starting to solve (4.3) with a small dmax and use the solution to warmstart (4.3) with αdmax where α ∈ R++ , α > 1. If we continue this procedure, we end up with the following algorithm, which we refer to as warmstart algorithm. 39 4.4. Warmstarting with mipgap check Algorithm 4: Warmstart Input: x,y,c,A,b,d0max ≥ 1,α > 1 Output: Solution z ∗ = cT x∗ for (4.1) begin k ← 0; Build Rd0max using (4.2); (x◦ , y ◦ ) ←FindFeasiblePoint(x,y,c,A,b,Rd0max ); Solve (4.3) for (x0 , y 0 ) with starting point (x◦ , y ◦ ); while dkmax < ms do k ← k + 1; k−1 ; dkmax ← αdmax Build Rdkmax using (4.2); Solve (4.3) for (xk , y k ) with starting point (xk−1 , y k−1 ); end (x∗ , y ∗ ) ← (xk , y k ) end Lemma 4.6. The warmstart method in Algorithm 4 converges to the optimal solution value z ∗ of the MILP (4.1). Proof. Recall that ms is the number of sections. Let 1 ≤ d0max < ms and α > 1. Hence there exists k ∈ Z+ , such that αk d0max ≥ ms . Therefore, at iteration k of Algorithm 4, we have k−1 dkmax = αdmax = αα · · · αd0max = αk d0max ≥ ms . For this k, the extended MILP (4.3) is equal to the original MILP (4.1). 4.4 Warmstarting with mipgap check ∗ | is the gap between the current best Definition 4.7. The mipgap |z k − zLP feasible value of the MILP, z k , and the optimal value of its LP relaxation, ∗ . zLP ∗ of the LP relaxation for to the original Lemma 4.8. The optimal value zLP MILP (4.1) gives a lower bound for the dmax -extended MILP (4.3). Proof. Combine Fact 4.2 and Lemma 4.5. 40 4.4. Warmstarting with mipgap check ∗ be the optimal value for the LP relaxation of (4.1) Lemma 4.9. Let zLP and z ∗ be the optimal value for the MILP (4.1). For any current best feasible value zdkmax of (4.3), ∗ ∗ |z ∗ − zLP | ≤ |zdkmax − zLP |. Proof. The result follows from Lemma 4.5 and Lemma 4.8. We can now extend Algorithm 4 with a stopping condition that verifies the mipgap. We call this the warmstart with ǫ algorithm. Algorithm 5: Warmstart with ǫ Input: x,y,c,A,b,d0max ≥ 1,α > 1,ǫ > 0 Output: Solution z ∗ = cT x∗ for (4.1), within ǫ begin ∗ of Solve for the LP relaxation solution zLP (4.1); k ← 0; Build Rd0max using (4.2); (x◦ , y ◦ ) ←FindFeasiblePoint(x,y,c,A,b,Rd0max ); Solve (4.3) for z 0 = cT x0 using starting point (x◦ , y ◦ ); ∗ | do while dkmax < ms and ǫ < |z k − zLP k ← k + 1; k−1 ; dkmax ← αdmax Build Rdkmax using (4.2); Solve (4.3) for (xk , y k ) with starting point (xk−1 , y k−1 ); end (x∗ , y ∗ ) ← (xk , y k ) end In the next chapter, we look at the numerical results for Algorithms 1 - 5 with a set of problems based on real road alignment designs. We analyze solving times and quality of solutions for the different algorithms and try to find good values for the parameters α and dmax . 41 Chapter 5 Numerical results Understanding is, after all, what science is all about-and science is a great deal more than mere mindless computation. Roger Penrose (1931–) There are two categories of questions that we seek to statistically examine in this section. First, which algorithm is the fastest, and which algorithm is the most robust to solve our MILP model? Second, we want to know how the solution value of our model compares to the solution value of a LP solution for a given non-optimal removal schedule. 5.1 5.1.1 Experimental setup Problems In order to compare the algorithms in Chapter 4, we based our test problems on a family of real road alignments. We used a road design software to generate 8 different road alignments with length between 1 and 40 kilometers. We denote the alignments with letters from A to H. From the 8 designs, we then created 44 problems with different block and access road configurations. For all of the problems we added two access roads, one at the beginning and one at the end of the road. For about half of all the problems we added two additional access roads in between the start and the end of the road. We further split each group into three subgroups of 5, 10, and 20 blocks. The problem configurations are given in Table A.1 of Appendix A. 5.1.2 Benchmark In order to evaluate the benefits of our model over current practices, we would need to conduct a field study on a sufficiently large set of real construction projects. Similar to our theoretical problem set, such projects 42 5.1. Experimental setup would need to have obstacles and access roads in different configurations. The effective earthwork costs could then be compared with the solution of our model. In our research, we did not have access to such data. We therefore tried to build a reference solution that, to our knowledge, would be as fair as possible to compare with our solution. Reference solution For our reference solution we assume that obstacles are removed in a sequential order. We assume that earthwork tasks can occur at the same time as the obstacle removal process, as long as the sections for these tasks are accessible. We therefore use the removal scheme in Procedure FindFeasiblePoint as reference block removal schedule. As in Procedure FindFeasiblePoint, we solve the LP for the fixed schedule. Note that we do not use any dmax restrictions in the reference solution. We call this reference solution the LP solution for mass removal schedule. We are aware that block removal schedules vary for different construction projects and they also depend heavily on the experience of the geotechnical engineer. A mass removal schedule, however, is perfectly valid in forest areas or in ground conditions where the operation of ordinary construction equipment is not possible. In very difficult ground settings, tractor mounted equipment is used to push large-sized granular material in front of the equipment to build a stable construction platform [Fwa06]. This corresponds to a sequential removal of ground obstacles. Parameters For some of the tests we changed parameters for the algorithms in Chapter 4, such as: − ε is the tolerance for the MILP solver; − εLP is the LP optimality tolerance to solve the LP relaxation; − ǫ is the mipgap defined in Definition 4.7; − dmax ] is the maximal allowed moving distance; − α is the scaling factor for dmax within the warmstarting scheme. 43 5.1. Experimental setup Note that except for the warmstart algorithms, ε = ǫ. For the warmstart algorithms, we use 0.05 if dkmax < ms , ε= ǫ if dkmax ≥ ms . ∗ (see Algorithm In some problems, the solving time for the LP relaxation zLP 5) can be quite long. To reduce this solving time, we change the LP solver optimality tolerance εLP . The LP solver in the Gurobi Optimizer has a default LP optimality tolerance of εLP = 10−6 . An increase in the LP optimality tolerance results in a faster solving time for an LP. Now if for the mipgap tolerance ǫ > 0, we have ∗ |z ∗ − zLP | ≤ ǫ, then for some δ > 0, by the triangle inequality, we have ∗ |z ∗ − zLP + δ| ≤ ǫ + δ. So if we let the new LP optimality tolerance ε˜LP > εLP , with δ = ε˜LP − εLP and the new mipgap tolerance ǫ˜ = ǫ+δ, then the convergence of Algorithm 5 within a mipgap tolerance of ǫ is still guaranteed. For our experiments with the warmstart and ǫ (d0max = 25, α = 10) and warmstart with ǫ (d0max = 5, α = 5) algorithm, we used ε˜LP = 10−5 to solve the LP relaxation and ǫ˜ = ǫ + δ for the mipgap tolerance. For all the other computations, we used the default LP optimality tolerance. Test configurations In our first set of experiments, our goal is to investigate the behavior of the solution and solving time with respect to dmax . Hence we run all the problems with the direct with x◦ and dmax algorithm (Algorithm 3). The configuration for these tests are given in Table 5.1. In our second set of experiments, our goal is to investigate the economical benefit and the performance of the algorithms in Chapter 4. The LP solution for mass removal schedule was used for an economical comparison to the exact solution to (2.1), as given by the algorithms in Chapter 4. Furthermore, a performance analysis was done. The LP with mass removal schedule was not used in the performance comparison since it is a pure LP, and it is of little value to compare the solving time of a single LP with the algorithms in Chapter 4. Moreover, the LP with mass removal schedule is not an optimal solution to (2.1). We therefore focus our performance analysis on the MILP 44 5.1. Experimental setup Table 5.1: Test configurations for dmax with the direct solve algorithm using x◦ (Algorithm 3) ǫ dmax ǫ dmax 0.01 5 25 100 200 0.05 5 25 100 200 Table 5.2: Economical test and performance test configurations ǫ 0.01 0.05 d0max α - - 25 10 warmstart (Algorithm 4) warmstart with mipgap (Algorithm 5) 5 5 warmstart with mipgap (Algorithm 5) - - direct solve (Algorithm 1) direct with x◦ (Algorithm 2) 25 10 warmstart (Algorithm 4) warmstart with mipgap (Algorithm 5) 5 5 warmstart with mipgap (Algorithm 5) Algorithm direct solve (Algorithm 1) direct with x◦ (Algorithm 2) 45 5.2. Results algorithms only. In the context of economical analysis and performance, we run the experiments with the configurations shown in Table 5.2. For all our tests, we used a CPU time limit of 12 hours. If an algorithm could not solve the problem within the time limit, we aborted the solving process and considered the problem as not solved. Remark 5.1. In practice, a MILP tolerance of 5% is sufficiently precise for an earth-moving problem. The data for earth quantities is based on approximations that mostly use triangular surface grids. The use of triangularization, measurement errors, and non-constant factors for earth shrinkage and swelling, introduce errors that are in the order of 5%. Test equipment and software All the experiments were solved with the Gurobi Optimizer 2.0.2 MILP solver (http://www.gurobi.com) on a standard 64-bit Linux (Ubuntu 9.10 Server) workstation with an Intel R 2.4GHz Quad-Core CPU and 8GB of memory. Since the solver makes use of multiple cores, the reported CPUtime may be greater than the wall-clock time in some problems. The algorithms were implemented in C++ and the code can be made available upon request. Please contact the author at valentin.koch@gmail.com. 5.2 5.2.1 Results Analysis of the dmax -heuristic In the top row in Figure 5.1 we show the boxplots of the relative error of the objective value with respect to the solution for dmax ≥ ms and the corresponding MILP solver tolerances of ε=1% and ε=5%. The bottom row contains the relative solving time with respect to the solving time for dmax ≥ ms . In the boxplots, the box represents the interquartile range (IQR) and contains 50% of the sample points. It gives therefore a relatively robust statistic. The vertical bar inside the box represents the median. The ends of the whiskers represent the lowest datum still within 1.5 IQR of the lower quartile, and the highest datum still within 1.5 IQR of the upper quartile. Outliers are depicted as small dots. For convenience, the mean is also plotted as a small cross. We only counted problems that could be solved for both, varepsilon=1% and ε=5%. The numbers for the plots are given in Table A.2 and A.3 in the Appendix. Looking at dmax = 5 in the 1% configuration, we see that the median relative error of 34.84% for the objective value is about 5% smaller than the 46 5.2. Results 5% mipgap 1% mipgap 1.4 Mean 1.2 Relative Error 1.4 1 1 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0 0 -0.2 -0.2 5 25 100 2 Relative Time Mean 1.2 5 ms 25 100 2 Mean 1.5 1.5 1 1 0.5 0.5 0 ms Mean 0 5 25 100 dmax ms 5 25 100 ms dmax Figure 5.1: Objective and solving time with increasing dmax and different tolerance 47 5.2. Results median relative error of 40.39% for the 5% configuration. Remember that the reference solution for the relative error is also solved within a tolerance of 1% or 5%, respectively. With regard to Remark 5.1, we conclude that, when using dmax = 5, the relative error is too large for a solution to be acceptable. At dmax = 25, however, the median error is 0% in both cases, with an IQR of 4.6% and 3.3%, respectively. At the same time, we can look at the relative solving time and we can see that for 1% tolerance, at dmax = 5, we need 10.2% of the time that the direct solve with x◦ algorithm uses to solve the MILP. The time for a 5% tolerance at dmax = 5 is slightly smaller, at 6.4% of the direct solve and x◦ , which can be explained with the earlier convergence. For dmax = 25 we ended up with 26.4% and 13.8%, respectively. A quick glance at the figure tells us that there is not much difference at dmax = 100 to the direct solve algorithm. Consider, however, that 80% of all problems have ms ≤ 100 and therefore have the same results with dmax = 100 and dmax ≥ ms . If we add the IQR to the solving tolerance, we see that a choice of dmax = 25 produces a reasonable result that lies usually within 5.6% of the real solution with 1% solving tolerance, and within 4.3% of the solution with 5% solving tolerance. The results with 1% solving tolerance do not lie within the 5% suggested by Remark 5.1. However, a solution from dmax = 25 is sufficiently close to the real solution in order to give a good upper bound for a subsequent MILP with a larger dmax or dmax ≥ ms . The median solving time with dmax = 5 for 1% tolerance is around 10 times faster, and with 1% tolerance around 16 times faster than the direct with x◦ algorithm. In more than half of all the problems, the solving time does not exceed 16% of the solving time with the direct and x◦ algorithm for both tolerances. In the dmax = 25 case, the solving time is about 4 and 7 times faster. We observe that the upper quartile for 1% is at about half and for 5% at about a third of the solving time with the direct with x◦ algorithm. We will use these observations in the analysis of the warmstart scheme in order to derive the effect of the starting point given by dmax = 5 and dmax = 25. 5.2.2 Performance In order to compare the performance of the algorithms, we use a performance profile as developed in [DM02]. For each algorithm, we plot the number of problems that are solved within a factor of the best CPU time. Let P be the set of all problems, and A be the set of algorithms. We compute the 48 5.2. Results performance ratio by rp,a = tp,a , min{tp,a : a ∈ A} where p ∈ P, and tp,a is the solving time of algorithm a for problem p. The percentage of problems that are solved within a factor τ ∈ R are computed by |{p ∈ P : rp,a ≤ τ }| ρa (τ ) = . |P| Hence for a particular problem p, an algorithm a needs τ times longer than the best algorithm for p, if τ > 1 when algorithm a is able to solve it. If algorithm a needs only τ = 1, the algorithm is considered as the best, or among the best algorithms for problem p. The percentage ρa of all problems solved by algorithm a at τ gives us an overall assessment of the performance of algorithm a. We conclude that the algorithm with the highest ρ value is likely to perform best on a given problem with a probability of ρa (τ ). For a good algorithm, the plot of ρ(τ ) should therefore be found above the other algorithms. If we increase τ to allow the algorithms to be within a larger factor of the best algorithm, then the algorithm with the highest ρ is the most robust algorithm. Hence an algorithm with the highest ρ value at a large τ , solves most of the problems. For a more detailed description of performance profiles, we refer the reader to [DM02]. Figure 5.2 shows the performance profiles for Algorithms 1 to 5 on a log2 scale, that is, |{p ∈ P : log2 (rp,a ) ≤ τ }| . ρa (τ ) = |P| Figure 5.2a shows the profile for ε=1% and Figure 5.2b shows the profiles with ε=5%. In Figure 5.3a and 5.3b, the profiles are shown on a normal scale with a short τ interval. Figures 5.4a and 5.4b show the profiles on a normal scale with a long τ interval. For all profiles, we have |A| = 5 and |P| = 44. For a ε=1%, we can see that the direct solve with x◦ algorithm outperforms the other algorithms with a probability of about 32%, and the probability that the warmstart with ǫ (d0max = 25, α = 10) algorithm outperforms all other algorithms is at about 27%. If we increase τ to let the algorithms be within a factor of 4 of the direct solve with x◦ algorithm, we can see that the warmstart with ǫ (d0max = 25, α = 10) algorithm wins with a high probability of about 78% and solves up to 80% of the problems when τ is extended to 5. We can see that the warmstart with ǫ (d0max = 25, α = 10) is the most robust 49 5.2. Results 1 ρa (log2 (τ )) 0.8 0.6 0.4 direct solve direct with x◦ warmstart (d0max = 25, α = 10) warmstart with ǫ (d0max = 25, α = 10) warmstart with ǫ (d0max = 5, α = 5) 0.2 0 0 1 2 3 4 5 6 7 log2 (τ ) (a) 1% MILP tolerance 1 ρa (log2 (τ )) 0.8 0.6 0.4 direct solve direct with x◦ warmstart (d0max = 25, α = 10) warmstart with ǫ (d0max = 25, α = 10) warmstart with ǫ (d0max = 5, α = 5) 0.2 0 0 1 2 3 4 5 6 7 log2 (τ ) (b) 5% MILP tolerance Figure 5.2: Algorithms performance profiles with log2 scale 50 5.2. Results algorithm (82% of all problems solved), followed by the warmstart with not mipgap check and the direct solve with x◦ algorithm (79.5% of all problems solved). We conclude that the warmstart with ǫ (d0max = 25, α = 10) algorithm is a good choice to solve most of the problems relatively fast. Clearly, the direct solve algorithm is not very competitive in these experiments. The warmstart with ǫ (d0max = 25, α = 10) converged with a dmax ≤ 25 on about 34% of all solved problems. Recall that in Section 5.2.1, for ε=1%, we found that the median and lower quartile of the relative solution error for dmax = 25 is 0% with an upper quartile of 4.6%. Therefore half of the solutions were above 4.6% of the real solution. From Lemma 4.9 we expect that the warmstart with ǫ algorithm will use at least one more iteration for these solutions. Furthermore, due to the IQR of 4.6% for the solution error, there is also a significant number of solutions within the IQR that are above the required ε=1% tolerance. Altogether, these problems account for the 66% where the warmstart with ǫ algorithms were required to increment dmax from 25 to 125 or 250, depending on α. For ε=5%, the warmstart with ǫ (d0max = 25, α = 10) algorithm dominates all other algorithms on 46% of all problems. The algorithm clearly stays ahead of the others until τ ≈ 35. At that point, either the warmstart with ǫ (d0max = 25, α = 10) or the warmstart without mipgap check would suffice to solve about 97.7% of all the problems. The warmstart with ǫ (d0max = 25, α = 10) algorithm is therefore the clear winner in the ε=5% experiments. An explanation why the warmstart with ǫ (d0max = 25, α = 10) algorithm performs always better or equal to the other algorithms when ε=5% lies in the mipgap check. While the median and lower quartile for the solution error at dmax = 25 are, similar to the ε=1% experiments, at 0%, the upper quartile is only at 3.3%. Therefore half of all the problems that lie in the IQR at dmax = 25 are already within ε=5%. From the results, we know that the algorithm converged on 61% of all solved problems with a dmax ≤ 25. This is a notable result that confirms our dmax -analysis. The warmstart with ǫ (d0max = 5, α = 5) algorithm shows a mixed performance. It starts better than the simple warmstart with no mipgap check algorithm and if τ ≈ 3, it also outperforms the direct with x◦ algorithm with a probability of over 70%. However, after τ = 6, the warmstart with ǫ (d0max = 5, α = 5) algorithm performs only better than the direct solve algorithm and finished by solving 88.6% of all the problems. The direct solve algorithm solved only 81%. In the next section, we will analyze the quality of the solutions with respect to this LP with mass removal schedule. 51 5.2. Results 0.8 0.7 0.6 ρa (τ ) 0.5 0.4 0.3 0.2 direct solve direct with x◦ warmstart (d0max = 25, α = 10) warmstart with ǫ (d0max = 25, α = 10) warmstart with ǫ (d0max = 5, α = 5) 0.1 0 0 2 4 6 8 10 12 τ (a) 1% MILP tolerance 1 0.8 ρa (τ ) 0.6 0.4 direct solve direct with x◦ warmstart (d0max = 25, α = 10) warmstart with ǫ (d0max = 25, α = 10) warmstart with ǫ (d0max = 5, α = 5) 0.2 0 0 2 4 6 8 10 12 τ (b) 5% MILP tolerance Figure 5.3: Algorithms performance profiles with normal scale and small intervale for τ 52 5.2. Results 1 0.8 ρa (τ ) 0.6 0.4 direct solve direct with x◦ warmstart (d0max = 25, α = 10) warmstart with ǫ (d0max = 25, α = 10) warmstart with ǫ (d0max = 5, α = 5) 0.2 0 0 20 40 60 80 100 τ (a) 1% MILP tolerance 1 0.8 ρa (τ ) 0.6 0.4 direct solve direct with x◦ warmstart (d0max = 25, α = 10) warmstart with ǫ (d0max = 25, α = 10) warmstart with ǫ (d0max = 5, α = 5) 0.2 0 0 20 40 60 80 100 τ (b) 5% MILP tolerance Figure 5.4: Algorithms performance profiles with normal scale and large interval for τ 53 5.2. Results 5.2.3 Economical analysis In Section 5.1.2, we defined the LP with mass removal schedule as our reference solution. Our reference solution is based on a LP model and does not necessarily represent the real cost of an actual earthwork construction project with current techniques. ∗ be the solution for problem p using the LP with mass removal Let zp,ref ∗ be the solution for problem p using algorithm a ∈ A. We schedule. Let zp,a define the gain gp,a by ∗ zp,a gp,a = 1 − ∗ . zp,ref For all problems that have been solved, we plot the gain in a histogram. The number of bars in the histogram is the square root of the number of problems that have been solved. We also plot a normal curve over the histogram that fits the data. Note that we observed that in some instances where the problems were not solved, the current solution was still better than the solution of the LP with mass removal schedule. It is not clear if one should count these solutions in the histogram plots. We decided to keep these solutions in our economical analysis on the following grounds. In practice, if a solver is interrupted with an intermediate solution that is better than a traditional LP solution, the intermediate solution will most likely be used. In Section 5.2.2 we identified the best three performing algorithms as the warmstart with ǫ (d0max = 25, α = 10) algorithm (Algorithm 5), the direct with x◦ algorithm (Algorithm 2), and the warmstart (d0max = 25, α = 10) algorithm (Algorithm 4). Since the warmstart with ǫ (d0max = 25, α = 10) algorithm solved more problems than the other two algorithms, we focus our economical analysis on the results of the warmstart with ǫ (d0max = 25, α = 10) algorithm. Figure 5.5 shows the histogram for the warmstart with ǫ (d0max = 25, α = 10) algorithm. Histograms for the other two algorithms are shown in Appendix B. The statistics for the histograms can be found in Appendix A, Tables A.4 and A.5. For the problems that have been solved, the use of the warmstart with ǫ (d0max = 25, α = 10) algorithm and ε=1% produced an average gain of 4.1%. The same algorithm delivered an average of 3.3% gain if ε=5%. From the plot we see that the negative gains in the ε=5% histogram, that are caused by the MILP tolerance, have a much greater influence on the average than in the ε=1% case, where almost no negative gains are encountered. The maximum gains are 17.9% in the ε=1% histogram, and 19% in the ε=5% histogram. 54 5.2. Results 5% mipgap Problems 1% mipgap 25 25 20 20 15 15 10 10 5 5 0 -0.2-0.15-0.1-0.05 0 0.05 0.1 0.15 0.2 0 -0.2-0.15-0.1-0.05 0 0.05 0.1 0.15 0.2 Figure 5.5: Savings with warmstart with ǫ (d0max = 25, α = 10), compared to mass block removal schedule Overall we can observe that the average gains for ε=1% are not significantly higher than for ε=5%. This is surprising since in general, the MILP solver stops improving the solution in the ε=5% case much earlier than in the ε=1% case. From the tables in the appendix we can see that the solving time for the ε=1% case is much higher than for ε=5%. The qualitative merit of our model can be assessed when looking at the best gains. For each of the two MILP tolerances, there were 8 problems with gains of 10% and more. This accounts for about 18% of all the problems. In the context of a large highway project, this is a significant amount of gain over the LP with mass schedule. We are interested in the problems that produced the biggest gains. In order to see if there are correlations between problem types and gains, we did a stratification analysis with the use of a scatter plot. The scatter plot is shown in Figure 5.6. The scatter plots for the other algorithms can be found in Appendix B. We first grouped the problems into the number of blocks. Each group is shown with a different color and symbol. Green stars show five block problems, blue circles show ten block problems, and red triangles show twenty block problems. Each group was internally sorted by the number of sections and plotted from left to right in ascending order. Hence, problem 1 in the scatter plot represents the five block problem with the least number of sections, whereas problem 17 represents a five block problems with the largest number of sections. Small points represent problems with two access roads, large points represent problems with four access roads. In the stratification analysis, we could not identify any visibly striking pattern. The ε=1% and the ε=5% results have a similar distribution. One might observe that in the groups of five block problems, there are less solu- 55 5.2. Results 1% mipgap 5% mipgap 0.2 0.2 0.15 0.15 0.1 0.1 0.05 0.05 0 0 -0.05 -0.05 0 10 20 30 40 50 0 10 20 30 40 50 Figure 5.6: Stratification of gains for warmstart with ǫ (d0max = 25, α = 10); * are 5, ⊙ are 10, are 20 blocks; small symbols are 2, big symbols are 4 access roads; within the block groups, points are sorted from left to right in ascending order of section numbers tions with significant gains than in the other two groups. For these problems, it seems, the LP with mass removal schedule provides often a good enough solution. If we group the results with the highest gains, we can observe a majority of problems with four access roads. Intuition tells us that with more access roads, there is more flexibility and hence a potential for cheaper solutions. We conclude that for a better stratification analysis, experiments with a larger set of problems would certainly give more insight. 56 Chapter 6 Conclusion Frustra fit per plura quod potest fieri per pauciora.4 William of Ockham (1287-1347) 6.1 Accomplishments We started the thesis with the description of traditional methods to balance earthwork and minimize hauling costs in the context of road construction. Two main methods were presented in Chapter 1. First, we explained the mass diagram, a simple tool that helps to balance the earthwork and to obtain an estimate on the directions of haul. Second, we showed a more detailed model of a LP by Mayer and Stark [MS81] that not only balances the earthwork, but also minimizes the hauling cost. We then pointed out that neither method captures the problem of blocks that restrict some of the earth movements to balance earth across the construction site. To address this problem, we presented an extension to the model of [MS81] in Chapter 2. In Chapter 3, we analyzed the model with respect to space complexity and we investigated the growth of the number of feasible block removal schedules. To improve the computation time for solving the MILP, we outlined a set of algorithms in Chapter 4. In Chapter 5, we presented the numerical results of the algorithms, together with an analysis of performance, robustness, and accuracy. 6.2 Research outcomes In the model analysis, we first identified a potential issue of memory growth. With a fixed number of blocks, we found that the theoretical upper bound on the growth of memory for the MILP model in (2.1) is O |M|4 . In 4 It is pointless to do with more what can be done with less. 57 6.2. Research outcomes our experiments, we use problems with |M| 200. However, models with |M| ≈ 200 required relatively large amounts of the available memory. The solver failed to solve some of these problems due to memory exhaustion. A second result of our model analysis was the enumeration of the feasible block removal schedules. We concluded that for a road with one access road at either end, there is no need to solve for the MILP, and one can just solve a LP for a fixed removal schedule that removes the blocks sequentially. If the access road is at some other place, we can quickly compute the number of feasible schedules with equation (3.4). If the number of feasible schedules is small, we can enumerate the schedules and solve a LP for every schedule. The schedule that produces the smallest objective value is the optimal removal schedule. We also showed that the worst case of number of schedules for one access road is bounded by O 3|I| . Therefore, if |I| is large, we recommend to solve the MILP instead of a schedule enumeration. The same observations were made for the typical configuration of two access roads, one on each end. In our numerical analysis, we showed that the MILP can be solved efficiently with the help of the algorithms from Chapter 4. In the experiments, we demonstrated the strength of the dmax -heuristic, and showed that there is a significant performance benefit when using the warmstart with ǫ, the simple warmstart, or the direct with x◦ algorithm, instead of solving the MILP directly. Furthermore, we showed that these three algorithms are considerably more robust than a method that solves directly for the MILP. As a last observation in our analysis, we showed the economic benefits of solving the MILP instead of using a LP with a fixed mass removal schedule. If we solve a problem to an accuracy of ε=1%, the MILP produces solutions that are about 4.1% cheaper in the average. If we aim for ε=5% accuracy, the MILP provides us with a solution that saves roughly 3.3% on a LP solution for a fixed mass removal schedule. At a first sight, these benefits may not seem very high. However, in the context of a highway project that involves earthwork costs of millions of dollars, 4.1% is a significant figure. In summary, we started by presenting a MILP model that overcomes the problems of the mass diagram and the LP model by Mayer and Stark. We pointed out potential issues with the model and efficient ways to solve the MILP. The result of the model presents an engineer with a good estimate on earthwork costs. The model is therefore a good tool to quantify a given vertical alignment, and therefore provides a valuable tool to find the most economic alignment. As already pointed out by Mayer and Stark [MS81], the objective of an earthwork contractor is to maximize its profit. A contract that is based on 58 6.3. Future research an engineer’s mass diagram, and that fixes the costs for excavation, borrow, and haul, may not reflect additional costs that result from a haul scheme that is not possible due to blocks. The presented MILP enables a contractor to estimate the real earthwork costs of a project, and at the same time provides the contractor with the cheapest haul scheme in order to maximize its profit. 6.3 6.3.1 Future research Warmstart schemes and εLP In Section 5.2.2, we explained the problems that we encountered when solving the LP relaxation with the default LP optimization tolerance of the Gurobi Optimizer. At the writing of this thesis, preliminary results of tests with εLP = 10−5 suggest an even better performance for the warmstart with ǫ (d0max = 25, α = 10) and warmstart with ǫ (d0max = 5, α = 5) algorithms. We also would like to investigate additional warmstart schemes with different d0max - and α-values. 6.3.2 Vertical alignment The above model can easily be extended to the vertical alignment. Based on a vertical optimization model by Easa [Eas87], Moreb introduced a LP model that combines the earthwork model by Mayer and Stark with roadway grade selection by using a piecewise linear approximation of the road profile [Mor96]. In [MA04], Moreb and Aljohani extended the model by using a quadratic approximation. The quadratic representation was then further extended in [Mor09] with a spline approximation, with improvements by Koch and Lucet [KL10]. The result was an improved and efficient linear model that avoids sharp connectivities and scaling problems of piecewise linear representations. The main difference between the LP model of Mayer and Stark, and the vertical alignment model by Moreb, is the earth volume. In the earthwork LP and in our MILP model, we treat the volume of a section i as a fixed constant Vi . We recall the balance constraint (1.3b) from the LP as (xij − xji ) = Vi j∈M j=i for all i ∈ S. However, in the vertical alignment LP by Moreb, the volume is a linear approximation. The volume is computed as the product of the surface area 59 6.3. Future research Ai of section i, with the average height of earth that is cut from the section or filled into that section. The cut or fill is represented with LP decision variables, namely ui for the height of the cut, and vi for the height of the fill in that section. The balance constraint then becomes j∈M j=i (xij − xji ) = Ai (ui − vi ) for all i ∈ S. The ground profile can then be modeled with spline segments. The spline segments do not need any additional decision variables besides ui and vi . We refer the reader to [Mor09, KL10] for the details. Since the model for vertical alignment needs only two additional continuous decision variables per section, the MILP model in (2.1) can readily be combined with the LP model in [Mor09, KL10]. 6.3.3 Schedule enumeration In Remark 3.25, we discussed the use of feasible schedule enumeration in the case of a small number of blocks. The use of a Depth-first search algorithm allows an enumeration in linear time. However, the number of feasible schedules itself is exponential. A natural question is what a small number of blocks means in this context. It would be interesting to compare the solving time for the enumeration of schedules, together with solving the LP for each schedule, and the solving time for the MILP. We believe that up to a certain number of blocks, the enumeration of all LP solutions would be significantly faster than solving the MILP. The decision process between the two methods would be simplified significantly with the identification of a threshold that indicates the number of blocks where the MILP will outperform an exhaustive enumeration. An additional topic of interest, with regard to the schedule enumeration, is the use of dynamic programming. In Section 3.2, we described a removal schedule in terms of binary strings, one binary string for each time step. An interesting observation is that the same binary string can be used in multiple schedules. In fact, after a certain time step, multiple schedules may even end with the same set of binary strings. We can therefore look at a binary string as a state in the removal process5 . Once two schedules share the same state, they need to solve the same subproblem to finish the removal process. Dynamic programming is particularly well suited for optimization problems with overlapping subproblems [Bel52]. It would be interesting to formulate 5 This observation was made by Jim Nastos. 60 6.3. Future research the feasible schedule enumeration as a dynamic program, and to compare the method with exhaustive enumeration, as well as the MILP method. If a threshold was found for the exhaustive enumeration and the MILP, it would be interesting to see if the threshold changes with dynamic programming. 6.3.4 Lagrangian relaxation The final topic for future research is the use of Lagrangian relaxation [Fis04, Lem07]. Let x ∈ Rn . Consider the MILP z = min cT x s.t. Ax = b Dx ≤ e x ≥ 0 and integral, where b ∈ Rm , e ∈ Rk , A ∈ Rm×n , and D ∈ Rk×n . The Lagrangian problem of the above MILP is zD (λ) = min cT x + λ(Ax − b) s.t. Dx ≤ e x ≥ 0 and integral, where λ = (λ1 , . . . , λm ) is a vector of Lagrange multipliers. In [Fis04], it is shown how to choose appropriately λ in order to approach z. From data of preliminary experiments, we believe that the MILP model could be solved much faster if some of the constraints are relaxed by moving them into the objective function. For example, we could readily relax the balance constraints (2.1b),(2.1c),(2.1d), and the stockpile constraints (2.1i),(2.1j), without affecting the binary variables. The relaxation of the feasible move constraints (2.1g) and (2.1h) would greatly reduce the problem size. 61 Bibliography [Aig07] Martin Aigner. A Course in Enumeration, volume 238 of Graduate Texts in Mathematics. Springer-Verlag, 2007. → pages 29, 30 [Bel52] Richard E Bellman. The theory of dynamic programming. Proceedings of the National Academy of Sciences of the United States of America, 38(8):716–719, 1952. → pages 60 [BS05] Cyril Banderier and Sylviane Schwer. Why delannoy numbers? Journal of Statistical Planning and Inference, 135(1):40–54, 2005. → pages 33 [CGF89] E P Chew, C J Goh, and T F Fwa. Simultaneous optimization of horizontal and vertical alignments for highways. Transportation Research Part B: Methodological, 23(5):315–329, 1989. → pages 2, 4 [CL06] Juey-Fu Cheng and Yusin Lee. Model for three-dimensional highway alignment. Journal of Transportation Engineering, 132(12):913–920, 2006. → pages 4 [CLRS09] Thomas H Cormen, Charles E Leiserson, Ronald L Rivest, and Clifford Stein. Introduction to Algorithms. The MIT Press, 3 edition, Jan 2009. → pages 22, 34 [Dan98] George B Dantzig. Linear Programming and Extensions. Princeton Landmarks in Mathematics. Princeton University Press, 1998. → pages 36 [DM02] Elizabeth D Dolan and Jorge J Mor´e. Benchmarking optimization software with performance profiles. Mathematical Programming, 91(2):201–213, 2002. → pages 48, 49 [Eas87] Said M Easa. Earthwork allocations with nonconstant unit costs. Journal of Construction Engineering and Management, 113(1):34–50, March 1987. → pages 7, 9, 59 62 Chapter 6. Bibliography [Eas88a] Said M Easa. Earthwork allocations with linear unit costs. Journal of Construction Engineering and Management, 114(4):641– 655, December 1988. → pages 7, 9 [Eas88b] Said M Easa. Selection of roadway grades that minimize earthwork cost using linear programming. Transportation Research. Part A: general, 22(2):121–136, 1988. → pages 4 [Fis04] M Fisher. The Lagrangian relaxation method for solving integer programming problems. Management science, 50(12):1861–1871, December 2004. → pages 61 [Fou05] R Fourer. Linear programming - software survey. Operations Research / Management Science Today, (3):46, 2005. → pages 35, 36 [Fwa06] T F. Fwa. The Handbook of Highway Engineering. CRC Press, 1 edition, September 2006. → pages 43 [GH08] Nicholas J. Garber and Lester A. Hoel. Traffic and Highway Engineering. CL-Engineering, 4 edition, June 2008. → pages 4 [Gom58] Ralph E Gomory. Outline of an algorithm for integer solutions to linear programs. Bulletin of the American Mathematical Society, 64:275–278, 1958. → pages 36 [JH90] A Jayawardane and F Harris. Further development of integer programming in earthwork optimization. Journal of Construction Engineering and Management, 116(1):18–34, Mar 1990. → pages 7 [JK06] Manoj Kumar Jha and Euncheol Kim. Highway route optimization based on accessibility, proximity, and land-use changes. Journal of Transportation Engineering, 132(5):435–439, May 2006. → pages 3 [JS03] J Jong and P Schonfeld. An evolutionary model for simultaneously optimizing three-dimensional highway alignments. Transportation Research Part B: Methodological, 37(2):107–128, February 2003. → pages 4 [JSJ06] Manoj Kumar Jha, Paul Schonfeld, and J-C Jong. Intelligent Road Design, volume 19. WIT Press, 2006. → pages 2, 3 63 Chapter 6. Bibliography [KL10] Valentin R Koch and Yves Lucet. A note on: Spline technique for modeling roadway profile to minimize earthwork cost. Journal of Industrial and Management Optimization, 6(2):393–400, May 2010. → pages 4, 59, 60 [Lay99] M G Lay. Ways of the World: A History of the World’s Roads and of the Vehicles that Used Them. Rutgers University Press, 1999. → pages 1 [Lem07] Claude Lemar´echal. The omnipresence of Lagrange. Annals of Operations Research, 153(1):9–27, September 2007. → pages 61 [LTL09] Yusin Lee, You-Ren Tsou, and Hsiao-Liang Liu. Optimization method for highway horizontal alignment design. Journal of Transportation Engineering, 135(4):217–224, Jan 2009. → pages 3 [MA04] Ahmad A Moreb and M Aljohani. Quadratic representation for roadway profile that minimizes earthwork cost. Journal of Systems Science and Systems Engineering, 13(2):245–252, April 2004. → pages 59 [Mak99] K Makanae. An application of parametric curves to highway alignment. Journal of Civil Engineering Information Processing System, 8:231–238, Jan 1999. → pages 4 [McA24] John Loudon McAdam. Remarks on the Present System of Road Making. London, Longman, Hurst, Rees, Orme, Brown, and Green, 1824. → pages 1 [MKV02] S M¨aki, R Kalliola, and K Vuorinen. Road construction in the Peruvian Amazon: process, causes and consequences. Environmental Conservation, 28(03):199–214, 2002. → pages 2 [Mor96] Ahmad A Moreb. Linear programming model for finding optimal roadway grades that minimize earthwork cost. European Journal of Operational Research, 93(1):148–154, August 1996. → pages 4, 59 [Mor09] Ahmad A Moreb. Spline technique for modeling roadway profile to minimize earthwork cost. Journal of Inustrual and Management Optimization, 5(2):275–283, 2009. → pages 4, 59, 60 64 [MS81] R Mayer and R Stark. Earthmoving logistics. Journal of the Construction Division, 107(2):297–312, Jun 1981. → pages ii, 6, 7, 11, 13, 57, 58 [PS98] Christos H Papadimitriou and Kenneth Steiglitz. Combinatorial Optimization: Algorithms and Complexity. Dover Publications, 1998. → pages 36 [SAV+ 03] Robert Sulanke, Jose Arregui, Danijela Vuji´c, Nelson Casta˜ neda, and David Callan. Another description of the central delannoy numbers: 10894. The American Mathematical Monthly, 110(5):443–444, May 2003. → pages 31 [Stu09] Les Stuart. Transportation in Canada 2008 - an overview. Transport Canada, pages 1–28, Jun 2009. → pages 1, 2 [Trb90] Trb. State of the art report 8: Guide to earthwork construction. Transportation Research Board, pages 1–119, Feb 1990. → pages 13 [Tri87] D Trietsch. A family of methods for preliminary highway alignment. Transportation Science, 21(1):17–25, February 1987. → pages 3 [Wei02] Eric W. Weisstein. CRC Concise Encyclopedia of Mathematics. CRC Press, December 2002. → pages 31, 33 65 Appendix A Tables Each design in Table A.1 is used with 2 access roads, and with 4 access roads for the same block configuration. An exception are the designs C and G with 20 blocks, where the use of 4 access roads was too hard to solve. Problems 23 and 24 are configurations were the blocks are mostly clustered together in adjacent sections. The data for the problems is available upon request. Please contact the author at valentin.koch@gmail.com. Table A.1: Problem configurations Problem Design 1,25 2,26 3,27 4,28 5,29 6,30 7,31 8,32 9,33 10,34 11,35 12,36 13,37 14,38 15,39 16,40 17,41 18 19,42 20,43 21,44 22 23 24 A B C D E F G H A B C D E F G H A C D E F G H F Length (km) Sections Access roads Blocks 4 40 1 1 1 1 9 9 4 40 1 1 1 1 9 9 4 1 1 1 1 9 9 1 101 201 51 100 52 103 100 201 101 201 51 100 52 103 100 201 101 51 100 52 103 100 201 103 2,4 2,4 2,4 2,4 2,4 2,4 2,4 2,4 2,4 2,4 2,4 2,4 2,4 2,4 2,4 2,4 2,4 2 2,4 2,4 2,4 2 2 2 5 5 5 5 5 5 5 5 10 10 10 10 10 10 10 10 20 20 20 20 20 20 5 10 66 Appendix A. Tables Table A.2: Statistics for objective function error with dmax ǫ dmax = 5 dmax = 25 dmax = 100 dmax = ms 0.01 Minimum Lower quartile Median Upper quartile Maximum Lower conf. for med. Upper conf. for med. 0.03344 0.21869 0.34841 0.60876 1.08125 0.23842 0.45841 -0.00907 0.00000 0.00000 0.04625 0.18926 0.00000 0.01304 -0.00351 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.05 Minimum Lower quartile Median Upper quartile Maximum Lower conf. for med. Upper conf. for med. -0.01434 0.22466 0.40394 0.60948 1.08494 0.30959 0.49830 -0.02768 0.00000 0.00000 0.03376 0.20007 0.00000 0.00828 -0.04701 0.00000 0.00000 0.00000 0.01690 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 Table A.3: Statistics for solving time with dmax ǫ dmax = 5 dmax = 25 dmax = 100 dmax = ms 0.01 Minimum Lower quartile Median Upper quartile Maximum Lower conf. for med. Upper conf. for med. 0.00000 0.03820 0.10172 0.16049 0.44444 0.06544 0.13801 0.00016 0.11630 0.26389 0.54474 1.42857 0.13677 0.39101 0.00051 0.59407 0.99488 1.00000 1.50000 0.87444 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 0.05 Minimum Lower quartile Median Upper quartile Maximum Lower conf. for med. Upper conf. for med. 0.00000 0.03800 0.06395 0.15395 0.33205 0.03800 0.09310 0.00968 0.09709 0.13793 0.28175 1.27652 0.09709 0.18435 0.21221 0.86429 1.00000 1.00000 1.74622 0.96588 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 67 Appendix A. Tables Table A.4: Histogram statistics for 1% tolerance Solved Min 1st qrt. Median 3rd qrt. Max Mean Std. dev Skew direct direct x◦ warmstart warmstart ǫ, 25 warmstart ǫ, 5 30 -0.00674 0.00041 0.02386 0.05308 0.17902 0.04341 0.05429 1.13166 35 0.00000 0.00000 0.01878 0.04827 0.17902 0.03950 0.05132 1.29910 35 -0.00674 0.00000 0.01878 0.05051 0.18987 0.04122 0.05519 1.35167 36 -0.00198 0.00000 0.02232 0.05079 0.17902 0.04088 0.05154 1.24205 33 -0.00674 0.00000 0.01696 0.04658 0.15804 0.03265 0.04227 1.34292 Table A.5: Histogram statistics for 5% tolerance Solved Min 1st qrt. Median 3rd qrt. Max Mean Std. dev Skew direct direct x◦ warmstart warmstart ǫ, 25 warmstart ǫ, 5 36 -0.01722 0.00000 0.01537 0.04733 0.17873 0.03865 0.05471 1.18703 42 0.00000 0.00000 0.00000 0.04055 0.17815 0.03176 0.05305 1.46810 43 -0.02907 0.00000 0.00000 0.04290 0.18987 0.03309 0.05527 1.35854 43 -0.02907 0.00000 0.00000 0.04290 0.18987 0.03309 0.05526 1.35773 39 -0.04251 -0.00115 0.00615 0.03903 0.15804 0.02331 0.04658 1.17794 68 Table A.6: CPU time in seconds for ε = 0.01, problems 1-22 direct success 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 79 911 21 26 14 55 14 316 7104 43200 603 2936 606 107860 376 15777 43200 43200 43324 43200 28 490 1 1 1 1 1 1 1 1 1 0 1 1 1 0 1 0 0 0 0 0 1 1 direct x◦ 77 1199 16 18 15 21 14 316 1603 22338 412 7522 134 43256 146 4864 84261 43200 43335 43200 22 692 success 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 0 0 0 0 1 1 warmstart 82 986 30 23 15 24 18 332 1507 43200 793 1528 100 468 206 5533 43200 43200 44453 43200 49 506 success 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 0 0 0 0 1 1 warmstart ǫ (25, 10) 73 2622 38 9 12 55 23 357 1380 17274 816 28593 16 540 285 5629 87130 43200 43224 43200 45 549 success 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 1 1 dkmax warmstart ǫ (5, 5) 250 250 250 25 25 250 250 250 250 250 250 250 25 250 250 250 250 25 25 25 250 250 68 3249 37 11 14 29 26 285 1171 18214 1472 3473 49 43231 223 11257 43337 43200 43176 43200 56 576 success 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 0 0 0 0 1 1 dkmax 125 625 125 25 25 125 125 125 125 625 125 125 25 125 125 625 5 125 5 125 125 625 Appendix A. Tables Problem 69 Table A.7: CPU time in seconds for ε = 0.01, problems 23-44 direct success 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 2 184 3 24 17 302 8531 43200 33 328 91 31175 197 43200 43200 2200 43200 43327 43200 56509 319 268 1 1 1 1 1 1 1 0 1 1 1 1 1 0 0 1 0 0 0 0 1 1 direct x◦ 2 17 3 25 16 298 15307 43200 20 188 31 22377 182 6486 43200 1425 9448 12576 44587 43208 270 173 success 1 1 1 1 1 1 1 0 1 1 1 1 1 1 0 1 1 1 0 0 1 1 warmstart 4 22 4 28 18 315 5369 43200 36 207 40 22458 196 7245 1536 2221 9628 43259 48538 43207 383 215 success 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 0 0 0 1 1 warmstart ǫ (25, 10) 5 9 2 9 8 337 10121 43200 15 39 35 8180 38 7362 7163 1430 11723 43202 49588 43200 353 221 success 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 0 0 0 1 1 dkmax warmstart ǫ (5, 5) 250 25 25 25 25 250 250 25 25 25 250 250 25 250 250 250 250 250 250 25 250 250 2 28 3 12 10 537 9506 43200 16 403 18 11565 54 3874 6891 43256 43342 43176 43200 43208 586 222 success 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 0 0 0 0 0 1 1 dkmax 25 125 25 25 25 625 125 625 25 125 25 125 25 125 125 5 25 5 125 125 625 125 Appendix A. Tables Problem 70 Table A.8: CPU time in seconds for ε = 0.05, problems 1-22 direct success 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 59 886 20 27 14 56 14 323 1232 43200 609 235 596 1862 383 43200 43200 43200 33900 43200 28 182 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 0 0 0 1 0 1 1 direct x◦ 19 601 14 18 15 19 14 322 182 2414 411 145 134 140 150 4972 5967 43200 43338 1187 22 181 success 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 1 1 1 warmstart 23 1013 29 23 13 23 18 338 150 3356 801 155 101 155 169 5635 1416 18899 43318 1295 27 189 success 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 warmstart ǫ (25, 10) 26 1044 30 9 13 8 8 367 163 1051 809 30 88 31 33 5755 1738 19198 43308 636 32 30 success 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 dkmax warmstart ǫ (5, 5) 250 250 250 25 25 25 25 250 250 25 250 25 25 25 25 250 250 25 250 25 250 25 40 107 38 15 15 12 14 386 242 885 1559 53 61 840 62 8530 43248 6564 42989 3106 39 69 success 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 0 1 1 1 dkmax 125 25 125 25 25 25 25 125 125 25 125 25 25 125 25 250 5 25 5 25 125 25 Appendix A. Tables Problem 71 Table A.9: CPU time in seconds for ε = 0.05, problems 23-44 direct success 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 2 154 3 25 17 309 186 5154 33 189 91 1536 200 14951 9991 1846 43200 8412 43200 57882 327 273 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 0 0 1 1 direct x◦ 2 18 3 26 16 305 169 5117 20 192 31 198 187 6628 1244 330 1130 1859 1732 6651 276 167 success 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 warmstart 3 21 4 28 19 323 204 5169 27 212 37 214 200 7412 1418 307 1289 2333 5447 7637 364 211 success 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 warmstart ǫ (25, 10) 2 8 2 8 8 343 228 151 15 37 24 41 37 7513 1533 302 2144 3782 409 36483 392 64 success 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 dkmax warmstart ǫ (5, 5) 25 25 25 25 25 250 250 25 25 25 25 25 25 250 250 25 250 250 25 25 250 25 3 12 3 14 12 686 351 319 18 55 22 60 64 5244 5606 43051 43200 42900 5617 38204 276 46 success 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 1 1 1 1 dkmax 25 25 25 25 25 625 125 25 25 25 25 25 25 125 125 5 25 5 25 25 125 25 Appendix A. Tables Problem 72 Appendix A. Tables Table A.10: Objective value in $ for direct solve, ε = 0.01 Problem Objective success 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 454426 9351412 316359 489574 60684 209271 77114005 56893108 455163 10000000 376427 493254 60687 212554 77114005 0 0 0 70262 0 451922 9351412 297389 482635 60253 209129 61940820 55096639 452486 9480000 297389 484506 60686 210088 61940820 56100000 454000 300248 0 63540 0 0 56487153 208258 1 1 1 1 1 1 1 1 1 0 1 1 1 0 1 0 0 0 0 0 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 0 0 1 0 0 0 0 1 1 ε LP mass value 0.00680 0.00549 0.00000 0.00446 0.00833 0.00540 0.00198 0.00714 0.00338 0.08160 0.00123 0.00982 0.00396 0.01895 0.00198 1.00000 1.00000 1.00000 0.02304 1.00000 0.00983 0.00966 -0.00000 0.00000 0.00124 0.00559 0.00000 0.00737 0.00905 0.04150 0.00000 0.00539 0.00395 0.01000 0.00000 0.02480 0.01480 0.00952 1.00000 0.04868 1.00000 1.00000 0.00000 0.00073 455181 9974094 329040 505297 68898 211057 76961263 56926279 472891 10204508 394819 519193 72079 212554 76961263 56975485 477309 516378 79597 217812 448897 9512758 300402 505297 67322 209244 61940820 55096639 468367 9620263 339047 482671 68391 211502 61940820 55132871 475342 360154 488042 77395 217519 64716813 56487153 214188 73 Appendix A. Tables Table A.11: Objective value in $ for direct solve with x◦ , ε = 0.01 Problem Objective success 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 454426 9351412 316359 488986 60684 209271 76961263 56926279 455163 9498097 376427 493254 60687 212554 76961263 56975485 462423 495000 70262 218000 448897 9351412 297389 482635 60253 209244 61940820 55096639 452486 9430000 297389 482671 60687 210088 61940820 55132871 455000 300248 485640 63540 217519 64716813 56487153 210165 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 0 0 0 0 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 0 1 1 1 0 0 1 1 ε LP mass value 0.00861 0.00366 0.00205 0.00327 0.00813 0.00627 0.00000 0.00771 0.00584 0.00866 0.00591 0.00995 0.00396 0.01170 0.00000 0.00857 0.03118 0.01470 0.10922 0.04520 0.00316 0.00000 0.00000 0.00092 0.00124 0.00614 0.00000 0.00737 0.00933 0.03620 0.00000 0.00160 0.00396 0.00993 0.00000 0.00802 0.01560 0.00952 0.00771 0.00821 0.04395 0.04289 0.00000 0.00975 455181 9974094 329040 505297 68898 211057 76961263 56926279 472891 10204508 394819 519193 72079 212554 76961263 56975485 477309 516378 79597 217812 448897 9512758 300402 505297 67322 209244 61940820 55096639 468367 9620263 339047 482671 68391 211502 61940820 55132871 475342 360154 488042 77395 217519 64716813 56487153 214188 74 Appendix A. Tables Table A.12: Objective value in $ for warmstart (d0max = 25, α = 10), ε = 0.01 Problem Objective success 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 454426 9351412 316359 489574 60684 209271 76961263 56982170 455163 9500000 376427 493254 60687 212554 77114005 56975485 459000 496000 64484 218000 451922 9351412 297389 483187 60253 209244 61940820 55096639 452486 9430000 297389 482671 60687 210088 61940820 55132871 451067 300248 485640 63540 217519 64716813 56487153 210165 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 0 0 0 0 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 0 0 0 1 1 ε LP mass value 0.00843 0.00000 0.00000 0.00446 0.00833 0.00627 0.00000 0.00869 0.00452 0.03160 0.00346 0.00982 0.00396 0.00996 0.00198 0.00857 0.02130 0.01730 0.05960 0.04510 0.00983 0.00984 0.00000 0.00206 0.00124 0.00614 0.00000 0.00737 0.00999 0.03620 0.00000 0.00160 0.00396 0.00993 0.00000 0.00802 0.00795 0.00952 0.00771 0.04868 0.04395 0.04289 0.00000 0.00975 455181 9974094 329040 505297 68898 211057 76961263 56926279 472891 10204508 394819 519193 72079 212554 76961263 56975485 477309 516378 79597 217812 448897 9512758 300402 505297 67322 209244 61940820 55096639 468367 9620263 339047 482671 68391 211502 61940820 55132871 475342 360154 488042 77395 217519 64716813 56487153 214188 75 Appendix A. Tables Table A.13: Objective value in $ for warmstart with ǫ (d0max = 25, α = 10), ε = 0.01 Problem Objective success 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 454426 9351412 316359 488986 60684 209271 76961263 56982170 455163 9498097 376427 493254 60687 212554 77114005 56975485 458003 496000 70234 218000 448897 9351412 297389 482635 60253 209244 61940820 55096639 452486 9480000 297389 482671 60748 210088 61940820 55132871 451067 300320 484188 63540 217519 0 56487153 208258 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 0 0 0 1 1 ε ǫ LP mass value 0.00929 0.00591 -0.00000 0.00327 0.00027 0.00542 0.00000 0.00869 0.00338 0.00859 0.00074 0.00989 0.00396 0.00996 0.00198 0.00857 0.02183 0.01730 0.10560 0.04510 0.00316 0.00966 0.00000 0.00092 0.00124 0.00614 0.00000 0.00737 0.00926 0.04150 0.00000 0.00160 0.00496 0.00996 0.00000 0.00802 0.00795 0.00976 0.00473 0.04868 0.04395 1.00000 0.00000 0.00068 0.01434 0.02938 0.06244 0.00328 0.00840 0.00631 0.00000 0.00876 0.01598 0.04553 0.26417 0.01203 0.00398 0.02209 0.00198 0.00857 0.02183 0.01730 0.33415 0.04510 0.00317 0.02938 0.00000 0.00092 0.00124 0.00618 0.00000 0.00742 0.01119 0.04150 0.00000 0.00161 0.00498 0.01024 0.00000 0.00809 0.00802 0.00986 0.00473 0.05154 0.04858 1.00000 0.00000 0.00106 455181 9974094 329040 505297 68898 211057 76961263 56926279 472891 10204508 394819 519193 72079 212554 76961263 56975485 477309 516378 79597 217812 448897 9512758 300402 505297 67322 209244 61940820 55096639 468367 9620263 339047 482671 68391 211502 61940820 55132871 475342 360154 488042 77395 217519 64716813 56487153 214188 76 Appendix A. Tables Table A.14: Objective value in $ for warmstart with ǫ (d0max = 5, α = 5), ε = 0.01 Problem Objective success 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 454426 9351412 316359 488986 60684 209271 76961263 56576215 455163 9498097 376427 493254 60687 212554 77114005 56902685 656851 496000 82236 218000 451922 9351412 297389 483187 60253 208491 61940820 54767023 452486 9430000 297389 483954 60735 210088 61940820 55231802 451067 403569 538301 77847 218000 64716813 56487153 210165 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 0 0 0 0 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 0 0 0 0 0 1 1 ε ǫ LP mass value 0.00843 0.00593 0.00000 0.00327 0.00833 0.00627 0.00000 0.00157 0.00600 0.00952 0.00154 0.00865 0.00396 0.01436 0.00198 0.00730 0.05385 0.01710 0.06211 0.04550 0.00983 0.00000 -0.00000 0.00206 0.00124 0.00255 0.00000 0.00139 0.00905 0.03620 0.00000 0.00425 0.00474 0.00974 0.00000 0.00980 0.00795 0.06654 0.10452 0.07685 0.04510 1.00000 -0.00000 0.00975 0.01434 0.02938 0.06244 0.00328 0.00840 0.00631 0.00000 0.00158 0.01598 0.04553 0.26417 0.01203 0.00398 0.05611 0.00198 0.00730 0.52950 0.01761 0.55954 0.04771 0.00993 0.02938 0.00000 0.00207 0.00124 0.00256 0.00000 0.00140 0.01119 0.03755 0.00000 0.00427 0.00477 0.01024 0.00000 0.00989 0.00802 0.60979 1.00000 0.45749 0.04801 0.04482 0.00000 0.01023 455181 9974094 329040 505297 68898 211057 76961263 56926279 472891 10204508 394819 519193 72079 212554 76961263 56975485 477309 516378 79597 217812 448897 9512758 300402 505297 67322 209244 61940820 55096639 468367 9620263 339047 482671 68391 211502 61940820 55132871 475342 360154 488042 77395 217519 64716813 56487153 214188 77 Appendix A. Tables Table A.15: Objective value in $ for direct solve, ε = 0.05 Problem Objective success 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 457354 9351412 316359 489574 60684 209271 77114005 56893108 458942 10000000 376427 494620 61213 212554 77114005 0 0 0 70262 0 451922 9351412 297389 503968 60253 209129 61940820 55096639 453208 9556950 297389 487393 60686 210088 61940820 56082250 468035 301565 0 63562 0 0 56487153 208258 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 0 0 0 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 0 0 1 1 ε LP mass value 0.01905 0.00549 0.04771 0.00446 0.00833 0.00540 0.00198 0.00714 0.02107 0.08160 0.02126 0.01462 0.01252 0.01895 0.00198 1.00000 1.00000 1.00000 0.02951 1.00000 0.00983 0.02854 -0.00000 0.04316 0.00124 0.00559 0.00000 0.00737 0.01264 0.04943 0.00000 0.01128 0.00395 0.01013 0.00000 0.02481 0.04392 0.01385 1.00000 0.04901 1.00000 1.00000 0.00000 0.00073 455181 9974094 329040 505297 68898 211057 76961263 56926279 472891 10204508 394819 519193 72079 212554 76961263 56975485 477309 516378 79597 217812 448897 9512758 300402 505297 67322 209244 61940820 55096639 468367 9620263 339047 482671 68391 211502 61940820 55132871 475342 360154 488042 77395 217519 64716813 56487153 214188 78 Appendix A. Tables Table A.16: Objective value in $ for direct solve with x◦ , ε = 0.05 Problem Objective success 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 455181 9351412 316359 505297 60684 211057 76961263 56926279 466699 9498097 376427 500476 60687 212554 76961263 56975485 466407 516000 70262 217812 448897 9512758 300402 505297 60253 209244 61940820 55096639 468367 9556950 297389 482671 60687 211502 61940820 55132871 463946 301010 488042 63607 217519 64716813 56487153 214188 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ε LP mass value 0.01577 0.02854 0.04304 0.03544 0.00813 0.01468 0.00000 0.00771 0.03908 0.04354 0.03718 0.02615 0.00396 0.02161 0.00000 0.00857 0.03946 0.05600 0.10922 0.04523 0.00316 0.04502 0.01003 0.04573 0.00124 0.00614 0.00000 0.00737 0.04460 0.04943 0.00000 0.00160 0.00396 0.01675 0.00000 0.00802 0.03549 0.01203 0.01259 0.04969 0.04395 0.04289 0.00000 0.02871 455181 9974094 329040 505297 68898 211057 76961263 56926279 472891 10204508 394819 519193 72079 212554 76961263 56975485 477309 516378 79597 217812 448897 9512758 300402 505297 67322 209244 61940820 55096639 468367 9620263 339047 482671 68391 211502 61940820 55132871 475342 360154 488042 77395 217519 64716813 56487153 214188 79 Appendix A. Tables Table A.17: Objective value in $ for warmstart (d0max = 25, α = 10), ε = 0.05 Problem Objective success 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 454426 9351412 316359 505297 60684 211057 78347224 56982170 461811 9498097 376427 497561 61213 212554 78347224 56975485 461947 495973 64484 218355 461946 9512758 300402 505297 60253 209244 61940820 55096639 456838 9425696 301197 482671 61213 211502 61940820 55132871 455873 312180 488042 63540 218062 64716813 57840191 210680 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ε LP mass value 0.01413 0.00000 0.04485 0.03544 0.00833 0.01468 0.01769 0.00869 0.02990 0.04354 0.00346 0.02044 0.01252 0.02161 0.01769 0.00857 0.03018 0.01731 0.05960 0.04761 0.03132 0.04502 0.01003 0.04573 0.00124 0.00614 0.00000 0.00737 0.02048 0.03620 0.01264 0.00160 0.01252 0.01675 0.00000 0.00802 0.01841 0.04738 0.01259 0.04868 0.04633 0.04289 0.02339 0.01254 455181 9974094 329040 505297 68898 211057 76961263 56926279 472891 10204508 394819 519193 72079 212554 76961263 56975485 477309 516378 79597 217812 448897 9512758 300402 505297 67322 209244 61940820 55096639 468367 9620263 339047 482671 68391 211502 61940820 55132871 475342 360154 488042 77395 217519 64716813 56487153 214188 80 Appendix A. Tables Table A.18: Objective value in $ for warmstart with ǫ (d0max = 25, α = 10), ε = 0.05 Problem Objective success 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 454426 9351412 316359 505297 60684 211057 78347224 56982170 461811 9498097 376427 497561 61213 212554 78347224 56975485 461947 495973 64484 218355 461946 9512758 300402 505297 60253 209244 61940820 55096639 456838 9425696 301197 482671 61213 211502 61940820 55132871 455873 312180 488042 63562 218062 64716813 57840191 210680 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ε ǫ LP mass value 0.01413 0.00000 0.04485 0.03544 0.00833 0.01468 0.00000 0.00869 0.02990 0.02020 0.00346 0.02044 0.01252 0.02161 0.00000 0.00857 0.03018 0.01709 0.05960 0.04761 0.03132 0.04502 0.01003 0.04573 0.00124 0.00614 0.00000 0.00737 0.02048 0.03620 0.01264 0.00160 0.01252 0.01675 0.00000 0.00802 0.01841 0.04738 0.01259 0.04901 0.04633 0.04289 0.02339 0.01216 0.01434 0.02938 0.06244 0.03674 0.00840 0.01489 0.01801 0.00876 0.03082 0.04553 0.26417 0.02087 0.01267 0.02209 0.01801 0.00857 0.03112 0.01761 0.06679 0.04999 0.03233 0.04714 0.01013 0.04792 0.00124 0.00618 0.00000 0.00742 0.02091 0.03756 0.01281 0.00161 0.01267 0.01703 0.00000 0.00809 0.01876 0.04974 0.01259 0.05154 0.04858 0.04482 0.02395 0.01270 455181 9974094 329040 505297 68898 211057 76961263 56926279 472891 10204508 394819 519193 72079 212554 76961263 56975485 477309 516378 79597 217812 448897 9512758 300402 505297 67322 209244 61940820 55096639 468367 9620263 339047 482671 68391 211502 61940820 55132871 475342 360154 488042 77395 217519 64716813 56487153 214188 81 Appendix A. Tables Table A.19: Objective value in $ for warmstart with ǫ (d0max = 5, α = 5), ε = 0.05 Problem Objective success 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 454426 9351412 316359 488986 60684 211057 78347224 56576215 461811 9498097 376427 501295 60687 219628 78499966 56975485 656851 495973 82236 217881 461946 9351412 297389 488396 60253 208491 61940820 57438725 465460 9425696 297389 490197 60735 211607 61940820 55231802 455042 403569 560000 77847 217944 64716813 57205059 218432 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 1 1 1 1 ε ǫ LP mass value 0.01413 0.02854 0.04655 0.00327 0.00833 0.01468 0.00000 0.00157 0.02990 0.04340 0.03643 0.02774 0.00396 0.04992 0.00195 0.00857 0.05385 0.01709 0.06379 0.04553 0.03132 0.02811 -0.00000 0.01271 0.00124 0.00255 0.00000 0.04784 0.03863 0.03620 0.00000 0.01693 0.00474 0.01724 0.00000 0.00980 0.01662 0.06654 1.00000 0.07685 0.04581 0.04289 0.01255 0.04759 0.01434 0.02938 0.06244 0.03674 0.00840 0.01489 0.01801 0.00876 0.03082 0.04553 0.26417 0.02087 0.01267 0.02209 0.01801 0.00857 0.03112 0.01761 0.06679 0.04999 0.03233 0.04714 0.01013 0.04792 0.00124 0.00618 0.00000 0.00742 0.02091 0.03756 0.01281 0.00161 0.01267 0.01703 0.00000 0.00809 0.01876 0.04974 0.01259 0.05154 0.04858 0.04482 0.02395 0.01270 455181 9974094 329040 505297 68898 211057 76961263 56926279 472891 10204508 394819 519193 72079 212554 76961263 56975485 477309 516378 79597 217812 448897 9512758 300402 505297 67322 209244 61940820 55096639 468367 9620263 339047 482671 68391 211502 61940820 55132871 475342 360154 488042 77395 217519 64716813 56487153 214188 82 Appendix B Figures 5% mipgap Problems 1% mipgap 20 20 15 15 10 10 5 5 0 -0.2 -0.1 0 0.1 0.2 0.3 0 -0.2 -0.1 0 0.1 0.2 0.3 Figure B.1: Savings with direct solve, compared to mass block removal schedule 83 Appendix B. Figures 5% mipgap 1% mipgap 35 25 30 Problems 20 25 15 20 10 15 10 5 5 0 -0.2-0.15-0.1-0.05 0 0.05 0.1 0.15 0.2 0 -0.2-0.15-0.1-0.05 0 0.05 0.1 0.15 0.2 Figure B.2: Savings with direct solve and x◦ , compared to mass block removal schedule 5% mipgap Problems 1% mipgap 25 25 20 20 15 15 10 10 5 5 0 -0.2 -0.1 0 0.1 0.2 0.3 0 -0.2-0.15-0.1-0.05 0 0.05 0.1 0.15 0.2 Figure B.3: Savings with warmstart (d0max = 25, α = 10), compared to mass block removal schedule 5% mipgap Problems 1% mipgap 20 20 15 15 10 10 5 5 0 -0.1 -0.05 0 0.05 0.1 0.15 0.2 0 -0.2-0.15-0.1-0.05 0 0.05 0.1 0.15 0.2 Figure B.4: Savings with warmstart and ǫ (d0max = 5, α = 5), compared to mass block removal schedule 84 Appendix B. Figures 1% mipgap 5% mipgap 0.2 0.2 0.15 0.15 0.1 0.1 0.05 0.05 0 0 -0.05 -0.05 0 5 10 15 20 25 30 35 40 0 5 10 15 20 25 30 35 40 Figure B.5: Stratification of gains for direct solve; * are 5, ⊙ are 10, are 20 blocks; small symbols are 2, big symbols are 4 access roads; within the block groups, points are sorted from left to right in ascending order of section numbers 1% mipgap 5% mipgap 0.2 0.2 0.15 0.15 0.1 0.1 0.05 0.05 0 0 0 10 20 30 40 50 0 10 20 30 40 50 Figure B.6: Stratification of gains for direct solve with x◦ ; green=5, blue=10, red=20 blocks; * are 5, ⊙ are 10, are 20 blocks; small symbols are 2, big symbols are 4 access roads; within the block groups, points are sorted from left to right in ascending order of section numbers 85 Appendix B. Figures 1% mipgap 5% mipgap 0.2 0.2 0.15 0.15 0.1 0.1 0.05 0.05 0 0 -0.05 -0.05 0 10 20 30 40 50 0 10 20 30 40 50 Figure B.7: Stratification of gains for warmstart (d0max = 25, α = 10); * are 5, ⊙ are 10, are 20 blocks; small symbols are 2, big symbols are 4 access roads; within the block groups, points are sorted from left to right in ascending order of section numbers 5% mipgap 1% mipgap 0.2 0.2 0.15 0.15 0.1 0.1 0.05 0.05 0 0 -0.05 -0.05 0 5 10 15 20 25 30 35 40 0 5 10 15 20 25 30 35 40 Figure B.8: Stratification of gains for warmstart with ǫ (d0max = 5, α = 5); * are 5, ⊙ are 10, are 20 blocks; small symbols are 2, big symbols are 4 access roads; within the block groups, points are sorted from left to right in ascending order of section numbers 86
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Optimizing earthwork block removal in road construction
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Optimizing earthwork block removal in road construction Koch, Valentin Raphael 2010
pdf
Page Metadata
Item Metadata
Title | Optimizing earthwork block removal in road construction |
Creator |
Koch, Valentin Raphael |
Publisher | University of British Columbia |
Date Issued | 2010 |
Description | In road construction, earthwork operations account for about 25% of the construction costs. Existing linear programming models for earthwork logistics optimization are designed to minimize the hauling costs and to balance the earth across the construction site. However, these models do not consider the removal of physical blocks that may influence the earthwork process. In this thesis, we extend the linear programming model of Mayer and Stark (1981) with the addition of a block removal schedule. The resulting model is a mixed-integer linear program. We analyze the model size and the schedule search space in order to make conclusion about the use of the model. Based on structural observations, we introduce a set of algorithms that significantly reduce the solving time of the model. Finally, we conduct numerical experiments to compare our solutions with the solutions of a traditional earthwork process that makes use of linear programming. From our numerical results, we conclude that an optimal removal schedule produces solutions that are 4.1% cheaper on average than a traditional method, with savings that can go as high as 19%. We conclude our discussion with possible extensions to the model, that can help an engineer to design roads that are more economical and ecological with respect to the earthwork operations. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-04-16 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivs 3.0 Unported |
DOI | 10.14288/1.0070953 |
URI | http://hdl.handle.net/2429/23722 |
Degree |
Master of Science - MSc |
Program |
Mathematics |
Affiliation |
Irving K. Barber School of Arts and Sciences (Okanagan) Mathematics, Department of (Okanagan) |
Degree Grantor | University of British Columbia |
Graduation Date | 2010-05 |
Campus |
UBCO |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/3.0/ |
Aggregated Source Repository | DSpace |
Download
- Media
- 24-ubc_2010_spring_koch_valentin.pdf [ 590.06kB ]
- Metadata
- JSON: 24-1.0070953.json
- JSON-LD: 24-1.0070953-ld.json
- RDF/XML (Pretty): 24-1.0070953-rdf.xml
- RDF/JSON: 24-1.0070953-rdf.json
- Turtle: 24-1.0070953-turtle.txt
- N-Triples: 24-1.0070953-rdf-ntriples.txt
- Original Record: 24-1.0070953-source.json
- Full Text
- 24-1.0070953-fulltext.txt
- Citation
- 24-1.0070953.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0070953/manifest