UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Simultaneous optimization of vertical and horizontal road alignments Hirpa, Dessalegn Amenu 2014

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


24-ubc_2014_september_hirpa_dessalegn.pdf [ 1.38MB ]
JSON: 24-1.0074354.json
JSON-LD: 24-1.0074354-ld.json
RDF/XML (Pretty): 24-1.0074354-rdf.xml
RDF/JSON: 24-1.0074354-rdf.json
Turtle: 24-1.0074354-turtle.txt
N-Triples: 24-1.0074354-rdf-ntriples.txt
Original Record: 24-1.0074354-source.json
Full Text

Full Text

Simultaneous optimization of verticaland horizontal road alignmentsbyDessalegn Amenu HirpaB.Sc., Bahir Dar University, 2004PGD., African Institute for Mathematical Sciences, 2011A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF SCIENCEinTHE COLLEGE OF GRADUATE STUDIES(Mathematics)THE UNIVERSITY OF BRITISH COLUMBIA(Okanagan)July 2014c© Dessalegn Amenu Hirpa, 2014AbstractOptimization of three-dimensional road alignments is recognized as anonlinear and nonconvex optimization problem. The development of mod-els that fully optimize a three-dimensional road alignment problem is notyet successful, because there are many factors involved and complexities inthe geometric specification of the alignment. At present, there are two op-timization approaches, the models that simultaneously optimize horizontaland vertical alignments, and those employing two or more stages of opti-mization processes.In this thesis, we develop a novel approach of solving a three dimen-sional road alignment problem where the optimal horizontal and verticalalignments are determined simultaneously. We develop the surrogate costmodel that approximate the earthwork cost and the pavement cost. Theproblem is modelled as a multiobjective optimization where the cost dueto the length of the road, and the cost due to the volume of earthwork arefound to be conflicting objectives. In order to study the proposed model,two case studies are tested and the numerical results are provided. Theexperimental results indicate that the problem is nonconvex, and that it is,indeed, a multiobjective optimization problem. Further developments andimprovements in the area of cost penalty parameters are recommended forfuture work.iiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . viiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiChapter 1: Introduction . . . . . . . . . . . . . . . . . . . . . . . 11.1 Background and motivation . . . . . . . . . . . . . . . . . . . 11.1.1 The horizontal alignment . . . . . . . . . . . . . . . 21.1.2 The vertical alignment . . . . . . . . . . . . . . . . . 41.1.3 The three-dimensional alignment . . . . . . . . . . . . 51.1.4 Surrogate models . . . . . . . . . . . . . . . . . . . . . 61.2 Problem statement . . . . . . . . . . . . . . . . . . . . . . . . 71.3 Research approach . . . . . . . . . . . . . . . . . . . . . . . . 81.4 Multiobjective Optimization . . . . . . . . . . . . . . . . . . . 91.5 Thesis outline . . . . . . . . . . . . . . . . . . . . . . . . . . . 10Chapter 2: Road alignment optimization . . . . . . . . . . . . 112.1 Road design system . . . . . . . . . . . . . . . . . . . . . . . 112.2 Horizontal alignment design . . . . . . . . . . . . . . . . . . . 112.2.1 Mathematical model of the horizontal alignment ge-ometry . . . . . . . . . . . . . . . . . . . . . . . . . . . 132.2.2 Horizontal alignment geometric constraints . . . . . . 162.3 Vertical alignment design . . . . . . . . . . . . . . . . . . . . 172.3.1 Vertical alignment geometric constraint . . . . . . . . 21iiiTABLE OF CONTENTS2.3.2 Other constraints . . . . . . . . . . . . . . . . . . . . . 21Chapter 3: Surrogate model for the cost functions . . . . . . 223.1 Surrogate cost model . . . . . . . . . . . . . . . . . . . . . . 233.2 Surrogate cost model for tangential road segment . . . . . . 243.2.1 Computing transition points . . . . . . . . . . . . . . 253.2.2 Length of the tangent road section . . . . . . . . . . . 273.2.3 Volume of ground cut . . . . . . . . . . . . . . . . . . 283.2.4 Volume of ground fill . . . . . . . . . . . . . . . . . . 303.2.5 Procedure for classifying a grid cell as a region of cut,fill, both cut and fill, or neither. . . . . . . . . . . . . 313.2.6 Cost calculation for the tangent road segment . . . . . 313.3 Surrogate cost model for circular road segment . . . . . . . . 323.3.1 Calculation of x-boundary crosses . . . . . . . . . . . 343.3.2 Length of circular road section . . . . . . . . . . . . . 363.3.3 Volume of ground cut . . . . . . . . . . . . . . . . . . 363.3.4 Volume of ground fill . . . . . . . . . . . . . . . . . . 373.3.5 Cost calculation for a circular road segment . . . . . 383.4 The overall surrogate cost function . . . . . . . . . . . . . . . 38Chapter 4: Multiobjective optimization model for 3D roadalignment . . . . . . . . . . . . . . . . . . . . . . . . 404.1 Variable definition . . . . . . . . . . . . . . . . . . . . . . . . 404.2 Problem formulation . . . . . . . . . . . . . . . . . . . . . . . 414.3 NOMAD: Nonlinear optimization with the MADS algorithm 414.4 Solution procedure . . . . . . . . . . . . . . . . . . . . . . . . 42Chapter 5: Numerical results . . . . . . . . . . . . . . . . . . . 435.1 Experimental setup . . . . . . . . . . . . . . . . . . . . . . . . 455.2 Case study 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . 455.2.1 Experimental results 1 . . . . . . . . . . . . . . . . . . 465.3 Case study 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . 485.3.1 Experimental results 2 . . . . . . . . . . . . . . . . . . 505.4 Implication of the results . . . . . . . . . . . . . . . . . . . . 52Chapter 6: Conclusion and Future work . . . . . . . . . . . . . 536.1 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54ivList of TablesTable 5.1 Design constants and parameter values. . . . . . . . . 43Table 5.2 Numerical results for case study 1. . . . . . . . . . . . 49Table 5.3 Numerical results for case study 2 when θ = 0.5. . . . 50Table 5.4 Numerical results for case study 2 when θ = 0. . . . . 52vList of FiguresFigure 1.1 An example a horizontal alignment. . . . . . . . . . 3Figure 1.2 An example vertical alignment, longitudinal profile. . 4Figure 1.3 An example 3D road alignment. . . . . . . . . . . . . 5Figure 2.1 An example horizontal alignment for the model . . . 12Figure 2.2 Horizontal alignment geometry . . . . . . . . . . . . . 14Figure 2.3 Example of horizontal alignment obtained using for-mulas (2.1), (2.2), and (2.3) . . . . . . . . . . . . . . 16Figure 2.4 An example section of vertical alignment. . . . . . . . 18Figure 2.5 Vertical alignment in hz-plane . . . . . . . . . . . . . 20Figure 3.1 An example, the projection of a tangent road segmentonto horizontal plane . . . . . . . . . . . . . . . . . . 26Figure 3.2 An example cut cross-section. . . . . . . . . . . . . . 28Figure 3.3 An example, fill cross-section. . . . . . . . . . . . . . 30Figure 3.4 An example of a pave only cross-section. . . . . . . . 31Figure 3.5 Example of a horizontal curve section . . . . . . . . . 34Figure 3.6 An example of road segments . . . . . . . . . . . . . . 39Figure 5.1 Three-dimensional view of the terrain. . . . . . . . . . 44Figure 5.2 The contour map. . . . . . . . . . . . . . . . . . . . . 44Figure 5.3 Initial horizontal alignment in case study 1. . . . . . 46Figure 5.4 The optimal horizontal alignments for case study 1. . 47Figure 5.5 The optimal vertical alignments for θ = 1. . . . . . . 48Figure 5.6 The optimal vertical alignments for θ = 0. . . . . . . 48Figure 5.7 The horizontal alignments for case study 2. . . . . . 49Figure 5.8 Two optimal horizontal alignments computed fromdifferent starting alignments. . . . . . . . . . . . . . . 50Figure 5.9 The optimal horizontal alignments. . . . . . . . . . . 51Figure 5.10 The optimal horizontal alignments. . . . . . . . . . . 51viAcknowledgementsI would like to express my sincere gratitude to my supervisors Dr. War-ren Hare and Dr. Solomon Tesfamariam for their continuous follow-up andunreserved support towards the completion of my thesis. Their guidanceand encouragements gave me the strength to persevere until the end. I, es-pecially, thank Dr. Hare for his critical comments, inspiration, and patience.I also thank Dr. Yves Lucet who is my graduate committee member. I amvery grateful for his guidance and valuable discussions during our weeklymeetings.I extend my thanks to the funding sources. This research is partiallyfunded by a Collaborative Research and Development (CRD) Grant fromthe Natural Sciences and Engineering Research Council (NSERC) sponsoredby Softree Technical System Inc., International scholarships from the Univer-sity of British Columbia, and by the NSERC Discovery Grants of Dr. Hareand Dr. Tesfamariam. The research was performed in the Computer-AidedConvex Analysis laboratory funded by the Canadian Foundation for Inno-vation (CFI) and by the British Columbia Knowledge Development Fund(BCKDF).I am grateful to Matiyas Ayalew for the good times we spent togetherat home. You are more than just a friend. I also thank all my colleagues formaking my graduate life truly enjoyable; particularly, Sukanto Mondal forhis valuable help with programming.My special thank is to my beloved wife, Genet Alemayehu. She mademe wonder at how people can be that patient. She sacrificed her emotionsand time just to make me feel comfortable and happy with my study. MyParents have a special place in my heart, thanks to them.Above all, I thank God for everything He provides me.viiDedicationThis thesis is dedicated to my wife, Genet Alemayehu and my son, SolenDessalegn.viiiChapter 1Introduction1.1 Background and motivationA good transport network is important in sustaining economic success,as it links people to jobs and delivers products to markets. It is also criticalto domestic and cross-border trades and tourism activities. Road trans-port continues to be the dominant mode of transportation for moving goodsbetween Canada and the U.S.A. [Rev]. In 2011, 56.5% of overall Canada-U.S.A. trade was transported by truck. These benefits come with cost re-lated to road construction. In road construction projects, the main factorsthat determine the location of a road are the construction and maintenancecosts, environmental requirements, and driver’s safety. The cost items tobe considered in the economic analysis of road projects are generally di-vided into two broad categories [JS04]: (1) supplier costs, which are directlyincurred by road (highway) administrators, and (2) user costs incurred byroad users. The supplier costs are further divided in three categories: (a)pavement, and other costs primarily depending upon the length of the align-ment, (b) right-of-way costs including those costs associated with land andenvironmental impacts, and (c) earthwork costs. The user costs are dividedinto three categories: (a) travel-time cost, (b) vehicle-operating cost, and(c) accident cost [CGF89, JS04, SJK06].Road design is a process for the determination of three-dimensional routelocation on the ground surface. It aims to connect two terminals (the startand end points) at minimum possible cost subject to the design, environ-mental, and social constraints [JM07].Since the number of alternative routes connecting two end points isunlimited, a traditional route location analysis, which has relied heavilyon human judgement and intuition, may overlook many good alternative[CGF89, LPZL13]. Thus, in order to reduce the workload for engineers, andto consider all possible route alternatives, an automated procedure that de-termines the road alignments and calculates the construction costs has beendeveloped. The automation of road design problem reduces the tedious anderror-prone manual tasks, most notably drafting [KJLS04]. In addition, this11.1. Background and motivationprocedure allows the use of optimization techniques in search of a good align-ment (an alignment with minimum cost) [OEC73]. Optimization techniquessave much of designer’s time and provide the decision maker with a powerfultool that searches an alignment with minimum cost from a large number ofalternative alignments. In fact, optimization of road alignment can yieldconsiderable savings in construction costs when compared with unoptimizeddesign procedures.The problem of road design can be broken down into three interconnectedstages: horizontal alignment, vertical alignment, and earthwork [HKL11].Optimization of a road alignment is often defined as the minimization of thecost of the project satisfying certain requirements. Optimization of roaddesign therefore implies the search for the optimum location of the verticaland horizontal alignments. Each stage is discussed further below.1.1.1 The horizontal alignmentThe Horizontal alignment is the projection of the three-dimensional roadon the horizontal plane. A typical horizontal alignment is made up of a se-quence of straights, circular curves, and transition curves. Transition curveshave the property that the radius of curvature changes progressively alongthem. The design of a horizontal alignment mainly involves the design ofhorizontal curves and tangents. Curves are provided for smooth transitionto avoid sudden changes in direction. The main considerations in horizon-tal alignment design are that, it should avoid lands which are expensive topurchase or restricted, obstacles which present engineering difficulties, andground which may involve large earthwork or structure costs. Therefore, thecost of road construction for the horizontal alignment problem depends onthe cost of acquiring land and on the output of the vertical alignment stage[HKL11].The horizontal alignment optimization seeks a low cost route while ad-hering to the design standards and reducing environmental impacts [AH11].In the literature, the following models have been developed for optimizinghorizontal alignments: calculus of variation [Nic73], network optimization[Tri87], dynamic programming [OEC73], and genetic algorithms [JJS00].Detailed discussion on the advantages and disadvantage of these methodscan be found in [SJK06].In optimization modelling using calculus of variation, the two end pointsare connected by a curve and integration of a cost function is minimized.This method can generate a smooth alignment and possibly yields a globaloptimum, but requires a continuously differentiable cost function, which21.1. Background and motivationmay not exist over different land-use patterns and geographic features suchas rivers and lakes [JJS00, SJK06]. In Network-optimization approaches, thealignment is represented by the arcs connecting the start and end points andthe modified alignment problem have well-developed solution algorithms,but the resulting alignments are not smooth [JJS00]. The principal assump-tion in dynamic programming approach is that a problem can be divided into a number of sub-problems (or stages) and that the contribution to the ob-jective function value from each sub-problem are independent and additive.For solving the horizontal alignment optimization problem, the stages aredefined to be the evenly spaced lines perpendicular to the axis connectingthe start and end points of the alignment [SJK06]. This approach may needless computer memory but also cannot yield a smooth alignment [JJS00]. Agenetic algorithm integrated with the geographic information systems (GIS)can be used to solve the horizontal road alignment between two given endpoints. First, the model can generate smooth alignments based on highwaydesign standards, and then the genetic algorithm approach can optimizevery complex cost functions including user costs, which have been ignoredin many existing models. Finally, the model directly exploits the informationin a GIS database, which reduces data preprocessing time and allows themodel to search through realistic and highly irregular spatial data [JJS00].An example of a horizontal road alignment section composed of straightsand circular curves is shown in Figure 1.3.Figure 1.1: An example a horizontal alignment.31.1. Background and motivation1.1.2 The vertical alignmentThe vertical alignment is the view of the centreline of the road whenseen along the longitudinal cross section of the road. Generally, the verti-cal alignment is composed of straight sections known as vertical tangentsor grades, and parabolic vertical curves, namely crest and sag curves, seeFigure 1.2. Usually, the vertical alignment design is based on a pre-selectedhorizontal alignment. The main purpose of the vertical alignment is to pro-vide point of elevations along the centreline of the horizontal alignment. Thedesign requirements and other constraints, many of which are conflicting innature, have made the design of a vertical alignment a rather complicatedproblem. Nonetheless, determination of the vertical alignment is a crucialstep in the road design problem since it has important implications on roadconstruction costs, traffic operations, vehicle fuel consumption, and safety[FCS02]. The need for optimization analysis in the selection of a desirablevertical alignment has long been recognized [FCS02].Figure 1.2: An example vertical alignment, longitudinal profile.In the vertical alignment optimization, one fits a road profile to theground profile while respecting various grade constraints and other roadspecifications. The objective is to minimize the cost of construction and thenegative impacts on the environment such as the natural landforms and soil[FCS02]. Many models are found for optimizing vertical alignments than forhorizontal alignments [SJK06]. The vertical alignment optimization modelsinclude linear programming [MS81, Eas88], Numerical search [Hay70], stateparametrization [GCF88], dynamic programming [GAA09, GCF88], genetic41.1. Background and motivationalgorithm [FCS02, AH11], and mixed integer linear programming [HKL11].See also [SJK06, GCF88, GAA09] for more references.1.1.3 The three-dimensional alignmentA three-dimensional road alignment is the superimposition of two di-mensional horizontal and vertical alignments. In essence, road alignmentdesign is a three-dimensioanl problem represented in X,Y , and Z coordi-nates. Figure 1.3 represents an example of a three-dimensional road.The development of models that fully optimize a three-dimensional roadalignment problem is not yet successful, because there are more factors in-volved and more complexities in the geometric specification of the alignment[SJK06, GAA09, ASAC05]. There are two basic approaches found in theliterature: models that simultaneously optimize the horizontal and verticalalignments [ASAC05, JS03, CGF89], and models that employ two or morestages of optimization to get an optimized three-dimensional road align-ment [Nic73, Par77]. Existing approaches for three-dimensional alignmentoptimization include genetic algorithms [JS03, KJSK07, KJS12], dynamicprogramming [LPZL13], neighbourhood search [CL06], and distance trans-form [DS06].Figure 1.3: An example 3D road alignment.Due to the rigorous geometric specifications of alignments and the com-51.1. Background and motivationplexity of the problem, none of the existing methods have completely solvedthe problem of alignments optimization [LPZL13]. The model developed in[CGF89] yields the simultaneously optimized three-dimensional road align-ment optimization problem. It uses a series of cubic spline functions toparametrize the road alignment, and the problem is formulated as a calcu-lus of variations problem. The method of constraint transcription is usedin optimal control theory is employed to transform constraints into a one-dimensional problem.A dynamic programming model is used in [Nic73] to optimize the three-dimensional road alignment problem in two stages. In this approach, atthe first stage, the model searches a relatively coarse grid of points for apreliminary alignment. Then a discrete variational calculus is adopted torefine the alignment so that the resulting alignment can deviate from thegrid points. The resulting solution is rough due to the storage requirementfor searching the initial grid.1.1.4 Surrogate modelsEngineering design heavily depends on computer simulations. The appli-cation of optimization to such systems, therefore, has to do with the objec-tive functions that may come from large scale computer operations or com-puter simulations. Unfortunately, accurate and high-fidelity simulations areoften computationally expensive, with evaluation times ranging from hoursto days [KL13]. In many cases, optimization of such objectives by direct ap-plication of optimization routines is impractical [KY11]. One reason is thatconventional optimization algorithms require tens, hundreds or even thou-sands of objective function calls per run, which makes the computational costof the whole optimization process intractable [KY11, WPR05, QHS+05]. Asecond reason is that simulation-based objective functions are often analyt-ically intractable, discontinuous, non-differentiable, noisy, or possibly fail toreturn a value [KY11]. Therefore, any technique that reduces the functionevaluation count is crucially important. Feasible handling of these objectivefunctions can be accomplished using surrogate models [KY11]. A surrogatesf of the function f is a function that shares some similarities with f , butis much cheaper to evaluate. A surrogate model is typically less accurateor has less quality than the true model, but it is cheaper to evaluate orconsumes fewer computing resources.For optimization problems, surrogate models can be regarded as approx-imation models for the objective function and the constraint set. The basicconcept of surrogate-based optimization is that the direct optimization of61.2. Problem statementthe computationally expensive model is replaced by an iterative process thatinvolves the creation, optimization, and updating of a fast and analyticallytractable surrogate model [KY11]. The surrogate should be a reasonablyaccurate representation of the true model, at least locally. The solutionobtained by optimizing the surrogate model is verified by evaluating thesolution of the true model which is used to update the surrogate modeliteratively until some termination criterion is met [KY11].The surrogate model generation is a key component of any surrogate-based optimization algorithm. There are two ways of generating the surro-gate model [KY11]a) a model constructed from physically-based low-fidelity models, andb) a model based on the function approximation of sampled data.Since the physical surrogate is based on a particular knowledge of the physi-cal system of interest, it can be formulated using analytical or semi-empiricalformulas. Once the surrogate model is built, an optimization algorithms canbe used to yield an approximation of the minimizer (or maximizer) of thetrue model.1.2 Problem statementRoad alignment optimization based on cost minimization requires com-prehensive formulation of costs and development of efficient solution algo-rithms. Among the major cost components that contributed in the con-struction of a road are the earthwork cost and pavement costs [CGF89].In this thesis, we develop a surrogate model for the cost of groundcut, ground fill, material waste, and pavement based on the ap-proximation of the physical terrain by planes.These cost items account for more than 50% of the construction cost [CGF89].Although the ground cut and fill costs depend on the varying nature of thesoil, we, regard it as a function of the volume of earth to be cut or filled.In the model, we have considered the volume of wasted material as thedifference between the volume of ground cut and ground fill. The average-end-area method is widely used to estimate the volume of earthwork betweentwo stations, however, it gives less accurate estimation when the terrain isirregular or mountainous [SJK06]. The pavement cost is computed usingroad length.71.3. Research approach1.3 Research approachThis thesis seeks to solve a three-dimensional road alignment optimiza-tion problem. This problem is presented as a complex optimization prob-lem in which both vertical and horizontal alignments are determined. Aswe have already discussed above, this problem can be solved in two ap-proaches. The first approach solves both vertical and horizontal alignmentssimultaneously, and the second approach uses multi-stage (or multi-level)based optimization methods. Both approaches, however, suffer from thehigh computational effort and large memory requirements [SJK06]. Three-dimensional road alignment optimization is a problem having a constrained,nonlinear, and non-differentiable structure that cannot be efficiently solvedby classic optimization techniques [ASAC05]. As a result, many heuristicbased methods, such as genetic algorithm, simulated annealing, and tabusearch, have been used to solve the problem.In our research, we modelled the problem as a simultaneous optimizationof vertical and horizontal alignment problems. In this approach the verticalalignment is composed of straight lines and the horizontal alignment is com-posed of straight lines and circular curves. The cost component consideredare the earthwork cost. A surrogate model that calculates the cost of groundcut, ground fill, ground waste, and paving is developed. The design vari-ables corresponding to the vertical and horizontal alignments are reflectedin the surrogate cost model formulation. The surrogate model calculates thecost quickly, but provides a less accurate earthwork cost than other mixedinteger linear programming (MILP) techniques.Furthermore, in mountainous regions the cost function for ground cutand fill tends to favor a circuitous, perhaps the longest distance, alignments,while the paving cost function favors the shortest distance alignment. Asa result, the two cost components conflict with each other. Therefore, theproblem is presented as a multiobjective optimization problem. This prob-lem is a nonlinear, nonconvex, non-differentiable constrained optimizationproblem. Consequently, the derivative-free optimization method, known asNOMAD (Nonlinear optimization with the MADS algorithm) is used tosolve the problem. By minimizing the cost function, NOMAD generates theradius of curvature and coordinates of the intersection points for the hori-zontal alignment, and the coordinates of vertical intersecting points for thevertical alignments.81.4. Multiobjective Optimization1.4 Multiobjective OptimizationIn real life, decision makers or engineers are frequently faced with deci-sion problems of several, perhaps conflicting, objectives. It is only naturalto want all of the choices or decisions to be as good as possible, in otherwords, optimal. Single objective optimization technology is not sufficientto deal with problems with conflicting objectives. Thus, multiobjective opti-mization is an appealing alternative. In multiobjective optimization, severalconflicting objective functions have to be optimized simultaneously over afeasible set determined by constraint functions. There is, usually, no uniquesolution that is simultaneously optimal for all objectives. As a result, onecan only consider a trade-off among the objectives, and the primary goalof multiobjective optimization is to seek the best trade-off to support thedecision maker in choosing a final preferred solution.Unlike the single objective optimization, in the multiobjective context,the interest is often more in the objective space. For this reason, the no-tion of optimality depends on how decision alternatives are compared andordered, that is, on the order relation in the objective space [Mie99]. Theorder relations in objective space can be defined by different dominance rela-tions such as Pareto dominance, Geoffrions dominance, and Lexicographicaldominance, of which the Pareto dominance is most commonly used [ASZ08].Although there is no universally accepted solution concept for decision prob-lems with multiple objectives, one would agree that a good solution mustnot be dominated by the other feasible alternatives [Yu74]. Generally, mul-tiobjective optimization problems (MOP) solution techniques are classifieddepending on the moment when the engineer or decision maker is able toestablish preferences relating the different objectives [CMVV11]. The pref-erence informations may be classified as: no preference, a priori preference,progressive preference, and a posteriori preference.Returning to horizontal alignment, in road design we have two conflictingobjectives: (1) the cost that depends on the length of a route and (2) the costthat depends on the volume of earthwork. Different types of costs will fa-vor different alignment configurations. For example, length-dependent costsand user costs tend to straighten the alignment, while location-dependentcosts tend to favor more indirect and circuitous alignment [SJK06]. In thissituation, there is a conflict between length of an alignment and the cost ofan alignment. In this thesis, the surrogate cost model is formulated for twocosts items, the cost that depends on the length of an alignment (pavingcost) and the cost that depends on the volume of earthwork (cost due toground cut, ground fill, and wasted material). Thus, a three-dimensional91.5. Thesis outlinealignment optimization problem is modelled as a multiobjective optimiza-tion problem.1.5 Thesis outlineIn Chapter 2, we present the mathematical models for the horizontaland vertical alignment geometries. The problem variables and constraintscorresponding to the horizontal alignment and vertical alignment are alsospecified. Chapter 3 deals with the formulation of the surrogate model forthe cost functions. The surrogate model corresponding to the tangent andcircular sections of the horizontal alignment are presented. The model cor-responding to the circular curve is more challenging. In Chapter 4, themultiobjective optimization model and NOMAD are discussed. In particu-lar, a solution procedure for the three-dimensional optimization problem ispresented. The numerical results are given in Chapter 5. Some concludingremarks are given in Chapter 6.10Chapter 2Road alignment optimizationThis chapter presents the mathematical formulation of the geometric de-sign for horizontal and vertical road alignments. The geometric constraintsrelated to each alignment is given. The mathematical formulation entails thespecification of the design elements that also reflect the cost of the individualalignment.2.1 Road design systemIn the road alignment development process, the first stage (the planningstage) is the selection of a corridor along which the road is to pass. Followingthis stage is the road alignment design stage, which deals with the location,alignment, and shape of a road. The selection of a specific alignment involvesthe determination of horizontal and vertical alignments of a road profile,subject to a set of constraints and requirements [CGF89].2.2 Horizontal alignment designIn designing a horizontal alignment geometry, a series of tangents andcurved sections are joined. Restrictions on the horizontal road alignment ge-ometry mainly come from code requirements and external limitations. Thesedesign codes require that the horizontal alignment of a road be composedof three types of design elements, namely line segments, circular arcs, andtransition curves. The circular and transition curves are typically combinedto form the curved section. The transition curves are often applied betweentangents and circular curves for reducing a sudden change in the curvature.In our model, the horizontal alignment is composed of the tangentialsegments and circular curves. In this approach the horizontal alignmentgeometry is required to satisfy two criteria: (1) the alignment should satisfythe orientation specified as tangent−circle− tangent, and (2) the start andend sections of the alignment should be tangent segment. The circular curvesare placed between two adjacent tangents to mitigate the sudden changes in112.2. Horizontal alignment designthe direction of the alignment. In this way the absence of transition curvesdoes not have a major effect on the applicability of the model. Thus, ahorizontal alignment for the model consists of two types of road segments:straight lines (tangents) and circular curves, see Figure 2.1.Figure 2.1: An example horizontal alignment for the modelThe exact shape of the horizontal alignment is determined by the set ofintersection points and the radius of a circular curve. The circular curvesare inserted between tangents as defined by the intersection points. Thereare three decision variables associated with each intersection point, namelythe x, y, and r, where x and y are the x−coordinate and y−coordinatesof intersection point, and r is the radius of the circular curve. Thus, if thenumber of intersection points is N , the horizontal alignment has to select 3Nvariables. The most important horizontal alignment constraints, which arerelated to the geometric design and code requirements, are given as follows:1. two adjacent circular segments should not overlap;2. the radius of each circular curve should be greater than the minimumturning radius (i.e. rmin ≤ r );3. the alignment should stay in the corridor (or design space).122.2. Horizontal alignment design2.2.1 Mathematical model of the horizontal alignmentgeometryWe design horizontal alignment as a curve that is described by a series ofintersection points, IP , and radius of curvature r that determine the designelements (the tangents and the circular curves).Determination of horizontal alignment design elementsWe assume we are given N intersection points of the horizontal alignment{IP1, IP2, · · · , IPN},and the corresponding set of radius of curvature{r1, r2, · · · , rN},between the start and end points Sh and Eh, where Sh and Eh are thehorizontal component of the start and end points of the three-dimensionalalignment. Then one can insert N + 1 tangential road segments and Ncircular curves to get the exact shape of a horizontal alignment. The hori-zontal alignment geometry is determined by design elements, the tangentsand circular curves. In this section, we calculate the points that determinethe tangent section and circular curves, such as the center of the curve (Ck)and transition points TCk (from tangential section to circular section) andCTk (from circular section to tangential), see Figure 2.2.We denote the vector from IPk to IPk−1 and IPk to IPk+1, respectively,by IPk(k−1) and IPk(k+1). Given intersection points IPk−1, IPk and IPk+1,and radius rk, let CIRk be the circle of radius rk that is tangent to IPk(k−1)and tangent to IPk(k+1) at TCk and CTk, respectively. We denote thecentre of CIRk by Ck = (xck , yck), we denote the points TCk = (xtck , ytck)and CTk = (xctk , yctk). Let θk be the angle between IPk(k−1) and IPk(k+1),then we calculate the coordinates of Ck, TCk, and CTk.132.2. Horizontal alignment designFigure 2.2: Horizontal alignment geometryCalculation of the transition points TCk and CTkFrom the dot product of vectors IPk(k−1) and IPk(k+1), we haveθk = arccos(IPk(k−1) · IPk(k+1)‖IPk(k−1)‖‖IPk(k+1)‖), (2.1)whereIPk(k−1) = IPk−1 − IPk = (xk−1 − xk, yk−1 − yk)IPk(k+1) = IPk+1 − IPk = (xk+1 − xk, yk+1 − yk).Also, we denoteDP = IPk(k−1) · IPk(k+1) and NP = ‖IPk(k−1)‖‖IPk(k+1)‖,where DP and NP means dot product and norm product, respectively. LetLk be the distances from IPk to TCk. Thentanβk2=Lkrk⇒ Lk = rk tanβk2,where βk is the central angle of the circular arc, see Figure Horizontal alignment designUsing the identitytanβk2=1− cosβksinβkand the fact that βk = pi − θk,we haveLk = rk1 + cos θksin θk.Notice that, because of the congruence of 4CkTCkIPk and 4CkCTkIPk,the length of line segments IPkTCk and IPkCTk are equal to Lk.Let Vk and Uk denote the vectors directed from IPk to TCk and IPkto CTk, respectively. That is,Vk = TCk − IPk and Uk = CTk − IPk.ThenVk = TCk − IPk = LkIPk(k−1)‖IPk(k−1)‖, and Uk = CTk − IPk = LkIPk(k+1)‖IPk(k+1)‖.Thus, we getTCk = IPk + LkIPk(k−1)‖IPk(k−1)‖, and CTk = IPk + LkIPk(k+1)‖IPk(k+1)‖. (2.2)Calculation of the center: CkLet Mk be the midpoint of the line segment that connects TCk and CTk.Clearly Mk lies on the line that passes through IPk and Ck. ThenMk =12(TCk + CTk)= IPk +Lk2(IPk(k−1)‖IPk(k−1)‖+IPk(k+1)‖IPk(k+1)‖).Now,Ck − IPk = (rk + dk)Mk − IPk‖Mk − IPk‖,where dk is the distance from IPk to the circle CIRk; calculated assecβk2=rk + dkrk⇒ dk = rk secβk2− rk = rk cscθk2− rk.Therefore,Ck = IPk + rk cscθk2(Mk − IPk‖Mk − IPk‖). (2.3)152.2. Horizontal alignment designExample 1. Given the following set of intersection points, radius of curva-ture, and the start and end points.Sh = (5, 2), Eh = (17, 17), IP1 = (10, 7), IP2 = (14, 5), IP3 = (16, 10), IP4 =(13, 13), r1 = 1.8, r2 = 2, r3 = 2, r4 = 1.8. Then, the horizontal design ele-ments (tangent section and circular arcs) that are generated using the aboveformulations are depicted in Figure 2.3.Figure 2.3: Example of horizontal alignment obtained using formulas(2.1), (2.2), and (2.3)2.2.2 Horizontal alignment geometric constraintsThe design of horizontal alignment geometry is restricted or constrained,mainly, by code requirements and external limitations such as control ar-eas or restricted areas. The set of constraints for horizontal alignment arediscussed below.Box constraints on the IP locationIn our model, the location of each intersection point is contained by abox. Let IPk = (xk, yk), and let xuk , xlk , yuk , and ylk be real numbers. Thenthe box constraint corresponding to IPk is given as follows.xlk ≤ xk ≤ xuk , and ylk ≤ yk ≤ yuk , k = 1, 2, · · · , N.162.3. Vertical alignment designThe circular curves should not overlapThis constraint depends on the length of the tangent section between twoadjacent circular curves. Two adjacent curves can meet only if the length ofthe tangent between them is zero. We note that, because of the requirementsdiscussed in Section 2.2 the length of the first and last tangential segmentcannot be zero. Thus, this constraint can be written, mathematically, asfollows.0 ≤ ‖TCk − IPk−1‖ − ‖CTk−1 − IPk−1‖, k = 1, 2, · · · , N. (2.4)The minimum turning radiusGiven the minimum radius (rmin) of the circular curve, the optimal radiusat each intersection point has to satisfy the minimal radius requirement set.Thusrmin ≤ rk, k = 1, 2, · · · , N. (2.5)2.3 Vertical alignment designVertical alignment is the view of three dimensional road alignment whenseen along the longitudinal cross section of the road. It is composed ofstraight lines known as vertical tangents and parabolic vertical curves. Thereare two forms of vertical curves: the crest vertical curve and sag verticalcurve. Curves are used to provide a gradual change in elevation betweensuccessive tangents for smooth traverse. Therefore, the standard approachfor designing vertical alignment is by the selection of proper grades for thetangent sections and proper curve lengths. However, since the length ofvertical curve relative to the length of tangent section is small, in our model,the design of vertical alignment is composed of only vertical tangents. Atypical section of vertical alignment is shown in Figure 2.4.In this section, we present the mathematical model of vertical alignmentdesign element, the tangent section. The vertical alignment is modelled asa continuous piecewise linear function. We begin the model by stretchinga surface orthogonal to the xy−plane along the horizontal alignment. Thesurface is so stretched to be flat, so we call it the hz−plane, where h is thedistance measured along the horizontal alignment. The projection onto thehz−plane of the three-dimensional road alignment is the vertical alignment.172.3. Vertical alignment design𝑙𝑒𝑛𝑔𝑡ℎ 𝑜𝑓 𝑎 𝑟𝑜𝑎𝑑 𝑒𝑙𝑒𝑣𝑎𝑡𝑖𝑜𝑛 𝑜𝑓 𝑎 𝑟𝑜𝑎𝑑 𝑣𝑒𝑟𝑡𝑖𝑐𝑎𝑙 𝑡𝑎𝑛𝑔𝑒𝑛𝑡𝑠 Figure 2.4: An example section of vertical alignment.Given the following sets of intersection points and the horizontal tangentpoints, calculated in Section 2.2.1,IP = {IP1, IP2, · · · , IPN}, (2.6)TC = {TC1, TC2, · · · , TCN}, (2.7)CT = {CT1, CT2, · · · , CTN}. (2.8)We partition the horizontal alignment between two adjacent tangent points(CTk−1 and TCk)into mk equally spaced points. Then, correspondingto these points, we design a set of vertical points whose projection onto the horizontal plane is required to stay on the line segment between(CTk−1 and TCk), see Figure 2.5. The sole criterion for determining thesevertical points is based on the minimization of the cost of earth cut, earthfill, earth balance, and the cost of paving. We define the set of the verti-cal points corresponding to horizontal tangent between CTk−1 and TCk as182.3. Vertical alignment designfollows:VPk ={V P1, V P2, · · · , V Pmk}, (2.9)where V Pj = (xkj , ykj , zkj ). The mathematical details of the design roadelevation zkj , and the xy−coordinates xkj and ykj , is given below.Note that since the length of the curved section of the horizontal align-ment is typically much smaller than the tangent section, we do not havevertical point corresponding to the circular curve.Next, we calculate the coordinates of the vertical point V Pj . Giventhe horizontal tangent points CTk−1 and TCk. Partition the line segmentbetween CTk−1 and TCk into mk equally spaced points, and let pj−1 =(xkj−1 , ykj−1) and pj = (xkj , ykj ) be consecutive points.For j = 1, 2, · · · ,mkxkj = xctk−1 +jmk(xtck − xctk−1), (2.10a)ykj = yctk−1 +jmk(ytck − yctk−1), (2.10b)xkj−1 = xctk−1 +j − 1mk(xtck − xctk−1), (2.10c)ykj−1 = yctk−1 +j − 1mk(ytck − yctk−1). (2.10d)Let zkj−1 and zkj be the design road elevations corresponding to pj−1 andpj , respectively. Now we parametrize the line segment between pj−1 and pjusing the parameter s as follows. For j = 1, 2, · · · ,mkj − 1mk≤ s ≤jmk.Then, the coordinates of a point p = (x, y, z) between V Pkj−1 and V Pkj areparametrized as follows,x(s) = xkj−1 + (xkj − xkj−1)(mks+ 1− j) = xctk−1 + s(xtck − xctk−1),(2.11a)y(s) = ykj−1 + (ykj − ykj−1)(mks+ 1− j) = yctk−1 + s(ytck − yctk−1),(2.11b)z(s) = zkj−1 + (zkj − zkj−1)(mks+ 1− j). (2.11c)192.3. Vertical alignment designFigure 2.5: Vertical alignment in hz-planeThe horizontal distance between V Pkj−1 and V Pkj , denoted as dkj isgiven bydkj =√(xkj − xkj−1)2 + (ykj − ykj−1)2.The calculation of the parameter s is discussed in the next chapter.The number of vertical alignment variables, and the precision of earth-work volume estimation depends on the magnitudes of mk. The number ofvertical alignment variables on each horizontal tangent section, except forthe first and last horizontal tangents, is mk+1. The first and last horizontal202.3. Vertical alignment designtangent section each have m0 and mN variables, respectively, because thestart and the end of a road are not variables.Suppose the number of horizontal intersection points is N . Then, thenumber of variables of the vertical alignment isM = m0 +mN +N−1∑k=1(mk + 1).2.3.1 Vertical alignment geometric constraintThe vertical tangent section is constrained by the standard design lim-its. Thus, the design of vertical alignment is constrained by the maximumgradient requirement and the vertical elevation boundaries.The maximum allowable slope (or gradient) of the roadOnce the vertical alignment variables are set, the design constraints haveto be checked before conducting cost calculations. In compliance with thedesign standards, the following requirements have to be satisfied in the de-sign of vertical alignment. Therefore, the vertical alignment design variableis required to respect the maximum grade constraint∣∣∣∣zkj − zkj−1dkj∣∣∣∣ ≤ Gmax.The maximum and minimum elevation of the roadWhile not strictly necessary, it is reasonable to place a maximum andminimum road elevation constraints. Given the vertical off-set z¯ from thecurrent ground elevation zg,zg − z¯ ≤ zkj ≤ zg + z¯, for all kj .2.3.2 Other constraintsIn addition to the horizontal and vertical alignments geometric con-straints, we require the circular section of the road to have the same el-evation (i.e. flat). This requirement is a constraint for the optimizationproblem which is given mathematically as follow,zkmk = z(k+1)1 , k = 1, 2, · · · , N.21Chapter 3Surrogate model for the costfunctionsIn this chapter we provide the mathematical formulation of the surrogatecost model for the three-dimensional road design problem. The surrogatemodel is developed based on the approximation of the ground surface by aseries of planar surfaces. The designed road is represented by a continuouspiecewise linear function which is expressed as parametric equations of aparameter t, where t is normalized to lie between 0 and 1. The horizontalplane in the region of interest is divided into grid cells of equal size, whichare small enough that the terrain above the grid is approximated by a singlelinear function. A matrix format is employed to store important groundinformation for the region of interest. Three different matrices are used tostore the data for the determination of the elevation of the ground at aspecified point in the xy−plane.The surrogate cost function is formulated as the sum of the costs due tothe paving and the earthwork. The cost associated with earthwork is thecost of ground cut, ground fill, and cost of waste material, which depend onthe volume of the material. The calculation of the volume of ground cut andground fill depends on the elevation difference between the approximatingplane and the design road. Since the volume of earthwork is approximated,the cost calculation is only an approximation of the true cost. The pavingcost is computed based on the length of the road. Hereafter, we refer to theapproximating plane elevation as the ground profile. Ground cutting costsare incurred when the design road profile is lower than the ground profile,while filling costs are calculated when the road profile is above the groundprofile.The three-dimensional road alignment optimization problem is based onthe minimization of a parametrized surrogate cost function subject to align-ment constraints. Penalty parameters for each of the cost components areintroduced. The optimal road alignment obtained by solving the surrogate-based optimization problem can be calibrated to different optimal align-ments by choosing different values of the penalty parameters. The objective223.1. Surrogate cost modelfunction, as the total cost C, may be expressed in a general form as:C = CcVc + CfVf + CpWL+ Cw‖Vf − Vc‖, (3.1)whereL is length of the roadVc is volume of ground cutVf is volume of ground fillW is width of the roadCc is the ground cut penalty parameterCf is the ground fill penalty parameterCp is the paving penalty parameterCw is the wasted material penalty parameter.Earthwork volume estimation is one of the most important componentsin formulating the road construction costs. In our approach, the earthworkvolume is estimated by continuous integration of the cross-sectional areabetween two end points along the road. The target road should minimizethe surrogate cost model and conform to the constraints on horizontal andvertical alignments.3.1 Surrogate cost modelConsider a three-dimensional road alignment between the starting pointSTART = (xs, ys, zs) and the ending point END = (xe, ye, ze). We regardthe projection onto the xy−plane (or horizontal plane) of a three dimensionalalignment as being composed of tangents and circular curves as discussed inChapter 2. We first divide the xy-plane into a square grid. We denote the(u, v)th grid cell by Guv. The terrain in each grid is then approximated bya linear functionz = Auvx+Buvy + Cuv,where x, y ∈ Guv, and Auv, Buv and Cuv are matrices containing the linearfunction data.We aim to calculate the cost of a road section (tangential section orcircular section) which is the sum of costs corresponding to each grid cellalong the section. To approximate the cost of a road segment, one requires233.2. Surrogate cost model for tangential road segmenti) an estimate of the volume of earth cut (Vc),ii) an estimate of the volume of earth filled (Vf ),iii) the total length of the road (L).We then apply the costing parameters Cc, Cf , Cw, Cp to create a parametrizedcost function (3.4). For the collection S of tangential road sections and cir-cular road sections, we define Vcξ is the total volume of ground cut on roadsegment ξ ∈ S, and Vfξ is the total volume of ground fill on road segmentξ ∈ S. The overall volume of cut, denoted Vc and volume of fill, denoted Vf ,is given as− Vc =∑ξ∈S Vcξ , and− Vf =∑ξ∈S Vfξ .The cost function formulations for the tangential segment and circularcurve are given in the following sections.3.2 Surrogate cost model for tangential roadsegmentSuppose ξk ∈ S corresponds to a tangential segment. That is, horizontalprojection is the tangent line that connects CTk−1 and TCk. We denote thishorizontal tangent by HTk. Suppose there are mk grade lines or verticaltangents along HTk. For j = 1, 2, 3, · · · ,mk, the vertical points, denoted byV Pkj and V Pkj−1 , are given asV Pkj−1 = (xkj−1 , ykj−1 , zkj−1), and V Pkj = (xkj , ykj , zkj ),where zkj−1 , zkj , xkj−1 , xkj , ykj−1 , and ykj are calculated using equations equa-tions (2.10a) to (2.10d). Then the parametric representation of a verticaltangent between V Pkj−1 and V Pkj is given asrt(s) =(x(s), y(s), z(s)),wherej − 1mk≤ s ≤jmk, and the parametric equations x(s), y(s), and z(s)are given in equations (2.11a) to (2.11c). Thus, we define the parametricroad profile and the parametric equation of a plane that approximates theterrain in grid Guv aszr(s) = z(s), and zg(s) = Auvx(s) +Buvy(s) + Cuv.243.2. Surrogate cost model for tangential road segmentFor x(s), y(s) ∈ Guv we havezr(s) = zkj−1 + (zkj − zkj−1)(mks+ 1− j), (3.2)andzg(s) = Auvx(s) +Buvy(s) + Cuv= Auv(xctk−1 + s(xtck − xctk−1))+Buv(yctk−1 + s(ytck − yctk−1))+ Cuv=(Auvxctk−1 +Buvyctk−1 + Cuv)+s(Auv(xtck − xctk−1) +Buv(ytck − yctk−1)). (3.3)3.2.1 Computing transition pointsSince it may be required to cut and fill the ground within a particular gridcell, we need to calculate the x and y values at the point of transition fromcut to fill or fill to cut. Next, we define the sets of parameters T jxk , Tjyk , Tjtkfor the x-boundary cross, y-boundary cross, and for the transition point.These parameters are used to calculate the road elevation in (3.2) and theground elevation in (3.3). Let xu, yv−1 be the x and y boundary points atwhich the horizontal tangent HTk crosses the grid cell Guv, see Figure 3.1.Then, forj − 1mk≤ s ≤jmkwe havex(s) = xctk−1 + s(xtck − xctk−1) = xu ⇒ s =xu − xctk−1xtck − xctk−1,y(s) = yctk−1 + s(ytck − yctk−1) = yv−1 ⇒ s =yv−1 − yctk−1ytck − yctk−1.Therefore,T jxk ={s | s =xu − xctk−1xtck − xctk−1,j − 1mk≤ s ≤jmk}T jyk ={s | s =yv−1 − yctk−1ytck − yctk−1,j − 1mk≤ s ≤jmk}.253.2. Surrogate cost model for tangential road segmentFigure 3.1: An example, the projection of a tangent road segment ontohorizontal planeIf there exist a transition from cut to fill (or fill to cut) within Guv, theparameter of transition s is calculated by solving the equation zg(s) = zr(s).Using equations (3.3) and (3.2), we solve for the transition point parameters.zr(s) = zg(s)⇔ zkj−1 + (zkj − zkj−1)(mks+ 1− j) =(Auvxctk−1+Buvyctk−1 + Cuv)+ s(Auv(xtck − xctk−1)+Buv(ytck − yctk−1))⇔(Auvxctk−1 +Buvyctk−1 + Cuv − zkj−1− (zkj − zkj−1)(1− j))= s(mk(zkj − zkj−1)−Auv(xtck − xctk−1)−Buv(ytck − yctk−1))⇔ s =ϕχ,whereϕ =(Auvxctk−1 +Buvyctk−1 − zkj + j(zkj − zkj−1) + Cuv),χ =(mk(zkj − zkj−1)−Auv(xtck − xctk−1)−Buv(ytck − yctk−1)).263.2. Surrogate cost model for tangential road segmentHence, forj − 1mk≤ s ≤jmk, the set T jtk of the transition parameters isgiven byT jtk ={s|s =(Auvxctk−1 +Buvyctk−1 − zkj + j(zkj − zkj−1) + Cuv)(mk(zkj − zkj−1)−Auv(xtck − xctk−1)−Buv(ytck − yctk−1))}.(3.4)We define the union of all sets of parameters asT jk = Tjxk ∪ Tjyk ∪ Tjtk .Let K be the cardinality of T jk . We sort Tjk in an increasing order to createT jks asT jks = {s1 < · · · < sK−1 < sK}.Next, we compute the volume of ground cut Vck , volume of ground fillVfk , and length of the road Ltk .3.2.2 Length of the tangent road sectionGiven the points V Pkj = (xkj , ykj , zkj ) and V Pkj−1 = (xkj−1 , ykj−1 , zkj−1),where xkj , xkj−1 , ykj , and ykj−1 are calculated using equations (2.10a) to (2.10d).The length of rt(s), denoted Ljtk , is given byLjtk =∫ jmkj−1mkdrt(s)=1mk√(xkj − xkj−1)2 + (ykj − ykj−1)2 + (zkj − zkj−1)2=1mk√1m2k(xtck − xctk−1)2 +1m2k(ytck − yctk−1)2 + (zkj − zkj−1)2=1m2k√(xtck − xctk−1)2 + (ytck − yctk−1)2 +m2k(zkj − zkj−1)2. (3.5)The length of tangential road segment is computed asLtk =mk∑j=1Ljtk . (3.6)273.2. Surrogate cost model for tangential road segment3.2.3 Volume of ground cutConsider two consecutive parameters si−1, si ∈ Tjks. If the consecutiveparameters bracket a cut, then the approximate volume of ground cut forthe tangential road segment over [si−1, si] is obtained by integrating thecross-section area over the effective length of the road in the specified gridcell. The elevation difference at s ∈ [si−1, si] is given byhc(s) = zg(s)− zr(s)=((Auvxctk−1 +Buvyctk−1 + Cuv)+ s(Auv(xtck − xctk−1)+Buv(ytck − yctk−1)))−(zkj−1 + (zkj − zkj−1)(mks+ 1− j))=(Auvxctk−1 +Buvyctk−1 − zkj + j(zkj − zkj−1) + Cuv)+s(Auv(xtck − xctk−1) +Buv(ytck − yctk−1)−mk(zkj − zkj−1)).Next, we calculate the cross-sectional area of the ground. We assumethe cross-section of the ground to be cut a trapezoid. An example is shownin Figures 3.2.Figure 3.2: An example cut cross-section.The cross-sectional area of the trapezoid is calculated asac =12h2c cot θ1 +Whc +12h2c cot θ2 = hc(W +12hcκ)(3.7)where W is the width of the road, θ1 and θ2 are side slope angles, andκ = cot θ1 + cot θ2.283.2. Surrogate cost model for tangential road segmentThe values of θ1 and θ2, and the value of hc are inter-dependent. More-over, we assume that 0 < θ1 <pi2and 0 < θ2 <pi2, which implies that0 < κ < ∞. Typically, for the surrogate model, we fix the values of θ1 andθ2. The approximate volume of ground cut over [si−1, si] is given asV jck =∫ sisi−1ac(s)ds, (3.8)where j = 1, 2, 3 · · · ,mk.Assuming si−1 and si bracket a cut, the cross-section area ac(s) is cal-culated as follows.ac(s) = hc(s)(W +12κhc(s))= Whc(s) +12κh2c(s)= W(zg(s)− zr(s))+12κ(zg(s)− zr(s))2= W[(Auvxctk−1 +Buvyctk−1 − zkj + j(zkj − zkj−1) + Cuv)+s(Auv(xtck − xctk−1) +Buv(ytck − yctk−1)−mk(zkj − zkj−1))]+12κ[(Auvxctk−1 +Buvyctk−1 − zkj + j(zkj − zkj−1) + Cuv)+s(Auv(xtck − xctk−1) +Buv(ytck − yctk−1)−mk(zkj − zkj−1))]2= W (Ω + sΘ) +12κ(Ω + sΘ)2=(WΩ +12κΩ2)+(WΘ + κΩΘ)s+12κΘ2s2 (3.9)whereΩ =(Auvxctk−1 +Buvyctk−1 − zkj + j(zkj − zkj−1) + Cuv),Θ =(Auv(xtck − xctk−1) +Buv(ytck − yctk−1)−mk(zkj − zkj−1)).If si−1 and si do not bracket a cut, then ac(s) is defined as 0. Integrating(3.8) using (3.9), we getCutvolume =(WΩ +12κΩ2)(si − si−1) +12(WΘ + κΩΘ)(s2i − s2i−1)+16κΘ2(s3i − s3i−1).293.2. Surrogate cost model for tangential road segmentTherefore, we compute the volume asV jck ={Cutvolume if si−1, si bracket a cut segment0 otherwise(3.10)We can now construct the volume of ground cut, Vck for the tangentialroad segment asVck =mk∑j=1V jck . (3.11)3.2.4 Volume of ground fillThe approximate volume of ground fill for the tangential road segmentover [si−1, si] is obtained by integrating the cross-section area over the effec-tive length of the road in the specified grid cell, an example fill cross-sectionis shown in Figure 3.3. The approximate volume calculation is given asV jfk ={ ∫ sisi−1af (s)ds if si−1, si bracket a fill segment0 otherwise,(3.12)where af (s) = −ac(s) and j = 1, 2, 3 · · · ,mk.Therefore, the ground fill volume approximation is calculated using thesame procedure as ground cut volume calculation.Figure 3.3: An example, fill cross-section.We can now construct the volume of ground fill, Vfk for the tangentialroad segment asVfk =mk∑j=1V jfk . (3.13)303.2. Surrogate cost model for tangential road segment3.2.5 Procedure for classifying a grid cell as a region of cut,fill, both cut and fill, or neither.In order to determine whether or not the ground corresponding to acertain grid cell has to be cut or filled, one can check the following situations.For i = 1, 2, · · · ,K − 1,1. if(zg(si)−zr(si))=(zg(si+1)−zr(si+1))= 0, then there is no groundcut or ground fill required. So V jck = Vjfk= 0;2. if(zg(si) − zr(si))(zg(si+1) − zr(si+1))≥ 0 but not both equal, thenthe grid cell is identified as cut region or fill region. Furthermore,− if(zg(si)− zr(si))> 0, it is cut region,− else if(zg(si)− zr(si))< 0, it is fill region,3. if(zg(si) − zr(si))(zg(si+1) − zr(si+1))< 0, then the grid cell hasboth cut region and fill region. In this case, we calculate the point oftransition from cut to fill or fill to cut using equation (3.4).Figure 3.4: An example of a pave only cross-section.3.2.6 Cost calculation for the tangent road segmentGiven the cost penalty parameters Cp, Cc, Cf and Cw, respectively, forthe costs of paving, ground cut, ground fill, and waste material, then, eachcost item is calculated asCjpk = CpWLjkCjck = CcVjckCjfk = CfVjfkCjwk = Cw|Vjck − Vjfk|,where Cjpk is the paving cost, Cjck is ground cut cost, Cjfkis the ground fillcost, and Cjwk is cost of waste material.313.3. Surrogate cost model for circular road segmentFor j = 1, 2, · · · ,mk the cost corresponding to the road segment thatconnects the points V Pkj−1 and V Pkj is the sum of all cost items related toearthwork and cost of paving. In particular, the earthwork cost is calculatedasCjEk = Cjck + Cjfk+ Cjwk . (3.14)Thus, the cost calculation of the tangential road segment is divided intotwo parts, the cost of paving and cost of earthwork,TcostEk =mk∑j=1CjEk (3.15a)TcostPk =mk∑j=1Cjpk . (3.15b)3.3 Surrogate cost model for circular roadsegmentSuppose segment ξ ∈ S corresponds to circular arc, that is, a circle ofradius rk that connects TCk and CTk, an example is shown in Figure 3.5.We denote the horizontal projection of this curve by HCk. The circular roadsegment is given in parametric form byrc(s) = (x(s), y(s), z(s)), (3.16)where,x(s) = xck + rk cos θck(s)y(s) = yck + rk sin θck(s)z(s) = zkmk + (z(k+1)1 − zkmk )s,where (xck , yck) is the center of HCk and θck(s) = θTCk + (θCTk − θTCk)s.Let xu−1 and yv−1 be the x boundary and y boundaries at which thehorizontal curve HCk crosses the grid cell Guv, see Figure 3.5. Note thatxu−1 = x0 + (u − 1)Dx and yv−1 = y0 + (v − 1)Dy, where Dx and Dy aredimensions of each grid cell. Next we calculate the set of parameters T kx andT ky where the circular road section crosses the x−boundary and y−boundary,323.3. Surrogate cost model for circular road segmentrespectively.xu−1 = xck + rk cos θk(s)⇒ θk(s) = arccos(xu−1 − xckrk)yv−1 = yck + rk sin θk(s)⇒ θk(s) = arcsin(yv−1 − yckrk).The parameters corresponding to the x boundary cross and the y boundarycross, denoted by sx and sy, respectively, are calculated assx =1(θCTk − θTCk)(arccos(xu−1 − xckrk)− θTCk)(3.17a)sy =1(θCTk − θTCk)(arcsin(yv−1 − yckrk)− θTCk). (3.17b)Therefore, the sets T kx and Tky are defined asT kx = {s∣∣ s = sx} (3.18)T ky = {s∣∣ s = sy}. (3.19)In order to calculate the earthwork cost for each grid cell, we need to de-termine the grid cells through which the circular curve passes. Given theintersection point IPk = (xk, yk), and suppose that the coordinates of TCkand CTk are calculated asTCk = (xtck , ytck), and CTk = (xctk , yctk).First, we determine the points at which the circular curve intersects the gridcells, and then we sort them. The sorting technique is discussed next.The Equation of a circle in rectangular coordinate system is defined as(x− xck)2 + (y − yck)2 = r2k. (3.20)We define xmin, xmax, ymin and ymax asxmin = min{xk, xtck , xctk}, xmax = max{xk, xtck , xctk}ymin = min{yk, ytck , yctk}, ymax = max{yk, ytck , yctk}We calculate the x boundary cross and y boundary cross in the ranges givenbelow.xrange = [xmin, xmax] and yrange = [ymin, ymax].333.3. Surrogate cost model for circular road segmentFigure 3.5: Example of a horizontal curve section3.3.1 Calculation of x-boundary crossesFor each x in the xrange we start by definingx¯ = Dx⌊xDx⌋,where b·c is the floor function which is defined below.Definition 1. Let x ∈ R. Then bxc is defined as:b·c : R→ Z : bxc = max{m ∈ Z : m ≤ x}.That is, bxc is the greatest integer less than or equal to x. Consider theline segment l that connects the points q1 = (x¯, ymin) and q2 = (x¯, ymax),which is given in parametric form asl(t) = q1 + (q2 − q1)t, 0 ≤ t ≤ 1.Essentially, l is a line segment which is parallel to the y−axis, and which lieson the grid lines. We calculate the intersection point(s) of l and the circle343.3. Surrogate cost model for circular road segmentin (3.20) by calculating the value(s) of the parameter t. The value(s) t canbe calculated by solving the equation(x¯− xck)2+(ymin + (ymax − ymin)t− yck)2− r2k = 0, 0 ≤ t ≤ 1, (3.21)There are three possibilities about the solution of t,1) no real solution is found,2) two distinct solutions are found,3) exactly one solution is found.If no real solution is found, then the line segment between q1 and q2 does notintersect the circular curve. However, if two distinct solutions are calculated,say t1 and t2 with 0 ≤ t1 ≤ 1 and 0 ≤ t2 ≤ 1, then we choose the rightparameter by calculating the angle corresponding to each parameter. Forexample, if bp1 =(x¯, ymin + (ymax − ymin)t1)is the point calculated usingthe parameter t1. Then we determine the angle θ1 for bp1, so that ifθTCk ≤ θ1 ≤ θCTk ,x¯ is accepted as x-boundary cross and bp1 is the boundary point, else it isrejected because it is not on the part of the curved road segment. Using thesame procedure, we calculate the y−boundary cross and the correspondingboundary point.Once the x-boundary cross x¯ and the y-boundary cross y¯ are deter-mined, we calculate the corresponding parameters sx and sy using (3.17a)and (3.17b), respectively. Let BPk be the set of distinct boundary pointsat which the circular curve crosses the grid cells, including the start andend points, and Tk be the set of parameters corresponding to BPk. Nowwe sort elements of BPk and Tk. However, since the circle, unlike the linesegment, is not convex set, the sorting is not by increasing or decreasingorder of elements of Tk. The distance measurement from the elements ofBPk to TCk is used for sorting.Let M be the cardinality of Tk, the sets after sorting are defined asTks = {0 = s1 < s2 < s3 < · · · , sM−1 < sM = 1} (3.22a)BPks = {bp1k, bp2k, · · · , bpMk }. (3.22b)353.3. Surrogate cost model for circular road segment3.3.2 Length of circular road sectionThe length of the circular road segment, denoted by Lck is given byLck = rkθck , (3.23)where θck = |θCTk − θTCk |, measured in radians, is the magnitude of centralangle of the curve.3.3.3 Volume of ground cutThe parametric equation of the ground profile for a circular road sectionis given aszg(s) = Auvx(s) +Buvy(s) + Cuv= Auv(xck + rk cos θck(s)) +Buv(yck + rk sin θck(s)) + Cuv=(Auvxck +Buvyck + Cuv)+ rk(Auv cos θck(s) +Buv sin θck(s)).Consider two consecutive parameters si−1, si ∈ Tks . If si−1 and si bracketa cut, then the approximate volume of ground cut for the circular roadsegment over [si−1, si] is obtained by integrating the cross-sectional areaover the effective length of the road in the specified grid cell. The elevationdifference at s ∈ [si−1, si] is calculated byhc(s) = zg(s)− zr(s)=(Auvxck +Buvyck + Cuv)+ rk(Auv cos θck(s) +Buv sin θck(s))−(zkmk + (z(k+1)1 − zkmk )s).Recalling the constraint stated in section (2.3.2), we obtainhc(s) =(Auvxck +Buvyck−zkmk +Cuv)+rk(Auv cos θck(s)+Buv sin θck(s)).The volume of ground cut for circular road segment is computed.Vck =M∑i=1∫ sisi−1ac(s)ds, (3.24)where, ac(s) = Whc(s) +12κh2c(s).363.3. Surrogate cost model for circular road segmentAssuming si−1 and si bracket a cut, the cross-sectional area ac(s) iscalculated as follows. We denote Λ =(Auvxck +Buvyck − zkmk + Cuv).ac(s) = Whc(s) +12κh2c(s)= W[Λ + rk(Auv cos θck(s) +Buv sin θck(s))]+12κ[Λ+rk(Auv cos θck(s) +Buv sin θck(s))]2=[WΛ +12κΛ2]+[WAuvrk + κAuvrkΛ]cos θck(s) +[WBuvrk+κBuvrkΛ]sin θck(s) +12κr2kA2uv cos2 θck(s) +12κr2kB2uv sin2 θck(s)+κAuvBuvr2k sin θck(s) cos θck(s) (3.25)If si−1 and si do not bracket a cut, then ac(s) is defined as 0.We denote ∆θk = θCTk − θTCk . Integrating (3.24) using (3.25), we getV ick =[WΛ +12κΛ2](si − si−1)+1∆θk[WAuvrk + κAuvrkΛ](sin θck(si)− sin θck(si−1))−1∆θk[WBuvrk + κBuvrkΛ](cos θck(si)− cos θck(si−1))+18∆θkκr2kA2uv[2(θck(si)− θck(si−1))+(sin 2θck(si)− sin 2θck(si−1))]+18∆θkκr2kB2uv[2(θck(si)− θck(si−1))−(sin 2θck(si)− sin 2θck(si−1))]+1∆θkκr2kAuvBuv[sin 2θck(si)− sin 2θck(si−1)].The volume of ground cut for a circular section is then computed asVck ={ ∑Mi=1 Vick if si−1, si bracket a cut segment0 otherwise(3.26)3.3.4 Volume of ground fillSimilar to the volume calculation for a ground cut, the ground fill volumeis calculated using the integration of the cross-sectional area of a ground fill.Vfk ={ ∑Mi=1∫ sisi−1af (s)ds if si−1, si bracket a cut segment0 otherwise(3.27)373.4. The overall surrogate cost functionwhere af (s) = −ac(s).3.3.5 Cost calculation for a circular road segmentGiven the cost penalty parameters Cp, Cc, Cf and Cw, respectively, forthe costs of paving, ground cut, fill, and waste material. The cost items forthe circular road section is given asCpk = CpWLckCck = CcVckCfk = CfVfkCwk = Cw‖Vck − Vfk‖.Therefore, the cost corresponding to the circular road section is relatedto the earthwork cost and paving cost. In particular, the earthwork cost iscalculated asCEk = Cck + Cfk + Cwk . (3.28)Thus, the cost calculation for the circular road segment HCk has twoparts, the cost of paving and the cost of earthwork.CcostEk = CEk (3.29a)CcostPk = Cpk . (3.29b)3.4 The overall surrogate cost functionIn this section, we present the surrogate cost model used to calculatethe cost of road construction from the start to the end of the road projectusing the cost formulas in (3.29a), (3.29b), (3.15a), and (3.15b). We notethat, because of the assumption that is made in Chapter 2, the start andend section of the road are required to be tangential road sections. Now,we divide the road alignment, between the start and the end, into segments,see Figure 3.6.383.4. The overall surrogate cost functionFigure 3.6: An example of road segmentsEach segment has a tangent part and a curve part, except for the lastsegment. The orientation of each segment is a tangent-curve. In this way,for k = 1, 2, 3, · · · , N , we have N + 1 segments including the last tangentsegment. Therefore, for i = 1, 2, 3, · · · , N + 1, the overall surrogate costmodel is given asCostE =N∑k=1(TcostEk + CcostEk)+ TcostEN+1 , (3.30)CostP =N∑k=1(TcostPk + CcostPk)+ TcostPN+1 , (3.31)where CostE and CostP are the overall earthwork cost and paving cost,respectively. It is easy to see that the surrogate cost calculation scaleslinearly with the number of intersection points and the number of grid cells.39Chapter 4Multiobjective optimizationmodel for 3D road alignmentIn this chapter, we formulate the three-dimensional road alignment op-timization as a multiobjective optimization problem. The surrogate costmodel developed in Chapter 3 has two components, the cost due to theearthwork and the cost related to the length of the road. There is, usually, aconflict between the two cost components. Thus, the solution of road align-ment optimization, based on cost minimization, should reflect the trade-offbetween the length of the road and the volume of earthwork. During theroad alignment design process, an engineer can have different objectives thatneed to be satisfied. Some of the objectives may favour the shortest roadpossible, while others might favour an indirect path with a lower earthworkcost. In this situation, there is a trade-off between the choice of cost penaltyparameters. The shortest path may incur more cost than al longer path.The multiobjective optimization formulation include the formulation ofthe individual cost components, and the set of constraints. The pavingcosts and earthwork costs are formulated in (3.30) and (3.31). The set ofconstraints are also specified in subsections (2.2.2) and (2.3.1).4.1 Variable definitionThe costs are formulated as functions of decision variables, such as ver-tical intersection points, radius of curvature and the road design elevationsat each vertical points. The parameters of the optimization model are thenumber of intersection points, the cost penalty parameters, and κ whichdepends on the side slope of the road.Let X =(x1, x2, · · · , xN)and Y =(y1, y2, · · · , yN)be the coordinatesof the intersection points, R =(r1, r2, · · · , rN)be the vector of radius ofcurvature, and Z =(z1, z2, · · · , zM)be the vector of elevations of the designroad. Define the objective functions CE and CP ,CE : R3N+M → R : (X,Y,R,Z) 7→ CostE(X,Y,R,Z),404.2. Problem formulationCP : R3N+M → R : (X,Y,R,Z) 7→ CostP (X,Y,R,Z),where CostE and CostP are given in equations (3.30) and (3.31).4.2 Problem formulationThe optimization model accepts the values of parameters and other con-stants as input and returns the values of decision variables. Given the startpoint START = (xs, ys, zs), the end point END = (xe, ye, ze), the maxi-mum allowable gradient Gmax, and the minimum radius of curvature rmin.If z¯ is the vertical off-set from the current ground elevation zg, the simulta-neous optimization of horizontal and vertical alignments is given as follows.Minimize{CostE(X,Y,R,Z) , CostP (X,Y,R,Z)}subject to : Horizontal alignment constraintsFor k = 1, 2, · · · , N,0 ≤ ‖TCk − IPk−1‖ − ‖CTk−1 − IPk−1‖,rmin ≤ rk,Vertical alignment constraintsFor k = 1, 2, · · · , N, and for j = 1, 2, · · ·mk,|zkj − zkj−1 | ≤ Gmax√(xkj − xkj−1)2 + (ykj − ykj−1)2,zg − z¯ ≤ zkj ≤ zg + z¯,Other constraintszkmk = z(k+1)1 , k = 1, 2, · · · , N.4.3 NOMAD: Nonlinear optimization with theMADS algorithmNOMAD [AAC+] is a black box optimization software that implementsthe Mesh Adaptive Direct Search (MADS) algorithm under general non-linear constraints [LD11]. When the objective functions and constraintfunctions defining the optimization problem are typically simulations whichpossess no exploitable properties such as derivatives, and they may fail toevaluate at some trial points, even feasible ones, the term blackbox prob-lem is used to denote this problem class [LD11]. Direct search methods are414.4. Solution proceduremeant for such a context since they use only function evaluations to drivetheir exploration in the space of variables and they can deal with missingfunction values. MADS is an iterative algorithm where at each iteration afinite number of trial points are generated, and the infeasible trial pointsare discarded. The objective function values at the feasible trial points arecompared with the current best feasible objective function value found sofar [AD06].4.4 Solution procedureThe solution procedure for solving the three-dimensional alignment op-timization problem is designed to use the derivative-free optimization algo-rithm. In our numerical tests, we use NOMAD [ALT09, LD11], but anyderivative-free optimization algorithm would work. Although the surrogatecost model in each grid cell is continuously differentiable, the overall costmodel is nondifferentiable because of two reasons. First, if the terrain corre-sponding to adjacent grid cells is approximated by different linear functions(or planes), the volume of earthwork is not necessarily equal. Thus, thecost function in each grid cell are not necessarily equal, and this results ina sharp turning on the boundary of the grid cells. Second, at a point wherethe ground changes from cut to fill or fill to cut, the cost function might havea sharp turning due to the difference in cost penalty parameters for groundcut and ground fill. Moreover, the surrogate cost model is nonconvex, thiscan be seen by noting the solution alignments of two initial roads, depictedin Figure 5.8. Thus, a derivative-free optimization algorithm is an appealingmethod to solve the problem.Cost =(θCE + (1− θ)CP), (4.1)where 0 ≤ θ ≤ 1.Given a feasible initial alignment, we run the solver for different values θand solved the problem. Thus, the optimization algorithm finds the optimaldecision vectors for both vertical and horizontal alignments simultaneously.42Chapter 5Numerical resultsThis chapter presents a numerical experiment that was employed to in-vestigate the proposed model. An actual ground profile in California withan area of 500 meters by 1000 meters is considered. The contours and three-dimensional view of the map are displayed in Figure 5.1 and Figure 5.2. (Wenote that the contour map scale is in the ratio of 1 m along the x-axis isequal to 2 m along the y−axis, and 1 m along the z−axis.) Based on theground profile, two example scenarios are generated and tested. In the firstscenario, a single alignment is generated to examine the effect of the costpenalty parameters on the solution found. In the second example, two dif-ferent alignments are generated to investigate whether the problem is convexor not. This is tested by looking at the optimal solutions of both alignments.In the second example, each alignment is tested on two different parameters.Based on discussions with our industrial partner (Softree Technical systemsInc.), the model parameters and design constants are fixed as presented inTable 5.1.Table 5.1: Design constants and parameter values.Parameter valueCp 0.3Cc 1Cf 0.5Cw 2Gmax (%) 15rmin (m) 20mk 5N 4Road width (m) 5κ 1Note that, for experimental purposes we set mk = 5 for all k, and wevaried the value of parameter θ that we introduced in Equation (4.1).43Chapter 5. Numerical resultsFigure 5.1: Three-dimensional view of the terrain. Note that the aspectratio for the x-axis, the y-axis, and the z-axis is, respectively, 10:10:1.104.767104.767106.6118106.6118106.6118106.6118106.6118106.6118106.6118106.6118106.6118106.6118106.6118106.6118106.6118108.4566 108.4566108.4566108.4566108.4566108.4566108.4566110.3014110.3014110.3014110.3014110.3014110.3014110.3014110.3014110.3014112.1462112.1462112.1462112.1462112.1462112.1462112.1462112.1462112.1462112.1462112.1462112.1462112.1462113.991113.991113.991113.991113.991113.991113.99113.991113.991113.991113.991113.991113.991115.8358115.8358115.8358115.8358115.8358115.8358115.8358115.8358115.8358115.8358115.8358117.6807117.6807117.6807117.6807117.6807117.6807117.6807117.6807119.5255119.5255119.5255119.5255119.5255119.5255119.5255121.3703121.3703121.3703121.3703121.3703123.2151123.2151123.2151123.2151125.0599125.0599125.0599126.9047xy5 10 15 20 25 30 35 40 45102030405060708090Figure 5.2: The contour map. The x-axis and y−axis are in tens of meters.445.1. Experimental setup5.1 Experimental setupThree artificial roads that satisfy the basic design and geometric con-straints were generated. All experiments were tested with a MATLABR2013b code performed on a Dell workstation equipped with an Intel(R)Core(TM) i7 2.8GHz processor, and a 16 GB of RAM using 64-bit Windowsoperating system. The two examples differ in that, in the first experiment, asingle alignment is solved for several parameters values, while in the secondexperiment, two different alignments are solved for the same parameter. Infact, we tested the two alignments on two different parameters, θ = 0 andθ = 0.5. In this second experiment, we used the same set of constraints andstopping criteria for both alignments. Each case study demonstrates theability of the model to solve a three-dimensional road alignment optimiza-tion problem in a reasonable amount of time.The MATLAB version of NOMAD [AAC+] is used to solve the prob-lems. The parameters of the solver that are set during the solution processinclude: Maximum time (maxtime), maximum number of function evalua-tion (maxfeval), and minimum mesh size.5.2 Case study 1An artificial three-dimensional road alignment that satisfies the con-straint set is generated. The road profile is calculated from the terraindata. The horizontal alignment is shown in Figure 5.3.Based on the surrogate cost model developed in Chapter 3, and themulti-objective optimization problem formulation in Chapter 4, we solvedthe problem for five different values of θ: θ ∈ {0, 0.25, 0.5, 0.75, 1}. We usedfour intersection points and maximum gradient of 12%. Because we fixedthe number of vertical grades between two consecutive intersection points,the number of variables depends on the number of intersection points only,and this makes the number of variables equal to forty.455.2. Case study 10 5 10 15 20 25 30 35 40 45 500102030405060708090100104.767104.767106.6118106.6118106.6118106.6118106.6118106.6118106.6118106.6118106.6118106.6118106.6118106.6118106.6118108.4566108.4566108.4566108.4566108.4566108.4566108.4566110.3014110.3014110.3014110.3014110.3014110.3014110.3014110.3014110.3014112.1462112.1462112.1462112.1462112.1462112.1462112.1462112.1462112.1462112.1462112.1462112.1462112.1462113.991113.991113.991113.991113.991113.991113.99113.991113.991113.991113.991113.991113.991115.8358115.8358115.8358115.8358115.8358115.8358115.8358115.8358115.8358115.8358115.8358117.6807117.6807117.6807117.6807117.6807117.6807117.6807117.6807119.5255119.5255119.5255119.5255119.5255119.5255119.5255121.3703121.3703121.3703121.3703121.3703123.2151123.2151123.2151123.2151125.0599125.0599125.0599126.9047xyFigure 5.3: Initial horizontal alignment in case study 1. Note that the x-axisand y−axis are in tens of meters.We used two stopping criteria, the number of black box evaluation andthe mesh size, which are set, for this case study, to 10,000 and 0.5, respec-tively.5.2.1 Experimental results 1The result of the experiment for case study 1 indicates that, when thevalue of θ is close to zero, the optimal road alignment, as expected, choosesthe shortest route. That is, the length of the road is minimized when θ valuesis zero. On the other hand, when θ values are close to one, the optimal routeminimizes the earthwork volume, which may result in a longer route. Thisresult indicates that, there are two conflicting objectives in the problem,the length of the road, and the amount of earthwork. Therefore, the roadoptimization problem is indeed multi-objective optimization problem. Theoptimal horizontal alignment corresponding to each of the parameters is465.2. Case study 1depicted in Figure 5.4.0 5 10 15 20 25 30 35 40 45 500102030405060708090100xy  104.767104.767106.6118106.6118106.6118106.6118106.6118106.6118106.6118106.6118106.6118106.6118106.6118106.6118106.6118108.4566108.4566108.4566108.4566108.4566108.4566108.4566110.3014110.3014110.3014110.3014110.3014110.3014110.3014110.3014110.3014112.1462112.1462112.1462112.1462112.1462112.1462112.1462112.1462112.1462112.1462112.1462112.1462112.1462113.991113.991113.991113.991113.991113.991113.99113.991113.991113.991113.991113.991113.991115.8358115.8358115.8358115.8358115.8358115.8358115.8358115.8358115.8358115.8358115.8358117.6807117.6807117.6807117.6807117.6807117.6807117.6807117.6807119.5255119.5255119.5255119.5255119.5255119.5255119.5255121.3703121.3703121.3703121.3703121.3703123.2151123.2151123.2151123.2151125.0599125.0599125.0599126.9047Int i al  al i gnme ntOpt i Al i gn1( t he t a=0;)Opt i Al i gn2( t he t a=0.25;)Opt i Al i gn3(t he t a=0.5;)Opt i Al i gn4(t he t a=0.75;)Opt i Al i gn5(t he t a=1;)Figure 5.4: The optimal horizontal alignments for case study 1. The x-axisand y−axis are in tens of meters.As it can be seen from Figure 5.4, the optimal horizontal alignment cor-responding to θ = 0 is not straight because of the box constraints at eachintersection point. The optimal vertical alignments for alignments corre-sponding to θ = 0 and θ = 1 are depicted in Figure 5.5 and Figure 5.6,respectively. The search stopped because of the maximum number of blackbox evaluations. Table 5.2 presents the numerical result of the experiment.475.3. Case study 20 200 400 600 800 1000 100 110 120 130 140 150 160 Horizontal distance (m)  Elevation (m)     Optimal vertical alignment  Ground profile  Figure 5.5: The optimal vertical alignments for θ = 1.0 200 400 600 800 1000100110120130140150160Ho r i zo nt al  di s t anc e  ( m)Elevation (m)  Opt i mal  ve r t i c al  al i g nme ntG r o und pr o fi l eFigure 5.6: The optimal vertical alignments for θ = 0.5.3 Case study 2In this experiment, we set the number of function evaluation to 70,000and the minimum mesh size to 0.05, in order to see if the two alignmentsconverge to the same or near the same solution. The inputs and modelparameters that we used in case study 1 are also used for this example. For485.3. Case study 2Table 5.2: Numerical results for case study 1.θ objective black-box mesh solution Cp CEvalue evaluations size time (s) ×103 ×1040 1105.1 10,000 0.5 358.62 1.19 334.80.25 1046.26 10,000 0.5 414.38 4.03 7.390.5 2165.04 10,000 0.5 429.43 4.24 6.530.75 2348.4 10,000 0.5 461.91 4.30 6.561 2885.06 10,000 0.5 436.36 4.40 1.42instance, we used four intersection points and maximum gradient of 12%.The horizontal alignment is shown in Figure 5.7. As mentioned earlier, thiscase study examines whether the problem is convex or not.0 10 20 30 40 500102030405060708090100xy  104.767104.767106.6118106.6118106.6118106.6118106.6118106.6118106.6118106.6118106.6118106.6118106.6118106.6118106.6118108.4566108.4566108.4566108.4566108.4566108.4566108.4566110.3014110.3014110.3014110.3014110.3014110.3014110.3014110.3014110.3014112.1462112.1462112.1462112.1462112.1462112.1462112.1462112.1462112.1462112.1462112.1462112.1462112.1462113.991113.991113.991113.991113.991113.991113.99113.991113.991113.991113.991113.991113.991115.8358115.8358115.8358115.8358115.835815.8358115.8358115.8358115.8358115.8358115.8358117.6807117.6807117.6807117.6807117.6807117.6807117.6807117.6807119.5255119.5255119.5255119.5255119.5255119.5255119.5255121.3703121.3703121.3703121.3703121.3703123.2151123.2151123.2151123.2151125.0599125.0599125.0599126.9047Int i al  al i g nme ntInt i al  al i g nme nt 2Figure 5.7: The horizontal alignments for case study 2. The x-axis andy−axis are in tens of meters.495.3. Case study 25.3.1 Experimental results 2The outcome of this experiment shows that, when θ = 0.5, both align-ments converge to different solutions. The optimal horizontal alignments areshown in Figure 5.8 . Even for the case when θ = 0, the vertical alignmentsare not the same, the optimal alignments are shown in Figure 5.9 and Figure5.10. Then it suffices to conclude that the problem is nonconvex. Table 5.3and Table 5.4 present the numerical results for both alignments.0 10 20 30 40 500102030405060708090100xy  104.767104.767106.6118106.6118106.6118106.6118106.6118106.6118106.6118106.6118106.6118106.6118106.6118106.6118106.6118108.4566108.4566108.4566108.4566108.4566108.4566108.4566110.314110.3014110.3014110.3014110.3014110.3014110.3014110.3014110.3014112.1462112.1462112.1462112.1462112.1462112.1462112.1462112.1462112.1462112.1462112.1462112.1462112.1462113.991113.991113.991113.991113.991113.991113.99113.991113.991113.991113.991113.991113.991115.8358115.8358115.8358115.8358115.8358115.8358115.8358115.8358115.8358115.8358115.8358117.6807117.6807117.6807117.6807117.6807117.6807117.6807117.6807119.5255119.5255119.5255119.5255119.5255119.5255119.5255121.3703121.3703121.3703121.3703121.3703123.2151123.2151123.2151123.2151125.0599125.0599125.0599126.9047Opt i Al i g n1 ( t he t a= 0 .5 ; )Opt i Al i g n2 ( t he t a= 0 .5 ; )Figure 5.8: Two optimal horizontal alignments computed from differentstarting alignments. Hence the problem is nonconvex. The x-axis andy−axis are in tens of meters.Table 5.3: Numerical results for case study 2 when θ = 0.5.Optimal obj. BBE min. mesh solution (s) Stoppedalignment value size time due to1 95.18 26,364 0.05 3,317 min. mesh size2 156.14 19,686 0.05 3,131 min. mesh size505.3. Case study 20 10 20 30 40 500102030405060708090100xy  104.767104.767106.6118106.6118106.6118106.6118106.6118106.6118106.6118106.6118106.6118106.6118106.6118106.6118106.6118108.4566 108.4566108.4566108.4566108.4566108.4566108.4566110.3014110.3014110.3014110.3014110.3014110.3014110.3014110.3014110.3014112.1462112.1462112.1462112.1462112.1462112.1462112.1462112.1462112.1462112.1462112.1462112.1462112.1462113.991113.991113.991113.991113.991113.991113.99113.991113.991113.991113.991113.991113.991115.8358115.8358115.8358115.8358115.8358115.8358115.8358115.8358115.8358115.8358115.8358117.6807117.6807117.6807117.6807117.6807117.6807117.6807117.6807119.5255119.5255119.5255119.5255119.5255119.5255119.5255121.3703121.3703121.3703121.3703121.3703123.2151123.2151123.2151123.2151125.0599125.0599125.0599126.9047Opt i Al i g n1 ( t he t a= 0 ; )Opt i Al i g n2 ( t he t a= 0 ; )Figure 5.9: The optimal horizontal alignments. Note that the x-axis andy−axis are in tens of meters.0 100 200 300 400 500 600 700 800 100 110 120 130 140 150 160 Horizontal distance (m)  Elevation (m)     Op timal vertical alignment1 (Th eta= 0)  Optimal vertical alignment2 (Theta= 0)  Figure 5.10: The optimal horizontal alignments.515.4. Implication of the resultsTable 5.4: Numerical results for case study 2 when θ = 0.Optimal obj. BBE min. mesh solution (s) Stoppedalignment value size time due to1 1222.61 63,771 0.05 4,692 min. mesh size2 1267.95 76,931 0.05 4,833 min. mesh size5.4 Implication of the resultsThe two case studies results show that, the problem, as multiobjectiveoptimization problems is dependent on the values of the cost penalty param-eters. The results, further, imply that role of the length dependent penaltyparameter is in conflict with other penalty parameters. For high valuesof paving cost penalty parameter and smaller values of other parameters,the problem tends to minimizes the cost (or objective) value by minimizingthe length of the road. On the other hand, for small values of paving costpenalty parameters, and higher values of the other parameters, the prob-lems minimizes the objective value by minimizing the volume of earthworkbetween the end points.Moreover, the numerical results of the experiment indicates that, theproblem is nonconvex. This is shown by the result obtained from Experi-ment 5.8. In this experiment, the number of black box evaluations and theminimum mesh size are set to let the program run until the two differentalignments converge. The results indicate that, even for the extreme case,when θ = 0, the solution alignments are different.Therefore, the optimal alignment is highly dependent on the choice ofthe cost penalty parameters.52Chapter 6Conclusion and Future workIn this thesis, a surrogate cost model is developed to solve a three-dimensional road alignment optimization problem. Optimal vertical andhorizontal alignments are determined simultaneously. The NOMAD solveris used to solve the problems, but any derivative-free optimization solvercould be applied. Results from the two case studies are presented whichshow that, the problem is multiobjective in terms of cost minimization. Asa result, the cost items can be classified into, those which depend on thelength of the road and those depending on the volume of the earthwork.They are, usually, conflicting because, during the optimization process, thelength dependent costs prefer the shortest route between the two end points,while the cost items that depend on the amount of earthwork choose a routewith the minimum amount of earthwork. Therefore, optimal solutions arevery dependent on the values of cost penalty parameters. Furthermore, thetime of solution for each experiment test shows that, the surrogate modelis solved in a reasonable amount of time, which is expected as it is only anapproximation of the numerically expensive true model.6.1 Future workWe have identified four cost penalty parameters, namely, the ground cutcost penalty parameter Cc, ground fill cost penalty parameter Cp, wastematerial cost penalty parameter Cw, and the paving cost penalty parameterCp. The surrogate can only approximate the cost of earthwork, waste mate-rial, and pavement. Any additional costs items such as transportation costcan be incorporated into the model. Even though the overall surrogate costis not diffenrentiable, it is almost differentiable, that is, the set of pointsat which the cost function is non-differentiable is finite. So, for future thegradient of the cost model can be calculated to make use of gradient basedalgorithms. In this thesis, the values of these parameters are obtained fromindustry partners and are fixed. However, we are also aware that the op-timal choice of this parameters should be done systematically, since theirvalues are highly affecting the optimal alignments.53Bibliography[AAC+] M.A. Abramson, C. Audet, G. Couture, J.E. Dennis, Jr.,S. Le Digabel, and C. Tribes. The NOMAD project. Softwareavailable at http://www.gerad.ca/nomad. → pages 41, 45[AD06] C. Audet and J. Dennis. Mesh adaptive direct search algorithmsfor constrained optimization. SIAM Journal on Optimization,17(1):188–217, 2006. → pages 42[AH11] B.A. Al-Hadad. An approach to the highway alignment develop-ment process using genetic algorithm based optimisation. PhDthesis, University of Nottingham, 2011. → pages 2, 5[ALT09] C. Audet, S. Le Digabel, and C. Tribes. NOMAD user guide.Technical Report G-2009-37, Les cahiers du GERAD, 2009. →pages 42[ASAC05] K. Aruga, J. Sessions, A. Akay, and W. Chung. Simultaneousoptimization of horizontal and vertical alignments of forest roadsusing tabu search. International Journal of Forest Engineering,16(2):137–151, 2005. → pages 5, 8[ASZ08] C. Audet, G. Savard, and W. Zghal. Multiobjective optimizationthrough a series of single-objective formulations. SIAM Journalon Optimization, 19(1):188–210, 2008. → pages 9[CGF89] E.P. Chew, C.J. Goh, and T.F. Fwa. Simultaneous optimizationof horizontal and vertical alignments for highways. Transporta-tion Research Part B: Methodological, 23(5):315 – 329, 1989. →pages 1, 5, 6, 7, 11[CL06] J. Cheng and Y. Lee. . Journal of Transportation Engineering,132(12):913–920, 2006. → pages 5[CMVV11] A. Custo´dio, J.F. Madeira, A.F. Vaz, and L.N. Vicente. Directmultisearch for multiobjective optimization. SIAM Journal onOptimization, 21(3):1109–1140, 2011. → pages 954Bibliography[DS06] M.J. De Smith. Determination of gradient and curvature con-strained optimal paths. Computer-Aided Civil and Infrastruc-ture Engineering, 21(1):24–38, 2006. → pages 5[Eas88] S.M. Easa. Selection of roadway grades that minimize earthworkcost using linear programming. Transportation Research Part A:General, 22(2):121–136, 1988. → pages 4[FCS02] T. Fwa, W. Chan, and Y. Sim. . Journal of TransportationEngineering, 128(5):395–402, 2002. → pages 4, 5[GAA09] A.B. Go¨ktepe, S. Altun, and P. Ahmedzade. Optimization ofvertical alignment of highways utilizing discrete dynamic pro-gramming and weighted ground line. Turkish Journal of Engi-neering and Environ-mental Sciences, 33(2):105–116, 2009. →pages 4, 5[GCF88] C.J. Goh, E.P. Chew, and T.F. Fwa. Discrete and continuousmodels for computation of optimal vertical highway alignment.Transportation Research Part B: Methodological, 22(6):399 –409, 1988. → pages 4, 5[Hay70] R.W. Hayman. Optimization of vertical alignment for highwaysthrough mathematical programming. Highway Research Record,1970. → pages 4[HKL11] W.L. Hare, V.R. Koch, and Y. Lucet. Models and algorithms toimprove earthwork operations in road design using mixed integerlinear programming. European Journal of Operational Research,215(2):470 – 480, 2011. → pages 2, 5[JJS00] J. Jong, M.K. Jha, and P. Schonfeld. Preliminary highway de-sign with genetic algorithms and geographic information sys-tems. Computer-Aided Civil and Infrastructure Engineering,15(4):261–271, 2000. → pages 2, 3[JM07] M.K. Jha and A. Maji. A multi-objective genetic algorithm foroptimizing highway alignments. In Computational Intelligencein Multicriteria Decision Making, IEEE Symposium on, pages261–266, April 2007. → pages 1[JS03] J. Jong and P. Schonfeld. An evolutionary model for simultane-ously optimizing three-dimensional highway alignments. Trans-55Bibliographyportation Research Part B: Methodological, 37(2):107 – 128,2003. → pages 5[JS04] M.K. Jha and P. Schonfeld. A highway alignment optimizationmodel using geographic information systems. TransportationResearch Part A: Policy and Practice, 38(6):455 – 481, 2004. →pages 1[KJLS04] E. Kim, M.K. Jha, D.J. Lovell, and P. Schonfeld. Intersectionmodeling for highway alignment optimization. Computer-AidedCivil and Infrastructure Engineering, 19(2):119–129, 2004. →pages 1[KJS12] M.W. Kang, M.K. Jha, and P. Schonfeld. Applicability of high-way alignment optimization models. Transportation ResearchPart C: Emerging Technologies, 21(1):257 – 286, 2012. → pages5[KJSK07] E. Kim, M. Jha, P. Schonfeld, and H. Kim. . Journal of Trans-portation Engineering, 133(2):71–81, 2007. → pages 5[KL13] S. Koziel and L. Leifsson. Surrogate-Based Modeling and Opti-mization. Springer, 2013. → pages 6[KY11] S. Koziel and X. Yang. Computational optimization, methodsand algorithms, volume 356. Springer, 2011. → pages 6, 7[LD11] S. Le Digabel. Algorithm 909: Nomad: Nonlinear optimizationwith the mads algorithm. ACM Transactions on MathematicalSoftware (TOMS), 37(4):44:1–44:15, February 2011. → pages41, 42[LPZL13] W. Li, H. Pu, H. Zhao, and W. Liu. Approach for optimizing 3dhighway alignments based on two-stage dynamic programming.Journal of Software, 8(11):2967–2973, 2013. → pages 1, 5, 6[Mie99] K. Miettinen. Nonlinear multiobjective optimization, volume 12.Springer, 1999. → pages 9[MS81] R.H. Mayer and R.M. Stark. Earthmoving logistics. Journal ofthe Construction Division, 107(2):297–312, 1981. → pages 4[Nic73] A.J. Nicholson. A variational approach to optimal route location.PhD thesis, University of Catenrbury, 1973. → pages 2, 5, 656Bibliography[OEC73] OECD. Optimization of road alignment by the use of computers.ProQuest, 1973. → pages 2[Par77] N.A. Parker. Rural highway route corridor selection. Trans-portation Planning and Technology, 3(4):247–256, 1977. →pages 5[QHS+05] N.V. Queipo, R.T. Haftka, W. Shyy, T. Goel, R. Vaidyanathan,and P.K. Tucker. Surrogate-based analysis and optimization.Progress in Aerospace Sciences, 41(1):1 – 28, 2005. → pages 6[Rev] Transportation In Canada 2011 Comprehensive Review. Minis-ter of public works and government services, canada, 2012. Cat.No. T1-23A/2011E-PDF, ISSN 1482-1311. → pages 1[SJK06] P.M. Schonfeld, J. Jong, and E. Kim. Intelligent road design,volume 19. WIT Press, 2006. → pages 1, 2, 3, 4, 5, 7, 8, 9[Tri87] D. Trietsch. Comprehensive design of highway networks. Trans-portation Science, 21(1):26–35, 1987. → pages 2[WPR05] X. Wan, J.F. Pekny, and G.V. Reklaitis. Simulation-based opti-mization with surrogate modelsapplication to supply chain man-agement. Computers & Chemical Engineering, 29(6):1317 –1328, 2005. → pages 6[Yu74] P.L. Yu. Cone convexity, cone extreme points, and nondomi-nated solutions in decision problems with multiobjectives. Jour-nal of Optimization Theory and Applications, 14(3):319–377,1974. → pages 957


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



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"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items