Optimizing Earthwork Block Removalin Road ConstructionbyValentin Raphael KochB.Sc. Hons., The University of British Columbia, 2008A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF SCIENCEinThe College of Graduate Studies(Mathematics)THE UNIVERSITY OF BRITISH COLUMBIA(Okanagan)April 2010c©Valentin Raphael Koch 2010AbstractIn road construction, earthwork operations account for about25% of thecon-struction costs. Existing linear programming models for earthwork logisticsoptimization are designed to minimize the hauling costs and to balance theearth across the construction site. However, these models do not considerthe removal of physical blocks that may influence the earthwork process.In this thesis, we extend the linear programming model of Mayer andStark [MS81] with the addition of a block removal schedule. The resultingmodel is a mixed-integer linear program. We analyze the model size andthe schedule search space in order to make conclusion about the use of themodel. Based on structural observations, we introduce a set of algorithmsthat significantly reduce the solving time of the model. Finally, we con-duct numerical experiments to compare our solutions with the solutions ofa traditional earthwork process that makes use of linear programming.From our numerical results, we conclude that an optimal removal sched-ule produces solutions that are 4.1% cheaper on average than a traditionalmethod, with savings that can go as high as 19%. We conclude our dis-cussion with possible extensions to the model, that can help an engineer todesign roads that are more economical and ecological with respect to theearthwork operations.iiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viGlossary of notation . . . . . . . . . . . . . . . . . . . . . . . . . . viiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Road design . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1.1 Road design and construction . . . . . . . . . . . . . 11.1.2 Road design optimization . . . . . . . . . . . . . . . . 21.1.3 Earthwork optimization . . . . . . . . . . . . . . . . . 41.2 Basic linear programming model . . . . . . . . . . . . . . . . 71.2.1 Notation . . . . . . . . . . . . . . . . . . . . . . . . . 71.2.2 Model parameters . . . . . . . . . . . . . . . . . . . . 71.2.3 Decision variables . . . . . . . . . . . . . . . . . . . . 91.2.4 Objective and constraints . . . . . . . . . . . . . . . . 102 Mixed integer programming model . . . . . . . . . . . . . . . 132.1 Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142.2 Model parameters . . . . . . . . . . . . . . . . . . . . . . . . 142.3 Decision variables . . . . . . . . . . . . . . . . . . . . . . . . 152.4 Objective and constraints . . . . . . . . . . . . . . . . . . . . 162.4.1 Objective function . . . . . . . . . . . . . . . . . . . . 172.4.2 Volume constraints . . . . . . . . . . . . . . . . . . . 172.4.3 Block removal indicator . . . . . . . . . . . . . . . . . 17iiiTable of Contents2.4.4 Move feasibility constraint . . . . . . . . . . . . . . . 182.4.5 Stockpile constraints . . . . . . . . . . . . . . . . . . 182.4.6 Block removal enforcement . . . . . . . . . . . . . . . 202.4.7 Observations . . . . . . . . . . . . . . . . . . . . . . . 203 Model analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 223.1 Space complexity . . . . . . . . . . . . . . . . . . . . . . . . 223.2 Schedule search . . . . . . . . . . . . . . . . . . . . . . . . . 264 Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 354.1 Mass schedule as starting point . . . . . . . . . . . . . . . . . 364.2 Maximal allowed moving distance dmax . . . . . . . . . . . . 374.3 Warmstarting with increasing dmax . . . . . . . . . . . . . . 394.4 Warmstarting with mipgap check . . . . . . . . . . . . . . . 405 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . 425.1 Experimental setup . . . . . . . . . . . . . . . . . . . . . . . 425.1.1 Problems . . . . . . . . . . . . . . . . . . . . . . . . . 425.1.2 Benchmark . . . . . . . . . . . . . . . . . . . . . . . . 425.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 465.2.1 Analysis of the dmax-heuristic . . . . . . . . . . . . . 465.2.2 Performance . . . . . . . . . . . . . . . . . . . . . . . 485.2.3 Economical analysis . . . . . . . . . . . . . . . . . . . 546 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 576.1 Accomplishments . . . . . . . . . . . . . . . . . . . . . . . . 576.2 Research outcomes . . . . . . . . . . . . . . . . . . . . . . . . 576.3 Future research . . . . . . . . . . . . . . . . . . . . . . . . . 596.3.1 Warmstart schemes and εLP . . . . . . . . . . . . . . 596.3.2 Vertical alignment . . . . . . . . . . . . . . . . . . . . 596.3.3 Schedule enumeration . . . . . . . . . . . . . . . . . . 606.3.4 Lagrangian relaxation . . . . . . . . . . . . . . . . . . 61Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62AppendicesA Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66B Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83ivList of Tables1.1 Major components of construction costs . . . . . . . . . . . . 32.1 Stockpile example data . . . . . . . . . . . . . . . . . . . . . . 183.1 Schedule example data . . . . . . . . . . . . . . . . . . . . . . 275.1 Test configurations for dmax . . . . . . . . . . . . . . . . . . . 455.2 Economical test and performance test configurations . . . . . 45A.1 Problem configurations . . . . . . . . . . . . . . . . . . . . . . 66A.2 Statistics for objective function error with dmax . . . . . . . . 67A.3 Statistics for solving time with dmax . . . . . . . . . . . . . . 67A.4 Histogram statistics for 1% tolerance . . . . . . . . . . . . . . 68A.5 Histogram statistics for 5% tolerance . . . . . . . . . . . . . . 68A.6 CPU time in seconds for ε = 0.01, problems 1-22 . . . . . . . 69A.7 CPU time in seconds for ε = 0.01, problems 23-44 . . . . . . 70A.8 CPU time in seconds for ε = 0.05, problems 1-22 . . . . . . . 71A.9 CPU time in seconds for ε = 0.05, problems 23-44 . . . . . . 72A.10 Objective value in $ for direct solve, ε = 0.01 . . . . . . . . . 73A.11 Objective value in $ for direct solve with x◦, ε = 0.01 . . . . . 74A.12 Objective value in $ for warmstart (d0max = 25,α = 10), ε = 0.01 75A.13 Objective value in $ for warmstart with ǫ (d0max = 25,α = 10),ε = 0.01 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76A.14 Objective value in $ for warmstart with ǫ (d0max = 5,α = 5),ε = 0.01 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77A.15 Objective value in $ for direct solve, ε = 0.05 . . . . . . . . . 78A.16 Objective value in $ for direct solve with x◦, ε = 0.05 . . . . . 79A.17 Objective value in $ for warmstart (d0max = 25,α = 10), ε = 0.05 80A.18 Objective value in $ for warmstart with ǫ (d0max = 25,α = 10),ε = 0.05 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81A.19 Objective value in $ for warmstart with ǫ (d0max = 5,α = 5),ε = 0.05 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82vList of Figures1.1 Example of vertical alignment with mass diagram . . . . . . . 51.2 Example of vertical alignment with sections . . . . . . . . . . 82.1 Earthwork scheme with removal schedule (black squares in-dicate blocks, gray squares indicate removed blocks). . . . . . 193.1 Potential schedules . . . . . . . . . . . . . . . . . . . . . . . . 273.2 Earthwork scheme with removal schedule. . . . . . . . . . . . 283.3 Schedule tree . . . . . . . . . . . . . . . . . . . . . . . . . . . 293.4 (6×4)-lattice in Z2 . . . . . . . . . . . . . . . . . . . . . . . 305.1 Objective and solving time with increasing dmax . . . . . . . . 475.2 Algorithms performance profiles with log2 scale . . . . . . . . 505.3 Algorithms performance profiles with small τ . . . . . . . . . 525.4 Algorithms performance profiles with large τ . . . . . . . . . 535.5 Savings with warmstart with ǫ (d0max = 25,α = 10) . . . . . . 555.6 Stratification for warmstart with ǫ (d0max = 25,α = 10) . . . . 56B.1 Savings with direct solve . . . . . . . . . . . . . . . . . . . . . 83B.2 Savings with direct solve and x◦ . . . . . . . . . . . . . . . . 84B.3 Savings with warmstart (d0max = 25,α = 10) . . . . . . . . . . 84B.4 Savings with warmstart and ǫ (d0max = 5,α = 5) . . . . . . . . 84B.5 Stratification for direct solve . . . . . . . . . . . . . . . . . . 85B.6 Stratification for direct solve with x◦ . . . . . . . . . . . . . . 85B.7 Stratification for warmstart (d0max = 25,α = 10) . . . . . . . . 86B.8 Stratification for warmstart with ǫ (d0max = 5,α = 5) . . . . . 86viGlossary of notationa∈A a is an element of setA. 4R Real numbers. 4R+ Nonnegative real numbers. 5φ :A→B The mappingAintoB by φ. 5[a,b] Interval{x∈R : a≤x≤b}. 5∞ Infinity. 5f′ Derivative of function f. 6integraltextIntegral. 5|·| Absolute value. 7|A| Cardinality of setA. 7{x : P(x)} The set of all x such that P(x) is satisfied. 8Z+ Nonnegative integers. 8A∪B The union of setsAandB. 8dist(i,j) Distance between section i and j. 8R++ Positive real numbers. 8minf The minimum of function f. 10summationtextSummation. 10A×B Cartesian product{(x,y) : x∈A,y∈B}. 12∧ Logical connective representing and. 14∨ Logical connective representing or. 14A\B Set difference{x : x∈Aand xnegationslash∈B}. 16∅ The empty set. 16Θ(·) Asymptotic tight bound. 22O(·) Asymptotic upper bound. 22Z++ Positive integers. 22A⊆B Ais a subset ofB. 22parenleftbignkparenrightbig Binomial coefficient (n choose k). 30Rn Real n-vectors (n×1 matrices). 35XT Transpose of matrix X. 35Rm×n Real m×n matrices. 38a←b Assign value b to a. 37viiAcknowledgementsMany people were involved in the creation of this thesis. First, I would liketo express my sincere gratitude to my supervisors, Dr. Yves Lucet and Dr.Warren Hare, for their help and support throughout this research. I greatlyenjoyed our lively discussions and conversations. It was a great privilege tohave supervisors that helped me to overcome so many challenges, whetherthey 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 guid-ance in my graduate studies, and for providing me insight into the beautifulsubject of convex analysis.It is a pleasure to thank my colleagues, especially Jim Nastos and MasonMacklem, for their patience to listen and the many ideas and input theycontributed.A very special thanks goes out to our industrial partner, for the profes-sional assistance and valuable help on engineering questions.I owe my deepest gratitude to my parents. Without their prayers andsupport, my academic career would have been impossible.Finally, I would like to thank my wife Delphine, and my children Jordanand Tabitha, for their love and patience during these intense years.This work was partially supported by the Ministry of Advanced Educa-tion of the Province of British Columbia through a Pacific Century GraduateScholarship, and by the Mathematics of Information Technology and Com-plex Systems network through a MITACS Accelerate internship.viiiDedicationTo Delphine, Jordan, and Tabitha.ixChapter 1IntroductionLa tendance la plus profonde detoute activit´e humaine est lamarche vers l’´equilibre.1Jean Piaget (1896–1980)1.1 Road design1.1.1 Road design and constructionTransportation of people and goods over land is, and has always been, animportant part of civilization. Construction of the first roadways dates backas early as 5000 BCE [Lay99]. The invention of the wheel created the needfor better roads and more sophisticated forms of road construction. Complexroad design, including the alignment, earthwork, and pavement was alreadyperformed during the Roman Empire. The first scientific approach to roaddesign was by Tr´esaguet in France at about 1764 [Lay99]; Tr´esaguet alsointroduced the first system of continuous road maintenance. The first com-prehensive guide on road design and construction was published by McAdamin 1816 [McA24], introducing the concept of the elevation of roads throughembankment. This resulted in much more stable roads, but also increasedthe cost of construction.Roadways are an essential part of the economy in most countries. For ex-ample, in Canada in 2008, trucking accounted for 54% of the value of tradewith the United States and commercial transportation services accountedfor 4.1% of Canada’s value-added GDP [Stu09]. At the same time, all levelsof Canadian government combined spent $20.9 billion, or 1.6% of Canada’sGDP, to build and maintain its 1 million kilometer road network. In viewof these numbers, it is a natural consequence that governments of coun-tries that depend heavily on road networks have a vast interest in keeping a1The most profound tendency of all Human activity is progression toward equilibrium.11.1. Road designhigh quality road infrastructure, while at the same time trying to lower theexpenditures on construction and maintenance. From an economic perspec-tive, it is therefore important to try to minimize the cost and maximize thequality and safety of the road infrastructure.The design of a road largely depends on safety requirements, which inturn depend on the traffic volume and vehicle speed. In 2007, in Canada,about 140,000 road collisions occurred, with about 16,000 seriously injuredpersons and about 3,000 fatalities. More than 100 million tonnes of dan-gerous 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 greatimpact on the environment. Maki [MKV02] shows that poorly plannedroad construction, especially in sensitive areas like the Amazon region, of-ten harms the local economy instead of helping it. According to Maki, it isimportant to optimize the land use through environmental and social impactanalysis. Appropriate classification and zoning of land should therefore betaken into account, and environmental impacts, like deforestation or earth-work, should be minimized when a road is built.It follows from the above discussion that the road infrastructure needs tobe optimized with respect to economic and social costs, safety requirements,and environmental impact.1.1.2 Road design optimizationMultiple factors determine the total transportation cost of a road. Jha et.al. [JSJ06] list five main types of costs that affect the total transportationcost:− 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 de-sign optimization. In [CGF89], Chew et. al. divide construction costs insix major categories, shown in Table 1.1 with an approximate associated21.1. Road designTable 1.1: Major components of construction costsPavement 30%Earthwork (half of this is haul cost) 25%Bridges 20%Drainage 10%Miscellaneous items 10%Land 5%percentage of the total construction cost. Operation and maintenance costsinclude roadway surface, roadside features, bridges, tunnels, snow and icecontrol, and traffic control devices and amount at about 5% of the construc-tion costs, over a period of 30 years. The user costs, which include vehicleoperating costs, travel time, and traffic accident costs, range from 300% to1000% of the construction costs over a period of 30 years [JSJ06]. Finally,social and environmental costs can be the most critical factor in highwaydesign. It is, however, difficult to quantify these costs and they are usuallyconsidered before the optimization of the design.The design of the road can roughly be divided into geometric design andpavement design. The geometric design determines the alignment of the roadand affects the construction costs, user costs, and social and environmentalcosts. The pavement design affects operation and maintenance, as well asuser costs. Optimization of the geometric design is traditionally done in twostages: the horizontal alignment design and the vertical design.The horizontal alignment is the road trajectory from a satellite’s eyeview. The primary objects that influence a horizontal alignment are forbid-den areas like lakes, ocean shores, steep mountains, nature reserves, and pri-vate property. Cost parameters in horizontal alignments might include landcosts and road distance. Early attempts at horizontal alignment optimiza-tion were based on dynamic programing [Tri87]. More recent approachesuse genetic algorithms [JK06] or local neighborhood heuristics [LTL09].Once the horizontal alignment is done, one can draw a line from thestarting point to the end point of the road and look at the vertical groundprofile along this line. The vertical alignment optimization tries to fit aroad to the vertical ground profile as close as possible, while respectinggrade constraints and other road specifications. The closer the road is tothe ground profile, the less earthwork needs to be done in order to cut orfill sections of the road. Usually, vertical alignment optimization includesearthwork allocation and earth-moving optimization. In civil engineering,31.1. Road designthe vertical alignment optimization is commonly done through a manualiteration of computations that estimate the earthwork volumes for a givenalignment. The alignment is then changed and the volumes are recomputeduntil a satisfying solution is found. Traditional mathematical models thatautomate the optimization of the vertical alignment use the method of linearprogramming (LP) [Eas88b, Mor96, Mor09, KL10].Recently, the road alignment optimization research community has fo-cused on the simultaneous optimization of horizontal and vertical alignment.In [CGF89], Chew et. al. modeled the road as a three-dimensional cubicspline. Jong and Schonfeld [JS03] use genetic algorithms to first solve for thethree-dimensional alignment as a three-dimensional piecewise linear spline,and then fit circular curves to the horizontal and vertical projections of thespline. Parametric curves were studied by Makanae [Mak99] to model theroad in three dimensions. Cheng and Lee [CL06] use a heuristic bi-levelapproach to iteratively optimize horizontal and vertical alignments.All of the above optimization methods depend on the cost evaluation ofthe earthwork. In this thesis, we focus on the aspect of earthwork optimiza-tion for a given alignment. The optimal cost of earthwork is essential forthe assessment of a given road alignment.1.1.3 Earthwork optimizationEarthwork depends significantly on the terrain and the geometric design ofthe road. A road in a mountainous area requires more cutting and fillingof earth than a road on the plain. An optimal vertical grade line minimizesthe earthwork and satisfies all the design constraints for the road. To deter-mine the optimality with respect to the earthwork, one needs to be able tocompute the earthwork volumes for a given alignment.In order to compute the earthwork volumes for a given grade line, cross-sections of the road profile are taken at regular intervals. The point alongthe horizontal axis where the cross-section is taken is called a station. Sta-tions are usually arranged at an equal distance from one another along thehorizontal axis of the grade line. The distance Lij is the actual length be-tween station i and j along the horizontal center line of the road. The areaAi of a cross-section at station i is determined and the volume V ∈R be-tween two stations i and j is computed as V = Lij(Ai+Aj)/2. A traditionalapproach uses the mass diagram to depict the net accumulation of cut orfill 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 accumulation41.1. Road designroad profile, g(¯x)ground profile, f(¯x)cumulative mass, h(¯x) =integraldisplay ¯x0f(u)−g(u)dubalance point balance point¯xFigure 1.1: Example of vertical alignment with mass diagramin cubic meter from an arbitrary starting point. Thus, the dif-ference in ordinates between any two stations represents the netaccumulation of cut or fill between these stations. If the firststation 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 cantherefore 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 theheight of the road profile at ¯x, where f : R+ → R and g : R+ → R arepiecewise continuous functions on [0,+∞[, then the mass diagram is thegraph of the Riemann integralh(¯x) =integraldisplay ¯x0(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 thecorresponding mass diagram beneath. From Definition 1.1 and Figure 1.1,we observe the following characteristics.51.1. Road design− If h′(¯x) < 0, then ¯x is the position of a fill. If h′(¯x) > 0, then ¯x is theposition of a cut.− Any extremum of h(¯x) represents a switch from cut to fill, or from fillto cut.− The net accumulation between a station ¯xi and a station ¯xj, where¯xi < ¯xj is h(¯xj)−h(¯xi).− For any two stations ¯xi and ¯xj, if h(¯xi) = h(¯xj), then the net accumu-lation between the two stations is zero. These points are called balancepoints, since they show points of balance between cut and fill volumes.− A consequence of the last point is that if h(0) = 0 and ¯xend is thelast station of the road, then an optimal alignment with respect to thebalance of earth quantities has h(¯xend) = 0.With the above observations, we can derive a rule of thumb to determinethe hauling directions.Remark 1.2. Let h(0) = 0. If h(¯x) > 0 and ¯xi < ¯x < ¯xj, where ¯xi and¯xj are balance points, then earth needs to be hauled ahead. The distancebetween ¯xi and ¯xj is the maximum push or haul distance. If h(¯x) < 0 and¯xi < ¯x < ¯xj, then earth needs to be hauled back. The values at the extremarepresent the volumes of material that must be hauled along the road.The mass diagram and the above rule of thumb have several shortcom-ings.− The mass diagram does not take into account haul costs that are notproportional to the haul distance.− The mass diagram does not take into account different soil character-istics.− The mass diagram does not take into account borrow pits or waste pitsthat can be used to borrow additional quantities of soil or dispose ofexcessive 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 also61.2. Basic linear programming modelminimizes the hauling cost and incorporates soil material distribution, bor-row pits, and waste pits. Based on [MS81], Easa developed a mixed integerlinear program (MILP) that uses non-constant unit costs [Eas87], and aquadratic program (QP) that uses linear unit costs [Eas88a]. Jayawardaneand Harris [JH90] extended the model in [MS81] by incorporating projectduration 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 toovercome the problem of blocks in the haul route.1.2 Basic linear programming modelIn order to better explain our model of earth-moving, it is helpful to beginby considering the LP model in Mayer and Stark [MS81]. An extensivedescription of the model is given in [MS81] and [JH90]. We give a shorteroutline in this section with some changes to the notation and refer the readerto the original paper for further details.1.2.1 NotationFor 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 meaningof the notation should be clear from the context.1.2.2 Model parametersIn this section we provide information on the basic LP model parameters.Let a section be the region between two fixed stations along the horizontalaxis. We divide the vertical alignment into ms sections of equal length.Consider Definition 1.1. If for a section i, ¯xi represents the station at thebeginning of the section, and ¯xi+1 represents the station at the end of thesection, thenh(¯x) =integraldisplay ¯xi+1¯xi(f(u)−g(u)) duis the required change in volume between the ground profile and the roadprofile in the section. The required change in volume for a section is positiveif the section represents a cut, negative if the section represents a fill. Asin [MS81, Eas87, JH90], we introduce mb borrow pits and mw waste pitsin order to obtain a feasible solution in the case of an unbalanced verticalalignment. A borrow pit or a waste pit has a volume that represents thecapacity. The following input parameters are therefore given:71.2. Basic linear programming modelsectionroad profileground profile1 2 3 4 5 6Figure 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 indexborrow pits;− W = {ms + mb + 1,...,ms + mb + mw} is the set of indices that isused to index waste pits;− M=S∪B∪W.Finally, we need the required change in volume and the moving costs, whichare defined by the following parameters:− Vi ∈R is the required change in volume of a section i ∈S in cubicunits, 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 materialfrom i to j, where i,j∈M,inegationslash= j.In Figure 1.2, we show the alignment from Figure 1.1 that is divided intoms = 6 sections of equal length. We denote the distance between section iand j as dist(i,j), and evaluate this usingdist(i,j) =|¯xi−¯xj|, (1.1)where ¯xi is the midpoint of section i and ¯xj is the midpoint of section j.81.2. Basic linear programming modelRemark 1.3. The distance dist(i,j) can be defined arbitrarily. A more so-phisticated definition than Equation (1.1) would incorporate the distancebetween the sections on the horizontal alignment as well as the difference ofaltitude. We leave this as a future research topic.Note that there are many different ways to define the cost cij of movingearth. 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 ofmaterial from section i to section j, where i,j∈Mand inegationslash= j, is defined ascij =cexc + dist(i,j)chau +cemb, if i,j∈S,cbor + dist(i,j)chau +cemb, if i∈B,j∈S,cexc + dist(i,j)chau +cwas, if i∈S,j∈W,+∞, otherwise,(1.2)where cexc > 0 is the unit cost of excavation, chau > 0 is the cost to haulone cubic unit of material over a unit distance, cemb > 0 is the unit cost forembankment, cbor > 0 is the unit cost to borrow material from a borrowpit, 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 negationslash= k, inegationslash= j, andj negationslash= k, then the strict triangle inequality cij < cik + ckj holds whenevercij < +∞.The hauling distance from or to borrow pits or waste pits needs to ac-count for dead-haul distances (the distance of the road from the pit-joiningsection to the pit). Excavation costs on a section are usually smaller thanborrowing costs from a pit. This is usually balanced by the lower embank-ment costs of borrowed material. Also material that is dumped into a wastepit does not need a lot of embankment work and the embankment cost todump material into a pit is much smaller than the embankment costs on asection.1.2.3 Decision variablesThe decision variable xij ∈R represents the number of cubic units of earththat will be hauled from section i to section j, where i,j∈Mand inegationslash= j.91.2. Basic linear programming model1.2.4 Objective and constraintsThe minimization of the earth-moving costs is modeled by LP (1.3) belowminsummationdisplayi,j∈Minegationslash=jcijxij (1.3a)s.t.summationdisplayj∈Mjnegationslash=i(xij−xji) = Vi for all i∈S, (1.3b)summationdisplayj∈Mjnegationslash=ixij ≤Vi, for all i∈B, (1.3c)summationdisplayj∈Mjnegationslash=ixji≤Vi, for all i∈W, (1.3d)xij ≥0, for all i,j∈M,inegationslash= j. (1.3e)Objective functionThe objective function in (1.3a) minimizes the total hauling cost. Illogicalmoves (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.ConstraintsFor every road section i that has either a demand or a supply, we want theamount of incoming and outgoing earth movements to balance the demandor supply of the section. Constraints (1.3b) balance the demand and supplybetween 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. Weenforce the capacity constraints with inequalities (1.3c) and (1.3d). Sinceone 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 tozero. This is ensured by Equation (1.3e).101.2. Basic linear programming modelObservationsRemark 1.6. Note that shrinkand swell factors that are mentioned in [MS81]are not included in this thesis since the model works with banked units andthe mentioned factors are only taken into account to compute haul andcompaction costs.Proposition 1.7. Let Vs = summationtexti∈S Vi, Vb = summationtexti∈B Vi, and Vw = summationtexti∈W Vi. IfVs≥0 and Vw−Vs > 0, then the LP (1.3) has always a solution. If Vs < 0and 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 feasiblepoint exists, as the waste or borrow pit provide sufficient slack. Since for alli,j∈M, cij ∈R++∪{+∞}and xij ≥0, we know thatminsummationdisplayi,j∈Minegationslash=jcijxij ≥0.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 > 0and x∗kj > 0, where inegationslash= k and knegationslash= j, then x∗ is not an optimal solution.Proof. If x∗ik ≤x∗kj, consider ˜xik = 0, ˜xkj = x∗kj−x∗ik, and ˜xij = x∗ij + x∗ik.ThenVi = x∗ij + x∗ik +summationdisplays∈Ssnegationslash=i,j,kx∗is−summationdisplays∈Ssnegationslash=ix∗si = ˜xij + ˜xik +summationdisplays∈Ssnegationslash=i,j,kx∗is−summationdisplays∈Ssnegationslash=ix∗si,(1.4)Vk = x∗kj−x∗ik +summationdisplays∈Ssnegationslash=j,kx∗ks−summationdisplays∈Ssnegationslash=i,kx∗sk = ˜xkj−˜xik +summationdisplays∈Ssnegationslash=j,kx∗ks−summationdisplays∈Ssnegationslash=i,kx∗sk,(1.5)Vj=−x∗ij−x∗kj +summationdisplays∈Ssnegationslash=jx∗js−summationdisplays∈Ssnegationslash=i,j,kx∗sj =−˜xij−˜xkj +summationdisplays∈Ssnegationslash=jx∗js−summationdisplays∈Ssnegationslash=i,j,kx∗sj,(1.6)and the volume constraints are satisfied. But from Lemma 1.5 we know thatcik +ckj > cij, and thus the total cost is cheaper, becausecikx∗ik+ckjx∗kj+cijx∗ij > (cij−ckj)x∗ik+ckjx∗kj+cijx∗ij = cik˜xik+ckj˜xkj+cij˜xij.111.2. Basic linear programming modelIf x∗ik > x∗kj, consider ˜xik = x∗ik−x∗kj, ˜xkj = 0, and ˜xij = x∗ij + x∗kj. Then,the same arguments show that equations (1.4), (1.5), and (1.6) hold, andthus 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 numberof moves from any section i∈Mto 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 productM×M. Now|M×M|=|M|2.Since we do not count moves where i = j, we have |M|2 −|M| decisionvariables. The constraint in (1.3b) appears for every section in S. Also,constraints (1.3c), and (1.3d) appear for each borrow pit inBand each wastepit inW, 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|2constraints.Remark 1.10. In practice, since we do not implement illogical moves, thenumber of decision variables and constraints are less than the number givenin Proposition 1.9. If we count only logical moves, we have|S|2−|S|potentialmoves between sections,|S||B|potential moves between sections and borrowpits, and |S||W| potential moves between sections and waste pits. Thusin the actual implementation of the model, we have |S|(|M|−1) decisionvariables. 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 theextension of the basic LP model (1.3).12Chapter 2Mixed integer programmingmodelHe has barred my way withblocks of stone; he has made mypaths crooked.Lamentations 3:9 (Holy Bible,NIV)The model presented in Section 1.2 is a major improvement over thetypical mass diagram due to the minimization of hauling distances and theincorporation of borrow and waste pits. However, in practice, a move fromone section to another might not be possible due to natural blocks betweenthe 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, arock-bed, or some other topographical feature that needs to be removed orotherwise dealt with before earth can be moved over the location of theblock.In a traditional earthwork construction process, clearing of blocks isusually 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, asthe one described in Section 1.2, is perfectly valid, since after the obstacleremoval, all the sections are accessible. This, however, does not accountfor earthwork that may occur during the obstacle removal (i.e. filling of aditch or blasting of rock beds). During obstacle removal, stockpiles may becreated and material may be borrowed elsewhere. This makes it difficult toassess the overall costs of obstacle removal, excavation, embankment, andhaul in a traditional setting. Moreover, the order of removal of blocks mightaffect the final cost of the earthwork. We therefore extend the model in[MS81] with the introduction of blocks.132.1. Notation2.1 NotationFor the remainder of this thesis, the symbol∧is a logical connective repre-senting and, and the symbol∨is a logical connective representing or. HenceA∧B means “A and B”, where A∨B means “A or B”2.2 Model parametersWe extend the basic LP formulation in Section 1.2 to develop a model thatincorporates nb blocks.Definition 2.2. Let I = {1,...,nb}. Each block is assigned to a roadsection. We letς:I→Sbe 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 theroad. Each access road is connected to one of the ms sections. Access roadsare typically attached to sections at the beginning and end of the road, butmay also occur at any point along the road and in between blocks. Thenumber and position of access roads influences the optimal schedule for theremoval of the blocks.Definition 2.3. Let R = {1,...,nr} be the set of all access roads. Eachaccess road is assigned to a road section. We let̺:R→Sbe the mapping that maps the access road index to the assigned sectionindex.Since each access road allows for transporting earth out to an externalwaste pit or importing earth from an external borrow site, we attach anexternal borrow and waste pit to every access road with infinite capacityeach. These pits are accounted for in setB andW.We introduce discrete time steps to remove the blocks. The removal ofa block requires one time step and the block is removed at the end of thetime step. Multiple blocks may be removed within the same time step. Ina worst-case scenario, where all the blocks need to be removed in sequentialorder, nb time steps are needed to remove the blocks. Once all the blocks areremoved, 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}.142.3. Decision variablesThe 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 setcijt = cij.Excavation and embankment costs on external pits should be higher thanfor sections.Remark 2.4. If more than one access road is present, then any movement ofearth between two sections could also go over the access roads. One possibleway to allow such moves is to represent each access road as a section. Eachsection that has an access road attached could be cloned and the clone beadded to the set S. In this thesis, however, we do not account for suchmoves and we leave it as a further research topic.In Section 1.2.2, we stated that the temporary stocking of earth betweensections is avoided due to Lemma 1.5. With the introduction of blocks, how-ever, this is not guaranteed anymore (as we demonstrate in Example 2.7).We therefore introduce for each section i a stockpile tolerance parameterǫstki ∈R+. The significance and usage of this parameter is shown in Section2.4.5.2.3 Decision variablesSimilar to the cost coefficients, we keep the same variables representing thequantities of earth that are moved between sections. However, since we nowhave the flexibility of moving earth between any two sections at any timestep, 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 fromsection 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 toremove the block during that time step, and a one indicates that the blockwas 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 viaykt∈{0,1},such that if block k is still present at time t, then ykt = 0.152.4. Objective and constraints2.4 Objective and constraintsThe minimization of the earth-moving costs, including the removal of blocks,can therefore be stated as MILP (2.1) belowminsummationdisplayt∈Tsummationdisplayi,j∈Minegationslash=jcijtxijt (2.1a)s.t.summationdisplayt∈Tsummationdisplayj∈Mjnegationslash=i(xijt−xjit) = Vi for all i∈S, (2.1b)summationdisplayt∈Tsummationdisplayj∈Mjnegationslash=ixijt≤Vi, for all i∈B, (2.1c)summationdisplayt∈Tsummationdisplayj∈Mjnegationslash=ixjit≤Vi, for all i∈W, (2.1d)usummationdisplayt=0summationdisplayj∈Mjnegationslash=ς(k)(xς(k)jt−xjς(k)t)≥ykuVς(k) for all k∈I,u∈T,Vς(k)≥0,(2.1e)usummationdisplayt=0summationdisplayj∈Mjnegationslash=ς(k)(xς(k)jt−xjς(k)t)≤ykuVς(k) for all k∈I,u∈T,Vς(k) < 0,(2.1f)xijt≤Myk,t−1,for all t∈T\{0},i,j∈M,k∈I,(i≤ς(k)≤j)∨(j≤ς(k)≤i)(2.1g)xijt≤Myk,t−1 + Myk+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)usummationdisplayt=0summationdisplayj∈Mjnegationslash=i(xijt−xjit)≤Vi + ǫstki , for all i∈S,u∈T,Vi ≥0(2.1i)162.4. Objective and constraintsusummationdisplayt=0summationdisplayj∈Mjnegationslash=i(xijt−xjit)≥Vi−ǫstki , for all i∈S,u∈T,Vi < 0(2.1j)usummationdisplayt=0summationdisplayk∈Iyki≥u, for all u∈T, (2.1k)xijt≥0, for all i,j∈M,inegationslash= 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 possibleearth move.2.4.1 Objective functionAs in the LP model, the MILP objective function (2.1a) minimizes the totalhauling costs, but this time over all the time steps. Note that the costcoefficients for the block removal variables are set to zero.2.4.2 Volume constraintsThe 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 overall the time steps.2.4.3 Block removal indicatorThe decision variable ykt, that indicates the removal of a block k at time stept, is binary. If we want to find an optimal removal schedule, we need to setsome of these removal indicators to one. Once we decided to remove blockk 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 theblock at time step t, we need to satisfy any potential demand or supply ofearth 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 indicatorconstraints relate the binary variables with the continuous variables.Remark 2.5. Note that once a block k is removed at time step t, there isno need to enforce the monotonicity ykt ≤yk,t+1 in order to represent thestatus of a block in the later time steps. If we want to move earth over the172.4. Objective and constraintsblock at some later time step u > t, we can just set the variables yk,u−1 to 1and the corresponding removal indicator constraint will already be satisfied.In fact, numerical tests showed that enforcing the monotonicity slows downthe solving process.2.4.4 Move feasibility constraintConsider 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 isnot possible. However, if yk,t−1 = 1,then any (reasonable) quantity can be moved. In order to indicate anyreasonable quantity, we use a very large constant M as upper bound forsuch a move. Equation (2.1g) models this relationship.Consider a move from section i to j. If there are blocks k and k+1, suchthat ς(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, beforethe 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 constraintsRemark 2.6. Proposition 1.8 from the LP model does not apply to the blockremoval MILP model in (2.1). An explanation follows from Example 2.7below.Example 2.7. Consider the example with sections and volumes given inTable 2.1. The road section 3 represents a block with demand of one cubicunit of earth. Road section 4 has a supply of one, but is not accessible dueto 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. SinceTable 2.1: Road sections and volumes (bold font represents a block, ≈ anaccess road).Section i ≈1 2 3 4Volume Vi 0 0 -1 1there is nb = 1 block, we have 2 time steps. It is easy to see that the optimalsolution to this problem is to first remove the block in section 3 by borrowing182.4. Objective and constraintsmaterial from the adjacent section 2, followed by filling the newly createddemand in section 2 with the cut of section 4. Hence the solution isx230 = 1x421 = 1y10 = 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. Ablack colored square indicates the presence of a block, a gray colored squareshows a block that was removed in a previous step. The arrows in betweenthe layers represent earth movements. With the given cost model, any other1 2 3 41 2 3 411 2 3 41Figure 2.1: Earthwork scheme with removal schedule (black squares indicateblocks, gray squares indicate removed blocks).solution would be more expensive.If we disregard the time steps and add the quantities for each move overtime, that is xij = xij0 + xij1, for all i,j∈{1,2,3,4}, then we havex42 = 1x23 = 1and the solution moves one cubic unit of earth from section 4 to 2, and fromsection 2 to 3. Hence, in this example, Proposition 1.8 does not apply tothe MILP in (2.1).192.4. Objective and constraintsFrom Remark 2.6, we conclude that a constraint is required to limit theamount of earth that can be stocked or borrowed from a road section. Thiscan 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 asection plus or minus the stockpile tolerance. However, since stockpiling canoccur 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 lowerbound.2.4.6 Block removal enforcementIn 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 wantto remove at least one block. This is guaranteed with constraint (2.1k). Thisconstraint also ensures that no more than nb + 1 time steps are needed.2.4.7 ObservationsProposition 2.8. Let Vs = summationtexti∈S Vi, Vb = summationtexti∈B Vi, and Vw = summationtexti∈W Vi.If Vs ≥0 and Vw−Vs > 0, then the MILP (2.1) has always a solution. IfVs < 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 borrowpit provide sufficient slack. Since for all i,j∈M, t∈T, cijt∈R++∪{+∞}and xijt≥0, we know thatminsummationdisplayt∈Tsummationdisplayi,j∈Minegationslash=jcijtxijt≥0.In Section 4.1, we show how to set the binary variables ykt in order to obtaina feasible point. Furthermore, in Section 3.2, we show that there are a finitenumber of ways to select the binary variables, and therefore a finite numberof feasible points to the MILP (2.1).Altogether, since the MILP (2.1) has a finite number of feasible pointsthat are bounded below, and there exists at least one such point, there existsa solution to the MILP (2.1).Proposition 2.9. The MILP (2.1) contains (|I|+1)(|M|2−|M|) contin-uous, and (|I|+ 1)|I| binary decision variables.202.4. Objective and constraintsProof. 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. Wehave 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 thepractical implementation of the model is smaller than the number given inProposition 2.9. If we multiply the number of continuous variables fromRemark 1.10 with|I|+1 time steps, we get|S|(|M|−1)(|I|+1) continuousvariables. The (|I|+1)|I|binary variables remain the same as in Proposition2.9.In the next section we examine some of the properties of Model (2.1).21Chapter 3Model analysisAut inveniam viam aut faciam.2Hannibal (247–183 BCE)3.1 Space complexityIn 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 pa-rameters, that is required to store the model in the memory of a computer.In this context, we employ the usual asymptotic notation for space complex-ity. The following definition is a reminder of the most common asymptoticbounds 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 that0≤c1g(n)≤f(n)≤c2g(n) for all n≥n0},and we say that g(n) is an asymptotic tight bound for f(n). We denote byO(g(n)) the set of functionsO(g(n)) ={f(n) : there exist c > 0, and n0∈Z++ such that0≤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 modelin (2.1) is of order Θparenleftbig|I||M|2 +|I|2parenrightbig.2I shall either find a way or make one.223.1. Space complexityProof. 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|,and12(|I||M|2 +|I|2)≤|I||M|2 +|I|2 +|I|−|I||M|−|M|2−|M|.We therefore conclude that the growth rate for the size of the total variablesis Θparenleftbig|I||M|2 +|I|2parenrightbig.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 capacityconstraint. Thus there are|B|+|W|capacity constraints. Altogether, thereare|S|+|B|+|W|=|M|balance and capacity constraints. The growth rate of the balance constraintsis therefore of order Θ(|M|).Lemma 3.4. The growth rate of the number of removal indicator constraints(2.1e) and (2.1f) is Θparenleftbig|I|2parenrightbig.Proof. We have one block removal indicator constraint for each indicatorvariable. From Proposition 2.9 we know that the number of indicator vari-ables is |I|(|I|+ 1). Hence the number of indicator variable constraints is|I|(|I|+ 1) which is in Θparenleftbig|I|2parenrightbig.Lemma 3.5. The growth rate of the number of feasible move constraints(2.1g) and (2.1h) is of order Oparenleftbig|M|2|I|2parenrightbig233.1. Space complexityProof. We consider the (|M|2−|M|) potential moves between all the sec-tions over|I|+1 time steps. Each move needs at most a total of|I|feasiblemove constraints. If the move occurs over one or more blocks, we use con-straint (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 nblocks 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 atotal of at most|I|constraints. We have therefore(|M|2−|M|)(|I|+ 1)|I|feasible move constraints, which is of orderOparenleftbig|M|2|I|2parenrightbig.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 have2|T||S|= 2(|I|+ 1)|S|≤2|I||M|+ 2|M|stockpile constraints which is inO(|M||I|).Lemma 3.7. The growth rate of the number of removal enforcement con-straints (2.1k) is of order Θ(|I|).Proof. At each time step, we need a constraint that sets the least numberof 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 Θparenleftbig|M|2|I|parenrightbig.Proof. We need a non-negativity constraint (2.1l) for each of the continuousvariables. From Proposition 2.9, we know that there are (|I|+1)(|M|2−|M|)continuous variables, which is in Θparenleftbig|M|2|I|parenrightbig.Lemma 3.9. The growth rate of the number of binary constraints (2.1m)is of order Θparenleftbig|I|2parenrightbig.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 Θparenleftbig|I|2parenrightbig.243.1. Space complexityProposition 3.10. The growth rate of the memory size for the MILP de-scribed in (2.1) is of order Oparenleftbig|M|7parenrightbig. If the number of blocks |I| is fixed,then the memory size for the MILP (2.1) is of order Oparenleftbig|M|4parenrightbig.Proof. We add the asymptotic bounds for all the constraints from Lemma3.3 to Lemma 3.9. The growth rate of the sum of all the constraints isthereforeΘ(|M|) + Θparenleftbig|I|2parenrightbig+Oparenleftbig|M|2|I|2parenrightbig+O(|M||I|) + Θparenleftbig|M|2|I|parenrightbig+ Θparenleftbig|I|2parenrightbig,(3.1)which is of orderOparenleftbig|M|2|I|2 +|I|2parenrightbig. From Lemma 3.2, we know that thenumber of variables is of order Θparenleftbig|I||M|2 +|I|2parenrightbig. Hence, the size of theproblem is bounded byOparenleftbig|M|2|I|2 +|I|2parenrightbig×Θparenleftbig|I||M|2 +|I|2parenrightbigwhich is of order Oparenleftbig|M|2|I|4 +|M|4|I|3parenrightbig. Since |B|≥ 1, and |W|≥ 1,we know that |I|<|M|, and hence |I|∈O(|M|). The asymptotic upperbound for the memory size is therefore Oparenleftbig|M|7parenrightbig. If we fix the number ofblocks|I|, the dimension of the problem is bounded byOparenleftbig|M|2parenrightbig×Θparenleftbig|M|2parenrightbig,which is of orderOparenleftbig|M|4parenrightbig.Remark 3.11. The asymptotic bound of Oparenleftbig|M|7parenrightbig in Proposition 3.10 is aworst-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 inChapter 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 withorderOparenleftbig|M|4parenrightbig.Regardless of the fact that we fix the number of blocks in practice, thememory growth rate inOparenleftbig|M|4parenrightbig remains a concern for roads with hundredsof sections. Most MILP solvers take advantage of the sparsity of (2.1) andthe bound may not be attained. However, in our numerical tests, we stillobserved a growth rate that quickly exhausted the memory of a standarddesktop computer. From the given bound and from practical experience us-ing a 64 bit desktop computer with 8 GB of RAM and the Gurobi Optimizer2.0 MILP solver, we recommend to keep the number of sections|S|≤200.253.2. Schedule search3.2 Schedule searchIn Section 3.1, we showed that the memory size grows in polynomial timewith the input. However, the difficult part in our model is the schedulingfor the block removal. In this section we analyze the scheduling problemmore closely. In particular, we look at the schedule search space and the setof feasible schedules in the schedule search space.For all time steps t∈T, we list the variables ykt as a binary string in theorder of block position. Since|T|=|I|+1, we have|I|+1 binary strings oflength|I|. At time t = 0, all the blocks are in place, thus yk0 = 0, and thebinary string contains zeros only. At a time step t > 0, the correspondingbinary 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 eachstring 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 newstring which is 10000. We continue the removal until the binary string is11111. Note that complete removal of all blocks may occur before time stept =|I|+ 1. In that case, the digits of the remaining|I|+ 1−t strings canbe considered all 1 for all remaining time steps (the constraints affected bythe binary variables for the remaining time steps are already satisfied).Definition 3.13. The schedule search space is the number of potentialschedules in the MILP (2.1) as defined by the binary selection. This in-cludes infeasible schedules.Lemma 3.14. The schedule search space is asymptotically bounded byOparenleftBig2|I|22|I|parenrightBig.Proof. Since there are|I|+ 1 binary strings of length|I|, we have2|I|(|I|+1)potential schedules. It follows that the schedule search space is inOparenleftBig2|I|22|I|parenrightBig.263.2. Schedule searchTable 3.1: Road sections and volumes (bold font represents a block, ≈ anaccess road).Section i 1 2 3 4 ≈5 6 7Volume Vi 2 -2 -1 -1 5 -1 -2000000010000110001110111111111(a) Feasible000001000010001110011101111111(b) InfeasibleFigure 3.1: Potential schedulesNot all schedules in the schedule search space are feasible. Blocks canonly be removed if they have an access road attached or if an adjacent blockthat is accessible is removed first. We are interested in the number of feasibleschedules that lead to a complete removal.Definition 3.15. The feasible schedules are the schedules that only removea block if the block is accessible at the appropriate time. Hence, if k∈I isa block at time step t∈T, then ykt = 1 implies that either there exists anaccess road r∈R, such thatς(k−1)≤̺(r)≤ς(k)∨ς(k)≤̺(r)≤ς(k + 1),or we haveyk−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 thatjoins section 5. Figure 3.1 shows two potential schedules. The schedule inFigure 3.1a is feasible, with a possible earthwork scheme shown in Figure3.2. However, the schedule in Figure 3.1b is not feasible, since the blocksat either end of the road are not accessible. There are other blocks betweenthe access road and the end of the road.Blocks can be removed sequentially or simultaneously. From an accessroad with blocks on either side, there are three possible ways to remove273.2. Schedule search1 2 3 4 5 6 71 2 3 4 5 6 711 2 3 4 5 6 711 2 3 4 5 6 721 2 3 4 5 6 711 2 3 4 5 6 72Figure 3.2: Earthwork scheme with removal schedule.283.2. Schedule search000001100100 001011111110011011101100 0111 0110 0111 0011111101111111111011111110 1111 1110 1111 0111 1111 01111111 1111 1111 1111 1111 1111Figure 3.3: Tree structure of all feasible block removal schedules for 4 blockswith an access road joining between the two middle blockadjacent blocks. We can remove the block to the left, to the right, or bothat once. If the access road is attached to the block, the block needs to beremoved first. Using edges to represent left, right, or simultaneous removal,we can draw a tree that contains all the feasible schedules for a given con-figuration of access points. An example for 4 blocks with an access roadattached 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 sin-gle access road from below. We therefore restrict the set of feasible schedulesby analyzing schedules that allow left and right removal only, using one sin-gle access road. Let the access road be at a section between block r andblock k =|I|−r. Consider the (k×r)-lattice of integral points in Z2 withcoordinate (0,0) at the lower left corner and coordinate (k,r) at the upperright corner. We associate a left block removal to a step upward and a rightblock removal to a step to the right. We call the vertical steps (1,0)-stepsand the horizontal steps (0,1)-steps to represent the increase of the y−andx−axis respectively. A lattice path consisting of (0,1)- and (1,0)-steps thatstarts at (0,0) and ends at (k,r) is in bijection with a schedule starting atthe root node with a zero binary string and ending at a leaf with a binarystring of ones. Figure 3.4a depicts such a path. The next Lemma providesa formula to count all the possible lattice paths starting at (0,0) and endingat (k,r). We provide the proof as outlined in [Aig07].Lemma 3.17 (Aigner [Aig07, p. 17]). Consider a (k×r)-lattice of integralpoints in Z2. The number of paths consisting of (1,0) and (0,1) steps from293.2. Schedule search(0,0)(6,4)rk(a) (1,0),(0,1)-lattice path(0,0)(6,4)(b) Delannoy pathFigure 3.4: (6×4)-lattice in Z2, representing 10 block configuration with asingle access road between block 6 and 4(0,0) to (k,r) isL(k,r) =parenleftbiggk + rrparenrightbigg. (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}containingexactly r L’s. This corresponds to parenleftbigk+rr parenrightbig possible combinations of (k + 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 haveL(k,r) =parenleftbigg|I|rparenrightbigg(3.3)feasible schedules for the restricted set of feasible schedules for a single accessroad 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 triangle11 11 2 11 3 3 11 4 6 4 1.303.2. Schedule searchLet the rows represent the number of blocks|I|, and the columns representblock r, where a single access road is between block k and r. The binomialcoefficients in the triangle represent the number of feasible left and rightremoval schedules. We notice that the central binomial coeffcients along thecentral axis of the triangle grow the fastest. They correspond to an accessroad between k and r, where k = r. If k = r, with the formula given in[Wei02, p. 369], thenL(r,r) =parenleftbigg2rrparenrightbigg=parenleftbigg|I|rparenrightbigg.Using Stirling’s approximation [Wei02, p. 2868] for the factorial function,we have parenleftbigg|I|rparenrightbigg≤ 4r√3r + 1, for all r≥1.Since|I|= 2r, we conclude that the number of feasible schedules for a singleaccess road and left and right removal only is inOparenleftbig2|I|parenrightbig.We now extend our restricted set of feasible schedules by adding anotherremoval step to the schedules with one access. So far we considered left-removal and right-removal steps only. But if the access road is in betweentwo blocks, we can also remove the left and right block at the same time, asshown in Figure 3.3. Remark 3.19 is therefore a lower bound on the worst-case. Instead of having only (0,1)- and (1,0)-steps in our lattice path, wenow add a diagonal (1,1)-step. This results in Delannoy paths as shown inFigure 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. Thenumber of Delannoy paths consisting of (1,0),(0,1) and (1,1) steps from(0,0) to (k,r) isD(k,r) =ksummationdisplayi=0parenleftbiggkiparenrightbiggparenleftbiggk + r−ikparenrightbigg. (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 astring 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 project3the (1,1)-steps onto the horizontal x-axis. The projected {L,R} string has3The projection of the diagonals onto a lattice is based on an idea by Jim Nastos.313.2. Schedule searchnow r−i L’s and k−i+i = k R’s. Using the formula for (1,0),(0,1)-latticepaths in Lemma 3.17, Equation (3.2), we haveL(k,r−i) =parenleftbiggk + r−ir−iparenrightbigg=parenleftbiggk + r−ikparenrightbigg(3.5)(1,0),(0,1)-paths. However, there are parenleftbigkiparenrightbig possible ways we can choose thediagonals that we project. Also, we have i∈{0,...,k}. Hence we multiplyEquation (3.5) by parenleftbigkiparenrightbig and sum i from 0 to k in order to get the formula inLemma 3.20.Proposition 3.21. The number of feasible schedules for a road with a singleaccess 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 fora given road with one access road. However, a lot of construction siteshave more than one access road. We are particularly interested in the verycommon case of two access roads, one at the start and one at the end of theroad.Lemma 3.22. The number of feasible schedules for a road with two accesspoints, one at the start and one at the end of the road, isS(|I|) =|I|summationdisplayj=1D(|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 stepat time t = 0. If we remove a block at the start of the road, we set the left-most digit from 0 to 1. If we remove a block at the end, we set the right-mostdigit to 1. Likewise, if we simultaneously remove on both sides, we set theouter-most digits to 1. Since there is no access point anywhere else, thisprocedure continues until all digits are set to 1. The last digit or the lasttwo digits that are set to 1 represent the meeting point of the constructioncrews, once all the blocks are removed. For any removal schedule with ameeting point j∈{1,...,|I|}, we can invert the schedule in the sense thatzeros are set to one, and ones to zero. A inverted schedule represents aschedule that has the meeting point as a single access road. For a fixed323.2. Schedule searchmeeting point, there is a bijection between all the schedules with two accessroads that meet at the meeting point and the schedules with a single accessroad attached to the meeting point. In order to get all feasible schedules, weconsider 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 num-bers along the diagonal from (0,0) to any (r,r) grow the fastest. To get aworst-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 approximationto compute an upper bound on the central Delannoy numbers. It is givenasD(k,k) = 5.82842709k(0.57268163k−1/2−0.06724283k−3/2+ 0.00625063k−5/2 +OparenleftBigk−7/2parenrightBig). (3.7)Proposition 3.23. The number of feasible schedules for a road with a singleaccess road between the center blocks of the road is inOparenleftbig3|I|parenrightbig.Proof. Consider Equation (3.4). If the single access road is between thecenter blocks, we have r = k. The dominating term in the sum of Lemma3.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. Since5.83|I|/2(|I|/2)−1/2 lessorsimilar 2.42|I|√2|I|−1/2,we conclude that the number of feasible schedules is inOparenleftbig3|I|parenrightbig.Proposition 3.24. The number of feasible schedules for a road with twoaccess roads, one on each end of the road, is in Oparenleftbig3|I|parenrightbig.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 formP(n) = (1 +√2)n−(1−√2)n2√2 .Since P(n) ≤ (1 +√2)n√2/4 ≈ 2.42n√2/4, we conclude that the numberof feasible schedules is inOparenleftbig3|I|parenrightbig.333.2. Schedule searchRemark 3.25. Consider Proposition 3.21, 3.23, and 3.24. If there is oneaccess road only and the number of blocks |I| is small, then the numberof feasible schedules is clearly tractable. In this case, it is reasonable toenumerate all the feasible schedules and to solve a LP for each of the feasibleschedules. The optimal removal schedule is the schedule for the LP solutionwith the lowest objective value. The enumeration of the feasible schedulescan be done with the use of a Depth-First search algorithm [CLRS09]. TheDepth-First algorithm runs linearly in the size of the tree. The solving timeis 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 anexponential number of feasible schedules. Whenever |I| is large, the enu-meration of the feasible schedules is not tractable and one should solve asingle MILP instance of the form in (2.1) instead.If there are two access roads, one on either side of the road, then thereis an exponential number of feasible schedules. Again, if |I| is large, oneshould solve the MILP in (2.1).34Chapter 4AlgorithmsAn algorithm must be seen to bebelieved.Donald Knuth (1938–)Model (2.1) is sufficient to solve the earth-moving problem. In this sec-tion we propose algorithmic techniques that improve the speed of conver-gence. 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 ofcontinuous variables and n = nc + nb denote the total number of variablesin (2.1). We rewrite the inequality constraints and group the left-hand sidesin 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 formmin cTxs.t. A(x,y)T ≥bx≥0y∈{0,1}nb,(4.1)where x ∈ Rnc, c ∈ Rnc++ is the cost vector, and x ≥ 0 means componentwise 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 algorithmin the common sense, we call it here the direct solve algorithm for futurecomparison.Algorithm 1: Direct solveInput: x,y,c,A,b (from (2.1))Output: Global minimum z∗ = cTx∗ for(4.1)beginBuild and solve MILP (4.1) for x∗;end354.1. Mass schedule as starting point4.1 Mass schedule as starting pointDefinition 4.1. The LP obtained by relaxing the constraints y ∈{0,1}nbwith 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 onthe 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 thecutting-plane algorithm of Gomory [Gom58]. A classical branch-and-boundalgorithm solves first the LP relaxation to the MILP. The algorithm thenpartitions the LP relaxation into subproblems which have one of the pre-viously relaxed integer variables fixed to closest integer lower- and upper-bounds. The continuation of this process results in a branch-and-boundtree with nodes that represent subproblems. The bounding part of the algo-rithm looks for the best complete solution zm in the tree. At each node k,where the solution zk to the subproblem is greater than zm, the algorithmcan discard any further subproblem created by this node, since all followingsolutions will be greater than zm. For a more detailed description of thebranch-and-bound algorithm, we refer to [PS98].A common technique to accelerate the solution process of a MILP is toprovide an initial upper bound z◦ on the solution value zm. This allows forearly bounds on the branch-and-bound tree. An even better technique is toprovide the solver with a feasible starting point (x◦,y◦), where z◦ = cTx◦.The number of nodes that need to be expanded can be reduced significantlywith a starting point (x◦,y◦) that produces a sufficiently low upper boundz◦. In this section, we describe a way to find a reasonable starting point.In a traditional mass diagram procedure, work is performed in a unidi-rectional way from the start of the road to the end (see Remark 1.2). Fora road with access attached to the starting section, a feasible block removalschedule is obtained by removing the blocks sequentially from the start ofthe road to the end, one block at each time step. If there is no access at thestarting section of the road, we select the section with an attached accessthat is closest to the starting section and first remove sequentially all theblocks on the left of the access, and then all the blocks on the right. Oncewe have the schedule, we can fix the binary variables in (4.1) and solve theproblem as a LP. This will give us feasible values for the continuous vari-ables, and together with the fixed binary variables, a feasible solution to(4.1).364.2. Maximal allowed moving distance dmaxProcedure FindFeasiblePoint(x,y,c,A,b)Input: x,y,c,A,b (from (2.1))Output: Feasible point (x∗,y∗) for (4.1)begina←FindFirstAccess(S);for k←a to 1 doy∗k,a−k+1←1;endfor k←a + 1 to|I|doy∗kk ←1;endfor k←1 to|I|dofor t←1 to|I|doif y∗k,t−1 = 1 then y∗kt←1;endendSolve (4.1) with fixed y∗kk as LP;endWe can now extend Algorithm 1 with procedure FindFeasiblePoint. We callthis the direct with x◦ algorithm.Algorithm 2: Direct with x◦Input: x,y,c,A,b (from (2.1))Output: Solution z∗ = cTx∗ for (4.1)beginBuild MILP (4.1);(x◦,y◦)←FindFeasiblePoint(x,y,c,A,b);Solve MILP (4.1) with (x◦,y◦) for (x∗,y∗);end4.2 Maximal allowed moving distance dmaxTo 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 thelargest number of sections over which earth is allowed to be moved. There-374.2. Maximal allowed moving distance dmaxfore, 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 makeslittle sense to move earth from a supply at one end of the road to a demandat the other end of the road. Intuition tells us that if there are other suppliesalong the road in between the two sections, then the closest supply to thedemanding section should be chosen. Hence if we disallow moves that goover too long a distance, then we restrict the search space to shorter, andthus probably cheaper moves. We enforce dmax with a set of constraints inthe form of Equation (4.2). If there are p such constraints, we can groupthem in a matrix Rdmax ∈{0,1}p×n and reformulate the MILP (4.1) in admax-extended formmin cTxs.t. A(x,y)T ≥bRdmaxx = 0x≥0y∈{0,1}nb.(4.3)Remark 4.4. It is clear from (4.2) that if dmax is small, then Rdmax becomeslarge (although very sparse). So instead of reducing the model size, weincrease it. However, with each constraint in Rdmax, we can remove a largenumber of feasible move constraints in A. This would be done before solvingthe model; potentially while building the model. Therefore, the model sizeis actually decreased. In fact, not only is the model size reduced, but alsothe problem becomes easier to solve due to the sparsity of Rdmax.Lemma 4.5. For any optimal solution value z∗ = cTx∗ to the MILP (4.1)and any optimal solution value z∗dmax = cTx∗dmax to the extended form (4.3),the inequalityz∗≤z∗dmaxis 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,Rdmaxx = 0,x≥0,y ∈{0,1}nb}. Clearly, Cdmax ⊆C.It follows that z∗≤z∗dmax.From Lemma 4.5 we can see that the optimality of the solution withrespect to (4.1) is not guaranteed when using dmax, unless dmax ≥|S|. Forthis 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 are384.3. Warmstarting with increasing dmaxapproximations. As such, a solution to (4.3) that is not optimal with respectto the original problem, may still be a good enough solution for an earth-moving problem. We have therefore an approximation algorithm that wecall direct with x◦ and dmax.Algorithm 3: Direct with x◦ and dmaxInput: x,y,c,A,b (from (2.1))Output: Solution z∗ = cTx∗ for (4.3)beginBuild Rd0max using (4.2);(x◦,y◦)←FindFeasiblePoint(x,y,c,A,b,Rd0max);Solve (4.3) for (x0,y0) with starting point(x◦,y◦);endNote that in this algorithm, we modified FindFeasiblePoint in order to takeRd0max into account. The change is implicit and we do not show the adaptedprocedure here.4.3 Warmstarting with increasing dmaxIn Section 4.2, we mention that the dmax-heuristic does not guarantee anoptimal solution for (4.1). However, the solution from (4.3) can be used asa starting point to solve (4.1) using a warmstart procedure. Further im-provement can be made by starting to solve (4.3) with a small dmax and usethe solution to warmstart (4.3) with αdmax where α ∈R++,α > 1. If wecontinue this procedure, we end up with the following algorithm, which werefer to as warmstart algorithm.394.4. Warmstarting with mipgap checkAlgorithm 4: WarmstartInput: x,y,c,A,b,d0max ≥1,α > 1Output: Solution z∗ = cTx∗ for (4.1)begink←0;Build Rd0max using (4.2);(x◦,y◦)←FindFeasiblePoint(x,y,c,A,b,Rd0max);Solve (4.3) for (x0,y0) with starting point(x◦,y◦);while dkmax < ms dok←k + 1;dkmax←αdk−1max;Build Rdkmax using (4.2);Solve (4.3) for (xk,yk) with starting point(xk−1,yk−1);end(x∗,y∗)←(xk,yk)endLemma 4.6. The warmstart method in Algorithm 4 converges to the optimalsolution 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 αkd0max ≥ms. Therefore, atiteration k of Algorithm 4, we havedkmax = αdk−1max = αα···αd0max = αkd0max≥ms.For this k, the extended MILP (4.3) is equal to the original MILP (4.1).4.4 Warmstarting with mipgap checkDefinition 4.7. The mipgap |zk−z∗LP|is the gap between the current bestfeasible value of the MILP, zk, and the optimal value of its LP relaxation,z∗LP.Lemma 4.8. The optimal value z∗LP of the LP relaxation for to the originalMILP (4.1) gives a lower bound for the dmax-extended MILP (4.3).Proof. Combine Fact 4.2 and Lemma 4.5.404.4. Warmstarting with mipgap checkLemma 4.9. Let z∗LP be the optimal value for the LP relaxation of (4.1)and z∗ be the optimal value for the MILP (4.1). For any current best feasiblevalue zkdmax of (4.3),|z∗−z∗LP|≤|zkdmax −z∗LP|.Proof. The result follows from Lemma 4.5 and Lemma 4.8.We can now extend Algorithm 4 with a stopping condition that verifiesthe mipgap. We call this the warmstart with ǫ algorithm.Algorithm 5: Warmstart with ǫInput: x,y,c,A,b,d0max ≥1,α > 1,ǫ > 0Output: Solution z∗ = cTx∗ for (4.1), within ǫbeginSolve for the LP relaxation solution z∗LP of(4.1);k←0;Build Rd0max using (4.2);(x◦,y◦)←FindFeasiblePoint(x,y,c,A,b,Rd0max);Solve (4.3) for z0 = cTx0 using starting point(x◦,y◦);while dkmax < ms and ǫ <|zk−z∗LP|dok←k + 1;dkmax←αdk−1max;Build Rdkmax using (4.2);Solve (4.3) for (xk,yk) with starting point(xk−1,yk−1);end(x∗,y∗)←(xk,yk)endIn the next chapter, we look at the numerical results for Algorithms 1 - 5with a set of problems based on real road alignment designs. We analyzesolving times and quality of solutions for the different algorithms and try tofind good values for the parameters α and dmax.41Chapter 5Numerical resultsUnderstanding is, after all, whatscience is all about-and scienceis a great deal more than meremindless computation.Roger Penrose (1931–)Thereare two categories of questionsthat we seek to statistically examinein this section. First, which algorithm is the fastest, and which algorithmis the most robust to solve our MILP model? Second, we want to knowhow the solution value of our model compares to the solution value of a LPsolution for a given non-optimal removal schedule.5.1 Experimental setup5.1.1 ProblemsIn order to compare the algorithms in Chapter 4, we based our test prob-lems on a family of real road alignments. We used a road design softwareto generate 8 different road alignments with length between 1 and 40 kilo-meters. We denote the alignments with letters from A to H. From the 8designs, we then created 44 problems with different block and access roadconfigurations. For all of the problems we added two access roads, one at thebeginning and one at the end of the road. For about half of all the problemswe added two additional access roads in between the start and the end ofthe road. We further split each group into three subgroups of 5, 10, and 20blocks. The problem configurations are given in Table A.1 of Appendix A.5.1.2 BenchmarkIn order to evaluate the benefits of our model over current practices, wewould need to conduct a field study on a sufficiently large set of real con-struction projects. Similar to our theoretical problem set, such projects425.1. Experimental setupwould need to have obstacles and access roads in different configurations.The effective earthwork costs could then be compared with the solution ofour model. In our research, we did not have access to such data. We there-fore tried to build a reference solution that, to our knowledge, would be asfair as possible to compare with our solution.Reference solutionFor our reference solution we assume that obstacles are removed in a sequen-tial order. We assume that earthwork tasks can occur at the same time asthe obstacle removal process, as long as the sections for these tasks are acces-sible. We therefore use the removal scheme in Procedure FindFeasiblePointas reference block removal schedule. As in Procedure FindFeasiblePoint,we solve the LP for the fixed schedule. Note that we do not use any dmaxrestrictions in the reference solution. We call this reference solution the LPsolution for mass removal schedule.We are aware that block removal schedules vary for different constructionprojects and they also depend heavily on the experience of the geotechni-cal engineer. A mass removal schedule, however, is perfectly valid in forestareas or in ground conditions where the operation of ordinary constructionequipment is not possible. In very difficult ground settings, tractor mountedequipment is used to push large-sized granular material in front of the equip-ment to build a stable construction platform [Fwa06]. This corresponds toa sequential removal of ground obstacles.ParametersFor some of the tests we changed parameters for the algorithms in Chapter4, 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.435.1. Experimental setupNote that except for the warmstart algorithms, ε = ǫ. For the warmstartalgorithms, we useε =braceleftBigg0.05 if dkmax < ms,ǫ if dkmax≥ms.In some problems, the solving time for the LP relaxation z∗LP (see Algorithm5) can be quite long. To reduce this solving time, we change the LP solveroptimality tolerance εLP. The LP solver in the Gurobi Optimizer has adefault LP optimality tolerance of εLP = 10−6. An increase in the LPoptimality tolerance results in a faster solving time for an LP. Now if for themipgap tolerance ǫ > 0, we have|z∗−z∗LP|≤ǫ,then for some δ > 0, by the triangle inequality, we have|z∗−z∗LP + δ|≤ǫ + δ.So if we let the new LP optimality tolerance ˜εLP > εLP, with δ = ˜εLP−εLPand the new mipgap tolerance ˜ǫ = ǫ+δ, then the convergence of Algorithm 5within a mipgap tolerance of ǫ is still guaranteed. For our experiments withthe 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 usedthe default LP optimality tolerance.Test configurationsIn our first set of experiments, our goal is to investigate the behavior ofthe solution and solving time with respect to dmax. Hence we run all theproblems with the direct with x◦ and dmax algorithm (Algorithm 3). Theconfiguration for these tests are given in Table 5.1. In our second set ofexperiments, our goal is to investigate the economical benefit and the per-formance of the algorithms in Chapter 4. The LP solution for mass removalschedule was used for an economical comparison to the exact solution to(2.1), as given by the algorithms in Chapter 4. Furthermore, a performanceanalysis was done. The LP with mass removal schedule was not used inthe performance comparison since it is a pure LP, and it is of little valueto compare the solving time of a single LP with the algorithms in Chapter4. Moreover, the LP with mass removal schedule is not an optimal solu-tion to (2.1). We therefore focus our performance analysis on the MILP445.1. Experimental setupTable 5.1: Test configurations for dmax with the direct solve algorithm usingx◦ (Algorithm 3)ǫ dmax0.01 525100200ǫ dmax0.05 525100200Table 5.2: Economical test and performance test configurationsǫ d0max α Algorithm0.01 - - 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)0.05 - - 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)455.2. Resultsalgorithms only. In the context of economical analysis and performance, werun the experiments with the configurations shown in Table 5.2. For allour tests, we used a CPU time limit of 12 hours. If an algorithm could notsolve the problem within the time limit, we aborted the solving process andconsidered the problem as not solved.Remark 5.1. In practice, a MILP tolerance of 5% is sufficiently precise foran earth-moving problem. The data for earth quantities is based on approx-imations that mostly use triangular surface grids. The use of triangulariza-tion, measurement errors, and non-constant factors for earth shrinkage andswelling, introduce errors that are in the order of 5%.Test equipment and softwareAll the experiments were solved with the Gurobi Optimizer 2.0.2 MILPsolver (http://www.gurobi.com) on a standard 64-bit Linux (Ubuntu 9.10Server) workstation with an IntelR© 2.4GHz Quad-Core CPU and 8GB ofmemory. Since the solver makes use of multiple cores, the reported CPU-time may be greater than the wall-clock time in some problems. The algo-rithms were implemented in C++ and the code can be made available uponrequest. Please contact the author at valentin.koch@gmail.com.5.2 Results5.2.1 Analysis of the dmax-heuristicIn the top row in Figure 5.1 we show the boxplots of the relative error ofthe objective value with respect to the solution for dmax ≥ ms and thecorresponding MILP solver tolerances of ε=1% and ε=5%. The bottomrow contains the relative solving time with respect to the solving time fordmax ≥ ms. In the boxplots, the box represents the interquartile range(IQR) and contains 50% of the sample points. It gives therefore a relativelyrobust statistic. The vertical bar inside the box represents the median. Theends of the whiskers represent the lowest datum still within 1.5 IQR of thelower quartile, and the highest datum still within 1.5 IQR of the upperquartile. Outliers are depicted as small dots. For convenience, the meanis also plotted as a small cross. We only counted problems that could besolved for both, varepsilon=1% and ε=5%. The numbers for the plots aregiven in Table A.2 and A.3 in the Appendix.Looking at dmax = 5 in the 1% configuration, we see that the medianrelative error of 34.84% for the objective value is about 5% smaller than the465.2. Results00.511.525 25 100 msdmax00.511.525 25 100 msRelativeTimedmax-0.200.20.40.60.811.21.45 25 100 ms5% mipgap-0.200.20.40.60.811.21.45 25 100 msRelativeError1% mipgapMeanMeanMeanMeanFigure 5.1: Objective and solving time with increasing dmax and differenttolerance475.2. Resultsmedian relative error of 40.39% for the 5% configuration. Remember thatthe reference solution for the relative error is also solved within a toleranceof 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 beacceptable. 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 canlook at the relative solving time and we can see that for 1% tolerance, atdmax = 5, we need 10.2% of the time that the direct solve with x◦ algorithmuses to solve the MILP. The time for a 5% tolerance at dmax = 5 is slightlysmaller, at 6.4% of the direct solve and x◦, which can be explained withthe 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 muchdifference at dmax = 100 to the direct solve algorithm. Consider, however,that 80% of all problems have ms≤100 and therefore have the same resultswith dmax = 100 and dmax≥ms.If we add the IQR to the solving tolerance, we see that a choice ofdmax = 25 produces a reasonable result that lies usually within 5.6% of thereal solution with 1% solving tolerance, and within 4.3% of the solution with5% solving tolerance. The results with 1% solving tolerance do not lie withinthe 5% suggested by Remark 5.1. However, a solution from dmax = 25 issufficiently close to the real solution in order to give a good upper boundfor a subsequent MILP with a larger dmax or dmax≥ms.The median solving time with dmax = 5 for 1% tolerance is around 10times faster, and with 1% tolerance around 16 times faster than the directwith x◦ algorithm. In more than half of all the problems, the solving timedoes not exceed 16% of the solving time with the direct and x◦ algorithmfor both tolerances. In the dmax = 25 case, the solving time is about 4 and7 times faster. We observe that the upper quartile for 1% is at about halfand 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 warmstartscheme in order to derive the effect of the starting point given by dmax = 5and dmax = 25.5.2.2 PerformanceIn order to compare the performance of the algorithms, we use a performanceprofile as developed in [DM02]. For each algorithm, we plot the number ofproblems that are solved within a factor of the best CPU time. Let P bethe set of all problems, and A be the set of algorithms. We compute the485.2. Resultsperformance ratio byrp,a = tp,amin{tp,a : a∈A},where p∈P, and tp,a is the solving time of algorithm a for problem p. Thepercentage of problems that are solved within a factor τ ∈R are computedbyρa(τ) = |{p∈P: rp,a≤τ}||P| .Hence for a particular problem p, an algorithm a needs τ times longer thanthe best algorithm for p, if τ > 1 when algorithm a is able to solve it. Ifalgorithm a needs only τ = 1, the algorithm is considered as the best, oramong the best algorithms for problem p. The percentage ρa of all problemssolved by algorithm a at τ gives us an overall assessment of the performanceof algorithm a. We conclude that the algorithm with the highest ρ valueis 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 theother algorithms. If we increase τ to allow the algorithms to be within alarger factor of the best algorithm, then the algorithm with the highest ρ isthe most robust algorithm. Hence an algorithm with the highest ρ value ata large τ, solves most of the problems. For a more detailed description ofperformance profiles, we refer the reader to [DM02].Figure 5.2 shows the performance profiles for Algorithms 1 to 5 on a log2scale, that is,ρa(τ) = |{p∈P: log2(rp,a)≤τ}||P| .Figure 5.2a shows the profile for ε=1% and Figure 5.2b shows the profileswith ε=5%. In Figure 5.3a and 5.3b, the profiles are shown on a normalscale with a short τ interval. Figures 5.4a and 5.4b show the profiles on anormal 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 outper-forms the other algorithms with a probability of about 32%, and the proba-bility that the warmstart with ǫ (d0max = 25,α = 10) algorithm outperformsall other algorithms is at about 27%. If we increase τ to let the algorithms bewithin a factor of 4 of the direct solve with x◦ algorithm, we can see that thewarmstart with ǫ (d0max = 25,α = 10) algorithm wins with a high probabilityof 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 robust495.2. Results00.20.40.60.810 1 2 3 4 5 6 7ρa(log 2(τ))log2(τ)direct solvedirect with x◦warmstart (d0max = 25,α = 10)warmstart with ǫ (d0max = 25,α = 10)warmstart with ǫ (d0max = 5,α = 5)(a) 1% MILP tolerance00.20.40.60.810 1 2 3 4 5 6 7ρa(log 2(τ))log2(τ)direct solvedirect with x◦warmstart (d0max = 25,α = 10)warmstart with ǫ (d0max = 25,α = 10)warmstart with ǫ (d0max = 5,α = 5)(b) 5% MILP toleranceFigure 5.2: Algorithms performance profiles with log2 scale505.2. Resultsalgorithm (82% of all problems solved), followed by the warmstart with notmipgap check and the direct solve with x◦ algorithm (79.5% of all problemssolved). We conclude that the warmstart with ǫ (d0max = 25,α = 10) algo-rithm 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 ≤25on 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 errorfor dmax = 25 is 0% with an upper quartile of 4.6%. Therefore half of thesolutions were above 4.6% of the real solution. From Lemma 4.9 we expectthat the warmstart with ǫ algorithm will use at least one more iterationfor these solutions. Furthermore, due to the IQR of 4.6% for the solutionerror, there is also a significant number of solutions within the IQR that areabove the required ε=1% tolerance. Altogether, these problems account forthe 66% where the warmstart with ǫ algorithms were required to incrementdmax from 25 to 125 or 250, depending on α.For ε=5%, the warmstart with ǫ (d0max = 25,α = 10) algorithm dom-inates all other algorithms on 46% of all problems. The algorithm clearlystays ahead of the others until τ ≈ 35. At that point, either the warm-start with ǫ (d0max = 25,α = 10) or the warmstart without mipgap checkwould 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) algorithmperforms always better or equal to the other algorithms when ε=5% lies inthe mipgap check. While the median and lower quartile for the solutionerror at dmax = 25 are, similar to the ε=1% experiments, at 0%, the upperquartile is only at 3.3%. Therefore half of all the problems that lie in theIQR at dmax = 25 are already within ε=5%. From the results, we know thatthe 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 per-formance. It starts better than the simple warmstart with no mipgap checkalgorithm and if τ ≈ 3, it also outperforms the direct with x◦ algorithmwith a probability of over 70%. However, after τ = 6, the warmstart withǫ (d0max = 5,α = 5) algorithm performs only better than the direct solvealgorithm and finished by solving 88.6% of all the problems. The directsolve algorithm solved only 81%.In the next section, we will analyze the quality of the solutions withrespect to this LP with mass removal schedule.515.2. Results00.10.20.30.40.50.60.70.80 2 4 6 8 10 12ρa(τ)τdirect solvedirect with x◦warmstart (d0max = 25,α = 10)warmstart with ǫ (d0max = 25,α = 10)warmstart with ǫ (d0max = 5,α = 5)(a) 1% MILP tolerance00.20.40.60.810 2 4 6 8 10 12ρa(τ)τdirect solvedirect with x◦warmstart (d0max = 25,α = 10)warmstart with ǫ (d0max = 25,α = 10)warmstart with ǫ (d0max = 5,α = 5)(b) 5% MILP toleranceFigure 5.3: Algorithms performance profiles with normal scale and smallintervale for τ525.2. Results00.20.40.60.810 20 40 60 80 100ρa(τ)τdirect solvedirect with x◦warmstart (d0max = 25,α = 10)warmstart with ǫ (d0max = 25,α = 10)warmstart with ǫ (d0max = 5,α = 5)(a) 1% MILP tolerance00.20.40.60.810 20 40 60 80 100ρa(τ)τdirect solvedirect with x◦warmstart (d0max = 25,α = 10)warmstart with ǫ (d0max = 25,α = 10)warmstart with ǫ (d0max = 5,α = 5)(b) 5% MILP toleranceFigure 5.4: Algorithms performance profiles with normal scale and largeinterval for τ535.2. Results5.2.3 Economical analysisIn Section 5.1.2, we defined the LP with mass removal schedule as our ref-erence solution. Our reference solution is based on a LP model and doesnot necessarily represent the real cost of an actual earthwork constructionproject with current techniques.Let z∗p,ref be the solution for problem p using the LP with mass removalschedule. Let z∗p,a be the solution for problem p using algorithm a∈A. Wedefine the gain gp,a bygp,a = 1− z∗p,az∗p,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 ofproblems that have been solved. We also plot a normal curve over thehistogram that fits the data. Note that we observed that in some instanceswhere the problems were not solved, the current solution was still betterthan the solution of the LP with mass removal schedule. It is not clearif one should count these solutions in the histogram plots. We decided tokeep these solutions in our economical analysis on the following grounds.In practice, if a solver is interrupted with an intermediate solution that isbetter than a traditional LP solution, the intermediate solution will mostlikely be used.In Section 5.2.2 we identified the best three performing algorithms as thewarmstart with ǫ (d0max = 25,α = 10) algorithm (Algorithm 5), the directwith 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 focusour economical analysis on the results of the warmstart with ǫ (d0max =25,α = 10) algorithm. Figure 5.5 shows the histogram for the warmstartwith ǫ (d0max = 25,α = 10) algorithm. Histograms for the other two algo-rithms are shown in Appendix B. The statistics for the histograms can befound 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 of4.1%. The same algorithm delivered an average of 3.3% gain if ε=5%. Fromthe plot we see that the negative gains in the ε=5% histogram, that arecaused by the MILP tolerance, have a much greater influence on the averagethan in the ε=1% case, where almost no negative gains are encountered. Themaximum gains are 17.9% in the ε=1% histogram, and 19% in the ε=5%histogram.545.2. Results0510152025-0.2-0.15-0.1-0.05 0 0.050.10.150.25% mipgap0510152025-0.2-0.15-0.1-0.05 0 0.050.10.150.2Problems1% mipgapFigure 5.5: Savings with warmstart with ǫ (d0max = 25,α = 10), comparedto mass block removal scheduleOverall we can observe that the average gains for ε=1% are not signifi-cantly higher than for ε=5%. This is surprising since in general, the MILPsolver stops improving the solution in the ε=5% case much earlier than inthe ε=1% case. From the tables in the appendix we can see that the solvingtime for the ε=1% case is much higher than for ε=5%. The qualitative meritof our model can be assessed when looking at the best gains. For each of thetwo 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 largehighway project, this is a significant amount of gain over the LP with massschedule.We are interested in the problems that produced the biggest gains. Inorder to see if there are correlations between problem types and gains, wedid a stratification analysis with the use of a scatter plot. The scatter plot isshown in Figure 5.6. The scatter plots for the other algorithms can be foundin 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 fiveblock problems, blue circles show ten block problems, and red triangles showtwenty block problems. Each group was internally sorted by the number ofsections and plotted from left to right in ascending order. Hence, problem 1in the scatter plot represents the five block problem with the least numberof sections, whereas problem 17 represents a five block problems with thelargest number of sections. Small points represent problems with two accessroads, large points represent problems with four access roads.In the stratification analysis, we could not identify any visibly strikingpattern. The ε=1% and the ε=5% results have a similar distribution. Onemight observe that in the groups of five block problems, there are less solu-555.2. Results-0.0500.050.10.150.20 10 20 30 40 505% mipgap-0.0500.050.10.150.20 10 20 30 40 501% mipgapFigure 5.6: Stratification of gains for warmstart with ǫ (d0max = 25,α = 10);* are 5, ⊙ are 10, trianglesolid are 20 blocks; small symbols are 2, big symbols are 4access roads; within the block groups, points are sorted from left to right inascending order of section numberstions with significant gains than in the other two groups. For these problems,it seems, the LP with mass removal schedule provides often a good enoughsolution. If we group the results with the highest gains, we can observea majority of problems with four access roads. Intuition tells us that withmore access roads, there is more flexibility and hence a potential for cheapersolutions. We conclude that for a better stratification analysis, experimentswith a larger set of problems would certainly give more insight.56Chapter 6ConclusionFrustra fit per plura quod potestfieri per pauciora.4William of Ockham (1287-1347)6.1 AccomplishmentsWe started the thesis with the description of traditional methods to balanceearthwork and minimize hauling costs in the context of road construction.Two main methods were presented in Chapter 1. First, we explained themass diagram, a simple tool that helps to balance the earthwork and toobtain an estimate on the directions of haul. Second, we showed a moredetailed model of a LP by Mayer and Stark [MS81] that not only balancesthe earthwork, but also minimizes the hauling cost. We then pointed outthat neither method captures the problem of blocks that restrict some of theearth movements to balance earth across the construction site. To addressthis problem, we presented an extension to the model of [MS81] in Chapter2. In Chapter 3, we analyzed the model with respect to space complexity andwe investigated the growth of the number of feasible block removal schedules.To improve the computation time for solving the MILP, we outlined a set ofalgorithms in Chapter 4. In Chapter 5, we presented the numerical resultsof the algorithms, together with an analysis of performance, robustness, andaccuracy.6.2 Research outcomesIn the model analysis, we first identified a potential issue of memory growth.With a fixed number of blocks, we found that the theoretical upper boundon the growth of memory for the MILP model in (2.1) is Oparenleftbig|M|4parenrightbig. In4It is pointless to do with more what can be done with less.576.2. Research outcomesour experiments, we use problems with |M|lessorsimilar 200. However, models with|M|≈200 required relatively large amounts of the available memory. Thesolver failed to solve some of these problems due to memory exhaustion.A second result of our model analysis was the enumeration of the feasibleblock removal schedules. We concluded that for a road with one access roadat either end, there is no need to solve for the MILP, and one can just solvea LP for a fixed removal schedule that removes the blocks sequentially. Ifthe access road is at some other place, we can quickly compute the numberof feasible schedules with equation (3.4). If the number of feasible schedulesis small, we can enumerate the schedules and solve a LP for every schedule.The schedule that produces the smallest objective value is the optimal re-moval schedule. We also showed that the worst case of number of schedulesfor one access road is bounded by Oparenleftbig3|I|parenrightbig. Therefore, if |I| is large, werecommend to solve the MILP instead of a schedule enumeration. The sameobservations 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 effi-ciently with the help of the algorithms from Chapter 4. In the experiments,we demonstrated the strength of the dmax-heuristic, and showed that thereis a significant performance benefit when using the warmstart with ǫ, thesimple warmstart, or the direct with x◦ algorithm, instead of solving theMILP directly. Furthermore, we showed that these three algorithms areconsiderably more robust than a method that solves directly for the MILP.As a last observation in our analysis, we showed the economic benefits ofsolving 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 solutionsthat 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 LPsolution for a fixed mass removal schedule. At a first sight, these benefitsmay not seem very high. However, in the context of a highway project thatinvolves earthwork costs of millions of dollars, 4.1% is a significant figure.In summary, we started by presenting a MILP model that overcomes theproblems of the mass diagram and the LP model by Mayer and Stark. Wepointed out potential issues with the model and efficient ways to solve theMILP. The result of the model presents an engineer with a good estimateon earthwork costs. The model is therefore a good tool to quantify a givenvertical alignment, and therefore provides a valuable tool to find the mosteconomic alignment.As already pointed out by Mayer and Stark [MS81], the objective of anearthwork contractor is to maximize its profit. A contract that is based on586.3. Future researchan 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 thatis not possible due to blocks. The presented MILP enables a contractor toestimate the real earthwork costs of a project, and at the same time providesthe contractor with the cheapest haul scheme in order to maximize its profit.6.3 Future research6.3.1 Warmstart schemes and εLPIn Section 5.2.2, we explained the problems that we encountered when solv-ing the LP relaxation with the default LP optimization tolerance of theGurobi Optimizer. At the writing of this thesis, preliminary results of testswith ε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 differ-ent d0max- and α-values.6.3.2 Vertical alignmentThe above model can easily be extended to the vertical alignment. Basedon a vertical optimization model by Easa [Eas87], Moreb introduced a LPmodel that combines the earthwork model by Mayer and Stark with roadwaygrade selection by using a piecewise linear approximation of the road profile[Mor96]. In [MA04], Moreb and Aljohani extended the model by using aquadratic approximation. The quadratic representation was then furtherextended in [Mor09] with a spline approximation, with improvements byKoch and Lucet [KL10]. The result was an improved and efficient linearmodel that avoids sharp connectivities and scaling problems of piecewiselinear representations.The main difference between the LP model of Mayer and Stark, and thevertical alignment model by Moreb, is the earth volume. In the earthworkLP and in our MILP model, we treat the volume of a section i as a fixedconstant Vi. We recall the balance constraint (1.3b) from the LP assummationdisplayj∈Mjnegationslash=i(xij−xji) = Vi for all i∈S.However, in the vertical alignment LP by Moreb, the volume is a linearapproximation. The volume is computed as the product of the surface area596.3. Future researchAi of section i, with the average height of earth that is cut from the sectionor filled into that section. The cut or fill is represented with LP decisionvariables, namely ui for the height of the cut, and vi for the height of thefill in that section. The balance constraint then becomessummationdisplayj∈Mjnegationslash=i(xij−xji) = Ai(ui−vi) for all i∈S.The ground profile can then be modeled with spline segments. The splinesegments do not need any additional decision variables besides ui and vi.We refer the reader to [Mor09, KL10] for the details. Since the model forvertical alignment needs only two additional continuous decision variablesper section, the MILP model in (2.1) can readily be combined with the LPmodel in [Mor09, KL10].6.3.3 Schedule enumerationIn Remark 3.25, we discussed the use of feasible schedule enumeration inthe case of a small number of blocks. The use of a Depth-first search algo-rithm allows an enumeration in linear time. However, the number of feasibleschedules itself is exponential. A natural question is what a small number ofblocks means in this context. It would be interesting to compare the solvingtime for the enumeration of schedules, together with solving the LP for eachschedule, and the solving time for the MILP. We believe that up to a certainnumber of blocks, the enumeration of all LP solutions would be significantlyfaster than solving the MILP. The decision process between the two meth-ods would be simplified significantly with the identification of a thresholdthat indicates the number of blocks where the MILP will outperform anexhaustive 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 removalschedule in terms of binary strings, one binary string for each time step. Aninteresting observation is that the same binary string can be used in multipleschedules. In fact, after a certain time step, multiple schedules may evenend with the same set of binary strings. We can therefore look at a binarystring as a state in the removal process5. Once two schedules share the samestate, they need to solve the same subproblem to finish the removal process.Dynamic programming is particularly well suited for optimization problemswith overlapping subproblems [Bel52]. It would be interesting to formulate5This observation was made by Jim Nastos.606.3. Future researchthe feasible schedule enumeration as a dynamic program, and to comparethe method with exhaustive enumeration, as well as the MILP method. If athreshold was found for the exhaustive enumeration and the MILP, it wouldbe interesting to see if the threshold changes with dynamic programming.6.3.4 Lagrangian relaxationThe final topic for future research is the use of Lagrangian relaxation [Fis04,Lem07]. Let x∈Rn. Consider the MILPz = min cTxs.t. Ax = bDx≤ex≥0 and integral,where b∈Rm,e∈Rk,A∈Rm×n, and D∈Rk×n. The Lagrangian problemof the above MILP iszD(λ) = min cTx + λ(Ax−b)s.t. Dx≤ex≥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. Fromdata of preliminary experiments, we believe that the MILP model could besolved much faster if some of the constraints are relaxed by moving theminto the objective function. For example, we could readily relax the balanceconstraints (2.1b),(2.1c),(2.1d), and the stockpile constraints (2.1i),(2.1j),without affecting the binary variables. The relaxation of the feasible moveconstraints (2.1g) and (2.1h) would greatly reduce the problem size.61Bibliography[Aig07] Martin Aigner. A Course in Enumeration, volume 238 of Grad-uate Texts in Mathematics. Springer-Verlag, 2007. → pages 29,30[Bel52] Richard E Bellman. The theory of dynamic programming. Pro-ceedings of the National Academy of Sciences of the United Statesof 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 ofhorizontal and vertical alignments for highways. TransportationResearch Part B: Methodological, 23(5):315–329, 1989. → pages2, 4[CL06] Juey-Fu Cheng and Yusin Lee. Model for three-dimensionalhighway alignment. Journal of Transportation Engineering,132(12):913–920, 2006. →pages 4[CLRS09] Thomas H Cormen, Charles E Leiserson, Ronald L Rivest, andClifford Stein. Introduction to Algorithms. The MIT Press, 3edition, Jan 2009. →pages 22, 34[Dan98] George B Dantzig. Linear Programming and Extensions. Prince-ton Landmarks in Mathematics. Princeton University Press,1998. →pages 36[DM02] Elizabeth D Dolan and Jorge J Mor´e. Benchmarking optimizationsoftware with performance profiles. Mathematical Programming,91(2):201–213, 2002. →pages 48, 49[Eas87] Said M Easa. Earthwork allocations with nonconstant unitcosts. Journal of Construction Engineering and Management,113(1):34–50, March 1987. →pages 7, 9, 5962Chapter 6. Bibliography[Eas88a] Said M Easa. Earthwork allocations with linear unit costs. Jour-nal of Construction Engineering and Management, 114(4):641–655, December 1988. →pages 7, 9[Eas88b] Said M Easa. Selection of roadway grades that minimize earth-work 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 integerprogramming problems. Management science, 50(12):1861–1871,December 2004. →pages 61[Fou05] R Fourer. Linear programming - software survey. OperationsResearch / Management Science Today, (3):46, 2005. → pages35, 36[Fwa06] T F. Fwa. The Handbook of Highway Engineering. CRC Press, 1edition, September 2006. →pages 43[GH08] Nicholas J. Garber and Lester A. Hoel. Traffic and HighwayEngineering. CL-Engineering, 4 edition, June 2008. →pages 4[Gom58] Ralph E Gomory. Outline of an algorithm for integer solutions tolinear programs. Bulletin of the American Mathematical Society,64:275–278, 1958. →pages 36[JH90] A Jayawardane and F Harris. Further development of integerprogramming in earthwork optimization. Journal of ConstructionEngineering and Management, 116(1):18–34, Mar 1990. →pages7[JK06] Manoj Kumar Jha and Euncheol Kim. Highway route optimiza-tion based on accessibility, proximity, and land-use changes. Jour-nal of Transportation Engineering, 132(5):435–439, May 2006. →pages 3[JS03] J Jong and P Schonfeld. An evolutionary model for simultane-ously optimizing three-dimensional highway alignments. Trans-portation Research Part B: Methodological, 37(2):107–128, Febru-ary 2003. →pages 4[JSJ06] Manoj Kumar Jha, Paul Schonfeld, and J-C Jong. IntelligentRoad Design, volume 19. WIT Press, 2006. →pages 2, 363Chapter 6. Bibliography[KL10] Valentin R Koch and Yves Lucet. A note on: Spline techniquefor modeling roadway profile to minimize earthwork cost. Journalof Industrial and Management Optimization, 6(2):393–400, May2010. →pages 4, 59, 60[Lay99] M G Lay. Ways of the World: A History of the World’s Roadsand of the Vehicles that Used Them. Rutgers University Press,1999. →pages 1[Lem07] Claude Lemar´echal. The omnipresence of Lagrange. Annals ofOperations Research, 153(1):9–27, September 2007. →pages 61[LTL09] Yusin Lee, You-Ren Tsou, and Hsiao-Liang Liu. Optimizationmethod for highway horizontal alignment design. Journal ofTransportation Engineering, 135(4):217–224, Jan 2009. →pages3[MA04] Ahmad A Moreb and M Aljohani. Quadratic representationfor roadway profile that minimizes earthwork cost. Journal ofSystems Science and Systems Engineering, 13(2):245–252, April2004. →pages 59[Mak99] K Makanae. An application of parametric curves to highwayalignment. Journal of Civil Engineering Information ProcessingSystem, 8:231–238, Jan 1999. →pages 4[McA24] John Loudon McAdam. Remarks on the Present System of RoadMaking. London, Longman, Hurst, Rees, Orme, Brown, andGreen, 1824. →pages 1[MKV02] S M¨aki, R Kalliola, and K Vuorinen. Road construction in thePeruvian Amazon: process, causes and consequences. Environ-mental Conservation, 28(03):199–214, 2002. →pages 2[Mor96] Ahmad A Moreb. Linear programming model for finding optimalroadway grades that minimize earthwork cost. European Journalof Operational Research, 93(1):148–154, August 1996. → pages4, 59[Mor09] Ahmad A Moreb. Spline technique for modeling roadway profileto minimize earthwork cost. Journal of Inustrual and Manage-ment Optimization, 5(2):275–283, 2009. →pages 4, 59, 6064[MS81] R Mayer and R Stark. Earthmoving logistics. Journal of theConstruction Division, 107(2):297–312, Jun 1981. →pages ii, 6,7, 11, 13, 57, 58[PS98] Christos H Papadimitriou and Kenneth Steiglitz. CombinatorialOptimization: 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 delan-noy numbers: 10894. The American Mathematical Monthly,110(5):443–444, May 2003. →pages 31[Stu09] Les Stuart. Transportation in Canada 2008 - an overview. Trans-port 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. →pages13[Tri87] D Trietsch. A family of methods for preliminary highway align-ment. 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, 3365Appendix ATablesEach design in Table A.1 is used with 2 access roads, and with 4 accessroads for the same block configuration. An exception are the designs C andG 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 clusteredtogether in adjacent sections.The data for the problems is available upon request. Please contact theauthor at valentin.koch@gmail.com.Table A.1: Problem configurationsProblem Design Length (km) Sections Access roads Blocks1,25 A 4 101 2,4 52,26 B 40 201 2,4 53,27 C 1 51 2,4 54,28 D 1 100 2,4 55,29 E 1 52 2,4 56,30 F 1 103 2,4 57,31 G 9 100 2,4 58,32 H 9 201 2,4 59,33 A 4 101 2,4 1010,34 B 40 201 2,4 1011,35 C 1 51 2,4 1012,36 D 1 100 2,4 1013,37 E 1 52 2,4 1014,38 F 1 103 2,4 1015,39 G 9 100 2,4 1016,40 H 9 201 2,4 1017,41 A 4 101 2,4 2018 C 1 51 2 2019,42 D 1 100 2,4 2020,43 E 1 52 2,4 2021,44 F 1 103 2,4 2022 G 9 100 2 2023 H 9 201 2 524 F 1 103 2 1066Appendix A. TablesTable A.2: Statistics for objective function error with dmaxǫ dmax = 5 dmax = 25 dmax = 100 dmax = ms0.01 Minimum 0.03344 -0.00907 -0.00351 0.00000Lower quartile 0.21869 0.00000 0.00000 0.00000Median 0.34841 0.00000 0.00000 0.00000Upper quartile 0.60876 0.04625 0.00000 0.00000Maximum 1.08125 0.18926 0.00000 0.00000Lower conf. for med. 0.23842 0.00000 0.00000 0.00000Upper conf. for med. 0.45841 0.01304 0.00000 0.000000.05 Minimum -0.01434 -0.02768 -0.04701 0.00000Lower quartile 0.22466 0.00000 0.00000 0.00000Median 0.40394 0.00000 0.00000 0.00000Upper quartile 0.60948 0.03376 0.00000 0.00000Maximum 1.08494 0.20007 0.01690 0.00000Lower conf. for med. 0.30959 0.00000 0.00000 0.00000Upper conf. for med. 0.49830 0.00828 0.00000 0.00000Table A.3: Statistics for solving time with dmaxǫ dmax = 5 dmax = 25 dmax = 100 dmax = ms0.01 Minimum 0.00000 0.00016 0.00051 1.00000Lower quartile 0.03820 0.11630 0.59407 1.00000Median 0.10172 0.26389 0.99488 1.00000Upper quartile 0.16049 0.54474 1.00000 1.00000Maximum 0.44444 1.42857 1.50000 1.00000Lower conf. for med. 0.06544 0.13677 0.87444 1.00000Upper conf. for med. 0.13801 0.39101 1.00000 1.000000.05 Minimum 0.00000 0.00968 0.21221 1.00000Lower quartile 0.03800 0.09709 0.86429 1.00000Median 0.06395 0.13793 1.00000 1.00000Upper quartile 0.15395 0.28175 1.00000 1.00000Maximum 0.33205 1.27652 1.74622 1.00000Lower conf. for med. 0.03800 0.09709 0.96588 1.00000Upper conf. for med. 0.09310 0.18435 1.00000 1.0000067Appendix A. TablesTable A.4: Histogram statistics for 1% tolerancedirect direct x◦ warmstart warmstart ǫ,25 warmstart ǫ,5Solved 30 35 35 36 33Min -0.00674 0.00000 -0.00674 -0.00198 -0.006741st qrt. 0.00041 0.00000 0.00000 0.00000 0.00000Median 0.02386 0.01878 0.01878 0.02232 0.016963rd qrt. 0.05308 0.04827 0.05051 0.05079 0.04658Max 0.17902 0.17902 0.18987 0.17902 0.15804Mean 0.04341 0.03950 0.04122 0.04088 0.03265Std. dev 0.05429 0.05132 0.05519 0.05154 0.04227Skew 1.13166 1.29910 1.35167 1.24205 1.34292Table A.5: Histogram statistics for 5% tolerancedirect direct x◦ warmstart warmstart ǫ,25 warmstart ǫ,5Solved 36 42 43 43 39Min -0.01722 0.00000 -0.02907 -0.02907 -0.042511st qrt. 0.00000 0.00000 0.00000 0.00000 -0.00115Median 0.01537 0.00000 0.00000 0.00000 0.006153rd qrt. 0.04733 0.04055 0.04290 0.04290 0.03903Max 0.17873 0.17815 0.18987 0.18987 0.15804Mean 0.03865 0.03176 0.03309 0.03309 0.02331Std. dev 0.05471 0.05305 0.05527 0.05526 0.04658Skew 1.18703 1.46810 1.35854 1.35773 1.1779468AppendixA.TablesTable A.6: CPU time in seconds for ε = 0.01, problems 1-22Problem direct suc- direct x◦ suc- warmstart suc- warmstart ǫ suc- dkmax warmstart ǫ suc- dkmaxcess cess cess (25,10) cess (5,5) cess1 79 1 77 1 82 1 73 1 250 68 1 1252 911 1 1199 1 986 1 2622 1 250 3249 1 6253 21 1 16 1 30 1 38 1 250 37 1 1254 26 1 18 1 23 1 9 1 25 11 1 255 14 1 15 1 15 1 12 1 25 14 1 256 55 1 21 1 24 1 55 1 250 29 1 1257 14 1 14 1 18 1 23 1 250 26 1 1258 316 1 316 1 332 1 357 1 250 285 1 1259 7104 1 1603 1 1507 1 1380 1 250 1171 1 12510 43200 0 22338 1 43200 0 17274 1 250 18214 1 62511 603 1 412 1 793 1 816 1 250 1472 1 12512 2936 1 7522 1 1528 1 28593 1 250 3473 1 12513 606 1 134 1 100 1 16 1 25 49 1 2514 107860 0 43256 0 468 1 540 1 250 43231 0 12515 376 1 146 1 206 1 285 1 250 223 1 12516 15777 0 4864 1 5533 1 5629 1 250 11257 1 62517 43200 0 84261 0 43200 0 87130 0 250 43337 0 518 43200 0 43200 0 43200 0 43200 0 25 43200 0 12519 43324 0 43335 0 44453 0 43224 0 25 43176 0 520 43200 0 43200 0 43200 0 43200 0 25 43200 0 12521 28 1 22 1 49 1 45 1 250 56 1 12522 490 1 692 1 506 1 549 1 250 576 1 62569AppendixA.TablesTable A.7: CPU time in seconds for ε = 0.01, problems 23-44Problem direct suc- direct x◦ suc- warmstart suc- warmstart ǫ suc- dkmax warmstart ǫ suc- dkmaxcess cess cess (25,10) cess (5,5) cess23 2 1 2 1 4 1 5 1 250 2 1 2524 184 1 17 1 22 1 9 1 25 28 1 12525 3 1 3 1 4 1 2 1 25 3 1 2526 24 1 25 1 28 1 9 1 25 12 1 2527 17 1 16 1 18 1 8 1 25 10 1 2528 302 1 298 1 315 1 337 1 250 537 1 62529 8531 1 15307 1 5369 1 10121 1 250 9506 1 12530 43200 0 43200 0 43200 0 43200 0 25 43200 0 62531 33 1 20 1 36 1 15 1 25 16 1 2532 328 1 188 1 207 1 39 1 25 403 1 12533 91 1 31 1 40 1 35 1 250 18 1 2534 31175 1 22377 1 22458 1 8180 1 250 11565 1 12535 197 1 182 1 196 1 38 1 25 54 1 2536 43200 0 6486 1 7245 1 7362 1 250 3874 1 12537 43200 0 43200 0 1536 1 7163 1 250 6891 1 12538 2200 1 1425 1 2221 1 1430 1 250 43256 0 539 43200 0 9448 1 9628 1 11723 1 250 43342 0 2540 43327 0 12576 1 43259 0 43202 0 250 43176 0 541 43200 0 44587 0 48538 0 49588 0 250 43200 0 12542 56509 0 43208 0 43207 0 43200 0 25 43208 0 12543 319 1 270 1 383 1 353 1 250 586 1 62544 268 1 173 1 215 1 221 1 250 222 1 12570AppendixA.TablesTable A.8: CPU time in seconds for ε = 0.05, problems 1-22Problem direct suc- direct x◦ suc- warmstart suc- warmstart ǫ suc- dkmax warmstart ǫ suc- dkmaxcess cess cess (25,10) cess (5,5) cess1 59 1 19 1 23 1 26 1 250 40 1 1252 886 1 601 1 1013 1 1044 1 250 107 1 253 20 1 14 1 29 1 30 1 250 38 1 1254 27 1 18 1 23 1 9 1 25 15 1 255 14 1 15 1 13 1 13 1 25 15 1 256 56 1 19 1 23 1 8 1 25 12 1 257 14 1 14 1 18 1 8 1 25 14 1 258 323 1 322 1 338 1 367 1 250 386 1 1259 1232 1 182 1 150 1 163 1 250 242 1 12510 43200 0 2414 1 3356 1 1051 1 25 885 1 2511 609 1 411 1 801 1 809 1 250 1559 1 12512 235 1 145 1 155 1 30 1 25 53 1 2513 596 1 134 1 101 1 88 1 25 61 1 2514 1862 1 140 1 155 1 31 1 25 840 1 12515 383 1 150 1 169 1 33 1 25 62 1 2516 43200 0 4972 1 5635 1 5755 1 250 8530 1 25017 43200 0 5967 1 1416 1 1738 1 250 43248 0 518 43200 0 43200 0 18899 1 19198 1 25 6564 1 2519 33900 1 43338 0 43318 0 43308 0 250 42989 0 520 43200 0 1187 1 1295 1 636 1 25 3106 1 2521 28 1 22 1 27 1 32 1 250 39 1 12522 182 1 181 1 189 1 30 1 25 69 1 2571AppendixA.TablesTable A.9: CPU time in seconds for ε = 0.05, problems 23-44Problem direct suc- direct x◦ suc- warmstart suc- warmstart ǫ suc- dkmax warmstart ǫ suc- dkmaxcess cess cess (25,10) cess (5,5) cess23 2 1 2 1 3 1 2 1 25 3 1 2524 154 1 18 1 21 1 8 1 25 12 1 2525 3 1 3 1 4 1 2 1 25 3 1 2526 25 1 26 1 28 1 8 1 25 14 1 2527 17 1 16 1 19 1 8 1 25 12 1 2528 309 1 305 1 323 1 343 1 250 686 1 62529 186 1 169 1 204 1 228 1 250 351 1 12530 5154 1 5117 1 5169 1 151 1 25 319 1 2531 33 1 20 1 27 1 15 1 25 18 1 2532 189 1 192 1 212 1 37 1 25 55 1 2533 91 1 31 1 37 1 24 1 25 22 1 2534 1536 1 198 1 214 1 41 1 25 60 1 2535 200 1 187 1 200 1 37 1 25 64 1 2536 14951 1 6628 1 7412 1 7513 1 250 5244 1 12537 9991 1 1244 1 1418 1 1533 1 250 5606 1 12538 1846 1 330 1 307 1 302 1 25 43051 0 539 43200 0 1130 1 1289 1 2144 1 250 43200 0 2540 8412 1 1859 1 2333 1 3782 1 250 42900 0 541 43200 0 1732 1 5447 1 409 1 25 5617 1 2542 57882 0 6651 1 7637 1 36483 1 25 38204 1 2543 327 1 276 1 364 1 392 1 250 276 1 12544 273 1 167 1 211 1 64 1 25 46 1 2572Appendix A. TablesTable A.10: Objective value in $ for direct solve, ε = 0.01Problem Objective success ε LP mass value1 454426 1 0.00680 4551812 9351412 1 0.00549 99740943 316359 1 0.00000 3290404 489574 1 0.00446 5052975 60684 1 0.00833 688986 209271 1 0.00540 2110577 77114005 1 0.00198 769612638 56893108 1 0.00714 569262799 455163 1 0.00338 47289110 10000000 0 0.08160 1020450811 376427 1 0.00123 39481912 493254 1 0.00982 51919313 60687 1 0.00396 7207914 212554 0 0.01895 21255415 77114005 1 0.00198 7696126316 0 0 1.00000 5697548517 0 0 1.00000 47730918 0 0 1.00000 51637819 70262 0 0.02304 7959720 0 0 1.00000 21781221 451922 1 0.00983 44889722 9351412 1 0.00966 951275823 297389 1 -0.00000 30040224 482635 1 0.00000 50529725 60253 1 0.00124 6732226 209129 1 0.00559 20924427 61940820 1 0.00000 6194082028 55096639 1 0.00737 5509663929 452486 1 0.00905 46836730 9480000 0 0.04150 962026331 297389 1 0.00000 33904732 484506 1 0.00539 48267133 60686 1 0.00395 6839134 210088 1 0.01000 21150235 61940820 1 0.00000 6194082036 56100000 0 0.02480 5513287137 454000 0 0.01480 47534238 300248 1 0.00952 36015439 0 0 1.00000 48804240 63540 0 0.04868 7739541 0 0 1.00000 21751942 0 0 1.00000 6471681343 56487153 1 0.00000 5648715344 208258 1 0.00073 21418873Appendix A. TablesTable A.11: Objective value in $ for direct solve with x◦, ε = 0.01Problem Objective success ε LP mass value1 454426 1 0.00861 4551812 9351412 1 0.00366 99740943 316359 1 0.00205 3290404 488986 1 0.00327 5052975 60684 1 0.00813 688986 209271 1 0.00627 2110577 76961263 1 0.00000 769612638 56926279 1 0.00771 569262799 455163 1 0.00584 47289110 9498097 1 0.00866 1020450811 376427 1 0.00591 39481912 493254 1 0.00995 51919313 60687 1 0.00396 7207914 212554 0 0.01170 21255415 76961263 1 0.00000 7696126316 56975485 1 0.00857 5697548517 462423 0 0.03118 47730918 495000 0 0.01470 51637819 70262 0 0.10922 7959720 218000 0 0.04520 21781221 448897 1 0.00316 44889722 9351412 1 0.00000 951275823 297389 1 0.00000 30040224 482635 1 0.00092 50529725 60253 1 0.00124 6732226 209244 1 0.00614 20924427 61940820 1 0.00000 6194082028 55096639 1 0.00737 5509663929 452486 1 0.00933 46836730 9430000 0 0.03620 962026331 297389 1 0.00000 33904732 482671 1 0.00160 48267133 60687 1 0.00396 6839134 210088 1 0.00993 21150235 61940820 1 0.00000 6194082036 55132871 1 0.00802 5513287137 455000 0 0.01560 47534238 300248 1 0.00952 36015439 485640 1 0.00771 48804240 63540 1 0.00821 7739541 217519 0 0.04395 21751942 64716813 0 0.04289 6471681343 56487153 1 0.00000 5648715344 210165 1 0.00975 21418874Appendix A. TablesTable A.12: Objective value in $ for warmstart (d0max = 25,α = 10), ε = 0.01Problem Objective success ε LP mass value1 454426 1 0.00843 4551812 9351412 1 0.00000 99740943 316359 1 0.00000 3290404 489574 1 0.00446 5052975 60684 1 0.00833 688986 209271 1 0.00627 2110577 76961263 1 0.00000 769612638 56982170 1 0.00869 569262799 455163 1 0.00452 47289110 9500000 0 0.03160 1020450811 376427 1 0.00346 39481912 493254 1 0.00982 51919313 60687 1 0.00396 7207914 212554 1 0.00996 21255415 77114005 1 0.00198 7696126316 56975485 1 0.00857 5697548517 459000 0 0.02130 47730918 496000 0 0.01730 51637819 64484 0 0.05960 7959720 218000 0 0.04510 21781221 451922 1 0.00983 44889722 9351412 1 0.00984 951275823 297389 1 0.00000 30040224 483187 1 0.00206 50529725 60253 1 0.00124 6732226 209244 1 0.00614 20924427 61940820 1 0.00000 6194082028 55096639 1 0.00737 5509663929 452486 1 0.00999 46836730 9430000 0 0.03620 962026331 297389 1 0.00000 33904732 482671 1 0.00160 48267133 60687 1 0.00396 6839134 210088 1 0.00993 21150235 61940820 1 0.00000 6194082036 55132871 1 0.00802 5513287137 451067 1 0.00795 47534238 300248 1 0.00952 36015439 485640 1 0.00771 48804240 63540 0 0.04868 7739541 217519 0 0.04395 21751942 64716813 0 0.04289 6471681343 56487153 1 0.00000 5648715344 210165 1 0.00975 21418875Appendix A. TablesTable A.13: Objective value in $ for warmstart with ǫ (d0max = 25,α = 10),ε = 0.01Problem Objective success ε ǫ LP mass value1 454426 1 0.00929 0.01434 4551812 9351412 1 0.00591 0.02938 99740943 316359 1 -0.00000 0.06244 3290404 488986 1 0.00327 0.00328 5052975 60684 1 0.00027 0.00840 688986 209271 1 0.00542 0.00631 2110577 76961263 1 0.00000 0.00000 769612638 56982170 1 0.00869 0.00876 569262799 455163 1 0.00338 0.01598 47289110 9498097 1 0.00859 0.04553 1020450811 376427 1 0.00074 0.26417 39481912 493254 1 0.00989 0.01203 51919313 60687 1 0.00396 0.00398 7207914 212554 1 0.00996 0.02209 21255415 77114005 1 0.00198 0.00198 7696126316 56975485 1 0.00857 0.00857 5697548517 458003 0 0.02183 0.02183 47730918 496000 0 0.01730 0.01730 51637819 70234 0 0.10560 0.33415 7959720 218000 0 0.04510 0.04510 21781221 448897 1 0.00316 0.00317 44889722 9351412 1 0.00966 0.02938 951275823 297389 1 0.00000 0.00000 30040224 482635 1 0.00092 0.00092 50529725 60253 1 0.00124 0.00124 6732226 209244 1 0.00614 0.00618 20924427 61940820 1 0.00000 0.00000 6194082028 55096639 1 0.00737 0.00742 5509663929 452486 1 0.00926 0.01119 46836730 9480000 0 0.04150 0.04150 962026331 297389 1 0.00000 0.00000 33904732 482671 1 0.00160 0.00161 48267133 60748 1 0.00496 0.00498 6839134 210088 1 0.00996 0.01024 21150235 61940820 1 0.00000 0.00000 6194082036 55132871 1 0.00802 0.00809 5513287137 451067 1 0.00795 0.00802 47534238 300320 1 0.00976 0.00986 36015439 484188 1 0.00473 0.00473 48804240 63540 0 0.04868 0.05154 7739541 217519 0 0.04395 0.04858 21751942 0 0 1.00000 1.00000 6471681343 56487153 1 0.00000 0.00000 5648715344 208258 1 0.00068 0.00106 21418876Appendix A. TablesTable A.14: Objective value in $ for warmstart with ǫ (d0max = 5,α = 5),ε = 0.01Problem Objective success ε ǫ LP mass value1 454426 1 0.00843 0.01434 4551812 9351412 1 0.00593 0.02938 99740943 316359 1 0.00000 0.06244 3290404 488986 1 0.00327 0.00328 5052975 60684 1 0.00833 0.00840 688986 209271 1 0.00627 0.00631 2110577 76961263 1 0.00000 0.00000 769612638 56576215 1 0.00157 0.00158 569262799 455163 1 0.00600 0.01598 47289110 9498097 1 0.00952 0.04553 1020450811 376427 1 0.00154 0.26417 39481912 493254 1 0.00865 0.01203 51919313 60687 1 0.00396 0.00398 7207914 212554 0 0.01436 0.05611 21255415 77114005 1 0.00198 0.00198 7696126316 56902685 1 0.00730 0.00730 5697548517 656851 0 0.05385 0.52950 47730918 496000 0 0.01710 0.01761 51637819 82236 0 0.06211 0.55954 7959720 218000 0 0.04550 0.04771 21781221 451922 1 0.00983 0.00993 44889722 9351412 1 0.00000 0.02938 951275823 297389 1 -0.00000 0.00000 30040224 483187 1 0.00206 0.00207 50529725 60253 1 0.00124 0.00124 6732226 208491 1 0.00255 0.00256 20924427 61940820 1 0.00000 0.00000 6194082028 54767023 1 0.00139 0.00140 5509663929 452486 1 0.00905 0.01119 46836730 9430000 0 0.03620 0.03755 962026331 297389 1 0.00000 0.00000 33904732 483954 1 0.00425 0.00427 48267133 60735 1 0.00474 0.00477 6839134 210088 1 0.00974 0.01024 21150235 61940820 1 0.00000 0.00000 6194082036 55231802 1 0.00980 0.00989 5513287137 451067 1 0.00795 0.00802 47534238 403569 0 0.06654 0.60979 36015439 538301 0 0.10452 1.00000 48804240 77847 0 0.07685 0.45749 7739541 218000 0 0.04510 0.04801 21751942 64716813 0 1.00000 0.04482 6471681343 56487153 1 -0.00000 0.00000 5648715344 210165 1 0.00975 0.01023 21418877Appendix A. TablesTable A.15: Objective value in $ for direct solve, ε = 0.05Problem Objective success ε LP mass value1 457354 1 0.01905 4551812 9351412 1 0.00549 99740943 316359 1 0.04771 3290404 489574 1 0.00446 5052975 60684 1 0.00833 688986 209271 1 0.00540 2110577 77114005 1 0.00198 769612638 56893108 1 0.00714 569262799 458942 1 0.02107 47289110 10000000 0 0.08160 1020450811 376427 1 0.02126 39481912 494620 1 0.01462 51919313 61213 1 0.01252 7207914 212554 1 0.01895 21255415 77114005 1 0.00198 7696126316 0 0 1.00000 5697548517 0 0 1.00000 47730918 0 0 1.00000 51637819 70262 1 0.02951 7959720 0 0 1.00000 21781221 451922 1 0.00983 44889722 9351412 1 0.02854 951275823 297389 1 -0.00000 30040224 503968 1 0.04316 50529725 60253 1 0.00124 6732226 209129 1 0.00559 20924427 61940820 1 0.00000 6194082028 55096639 1 0.00737 5509663929 453208 1 0.01264 46836730 9556950 1 0.04943 962026331 297389 1 0.00000 33904732 487393 1 0.01128 48267133 60686 1 0.00395 6839134 210088 1 0.01013 21150235 61940820 1 0.00000 6194082036 56082250 1 0.02481 5513287137 468035 1 0.04392 47534238 301565 1 0.01385 36015439 0 0 1.00000 48804240 63562 1 0.04901 7739541 0 0 1.00000 21751942 0 0 1.00000 6471681343 56487153 1 0.00000 5648715344 208258 1 0.00073 21418878Appendix A. TablesTable A.16: Objective value in $ for direct solve with x◦, ε = 0.05Problem Objective success ε LP mass value1 455181 1 0.01577 4551812 9351412 1 0.02854 99740943 316359 1 0.04304 3290404 505297 1 0.03544 5052975 60684 1 0.00813 688986 211057 1 0.01468 2110577 76961263 1 0.00000 769612638 56926279 1 0.00771 569262799 466699 1 0.03908 47289110 9498097 1 0.04354 1020450811 376427 1 0.03718 39481912 500476 1 0.02615 51919313 60687 1 0.00396 7207914 212554 1 0.02161 21255415 76961263 1 0.00000 7696126316 56975485 1 0.00857 5697548517 466407 1 0.03946 47730918 516000 0 0.05600 51637819 70262 0 0.10922 7959720 217812 1 0.04523 21781221 448897 1 0.00316 44889722 9512758 1 0.04502 951275823 300402 1 0.01003 30040224 505297 1 0.04573 50529725 60253 1 0.00124 6732226 209244 1 0.00614 20924427 61940820 1 0.00000 6194082028 55096639 1 0.00737 5509663929 468367 1 0.04460 46836730 9556950 1 0.04943 962026331 297389 1 0.00000 33904732 482671 1 0.00160 48267133 60687 1 0.00396 6839134 211502 1 0.01675 21150235 61940820 1 0.00000 6194082036 55132871 1 0.00802 5513287137 463946 1 0.03549 47534238 301010 1 0.01203 36015439 488042 1 0.01259 48804240 63607 1 0.04969 7739541 217519 1 0.04395 21751942 64716813 1 0.04289 6471681343 56487153 1 0.00000 5648715344 214188 1 0.02871 21418879Appendix A. TablesTable A.17: Objective value in $ for warmstart (d0max = 25,α = 10), ε = 0.05Problem Objective success ε LP mass value1 454426 1 0.01413 4551812 9351412 1 0.00000 99740943 316359 1 0.04485 3290404 505297 1 0.03544 5052975 60684 1 0.00833 688986 211057 1 0.01468 2110577 78347224 1 0.01769 769612638 56982170 1 0.00869 569262799 461811 1 0.02990 47289110 9498097 1 0.04354 1020450811 376427 1 0.00346 39481912 497561 1 0.02044 51919313 61213 1 0.01252 7207914 212554 1 0.02161 21255415 78347224 1 0.01769 7696126316 56975485 1 0.00857 5697548517 461947 1 0.03018 47730918 495973 1 0.01731 51637819 64484 0 0.05960 7959720 218355 1 0.04761 21781221 461946 1 0.03132 44889722 9512758 1 0.04502 951275823 300402 1 0.01003 30040224 505297 1 0.04573 50529725 60253 1 0.00124 6732226 209244 1 0.00614 20924427 61940820 1 0.00000 6194082028 55096639 1 0.00737 5509663929 456838 1 0.02048 46836730 9425696 1 0.03620 962026331 301197 1 0.01264 33904732 482671 1 0.00160 48267133 61213 1 0.01252 6839134 211502 1 0.01675 21150235 61940820 1 0.00000 6194082036 55132871 1 0.00802 5513287137 455873 1 0.01841 47534238 312180 1 0.04738 36015439 488042 1 0.01259 48804240 63540 1 0.04868 7739541 218062 1 0.04633 21751942 64716813 1 0.04289 6471681343 57840191 1 0.02339 5648715344 210680 1 0.01254 21418880Appendix A. TablesTable A.18: Objective value in $ for warmstart with ǫ (d0max = 25,α = 10),ε = 0.05Problem Objective success ε ǫ LP mass value1 454426 1 0.01413 0.01434 4551812 9351412 1 0.00000 0.02938 99740943 316359 1 0.04485 0.06244 3290404 505297 1 0.03544 0.03674 5052975 60684 1 0.00833 0.00840 688986 211057 1 0.01468 0.01489 2110577 78347224 1 0.00000 0.01801 769612638 56982170 1 0.00869 0.00876 569262799 461811 1 0.02990 0.03082 47289110 9498097 1 0.02020 0.04553 1020450811 376427 1 0.00346 0.26417 39481912 497561 1 0.02044 0.02087 51919313 61213 1 0.01252 0.01267 7207914 212554 1 0.02161 0.02209 21255415 78347224 1 0.00000 0.01801 7696126316 56975485 1 0.00857 0.00857 5697548517 461947 1 0.03018 0.03112 47730918 495973 1 0.01709 0.01761 51637819 64484 0 0.05960 0.06679 7959720 218355 1 0.04761 0.04999 21781221 461946 1 0.03132 0.03233 44889722 9512758 1 0.04502 0.04714 951275823 300402 1 0.01003 0.01013 30040224 505297 1 0.04573 0.04792 50529725 60253 1 0.00124 0.00124 6732226 209244 1 0.00614 0.00618 20924427 61940820 1 0.00000 0.00000 6194082028 55096639 1 0.00737 0.00742 5509663929 456838 1 0.02048 0.02091 46836730 9425696 1 0.03620 0.03756 962026331 301197 1 0.01264 0.01281 33904732 482671 1 0.00160 0.00161 48267133 61213 1 0.01252 0.01267 6839134 211502 1 0.01675 0.01703 21150235 61940820 1 0.00000 0.00000 6194082036 55132871 1 0.00802 0.00809 5513287137 455873 1 0.01841 0.01876 47534238 312180 1 0.04738 0.04974 36015439 488042 1 0.01259 0.01259 48804240 63562 1 0.04901 0.05154 7739541 218062 1 0.04633 0.04858 21751942 64716813 1 0.04289 0.04482 6471681343 57840191 1 0.02339 0.02395 5648715344 210680 1 0.01216 0.01270 21418881Appendix A. TablesTable A.19: Objective value in $ for warmstart with ǫ (d0max = 5,α = 5),ε = 0.05Problem Objective success ε ǫ LP mass value1 454426 1 0.01413 0.01434 4551812 9351412 1 0.02854 0.02938 99740943 316359 1 0.04655 0.06244 3290404 488986 1 0.00327 0.03674 5052975 60684 1 0.00833 0.00840 688986 211057 1 0.01468 0.01489 2110577 78347224 1 0.00000 0.01801 769612638 56576215 1 0.00157 0.00876 569262799 461811 1 0.02990 0.03082 47289110 9498097 1 0.04340 0.04553 1020450811 376427 1 0.03643 0.26417 39481912 501295 1 0.02774 0.02087 51919313 60687 1 0.00396 0.01267 7207914 219628 1 0.04992 0.02209 21255415 78499966 1 0.00195 0.01801 7696126316 56975485 1 0.00857 0.00857 5697548517 656851 0 0.05385 0.03112 47730918 495973 1 0.01709 0.01761 51637819 82236 0 0.06379 0.06679 7959720 217881 1 0.04553 0.04999 21781221 461946 1 0.03132 0.03233 44889722 9351412 1 0.02811 0.04714 951275823 297389 1 -0.00000 0.01013 30040224 488396 1 0.01271 0.04792 50529725 60253 1 0.00124 0.00124 6732226 208491 1 0.00255 0.00618 20924427 61940820 1 0.00000 0.00000 6194082028 57438725 1 0.04784 0.00742 5509663929 465460 1 0.03863 0.02091 46836730 9425696 1 0.03620 0.03756 962026331 297389 1 0.00000 0.01281 33904732 490197 1 0.01693 0.00161 48267133 60735 1 0.00474 0.01267 6839134 211607 1 0.01724 0.01703 21150235 61940820 1 0.00000 0.00000 6194082036 55231802 1 0.00980 0.00809 5513287137 455042 1 0.01662 0.01876 47534238 403569 0 0.06654 0.04974 36015439 560000 0 1.00000 0.01259 48804240 77847 0 0.07685 0.05154 7739541 217944 1 0.04581 0.04858 21751942 64716813 1 0.04289 0.04482 6471681343 57205059 1 0.01255 0.02395 5648715344 218432 1 0.04759 0.01270 21418882Appendix BFigures05101520-0.2 -0.1 0 0.1 0.2 0.35% mipgap05101520-0.2 -0.1 0 0.1 0.2 0.3Problems1% mipgapFigure B.1: Savings with direct solve, compared to mass block removalschedule83Appendix B. Figures05101520253035-0.2-0.15-0.1-0.05 0 0.050.10.150.25% mipgap0510152025-0.2-0.15-0.1-0.05 0 0.050.10.150.2Problems1% mipgapFigure B.2: Savings with direct solve and x◦, compared to mass block re-moval schedule0510152025-0.2-0.15-0.1-0.05 0 0.050.10.150.25% mipgap0510152025-0.2 -0.1 0 0.1 0.2 0.3Problems1% mipgapFigure B.3: Savings with warmstart (d0max = 25,α = 10), compared to massblock removal schedule05101520-0.2-0.15-0.1-0.05 0 0.050.10.150.25% mipgap05101520-0.1 -0.05 0 0.05 0.1 0.15 0.2Problems1% mipgapFigure B.4: Savings with warmstart and ǫ (d0max = 5,α = 5), compared tomass block removal schedule84Appendix B. Figures-0.0500.050.10.150.20 5 10 15 20 25 30 35 405% mipgap-0.0500.050.10.150.20 5 10 15 20 25 30 35 401% mipgapFigure B.5: Stratification of gains for direct solve; * are 5, ⊙ are 10, trianglesolidare 20 blocks; small symbols are 2, big symbols are 4 access roads; withinthe block groups, points are sorted from left to right in ascending order ofsection numbers00.050.10.150.20 10 20 30 40 505% mipgap00.050.10.150.20 10 20 30 40 501% mipgapFigure B.6: Stratification of gains for direct solve with x◦; green=5, blue=10,red=20 blocks; * are 5,⊙are 10, trianglesolid are 20 blocks; small symbols are 2, bigsymbols are 4 access roads; within the block groups, points are sorted fromleft to right in ascending order of section numbers85Appendix B. Figures-0.0500.050.10.150.20 10 20 30 40 505% mipgap-0.0500.050.10.150.20 10 20 30 40 501% mipgapFigure B.7: Stratification of gains for warmstart (d0max = 25,α = 10); *are 5, ⊙ are 10, trianglesolid are 20 blocks; small symbols are 2, big symbols are 4access roads; within the block groups, points are sorted from left to right inascending order of section numbers-0.0500.050.10.150.20 5 10 15 20 25 30 35 405% mipgap-0.0500.050.10.150.20 5 10 15 20 25 30 35 401% mipgapFigure B.8: Stratification of gains for warmstart with ǫ (d0max = 5,α = 5);* are 5, ⊙ are 10, trianglesolid are 20 blocks; small symbols are 2, big symbols are 4access roads; within the block groups, points are sorted from left to right inascending order of section numbers86
- 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-04-16
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) |
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