Road Design Optimization with aSurrogate FunctionbyYasha PushakA THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFB.SC. HONOURS IN COMPUTER SCIENCE HONOURS IN MATHEMATICSinIrving K. Barber School of Arts and Sciences(Computer Science and Mathematics)THE UNIVERSITY OF BRITISH COLUMBIA(Okanagan)April 2015c© Yasha Pushak, 2015AbstractPlanning new highways is a complex process requiring consideration ofseveral non-trivial cost factors. An initial path is first refined by adjustingits vertical alignment to reduce earthwork cutting, filling, and hauling costs.This optimization problem is then used as a sub-problem that needs to berepeatedly solved with different horizontal alignments, to find an optimalpath. The complete optimization process requires significant computationtime to solve real-world problems. A multi objective surrogate function,which approximates the earthwork and paving costs and can be quicklyevaluated, can be used to reduce the running time of the full optimizationalgorithm. This thesis builds upon an existing surrogate function, using a“Mass Haul” to approximate the earthwork hauling costs. Other constraintsmade to simplify the former surrogate are relaxed to improve the quality ofthe solutions returned by the surrogate. Numerical results compare solutionsreturned by the two versions of the surrogate and investigate the Pareto frontbetween the earthwork and paving costs.iiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiListings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . ixChapter 1: Introduction . . . . . . . . . . . . . . . . . . . . . . . 11.1 Road Design in the Literature . . . . . . . . . . . . . . . . . . 11.2 Research Goal . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.3 Thesis Outline . . . . . . . . . . . . . . . . . . . . . . . . . . 3Chapter 2: Existing Models . . . . . . . . . . . . . . . . . . . . 42.1 The Three-Dimensional Model . . . . . . . . . . . . . . . . . 42.1.1 The Earthwork Problem . . . . . . . . . . . . . . . . . 42.1.2 The Vertical Alignment Problem . . . . . . . . . . . . 52.1.3 The Horizontal Alignment Problem . . . . . . . . . . . 72.2 The Surrogate Function Hirpa-2014 . . . . . . . . . . . . . . . 92.2.1 Road Alignment Model . . . . . . . . . . . . . . . . . 92.2.2 Terrain Model . . . . . . . . . . . . . . . . . . . . . . 112.2.3 Cost Formulation . . . . . . . . . . . . . . . . . . . . . 12Chapter 3: Implementing the Surrogate . . . . . . . . . . . . . 153.1 Improving the Surrogate . . . . . . . . . . . . . . . . . . . . . 153.1.1 The Jacobian . . . . . . . . . . . . . . . . . . . . . . . 153.1.2 Computing the Integral . . . . . . . . . . . . . . . . . 173.1.3 Relaxing the Constant Curve Elevation Constraint . . 193.1.4 Mass-Haul . . . . . . . . . . . . . . . . . . . . . . . . . 22iiiTABLE OF CONTENTS3.2 Converting between the Full Model and the Surrogate . . . . 233.2.1 Converting the Terrain . . . . . . . . . . . . . . . . . . 243.2.2 Converting the Road Alignment . . . . . . . . . . . . 25Chapter 4: Numerical Tests . . . . . . . . . . . . . . . . . . . . 274.1 Comparing the Mass-Haul Algorithm . . . . . . . . . . . . . . 274.2 Pareto Front . . . . . . . . . . . . . . . . . . . . . . . . . . . 284.3 Cost Parameter Stability . . . . . . . . . . . . . . . . . . . . . 31Chapter 5: Conclusion . . . . . . . . . . . . . . . . . . . . . . . . 33Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38Appendix A: Source Code . . . . . . . . . . . . . . . . . . . . . . . 39A.1 Mass Haul Algorithm . . . . . . . . . . . . . . . . . . . . . . 39A.2 Mass Haul Costs . . . . . . . . . . . . . . . . . . . . . . . . . 42ivListings3.1 Source code written to calculate the closed form integrals forthe volume of cutting and filling of tangent sections in thesurrogate function Hirpa-2014. . . . . . . . . . . . . . . . . . 183.2 Source code written to calculate the closed form integrals forthe volume of cutting and filling of circular sections in thesurrogate function Hirpa-2014. . . . . . . . . . . . . . . . . . 183.3 Source code written to calculate the closed form integrals forthe volume of cutting and filling of curve sections in the sur-rogate function Pushak-2015. . . . . . . . . . . . . . . . . . . 20A.1 An Implementation of the Mass Haul algorithm. An emphasiswas made on code readability and documentation for futureuse. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39A.2 A helper function used by the Mass Haull Algorithm in list-ing A.1 to determine the costs of hauling some volume ofearth over a set distance. . . . . . . . . . . . . . . . . . . . . 42vList of FiguresFigure 2.1 A sample solution for the Earthwork Problem. Theelevation profile of the ground is the solid black line.The dotted black line is the elevation profile of theroad. . . . . . . . . . . . . . . . . . . . . . . . . . . . 5Figure 2.2 A summary of the Earthwork Problem. . . . . . . . . 5Figure 2.3 A sample solution for the Vertical Alignment Prob-lem. The elevation profile of the ground is the solidblack line. The dotted black line is the elevationprofile of the road. We see that in comparison toFigure 2.1, the modified vertical alignment has re-duced the total volume of cutting and filling, whichdecreases the total cost of the vertical alignment. . . 6Figure 2.4 A summary of the Vertical Alignment Problem. TheEarthwork Problem is used as a component of theVertical Alignment Problem. . . . . . . . . . . . . . . 7Figure 2.5 A sample solution for the Horizontal Alignment Prob-lem. The solid black curvy-linear line represents theinitial alignment selected by engineers. The gray linesare the cross sections which define the horizontal cor-ridor. The elevation of the ground is given at each ofthe points along the cross sections. The dotted blackline represents a hypothetical optimal solution to theHorizontal Alignment Problem. . . . . . . . . . . . . 7Figure 2.6 A summary of the Horizontal Alignment Problem.The Vertical Alignment Problem is treated as a black-box, for which the optimal solution is used as the ob-jective function for the Horizontal Alignment problem. 8Figure 2.7 The horizontal alignment model used by the surrogatefunction. . . . . . . . . . . . . . . . . . . . . . . . . . 9viLIST OF FIGURESFigure 2.8 The vertical alignment model used by the surrogatefunction. Here, the number of vertical segments perhorizontal tangent is mi = 5 for 1 ≤ i ≤ n. Notethat, when projected onto the xy-plane, the lengthof each vertical segment is exactly 1/mi the length ofthe corresponding horizontal tangent. The height ofthe road is constant over the length of each circularcurve, hence z5 = z6. . . . . . . . . . . . . . . . . . . 10Figure 2.9 Visualization of the piecewise-linear function used tomodel the elevation of the terrain. Here δx = δy = 10(meters). . . . . . . . . . . . . . . . . . . . . . . . . . 12Figure 2.10 On the left, the cross-sectional area of cut for a points1 that is in a cut region. On the right, the cross-sectional area of fill for a point s2 that is in a fillregion. . . . . . . . . . . . . . . . . . . . . . . . . . . 14Figure 3.1 A sample of an optimal solution found using the sur-rogate function with Nomad [AAC+]. The test wasrun with and without the constraint which imposesa constant elevation over the circular sections of theroad. The black line is the original surrogate with theconstraint and the white line is the relaxed one. Thex and y axes are in meters, the z axis is in decameters. 20Figure 3.2 The volume of earthwork that needs to be done alongthe length of the roads in Figure 3.1. The originalsurrogate is in black and the version with the relaxedcurve elevation constraint is in gray. A value abovezero indicates a region of cut and a value below zeroindicates a region of fill. . . . . . . . . . . . . . . . . 21Figure 3.3 The visual representation of using the distance andoffset coordinates of the terrain coordinate as the in-put to the surrogate function. . . . . . . . . . . . . . 24Figure 4.1 A comparison of the surrogate function with and with-out the Mass-Haul algorithm. The original optimalsurrogate road alignment is in white, the alignmentfound with the surrogate function using the Mass-Haul algorithm is in black. The x and y axes are inmeters, the z axis is in decameters. . . . . . . . . . . 28viiLIST OF FIGURESFigure 4.2 The volume of earthwork that needs to be done alongthe road alignments in Figure 4.1. A positive valueindicates a region of cut, while a negative value indi-cates the region of fill. The original optimal surrogateroad alignment is in gray, the alignment found withthe surrogate function using the Mass-Haul algorithmis in black. . . . . . . . . . . . . . . . . . . . . . . . . 29Figure 4.3 The plot of all solutions found by the weighted summethod. The boxed points are the ones that belongto the pareto front. Note the vertical scale does notstart at zero and all costs are within 10%. . . . . . . 30Figure 4.4 The plot of the pareto fronts found by the variousalgorithms with their respective parameters. . . . . . 31Figure 4.5 Five different optimal road alignments found with dif-ferent weights for the cost parameters. The white linerepresents a value of θ = 0 (all user costs) and theblack line represents a value of θ = 1 (all earthworkcosts). The three shades of the grey lines representsthe values 0.25, 0.5 and 0.75. The x and y axes arein meters, the z axis is in decameters. . . . . . . . . . 32viiiAcknowledgmentsI would like to thank my supervisors Dr. Yves Lucet and Dr. WarrenHare, as well as Softree Technical Systems Inc. for the input data used togather results and engineering advise on road design. Thank you for thefunding received through CFI for the CA2 lab, and the NSERC funding forthe computer used to gather the results.ixChapter 1IntroductionModern road design is a complex process requiring consideration of sev-eral different cost factors, be them political, environmental or financial. De-termining the optimal road alignment for new highways or forestry roads isa non-trivial problem which cannot be solved by hand. In order to accu-rately model and capture the different costs and constraints computers arerequired to aid the process and determine the optimal solution.Current numerical optimization techniques exist to find locally optimalsolutions to the road design problem. These methods yield satisfactory solu-tions, however require a significant amount of computation time to converge.This thesis is primarily focused with developing a surrogate function whichapproximates the costs of the full optimization problem, yet is quick to eval-uate. The intended use for the surrogate function is to be combined withexisting methods to speed up the overall optimization process.1.1 Road Design in the LiteratureThe first problem which needs to be solved when designing a new high-way or forestry road, is planning an initial alignment. When connecting astart and destination point with a road, it can be a nontrivial choice of de-ciding if the road should go around one side or another of a mountain, lake,or other obstacle. Currently this problem is done by hand, with engineerssketching an initial plan for a road alignment. Most existing work in roaddesign optimization has been done in refining the road alignments plannedby engineers. There has been some work done [PLH14] which addresses thisproblem by finding a set of locally optimal paths which serve as candidatescorridors for further refinement. Similar methods for finding spatially dis-tinct corridors have been done by Church et al. [CLL92, LC93], which waslater improved by Scaparra et al. [SCM14]. Another approach for this typeof problem was used by Zhang et al. [ZA05] using genetic algorithms.Once an initial corridor is selected, the alignment needs to be refined.Typically this is treated only as a 2-dimensional problem, seeking to findthe optimal horizontal alignment within a corridor. There have been sev-11.2. Research Goaleral different approaches used to solve the horizontal alignment problem.The typical formulation begins with some form of a locally optimal corri-dor and aims to find the optimal road alignment within the corridor. Thishas been done using, for example, genetic algorithms as done by Kang etal. [Kan08, KJS12] or Jha et al. [Jha03, JS04, JK06]. Mondal [Mon14] solvedthe problem using derivative-free optimization solvers as NOMAD [AAC+]and HOPSPACK [Lib]. The cost formulation for the road varies, but a com-mon approach is to separate the horizontal and vertical alignment problems.One way of doing this is as a bi-level optimization problem using the ver-tical alignment as a sub-optimization problem to the horizontal alignmentoptimization problem. One example of a bi-level approach was done byMondal [Mon14]. In a similar method to ones describes in this thesis, Lee etal. [LTL09] modeled the horizontal alignment as a piecewise-linear line andapplied safety constraints after finding the optimal piecewise-linear align-ment. Another take to improving the efficiency of the horizontal alignmentproblem was done by Kang et al. [KSJ07] using feasible gates to reduce thesearch space.Closely linked to the cost of the horizontal alignment, is the road’s ver-tical alignment, which is needed to determine the earthwork costs. In orderto simplify models, the vertical alignment problem is typically treated sepa-rately from the horizontal problem, or as a subproblem of the full alignmentproblem. It is primarily concerned with minimizing the earthwork costs as-sociated with constructing new roads, which are a significant cost factor inroad construction [JSJ06]. There are three components to the earthworkcosts: cutting, filling and hauling. The first two are static for a given roadalignment and are proportional to the volume of earthwork that needs to bedone along the road. The third cost is determined by the strategy used tomove the earth along the road. The goal is to minimize the cost of excavationand embankment, while finding the earth-hauling plan which minimizes thecost of transporting the material. Some work on minimizing the earthworkcosts was done by Hare et al. [HKL11], which was later expanded upon in amethod for solving the vertical alignment [HHLR14].1.2 Research GoalThis thesis builds upon an existing surrogate model developed by Hirpa[Hir14]. The surrogate function is designed to serve as a quick approximationof the construction costs for a given three-dimensional road alignment. Theultimate goal for the surrogate is to be combined with a bi-level optimization21.3. Thesis Outlineapproach developed by Mondal [Mon14]. Since the full three-dimensionalproblem is structured as a bi-level optimization problem, it requires fullysolving the vertical alignment for every horizontal alignment, which yieldsa slow algorithm. The goal is to combine the surrogate with this methodin a way that uses the surrogate function to quickly make progress, whileoccasionally using the full model to validate that the surrogate function ismaking steps in the correct direction.1.3 Thesis OutlineIn Chapter 1 we introduce the road design problems and provide moti-vation for construction and use of the surrogate function. Chapter 2 beginswith an overview of the full three-dimensional alignment problem and onemethod for how it has been solved in literature. Chapter 2.1 provides somebackground into how the earthwork, vertical and horizontal problems areused to form the full, three-dimensional model. In Chapter 2.2 we discussthe formulation for an existing surrogate function by Hirpa [Hir14]. Chap-ter 3 describes the work done for this thesis to implement the surrogate func-tion described in Chapter 2.2. Chapter 3.1 contains several modificationsand improvements made to the existing surrogate. Chapter 3.2 discussespossible methods for converting between the full model and the surrogatemodel. In Chapter 4 we briefly discuss some preliminary results obtainedwith the surrogate and provide concluding remarks and ideas for future workin Chapter 5.3Chapter 2Existing ModelsThis chapter provides an overview of the specific models that the lab isusing. We describe how the road and terrain are modeled and briefly coverthe cost formulation.2.1 The Three-Dimensional ModelThis model is used to refine an initial alignment provided by an engi-neer. An initial corridor is selected manually and the elevation data withinthe corridor is used as input for this model. A derivative-free solver isused to find the optimal road alignment within the corridor provided byengineers. This model has been developed as part of an ongoing researchproject [HKL11, HHLR14]. Most recently the horizontal component wascompleted by Mondal [Mon14].2.1.1 The Earthwork ProblemThe road design problem requires consideration of several different costfactors. Once a horizontal and vertical alignment have been chosen, we thenneed to find the optimal earthwork solution. Suppose we begin with thevertical alignment seen in Figure 2.1. The goal of the Earthwork Problemis to find the optimal strategy to move the earth along the road.Some of the material that needs to be cut can be used to fill in groundto bring it to the level of the road elsewhere. How far and how much of thismaterial is moved affects the cost of the earthwork solution. In some cases,the type of material that exists along the road needs to be replaced withbetter quality earth, which will increase the cost of the earthwork solutionas well. In order to build the final road, some material is required to bebrought to the road, while other material must be removed and depositedat waste pits. Depending on the distance a given volume of material needsto be moved, different types of trucks or bulldozers will be used, which againalters the cost of the final earthwork solution.42.1. The Three-Dimensional ModelCut RegionCut RegionFill Region Fill RegionWaste PitWaste PitWaste Pit𝑑𝑖𝑠𝑡𝑎𝑛𝑐𝑒ℎ𝑒𝑖𝑔ℎ𝑡Figure 2.1: A sample solution for the Earthwork Problem. The elevationprofile of the ground is the solid black line. The dotted black line is theelevation profile of the road.A summary of the earthwork problem is given in Figure 2.2. For moredetails on its formulation and how it can be solved the reader is referredto [HKL11].The Earthwork Problem Input: A vertical alignment and ground profile Minimize: The earth hauling costs Output: The cost of cutting, filling, and haulingFigure 2.2: A summary of the Earthwork Problem.2.1.2 The Vertical Alignment ProblemThe model used for the vertical alignment is a piecewise-quadratic func-tion. In order to meet minimum safety criteria, various constraints areimposed upon the vertical alignment. These constraints include a restric-tion on the maximum curvature of the vertical alignment and require thatthe alignment is both continuous and smooth.52.1. The Three-Dimensional Model𝑑𝑖𝑠𝑡𝑎𝑛𝑐𝑒ℎ𝑒𝑖𝑔ℎ𝑡Waste PitWaste PitWaste PitCut RegionCut RegionFill Region Fill RegionFigure 2.3: A sample solution for the Vertical Alignment Problem. Theelevation profile of the ground is the solid black line. The dotted black lineis the elevation profile of the road. We see that in comparison to Figure 2.1,the modified vertical alignment has reduced the total volume of cutting andfilling, which decreases the total cost of the vertical alignment.The cost formulation for the Vertical Alignment Problem is closely linkedwith the Earthwork Problem. By altering the vertical alignment, we canadjust the volume of cutting and filling that is required along the road.Minimizing these volumes reduces the excavation and embankment costs,as well as the volume of earthwork that needs to be hauled, see Figure 2.3.Hence, the cost of a given vertical alignment can be found by solving theEarthwork Problem.While the Earthwork Problem could then be treated as a black box andsolved as a subproblem of the Vertical Alignment problem, it is faster inpractice to solve the two by simultaneously optimizing the vertical alignmentand the earth hauling costs. The final output of the Vertical Alignmentproblem yields the vertical alignment with the optimal earthwork costs. Avisual summary of the Vertical Alignment Problem is presented in Figure 2.4.For more details on how the Vertical Alignment Problem is modeled andsolved, the reader is referred to [HHLR14].62.1. The Three-Dimensional ModelCost of theAlignmentThe Vertical Alignment Problem Input: A ground profile Minimize: The Earthwork costs Output: The vertical alignment with minimal costUpdateA vertical AlignmentThe Earthwork Problem Input: A vertical alignment and ground profile Minimize: The earth hauling costs Output: The cost of cutting, filling, and haulingFigure 2.4: A summary of the Vertical Alignment Problem. The EarthworkProblem is used as a component of the Vertical Alignment Problem.2.1.3 The Horizontal Alignment ProblemThe goal for the Horizontal Alignment Problem is to refine an initialalignment provided by engineers. To this end, a corridor around the align-ment is given which provides the feasible region for the road and the terrainmodel for the Horizontal Alignment Problem.𝑥𝑦Figure 2.5: A sample solution for the Horizontal Alignment Problem. Thesolid black curvy-linear line represents the initial alignment selected by en-gineers. The gray lines are the cross sections which define the horizontalcorridor. The elevation of the ground is given at each of the points alongthe cross sections. The dotted black line represents a hypothetical optimalsolution to the Horizontal Alignment Problem.72.1. The Three-Dimensional ModelThe Horizontal Alignment Problem Input: A horizontal corridor and terrain data Minimize: The Earthwork and Paving Costs Output: The horizontal and corresponding vertical alignment with minimal costCost of theAlignmentUpdateGround Profileof Horizontal AlignmentCost of theAlignmentThe Vertical Alignment Problem Input: A ground profile Minimize: The Earthwork costs Output: The vertical alignment with minimal costUpdateA vertical AlignmentThe Earthwork ProblemFigure 2.6: A summary of the Horizontal Alignment Problem. The VerticalAlignment Problem is treated as a black-box, for which the optimal solutionis used as the objective function for the Horizontal Alignment problem.The initial, horizontal road alignment is modeled as a piecewise curvy-linear line. At equally spaced intervals along this line we are given a cross-section which is perpendicular to the road. The feasible corridor for theHorizontal Alignment Problem is then defined by these cross-sections. Ad-ditionally, each cross-section is given a set of equally spaced points, whichare used to define the terrain model, see Figure 2.5. At each point, in-formation such as the elevation of the ground, the type and percentage ofmaterials present are specified.Given a horizontal alignment, the terrain data along the alignment canbe used as an input to the Vertical Alignment Problem. The cost of theHorizontal Alignment Problem includes both the earthwork costs from thevertical solution and the cost of paving the road. The Vertical AlignmentProblem can then be treated as a black box used in the objective functionof the Horizontal Alignment Problem. A derivative-free method is appliedto the Horizontal Alignment Problem to find the optimal solution. A visualsummary of the Horizontal Alignment Problem is presented in Figure 2.6.For more details on the model and implementation of the HorizontalAlignment problem, the reader is referred to Mondal [Mon14].82.2. The Surrogate Function Hirpa-20142.2 The Surrogate Function Hirpa-2014The surrogate function described in this section was designed by Hirpa [Hir14]to provide a quick approximation of the costs associated with building a newroad.2.2.1 Road Alignment ModelSimilar to the full model, the surrogate splits the road’s alignment intotwo parts: the horizontal alignment and the vertical alignment.The horizontal alignment is a piecewise circular-linear curve. The coor-dinates of the start and end of the road are fixed and a set of n intersectionpoints and radii define the road’s horizontal alignment. Let (xs, ys) and(xe, ye) be the coordinates for the start and end of the road respectively andlet {(x1, y1, r1), (x2, y2, r2)...(xn, yn, rn)} be the set of intersection points andtheir radii, as seen in Figure 2.7.Horizontal Road Alignment(𝑥1, 𝑦1)(𝑥𝑛, 𝑦𝑛)(𝑥𝑒 , 𝑦𝑒)(𝑥𝑠 , 𝑦𝑠)𝑟1𝑟𝑛𝑥𝑦 𝑇𝐶1𝑇𝐶𝑛 𝐶𝑇𝑛𝐶𝑇1Figure 2.7: The horizontal alignment model used by the surrogate function.We denote the segments between circular pieces as tangents. The tran-sition points from tangents to circular pieces form the set{TC1, TC2 , ..., TCn}92.2. The Surrogate Function Hirpa-2014and the transition points from circular pieces to tangents the set{CT1, CT2, ..., CTn}.These sets are not inputs to the model, as they can be computed from theintersection points and their radii.Also included in the horizontal model is the width of the road, w inmeters, which is assumed to be constant for the length of the road.The vertical alignment model adds elevation data to the horizontal align-ment to create the 3 dimensional alignment. The ith horizontal tangentsection is split into mi segments of equal length. These mi segments de-fine a piecewise-linear vertical alignment for each horizontal tangent, as inFigure 2.8.Vertical Road Alignment𝑧𝑧𝑠𝑧1𝑧2𝑧3𝑧4𝑧5 𝑧6…… 𝑧𝑒Figure 2.8: The vertical alignment model used by the surrogate function.Here, the number of vertical segments per horizontal tangent is mi = 5 for1 ≤ i ≤ n. Note that, when projected onto the xy-plane, the length of eachvertical segment is exactly 1/mi the length of the corresponding horizontaltangent. The height of the road is constant over the length of each circularcurve, hence z5 = z6.Alignment ConstraintsThere are several constraints imposed on the surrogate alignment, eitherfor simplicity, or to ensure that the road produced meets safety criteria.An implicit constraint enforced by the alignment’s formulation ensures that102.2. The Surrogate Function Hirpa-2014the road alignment is continuous. In practice a full road alignment mustalso be smooth, however this constraint is relaxed in the vertical model forsimplicity.In the horizontal alignment the transition point CTi must not pass thepoint TCi+1, in order to ensure a smooth transition between the curves.A minimum radius constraint is imposed on the horizontal alignment forsafety. To reduce the search space in the optimization process, a maximumcurvature radius is also specified.In the vertical alignment, the elevation distance between two consecutivepoints is bounded to ensure a maximum road grade. For simplicity in im-plementation, each horizontal tangent section is assumed to have the samemi = m number of vertical segments. One constraint of note was made tosimplify the cost calculations, which was to make the road elevation constantover the circular segments.2.2.2 Terrain ModelThe terrain model used for the surrogate function is a piecewise-linearfunction. This model was chosen for its simplicity and because a piecewisecontinuous function was needed to calculate the integrals in Chapter 2.2.3.Each linear plane of the model is a rectangle with fixed dimensions, δx andδy.The model, M(x, y) = z, is stored using three matrices, A, B and C,whereA =∂M∂x,B =∂M∂y,C = z(x, y)− x∂M∂x− y∂M∂y.We introduce the variablesix := floor(xδx) + 1,jy := floor(yδy) + 1,such thatA(ix, jy) =∂M∂x∣∣∣∣(x,y),112.2. The Surrogate Function Hirpa-2014Figure 2.9: Visualization of the piecewise-linear function used to model theelevation of the terrain. Here δx = δy = 10 (meters).B(ix, jy) =∂M∂y∣∣∣∣(x,y).The value of z for the model can then be found usingz = M(x, y) = xA(ix, jy) + yB(ix, jy) + C(ix, jy).2.2.3 Cost FormulationRoad construction contains many different factors contributing to thefinal cost of the road. For a typical new highway, two of the dominatingcosts are the earthwork and paving costs [JSJ06]. In order to capture bothof these cost factors, the surrogate is a multi-objective function that returnstwo cost components, Costearth and Costpaving.122.2. The Surrogate Function Hirpa-2014The cost calculations uses a parametric form of the road alignment givenby r(s), wherer(s) = (x(s), y(s), z(s)).Earthwork CostsThe surrogate earthwork costs include three parts: the cost of cutting,filling and the waste cost. The total earthwork costs are calculated asCostearth = VcCc + VfCf + VwCw,where Vc, Vf and Vw are the volume of cutting, filling and earth wasted incubic meters, respectively, and the user specified parameters Cc, Cf and Cware the costs of cutting, filling and waste costs per cubic meter, respectively.The waste costs are a simplification used in place of the earth haulingcosts. The waste volume is simply the difference of the volume of cuttingand filling, i.e.,Vw = |Vc − Vf |.This waste cost then applies a penalty to any excess earth that needs to bebrought to or removed from the road. We note that in this cost formulationthere are no earthwork hauling costs.In order to calculate the volume of cut and fill, we need the cross sectionalarea and length of the road. Since the road is defined as a piecewise circular-linear curve and the terrain is modeled as a piecewise-linear function, thesecalculations need to be split up into several sections. Hence we haveVc =∑i∈CR∫ biaiareac(s)ds,where areac(s) is the cross-sectional area of the road that needs to be cutat s and i ∈ CR if [ai, bi] brackets a cut region of the road, i.e., one wherethe ground needs to be cut. Similarly, we haveVf =∑i∈FR∫ biaiareaf (s)ds,where areaf (s) is the cross-sectional area of the road that needs to be filledat s and i ∈ FR if [ai, bi] brackets a fill region of the road.The cross-sectional area of the road that needs to be cut is calculated asa trapezoid defined by the two angles θ1 and θ2, the width of the road, w,132.2. The Surrogate Function Hirpa-2014and the distance of the road from the ground, hc(s), as in Figure 2.10. Wedenote zr(s) as the height of the road at s and zg(s) = M(x(s), y(s)) as theheight of the ground at (x(s), y(s)). We can find the distance of the roadfrom the ground as hc(s) = zg(s)− zr(s). Thenareac = hc(s)(w +12κhc(s)),where κ = cot(θ1) + cot(θ2).𝑤ℎ𝑐(𝑠1) 𝜃2𝜃1𝑎𝑟𝑒𝑎𝑐(𝑠1)𝑧𝑔(𝑠1)𝑧𝑟(𝑠1)𝑤ℎ𝑓(𝑠2)𝜃2𝜃1𝑎𝑟𝑒𝑎𝑓(𝑠2)𝑧𝑟(𝑠2)𝑧𝑔(𝑠2)Figure 2.10: On the left, the cross-sectional area of cut for a point s1 thatis in a cut region. On the right, the cross-sectional area of fill for a point s2that is in a fill region.Symmetrically we let hf (s) = zr(s) − zg(s), and calculate the cross-sectional area of road that needs to be filled asareaf = hf (s)(w +12κhf (s)).Paving CostsThe paving costs are simply calculated using the length and width of theroadCostpaving =∫wCpds,where Cp is the cost of paving the road per meter squared and w is thewidth of the road.14Chapter 3Implementing the SurrogateAn initial implementation of the Surrogate as described in Chapter 2.2was written by Hirpa [Hir14] in MATLAB [Inc]. This code was used asa proof of concept for the road design surrogate function. Beginning withthis prototype implementation as a reference, a new surrogate function wasdeveloped for this thesis.3.1 Improving the SurrogateIn order to use the surrogate function it needed to be converted fromMATLAB to C code. This change was required to reduce the running timeof the algorithm for use in the full algorithm. This process was done usinga MATLAB to C code conversion tool. In order to use this tool significantrefactoring of the original code was required.During refactoring, some implementation and design problems were un-covered and fixed. In particular there were two errors in how the calculationsfor the volume of cutting and filling were performed that were discovered.In order to ensure the quality of the final surrogate twenty-four unit testswere implemented. During the refactoring and testing some undesirablebehaviours of the solutions returned by the surrogate were observed. Mod-ifications to the surrogate were made to remove these defects and improvethe quality of solutions found with the surrogate. We denote the originalsurrogate design by Hirpa-2014 and the modified one with Pushak-2015.3.1.1 The JacobianThe original coordinates used to define the road alignment model arestandard Euclidean coordinates. As described in Chapter 2.2.3, these co-ordinates are transformed to parametric form to calculate the cost of thesurrogate alignment. Recall the parametric formr(s) = (x(s), y(s), z(s)),153.1. Improving the Surrogateand the calculation for the volume of cutVc =∑i∈CR∫ biaiareac(s)ds.Since this integral is computed using the parameter r(s) in place of (x, y, z),the conversion from dxdydz to ds requires the use of a Jacobian, which wasomitted in the original design. In this case, the Jacobian is computed as||r′(s)||, where || · || is the Euclidean norm, sodxdydz = ||r′(s)||ds = ||(x′(s), y′(s), z′(s))||ds.To simplify the implementation, we vary the parameter s from 0 to 1 foreach continuous region [ai, bi] that can be computed with a single integral.The parametric form of the equations for the tangent sections connectingpoints (xj , yj , zj) and (xj+1, yj+1, zj+1) can be written asx(s) = xj + (xj+1 − xj)s,y(s) = yj + (yj+1 − yj)s,z(s) = zj + (zj+1 − zj)s.So the Jacobian is||r′(s)|| = ||((xj+1 − xj), (yj+1 − yj), (zj+1 − zj))||= ||(xj+1, yj+1, zj+1)− (xj , yj , zj)||,which is simply the distance between the two points.For the circular sections of the road, the parametric equation is writtenasx(s) = Cj + rj cos(θ(s)),y(s) = Cj + rj sin(θ(s)),z(s) = zj = zj+1,where Cj is the center of the circle, rj is the radius of the circle and θ(s) =θj + (θj+1 − θj)s, such that r(0) = (xj , yj , zj) and r(1) = (xj+1, yj+1, zj+1).163.1. Improving the SurrogateThe Jacobian is then||r′(s)|| = ||(−rj sin(θ(s))θ′(s), rj cos(θ(s))θ′(s), 0)||=√(rjθ′(s) sin(θ(s)))2 + (rjθ′(s) cos(θ(s)))2=√(rjθ′(s))2(sin2(θ(s)) + cos2(θ(s)))=√(rjθ′(s))2= |rjθ′(s)|= |rj(θj+1 − θj)|.3.1.2 Computing the IntegralA closed form version of the volume of cut and fill integrals was usedin the prototype implementation of the surrogate function. These integralswere originally solved by hand and typed into the code. This process requirescareful attention to notation and is prone to human error.As an example, an easy error that was made stated that the volume offill was the negative of the volume of cut. While intuitively correct, thedirect implementation of the integral could not simply be negated. Recallfrom Chapter 2.2.3areac = hc(s)(w +12κhc(s)),hc(s) = zg(s)− zr(s)andareaf = hf (s)(w +12κhf (s)),hf (s) = zr(s)− zg(s).The correct relationship ishf (s) = −hc(s),which yieldsareaf = −hc(s)(w − 12κhc(s)) 6= −hc(s)(w + 12κhc(s)) = areac.To avoid human error, the symbolic MATLAB code in Listing 3.1 wasused to compute the closed form integrals that were used in the surrogate.173.1. Improving the Surrogate1 %Set up the symbol ic code .syms x1 x2 y1 y2 z1 z2 A B C kappa w j m s s1 s2 ;3%The implementation f o r the su r roga t e uses ( x1 , y1 ) and ( x2 , y2 )to r ep r e s en t CT( i ) and TC( i +1) . m i s the number o f v e r t i c a ltangent s e c t i on s , so z1 and z2 . Here A, B and C are the threematrix e n t r i e s that d e f i n e the t e r r a i n model in t h i s r eg i on .5 zg = ( (A∗x1 + B∗y1 + C) + s ∗ (A∗ ( x2−x1 ) + B∗ ( y2−y1 ) ) ) ;z r = ( z1 + ( z2−z1 ) ∗ (m∗ s+1− j ) ) ;7%Calcu la t i on o f the Jacobian . We d iv id e x2 and y2 by m to f i ndthe l ength o f the cur rent v e r t i c a l segment .9 Jacobian = norm ( [ x2/m y2/m z2 ] − [ x1 y1 z1 ] ) ;%Ca lcu la t e the tangent volume o f cut .11 hc = zg − zr ;ac = hc∗ (w+(1/2)∗hc∗kappa ) ;13 Vc = in t ( ac∗Jacobian , s , s1 , s2 ) ;%Ca lcu la t e the tangent volume o f f i l l .15 hf = zr − zg ;a f = hf ∗ (w+(1/2)∗hf ∗kappa ) ;17 Vf = in t ( a f ∗Jacobian , s , s1 , s2 ) ;Listing 3.1: Source code written to calculate the closed form integrals for thevolume of cutting and filling of tangent sections in the surrogate functionHirpa-2014.Similarly, the closed form integral for volume of cut and fill integralsfor the circular sections of the road can be calculated using the symbolicMATLAB code in Listing 3.2.1 %Set up the symbol ic code .syms x1 x2 y1 y2 z1 z2 r A B C kappa w j m s s1 s2 theta1 theta23%Set up theta ( s ) .5 angMag = theta2−theta1 ;theta = theta1 + angMag∗ s ;7%Calcu la t e the e l e v a t i o n o f the road and ground . Theimplementation f o r the su r roga t e uses ( x1 , y1 ) and ( x2 , y2 ) tor ep r e s en t TC( i ) and CT( i +1) . Here A, B and C are the threematrix e n t r i e s that d e f i n e the t e r r a i n model in t h i s r eg i on .9 zg = (A∗x1 + B∗y1 + C) + r ∗ (A∗ cos ( theta ) + B∗ s i n ( theta ) ) ;z r = z1 ;11%Calcu la t e the Jacobian .13 Jacobian = abs ( r ∗angMag) ;%Calcuate the volume o f cut f o r curves .15 hc = zg − zr ;183.1. Improving the Surrogateac = w∗hc + (1/2) ∗kappa∗hc ˆ2 ;17 Vc = in t ( ac∗Jacobian , s , s1 , s2 ) ;%Calcuate the volume o f f i l l f o r curves .19 hf = zr − zg ;a f = w∗hf + (1/2) ∗kappa∗hf ˆ2 ;21 Vf = in t ( a f ∗Jacobian , s , s1 , s2 ) ;Listing 3.2: Source code written to calculate the closed form integrals for thevolume of cutting and filling of circular sections in the surrogate functionHirpa-2014.3.1.3 Relaxing the Constant Curve Elevation ConstraintPreliminary tests using the new surrogate were performed by using NO-MAD [AAC+], a derivative-free optimization solver. An example of oneof the graphical solutions is displayed in Figure 3.1, which highlights anundesirable behaviour of the circular road sections. We see that in the orig-inal surrogate model the constraint requiring the elevation to be constantover the curve located on a steep hill causes the road alignment to extendfar underneath of the level of the ground. Figure 3.2 shows the volume ofearthwork that needs to be done along the road alignment in Figure 3.1.The two large bumps occur during the circular sections of the road, addingsignificant cost to the optimal solutions. The second curve in this imagecorresponds to the one noted in the first figure, while the other one is notvisible in the other image due to the angle.This behaviour of requiring increased earthwork around the circular sec-tions of the road was hypothesized to be dominating the other costs andhence was detracting from the quality of the solutions found with the sur-rogate. With the earthwork volume integrals being calculated symbolicallythere was no longer a need for a simplified model and so the constant ele-vation constraint was relaxed.For the circular sections of the road, the new parametric equation iswritten asx(s) = Cj + rj cos(θ(s)),y(s) = Cj + rj sin(θ(s)),z(s) = zj + (zj+1 − zj)s,where Cj is the center of the circle, rj is the radius of the circle and θ(s) =θj + (θj+1 − θj)s, such that r(0) = (xj , yj , zj) and r(1) = (xj+1, yj+1, zj+1).193.1. Improving the SurrogateFigure 3.1: A sample of an optimal solution found using the surrogate func-tion with Nomad [AAC+]. The test was run with and without the constraintwhich imposes a constant elevation over the circular sections of the road.The black line is the original surrogate with the constraint and the whiteline is the relaxed one. The x and y axes are in meters, the z axis is indecameters.The Jacobian is then||r′(s)|| = ||(−rj sin(θ(s))θ′(s), rj cos(θ(s))θ′(s), zj+1 − zj)||=√(rjθ′(s) sin(θ(s)))2 + (rjθ′(s) cos(θ(s)))2 + (zj+1 − zj)2=√(rj(θj+1 − θj))2 + (zj+1 − zj)2.Unlike the tangent sections of the road, the curves were only given anelevation at their beginning and end, defining a constant slope over thecurve. The results seen in figures 3.1 and 3.2 show that the single elevationchange was still able to yield solutions where the road was much better ableto match the elevation of the ground. Recalling that the road alignments arerestricted within corridors, we see in Figure 3.1 that the new road alignmentfound has moved substantially to the left.The updated Jacobian and parametric equations are solved with thesymbolic MATLAB code in Listing 3.3 to produce the final, closed formintegrals used in the new version of the surrogate function.1 %Set up the symbol ic code .syms x1 x2 y1 y2 z1 z2 r A B C kappa w j m s s1 s2 theta1 theta2203.1. Improving the SurrogateFigure 3.2: The volume of earthwork that needs to be done along the lengthof the roads in Figure 3.1. The original surrogate is in black and the versionwith the relaxed curve elevation constraint is in gray. A value above zeroindicates a region of cut and a value below zero indicates a region of fill.3%Set up theta ( s ) .5 angMag = theta2−theta1 ;theta = theta1 + angMag∗ s ;7%Calcu la t e the e l e v a t i o n o f the road and ground .9 zg = (A∗x1 + B∗y1 + C) + r ∗ (A∗ cos ( theta ) + B∗ s i n ( theta ) ) ;z r = z1 + ( z2−z1 ) ∗ s ;11%Calcu la t e the Jacobian .13 Jacobian = sq r t ( r ˆ2∗angMagˆ2 + ( z2−z1 ) ˆ2) ;%Calcuate the volume o f cut f o r curves with vary ing e l e v a t i o n .15 hc = zg − zr ;ac = w∗hc + (1/2) ∗kappa∗hc ˆ2 ;17 Vc = in t ( ac∗Jacobian , s , s1 , s2 ) ;%Calcuate the volume o f f i l l f o r curves with vary ing e l e v a t i o n .19 hf = zr − zg ;a f = w∗hf + (1/2) ∗kappa∗hf ˆ2 ;21 Vf = in t ( a f ∗Jacobian , s , s1 , s2 ) ;Listing 3.3: Source code written to calculate the closed form integrals forthe volume of cutting and filling of curve sections in the surrogate functionPushak-2015.213.1. Improving the Surrogate3.1.4 Mass-HaulIn some roads the earth hauling costs comprise roughly half of the totalearthwork costs [JSJ06]. Given this fact, we hypothesized that assumingthe earth hauling costs are negligible will significantly affect the solutionsreturned by the surrogate. A Mass-Haul algorithm was used to quicklyapproximate the cost of the optimal earth hauling solution. The Mass-Haul algorithm was chosen for its relatively low runtime complexity, despiteproviding a quick cost approximation.Algorithm 1: The Mass Haul Algorithm.input : An array of volumes, vol, and an array of distances, dist,each of size n.output: An estimate of the cost of the optimal earthwork solution./* Preprocess so that the earthwork balances out in theend. */if totalV olumeCut > totalV olumeFill thenFind largest vol(i) > 0;Move excess mass from vol(i) to nearest waste pit;elseFind largest vol(i) < 0;Move extra mass from nearest waste pit to vol(i);end/* The main loop for the Mass-Haul Algorithm. */j = 1;for i = 1 to n dowhile vol(i) > 0 dowhile vol(j) <= 0 doj++;end/* Guaranteed to find j > 0 because ofpreprocessing. */Move vol(i) to vol(j);endendThe Mass-Haul algorithm takes as input two arrays of length n. The firststores the volume of earthwork needing to be done along different segmentsof the road and the second stores the length of each segment of the road.223.2. Converting between the Full Model and the SurrogateThe cost of moving material depends on the volume of material and thedistance it needs to be moved. Earth that is moved a longer distance can bemoved with different kinds of trucks, while short movements can simply bedone with a bulldozer. Preprocessing can be done to determine the optimaltype of truck for a given distance. Additional input parameters are used inthe Mass-Haul algorithm to model these costs.An initial pass is done by the algorithm to determine whether or notthere is a greater volume of cutting or filling that needs to be done alongthe road. In the case where there is more cutting, the segment of the roadwith the most volume of cut is found and sufficient material is moved tothe nearest waste pit to ensure that the rest of the road’s earthwork will bebalanced, that is, the road can be made without using another waste pit.Symmetrically, when there is more filling required, the segment with thelargest volume of fill is found and extra material is brought to the road.A second pass is then done, which scans along the road until it findsthe first segment which requires earthwork. If cutting is required, then thenearest available fill segment is used to move the earth. The filling segmentsare handled symmetrically.The time complexity of the Mass-Haul algorithm is found by observingthat the both the preprocessing and main body of the algorithm requirea constant number of passes over the road. The preprocessing performs asingle pass to determine the position with the most mass and whether ornot the total volume of cutting or filling is greater. The main body of thealgorithm, while it appears to have multiple nested loops, is in fact also donein linear time. This is because it uses two variables, i and j which each passfrom 1 to n, tracking the positions of cut and fill respectively. This yields atime complexity of O(n), where n is the number of segments along the road.Pseudo-code for the mass-haul algorithm is included in Algorithm 1. Thefull source code available in appendix A.1.3.2 Converting between the Full Model and theSurrogateUsing the surrogate in combination with the full model requires defininga mapping between the full model’s road alignment and terrain data. Whilethe full optimization process has not yet been adapted to use this surro-gate, this section describes a two hypothetical conversion techniques thatcould be used. The sketch of the conversion techniques and their relativecomplexity are discussed in Chapter 3.2.1. In Chapter 3.2.2, we examine233.2. Converting between the Full Model and the Surrogatethe challenges that arise from the two methods when converting the roadalignments’ coordinates between the two models and present methods forthe conversion.3.2.1 Converting the TerrainThis section summarizes two different methods for converting the initialterrain corridor data provided by engineers into the piecewise linear stepfunction used by the surrogate.Option AThe first approach is the simplest method for converting the terrain. Inthe initial input data, the distance between each cross-section is always aconstant size. Additionally, the points are always equally spaced along thecross sections. Instead of using the proper Euclidean coordinates, we cansimply use the distance along the road as the x coordinate in the surrogate,and the offset from the initial alignment as the y coordinate. We denote thissystem using (d, o), where d is the distance and o is the offset. The visualinterpretation for this coordinate system is analogous to simply stretchingthe original coordinate system as seen in Figure 3.3.Figure 3.3: The visual representation of using the distance and offset coor-dinates of the terrain coordinate as the input to the surrogate function.One possible flaw with this model, is that it may not provide accuratesolutions for the curved portions of the corridor. Numerical testing needsto be done to determine whether or not this simplification will significantlyaffect the quality of the solutions found by the surrogate.Option BIn this section, we describe the sketch of an algorithm for another ap-proach to the terrain conversion problem. Here, we use the euclidean coor-243.2. Converting between the Full Model and the Surrogatedinates provided by the initial corridor instead of the distance and offset.Since the surrogate function still requires the evenly spaced step function,this could be generated using interpolation to find the planes of best fit.We begin by drawing a Euclidean rectangle that contains the entire cor-ridor. This rectangle is filled by a two-dimensional array, which will beused to form the new terrain model. A first pass is done, scanning each ofthe original terrain points and putting them into the positions of the two-dimensional array that correspond to the nearest Euclidean coordinates.Next, each position in the two dimensional array is scanned and the planeof best-fit is found to form the piecewise-linear function. If there are no datapoints near to one of the array entries, an elevation level of infinity can beused to indicate that this is outside of the corridor, and hence an infeasibleregion of the map.While this method does require more time to set up the terrain, it shouldprovide a more accurate representation of the terrain than the one presentedin Option A.3.2.2 Converting the Road AlignmentThe methods for converting the road alignment coordinates between theterrain models are described in Chapter 3.2.1. In this section we see thesimple method in Option A above requires additional work when convertingthe road alignments’ coordinates between the two models. In comparison,Option B yields a trivial conversion between the two coordinate systems.Option ATo convert a horizontal alignment in the full model to the surrogate, weconstruct a distance and offset coordinate system in the full model, whichtranslates to x and y coordinates in the surrogate. Since the full modeluses x and y coordinates, rather than distance and offset, it is non-trivialto convert a Euclidean (x, y) coordinate into a distance and offset (d, o)coordinate system.We begin by forming the set of boxes with smallest area that contain eachpair of adjacent, cross-sections for the initial intersection points in the fullmodel. To convert a given point from the full model, we find the boxes whichcontain the point. This narrows down the distance of the point to being onat most two tangent sections of your road. Using a similar strategy, we candetermine which pair of cross-sections (within the box) the point is between.The point can then be projected onto the two nearest cross-sections. We253.2. Converting between the Full Model and the SurrogateInterpolate the projected point’s d and o coordinates on each cross sectionby using the (d, o) coordinates of the two nearest data points on each of thecross sections. Using this information, we can interpolate the final distanceand offset coordinates used by the surrogate.Converting a horizontal alignment in the surrogate model to the fullmodel is comparatively simple. For each plane in the surrogate’s terrain,we keep track of the corresponding point in the original model. We roundthe (d, o) coordinates of the road to the nearest terrain data points andinterpolate.Option BIn this model, the coordinate systems used in the surrogate and thefull model are both Euclidean Coordinates. Hence, we do not require anyadditional techniques to convert the road alignments’ coordinates.26Chapter 4Numerical TestsSample numerical tests were run and are presented in this Chapter. Dueto time constraints comparing the solutions returned by the surrogate withthe full algorithm is left as future work.All numerical tests were done using input parameters for the constructioncosts provided by Softree Technical Systems Inc. The terrain data wasacquired from the USGS National Map Viewer [USG] and constructed usinga digital elevation model from California. Numerical tests were run on a64-bit desktop computer running Windows 7 with 32 GB of RAM. Theprocessor used was an Intel(R) Xeon(R) CPU E5-1620 v2 running at 3.70GHz. The surrogate was prototyped and tested in MATLAB R2014a, usingthe Optimization Toolbox.4.1 Comparing the Mass-Haul AlgorithmThe results from a sample solution found with two versions of the sur-rogate give indication for an improvement in the quality of solutions foundwith the surrogate using the Mass-Haul algorithm. Two versions of the sur-rogate were compared using Nomad [AAC+], a derivative-free optimizationsolver to find the optimal solutions. The first used the Mass-Haul algorithm,while the other used the original waste cost calculations instead. Both usedthe relaxed constraint for the circular sections of the road.The two solutions found with the surrogate are presented in Figure 4.1.Recall that the intersection points of the road are restricted within boxes tocapture the corridor-like constraint of the full algorithm. In this case, theboxes are squares with side lengths of 100 meters. From this observation wecan see that the optimal road alignments found using the two version of thesurrogate function are spatially very different.In Figure 4.2 we see the plot of earthwork that needs to be done alongthe road alignments in Figure 4.1. In particular we have labeled three pointsof note along the road alignments. In the first and third, we see that thereare two downward spikes, corresponding to regions of fill in the solutionreturned by the original surrogate. In the original alignment, we note that274.2. Pareto FrontFigure 4.1: A comparison of the surrogate function with and without theMass-Haul algorithm. The original optimal surrogate road alignment is inwhite, the alignment found with the surrogate function using the Mass-Haul algorithm is in black. The x and y axes are in meters, the z axis is indecameters.there is nowhere nearby that needs to be cut and hence the material neededto fill Points 1 and 3 must be brought from relatively far away. In thesecond alignment, which uses the Mass-Haul algorithm, such alignments arepenalized. Indeed, we see these downward spikes are no longer present forthe new alignment. On the other hand, we note that at Point 2 we seein both alignments a spikes going both up and down. In this case, whilethe total volume of earthwork required is larger than at Points 1 or 3, wesee that the same general trend is still present in the new alignment. TheMass-Haul algorithm does not penalize such alignments as heavily, since theearth in these regions only needs to be moved a short distance to level outthe road.4.2 Pareto FrontBuilding upon the original surrogate function Hirpa et al. [HHLT14]determined the pareto front between the user costs and the earthwork costsfor building new roads. The pareto front is defined to be the set of non-dominated points of a multi-objective function. That is, given a set ofcandidate road alignments the surrogate function will map these alignment284.2. Pareto Front0 100 200 300 400 500 600 700 800 900−4−20246810123Figure 4.2: The volume of earthwork that needs to be done along the roadalignments in Figure 4.1. A positive value indicates a region of cut, whilea negative value indicates the region of fill. The original optimal surrogateroad alignment is in gray, the alignment found with the surrogate functionusing the Mass-Haul algorithm is in black.to a vector of cost components. A point belongs to the set of pareto-optimalpoints if it is a non-dominated point, i.e., there does not exist an alignmentwhose cost vector has each component strictly less than that of the non-dominated point.The original version of the surrogate function returned a vector of twocosts, the earthwork costs and the paving costs. Which this was used todetermine the pareto front, the paving costs were later relabeled as the usercosts. The user costs are determined by the amount of time required to travelalong a road and hence, similar to the paving costs, are length dependent.Using a similar strategy, an updated version of the pareto front has beenfound using the new surrogate function. Numerical results were done withthree different solvers. The first was using a Multi-Objective Genetic Algo-rithm (MOGA), which is available as a solver in the MATLAB software. Thesecond used a derivative-free solver known as a Direct-Multisearch [CMVV11],or DMS. Finally, a weighted-sum method (WS) was used with the fminconsolver available in MATLAB. Since the fmincon solver optimizes a single-valued objective function, a parameter θ was varied from 0 to 1 at intervalsof 0.01. The objective value used was then cost = θcostearthwork + (1 −θ)costuser. The MOGA solver used a stopping criterion based on the toler-294.2. Pareto Frontance for the changes in the objective function. The DMS and WS methodsused as a stopping criterion the maximum number of function evaluations.Two choices of parameters were chosen for the MOGA and DMS solversin finding the pareto front. The results are summarized in table 4.1. Formore information on the parameter selection the reader is referred to [HHLT14].Table 4.1: A summary of the parameter choices and results for the solversused to find the pareto front.Solver Tolerance function evaluations time (seconds)MOGA 10−4 21,833 386MOGA 10−8 960,115 16,821DMS N/A 24,000 444DMS N/A 800,000 54,306WS (100 runs) N/A 24,000 (each) 8463 (total)103 104 1054200440046004800500052005400Earthwork Costs User Costs WSWS (pareto)Figure 4.3: The plot of all solutions found by the weighted sum method.The boxed points are the ones that belong to the pareto front. Note thevertical scale does not start at zero and all costs are within 10%.In Figure 4.3 we present the plot of all of the solutions along with thepareto front found using the weighted sum method. Figure 4.4 contains thepareto fronts found by all three algorithms, including the two parameter304.3. Cost Parameter Stabilitychoices for the MOGA and DMS algorithm. We see that MOGA with atolerance of 10−8, DMS with 800,000 function evaluations and WS return arelatively similar pareto front, indicating that the methods are convergingto the same solution.103 104 10543004400450046004700480049005000Earthwork CostsUser Costs WS (pareto)DMS (f−evals = 24,000)MOGA (Tolerance=1e−4)DMS (f−evals = 800,000)MOGA (Tolerance=1e−8)Figure 4.4: The plot of the pareto fronts found by the various algorithmswith their respective parameters.4.3 Cost Parameter StabilityA similar approach to the Weighted Sum method in Chapter 4.2 was usedto find the optimal solutions for roads with five different weights to the costparameters. These solutions were found using NOMAD with the surrogateformulation Pushak-2015. The objective function used for the solutions wassimilar to the weighted sum method, however only five different values forθ were used, i.e., 0, 0.25, 0.5, 0.75 and 1.We see that when only the user costs are considered, the road align-ment takes the straightest possible path within the corridor constraints, asexpected. The other four alignments show a variety of different solutions,indicating that the exact choice of cost parameters significantly affects theoptimal road alignment found.314.3. Cost Parameter StabilityFigure 4.5: Five different optimal road alignments found with differentweights for the cost parameters. The white line represents a value of θ = 0(all user costs) and the black line represents a value of θ = 1 (all earthworkcosts). The three shades of the grey lines represents the values 0.25, 0.5 and0.75. The x and y axes are in meters, the z axis is in decameters.32Chapter 5ConclusionThe next step for the surrogate function is more rigorous testing to com-pare the quality of its optimal solutions to those of the full three-dimensionalmodel. To isolate the components of the surrogate function for comparison,a first test would compare only the vertical alignment solution quality. Thiscould be done by running the full model’s vertical alignment optimizationproblem, and comparing it to the surrogate’s optimal solution using a fixedhorizontal alignment. Such testing would help to reveal the impact of theMass-Haul algorithm compared to the original design of the surrogate func-tion.Provided that the surrogate function’s vertical solutions are promising,a second set of tests need to be done comparing the solutions for the fullthree-dimensional road alignment. The results of which, while perhaps moreinteresting for the final intent of the surrogate, may prove less telling thanthe vertical alignment tests. The added complexity and increased numberof variables in the full, three-dimensional problem is likely to cause thesurrogate function and the full problem to converge to different solutions.For proof of the fact that the function can converge to many different localminima, we only need to consider the example in Chapter 4.3. With thisexample in mind, it is a likely assumption that the surrogate and full problemwill find different local minima. A metric that may prove more telling, isto compare the price of the optimal solutions found with the surrogate andthe full algorithm, when both are evaluated using the full model.We propose an alternative test which may help to determine whether thesurrogate captures similar local minima to the full algorithm. The surrogatefunction could be given the optimal solution found by the full algorithm asits initial alignment. If this alignment were to be significantly altered whenoptimized by the surrogate function, it would indicate that the two modelsdo not contain the same argminima. If this were to be determined to be thecase, the design of the surrogate model would need to be re-examined andimproved before integration with the full model.Assuming further testing yields positive results, there are a few possi-ble directions for the surrogate. Future plans for the full model include33Chapter 5. Conclusiondeveloping a three-dimensional model which simultaneously optimizes thehorizontal and vertical alignment. Should this be done, the surrogate func-tion may prove easy to integrate. However, since the surrogate model isable to optimize both the horizontal and vertical components, while thecurrent full model performs these steps separately, combining the two is notas simple. A simple, but likely inefficient option, would be to separate thetwo problems with the surrogate similar to the full model. First taking astep in the horizontal direction, and then quickly approximating the optimalvertical solution with the surrogate. While this method would prove easi-est in implementation, it does not take advantage of the potential availablethrough the surrogate function to simultaneously optimize both the verticaland horizontal alignments.Alternatively, the reverse type of approach could be done, by simultane-ously optimizing the vertical and horizontal alignments with the surrogatefunction. When the full model is then used as the objective function, thevertical component of the alignment could either be simply ignored, or usedas a warm start in solving the vertical problem. In either case, the currentiterate’s vertical alignment would not be the one at which the full model’sobjective function is evaluated. With this method, should the vertical align-ments produce by the surrogate be far from the optimal vertical alignments,then the step with the full model should indicate that the surrogate is notmaking progress, alerting the algorithm that a drastic step is necessary.Finally, the precision of the surrogate function could be further improvedby evaluating earthwork using a piecewise quadratic vertical alignment in-stead of a piecewise linear one. While the computation is more complicated,symbolic computation should help and result in a comparable computationtime.34Bibliography[AAC+] M.A. Abramson, C. Audet, G. Couture, J.E. Dennis, Jr.,S. Le Digabel, and C. Tribes. The NOMAD project. Softwareavailable at https://www.gerad.ca/nomad/. → pages vii, 2,19, 20, 27[CLL92] Richard L. Church, Scott R. Loban, and Kristi Lombard. Aninterface for exploring spatial alternatives for a corridor locationproblem. Computers & Geosciences, 18:1095–1105, 1992. →pages 1[CMVV11] A.L. Custo´dioo, J. F. A. Madeira, A. I. F. Vax, and L.N. Vicente.Direct multisearch for multiobjective optimization. SIAM, 2011.→ pages 29[HHLR14] Warren Hare, Shahadat Hossain, Yves Lucet, and Faisal Rah-man. Models and strategies for efficiently determining an opti-mal vertical alignment of roads. Computers & Operations Re-search, 44(0):161 – 173, 2014. → pages 2, 4, 6[HHLT14] Dessalegn Hirpa, Warren Hare, Yves Lucet, and SolomonTesfamariam. A multiobjective optimization framework forthree-dimensional road alignment design. Masters Thesis, 2014.→ pages 28, 30[Hir14] D. Hirpa. Simultaneous optimization of vertical and horizon-tal road alignments. Master’s thesis, The University of BritishColumbia, 2014. → pages 2, 3, 9, 15[HKL11] Warren L. Hare, Valentin R. Koch, and Yves Lucet. Modelsand algorithms to improve earthwork operations in road designusing mixed integer linear programming. European Journal ofOperational Research, 215(2):470 – 480, 2011. → pages 2, 4, 5[Inc] The Mathworks Inc. Software available at https://www.mathworks.com. [link]. → pages 1535Bibliography[Jha03] Manoj K. Jha. Criteria-based decision support system for select-ing highway alignments. Journal of Transportation Engineering,129(1):33–41, 2003. → pages 2[JK06] Manoj K. Jha and Eungcheol Kim. Highway route optimizationbased on accessibility, proximity, and land-use changes. Journalof Transportation Engineering, 132(5):435–439, 2006. → pages2[JS04] MK Jha and P Schonfeld. A highway alignment optimizationmodel using geographic information systems. TransportationResearch Part A, 38(6):455–481, 2004. → pages 2[JSJ06] Manoj Kumar Jha, Paul Schonfeld, and J-C Jong. IntelligentRoad Design, volume 19. WIT Press, 2006. → pages 2, 12, 22[Kan08] Min Wook Kang. An alignment optimization model for a simplehighway network. PhD thesis, University of Maryland, 2008. →pages 2[KJS12] Min-Wook Kang, Manoj K. Jha, and Paul Schonfeld. Applica-bility of highway alignment optimization models. TransportationResearch Part C: Emerging Technologies, 21(1):257 – 286, 2012.→ pages 2[KSJ07] Min Wook Kang, Paul Schonfeld, and Jyh-Cherng Jong. High-way alignment optimization through feasible gates. Journal ofAdvanced Transportation, 41(2):115–144, 2007. → pages 2[LC93] Kristi Lombard and Richard L. Church. The gateway shortestpath problem: Generating alternative routes for a corridor lo-cation problem. Geographical Systems, 1:22–45, 1993. → pages1[Lib] Sandia National Libraries. Hybrid Optimization Parallel SearchPACKage. Software available at https://software.sandia.gov/trac/hopspack/. → pages 2[LTL09] Yusin Lee, You-Ren Tsou, and Hsiao-Liang Liu. Optimizationmethod for highway horizontal alignment design. Journal ofTransportation Engineering, 135(4):217–224, 2009. → pages 236Bibliography[Mon14] Sukanto Mondal. Horizontal alignment optimization in road de-sign. Master’s thesis, University of British Columbia, 2014. →pages 2, 3, 4, 8[PLH14] Y. Pushak, Y. Lucet, and W. Hare. Multiple-path selection fornew highway alignments using discrete algorithms. EuropeanJournal of Operational Research, Submitted on Oct 23, 2014.→ pages 1[SCM14] Maria P. Scaparra, Richard L. Church, and F. Antonio Medrano.Corridor location: the multi-gateway shortest path model. Jour-nal of Geographical Systems, 16:287–309, 2014. → pages 1[USG] USGS. The national map viewer.http://viewer.nationalmap.gov/. → pages 27[ZA05] Xingdong Zhang and Marc P Armstrong. Using a genetic algo-rithm to generate alternatives for multi-objective corridor loca-tion problems. Proceedings of GeoComputation 2005, 2005. →pages 137Appendix38Appendix ASource CodeA.1 Mass Haul Algorithm1 f unc t i on [ cost , volume ] = massHaul ( volume , d i s tance , c l f , c lo ,c l e , chf , cho , che , l f , l o )%Mass haul func t i on c rea ted to be t t e r approximate the co s t o fthe earthwork f o r the su r roga t e func t i on .3 %Author : Yasha Pushak%Last Edited : Apr i l 26 , 20155 %volume − An array conta in ing the volume o f earthwork beingmoved at var i ous po in t s a long the road . The array i s paddedwith z e ro s at the end f o r speed .%d i s t anc e − The length o f each s e c t i o n correspond ing to thevolume array . The array i s padded with z e ro s at the end f o rspeed .7 %c l f − The co s t o f l oad ing f o r a Free Haul%c l o − The co s t o f l oad ing f o r an Over Haul9 %c l e − The co s t o f l oad ing f o r an End Haul%chf − The co s t o f hau l ing f o r a Free Haul11 %cho − The co s t o f hau l ing f o r an Over Haul%che − The co s t o f hau l ing f o r and End Haul13 %l f − The maximum length f o r which a Free Haul i s cheapest .%l o − The maximum length f o r which an Over Haul i s cheaper thanan End Haul .15%Prepros s ing17 %Find out how many non−zero e n t r i e s are in the ar rays .n = sum( d i s t ance ˜=0) ;19%CumulativeDistance i s used in the p r ep ro c e s s i ng f o r speed .21 cumulat iveDistance = ze ro s ( s i z e ( d i s t ance ) ) ;cumulat iveDistance (1 ) = d i s t ance (1 ) ;23 f o r i = 2 : numel ( d i s t ance )cumulat iveDistance ( i ) = cumulat iveDistance ( i −1) + d i s t anc e ( i ) ;25 end27 %Preproces s to remove any exce s s mass . Move the exce s s mass fromthe po int with the most mass to the nea r e s t p i t ( e i t h e r thes t a r t o f the end ) .39A.1. Mass Haul Algorithm29 %Find how much we have l e f t over .ex c e s s = sum( volume ) ;31 %Find where most o f the mate r i a l i s coming from .i f ( exc e s s > 0)33 maxv = −1;argmaxv = −1;35 f o r i = 1 : ni f (maxv < volume ( i ) )37 maxv = volume ( i ) ;argmaxv = i ;39 endend41 argmost = argmaxv ;e l s e i f ( exc e s s < 0)43 minv = 1 ;argminv = −1;45 f o r i = 1 : ni f (minv > volume ( i ) )47 minv = volume ( i ) ;argminv = i ;49 endend51 argmost = argminv ;end53 volume ( argmost ) = volume ( argmost ) − exce s s ;%We assume we are moving the earth from the middle o f the bucket.55 haulDis t = d i s t ance ( argmost ) /2 ;57 i f ( cumulat iveDistance ( argmost ) < cumulat iveDistance (n)−cumulat iveDistance ( argmost ) )%The c l o s e s t p i t i s as the s t a r t .59 f o r i = 1 : argmost−1%This loop could be r ep laced us ing the in fo rmat ion incumulat iveDistance , but s i n c e t h i s code i s intended to beused l a t e r , I have l e f t i t t h i s way f o r r e a d ab i l i t y .61 haulDis t = haulDis t + d i s t ance ( i ) ;end63 e l s e%The c l o s e s t p i t i s at the end .65 f o r i = argmost+1:nhaulDis t = haulDis t + d i s t ance ( i ) ;67 endend69 %Find the co s t o f hau l ing the earth to the waste p i t .d epo s i t = abs ( exce s s ) ;71 co s t = haulCost ( depos i t , haulDist , c l f , c lo , c l e , chf , cho , che , l f , l o ) ;40A.1. Mass Haul Algorithm73 %The main body o f the a lgor i thm .j = 1 ;75 done = 0 ;f o r i = 1 : n77 i f ( done )break ;79 endwhi l e ( volume ( i ) > 0)81 whi le ( volume ( j ) >= 0)j = j + 1 ;83 i f ( j == n+1)%The volume ar rays are padded with extra z e r o s f o re f f i c i e n c y , we want to stop when we have reached the end .85 done = 1 ;break ;87 endend89%I f j reached the end o f the array we are done .91 i f ( done )break ;93 end95 %Assume we are moving the mass from the cente r o f eachs e c t i o nhaulDis t = d i s t ance ( i ) /2 + d i s t ance ( j ) /2 ;97 i f ( j > i )hau lDis t = haulDis t + cumulat iveDistance ( j−1) −cumulat iveDistance ( i ) ;99 e l s ehau lDis t = haulDis t + cumulat iveDistance ( i −1) −cumulat iveDistance ( j ) ;101 end%we can ’ t depo s i t more than there i s space for , nor can wedepos i t more than we are hau l ing .103 depos i t = min ( abs ( volume ( j ) ) , abs ( volume ( i ) ) ) ;%Move the earth .105 volume ( i ) = volume ( i ) − depos i t ;volume ( j ) = volume ( j ) + depos i t ;107 %Find the co s t o f the movement .co s t = cos t + haulCost ( depos i t , haulDist , c l f , c lo , c l e , chf , cho ,che , l f , l o ) ;109 endend111 endListing A.1: An Implementation of the Mass Haul algorithm. An emphasiswas made on code readability and documentation for future use.41A.2. Mass Haul CostsA.2 Mass Haul Costs1 f unc t i on [ co s t ] = haulCost ( depos i t , haulDist , c l f , c lo , c l e , chf , cho ,che , l f , l o )%Helper func t i on to c a l c u l a t e the hau l ing co s t g iven the amountbeing3 %hauled and the d i f f e r e n t c o s t s f o r haul ing , loading , and t h e i r%correspond ing l eng th s .5 %Author : Yasha Pushak%Last Edited : Apr i l 5 , 20157 %depos i t − The volume o f earthwork being moved%haulDis t − The d i s t anc e the earth i s be ing moved9 %c l f − The co s t o f l oad ing f o r a Free Haul%c l o − The co s t o f l oad ing f o r an Over Haul11 %c l e − The co s t o f l oad ing f o r an End Haul%chf − The co s t o f hau l ing f o r a Free Haul13 %cho − The co s t o f hau l ing f o r an Over Haul%che − The co s t o f hau l ing f o r and End Haul15 %l f − The maximum length f o r which a Free Haul i s cheapest .%l o − The maximum length f o r which an Over Haul i s cheaper thanan End Haul .17i f ( hau lDis t < l f )19 %Free Haulco s t = c l f ∗ depos i t + chf ∗ depos i t ∗ haulDis t ;21 e l s e i f ( hau lDis t < l o )%Over Haul23 co s t = c l o ∗ depos i t + cho∗ depos i t ∗ haulDis t ;e l s e25 %End Haulco s t = c l e ∗ depos i t + che∗ depos i t ∗ haulDis t ;27 endendListing A.2: A helper function used by the Mass Haull Algorithm inlisting A.1 to determine the costs of hauling some volume of earth overa set distance.42
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Undergraduate Research /
- Road design optimization with a surrogate function
Open Collections
UBC Undergraduate Research
Road design optimization with a surrogate function Pushak, Yasha Apr 30, 2015
pdf
Page Metadata
Item Metadata
Title | Road design optimization with a surrogate function |
Creator |
Pushak, Yasha |
Date Issued | 2015-04 |
Description | Planning new highways is a complex process requiring consideration of several non-trivial cost factors. An initial path is first refined by adjusting its vertical alignment to reduce earthwork cutting, filing, and hauling costs. This optimization problem is then used as a sub-problem that needs to be repeatedly solved with different horizontal alignments, to find an optimal path. The complete optimization process requires significant computation time to solve real-world problems. A multi objective surrogate function, which approximates the earthwork and paving costs and can be quickly evaluated, can be used to reduce the running time of the full optimization algorithm. This thesis builds upon an existing surrogate function, using a \Mass Haul" to approximate the earthwork hauling costs. Other constraints made to simplify the former surrogate are relaxed to improve the quality of the solutions returned by the surrogate. Numerical results compare solutions returned by the two versions of the surrogate and investigate the Pareto front between the earthwork and paving costs. |
Type |
Text |
Language | eng |
Series | University of British Columbia. COSC 449 |
Date Available | 2017-02-01 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivs 2.5 Canada |
DOI | 10.14288/1.0224104 |
URI | http://hdl.handle.net/2429/56974 |
Affiliation |
Other UBC |
Peer Review Status | Unreviewed |
Scholarly Level | Undergraduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/2.5/ca/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 52966-pushak_yasha_COSC_449_Spring2015.pdf [ 1.87MB ]
- Metadata
- JSON: 52966-1.0224104.json
- JSON-LD: 52966-1.0224104-ld.json
- RDF/XML (Pretty): 52966-1.0224104-rdf.xml
- RDF/JSON: 52966-1.0224104-rdf.json
- Turtle: 52966-1.0224104-turtle.txt
- N-Triples: 52966-1.0224104-rdf-ntriples.txt
- Original Record: 52966-1.0224104-source.json
- Full Text
- 52966-1.0224104-fulltext.txt
- Citation
- 52966-1.0224104.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:
https://iiif.library.ubc.ca/presentation/dsp.52966.1-0224104/manifest