Science, Irving K. Barber Faculty of (Okanagan)
Computer Science, Mathematics, Physics and Statistics, Department of (Okanagan)
DSpace
UBCO
Akhmet, Ayazhan
Master of Science - MSc
University of British Columbia
This thesis defines a multi-objective optimization model that seeks to find road profiles that would be optimal for the manufacturers in terms of the road construction cost and at the same time for the users in terms of the vehicle operating costs, specifically in terms of the fuel consumption cost. The research implements and validates the formula for the fuel consumption cost. It further presents and examines a variety of well-known methods: three classical scalarization techniques (the ε-constraint method, weighted sum method, and weighted metric methods) and two widely-used evolutionary methods (NSGA-II and FP-NSGA-II). Moreover, to accelerate the performance of the chosen scalarization approaches, a warm start strategy is proposed. Numerical experiments are performed on 30 road samples for Caterpillar 793D off-highway trucks to determine the most robust approach for the proposed multi-objective optimization problem. The results are analyzed using the commonly-used performance indicators, namely, hypervolume (to assess the convergence of solutions), spacing (to assess the diversity of solutions), and CPU time (to assess the speed). The research determines that the most promising and recommended method for the proposed problem is the ε-constraint method (successfully solved approximately 75% of test problems) followed by the weighted sum method (successfully solved approximately 50% of test problems). Moreover, the research shows that the warm start strategy improves the performance of the scalarization techniques.
Multi-Objective Optimization forRoad Vertical Alignment DesignbyAyazhan AkhmetB.Sc., Nazarbayev University, 2019A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF SCIENCEinTHE COLLEGE OF GRADUATE STUDIES(Computer Science)THE UNIVERSITY OF BRITISH COLUMBIA(Okanagan)June 2021© Ayazhan Akhmet, 2021The following individuals certify that they have read, and recommendto the College of Graduate Studies for acceptance, a thesis/dissertation en-titled:Multi-Objective Optimization for Road Vertical AlignmentDesignsubmitted by Ayazhan Akhmet in partial fulfilment of the requirementsof the degree of Master of ScienceDr. Warren Hare, Irving K. Barber School of Arts and ScienceSupervisorDr. Yves Lucet, Irving K. Barber School of Arts and ScienceSupervisory Committee MemberDr. Solomon Tesfamariam, School of EngineeringSupervisory Committee MemberDr. Amir Ardestani-Jaafari, Faculty of ManagementUniversity ExamineriiAbstractThis thesis defines a multi-objective optimization model that seeks tofind road profiles that would be optimal for the manufacturers in terms ofthe road construction cost and at the same time for the users in terms ofthe vehicle operating costs, specifically in terms of the fuel consumptioncost. The research implements and validates the formula for the fuel con-sumption cost. It further presents and examines a variety of well-knownmethods: three classical scalarization techniques (the ε-constraint method,weighted sum method, and weighted metric methods) and two widely-usedevolutionary methods (NSGA-II and FP-NSGA-II). Moreover, to acceler-ate the performance of the chosen scalarization approaches, a warm startstrategy is proposed.Numerical experiments are performed on 30 road samples for Caterpillar793D off-highway trucks to determine the most robust approach for the pro-posed multi-objective optimization problem. The results are analyzed usingthe commonly-used performance indicators, namely, hypervolume (to assessthe convergence of solutions), spacing (to assess the diversity of solutions),and CPU time (to assess the speed).The research determines that the most promising and recommendedmethod for the proposed problem is the ε-constraint method (successfullysolved approximately 75% of test problems) followed by the weighted summethod (successfully solved approximately 50% of test problems). Moreover,the research shows that the warm start strategy improves the performanceof the scalarization techniques.iiiLay SummaryA high-quality road infrastructure is important in ensuring the economicgrowth of the country. But the expenses for designing new roads are notcheap. Along with large capital investments for road construction, roads alsoincur future user costs. In this thesis, we present an optimization framework,that takes into account the cost of building the road and the cost of usingthe road, as well as strategies for solving it. Our experiments determine thatthe most recommended method for our problem is the ε-constraint methodwith a warm start.ivTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . xDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiChapter 1: Introduction . . . . . . . . . . . . . . . . . . . . . . . 11.1 Motivation and objective . . . . . . . . . . . . . . . . . . . . . 11.2 Research approach . . . . . . . . . . . . . . . . . . . . . . . . 21.3 Thesis outline . . . . . . . . . . . . . . . . . . . . . . . . . . . 3Chapter 2: Preliminaries . . . . . . . . . . . . . . . . . . . . . . 42.1 Vertical alignment design . . . . . . . . . . . . . . . . . . . . 42.2 Multi-objective optimization . . . . . . . . . . . . . . . . . . . 52.2.1 Problem formulation . . . . . . . . . . . . . . . . . . . 62.2.2 Pareto optimality . . . . . . . . . . . . . . . . . . . . . 62.2.3 Pareto front filtering . . . . . . . . . . . . . . . . . . . 8Chapter 3: Model . . . . . . . . . . . . . . . . . . . . . . . . . . 113.1 Road construction cost . . . . . . . . . . . . . . . . . . . . . . 113.1.1 Variables and parameters . . . . . . . . . . . . . . . . 113.1.2 Formula . . . . . . . . . . . . . . . . . . . . . . . . . . 123.2 Vehicle operating cost . . . . . . . . . . . . . . . . . . . . . . 12vTABLE OF CONTENTS3.2.1 Literature review . . . . . . . . . . . . . . . . . . . . . 123.2.2 Variables and parameters . . . . . . . . . . . . . . . . 143.2.3 Formula . . . . . . . . . . . . . . . . . . . . . . . . . . 15Chapter 4: Methods . . . . . . . . . . . . . . . . . . . . . . . . . 214.1 Literature review . . . . . . . . . . . . . . . . . . . . . . . . . 214.2 ε-constraint method . . . . . . . . . . . . . . . . . . . . . . . 244.2.1 Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . 244.2.2 Dynamic ε selection . . . . . . . . . . . . . . . . . . . 254.3 Weighted sum method . . . . . . . . . . . . . . . . . . . . . . 274.3.1 Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . 284.3.2 Dynamic weight selection . . . . . . . . . . . . . . . . 284.4 Weighted metric method . . . . . . . . . . . . . . . . . . . . . 304.4.1 Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . 314.4.2 Dynamic weight selection . . . . . . . . . . . . . . . . 324.5 Genetic algorithms . . . . . . . . . . . . . . . . . . . . . . . . 324.5.1 NSGA-II . . . . . . . . . . . . . . . . . . . . . . . . . 344.5.2 FP-NSGA-II . . . . . . . . . . . . . . . . . . . . . . . 354.6 Warm start strategy . . . . . . . . . . . . . . . . . . . . . . . 38Chapter 5: Numerical results . . . . . . . . . . . . . . . . . . . 405.1 Performance indicators . . . . . . . . . . . . . . . . . . . . . . 405.1.1 Literature review . . . . . . . . . . . . . . . . . . . . . 405.1.2 Hypervolume . . . . . . . . . . . . . . . . . . . . . . . 425.1.3 Spacing . . . . . . . . . . . . . . . . . . . . . . . . . . 445.1.4 CPU time . . . . . . . . . . . . . . . . . . . . . . . . . 455.2 Experiment description . . . . . . . . . . . . . . . . . . . . . 465.3 Experimental results . . . . . . . . . . . . . . . . . . . . . . . 485.3.1 Effect of warm start strategy on speed . . . . . . . . . 485.3.2 Examination of performance indicators . . . . . . . . . 495.3.3 Comparison of time . . . . . . . . . . . . . . . . . . . 515.3.4 Summary of results . . . . . . . . . . . . . . . . . . . . 55Chapter 6: Conclusion . . . . . . . . . . . . . . . . . . . . . . . . 566.1 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 566.2 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . 576.2.1 Model extension . . . . . . . . . . . . . . . . . . . . . 576.2.2 Methods extension . . . . . . . . . . . . . . . . . . . . 57Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58viTABLE OF CONTENTSAppendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65Appendix A: Tables . . . . . . . . . . . . . . . . . . . . . . . . . . 65Appendix B: Figures . . . . . . . . . . . . . . . . . . . . . . . . . . 66viiList of TablesTable 5.1 Timeout limits . . . . . . . . . . . . . . . . . . . . . . 46Table 5.2 Parameters for NSGA-II and FP-NSGA-II (Part 1) . . 47Table 5.3 Parameters for NSGA-II and FP-NSGA-II (Part 2) . . 48Table 5.4 Average speed improvement due to warm start . . . . 49Table A.1 Parameters α, β1 and β2 . . . . . . . . . . . . . . . . . 65Table A.2 Idle fuel rates for various vehicle types . . . . . . . . . 65viiiList of FiguresFigure 1.1 Annual vehicle operating cost due to poor roads . . . 2Figure 2.1 Graphical illustration of vertical alignment . . . . . . 5Figure 2.2 Pareto dominance example . . . . . . . . . . . . . . . 6Figure 2.3 Graphical illustration of ideal and nadir points . . . . 8Figure 2.4 Graphical illustration of non-dominated sorting . . . 9Figure 3.1 Graphical illustration of tractive and resistance forces 18Figure 3.2 Graphical illustration of traffic flows . . . . . . . . . . 19Figure 4.1 Graphical illustration of ε-constraint method . . . . . 27Figure 4.2 Graphical illustration of weighted sum method . . . . 28Figure 4.3 Graphical illustration of weighted metric method . . . 31Figure 4.4 Graphical illustration of crowding distance . . . . . . 35Figure 4.5 Graphical illustration of NSGA-II and FP-NSGA-II . 37Figure 4.6 Graphical illustration of warm start strategy . . . . . 38Figure 5.1 Graphical illustration of hypervolume . . . . . . . . . 43Figure 5.2 Pareto front examples . . . . . . . . . . . . . . . . . . 44Figure 5.3 Graphical illustration of spacing . . . . . . . . . . . . 45Figure 5.4 Results for hypervolume indicator . . . . . . . . . . . 50Figure 5.5 Results for spacing indicator . . . . . . . . . . . . . . 51Figure 5.6 Comparison of time . . . . . . . . . . . . . . . . . . . 53Figure 5.7 Extra comparisons of time . . . . . . . . . . . . . . . 54Figure B.1 Resulting Pareto fronts for 30 test problems . . . . . 66ixAcknowledgementsI take this opportunity to sincerely thank everyone who supported andguided me in this journey.First and foremost, I express a deep sense of gratitude and appreciationto my supervisor, Dr. Warren Hare, for leading me throughout the courseof this research. This work has not been realized without his constant en-couragement, invaluable advice, and cordial support.Acknowledgment is extended to my graduate committee member,Dr. Yves Lucet, for giving detailed feedback and wholesome insights. Thecourse of numerical optimization that was taught by him broadened myknowledge and established a solid foundation for my research.I am very grateful to my graduate committee member, Dr. SolomonTesfamariam, for giving constructive reviews that put me on a right trackat the start of my thesis.Special thanks to our industrial partner, Softree Technical Systems Inc.,for supplying us with practical datasets, as well as with insightful feedback.I gratefully acknowledge the funding provided by The University ofBritish Columbia (UBC) and the Natural Sciences and Engineering ResearchCouncil (NSERC).Deserving a special mention, I express my sincere appreciation to theComputer-Aided Convex Analysis (CA2) laboratory for providing persistentaccess to computers and services.I owe a great debt of gratitude to my family and my friends for their love,faith, and understanding. Especially, I would like to thank my labmates,Nusrat Suzana Momo and Khandoker Md Ayman, for sharing great ideasand having happy times together as graduate students.Above all, I thank God for giving all grace, good health, and strengthto complete this work.xDedicationDedicated to the individuals whom you can tell your night dreams.xiChapter 1Introduction1.1 Motivation and objectiveThe era of road design began about 2000 years ago when the Romansinitiated the development of techniques for building durable roads to easethe mobility of armies within their empire [Abr13]. Since then, the roaddesign process has been influenced by the idea that better roads result inlower road user costs. Well-designed roads increased the probability that theRoman armies would reach their destinations with lower fuel (analogously,food and water supplies), lower equipment repair costs, fewer accidents, aswell as in less time. Today, the goal of road design is still analogous.Even though we attempt to serve several objectives while constructing anew road, we often consider only the minimization of the road constructioncost and neglect the basic importance of user costs affecting the country’seconomy. This is similar to how people often choose what to eat and drinkwithout considering the simple nutritional benefits of food options. There-fore, we should increase consideration of the road user costs while planninga road.According to the most recent report for the Canadian Automobile Asso-ciation (CAA), the poor roads for Canada results in higher vehicle operatingcosts of about $3 billion per year [CPC21]. The statistics by the provincesare given in Figure 1.1. These expenses encompass several cost componentssuch as fuel consumption, tire usage, depreciation, repair, and maintenanceof the vehicle. Among all these costs, the fuel consumption cost is the majorcomponent that is sensitive to road slope, rolling resistance, and air drag.This cost is very significant for heavy trucks that are used in forestry ormining, which is our main consideration in this research.This thesis defines a multi-objective optimization model for designinga road that would be optimal for the manufacturers in terms of the roadconstruction cost and at the same time for the users in terms of the vehicleoperating costs, specifically in terms of the fuel consumption cost. It furtherpresents and examines a variety of methods for solving this multi-objectiveoptimization problem.11.2. Research approach$203 $266$82 $44 $44 $82$747$17$1,395$79$2,959$0$500$1,000$1,500$2,000$2,500$3,000AB BC MB NB NL NS ON PE QC SK CanadaVehicle operating cost ($ millions per year)Figure 1.1: Statistics showing the annual vehicle operating cost due to poorroads, by Canadian provinces and territories ($ millions per year) [CPC21].Note that these costs are not normalized by the population.One may wonder why we focus on the multi-objective optimization modelinstead of defining a single-objective optimization problem that aggregatesboth objective functions into one scalar function. Firstly, as the fuel pricevaries from region to region and varies over the lifetime of the mine, we decideto focus on the fuel amount used in mL, whereas the road construction costis estimated in dollars because it is a short-term event. Secondly, the totalfuel consumption cost also depends on the lifetime usage of the road whichis information that may not be available, for example, in a new mine wherethe exact amount of minerals are unknown. Therefore, for our problem,these two cost functions cannot be added together into a single-objectiveproblem.1.2 Research approachGenerally, the road design procedure can be considered as finding anoptimal road geometry for a given ground profile. Due to computationalbarriers, many studies split this process into horizontal alignment and verti-cal alignment [MR12, BHLH17]. In this thesis, we assume that the horizon-tal alignment is predetermined, and we consider only the vertical alignmentdesign that follows this fixed horizontal alignment. We optimize the verti-cal alignment by performing optimal earthwork operations to construct anoptimal road under some design requirements.In our research, we examine how road vertical alignment affects the road21.3. Thesis outlineuser costs, particularly the fuel consumption cost, and we model a multi-objective approach that takes in both the cost of road construction and thecost of road usage.In any optimization problem, a decision-maker first needs to identify thecore decision variables that will impact the overall quality of decisions. Theobjective functions represent this quality and are then used to make themost suitable decision. Thus, we first present the functional form of theroad construction cost and road user cost in the context of road verticalalignment design.Then, we review and examine several well-known multi-objective opti-mization methods. We also present an approach (called warm start strategy)for our problem that can improve the performance of some chosen methods.We test the performance of these selected methods on 30 road samples andattempt to determine the most promising approach for our model.1.3 Thesis outlineChapter 2 serves as a brief introduction to the road vertical alignmentdesign followed by the key concepts of multi-objective optimization necessaryfor this research.Chapter 3 formulates the functional forms of the road construction costand road user costs for our multi-objective optimization model in the contextof road vertical alignment design.Chapter 4 reviews several approaches for solving the proposed multi-objective problem. It provides brief descriptions of the algorithms alongwith their pseudocodes.Chapter 5 describes the experimental framework performed to assessthe performance of the selected methods and presents an analysis of thenumerical results. It examines the most promising and appropriate methodfor our proposed model.The final chapter, Chapter 6, summarizes all the main points of thestudy and highlights some suggestions for possible future work.3Chapter 2PreliminariesIn this chapter, we review some essential preliminary material to set thefoundation of this research.In Section 2.1, we present a brief introduction to the road vertical align-ment design, and in Section 2.2, we summarize the basic concepts of multi-objective optimization.2.1 Vertical alignment designDefinition 2.1 (Ground profile and road profile [MR12]). The ground pro-file is the vertical profile of the ground prior to the construction of the road.The road profile is the vertical profile of the ground after the constructionof the road.Whereas the horizontal alignment design examines the projection of thethree-dimensional ground profile on a horizontal plane, in the vertical align-ment design we focus on the projection on a vertical plane. In other words,the vertical alignment is a longitudinal cross-section of the horizontal align-ment. It is used to determine the pavement elevation, road grade, and otherroad characteristics.To model the vertical alignment, we split a road into smaller units knownas sections. We index them by the set S = {1, 2, . . . ,m}. We represent theroad with a spline which is defined as a piecewise quadratic function givenin the following formPi(x) = ai,1 + ai,2x+ ai,3x2, (2.1)where ai,1, ai,2, ai,3 are the coefficients of the quadratic polynomial spanningalong the section i ∈ S. Each polynomial piece is referred as a segment (alsoknown as a spline segment).For instance, Figure 2.1 illustrates a graphical example of the verticalalignment in profile view. The y-axis shows the ground elevation (or heightof the ground), while the x-axis represents the horizontal distance of theroad in meters. In this example, the road is divided into 5 sections and is42.2. Multi-objective optimizationSegment 1 Segment 2 Segment 3 Segment 4Ground profileRoad profileSection 1 Section 2 Section 3 Section 4 Section 5Horizontal distanceGround elevationFigure 2.1: Graphical illustration of vertical alignment.approximated using 4 (spline) segments. The ground profile and road profileare given with the blue and orange lines, respectively. The road profile isconstructed by changing the height of the ground elevation at different pointsand by moving the earth between sections.In other words, the input of the vertical alignment design is the groundprofile, and the output is the road profile provided as a quadratic spline.Note that the vertical alignment is determined for the fixed horizontal align-ment.For further details on the vertical alignment design, we refer to [MR12,BHLH17].2.2 Multi-objective optimizationMany real-life optimization problems in the agricultural, construction,and manufacturing sectors involve multiple objectives. These objectives areoften conflicting in nature in a way that a solution that minimizes one ofthe objectives might worsen another one. Therefore, all objectives needto be handled at the same time. This creates the field of Multi-ObjectiveOptimization (MOO).52.2. Multi-objective optimization2.2.1 Problem formulationWe consider a MOO problem given in the formminimize {f1(x), f2(x), . . . , fk(x)}subject to x ∈ X ,(2.2)where x ∈ Rn is a decision variable, fj : Rn → R for j ∈ J = {1, ..., k}with k ≥ 2 are objective functions, and X ⊂ Rn is called a feasible set.The image of the feasible set is called the objective function space.2.2.2 Pareto optimalityIn MOO, we need to select the ‘best’ solutions from the domain of feasiblesolutions. Thus, we need the concept of Pareto dominance to say whichsolution is better than others.Definition 2.2 (Pareto dominance [JMC09]). Consider the problem (2.2).Suppose that x,y ∈ X . Then we say that x Pareto dominates y, denotedby x ≺ y, if1. fj(x) ≤ fj(y) for all j ∈ J , and2. fj(x) < fj(y) for at least one j ∈ J .Example 2.3. Suppose that we have f = [f1, f2] where f1 and f2 are theobjective functions that are being minimized. And suppose that we havefour solutions a, b, c, and d, as illustrated in Figure 2.2. Here, c minimizesFigure 2.2: Pareto dominance example with four solutions a, b, c, and d,where a is a non-dominated solution (red point) and a, b, and d are domi-nated solutions (blue points).62.2. Multi-objective optimizationboth f1 and f2 better than b. Thus, both cases in Definition 2.2 hold, andc ≺ b. Also, c ≺ a, since c minimizes f2 better than a and minimizes f1as good as a. In similar way, c ≺ d, a ≺ b, and d ≺ b, while a and d areincomparable according to the definition. Among all these four solutions,there is no solution that Pareto dominates the solution c. Therefore, werefer to c as a non-dominated solution, whereas a, b, and d are dominatedones.We can now define the concept of solution in the multi-objective opti-mization context using the Pareto dominance relation.Definition 2.4 (Pareto optimal set [JMC09, ABC+20]). For a set of solu-tions P ⊂ X of the problem (2.2), Pareto optimal set P∗ is a set of solutionsthat are not Pareto dominated by any other solution of P, that is,P∗ = {x ∈ P | 6 ∃y ∈ P : y ≺ x}.Definition 2.5 (Pareto front [JMC09, ABC+20]). For the Pareto optimalset P∗ of the problem (2.2), Pareto front PF∗ is the image of the solutionsof the set P∗, that is,PF∗ ={(f1(x), f2(x), ..., fk(x)) | x ∈ P∗}.Instead of finding a single optimal solution as in a single-objective case,the aim of MOO is to estimate a Pareto optimal set that makes all theobjectives have as good outcomes as possible. It attempts to converge closerto the true Pareto optimal set and to discover a diverse set of solutions thatare well-distributed among themselves.In the next definition, we define some bounds of the Pareto front PF∗.Definition 2.6 (Ideal point and nadir point [JMC09]). For the Pareto op-timal set P∗ of the problem (2.2), the ideal point z∗ = (z∗1 , . . . , z∗k)T of thePareto front PF∗ is defined asz∗i = minx∈P∗fi(x)for i = {1, 2, . . . , k}. Similarly, the nadir point znad = (znad1 , . . . , znadk )T ofthe Pareto front PF∗ is defined asznadi = maxx∈P∗fi(x)for i = {1, 2, . . . , k}.72.2. Multi-objective optimizationFor instance, Figure 2.3 illustrates the Pareto front for the case of twoobjectives f1 and f2. The ideal point and nadir point are given with blueand red dots, respectively. The ideal point represents the simplest point thatPareto dominates the entire Pareto front, whereas the nadir point representsthe simplest point that is Pareto dominated by the entire Pareto front. Thesetwo points give a rectangular region that contains the entire Pareto front.𝑧∗𝑧!"#Pareto front𝑚𝑖𝑛𝑚𝑖𝑛𝑚𝑎𝑥𝑚𝑎𝑥Figure 2.3: Graphical illustration of the ideal and nadir points of the Paretofront for bi-objective problem.The notion of the ideal point is used in the descriptions of some se-lected MOO methods in Chapter 4, whereas the nadir point is used whencomputing the hypervolume indicator in Chapter 5.2.2.3 Pareto front filteringAlgorithm 2.1 outlines the pseudocode for identifying the Pareto optimalsolutions for the given solution set P.The main idea of the algorithm is to classify the solutions into two cat-egories, namely, dominated and non-dominated. The set I saves the labelsfor each solution of the set P in a way that it labels the dominated solutionswith 0 and the non-dominated ones with 1.At the beginning of the algorithm, we assume that all the elements of theset P are non-dominated. It starts with the first element of the set P, andat Steps 5-6 it identifies the solutions that are dominated by the selectedelement and marks them with a label 0. Then, at Step 7, it moves to the nextelement of the set P that has a label 1 and repeats the process. In otherwords, the algorithm selects the next element which is a ‘potential’ non-dominated solution and identifies the solutions that it Pareto dominates. Inthis way, in each iteration, we ignore the already found dominated solutions82.2. Multi-objective optimizationAlgorithm 2.1: Pareto optimal set filteringInput : Set of solutions P = {p1, p2, . . . , pn}Output: Pareto optimal set P∗1) Declare a set I of size n2) I = 13) i = 14) while i ≤ n do5) for j such that I(j) = 1 do6) if j 6= i and pi ≺ pj then I(j) = 0end7) Increment i until I(i) 6= 0end8) P∗ = P(I)and perform comparisons among the ‘potential’ non-dominated solutions.The algorithm stops when it reaches the last element of the set P, and itreturns the set of non-dominated solutions, that is, Pareto optimal set P∗.By applying Algorithm 2.1 on the set P, we find the Pareto optimal setcomposing the so-called Front 1. By removing the solutions of Front 1 fromthe set P and applying Algorithm 2.1 on the new set, we identify Front 2,and so on. This strategy is called a non-dominated sorting and is illustratedin Figure 2.4(a) [DPAM02]. Figure 2.4(b) shows an example of fronts forDominatedDominatedNon-dominatedNon-dominatedNon-dominatedNon-dominatedNon-dominatedNon-dominatedFront 1Front 2Front NFronts(a) Graphical illustration of non-dominated sorting.Front 1 (Pareto front)Front 3Front 4Front 2(b) Graphical illustration ofthe fronts.Figure 2.4: Graphical illustration of non-dominated sorting and fronts.92.2. Multi-objective optimizationbi-objective problem.The concept of the non-dominated sorting is used later in Chapter 4while describing the genetic algorithms selected for this thesis.10Chapter 3ModelIn this chapter, we present the objective functions of our MOO model,namely, road construction cost in Section 3.1 and vehicle operating cost inSection 3.2.3.1 Road construction costIn the road vertical alignment design, the aim of the road constructionprocess is to transform the given ground profile into a road profile usingearthwork operations. Thus, the effect of earthwork cost on the road con-struction cost is significant [Fwa89]. The earthwork cost consists of thefollowing three components.1. Excavation cost, denoted by Ccut, is the cost of cutting earth from asection;2. Hauling cost, denoted by Chaul, is the cost of transporting earth be-tween sections;3. Embankment cost, denoted by Cfill, is the cost of filling a section withearth.In this thesis, the road construction cost is retrieved from the Mixed-Integer Linear Programming (MILP) model with a single material pavementproposed by [MR12, BHLH17]. This cost function is considered as a black-box, that is, the internal structure and code implementation of the costfunction is not known. Below we present the general functional form for theroad construction cost.3.1.1 Variables and parametersDecision variablesA decision variable is a variable whose optimal value is determined dur-ing the optimization process. The decision variables in the road constructionproblem that are relevant to this research are listed below.113.2. Vehicle operating cost− ai,1, ai,2, ai,3 ∈ R: the coefficients of the quadratic polynomial definingthe spline spanning along section i ∈ S.Design parametersA design parameter is a predefined parameter used as an input for theoptimization problem. The design parameters are listed below.− pi ∈ R+: the cost per cubic unit of earth that is excavated fromsection i ∈ S.− cij ∈ R+: the cost per cubic unit of earth that is hauled fromsection i ∈ S to section j ∈ S (j 6= i).− qi ∈ R+: the cost per cubic unit of earth that is embanked tosection i ∈ S.3.1.2 FormulaThe road construction cost, denoted by Cconstruct, includes the total costof excavation, hauling, and embankment, as follows.Cconstruct = Ccut + Chaul + Cfill, (3.1)It estimates the total earthwork cost performed to construct a new road.For more details on the road construction cost, we refer to [MR12, BHLH17].3.2 Vehicle operating costVehicle operating costs consist of the daily costs spent by the users whileoperating vehicles on a given road profile. Fuel consumption, oil usage, tirewear, repairs, and maintenance of the vehicle are all examples of vehicleoperating costs that differ with vehicle use [TV03]. Among them, fuel con-sumption is the component that is most impacted by the vertical alignmentof the road, thus, is chosen as the second objective function of the model[JSJ06, TV03].3.2.1 Literature reviewOur goal is to derive an equation to quantify the fuel consumption costbased on road alignment. However, there is no clear relationship between123.2. Vehicle operating costfuel consumption and road alignment. A number of studies have been con-ducted to identify the effect of various parameters on fuel consumption. Webegin by discussing some of the previous methods to estimate fuel consump-tion.In 1987, it was proposed to estimate fuel consumption using only onevariable, namely, vehicle running speed [JSJ06]. After conducting severalexperiments, Jha [JSJ06] advanced the equation by including the road gradein calculations. Moreover, the author introduced a model for predicting fuelusage for the entire year taking into account the interest rate affecting thefuel price.Akc¸elik et. al. [ASB12] proposed the model for estimating the amountof fuel consumed by taking into consideration the physics laws. It relies onthe engine-generated tractive force, which is required to drive the vehicle,and considers the main external resistance forces affecting the vehicle. Themodel is analyzed using the software called SIDRA INTERSECTION and iscalibrated for different vehicle parameters [ASB12]. The user can estimatethe fuel consumption cost for any vehicle type by providing the model withthe needed vehicle parameters.In 2013, Swedish National Road and Transport Research Institute pro-posed a model for estimating the fuel consumption cost generated with ahelp of simulation software called VETO [CHE13]. Carlson et. al. [CHE13]simulated the traffic flow of different vehicle types for a systematic variationof roads and speed conditions. By generating a dataset, the influence ofroad variables on fuel usage was analyzed to produce the regression modelfor the fuel consumption cost. This model has a high degree of explanation,with R2 = 0.99. However, this model is overly complex and is only valid forspecific vehicle categories [CHE13].Another model created based on the collected data is introduced byKubler [Kub15]. Kubler [Kub15] collected the data for fuel usage in variousenvironmental settings for Caterpillar 785D off-highway trucks. The unique-ness of this model is that it depends only on one variable, namely on roadgrade. Thus, the model is very simple with the high coefficient of determi-nation of R2 = 0.89, and it considers the fuel usage for loaded trucks andunloaded trucks separately. Nevertheless, as the authors mentioned, thisregression model is valid only for the Caterpillar 785D off-highway trucks[Kub15].Soofastaei et. al. [SAKK16] outlines the model to evaluate an approx-imate amount of fuel burnt by a haul truck to move one tonne of minedmaterial in an hour. According to their model, the fuel consumption costdepends on the vehicle rimpull, which is defined as an amount of force ex-133.2. Vehicle operating costerted by tires to the road surface. As the authors mentioned, the vehiclerimpull can be estimated using the so-called Rimpull Table which is uniquefor all vehicle types [SAKK16].Another regression model with a coefficient of determination R2 = 0.84is outlined by Svenson and Fjeld [SF16]. They conducted an experimentwith a test truck to collect data for a different road grade, road curvature,surface roughness, and vehicle payloads in dry summer and wet fall weath-ers [SF16]. Then, the regression analysis was performed on the dataset toidentify the relationship between the fuel consumption cost and road param-eters. Nevertheless, the model is highly dependent on the surface roughnesswhich is not easily accessed in real life.Based on the investigations of different approaches to estimate fuel con-sumption cost, this thesis uses the model proposed by [ASB12] that can beadjusted for different environmental settings and can be applied for differ-ent vehicle types. Using this model, the user can estimate the fuel usage byproviding all the needed vehicle parameters with the environmental factorsaffecting the vehicle.3.2.2 Variables and parametersDecision variablesThe decision variables defining the fuel consumption cost are listed be-low.− ai ∈ R: the engine-induced acceleration (m/s2) for section i ∈ S.− si ∈ R: the road grade or slope (m/100 m) of the spline segmentspanning along section i ∈ S.− vi ∈ R+: the vehicle running speed (m/s) for section i ∈ S.Design parametersThe design parameters required for the vehicle operating cost can beclassified into two categories: vehicle parameters and environment parame-ters.Vehicle parameters are the predefined characteristics of the vehicle, suchas vehicle mass, horsepower, idle fuel rate, and others. The vehicle param-eters are listed below.− A ∈ R+: the projected frontal area (m2) of the vehicle.143.2. Vehicle operating cost− M ∈ R+: the vehicle total mass (kg).− Pe ∈ R+: the vehicle engine power (kW).− ptract ∈ R+: the percent mass acting on tractive axle (%) of the vehicle.− α ∈ R+: the constant idle fuel rate (mL/s) of the vehicle.− β1 ∈ R+: the fuel usage per unit of energy (mL/kJ) (also called as theenergy efficiency parameter).− β2 ∈ R+: the fuel usage per unit of energy-acceleration (also called asthe energy-acceleration efficiency parameter) (mL/(kJ.m/s2)).− η ∈ R+: the power transmission efficiency (%) of the vehicle.Environment parameters refer to the road resistance, gravity, air density,and other environmental factors affecting the vehicle motion. The environ-ment parameters are listed below.− cair ∈ R+: the air drag coefficient.− croll ∈ R+: the rolling coefficient depending on the road surface type.− ctire1, ctire2 ∈ R+: the rolling coefficient depending on vehicle’s tiretype such as bias ply or radial.− g ∈ R+: the acceleration due to gravity (m/s2).− r ∈ R+: the road resistance (%).− ∆di ∈ R+: the distance or length (m) of section i ∈ S.− µ ∈ R+: the coefficient of friction between tires and pavement.− ρair ∈ R+: the air density (kg/m3).3.2.3 FormulaThe paper by [ASB12] proposed the model for estimating the value offuel consumed (mL) for a simulation distance ∆di for section i ∈ S. Themodel has the following three components.1. Fuel consumption due to idling, denoted by Cidle,i, is an estimated fuelamount used to maintain engine operation at constant idle fuel rate αfor section i ∈ S;153.2. Vehicle operating cost2. Fuel consumption due to tractive force, denoted by Ctract,i, is an es-timated fuel usage to increase tractive force Rtract,i that moves thevehicle along section i ∈ S;3. Fuel consumption due to engine-induced acceleration, denoted byCaccel,i, is an estimated fuel amount used by the engine to acceler-ate the vehicle along section i ∈ S. Note that the engine-inducedacceleration refers to the acceleration produced while pushing the gaspedal of the vehicle, and it ignores the acceleration produced by theenvironmental factors such as gravitational force.For each section i ∈ S, the cost function Cfuel,i is defined asCfuel,i ={Cidle,i + Ctract,i + Caccel,i for ai > 0Cidle,i + Ctract,i for ai ≤ 0(3.2)whereCidle,i = α∆di/vi, (3.3)Ctract,i = β1 ∆diRtract,i/1000, (3.4)Caccel,i = β2M ∆di a2i /1000. (3.5)Note that the vehicle parameters α, β1, and β2 can be retrieved from[AB03, ASB12]. We provide sample values in Table A.1 and A.2 (in Ap-pendix A).For each section i ∈ S, we first need to estimate four parameters, namely,vehicle running speed vi, vehicle tractive force Rtract,i, engine-induced ac-celeration ai, and road grade si.Vehicle running speedThe vehicle running speed vi is a rate at which the vehicle travels thesection i ∈ S. With the purpose of simplicity of the model, we assume thatthe vehicle running speed is constant and predefined by the user.Vehicle tractive forceThe vehicle tractive force Rtract,i (N) is an amount of force required todrive the vehicle along the section i ∈ S. It is directly proportional to theengine power Pe generated by the vehicle and is inversely proportional tothe vehicle running speed vi, as follows.Rtract,i = min( 3600 η Pevi, Rmaxtract), (3.6)163.2. Vehicle operating costwhere Rmaxtractis the constant maximum tractive force (N) which can be pro-duced with no harm to the vehicle. Note that we include Rmaxtractin Equation(3.6) to ensure that the vehicle tractive force does not approach infinity atlow vehicle speeds [RLD+01]. Moreover, it can be quantified asRmaxtract= M g µptract.According to Equation (3.6), the vehicle tractive force is always posi-tive. This implies that the vehicle is moving without any stops. So, thismodel does not consider the case when the vehicle is parked with the enginerunning. We ignore this case because it indeed does not affect the roadgeometry.Engine-induced accelerationTo estimate the engine-induced acceleration, we assume that the vehicletractive force Rtract,i (N) is the sum of inertia force (N), air drag force Rair,i(N), grade resistance force Rgrade,i (N), and rolling resistance Rroll,i (N), asfollows [RLD+01]Rtract,i = M ai +Rair,i +Rgrade,i +Rroll,i, (3.7)where Rair,i = cair Aρair v2i /2,Rgrade,i = M g si,Rroll,i = M g croll (ctire1 + 3.6 ctire2 vi)/1000.For more details on the external resistance forces Rair,i, Rgrade,i andRroll,i, see [RLD+01].By rearranging Equation (3.7), the engine-induced acceleration ai is de-fined asai =Rtract,i −Rair,i −Rgrade,i −Rroll,iM.The engine-induced acceleration ai is positive when the vehicle tractiveforce Rtract,i is more than the external resistance forces. In other words,we increase the tractive force to overcome all external resistance forces ex-posed by nature to accelerate the vehicle. As a result, the engine inducesan acceleration and uses more fuel. However, when the produced tractiveforce by the vehicle is less than the external resistance forces, the vehicle isdecelerating, and so, no extra fuel is consumed by the engine. This processis visualized in Figure 3.1.173.2. Vehicle operating costFigure 3.1: Graphical illustration of tractive force and external resistanceforces.Road gradeThe road grade si is a slope of the spline segment spanning along thesection i ∈ S. Using Equation (2.1), the road grade si issi(x) =dPi(x)dx= ai,2 + ai,3x. (3.8)Total vehicle operating costAlong the given vertical alignment, the vehicle can move in two directionsas follows.1. Left-to-right traffic flow is the vehicle motion from left to right alongthe given vertical alignment;2. Right-to-left traffic flow is the vehicle motion from right to left alongthe given vertical alignment.The road vertical alignment, which is optimal for one of the directionsof the traffic flows, might be inefficient for the other direction of the trafficflow by consuming significantly more fuel. With this in mind, we need toconsider both directions of the traffic flow at the same time while optimizingthe total vehicle operating cost.Suppose that for the given vertical alignment, we label the left-to-righttraffic flow with sections S = {1, 2, . . . , n}, while for the right-to-left trafficflow we use its reverse order S ′ = {n, n − 1, . . . , 1}. In addition, suppose183.2. Vehicle operating costthat the road grade is si for section i ∈ S in the left-to-right traffic flow.Then, in the right-to-left traffic flow, the road grade will be −si for sectioni ∈ S.For instance, in Figure 3.2, we present a road profile with 5 sections andillustrate the road grades si for each section i ∈ {1, 2, . . . , 5} in two trafficflow directions.𝒔𝟏 𝒔𝟐 𝒔𝟑 𝒔𝟒 𝒔𝟓Left-to-right traffic flow (Loaded vehicle)Right-to-left traffic flow (Unloaded vehicle)Grade-𝒔𝟏-𝒔𝟐-𝒔𝟑-𝒔𝟒 -𝒔𝟓Horizontal distanceGround elevationFigure 3.2: Graphical illustration of traffic flows for the road with 5 sectionsand grades s1, s2, s3, s4, s5.In this model, we assume that the leftmost endpoint of the vertical align-ment (corresponding to Station 0) is a mining site or a forestry loading site.Thus, we force the ground profile and road profile to have the same groundelevations at the point associated with Station 0. With this in mind, wealso assume that the vehicle is loaded when moving from left to right, butis unloaded while moving from right to left.Ultimately, the total vehicle operating cost spent for traveling the wholealignment in both directions of the traffic flow isCfuel =∑i∈SCfuel,i +∑i∈S′Cfuel,i. (3.9)193.2. Vehicle operating costUnit testingIn order to ensure that the code for our model was high-quality, the unittesting was done along the way. The unit testing included the followingsmall roads:1. Linear flat ground with 5 sections;2. Linear inclined ground with 5 sections;3. Linear declined ground with 5 sections;4. Quadratic ground with a sag with 5 sections;5. Quadratic ground with a crest with 5 sections;6. Quadratic ground with a sag and crest with 5 sections;7. Quadratic ground with a crest and sag with 5 sections.ValidationThe current model is validated on the dataset found in [SAKK16]. Thedataset was collected from a surface coal mine with Caterpillar 793D off-highway trucks (CAT 793D). In the validation test, the vehicle parametersfor CAT 793D are retrieved from the official website of Caterpillar manufac-turer. As a result, the coefficient of determination for the proposed modelis R2 = 0.93 for this dataset.20Chapter 4MethodsIn this chapter, we discuss the MOO methods used in this thesis.In Section 4.1, we give a brief survey on some MOO techniques. InSections 4.2-4.5, we present the algorithm descriptions of the selected MOOmethods followed by the warm start strategy in Section 4.6 that can improvethe performance of some MOO methods.4.1 Literature reviewIn the literature, there are plenty of MOO techniques, and we start bydiscussing some of the well-known approaches for solving MOO problems.One of the widely used approaches is scalarization. Therein, the objec-tive functions are merged (or reformulated as constraints) to form a single-objective optimization problem which is then solved by known methods.Below we first discuss some of the classical and commonly used scalariza-tion techniques.The weighted sum method is one of the well-known and classical scalar-ization approaches. The method combines all the objectives into a scalarsingle objective by pre-multiplying each objective function with predefinedweights. A positive fact for this technique, which was first proposed by Gassand Satty [GS55], is its simplicity and easiness of implementation. However,the main limitation of the weighted sum method is the incapability of findingsolutions in nonconvex regions of the Pareto front.Goal programming is also one of the classical scalarization approaches,which was extended from linear programming and was firstly invented byCharnes and Cooper [CCF55, CC61]. The method does not attempt tooptimize the objectives directly, instead, it declares the objective functionsas the problem constraints restricted with the specific goals, and then itoptimizes deviations from these goals. Even though the method is verysimple and easy to implement, it highly relies on the choice of the goals whichare predefined by the user. Similar approaches that use this idea are thegoal attainment method, interactive goal programming, and lexicographicalgoal programming [JMC09].214.1. Literature reviewThe weighted metric method is another scalarization approach that issimilar to the weighted sum method but uses another means of aggregation ofthe objective functions [Asg98, Deb01]. It seeks the nearest possible solutionto a user-defined reference point, which is generally chosen as the ideal point,using either Manhattan, Euclidean, or Chebyshev metric. Depending on themetric type, the method can succeed in finding solutions in nonconvex partsof the Pareto front.Benson’s method is another example of a scalarization approach. It issimilar to the weighted metric method but it chooses a reference point to bea feasible non-Pareto optimal solution instead of the ideal point, and then itaims to find the farthest solution from the selected reference point [Deb01].This method can also uncover the nonconvex portion of the Pareto front. Inaddition, the method does not require the users to choose the reference point.Instead, it searches and selects the reference point on its own by definingsome extra constraints that restrict the search region. Nevertheless, thoseextra constraints may complicate the optimization process [Deb01].The ε-constraint method is also one of the commonly used scalarizationtechniques. This method was firstly proposed by Chankong and Haimes[CH83]. In this approach, one of the objectives is optimized while the othersare converted into the constraints which are restricted by some allowablethresholds ε. The main advantage of this method is that it works for bothconvex and nonconvex Pareto fronts, and it can better approximate thePareto front by varying ε values. On the other hand, in practice, it may bedifficult to specify the values of ε. We purpose one method to select ε valuesin Subsection 4.2.2 herein.The key benefit of all scalarization techniques is that once the given MOOproblem has been reformulated as a single-objective optimization problem,we can apply all known theoretical results and existing numerical algorithmsdedicated to this single-objective case. On the other hand, the scalarizationapproaches generally tend to find a single solution for each scalarization,thus, the quality of approximation of the Pareto front relies on the numberof scalarizations performed.The MOO approaches that can handle this drawback of scalarizationare multi-objective evolutionary approaches. For instance, NSGA-II uses agenetic algorithm for producing a set of solutions in each iteration, instead offinding a single solution. In this way, the method attempts to approximatethe entire Pareto front in a single run. There are a number of MOO methodsintegrated with evolutionary strategies, such as Vector Evaluated GeneticAlgorithm (VEGA) [Sch85], Strength Pareto Evolutionary Algorithm 2ndversion (SPEA-II) [EML01], Non-Dominated Sorting Genetic Algorithm 2nd224.1. Literature reviewversion (NSGA-II) [DPAM02], Non-Dominated Sorting Genetic Algorithm3rd version (NSGA-III) [DJ13], and Front Prediction based Non-DominatedSorting Genetic Algorithm 2nd version (FP-NSGA-II) [FTG15].Another drawback of the scalarization approaches is related to the het-erogeneous structure of the optimization problem. Consider the case whenone of the objectives is a computationally expensive function that requires alot of time to evaluate, while the other objectives are given analytically andare easy to assess. Then it would be beneficial to use this heterogeneousstructure of the optimization problem. However, scalarization techniquesmay break the heterogeneous structure of the given problem, leaving uswith a high computational effort that is mainly exposed by one expensiveobjective function.To take advantage of the heterogeneity, Thomann and Eichfelder [TE19]introduced a trust-region algorithm for a multi-objective setting. The methodemploys a basic trust-region approach by defining a local region within whichthe model can be trusted and the objectives can be approximated with suit-able models. In each iteration, the search direction is determined with ahelp of local ideal points [TE19]. Nevertheless, as the authors stated, thisapproach is developed for smooth MOO problems. Thus, it is not applica-ble for our proposed problem since the vehicle operating cost presented inChapter 3 is nonsmooth. A similar approach for the case when all of theobjectives are expensive functions is presented in [RK14], and another onefor non-expensive objectives with Taylor models is proposed in [VOS14].To handle the nonsmooth MOO problems, one can employ the multi-objective proximal bundle method introduced by Ma¨kela¨ et. al. [MKO14].The main idea is to transform all the objectives and constraints into oneobjective using a so-called improvement function. Then, it solves a nons-mooth unconstrained single-objective optimization problem using the prox-imal bundle method [MKO14]. It finds a descent direction that improvesall the objectives of the original problem at each iteration, and it gath-ers subgradients from the previous iterations into a bundle. In that way,it approximates the whole subdifferential instead of using only one arbi-trary subgradient at each point. Recently, Ma¨kela¨ and Montonen [MM20]advanced the multi-objective proximal bundle method with a new versionof the improvement function that can improve the algorithm performance.Nevertheless, as the authors stated, more numerical testing is needed.In this thesis, as a starting point of the MOO for vertical alignment de-sign, we focus on three classical scalarization methods (ε-constraint method,weighted sum method, and weighted metric method) and two multi-objectiveevolutionary methods (NSGA-II and FP-NSGA-II).234.2. ε-constraint method4.2 ε-constraint methodIn the ε-constraint method, we optimize one of the objectives while weuse the others as constraints that are bounded by some threshold ε. Forsimplicity, we consider a bi-objective case of the problem (2.2). We refor-mulate the given MOO problem to a single-objective optimization problemof the formminimize f2(x)subject to f1(x) ≤ ε,x ∈ X .(4.1)In the above formulation, the value of ε needs to be in the range betweensome user-defined lower bound f and upper bound f¯ of the objective f1, thatis, ε ∈ [f, f¯ ]. Since we aim to approximate the whole Pareto front, we choosef = minx∈Xf1(x) and f¯ = f1(x∗) where x∗ = argminx∈Xf2(x).In order to find a set of Pareto optimal solutions, we need to solve theproblem (4.1) for different values of ε ∈ [f, f¯ ]. While the idea is straight-forward, it raises a complex query. What values of ε should we consider tosucceed in approximating the entire Pareto front? In the following subsec-tions, we address this question.4.2.1 AlgorithmAlgorithm 4.1, inspired by [JMC09], outlines the pseudocode of the ε-constraint method for bi-objective optimization problem.Algorithm 4.1: Bi-objective ε-constraint methodInput : Objectives f1, f2 : Rn → R, objective bounds f, f¯ ∈ R,and increment δ ∈ R+Output: Set of solutions P1) Declare a matrix D with 4 columns2) ε = f3) while ε ≤ f¯ do4) Approximately solve the problem (4.1) with ε to find xε ∈ Rn5) D(i, : ) = [ ε, xε, f1(xε), f2(xε) ]6) ε = ε+ δend7) P = D( : , 2)244.2. ε-constraint methodThe general idea of the method is to solve the problem (4.1) for variousvalues of ε which are iteratively increased by a predefined δ. Note that thechoice of δ highly affects the performance of the algorithm. For instance,if δ is too small, the algorithm may converge to the same optimal solutionmultiple times, and thus, it may lead to a computational burden. On theother hand, if δ is too large, the algorithm may miss some Pareto-optimalsolutions. The algorithm presented in the next subsection seeks to addressthis drawback.4.2.2 Dynamic ε selectionIn Algorithm 4.2, we provide the pseudocode for the bi-objectiveε-constraint method with dynamic ε selection. The main idea of the dynamicε selection is to update the ε value by analyzing the obtained solutions. Thisis done in the following way.We create a matrix D to save all results found during the algorithmAlgorithm 4.2: Bi-objective ε-constraint method with dynamic εselectionInput : Objectives f1, f2 : Rn → R, objective bounds f, f¯ ∈ R,and tolerance tolε ∈ R+Output: Set of solutions P1) Declare a matrix D with 4 columns2) Approximately solve the problem (4.1) with ε = f to find xε ∈ Rn3) D(1, : ) = [ ε, xε, f1(xε), f2(xε) ]4) Approximately solve the problem (4.1) with ε = f¯ to find xε ∈ Rn5) D(2, : ) = [ ε, xε, f1(xε), f2(xε) ]6) i = 17) while not stopping criteria do8) ε = (D(i, 1) +D(i+ 1, 1))/29) Approximately solve the problem (4.1) with ε to find xε ∈ Rn10) D(end+ 1, : ) = [ ε, xε, f1(xε), f2(xε) ]11) Sort D in ascending order of D( : , 1)12) Find i such that norm(D(i, 3 : 4)−D(i+ 1, 3 : 4)) is maximumand D(i, 1)−D(i+ 1, 1) > tolε13) if i = ∅ then stopend14) P = D( : , 2)254.2. ε-constraint methodprocess: ε values and the corresponding solutions xε with the objectivefunction values f1(xε) and f2(xε).First, we need to determine the region where the true Pareto front liesusing the objective bounds f and f¯ . We solve the problem (4.1) with thevalue of ε set to these user-defined objective bounds and record the resultsinto the matrix D. We set i = 1 and define our search region to be the spacebetween D(i, 1) and D(i+ 1, 1).At each iteration, we select ε in the middle of the search region. Then,we optimize the problem (4.1) with the selected ε and record the outputinto D. By performing these steps, in each iteration, we divide one searchregion into two ‘potential’ search regions from which we select the next ε.To identify the next search region, we sort the matrix D in ascending or-der of its column that stores the values of ε. Then, we find the biggest regioni where no solution is found yet and which is greater than the predefinedtolerance tolε. If no such region is found, then the algorithm terminates. Forour research, we choose the tolerance tolε to be 5% of the region betweenthe objective bounds f and f¯ .We repeat this process until some stopping criteria are met. These stop-ping criteria are defined by the user. For example, it might be the maximumnumber of found solutions or some other tolerance restrictions.To show the effect of the dynamic selection approach, we present anexample in Figure 4.1 that illustrates the approximated Pareto fronts for theTNK test problem that is presented in [TWFT95]. It visualizes the resultsfound at each iteration by the ε-constraint method (Figure 4.1(a)) and theε-constraint method with dynamic ε selection (Figure 4.1(b)). The lightgreen dots are the solutions of the true Pareto front, the red dot representsthe value of ε from the current iteration, and the blue dots correspond tothe values of ε from the previous iterations.Watching Figure 4.1, the ε-constraint method spends several numberof iterations searching in the region where no solutions of the true Paretofront exist. To put it in another way, the dynamic ε selection approachhelps the algorithm to identify the regions where predominantly to searchfor solutions.In this thesis, we use the ε-constraint method with the dynamic ε se-lection approach due to its effectiveness in the ε selection. In the comingchapters, for simplicity, we refer to this method as the ε-constraint method.264.3. Weighted sum method0 0.2 0.4 0.6 0.8 100.20.40.60.8100.20.40.60.8100.20.40.60.8100.20.40.60.8100.20.40.60.8100.20.40.60.8100.20.40.60.81Number of iterations: 27Number of solutions: 21𝜀-constraint method(a) ε-constraint method.0 0.2 0.4 0.6 0.8 100.20.40.60.810 0.2 0.4 0.6 0.8 100.20.40.60.8100.20.40.60.8100.20.40.60.8100.20.40.60.8100.20.40.60.8100.20.40.60.8100.20.40.60.810 0.2 0.4 0.6 0.8 100.20.40.60.81Number of iterations: 26Number of solutions: 24𝜀-constraint method with dynamic 𝜀 selection(b) ε-constraint method with dynamicε selection.Figure 4.1: TNK problem solved by ε-constraint method and ε-constraintmethod with dynamic ε selection.4.3 Weighted sum methodWe consider the problem given in (2.2) with k objectives. In the weightedsum approach, we aggregate all the objectives by pre-multiplying each ob-jective fi by user-defined weights wi and form the following problemminimizek∑i=1wi fi(x)subject to x ∈ X ,(4.2)where wi ≥ 0 for all i = {1, . . . , k} and∑ki=1wi = 1.As the ε-constraint method, this approach is very simple. To show howthe weighted sum method performs to find a solution on the Pareto front,we consider a bi-objective case of the problem (2.2). Then, the weightedsum method aims to minimizey = w1 f1(x) + w2 f2(x)subject to some given constraints. As y is a linear combination of the ob-jectives f1 and f2, we can represent it in the objective function space asa straight line which has the slope defined by the weights w1 and w2 as274.3. Weighted sum method(a) Convex case. (b) Nonconvex case.Figure 4.2: Graphical illustration of the weighted sum method for convexand nonconvex Pareto fronts.shown in Figure 4.2(a). By varying the weight combinations, we can repre-sent different y lines that touch the Pareto front, and for the convex Paretofront, we have enough room to quantify such points with different weights.However, for the nonconvex case presented in Figure 4.2(b), we can see thatthere are points of the nonconvex region of the Pareto front that cannot bereached for any combinations of the weight.4.3.1 AlgorithmAlgorithm 4.3 presents the pseudocode for the weighted metric methodfor the bi-objective optimization problem.The algorithm idea is similar to the ε-constraint method. It finds aset of Pareto optimal solutions by iteratively updating the weights by afixed increment δ. As before, this approach of choosing the weights mayincrease the computational time. To alleviate this drawback as in Subsection4.2.2, we introduce the dynamic weight selection which assists the algorithmwith the choice of weights. We discuss the proposed approach in the nextsubsection below.4.3.2 Dynamic weight selectionAlgorithm 4.4 outlines the pseudocode of the weighted sum method in-tegrated with the dynamic weight selection approach. The main idea of thealgorithm is the same as in Subsection 4.2.2.284.3. Weighted sum methodAlgorithm 4.3: Bi-objective weighted sum methodInput : Objectives f1, f2 : Rn → R and increment δ ∈ (0, 1)Output: Set of solutions P1) Declare a matrix D with 4 columns2) w1 = 0 and w2 = 1− w13) w = (w1, w2)T4) while w1 ≤ 1 do5) Approximately solve the problem (4.2) with w to find xw ∈ Rn6) D(i, : ) = [ w, xw, f1(xw), f2(xw) ]7) w1 = w1 + δ and w2 = 1− w18) w = (w1, w2)Tend9) P = D( : , 2)The procedure determines which weight combinations should be con-sidered in order to cover a representative portion of the Pareto front. Aspreviously, this is accomplished by archiving and analyzing all of the resultsobtained during the algorithm’s execution. We record the weight vector w,the corresponding solutions xw, and the objective function values f1(xw)and f2(xw).Firstly, the method determines the bounds where each of the objectivesis optimal using the weights w = (0, 1)T and w = (1, 0)T . In this way, weidentify the region where the true Pareto front lies. Then, it follows the sameprocedure presented in Subsection 4.2.2. It chooses the weights dependingon the region where no solution is found yet.The algorithm stops either when it completes the search for solutions orwhen it satisfies some stopping condition defined by the user. Note that forthis thesis, we select the tolerance tolw to be 5% of the allowable range forthe weight w, which are 0 and 1. In other words, we choose tolw to be equalto 0.05.Similar to the dynamic ε selection approach presented in Subsection4.2.2, the proposed method can also perform better than the weighted summethod outlined in Algorithm 4.3. Thus, as before, we choose to use theweighted sum method with the dynamic weight selection due to its effec-tiveness in weight adjustments. In the coming chapters, we simply refer tothis method as the weighted sum method.294.4. Weighted metric methodAlgorithm 4.4: Bi-objective weighted sum method with dynamicweight selectionInput : Objectives f1, f2 : Rn → R and tolerance tolw ∈ (0, 1)Output: Set of solutions P1) Declare a matrix D with 4 columns2) Approximately solve the problem (4.2) with w = (0, 1)T to findxw ∈ Rn3) D(1, : ) = [ w, xw, f1(xw), f2(xw) ]4) Approximately solve the problem (4.2) with w = (1, 0)T to findxw ∈ Rn5) D(2, : ) = [ w, xw, f1(xw), f2(xw) ]6) i = 17) while not stopping criteria do8) w = (D(i, 1) +D(i+ 1, 1))/29) Approximately solve the problem (4.2) with w to find xw ∈ Rn10) D(end+ 1, : ) = [ w, xw, f1(xw), f2(xw) ]11) Sort D in ascending order of D( : , 1)12) Find i such that norm(D(i, 3 : 4)−D(i+ 1, 3 : 4)) is maximumand D(i, 1)−D(i+ 1, 1) > tolw13) if i = ∅ then stopend14) P = D( : , 2)4.4 Weighted metric methodWe consider the problem (2.2) with k objectives. In the weighted metricmethod, we aim to find the closest feasible solution to an ideal point z∗using the lp-metric for p ∈ [ 1, ∞ ]. For p ∈ [ 1, ∞ ), this is done by solvinga single-objective optimization problem of the formminimize( k∑i=1wi |fi(x)− z∗i |p)1/psubject to x ∈ X ,(4.3)and for p =∞, by solving the problem of the formminimize maxi={1,...,k}{wi |fi(x)− z∗i |}subject to x ∈ X ,(4.4)304.4. Weighted metric methodwhere wi ≥ 0 for all i = {1, . . . , k} and∑ki=1wi = 1.The method can perform differently depending on the lp-metric. Thewell-known metric types are Manhattan metric (p = 1), Euclidean metric(p = 2), and Chebyshev metric (p =∞).To show how the metric types differ from each other and how they affectthe performance of the weighted metric method, we consider the bi-objectiveoptimization problem. Figure 4.3 visualizes how the weighted metric methodwith different metric types works.(a) Manhattan metric. (b) Euclidean metric. (c) Chebyshev metric.Figure 4.3: Graphical illustration of the weighted metric method with dif-ferent metric typesAs can be seen from Figure 4.3(a), the weighted metric method withp = 1 is similar to the weighted sum method discussed in Section 4.3, andthus, it also cannot approximate the points of the nonconvex region of thePareto front [JMC09]. With p = 2, the method can perform better comparedto p = 1 as presented in Figure 4.3(b). On the other hand, it increases thedegree of the polynomial that is being optimized in the problem (4.3), andthis may complicate the optimization process. Lastly, the weighted metricmethod with p =∞ can reach the solutions lying on the nonconvex portionof the Pareto front as illustrated in Figure 4.3(c). However, note that thereformulated problem given in (4.4) is nonsmooth.4.4.1 AlgorithmAlgorithm 4.5 presents the pseudocode of the weighted metric methodwith the lp-metric. The procedure is very similar to the weighted summethod outlined in Algorithm 4.3. The only difference is that the userneeds to input the metric type with which the method approximates thePareto front.314.5. Genetic algorithmsAlgorithm 4.5: Bi-objective weighted metric methodInput : Objectives f1, f2 : Rn → R, increment δ ∈ (0, 1), andmetric type p ∈ [ 1, ∞ ]Output: Set of solutions P1) Declare a matrix D with 4 columns2) w1 = 0 and w2 = 1− w13) w = (w1, w2)T4) while w1 ≤ 1 do5) Approximately solve the problem (4.3) with w and lp-metric tofind xw ∈ Rn6) D(i, : ) = [ w, xw, f1(xw), f2(xw) ]7) w1 = w1 + δ and w2 = 1− w18) w = (w1, w2)Tend9) P = D( : , 2)4.4.2 Dynamic weight selectionLike other scalarization techniques, in order to find several solutions ofthe Pareto optimal set, we need to run the weighted metric method withdifferent weight combinations. Unfortunately, this can produce some compu-tational effort. Therefore, as before, we apply the dynamic weight selectionapproach to the weighted metric method. This approach follows the samelogic discussed in Subsection 4.3.2, and thus, we skip the pseudocode forthis method.In this thesis, we use the weighted metric method with dynamic weightselection to which in the coming chapters we simply refer as the weightedmetric method.4.5 Genetic algorithmsThe genetic algorithm (GA) is an optimization method based on theconcepts of natural selection and genetics. The main idea is to imitatean evolutionary process on a set of feasible solutions. It removes the poorsolutions and generates a set of better solutions. The pseudocode is retrievedfrom [AH17] and is presented in Algorithm 4.6.In GA, solutions are referred as individuals, while we call a set of solu-tions as a population. Moreover, we call new individuals created during the324.5. Genetic algorithmsAlgorithm 4.6: GAInput : Objective f : Rn → R, a feasible region X ∈ Rn, an initialpopulation P 0 = {x1,x2, . . . ,xN}, and mutationprobability γ ∈ (0, 1)Output: Final population P final1) t = 02) while not stopping criteria do3) Evaluate the fitness of all x ∈ P t using f(x) . Fitness4) Select 2 parents from P t and go to 5) . Selectionor select 1 survivor from P t and go to 7)5) Generate a new offspring xnew from 2 parents . Crossover6) Mutate xnew with the probability γ . Mutation7) if xnew 6∈ X thenReturn to 4)elseAdd xnew to P t+1if |P t+1| < N then return to 4)end8) t = t+ 1end9) P final = P tevolutionary process as offspring.Given an initial population P 0, GA begins by assigning each individuala fitness score that is used to select survivors and parents. The fitnessdetermines how strong an individual x is based on the function value f(x).The likelihood that an individual will be selected to reproduce the nextpopulation depends on the individual’s fitness score.The randomly selected parents create offspring using the crossover. Asthe offspring in nature are never exact copies of their parents, they go underthe mutation step to create more variation. Then, the feasible offspring areadded to the next population P t+1, while the infeasible ones are rejected andwe go back to the selection step. Note that the child population P t+1 musthave as many individuals as its parent population P t. This evolutionaryprocess stops when some stopping condition is satisfied, and the algorithmreturns the final population with better solutions. A common stoppingcriterion is to set the maximum number of generations or to set the maximumfunction evaluations [AH17].334.5. Genetic algorithmsThe GA is very flexible so that the user can choose crossover and muta-tion techniques. Some of them can be found in [AH17]. In this thesis, theversion of GA that we used was NSGA-II taken from [DPAM02], and allparameters were left at their default values. The details of the parameterconfigurations are given in Section 5.2.In the GA, the definition of fitness is indeed very flexible. For a single-objective optimization problem, the fitness of the given population can berepresented by the function values evaluated for this population. In themulti-objective case, we seek to quantify the fitness score in a way thatfor the given population it shows the goodness of an individual for all theobjectives at the same time. Therefore, in the coming GA methods, we usethe non-dominated sorting approach presented in Chapter 2 for quantifyingthe fitness score of the individuals.4.5.1 NSGA-IIOne of the well-known benchmark GA methods developed for the multi-objective setting is NSGA-II which has succeeded in achieving the uniformlydistributed and well-approximated Pareto front for bi-objective problems[FTG15]. The main goal of NSGA-II is to find a set of Pareto-optimalsolutions using the GA and non-dominated sorting. In Algorithm 4.7, weAlgorithm 4.7: NSGA-IIInput : Objectives f = {f1, f2, . . . , fk}, a feasible region X ∈ Rn,an initial population P 0 := {x1,x2, . . . ,xN}, andmutation probability γ ∈ (0, 1)Output: Final population P final of non-dominated solutions1) t = 02) while not stopping criteria do3) Generate a new population Qt using GA(f,X , P t, γ)4) Rt = P t ∪Qt5) Apply non-dominated sorting on Rt6) Calculate crowding distances for all x ∈ Rt7) Generate P t+1 by choosing N solutions of Rt based on thefitness and crowding distance8) t = t+ 1end9) P final = P t344.5. Genetic algorithmsprovide the pseudocode for NSGA-II [DPAM02].To better understand the algorithm, we discuss Steps 3-7 in more detailbelow.At Step 3, we create a child population Qt from the parent populationP t using the GA. Then, at Step 4, the parent population P t and childpopulation Qt are combined to produce a population Rt of size 2N . In thisway, we preserve the best solutions of the current population in the nextgeneration.At Step 5, we evaluate the fitness of all individuals of the population Rtby ranking them based on the Pareto dominance as stated above. Then, atStep 6, the crowding distance is calculated for each individual of the pop-ulation Rt. As shown in Figure 4.4, the crowding distance of an individualpi is a perimeter of the rectangle formed by two closest neighbors pi−1 andpi+1 in the objective function space. It shows how close an individual is toits neighbors in the objective function space. Note that for the less crowdedsolutions, the crowding distance is higher. Therefore, in order to reach abetter diversity in solution, the solutions with higher crowding distances arepreferable.𝑝!"#𝑝! 𝑝!$#Figure 4.4: Graphical illustration of the crowding distance.Lastly, at Step 7, we generate the next population P t+1 by selecting onlythe best N solutions of Rt based on the rank and highest crowding distances.NSGA-II terminates when some stopping condition is satisfied and returnsthe final population of non-dominated solutions.4.5.2 FP-NSGA-IIThe paper by [PF03] showed that the performance of NSGA-II decreasesfor the optimization problems with higher dimensions and more variables.Thus, Fettaka et. al. [FTG15] proposed a new hybrid algorithm that applies354.5. Genetic algorithmsNSGA-II with the prediction algorithm FP-NSGA-II to improve convergenceto the true Pareto optimal set. The pseudocode is given in Algorithm 4.8[FTG15].The method is similar to NSGA-II, except it has an additional step forthe predictive operator to produce better offspring. This step is called frontprediction. The general idea of the front prediction is to identify the direc-tion in the decision variable space toward the true Pareto optimal solutions[FTG15]. This is done in the following way.Given the population P t, we apply the Pareto front filtering, which wasdiscussed in Chapter 2, on P t and find the Pareto optimal set, which formsthe solutions of the first front F t1. In the same way, by applying the Paretofront filtering once again on P t \F t1, we identify the solutions on the secondfront F t2.Algorithm 4.8: FP-NSGA-IIInput : Objectives f = {f1, f2, . . . , fk}, a feasible region X ∈ Rn,an initial population P 0 := {x1,x2, . . . ,xN}, andmutation probability γ ∈ (0, 1)Output: Final population P final of non-dominated solutions1) t = 02) while not stopping criteria do3) Determine Front in P t. Label this F t14) Create P˜ t = P t \ F t15) Determine Front in P˜ t. Label this F t2. (Note, if P˜t = ∅, thenF t2 = ∅)6) if F t2 6= ∅ then7) foreach x ∈ F t1 do8) Determine y∗ = argminy∈F t2 norm(f(x)− f(y))9) Calculate step = x− y∗10) xpred = x + step11) Add xpred to Ptendend12) P t+1 = NSGA-II(f ,X , P t, γ)13) t = t+ 1end14) P final = P t364.5. Genetic algorithmsIf the second front F t2 is not empty, then at Step 8, for each solutionx in F t1, we determine the solution y∗ in F t2 which is the nearest to xin the objective function space. Then, at Steps 9-10, we calculate a vectordifference between x and y∗ and refer to it as a step direction with which weextrapolate the predicted solutions xpred. At Step 11, we add the predictedsolution xpred to the population Pt.Then, we continue with the NSGA-II applied to the population P t.Note that, Steps 6-12 of the algorithm increase the population size of P t.In other words, suppose that for each t, we generate N tpred number of thenew predicted points. Then, the size of the population P t that is passedto NSGA-II is N + N tpred. This implies that FP-NSGA-II has Ntpred moref-calls than NSGA-II in each iteration.Lastly, Algorithm 4.8 terminates when some predefined stopping condi-tion is satisfied and returns the final population of non-dominated solutions.To give some idea of how NSGA-II and FP-NSGA-II are performing, weimitate one simulation from the paper by [FTG15] and present in Figure4.5 that illustrates the Pareto fronts found by the methods for the TNKproblem. The light green dots are the true Pareto front, while the resultsfor NSGA-II and FP-NSGA-II are given by a red cross symbol and bluecircles, respectively.Watching Figure 4.5, we can see that FP-NSGA-II performs better thanNSGA-II in every generation. This example supports the results presentedin the paper by [FTG15] that compares both methods for a fixed number ofgenerations.Figure 4.5: TNK problem solved by NSGA-II and FP-NSGA-II.374.6. Warm start strategy4.6 Warm start strategyReferring to the scalarization techniques discussed in Sections 4.2-4.4,we notice that all of them solve a new single-objective optimization problemin each iteration. Depending on the original problem, these single-objectiveproblems can be classified as smooth/nonsmooth, convex/nonconvex, lin-ear/quadratic/nonlinear, differentiable/nondifferentiable, etc. Based on thisinformation, the user can use any appropriate optimization solver.Generally, before the optimization solver starts an optimization process,the user can set a starting point that indicates the initial values for thedecision variables. This starting point represents the best guess where weexpect to find an optimal point. If no starting point is provided by the user,then the optimization solver takes a default starting point.For most optimization methods, the better the starting point, the fasterthe convergence. With this in mind, we implement the warm start strategy.The main idea of the warm start strategy is to provide the optimizationsolver with ‘good’ starting points so that it can converge to the optimalsolution faster. Figure 4.6 visualizes the warm start strategy for the verticalalignment setting.Consider the bi-objective optimization problem with objectives Cconstructand Cfuel proposed in Chapter 3. Recall that in the vertical alignmentdesign we input a ground profile to construct the road profile as an output.Moreover, as mentioned in the previous sections, the first step of all of theselected scalarization techniques is to identify the region where the truePareto front lies by optimizing each objective separately. So, as shownin Figure 4.6, we start by minimizing one of the objectives, Cconstruct, fora given ground profile and some default starting point that is set by theGround profileRoad profileDefault starting pointRoad profileYES NOStopping criteria?ARCHIVEof starting pointsRetrieve a starting pointOutput:ARCHIVE of starting pointsmin 𝐶!"#$%&'!% min {𝐶!"#$%&'!%, 𝐶(')*}Starting pointGround profileFigure 4.6: Graphical illustration of the warm start strategy.384.6. Warm start strategyoptimization solver. Then, we save the resulting road profile to use it lateras a starting point for the next optimization problem.We record all road profiles found during the algorithm process into anarchive of starting points. If the stopping criteria defined by the user arenot met, then we retrieve one starting point from this archive in order to useit as a warm start in the next optimization procedure. In particular, for theε-constraint method, our choice of the starting point depends on the valueof the current ε. We select the road profile that was obtained for the ε valuethat is closest to the current ε. Similarly, for the weighted sum methodand weighted metric methods, we use the weight values to determine thepotential starting point.With the selected starting point and the given ground profile, we moveto the next optimization problem and save the resulting road profile into thearchive. In this way, we use a warm start in each iteration of the scalarizationtechniques to select the potentially ‘good’ starting points that can acceleratethe convergence process.39Chapter 5Numerical resultsIn this chapter, we describe the numerical experiments that were per-formed to analyze the efficiency of the chosen MOO methods for solving theproposed bi-objective optimization problem.In Section 5.1, we discuss the well-known performance indicators that areused to compare the results found by the MOO methods. Then, in Section5.2, we outline how the experiment is conducted and give some overview onthe experimental setup and test problems. In Section 5.3, we present all thenumerical results analyzed by the selected performance indicators.5.1 Performance indicatorsTo identify the most promising technique for solving the proposed MOOproblem, we want to compare the performance of the selected methods in aquantitative manner. In the single-objective minimization problem, we saythat the method with the smaller objective function value is better. In themulti-objective case, it is not straightforward what the appropriate qual-ity measures are for the approximated Pareto fronts. In order to compareseveral sets of the solutions found by different MOO methods, we use perfor-mance indicators that attempt to quantify the quality of the found solutionset using a scalar value. In the coming subsections, we review several perfor-mance indicators found in the literature and then discuss the selected onesmore closely.5.1.1 Literature reviewThere are numerous studies conducted to develop the performance in-dicators that aim to assess the quality of Pareto front estimates. Theyattempt to examine how close the approximation of the Pareto front is tothe true Pareto front and how diverse the found solutions are. As the au-thors [ZTL+03, CSQ12] mentioned, there is no performance indicator thatsurpasses other indicators and is universal for all the MOO problems. Below405.1. Performance indicatorswe discuss some of the well-known and widely-used performance indicatorsalong with their strengths and weaknesses.Hypervolume is the most well-known performance indicator introducedby Zitzler [Zit99]. It measures the volume of space covered by the points ofthe obtained Pareto front and bounded by some predefined reference pointin the objective function space. A positive fact for the hypervolume is thatit can capture the information about the convergence and diversity of thefound solution set without knowing the true Pareto optimal set. On theother hand, the measure is sensitive to the choice of the reference point[CSQ12]. We attempt to address this for our problem in Subsection 5.1.2.Generational distance [VV99] is another commonly used performanceindicator that evaluates how close the found Pareto optimal solutions are tothe true Pareto optimal set. It calculates the average Euclidean distance inthe objective function space from each obtained solution to the closest pointon the true Pareto front which must be known beforehand [ABC+20]. Aweakness of this indicator becomes apparent by considering the case whenthe method finds only one solution and this solution lies on the true Paretofront. Then, the generational distance will be 0 implying that the foundresult is an ideal Pareto front approximation, which is clearly not correct.Inverted generational distance [CC05], averaged Hausdorff distance[SELC12], and modified inverted generational distance [IMTN15] are similarvariations of the generational distance that can handle its drawback but stillrequire the true Pareto front to be known.Spread, also known as the ∆-metric [DPAM02], is also a common per-formance indicator. It is a measure showing the diversity and uniformityof the solutions along the true Pareto front. However, the spread’s equa-tion proposed by [DPAM02] is applicable only for bi-objective problems. In[ZJZ+06], the spread indicator is advanced for problems with more than twoobjectives. The main disadvantage of this performance indicator is that itneeds the true Pareto optimal set to be predetermined before the assessment.Similar to the spread, one can use another performance indicator calledspacing [Sch95, DPAM02]. Spacing also analyzes the diversity and unifor-mity of the found Pareto optimal set in the objective function space butdoes not need any information about the true Pareto front. Nevertheless,as it only considers the distance between a point and its neighbor, it doesnot assess the convergence aspect, that is, it does not examine how close thefound solutions are to the true Pareto optimal set.Another common performance indicator is the final number of non-dominated solutions. As its name suggests, it counts the number of non-dominated solutions of the found solution set and does not need the infor-415.1. Performance indicatorsmation about the true Pareto optimal set. However, this indicator mightbe inapplicable in general [ABC+20]. For instance, suppose that we havemethod A that determines an approximated Pareto optimal set with thehundred non-dominated solutions, whereas another method B finds onlyone solution that can Pareto dominate all the solutions found by A. Then,according to the final number of non-dominated solutions, we reach an in-correct conclusion that method A is better than B.Error ratio is another performance indicator that deals with the car-dinality of the found Pareto optimal solutions [VV99]. It identifies whatportion of the found solutions belongs to the true Pareto optimal set. Thisindicator may suffer from rounding errors during the classification process.Moreover, the true Pareto front must be available in advance.Lastly, the CPU time plays a crucial role while comparing the perfor-mance of the optimization methods. It records the time required to find theapproximation of the Pareto optimal set with appropriate quality. Note thatthe quality of the solutions can be defined by the above-stated performanceindicators.More information on other performance indicators can be found in theexisting surveys [ZTL+03, CSQ12, ABC+20].In this thesis, we employ the performance indicators that can be evalu-ated for unknown true Pareto front: hypervolume, spacing, and CPU time.The reason why we choose these indicators is that using the CPU time wecan identify the fastest MOO method, and using the other two indicatorswe can define the conditions for good quality results. While the spacingindicator can purely examine the diversity and uniformity of the obtainedsolutions, the hypervolume indicator can assess both convergence and diver-sity of the results. In this way, we attempt to determine the most promisingMOO method that can find relatively good results for the proposed MOOproblem in road vertical alignment design.In the next subsections, we discuss the selected performance indicatorsin more detail.5.1.2 HypervolumeHypervolume is a measure of the region in the objective function spacedominated by the Pareto optimal solutions pi and bounded from top bya reference point as shown in Figure 5.1 [FTG15, ABC+20]. We choosethe reference point to be the point slightly worse than the nadir point ofthe obtained Pareto front. For consistency in calculations and fairness incomparisons, the reference point must be the same for all the compared425.1. Performance indicators𝑝!Reference pointHypervolume𝑝"𝑝#𝑝$Figure 5.1: Graphical illustration of the hypervolume indicator for bi-objective problem. Hypervolume is the area of the shaded region.MOO methods. Also, note that the reference point is not the nadir pointitself, but is related to it. If we choose the nadir point as the referencepoint for all the compared methods, then the space covered by the foundedge points would be 0, that is the hypervolume value for those edge pointswould be 0. This implies that the method could not find those solutionsat the edges of the Pareto front, which is actually not true. This is themain reason why we do not choose the nadir point as the reference point forcalculating the hypervolume indicator.Hypervolume indicator can tell the user about the closeness of the foundsolutions to the true Pareto front and, at the same time, about the diversityof the solutions [FTG15, ABC+20]. Moreover, there is no need to know thetrue Pareto front to calculate this metric. The MOO method with a higherhypervolume is considered to be better than the others.Suppose that we optimize a small road with 22 sections using theε-constraint method and weighted sum method. As a result, we obtainthe Pareto fronts shown in Figure 5.2(a). Choosing the reference point asthe point that is 5% worse than the nadir point, the hypervolume indicatorsare 3.33 × 105 for the ε-constraint method and 3.16 × 105 for the weightedsum method. Visually and numerically, we can conclude which method per-forms better. Nonetheless, these values for hypervolume are quite large.This indeed may complicate the analysis process later.To alleviate these issues, we normalize the obtained solutions into1-to-1 box as shown in Figure 5.2(b). Then, by selecting the reference pointat a point (1, 1)T , we find that the hypervolume indicators are 0.80 for theε-constraint method and 0.75 for the weighted sum method. Note that thehypervolume for an ideal case will be 1. In this way, we can easily identify435.1. Performance indicators(a) Original objective function space. (b) Normalized objective functionspace.Figure 5.2: Example of the found Pareto fronts for the road with 22 sections.what portion of the 1-to-1 box is covered without visualizing the results,and we can come to the same conclusion on which method performs betterfor the given test problem.Note that while normalizing the results into the 1-to-1 box, we do not putthe solutions located at the edges at the points (1, 0)T and (0, 1)T , instead,we keep 5% of space farther from the edges of the plot. This is done in orderto track whether the MOO method could identify those solutions or not.If we put them at the corners and the method indeed finds them, then thespace covered by those solutions would be still 0.There are several approaches to estimate the hypervolume indicator forthe given Pareto front. In this thesis, we calculate the hypervolume indicatorusing the algorithm proposed by [Fle03] and coded by [Kru11].5.1.3 SpacingSpacing is a measure showing how evenly and uniformly the resultingPareto optimal solutions are distributed among themselves. For the approx-imated Pareto front PF∗, the spacing indicator ∆′ is calculated as follows[JOZF14, ABC+20].445.1. Performance indicators∆′ =|PF∗|−1∑i=1|di − d¯||PF∗| − 1 , (5.1)where di are the Euclidean distances between consecutive solutions in theobjective function space as shown in Figure 5.3, and d¯ is the average of alldi for i ∈ [ 1, |PF∗| − 1 ]. Equation (5.1) is applicable only for bi-objectiveoptimization problems. For more than two objectives, refer to [ZJZ+06].𝑑! 𝑑" 𝑑# 𝑑$Figure 5.3: Graphical illustration of the spacing indicator for bi-objectiveproblem. d1, d2, d3, and d4 are the Euclidean distances between consecutivesolutions.Note that the spread indicator for the ideally distributed dataset is 0,while for the clustered distribution the value increases. Moreover, this indi-cator does not give any information about how good the found results arecompared to the true Pareto front. It just assesses how well the obtainedresults are distributed among themselves.5.1.4 CPU timeCPU time (sec.) is a measure indicating an approximate performancetime required to obtain an estimate for the Pareto optimal set with somegood quality under some timeout limit. We define good quality results withsome successful test condition(s). In this way, we aim to know whetherthe chosen MOO method works well for the proposed MOO model. Inthis thesis, the successful test conditions are defined using the hypervolumeand/or spacing indicators.Using the CPU time, we implement the so-called performance profilesdeveloped by Dolan and More´ [DM02] to identify which method performsbetter and finds good quality results. The performance profile considers the455.2. Experiment descriptionnumber of problems resolved by the method for the given successful testcondition(s), as well as the CPU time required to solve them.5.2 Experiment descriptionThe experiment is performed on 30 test problems, that is, on 30 roadsamples with distinct ground profiles. The road samples vary in length (num-ber of stations). The data are provided by our industrial partner SoftreeTechnical Systems Inc. For each problem, we test the following 10 methods.1. ε-constraint method;2. ε-constraint method with warm start;3. Weighted sum method;4. Weighted sum method with warm start;5. Weighted metric method (p = 1);6. Weighted metric method (p = 1) with warm start;7. Weighted metric method (p =∞);8. Weighted metric method (p =∞) with warm start;9. NSGA-II;10. FP-NSGA-II.Each method is run until either the problem is solved, or the allocatedtimeout limit is hit. The timeout limit is presented in Table 5.1.The visualizations of the found Pareto fronts for all 30 test problemsappear in Appendix B.Table 5.1: Timeout limits that are set for the experiment.Stations Roads Time (sec.)10-20 stations 8 30020-30 stations 7 30030-50 stations 4 180050-80 stations 4 360080-100 stations 7 10800465.2. Experiment descriptionNote that for the proposed MOO problem we skip the weighted metricmethod with p = 2 since the reformulated problem given in (4.3) with theEuclidean metric becomes Mixed-Integer Polynomial Program (MIPP) ofdegree higher than 2 that is also nonsmooth and nonconvex. For our MOOmodel, this MIPP problem was not solved by the tested solvers (SCIP,BARON, and BMIBNB) in the given timeout limits.Before conducting the main experiment on the 30 test problems, we im-plemented all the chosen scalarization approaches with/without the warmstart strategy as described in the pseudocodes provided in Chapter 4. Forthe NSGA-II method, we used the free-accessed code written in MATLABby [Abr19], while for the FP-NSGA-II method, we have accordingly modi-fied the code for NSGA-II using the pseudocode given in [FTG15]. All theselected MOO algorithms were unit tested on the 14 benchmark test prob-lems with two objectives (ZDT1-ZDT6, SRN, TNK, and other test problemsgiven in [iJ15]).Even though all the selected scalarization methods are deterministic, werun each of them 3 times on all 30 test problems to record the minimum CPUtime spent. In this way, we aim to determine the best CPU time recordedfor the least number of external operations performed in the background bythe computer.The methods such as NSGA-II and FP-NSGA-II are stochastic, andtherefore, they tend to find a different set of solutions in each run. With thisin mind, we run these methods 5 times each and record all the results found.The parameter settings used for NSGA-II and FP-NSGA-II are similar tothose used by [DPAM02, FTG15] and are summarized in Tables 5.2-5.3.Table 5.2: Parameter configurations for NSGA-II and FP-NSGA-II similarto those used by [DPAM02, FTG15]. Any other parameters were left attheir default values.Setting ValueSelection Binary tournament selectionCrossover method Simulated binary crossoverCrossover rate 0.9Distribution index for crossover 20Mutation method Polynomial mutationMutation rate 1/Number of decision variablesDistribution index for mutation 20Constraint handling strategy Constrained tournament method475.3. Experimental resultsTable 5.3: Population size used for NSGA-II and FP-NSGA-II.Stations Population size10-20 stations 10020-30 stations 10030-50 stations 20050-80 stations 50080-100 stations 1000We perform the experiment using the vehicle characteristics for theCAT 793D off-highway truck (as in Subsection 3.2.3) that can be retrievedfrom the official website of Caterpillar manufacturer. Note that the model isapplicable for any heavy truck type. Moreover, we assume that the vehiclerunning speed is constant and is 10 m/s.In this thesis, the optimization process is carried out using the well-known CPLEX solver retrieved from the academic edition of the IBM ILOGCPLEX Optimizer 12.10.0 [Cpl09]. This solver is accessed through a freetoolbox YALMIP [Lo¨f04] used for modeling and optimization in the pro-gramming platform MATLAB. For this thesis, we use the R2019b versionof the MATLAB [MAT19]. The technical characteristics of the computeron which the experiment was performed are a Dell workstation with an In-tel(R) Core(TM) i7 2.8GHz processor, a 16 GB of Random Access Memory(RAM), and a 64-bit Windows 10 operating system.5.3 Experimental resultsWe first focus on the effect of the warm start strategy on the convergencespeed of the selected scalarization techniques. Then, we aim to determinewhich is the most promising MOO method for solving the proposed opti-mization problem in the vertical alignment design.5.3.1 Effect of warm start strategy on speedAs mentioned in Section 4.6, the main idea of implementing the warmstart strategy on the scalarization techniques is to speed up the optimizationprocess by converging to the optimal solution faster.Table 5.4 presents the approximate percentage value by which the warmstart strategy accelerates the convergence speed of the chosen scalarizationtechniques. In the brackets, we state the standard deviations of the results.485.3. Experimental resultsTable 5.4: Average speed improvement (stdev) due to the warm start strat-egy.Stations ε-constraintmethodWeightedsummethodWeightedmetricmethod(p = 1)Weightedmetricmethod(p =∞)10-20 stations 15.1% (2.0) 6.8% (2.2) 15.2% (1.6) 6.8% (0.1)20-30 stations 17.3% (1.5) 18.9% (0.7) 21.7% (0.8) 6.2% (0.5)30-50 stations 18.0% (0.4) 14.0% (2.8) 11.3% (3.6) *50-80 stations 13.9% (0.1) 7.6% (0.2) * *80-100 stations 13.9% (0.5) 7.6% (0.3) 8.2% (0.8) *Average 15.7% (0.9) 11.0% (1.2) 14.1% (1.7) 6.5% (0.3)* Terminated due to the timeout limit.Moreover, we ignore the test problems that terminated with the timeoutlimit. If every problem belonging to the same entry terminated due to thetimeout limit, then we mark such entries with a star symbol. As a result,we can see that the warm start strategy indeed improves the CPU time ofthe selected MOO methods by approximately 5-15%.In the next subsection, we present results for the other performanceindicators and discuss how the warm start strategy has affected the qualityof the results.5.3.2 Examination of performance indicatorsTo present the results obtained for the hypervolume and spacing indica-tors for all 30 test problems and all 10 MOO methods, we use the box-and-whisker plots given in Figures 5.4-5.5.Each box-and-whisker plot shows the distribution of numerical results foreach MOO method by visualizing the minimum, lower quartile, median, up-per quartile, and maximum of the selected indicator. The box-and-whiskerplots label the mean values with a cross symbol and the outliers with thedots outside the box. The ends of the box are the upper and lower quar-tiles, whereas a line inside the box is the median. The whiskers are the twovertical lines that extend outside the box to the maximum and minimumobservations.For Figure 5.4, the higher hypervolume values are preferable, and theideal result would be 1. For Figure 5.5, the lower spacing values are desired,and the ideal result would be 0.495.3. Experimental resultsHypervolume0.00.10.20.30.40.50.60.70.80.91.0Spread 0.00.10.20.30.40.50.60.70.80.91.0Hypervolume (higher is better)Eps-constraint method Eps-constraint method with warm startWeighted sum methodWeighted sum method with warm startWeighted metric method (p=1)Weighted metric method (p=1) with warm startWeighted metric method (p=∞)Weighted metric method (p=∞) with warm startNSGA-II (5 runs)FP-NSGA-II (5 runs)Spread 0.00.10.20.30.40.50.60.70.80.91.0Hypervolume (higher is better)Eps-constraint method Eps-constraint method with warm startWeighted sum methodWeighted sum method with warm startWeighted metric method (p=1)Weighted metric method (p=1) with warm startWeighted metric method (p=∞)Weighted metric method (p=∞) with warm startNSGA-II (5 runs)FP-NSGA-II (5 runs)Spread 0.00.10.20.30.40.50.60.70.80.91.0Hypervolume (higher is better)Eps-constraint method Eps-constraint method with warm startWeighted sum methodWeighted sum method with warm startWeighted metric method (p=1)Weighted metric method (p=1) with warm startWeighted metric method (p=∞)Weighted metric method (p=∞) with warm startNSGA-II (5 runs)FP-NSGA-II (5 runs)!-constraint method!-constraint method with warm startighted sum ethodighted sum ethod with warm startighted metric method (# =1)eighted metric method (# =1) with warm starteighted metric method (# =∞)eighted metric method (# =∞) with warm startI ( r ns)-II (5 runs)Figure 5.4: Results for hypervolume indicator (higher is better). In a box-and-whisker plot: the ends of the box are the upper and lower quartiles; aline inside the box is the median; a cross is the mean; the whiskers are thetwo lines that extend outside the box to the highest and lowest observations.First of all, we focus on scalarization techniques. As we can see from bothplots, the scalarization techniques with the warm start strategy perform aswell as (or slightly better in some cases) the versions without the warm start.Combining this information with the timing results in Subsection 5.3.1, itis clear that using the warm start is a good choice.The cases where the warm start strategy slightly improves the perfor-mance of some of the MOO methods can be explained as follows. Among 30test problems, the scalarization techniques solved the bulk of them within atime range that is less than the timeout limit. Nonetheless, there were stillseveral test problems that required the entire timeout limit to terminate. So,the warm start strategy showed the improvement only in those instances byaccelerating the optimization process and performing more iterations to findmore solutions before reaching the given timeout limit.Now, we focus our attention on the numerical results found by theNSGA-II and FP-NSGA-II methods given in Figures 5.4-5.5. For each ofthem, we illustrate the results for all 5 runs. In all runs, both methods obtaina wide spread of results for the hypervolume and spacing indicators. This505.3. Experimental resultsSpacing0.00.10.20.30.40.50.60.70.80.91.0Spread 0.00.10.20.30.40.50.60.70.80.91.0Hypervolume (higher is better)Eps-constraint method Eps-constraint method with warm startWeighted sum methodWeighted sum method with warm startWeighted metric method (p=1)Weighted metric method (p=1) with warm startWeighted metric method (p=∞)Weighted metric method (p=∞) with warm startNSGA-II (5 runs)FP-NSGA-II (5 runs)Spread 0.00.10.20.30.40.50.60.70.80.91.0Hypervolume (higher is better)Eps-constraint method Eps-constraint method with warm startWeighted sum methodWeighted sum method with warm startWeighted metric method (p=1)Weighted metric method (p=1) with warm startWeighted metric method (p=∞)Weighted metric method (p=∞) with warm startNSGA-II (5 runs)FP-NSGA-II (5 runs)Spread 0.00.10.20.30.40.50.60.70.80.91.0Hypervolume (higher is better)Eps-constraint method Eps-constraint method with warm startWeighted sum methodWeighted sum method with warm startWeighted metric method (p=1)Weighted metric method (p=1) with warm startWeighted metric method (p=∞)Weighted metric method (p=∞) with warm startNSGA-II (5 runs)FP-NSGA-II (5 runs)!-constraint method!-constraint method with warm startighted sum ethodighted sum ethod with warm startighted metric method (# =1)ighted metric method (# =1) with warm starteighted metric method (# =∞)eighted metric method (# =∞) with warm startI ( r ns)-II (5 runs)Figure 5.5: Results for spacing indicator (lower is better). In a box-and-whisker plot: the ends of the box are the upper and lower quartiles; a lineinside the box is the median; a cross is the mean; the whiskers are the twolines that extend outside the box to the highest and lowest observations.implies that these selected GA methods are unstable for the proposed MOOproblem. A possible reason for their non-stable performance can be thevast amount of equality constraints required to define the objective func-tion for the road construction cost. These constraints are retrieved fromthe blackbox model described in [MR12, BHLH17]. Generally, the equalityconstraints reduce the feasible search space which, in turn, complicates theselection of feasible solutions for the chosen GA methods.Overall, based on the results shown in Figures 5.4-5.5, the ε-constraintmethod with the warm start is a candidate for the promising MOO methodfor our proposed problem followed by the weighted metric method (p =∞)with the warm start.5.3.3 Comparison of timeNow, we analyze the CPU time performance of the chosen MOO meth-ods by visualizing the performance profiles for different successful test con-dition(s).Before discussing the results, we briefly present an interpretation of the515.3. Experimental resultsperformance profiles. On its y-axis, we have the percentage of test problemsthat passed the specified successful test condition(s), whereas, on its x-axis,we have the performance ratio showing the ratio of time cost. For example,an (x, y) point of (5, 0.3) for an examined solver means that given 5 timesas long as the fastest solver, the examined solver successfully solved 30% ofproblems. As such, the method with the highest speed of convergence wouldbe found at τ = 0 with the highest value of the y-axis, while the methodsolving the most number of test problems with good quality would be foundclose to the top part of the plot. So, to determine the fastest and efficientmethod, we need to look for the graph closest to the top-left region of theperformance profile.Moving on to the discussion of the results, we make the following impor-tant notes. For the scalarization techniques, we focus only on their warmstart versions, because, as discussed before, they obtain as good (or slightlybetter) results as their original versions within a shorter amount of time.For the GA methods, we consider only the best run out of 5 runs for eachtest problem.Based on the box-and-whisker plots presented in Figures 5.4-5.5, wefirst choose to examine the performance of the methods for some non-strictsuccessful test conditions: hypervolume ≥ 0.75 and spacing ≤ 0.25 (bothof them implied at the same time). In this way, we aim to identify themost suitable and promising method for our MOO problem that determinesaverage quality results for most of the test problems.Figure 5.6 shows the performance profile for the selected successful testconditions. As can be seen from Figure 5.6, the ε-constraint method withthe warm start successfully solves the most number of test problems (23problems out of 30) with relatively good time performance. Nonetheless,the fastest method is the weighted sum method with the warm start, as itsolves half of the test problems the fastest. Moreover, we did an additionaltest by removing the ‘best’ method (that is, ε-constraint method with thewarm start) from the performance profile to verify that the weighted summethod is the second ‘best’ method. This was done to check that there isno switching phenomenon (see [GS16]). Lastly, for Figure 5.6, note thateven the best runs of the GA methods show poor results for the proposedoptimization model since they use the entire timeout limit to stop.In Figure 5.7, we examine the CPU time performance for other differentsuccessful test conditions.For the hypervolume ≥ 0.75 (Figure 5.7(a)), we obtain almost similarresults as in Figure 5.6, while for the spacing ≤ 0.25 (Figure 5.7(b)), wereveal the similar pattern as in Figure 5.6 but with more number of the525.3. Experimental results0 5 10 15 20 2500.10.20.30.40.50.60.70.80.910 5 10 15 20 2500.10.20.30.40.50.60.70.80.91Figure 5.6: Performance profile with the successful test conditions of hyper-volume ≥ 0.75 and spacing ≤ 0.25.successfully solved test problems. This implies that the compared MOOmethods can reach diversity easier than convergence. For the hypervolume≥ 0.90 (Figure 5.7(c)), only a small portion of the test problems passed thisstrict condition, and most of the methods could not reach the high levelof convergence. For the spacing ≤ 0.1 (Figure 5.7(d)), we reveal that theε-constraint method with the warm start performs best by solving almosthalf of the problems with a high level of diversity. Moreover, we can noticethat the GA methods obtain better results for diversity than the weightedsum method with the warm start and the weighed metric method (p = 1)with the warm start. From these two plots, we can again conclude that forthe given MOO model the diversity and uniformity of the solutions can beobtained easier than the convergence.Overall, Figure 5.7 supports the conclusions of Figure 5.6 that the ε-constraint method with the warm start provides the most robust conver-gence.535.3. Experimental results0 5 10 15 20 2500.10.20.30.40.50.60.70.80.91(a) Hypervolume ≥ 0.75.0 5 10 15 20 2500.10.20.30.40.50.60.70.80.91(b) Spacing ≤ 0.25.0 2 4 6 8 10 1200.10.20.30.40.50.60.70.80.91(c) Hypervolume ≥ 0.90.0 5 10 15 2000.10.20.30.40.50.60.70.80.91(d) Spacing ≤ 0.10.0 5 10 15 20 2500.10.20.30.40.50.60.70.80.91Figure 5.7: Performance profiles with different successful test conditions.545.3. Experimental results5.3.4 Summary of resultsWe tested 10 MOO methods on 30 road samples and analyzed the ob-tained results using the performance indicators. We find the following re-sults.− The warm start strategy accelerates the performance of the scalariza-tion techniques by approximately 5% to 15%.− The warm start versions of the scalarization methods obtain similar(or slightly better) results as their original versions.− The most promising and recommended MOO method for the proposedoptimization model is the ε-constraint method with the warm startstrategy.− The second best MOO method is the weighted sum method with thewarm start strategy. It shows the best performance in time but has anotably lower success rate than the ε-constraint method.− The GA methods perform unstably and poorly for our MOO problem.− The diversity and uniformity of the solutions can be reached easierthan the convergence.55Chapter 6Conclusion6.1 ConclusionThe cost of road design is significant as it needs large capital investmentsfor building new roads and incurs future user costs as well. This thesisproposes a multi-objective approach that seeks to find road profiles that areoptimal for the manufacturers in terms of the road construction cost and atthe same time for the users in terms of the vehicle operating costs.In our model, we use the road construction cost defined by the MILPmodel given in [MR12, BHLH17]. As the vehicle operating cost, we employthe fuel consumption cost proposed by [ASB12]. It estimates fuel usagebased on the vehicle parameters and environmental factors that affect vehiclemotion. Moreover, the fuel consumption model was validated on the datasetgiven in [SAKK16] with the coefficient of determination R2 = 0.93.In this thesis, we examine and test several well-known MOO methodson 30 road samples to determine the most promising approach for optimiz-ing our model. We focus on three classical scalarization techniques (theε-constraint method, weighted sum method, and weighted metric methods)and two commonly-used GA methods (NSGA-II and FP-NSGA-II). More-over, we propose the warm start strategy to accelerate the performance ofthe chosen scalarization techniques. The results were assessed with hyper-volume indicator (to test the convergence of solutions), spacing indicator (totest the diversity of solutions), and CPU time (to test the speed).As a result, we determine that the warm start versions of the scalariza-tion techniques obtain similar (or slightly better) results as their originalversions, as well as in less time. Moreover, the most robust and recom-mended MOO method for our proposed problem is the ε-constraint methodwith the warm start strategy. We also identify that the compared GA meth-ods perform poorly and unstably for our problem.566.2. Future work6.2 Future workBelow we outline several suggestions for future work as an improve-ment/extension of this research.6.2.1 Model extension− Our model considers only bi-objective optimization problem. Onecould include more objectives such as tire wear, travel time, etc.− Our model assumes that the vehicle running speed is constant anduser-defined. One could advance the model by introducing a new de-cision variable that can control the vehicle running speed.− Our model focuses only on the vertical alignment design. One couldextend the existing model to the horizontal alignment design and couldoptimize the horizontal and vertical alignments at the same time.6.2.2 Methods extension− Our dynamic ε/weight selection approach for the scalarization tech-niques uses Euclidean distance to calculate the range of search regions.One could improve this approach by implementing different distancemetrics so that it would select the search regions more efficiently.− Our tested GA methods suffer from the vast amount of equality con-straints that significantly reduce the feasible search space. One couldadvance them by introducing an efficient way of handling the equalityconstraints or could implement some other GA methods.− Our experiment tests only several commonly-used MOO methods. Onecould choose to employ other MOO methods such as the derivative-freeMOO methods, heterogeneous MOO methods, or some other MOOmethods that can address nonsmooth/nonconvex/nonlinear problems.57Bibliography[AB03] R. Akc¸elik and M. Besley. Operating cost, fuel consumption,and emission models in aaSIDRA and aaMOTION. In 25th Con-ference of Australian Institutes of Transport Research (CAITR2003), pages 1–15. University of South Australia Adelaide, Aus-tralia, 2003. → pages 16, 65[ABC+20] C. Audet, J. Bigeon, D. Cartier, S. L. Digabel, and L. Salomon.Performance indicators in multiobjective optimization. Euro-pean Journal of Operational Research, 2020. → pages 7, 41, 42,43, 44[Abr13] S. Abrams. The unseen history of our roads. Road &Track. https://www.roadandtrack.com/car-culture/a4447/the-road-ahead-road-evolution.html, 2013. Accessed:2021-05-15. → pages 1[Abr19] A. F. Abril. Real coded (integer handling) NSGA-II. MATLAB Central File Exchange. https://www.mathworks.com/matlabcentral/fileexchange/71609-real-coded-integer-handling-nsga-ii.html, 2019.Accessed: 2021-01-05. → pages 47[AH17] C. Audet and W. Hare. Derivative-free and blackbox optimiza-tion. Springer, 2017. → pages 32, 33, 34[ASB12] R. Akc¸elik, R. Smit, and M. Besley. Calibrating fuel consump-tion and emission models for modern vehicles. In IPENZ Trans-portation Group Conference, Rotorua, New Zealand, pages 1–13,2012. → pages 13, 14, 15, 16, 56[Asg98] M. J. Asgharpour. Multiple criteria decision making, volume 1.Tehran University Publishing, Tehran, 1998. → pages 22[BHLH17] V. Beiranvand, W. Hare, Y. Lucet, and S. Hossain. Multi-haulquasi network flow model for vertical alignment optimization.58BibliographyEngineering Optimization, 49(10):1777–1795, 2017. → pages 2,5, 11, 12, 51, 56[CC61] A. Charnes and W. W. Cooper. Management models and indus-trial applications of linear programming. John Wiley & Sons,LTD, New York, 1961. → pages 21[CC05] C. A. C. Coello and N. C. Corte´s. Solving multiobjective opti-mization problems using an artificial immune system. GeneticProgramming and Evolvable Machines, 6(2):163–190, 2005. →pages 41[CCF55] A. Charnes, W. W. Cooper, and R. O. Ferguson. Optimal esti-mation of executive compensation by linear programming. Man-agement Science, 1(2):138–151, 1955. → pages 21[CH83] V. Chankong and Y. Y. Haimes. Multiobjective decision makingtheory and methodology. Elsevier Science, New York, 1983. →pages 22[CHE13] A. Carlson, U. Hammarstro¨m, and O. Eriksson. Models andmethods for the estimation of fuel consumption due to infras-tructure parameters. MIRAVEC Deliverable D2.1, 2013. →pages 13[CPC21] CPCS. Cost of poor roads in Canada - final report for Cana-dian Automobile Association (CAA). In Solutions for GrowingEconomies, pages 1–16, 2021. → pages 1, 2[Cpl09] IBM ILOG Cplex. V12. 1: User’s manual for CPLEX. Inter-national Business Machines Corporation, 46(53):157, 2009. →pages 48[CSQ12] S. Cheng, Y. Shi, and Q. Qin. On the performance metricsof multiobjective optimization. In International Conference inSwarm Intelligence, pages 504–512. Springer, 2012. → pages 40,41, 42[Deb01] K. Deb. Multi-objective optimization using evolutionary algo-rithms. John Wiley & Sons, LTD, New York, 2001. → pages2259Bibliography[DJ13] K. Deb and H. Jain. An evolutionary many-objective optimiza-tion algorithm using reference-point-based nondominated sort-ing approach, part I: solving problems with box constraints.IEEE Transactions on Evolutionary Computation, 18(4):577–601, 2013. → pages 23[DM02] E. D. Dolan and J. J. More´. Benchmarking optimization soft-ware with performance profiles. Mathematical Programming,91(2):201–213, 2002. → pages 45[DPAM02] K. Deb, A. Pratap, S. Agarwal, and T. Meyarivan. A fast andelitist multiobjective genetic algorithm: NSGA-II. IEEE Trans-actions on Evolutionary Computation, 6(2):182–197, 2002. →pages 9, 23, 34, 35, 41, 47[EML01] Z. Eckart, L. Marco, and T. Lothar. Improving the strengthpareto evolutionary algorithm for multiobjective optimization.EUROGEN, Evolutionary Methods for Design, Optimizationand Control with Applications to Industrial Problems, pages 1–21, 2001. → pages 22[Fle03] M. Fleischer. The measure of pareto optima applications tomulti-objective metaheuristics. In International Conferenceon Evolutionary Multi-Criterion Optimization, pages 519–533.Springer, 2003. → pages 44[FTG15] S. Fettaka, J. Thibault, and Y. Gupta. A new algorithm us-ing front prediction and NSGA-II for solving two and three-objective optimization problems. Optimization and Engineering,16(4):713–736, 2015. → pages 23, 34, 35, 36, 37, 42, 43, 47[Fwa89] T. F. Fwa. Highway vertical alignment analysis by dynamicprogramming. Transportation Research Record, 1239(1-9):2–3,1989. → pages 11[GS55] S. Gass and T. Saaty. The computational algorithm for the para-metric objective function. Naval Research Logistics Quarterly,2(1-2):39–45, 1955. → pages 21[GS16] N. Gould and J. Scott. A note on performance profiles for bench-marking software. ACM Transactions on Mathematical Software(TOMS), 43(2):1–5, 2016. → pages 5260Bibliography[iJ15] Metaheuristic Algorithms in Java. Problems included in jMetal.http://www.jmetal.sourceforge.net/problems.html, 2015.Accessed: 2021-02-15. → pages 47[IMTN15] H. Ishibuchi, H. Masuda, Y. Tanigaki, and Y. Nojima. Modifieddistance calculation in generational distance and inverted gen-erational distance. In International Conference on EvolutionaryMulti-Criterion Optimization, pages 110–125. Springer, 2015. →pages 41[JMC09] A. L. Jaimes, S. Z. Martınez, and C. A. C. Coello. An introduc-tion to multiobjective optimization techniques. Optimization inPolymer Processing, pages 29–57, 2009. → pages 6, 7, 21, 24, 31[JOZF14] S. Jiang, Y. S. Ong, J. Zhang, and L. Feng. Consistencies andcontradictions of performance metrics in multiobjective opti-mization. IEEE Transactions on Cybernetics, 44(12):2391–2404,2014. → pages 44[JSJ06] M. K. Jha, P. M. Schonfeld, and J. C. Jong. Intelligent roaddesign, volume 19. WIT press, 2006. → pages 12, 13[Kru11] J. W. Kruisselbrink. Hypervolume computa-tion. MATLAB Central File Exchange. https://www.mathworks.com/matlabcentral/fileexchange/30785-hypervolume-computation.html, 2011. Accessed:2021-02-14. → pages 44[Kub15] K. A. C. Kubler. Optimisation of off-highway truck fuel con-sumption through mine haul road design. Master’s thesis, Uni-versity of Southern Queensland, Queensland, Australia, 2015.→ pages 13[Lo¨f04] J. Lo¨fberg. YALMIP: A toolbox for modeling and optimizationin MATLAB. In In Proceedings of the CACSD Conference, pages284–289, Taipei, Taiwan, 2004. → pages 48[MAT19] MATLAB. 9.7.0.1190202 (R2019b). The MathWorks Inc., Nat-ick, Massachusetts, 2019. → pages 48[MKO14] M. Ma¨kela¨, N. Karmitsa, and W. Outi. Multiobjective proximalbundle method for nonsmooth optimization. Technical report,TUCS, 2014. → pages 2361Bibliography[MM20] M. Ma¨kela¨ and O. Montonen. New multiobjective proximal bun-dle method with scaled improvement function. In NumericalNonsmooth Optimization, pages 461–479. Springer, 2020. →pages 23[MR12] F. Md Rahman. Optimizing the vertical alignment under earth-work block removal constraints in road construction. Master’sthesis, University of British Columbia (Okanagan), Kelowna BC,Canada, 2012. → pages 2, 4, 5, 11, 12, 51, 56[PF03] R. C. Purshouse and P. J. Fleming. Evolutionary many-objectiveoptimisation: an exploratory analysis. In The 2003 Congresson Evolutionary Computation, 2003 (CEC’03), volume 3, pages2066–2073. IEEE, 2003. → pages 35[RK14] J. H. Ryu and S. Kim. A derivative-free trust-region methodfor biobjective optimization. SIAM Journal on Optimization,24(1):334–362, 2014. → pages 23[RLD+01] H. Rakha, I. Lucic, S. H. Demarchi, J. R. Setti, and M. V. Aerde.Vehicle dynamics model for predicting maximum truck accelera-tion levels. Journal of Transportation Engineering, 127(5):418–425, 2001. → pages 17[SAKK16] A. Soofastaei, S. M. Aminossadati, M. S. Kizil, and P. Knights.The influence of rolling resistance on haul truck fuel consump-tion in surface mines. Tribology International Journal, 2(1):215–228, 2016. → pages 13, 14, 20, 56[Sch85] J. D. Schaffer. Multiple objective optimization with vector eval-uated genetic algorithms. In Proceedings of the First Interna-tional Conference on Genetic Algorithms and Their Applica-tions, pages 93–100. Lawrence Erlbaum Associates. Inc., Pub-lishers, 1985. → pages 22[Sch95] J. R. Schott. Fault tolerant design using single and multicrite-ria genetic algorithm optimization. PhD thesis, MassachusettsInstitute of Technology, Cambridge, Massachusetts, 1995. →pages 41[SELC12] O. Schutze, X. Esquivel, A. Lara, and C. A. C. Coello. Usingthe averaged hausdorff distance as a performance measure in62Bibliographyevolutionary multiobjective optimization. IEEE Transactionson Evolutionary Computation, 16(4):504–522, 2012. → pages 41[SF16] G. Svenson and D. Fjeld. The impact of road geometry andsurface roughness on fuel consumption of logging trucks. Scan-dinavian Journal of Forest Research, 31(5):526–536, 2016. →pages 14[TE19] J. Thomann and G. Eichfelder. A trust-region algorithm forheterogeneous multiobjective optimization. SIAM Journal onOptimization, 29(2):1017–1047, 2019. → pages 23[TV03] R. J. Thompson and A. T. Visser. Mine haul road maintenancemanagement systems. Journal of the Southern African Instituteof Mining and Metallurgy, 103(5):303–312, 2003. → pages 12[TWFT95] M. Tanaka, H. Watanabe, Y. Furukawa, and T. Tanino. Ga-based decision support system for multicriteria optimization. In1995 IEEE International Conference on Systems, Man and Cy-bernetics. Intelligent Systems for the 21st Century, volume 2,pages 1556–1561. IEEE, 1995. → pages 26[VOS14] K. D. V. Villacorta, P. R. Oliveira, and A. Soubeyran. A trust-region method for unconstrained multiobjective problems withapplications in satisficing processes. Journal of OptimizationTheory and Applications, 160(3):865–889, 2014. → pages 23[VV99] D. A. Van Veldhuizen. Multiobjective evolutionary algorithms:classifications, analyses, and new innovations. PhD thesis,School of Engineering of the Air Force Institute of Technology,Dayton, Ohio, 1999. → pages 41, 42[Zit99] E. Zitzler. Evolutionary algorithms for multiobjective optimiza-tion: methods and applications. PhD thesis, Swiss Federal Insti-tute of Technology Zurich, Zurich, Switzerland, 1999. → pages41[ZJZ+06] A. Zhou, Y. Jin, Q. Zhang, B. Sendhoff, and E. Tsang. Com-bining model-based and genetics-based offspring generation formulti-objective optimization using a convergence criterion. In2006 IEEE International Conference on Evolutionary Compu-tation, pages 892–899. IEEE, 2006. → pages 41, 4563Bibliography[ZTL+03] E. Zitzler, L. Thiele, M. Laumanns, C. M. Fonseca, and V. G.Da Fonseca. Performance assessment of multiobjective optimiz-ers: an analysis and review. IEEE Transactions on EvolutionaryComputation, 7(2):117–132, 2003. → pages 40, 4264AppendicesAppendix ATablesTable A.1: Parameters α, β1 and β2 retrieved from [AB03].Parameter Symbol Unit LightvehicleHeavyvehicleIdle fuel rate α mL/s 0.375 0.556Energy efficiency parameter β1 mL/kJ 0.009 0.008Energy-acceleration efficiencyparameterβ2 mL/(kJ.m/s2) 0.003 0.002Table A.2: Idle fuel rates for various vehicle types. Source: Argonne Na-tional Laboratory, www.anl.gov, 2014.Vehicle type Class Fuel typeIdle fuel rate (gal/hr)No load With loadMedium heavy truck 6 Gasoline 0.84 -Delivery truck 5 Diesel 0.84 1.10Tow truck 6 Diesel 0.59 1.14Medium heavy truck 6-7 Diesel 0.44 -Transit bus 7 Diesel 0.97 -Combination truck 7 Diesel 0.49 -Bucket truck 8 Diesel 0.90 1.50Tractor-semitrailer 8 Diesel 0.64 1.1565Appendix BFiguresIn this appendix, we illustrate the Pareto fronts for all 30 test problems.In each problem, we plot only the non-dominated solutions found by thescalarization methods with the warm start and the best run of GA methods.(a) Test problem 1. (b) Test problem 2. (c) Test problem 3.(d) Test problem 4. (e) Test problem 5. (f) Test problem 6.Figure B.1: Resulting Pareto fronts for 30 test problems.66Appendix B. Figures(g) Test problem 7. (h) Test problem 8. (i) Test problem 9.(j) Test problem 10. (k) Test problem 11. (l) Test problem 12.Figure B.2: Resulting Pareto fronts for 30 test problems (cont.).67Appendix B. Figures(m) Test problem 13. (n) Test problem 14. (o) Test problem 15.(p) Test problem 16. (q) Test problem 17. (r) Test problem 18.Figure B.3: Resulting Pareto fronts for 30 test problems (cont.).68Appendix B. Figures(s) Test problem 19. (t) Test problem 20. (u) Test problem 21.(v) Test problem 22. (w) Test problem 23. (x) Test problem 24.Figure B.4: Resulting Pareto fronts for 30 test problems (cont.).69Appendix B. Figures(y) Test problem 25. (z) Test problem 26. (aa) Test problem 27.(ab) Test problem 28. (ac) Test problem 29. (ad) Test problem 30.Figure B.5: Resulting Pareto fronts for 30 test problems (cont.).70
Thesis/Dissertation
eng
Computer Science
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Hare, Warren
Multi-objective optimization for road vertical alignment design
Text