A note on: Spline technique for modeling roadway profile to minimize earthwork cost

The gap constraint used in A. A. Moreb, Spline technique for modeling 
 roadway profile to minimize earthwork cost, J. Ind. Man. & Opt., 5 (2) (2009), 275-283 introduces unnecessary errors, while the slope 
constraint may be violated for second- and higher-order splines. In this note we amend the gap constraint, while 
maintaining the linearity of the model. We also present an improved slope constraint for 
linear and quadratic splines, and show that it becomes nonlinear for cubic and 
higher order splines. The improvements also apply to A. Moreb, M. Aljohani, 
 Quadratic representation for roadway profile that minimizes earthwork cost, 
 J. Sys. Sci. & Sys. Eng., 13 (2) (2004), 245-252.


1.
Introduction. Based on Easa [2], Moreb [6] introduced a linear programming (LP) model that combines grade selection, earthwork, and transportation using a piecewise linear approximation of the road profile. In [8], Moreb and Aljohani extended the model by using a quadratic approximation. The quadratic representation was then further extended in [7] with a spline approximation and tested using second and third degree polynomials. The result was an improved and efficient linear model that avoids sharp connectivities and scaling problems of piecewise linear representations. However, most of the constraints from [6] were adopted in the models in [8] and [7]. In this process, errors were introduced in the gap constraint and the slope constraint. In this paper we will discuss the affected constraints and suggest some improvements to avoid the errors.
Valentin Koch was partially supported by the Ministry of Advanced Education of the Province of British Columbia through a Pacific Century Graduate Scholarship. Yves Lucet was partially supported by a Discovery grant from the Natural Sciences and Engineering Research Council of Canada.

VALENTIN R. KOCH AND YVES LUCET
2. Current model. For consistency, we use the same notation as in [7] with minor adjustments. The model divides a road into sections. The vertical ground profile in each section i is represented as a point (x i , h i ), where x i is the midpoint of section i on the x-axis and h i is the average height of the ground profile in section i. The road profile is modeled by a spline function that is divided into m spline segments. Each spline segment spans one or multiple sections. A section is indexed either directly with i ∈ {1, . . . , n} or as g, i through its segment g ∈ {1, . . . , m} and section within that segment i ∈ {1, . . . , n g }. The parameter n ∈ Z + is the number of sections in the road; m ∈ Z + is the number of spline segments in the road; n g ∈ Z + is the number of sections in segment g, g ∈ {1, 2, . . . , m}; A i ∈ R + is the surface area of section i, i ∈ {1, 2, . . . , n}; x i ∈ R + is the midpoint of section i (also indexed x g,i ); h i ∈ R + is the average height of the ground profile in section i; m g ∈ Z + is the order of the polynomial in segment g; t g ∈ R + is the x-coordinate at the end of segment g; s g,i ∈ R + is the x-coordinate at the beginning of section i in segment g; c ij ∈ R + is the cost per cubic unit for excavation at section i, hauling from section i to j, and embankment at section j where i, j ∈ {1, 2, . . . , n}, i = j.
The decision variable u i ∈ R + is the height of earth to be removed from section i; v i ∈ R + is the height of earth to be added to section i; T ij ∈ R + is the cubic units of earth to be transported from section i to section j where i, j ∈ {1, . . . , n}, i = j; a gk ∈ R is the coefficient of x k in the polynomial of segment g, k ∈ {0, . . . , m g }.
The following constraint balances the volume transported to or from a section with the work done to that section.
The gap constraint ensures that the gap between ground profile and road profile in a section matches the work done to that section. It is where P g (x) = mg k=0 a gk x k , and x gi represents the midpoint of the i-th section in segment g.
In order to comply with road specifications, upper bounds U and lower bounds L are imposed to the slope with While the constraints for zero order continuity of the spline are not listed explicitly in [8,7], they are implicit. For completeness, we state them here as P g (t g ) = P g+1 (t g ), g = 1, 2, . . . , (s − 1).
At last, C 1 smoothness is required for the spline.
Additional constraints mentioned in the paper are not relevant for the present note.
3.1. The gap constraint. In Equation (3), the difference between the height of the averaged ground profile in section i of segment g and the polynomial of the spline at the midpoint x gi of that section is computed. This distance is required to equal the height of earth that was removed or added to that section. This constraint works if the spline P g is linear on segment g, as in [6]. If the spline is linear, the line that represents the average of the polynomial in section i of segment g will cut the polynomial exactly at the midpoint x gi . So when m g = 1, the constraint will balance the volume between the average of the ground profile and the average of the polynomial with the volume that is excavated or embanked. If, however, m g > 1, the evaluation of the polynomial at x gi may not represent the actual height of the average of the spline in that section since the average of the polynomial may intersect the polynomial at some other point than x gi . In Table 1, we show the error that is introduced in road section 1 for case 2 and 3 in the numerical example in [7]. We use the coefficients given in the paper. Let P denote the spline andP 1 denote the average of the spline over section 1. From the table we see that the errors in these examples are small because the spline is almost piecewise linear (the coefficient of quadratic parts are around 10 −5 , and the cubic part around 10 −8 ). Since most roads require smooth and gradual change in grade along the vertical curve [1], we do not expect a significant error in most cases.
However, there is no need to introduce this error: if we let s gi denote the xcoordinate of the beginning at the i-th section in segment g, we want which results in As we can see, equation (8) is linear in the decision variables a gk , u gi , and v gi . We can therefore replace the constraints in (3) with (8) and still solve the problem as a linear program without introducing the error mentioned above.
3.2. The slope constraint. One of the most important factors in road design is the stopping sight distance (SSD) [1,4,5]. In the case of parabolic vertical curves (crest or sag), the minimum length of the curve depends on the SSD and the change of road grade [5]. A violation of the maximal allowed slope may therefore result in a violation of the SSD and can greatly impact the safety. In [7], road specifications may include an upper and lower bound on the slope of the road, which is modeled with constraint (4). In inequality (4), the bounds are applied to the derivative of the spline evaluated at the midpoint x gi of the i-th section in segment g. However, this constraint will not ensure that the road complies with the specifications. If m g > 1, the maximal slope of the polynomial within a section can be attained at the endpoints or at inflection points inside the section, not just at the midpoint x gi of that section (this was already pointed out by [5] on a similar model with a 5 th order polynomial).  Table 2, we show the violation that occurs in the example in [7] for case 1 without the mid-road constraint from [2] and grade specifications of ±1%. All other data is taken from [7] where P denotes the spline and x * ∈ [0, 4000] the extrema.
We therefore suggest an exact approach using minimization and maximization over the full spline: where t 0 = 0. Consider section (g, i) and name I = [t g−1 , t g ]. If x * ∈ I is an extrema of P g (x) over I, then either x * is an endpoint of I or P ′ g (x * ) = 0. If x * is an endpoint, we can enforce the slope constraint with the inequalities L ≤ mg k=1 ka gk t k−1 g ≤ U, g = 1, 2, . . . , s, which are linear in the decision variables a gk . Note that one additional constraint is needed for the first spline segment that is evaluated at t 0 . In order to enforce the slope constraint when x * is not an endpoint, we want x * ∈ [s gi , s gi+1 ] and If the polynomial over I is of degree 2, x * is one of the endpoints of I and we can use the linear constraints in (11) and solve the problem as a linear programming problem. However, if we have a polynomial of degree 3, the second derivative becomes 2a g2 + 6a g3 (x * ) = 0 and x * = −a g2 /(3a g3 ), resulting in a nonlinear problem. Consequently, the slope constraint in [7] could be violated for cubic (and higher order splines) with potentially large errors for polynomials that are very nonlinear.

Experimental setup.
In this section, we use a quadratic spline to show the differences between the solution of the original model and the solutions with the corrected constraints (8) and (11). We use a piecewise linear ground profile that has a similar shape as the one in [7]. The road used is a rural collector with a design speed of 100 km/h and four lanes. The data for the ground profile is given in Table 3. All values are given in meters. As in [7], we require the initial and final road height to be y init = 32 m and y final = 31 m respectively. We do not require a road height anywhere else in our example. We let A i = 2880 m 2 , i = 1, . . . , n. The parameters x i , h i , i = 1, . . . , n can be computed with the given values above. We use the same cost model as in [7] converted to metric values. Thus, using different costs for excavation, haul, and embankment, we have c ij = 3.139 + 1.373|i − j| + 7.063, i, j = 1, . . . , n, i = j.
The spline contains 3 segments. The data for the spline is given in Table 4. Let S denote the minimum required stopping sight distance. According to the guidelines of the American Association of State Highway and Transportation Officials (AASHTO) and based on our design speed, we let S = 205 meters [3]. With the given S, we can now compute the maximum slopes of the spline. For simplicity, we only compute the values for the first segment in the spline. In [1], the formula for the length L of a vertical crest curve, where L ≥ S, is given as where η 1 is the height of the drivers eye, η 2 the height of the sighted object, and Λ the absolute change of gradient before and after the crest curve. Since in the first spline segment, L = 400 meters, we have indeed L ≥ S. Recommended values by AASHTO for the heights are η 1 = 1.07 m and η 2 = 0.15 m [3]. We solve (13) for Λ and use the recommended values for η 1 and η 2 and end up with Λ = 0.0384 or 3.84%. Now Λ = |ζ 1 − ζ 2 |, where ζ 1 is the gradient before the crest curve and ζ 2 is the gradient after the crest curve. If, for simplicity, we let |ζ 1 | = |ζ 2 |, and |U | = |L|, we have Similar computations can be done for the sag curve in the second segment and the crest curve in the third segment. For the numerical experiment, however, we use L = −0.019 and U = 0.019 for all segments.

4.2.
Results. We run three different experiments with the following configurations.
1. The original model from [7]; 2. The model from [7] with the corrected gap constraint (8) applied; 3. The model from [7] with the corrected gap constraint (8) and the corrected slope constraint (11) applied. The results for the three cases with the corresponding polynomial for the first spline segment are given in Table 5. The different road profiles are shown in Figure 1. In Table 6, we show the absolute and relative error that is introduced into the objective value with case 1 with respect to case 2 and 3. The relative error that is introduced by the gap constraint is about 9%, which, in the context, is not unacceptable. Remember that we approximate the change in volume of earth as a linear function of height, where as in reality, it is a quadratic function of height [5]. However, the relative error of 43%, that is introduced with the wrong slope constraint,  is of much greater concern. Thus, if the road is built according to the correct specifications, the price estimate will change significantly and the original model is not a good approximation anymore.
The last part of our analysis concerns the violation of the stopping sight distance. With the given design constraints in Section 4.1, we computed a maximal grade change Λ = 0.0384. Table 7 shows the violation of the maximal grade change. The original model violates the maximal grade change Λ by 25% and therefore clearly violates the SSD requirements. Note that we did not show an example for a cubic spline. As explained in Section 3.2, the correct application of the slope constraints to a cubic spline results in a nonlinear program. Further research is needed to determine if it can be solved efficiently.

5.
Conclusion. In this note we improve the modeling of the gap and slope constraints proposed in [7]. Earthwork in roadway construction is not a precise task. It is therefore clear that whatever model we use, it will always be an approximation to the real process and carries uncertainties. In this context, the error that is introduced in (3) might not have a dramatic effect in practice. However, since we can avoid this additional error, we suggest to use the constraint (8) instead.
The grade constraint for the road has a bigger impact, especially if the model is used for highway alignments. A small violation of the slope can result in a violation of the required minimum stopping sight distance on the vertical curve, which may greatly impact the safety. As shown in our discussion, the potential slope violations that may occur in constraint (4) can be avoided for linear and quadratic splines with equation (11). If we want to use a cubic spline, we cannot enforce the constraint with linear programming. In that case, one could use nonlinear programming. Another approach could be to solve the problem with the original constraint (4) and then check for violations afterwards. If a violation occurs in a section, the section would be split in two so that the midpoint x i of one of the new sections falls onto the point of violation. Constraints (3) and (4) would then be applied to the newly created sections and constraints (5) and (6) would be applied to the end slope of the first of the newly created sections. This process could then be repeated iteratively. However, it is not clear how many iterations it would take, and under what conditions this process would end with a solution where no violation occurs. We leave it as a further research topic.