Slope Optimization (or “Sloop”): Customized Optimization for Road Longitudinal Proﬁle Eco-Design

: Current transportation systems contributed to one quarter to the Greenhouse Gases (GhG) emissions and the mobility demand increases continuously. The transportation infrastructure design has to be optimized to mitigate these emissions. The methodology “Sloop” (acronym for Slope Optimization) has recently been set up to optimize the longitudinal road proﬁle with respect to a Global Warming Potential (GWP) criterion calculated for both the construction and operational phases while incorporating accurate vehicle models. This paper proposes a customized optimization strategy that signiﬁcantly improves the Sloop methodology. From an initial longitudinal proﬁle generated by road designers, our algorithm identiﬁes the optimized proﬁle in terms of GWP. This optimization step is complex due to the large number of degrees of freedom, along with various constraints to obtaining a feasible solution and the computational cost of the operational phase assessment. The most stringent constraint entails connecting the proﬁle with the existing road (i.e., a constraint on the ﬁnal proﬁle altitude). We have developed a speciﬁc algorithm based on the Sequential Quadratic Programming (SQP) method; this algorithm takes into account the quasi-quadratic nature of the constraint on the ﬁnal altitude. In a real case study, our algorithm outputs a proﬁle whose GWP assessment saves 4% compared to the GWP assessment of the initial proﬁle. The beneﬁt of our speciﬁc protocol for treating the ﬁnal altitude constraint is demonstrated in this example. By means of this new efﬁcient and open-box algorithm, proﬁle evolution at each iteration is analyzed to determine the most sensitive degrees of freedom, and the sensitivity of the optimization with respect to the main construction parameters and trafﬁc assumptions are conducted. The results are consistent and stable.


Introduction
According to [1], "only an acceleration in structural changes to the way the world produces and consumes energy can break the Greenhouse Gases (GhG) emissions trend for good". Knowing that emissions from transport amounted to one quarter of the GhG emissions in 2016 and that the mobility demand increases year after year [2], all new transportation infrastructure project have to be optimized to mitigate these emissions.
Within the eco-design framework, a methodology called PEAM (for Project Energy Assessment Method) was specifically developed to assess transportation projects according to the geometry of the infrastructure taking into account construction and operation phases i.e., the planned traffic impact. PEAM was then applied to either evaluate High-Speed Line projects [3], road speed-sectioning [4] or to assess the longitudinal profile i.e., slopes of road project [5].
Researchers specialized in vehicle fuel consumption have pointed out the importance of slopes traversed on fuel consumption, from the seminal investigations of Webb in 1952 [6] to recent papers of [7][8][9], yet their logic still treats slopes as constant. The construction phase has long been a subject of research, dating back to the mathematician Monge in 1781 [10]. Most research conducted on such assessments has focused on the construction and maintenance phases since these steps are monitored by road builder [11,12]. Most of the advanced research in this field of road design, e.g., [13,14], account for the influence of slopes on traffic-related fuel consumption by use of an average slope in a statistical model. This modeling set-up only provides a rough approximation of reality, depending on the domain of validity inherent in the statistical model. In a recent work [15], the authors emphasize the necessity of an accurate model of vehicles to assess the effect of slopes on fuel consumption. An optimization of the road profile is proposed based only on fuel consumption with a profile finally difficult to build.
In [5], the authors implement a local optimization called Sloop for Slope Optimization. The starting point is a profile proposed on the drawing board. The objective function is the global warming potential (GWP) of both the construction and operation phases i.e., the traffic GhG emission. The degrees of freedom in this optimization process pertain to the longitudinal road profile. The operational phase is evaluated by simulating traffic over a 10-years period. The traffic simulation is implemented by means of vehicle dynamic models associated with an engine model in a free-flow regime. The construction phase is evaluated by examining the earthworks, which are estimated by computing 3D geometric variations between the natural terrain and the longitudinal profile. The optimization was carried out with a black-box algorithm (i.e., the toolbox optimization from Matlab c ) as a proof of concept of the methodology.
Sloop faces a major limitation since the optimization step becomes highly complicated by the use of such precise models. This paper proposes an efficient open-box optimization algorithm specially tailored for this problem which paves the way to the development of an efficient tool for the practitioners.
This paper focuses on the optimization step of the methodology, in ensuring: • the optimization problem is fully described; • the choice of the algorithm is justified; • the algorithm is described and specially tailored for the problem.
By means of this new efficient and open-box algorithm, we analyze the sensitivity of • the evaluation model with respect to the parameters of the longitudinal profiles; • the optimization with respect to the main construction parameters and traffic assumptions.
The optimization problem will be completely described, in terms of degrees of freedom, constraints, and objective function in Sections 2 and 3. Section 4 constitutes the core of the paper by justifying, describing and adapting the optimization algorithm. Section 5 will then present a case study whose highlights and limitations will be discussed in Section 6. The last section concludes this paper by offering some leads for handling the limitations of our approach.

Degrees of Freedom: Parameters of the Longitudinal Profile Model
Our algorithm minimizes the GWP of a road by choosing a feasible longitudinal profile. This section displays the equations of such a profile. Equations of this type depend on the degrees of freedom parameters, which also serve as the optimization problem input. To obtain a feasible solution, the degrees of freedom are constrained. The corresponding constraints are listed in the second part of this section.

Road Longitudinal Profile
The longitudinal profile Ps is modeled as illustrated in Figure 1. It is constituted of a sequence of straight lines connected by circular arcs. The following set of road-specific equations describes this model. (1) where α k denotes the longitudinal slope of the straight line k, x k , h k are the abscissa and altitude of the point of the intersection of straight line k and straight line k + 1. b k−1 , a k are internal variables coming from the computation of the circle connecting this straight line k to the next straight line.
y ck parallel to line k straigth line k + 1 with slope ® k+1 x ck Figure 1. Longitudinal profile (in red) modeled as a sequence of straight lines connected by circular arcs (adapted from [5]). C x k ,h k ,α k ,R k ,x k+1 ,h k+1 ,α k+1 is the equation of the circle, with radius R k , that connects straight line k to straight line k + 1. This circle is tangent to straight line k at the point with abscissa a k and it is tangent to the straight line k + 1 at the point with abscissa b k . By virtue of these equations, the derivative of Ps(x) is continuous. Road designers often use circular arc to "smooth" the profile: large radii (R k > 10,000 m) tend to be chosen: circular arcs cannot be neglected in the longitudinal profile model.
The center of the circle, whose coordinates are x ck , y ck , lies at the intersection of the parallels at a distance R k of the straight lines k and k + 1. a k is the projection of the center's circle on the straight lines k. b k is the projection of the center on the straight line k + 1. The following equations specify this construction, where x ck , y ck is the solution to the following system of equations.
(2) t k is the offset of the straight line k (t k+1 for k + 1). These equations have been written under the assumption that the circle is above the straight lines i.e., the curve is locally convex as it is the case presented on Figure 1. In the concave case, R k 1 + α 2 k and R k+1 1 + α 2 k have to be subtracted . After solving this system of equations, a k (b k ) can be computed as intersection of line α k (α k+1 ) and the perpendicular to this line passing through the center of the circle. Solutions to this double system of equations are: The sign of the last equation must be adapted when the curve is locally concave. The case when one line is vertical has not been treated since it is unrealistic for a road context (the road is at least continuous i.e., without step). The boundary conditions model the connection with the existing road: Ps(x n ) = Ps e . α b is the longitudinal slope of the existing road before x b , which is the beginning of the profile to optimize, α e is the longitudinal slope of the existing road after x e , which is the end of the profile to optimize.
The degrees of freedom are: x k , α k , R k , with k = 1 to n − 1 and R 0 , R n , α n . The other variables are dependent from these degrees of freedom. h k is calculated from x k−1 , x k and α k . a k , b k are computed with the equation of circle C(x k , h k , α k , R k , x k+1 , h k+1 , α k+1 ) as mentioned above.

Constraints
The degrees of freedom cannot evolve freely. Each degree of freedom is first bounded by rules expressed in the road standards. For example, the slope cannot exceed a certain threshold for safety reasons: a loaded truck should not decelerate too strongly [16]. The degrees of freedom must also be restricted in order to model a road that can actually be built: feasible slopes depends on the type of materials due to the stability of embankment, the road project must be tied to the existing road. To specify these before the construction, the road is built virtually in a sequence from the beginning to the end of the profile i.e., from k = 0 to n. Some of these constraints are linear, which can be taken into account relatively easily, while others are nonlinear.

Linear Constraints
The intersection of both straight lines k and k + 1 has to be located after the previous intersection, as ensured by the following inequality: with k varying from 1 to n. Moreover, the last intersection must be far enough from the beginning of the actual road to link both straight lines by a circle with a compatible radius. This is expressed by the following inequality where x end is set by the optimization process user (x end has to be set before the end of the studied profile, x end < x e ).
x n < x end (6) Two consecutive straight lines have to be different in order to generate an intersection point.
d α is chosen to ensure a reasonable intersection point. Activation of this constraint serves as a signal to stakeholders that the road profile can be simplified by considering one ramp with the same slope instead of two ramps with two different slopes. Activation of an inequality means that this inequality becomes a strict equality during an iteration of the optimization algorithm.

Non-Linear Constrains
In order to build the road, it is essential that (see Figure 1) : These inequalities are nonlinear because a k and b k depend non linearly on α k and R k . The profile needs to be connected to the existing road, which implies that the profile assumes prescribed values on the boundary. According to the previous subsection, these boundary conditions are expressed as: Ps These constraints penalize heavily the optimization problem as • they are nonlinear ; • they are equality constraints, hence they must be taken into account at each iteration of the optimization process.
Since the constraints are specified as if the road were being built sequentially, the constraint Ps(x 0 ) = Ps o is automatically satisfied. According to our model, the derivative is also automatically continuous. On the other hand, the constraint Ps(x f ) = Ps f must be added to the problem. Given that the optimization process begins with a compatible solution, respecting this constraint merely entails considering that the derivative of this constraint equal zero in the direction of the descent, which means that the gradient of this constraint has to be computed. Our model (see Figure 2 to support this analysis), therefore yields: ∂Ps Except for the derivative according to the last radius R n , which is computed with Equation (4).
It turns out that the constraint is quadratic as the second derivatives are constant except for the degree of freedom R n . This " almost quadratic " property will be utilized in order to simplify the optimization problem.

Objective Function
Starting from a feasible road longitudinal profile, the objective function outputs the GWP of the road. This function is the sum of two independent assessments: • the assessment of the construction phase; • the assessment of the operation phase.

Construction
The GWP assessment of the construction phase has been fully described in [5]. Let's recall herein the main features to provide an overview.
Knowing a longitudinal profile, a 3D model of the road sub-grade, which supports the pavement layers, is calculated. Earthworks are required to transform the natural terrain into this sub-grade road. Different longitudinal profiles imply significant differences in excavation quantities. On the other hand, these longitudinal profiles do not lead to significant changes to the pavement layers. Their widths and lengths mainly depend on traffic levels. Therefore, road pavement layers have not been taken into account in the optimization process.
Earthworks include the following tasks: raw material extraction, transportation, and deposit as fill (in the case of reuse) or storage as a final deposit [17]. Assessment of Earthworks includes an estimate of the fuel consumed by earthworks equipment (bulldozer, scraper, dump truck) during each task. Soil treatment also needs to be quantified. The consumed fuel depends on the earth movements which are estimated from volume differences between the natural ground level and the road project longitudinal profile. The GhG emissions is assessed using the software ECORCE M, which features an extensive database relative to earthworks equipment [18].

Volume Computation
The longitudinal profile of the project, denoted Ps, and the natural soil profile, Pn, are altitudes relative to the metric point, i.e., distance x; they are discretized every d meters. For each point on these discretized profiles, the difference in altitude is computed and multiplied by the road infrastructure width, w ph , to yield a surface area measure. In turn, this surface area is multiplied by the segment length, d, to obtain a volume. Embankments with a given transverse slope, s b , are added along each side of this volume as displayed in Figure 3.
An elementary volume, V is computed for each point on the discretized profiles.

Assessment of Ghg Emissions from the Volume Computation
Starting from these volumes, the GWP are estimated by taking into account the raw material extraction, the transportation between excavation and embankment, the laying of materials, and the quicklime treatment of natural soil.
GWP is calculated with the following equations: gwp f (gwp cnu ) is the GWP per cubic meter for the fill (the cutting spread in a final deposit).
gwp cnu = gwp extrac + gwp trans × Dist D + gwp spread + gwp lime × r limeC (15) where gwp extrac is the GWP per cubic meter to extract the materials. gwp trans is the GWP to transport one cubic meter of materials over one kilometer. Dist FC is the average distance between cut and fill, gwp spread is the GWP to spread one cubic meter of materials, gwp lime is the GWP associated with one cubic meter of quicklime, r limeF is the quicklime rate of treated embankment materials, r limeC is the quicklime rate of treated materials stored in the final deposit and lastly Dist D is the average distance between cut and final deposit. GWP AddMat is the GWP required by additional material in the case of a negative balance between fills and cuts (i.e., for a work site with a material deficit).
The optimization process includes GWP c as part of the objective function to be minimized. For organizational purposes, a work site is in a state of either materials deficit or surplus. It becomes too complicated to switch the organization from one type of work site to another. For this reason, if the work site is in a state of surplus for the proposed profile, the following nonlinear inequality gets added to the optimization problem: The inverse inequality would be added in the case of a work site with a materials deficit.

Operations
A traffic simulation model has been used to assess the operational cost (see [5] for more details). Different classes of vehicles are modeled from small passenger cars to heavy-duty vehicles. Both diesel and gasoline-powered engines are considered for passenger cars. An average vehicle model has been developed to represent the fuel consumption of the entire class of like vehicles. Fuel consumption is obviously dependent on the slope and the speed, which is limited by the applicable speed limit or the vehicle's capabilities. The vehicle models have been developed under the customized Matlab/Simulink c library, called VEHLIB [19]. Fuel consumption is derived working backwards from wheel forces to the engine through component sub-models. The chassis model describes the aerodynamic and rolling resistance forces used to solve the fundamental dynamic principle. Wheel forces are then known and relayed to the engine shaft using a power model composed of ratios and efficiencies of the final gear and gearbox. If needed, a clutch model may also be involved (at low speed). A specific fuel consumption map is then used to assess engine fuel consumption, which then gets transposed into GWP (CO 2 emissions).

Optimization
The previous section presented the methodology employed to assess the cost, in GWP, of a road longitudinal profile. This section will demonstrate how to derive an optimal profile that minimizes GWP. The mathematical formulation is as follows: where J(Ps) is the objective function of the longitudinal profile Ps, in accordance with the previous section. C(Ps) ≥ 0 summarizes the non linear inequalities described in Section 2.2.1 and the inequality added in the construction phase (Equation (16))

Choice of Algorithm
The following discussion presents the features of this optimization problem. The choice of a suitable algorithm will then be justified.
Optimization type The problem is a finite dimensional nonlinear optimization with various linear and nonlinear inequalities and one nonlinear equality. Global/Local optimization The optimized profile cannot be too dissimilar from the initial profile for economic and mechanical reasons: earth movements are costly, embankments too large are not stable. The optimization is rather local. Initial solution An initial designed profile is available; it is considered feasible to the extent that it satisfies the constraints. Smoothness We analyze this mathematical property for the objective function and for the constraints.

•
The objective function is relatively smooth as the sum of two smooth functions, i.e.: 1. The construction assessment is based on a volume computation. This function has continuous derivatives of all orders (C ∞ ). 2. The operational assessment is the sum of the fuel consumption of various vehicles rolling virtually on the road. This sum numerically regularizes the integrand. However, the integrand is not as smooth as the construction cost. For example, a slight change in slope can lead to a drop in engine efficiency that locally alters the consumption or can lead to mechanical braking, which changes overall fuel consumption on the road. This function is at least C 1 .
• The nonlinear constraints stem from geometric considerations; they are C ∞ .

Computational burden
The construction assessment and nonlinear constraints do not require a large amount of computational operations since they rely on simple geometric analysis. On the other hand, the operational assessment requires the simulation of vehicle dynamics, which is a computation time-consuming task. Nearly quadratic equality According to Section 2.2.2, the equality constraint on altitude is nearly quadratic. Given that such a constraint must be verified at each iteration, an efficient algorithm needs to take advantage of this property.
This type of optimization is classically solved by sequential quadratic programming (SQP) [20][21][22]. SQP is a local optimization solver starting from an initial solution; however, SQP exhibits certain weaknesses as regards this specific problem. The following discussion will list these weaknesses and our solutions to mitigate them.

•
Generic SQP algorithms do not take the quadratic equality constraints into account. We have developed a dedicated solution to take advantage of this feature. • SQP imposes that the cost function and constraints can be continuously differentiated twice (C 2 ), whereas the operational assessment is at least C 1 . At each iteration, we verify that the quadratic approximation of the Lagrangian is validated by comparing this approximation with the actual value of the objective function. • SQP requires a computation of the Hessian of the Lagrangian and then of the Hessian of the objective function. The numerical computation of this Hessian is not straightforward because: the number of inputs may be high; -the computation of the operational assessment is time-consuming.

Description of the Algorithm
Our algorithm is an SQP with BFGS approximation adapted to our specific equality constraint.

Sqp Algorithm
We have followed [23] for a practical implementation of SQP, including the BFGS approximation of the Hessian. To simplify the presentation, optimization problem (17) is reformulated: SQP is an iterative algorithm starting from an initial profile: with d k is the descent increment; it is computed according to a two-step process.
1. by the descent direction p k solving the Karush Kuhn Tucker (KKT) conditions [24] of the quadratic optimization with linear equality constraints. The quadratic program is as follows: with is the Lagrangian of the non linear optimization problem. λ k is the Lagrange's multiplier at step k. ∇ 2 xx is the Hessian according to x and ∇ is the gradient operator. The corresponding KKT conditions are: The linear equalities ∇ h(x k )p k = −h(x k ) include: • the linear approximation of the non linear equality constraint on the altitude; • the linear approximation of the active inequalities, i.e., the inequalities, numbered l, which are saturated and hence become equalities, g l (x k ) = 0.
This step requires to compute the gradient and hessian of the objective function and constitutes the main computational cost of the algorithm. The gradient is computed numerically by centered finite difference; this requires (n + 1) simulations of traffic. We have applied the BFGS approximation to compute the hessian according to the equations given in [23] so as to avoid to carrying out n × (n + 1)/2 traffic simulations in order to compute the Hessian by finite difference. 2. by reducing the length of the increment d k = µp k with 0 < µ ≤ 1 being accommodated: • activation of one non active linear inequality; • activation of one non active linearized nonlinear inequality ; • an increment can respect linearized nonlinear inequalities while not respecting the nonlinear inequalities themselves. In this case, the increment's length is reduced until the inequalities are validated, with this reduction factor being due to the linear approximation of non linear constrains; • along the same line, a reduced increment can output a solution that fails to respect the linear approximation of the non linear equality on altitude. In this case, a reduced factor is computed, as will be described in the next subsection; • if the reduced increment yields a solution which increases the objective, the increment is again reduced until the objective decreases with this reduction factor being due to the quadratic approximation of the cost;

Adaptation of Sqp to the Nonlinear Equality on Altitude
In our particular study, the nonlinear equality constraint of the optimization problem (17) is a sequence of straight lines and circle arcs, which means that this constraint is in fact quadratic equality constraint as described in Section 2.2.1 except for the last parameter R n .
Numerically, the constraint on altitude is written as follows: where tol is the tolerance and x k+1 f = x k f + µ alt p k . µ alt is the reduction factor of the descent p k which yields the profile x k+1 f to be rendered consistent with the altitude constraint.
Let's consider that Ps(x) is quadratic, thus yielding: if this last inequality is not true when µ alt = 1, then µ alt would have to be reduced. Let's now suppose that the +tol has been exceeded, µ alt then provides the solution to the second-order equation: The calculation is straightforward. The equation is symmetrical if −tol has been exceeded. This equation thus directly yields the reduction factor µ alt of the increment. This solution is much more effective than when ignoring the quadratic property of this constraint, as it is the case with a generic algorithm. A generic algorithm actually guesses and approximates µ alt iteratively, in requiring at each iteration the longitudinal profile computation. Since the profile is nearly quadratic and as a result of numerical error, this factor might need to be reduced as well by another reduction factor in order to numerically respect the altitude constraint (this second reduction factor is denoted µ altnol ).

Case Study
Our case study is a road project built prior to our study. Its construction and traffic data have been fully described in [5]. We will recall herein the key points necessary to understanding the article.
The project is an 8.7-km long, 2 × 2-lane road within the French national highway network; the optimized portion of the project extends 4.75 km in length. The entire section was covered by the same asphalt pavement mix. All technical data (earthworks, pavement structure) required to assess GWP during construction work were collected from the relevant road authorities.

Construction Phase Parameters
During the construction phase, the GhG emitted to extract, spread, transport and treat are indicated, according to Equations (12) and (13) plus the parameters given in Table 1.
The site specific parameters are: average distance between cut and fill (Dist FC = 1.132 km (0.621 mi)), average distance between cut and final deposit (Dist D = 0.869 km (0.534 mi)), average quicklime rate of the material treated for fill (r limeF = 1.04%) and the quicklime rate for the final deposit, set at 0% given the absence of treatment (r limeC = 0%). These numbers are accurate because dedicated experiments were carried out in 2012 as part of the national TERDOUEST project (funded by the French National Research Agency). This site had a material surplus: inequality constraint (16) is added to the optimization problem.

Simulation Parameters of the Operations Phase
Traffic was measured on a daily basis. Table 2 displays its composition. The average daily traffic equaled approx. 7000 vehicles in both directions. We decided to consider the traffic as stable, in terms of trend in the number of vehicles over 25 years. Different traffic assumptions can be tested by using the new algorithm presented in this paper. The variation of the optimization with respect to traffic assumptions is analyzed in Section 5.5. To simulate this traffic, we divided it into vehicle classes. In each class, dynamic models were introduced to represent the vehicle behavior. In France, diesel engines represent roughly 70% of all passenger cars, and nearly 100% in trucks and vans. Moreover, the "diesel personal car" class was split into three categories, modeled by a given vehicle type (small, medium or large), as available in the VEHLIB base. These models were calibrated using data from vehicles that had been tested and characterized at Ifsttar. The heavy-duty vehicle characteristics were provided by the CEREMA-Lyon Institute. The number of vehicles was determined by calibration with the data obtained from CEREMA-Lyon during the traffic measurement campaign. Table 3 lists the technical characteristics of the various vehicles used for this assessment : • number of vehicles/day in each direction over the studied part of the road; • mass of each type of vehicle (including load and passengers); • product of the drag coefficient and the front surface area, denoted S.C x ; • tire rolling resistance coefficients, (C rr0 ) represents the constant part (C rr1 ) the part proportional to vehicle speed.  Figure 4 displays the initial and optimized profiles, along with the natural terrain. The initial profile has been designed manually by specialists in eco-design in adopting the same objective as for the algorithm, i.e.: mitigation of GWP during both the construction and operational phases. For this reason, the profile is smooth, meaning that slopes are limited. Let's begin with this profile in order to verify whether the algorithm can find a better profile than that optimized by specialists.

Optimized Profiles
This profile is denoted Ps initial (x); its length equals 4.75 km. Since the actual road is divided into six subsections the profile is also divided into six subsections (n = 6). The number of degrees of freedom is therefore 19, according to the formulation of the simplified optimization problem (see Equation (17)), i.e.: • 6 abscissas defining the intersection points between two consecutive straight lines; the displacements of abscissas are limited to +/− 100 m in order to obtain a profile not to far from the initial. • 6 slopes, one per section; Slopes can vary from −5.3% to 5.3%. These bounds are set in order to ensure stability of the embankment. • 7 radii for the circle linking the sections, since the two radii of the two extremities are also degrees of freedom; Radii vary from +/− 2000 m around the initial radii in order to keep a profile relatively smooth.

•
The tolerance on the altitude constrain is +/− 0.1 m. This step with a maximum height of 10 cm is easily flattened during construction.
According to Figure 4, the optimized profile has a simpler shape and is significantly below the initial profile. Figure 5 presents the profile derived by the algorithm at various iterations. For the sake of clarity, we have not drawn any iterations that fail to provide graphical information on the algorithm's evolution. Table 4 shows the evolution in emissions stemming from the construction and operation phases. Each iteration gives the reduction factors limiting the length of the descent increment. The final reduction factor is the product of the following reduction factors, which are to be applied sequentially (the explanation of these reduction factors is given in the Section 4.2): 1. µ lin : activation of one inactive linear inequality; 2. µ nolin : activation of one inactive linearized non linear inequality; 3. µ nol : linear approximation of the non linear equality constraints; 4. µ alt : compatibility of the quadratic approximation of the altitude constraint (µ alt is the solution of Equation (24)); 5. µ altnol : quadratic approximation of the altitude constraint; 6. µ obj : quadratic approximation of the objective function.
These reduction factors are set to 1 when inactive. For example, if multiplying the 5 first reduction factors yields a solution that decreases the objective function then µ obj = 1. For the sake of clarity, reduction factors µ lin , µ nolin , µ nol have already multiplied to display a single synthetic factor µ ineq = µ lin × µ nolin × µ nol in the table. c alt is the difference in altitude between the initial profile and the profile computed by the algorithm. We have set a tolerance of 0.1 m, which implies that −0.1 ≤ c alt ≤ 0.1.
The first row of the table provides the assessment of the initial solution. The construction cost accounts for less than 8% of the overall cost. The bottom row of the table displays the difference in the optimized solution relative to the initial solution. A GWP savings of about 4% has been gained: the GWP of the operational phase is stable, while the construction phase GWP has been halved.
According to this table, the algorithm first finds a solution that decreases the construction cost of the initial solution and then refines this solution (Iterations 2 to 4). After this first exploratory step, distinct profiles are identified whose construction phase assessment shows a significantly smaller GWP but whose operational phase emissions exert a greater impact (Iterations 5 and 6). Next, the algorithm proposes distinct profiles whose the operational phase emissions is lowered by smoothing the slopes (Iterations 7 to 9) but the emissions due to the construction increase. Iteration 10 refines the output of the last iteration to decrease the construction cost. Iterations 11 through 16 do not substantially improve the objective function, as the constraint on altitude is bounded by the upper limit of the tolerance interval (c alt = 0.1), which serves to limit the descent (µ alt × µ altnol is small). Then, the algorithm determines other profiles that decrease the construction emissions (Iterations 17 to 18) then decrease the traffic emissions (Iteration 19) before ultimately refining the solution (Iterations 20 to 25).

With Optimization Variables
One of the interests of an open-box algorithm is to be able to analyze each iteration of the algorithm. We can thus relate the variations of the objective function to the variations of the degrees of freedom. A sensitivity analysis of the objective function with respect to the degrees of freedom of the optimization is performed by using this method. Table 5 shows the iterations whose the variations in traffic and construction emissions were the largest: these are iterations 2, 5 and 7. The first row of the table provides the assessment of the initial solution. The first 4 lines present the 6 variations of the intersection points between two consecutive straight lines. The following 4 lines shows the variations for the 6 slopes. The last lines displays the variations for the 7 radii.
Major variations in GhG emissions occur at iteration 5 (variation of 851 and −2838 Mteq CO 2 for the traffic and construction, respectively). They correspond to significant variations for the slopes (variation from 0.18 to 0.54%). The variation of the other parameters are not significant (maximum variation of 2 m for the intersection points and 5 m for the radii). If this result is confirmed on various case studies, the optimization could be simplified by considering only slopes.

On Construction Parameters
The purpose of this section is to estimate the influence of the variation of the construction phase parameters on the Sloop solution. Table 1 displays the main parameters. According to Equations (12) to (15), the influence of these parameters are linear on the cost construction. If this influence were important, Sloop solution could be totally different. To test this assumption, we performed 4 optimizations, one by parameter, by increasing each parameter by 5%. Table 6 displays this sensitivity analysis. The column Parameter specifies the parameter that is modified. The first two lines correspond to the optimization with the initial values of the parameters. Profile init is the starting point of the optimization. Profile opt is the profile optimized. By subtracting the cost of initial profile of each change of parameter (101,223; 101,232; 101,226; 101,552 Mteq CO 2 ) to the cost of initial profile when the parameters are fixed to their initial values (101,210 Mteq CO 2 ), we observe that only the change of gwp lime has a significant impact on the cost function with an increase of 0.4%. This is consistent with the practitioners of Life Cycle Analysis who underline the influence of Lime on GWP [5]. A comparison of the cost of the gwp lime profile optimized (97,848 Mteq CO 2 ) with the profile optimized computed with the initial value (97,505 Mteq CO 2 ) shows a small difference (0.35%) consistent with our expectations: as the cost of construction increased relatively to the operation cost, the optimization process focuses more its effort on the construction cost (3717 Mteq CO 2 when the parameter changes and 3917 Mteq CO 2 in the initial scenario). The costs of the profile optimized when the other parameters change do not differ significantly with the profile optimized with the initial values (difference < 0.06%).

Evolution of Vehicle Fleet Scenario
The scenario proposed to develop and test our model is based on the fleet measured in 2018 in France. At that time, the part of electrical vehicle were negligible. This scenario called scenario 0 assumes that the share of electrical vehicle will remain negligible in the next 25 years. This scenario is not realistic except to demonstrate the feasibility of Sloop. With the algorithm presented in this paper, practitioners will be able to test different scenarios and to analyze the results. This subsection is an example of this approach and gives a first idea of the sensitivity of the fleet composition on the results.
Two scenarios including electrical vehicles have been added. As for the internal combustion engine vehicle, the electrical vehicles have been separated in different classes: small, medium and large classes for the passengers cars and Mini-truck and van for the duty vehicle. Their characteristics are presented by Table 7. The first scenario (scenario 1) is optimistic, the second scenario (scenario 2) is a medium way between scenario 0 and scenario 1. In these scenarios, the global number of vehicle remains the same along the next 25 years. The number of electrical vehicle increases with a linear rate from 0 (the first year) to a maximum which is different according to the scenario 1 or 2. Conversely, the number of diesel vehicle decreases following a linear rate to reach a certain part of the passengers and duty vehicle in 25 years. In scenario 1, the diesel passengers cars completely disappear in 25 years, and 50% of passengers cars will be electric. Thus 50% remain gasoline cars. For duty vehicles (Mini-truck and van), 50% become electric and replace the diesel cars. In scenario 2, a quater of passengers cars and duty vehicles are electrified in 25 years. 50% of diesel passengers cars disappear. Table 8 presents the number of vehicles per days for the first year and in 25 years for the two scenarios.
A key point, when considering electrical vehicle, is the energetic mix of the country which is highly influenced by the mode of production of electricity. We consider in this study a GWP of 75 gCO 2 per kW·h, mean value for the French electrical mix highly composed of nuclear electricity [25]. It is clear that the results are influenced by this value for others countries. Figure 6 shows the optimized profiles according to scenarios 0, 1 and 2. The curves of the 3 scenarios are close to one another, but not identical. Globally, Figure 6 shows that the profile is not highly influenced par the introduction of electrical vehicle for this particular case study. This is confirmed by Table 9 that presents the gains between initial and optimized solution according to the 3 scenarios.

Discussion
• Since vehicle consumption accounts for roughly 90% of GHG emissions, we anticipated finding an optimized profile with milder slopes than the initial profile because, in this case, consumption would be lower. However, the optimized profile proved to have steeper slopes. This finding can be explained by the fact that the ramps remain sufficiently mildly inclined so as to not significantly modify engine efficiency and not trigger braking. If consumption remains rather stable, the algorithm will minimize construction-related emissions, which are largely due to the lime treatment of embankments as underlined by the sensitivity analysis. For this reason, the optimized profile minimizes the embankment while holding vehicle consumption stable. Yet an optimization that only takes construction into account would output a profile with a steep ramp, as is the case for Iteration 5 in Table 4 in order to avoid the embankment. In this case, the traffic-related emissions increase significantly.

•
As expected, according to Table 4, the nearly quadratic constraint on the altitude often limits the descent during several iterations: µ alt lies below 1 for two thirds of the iterations and in half of the iterations, µ alt is the appropriate reduction factor of the increment to validate the altitude constraint (µ alt < 1 and µ altnol = 1), a finding that underscores the benefit of taking this property into account by the algorithm under the real-world scenario.

•
This algorithm explores various feasible profiles with distinct features, thus providing stakeholders with different possibilities. They can refine one of the proposed profiles if deemed relevant. For example, given that the optimized solution includes two successive slopes very close to one another, stakeholders can propose a simpler profile to replace these two successive ramps by just one ramp. In its practical implementation, our algorithm takes this special case into account by automatically equalizing successive slopes when they are nearly equal. This type of ad hoc procedure is not robust (i.e., a discontinuity is artificially created) and difficult to tune (e.g., how to numerically define what constitutes "nearly equal" successive slopes?).
• This last point suggests a limitation of the algorithm, which depends on the number of slopes defined in the initial profile by road designers. An alternative would be to directly optimize the profile in an infinite dimensional space through the use of variational methods. The optimized solution would then be approximated by lines and circles in order to transform this optimized, but not necessarily feasible, solution into a longitudinal road profile that can indeed be built. Research along these lines is currently underway; this approach does not guarantee a feasible solution.

•
The degrees of freedom are different in nature and are expressed in different units: angles, abscissa, radii. The implementation of the optimization algorithm does not directly use them: a scaling is first performed.

•
The optimization with just a single criterion oversimplifies the problem. One perspective would call for a multi-criteria optimization in order to propose to stakeholders a set of efficient solutions (a Pareto front). It should be noted that stakeholders already have a feasible profile at each iteration in using our algorithm.

•
Our vehicle consumption models are more sensitive to slope than the usual models used in environmental assessment. However, they do not take into account the specific deformation of the tire on uphill slopes as described in the article [26]. This phenomenon depends on the tire pressure which is not taken into account in our model. This is a point of improvement in our modeling.

•
Our case study does not include civil engineering structures such as bridges or tunnels. Our methodology is based on Project Energy Assessment Method (PEAM) and PEAM also models civil engineering structures as mentioned in [3], so there are no theoretical limits to include civil engineering structures in our methodology. Implementation can be complicated because the construction cost function may not be continuous. We plan to implement them according to users' requests.

Conclusions
Sloop is a methodology that optimizes a road's longitudinal profile according to GWP criteria. Its solution entails a customized sequential quadratic program starting from an initially designed road profile. At each iteration, a quadratic optimization problem is solved. The length of the descent increment often depends on a linear approximation of the nonlinear constraint that models the project's connection to the existing road. We have utilized the quasi-quadratic property of this constraint to approximate this length. During an application of our algorithm to a real-world example, this approximation was activated in two thirds of iterations. In comparison with a black-box optimization software, stakeholders access a feasible solution at each iteration of the algorithm and are then aware of the cause of descent increment limitation, allowing them to ultimately simplify the road model if deemed relevant. By using this property, we have been able to demonstrate that slopes were the most sensitive degrees of freedom for the road project studied. This finding will be used to simplify the optimization on the next studied cases. Using this new efficient algorithm, different scenarios of traffic evolution and a sensitivity analysis with respect to the main parameters of constructions were carried out. Results are consistent and demonstrate the stability of the optimization. One limitation of the proposed algorithm is its dependence on the initial number of ramps established in the initial design. Two solutions are envisaged to overcome this limitation. The simplest way would be to present a meta level analysis of the algorithm to the stakeholders that exposes the cause of this descent step limitation, in offering them insights into modifying this initial number of ramps. The second, more elegant solution would entail reformulating the problem within the variational analysis framework without the guarantee to obtain a feasible solution. An open-source software based on this paper is currently under development in order to popularize Sloop among practitioners.