Multi-Haul Quasi Network Flow Model for Vertical Alignment Optimization

TL;DR

This paper introduces the multi-haul quasi network flow (MH-QNF) model for vertical road alignment optimization, significantly improving accuracy, robustness, and computational efficiency over existing models.

Contribution

The paper presents a novel MH-QNF model that enhances the accuracy and speed of vertical alignment optimization compared to previous models.

Findings

MH-QNF solves over 93% of test problems with less than 1% error.

MH-QNF is approximately 8 times faster than the CTG model.

The model outperforms existing state-of-the-art models in robustness and efficiency.

Abstract

The vertical alignment optimization problem for road design aims to generate a vertical alignment of a new road with a minimum cost, while satisfying safety and design constraints. We present a new model called multi-haul quasi network flow (MH-QNF) for vertical alignment optimization that improves the accuracy and reliability of previous mixed integer linear programming models. We evaluate the performance of the new model compared to two state-of-the-art models in the field: the complete transportation graph (CTG) and the quasi network flow (QNF) models. The numerical results show that, within a 1% relative error, the proposed model is robust and solves more than 93% of test problems compared to 82% for the CTG and none for the QNF. Moreover, the MH-QNF model solves the problems approximately 8 times faster than the CTG model.

Tables

Table 1. Table 1: The specifications of the basic roads used to assemble the test collection.

Road	Length (km)	Section Length (m)	Number of sections
A	1	20	50
B	5	100	50
C	2	20	100
D	3	20	150
E	15	100	150
F	20	100	200
G	9	20	450

Table 2. Table 2: Cost components in $ / c ​ u . m formulae-sequence currency-dollar 𝑐 𝑢 𝑚 \$/cu.m .

Cost component	Cost Type	Value
Excavation cost	M1	4.000
	M2	4.000
	M3	20.000
	M4	4.000
Embankment cost	M1	2.000
	M2	2.000
	M3	1.800
	M4	2.000
Cost component	Cost Type	Value	Avg Value
Loading cost	Short	0.000	1.067
	Middle	0.600
	Long	2.600
Hauling cost	Short	0.008	0.005
	Middle	0.004
	Long	0.002

Table 3. Table 3: Summary of model accuracy for each model.

#	Model	opt. found	min. error (%)	mean error (%)	max. error (%)
1	CTG-B	51	–	–	–
2	MQN-B	56	0.00	0.47	2.20
3	QNS-B	56	27.15	57.10	100.92
4	QNM-B	56	12.33	25.65	40.74
5	QNL-B	56	18.93	25.81	37.95
6	QNA-B	56	22.19	40.22	62.98
7	CTG-S1	45	–	–	–
8	MQN-S1	57	0.00	0.37	1.43
9	QNS-S1	56	27.15	57.56	100.92
10	QNM-S1	56	12.45	25.78	40.74
11	QNL-S1	56	18.93	25.81	37.39
12	QNA-S1	56	22.36	40.44	62.98

Table 4. Table 4: Solution times (in second) for 1% relative error.

Test #	MQN-B	CTG-B	MQN-S1	CTG-S1
1	23.84	19.02	15.21	33.71
2	11.00	27.18	15.28	35.41
3	26.51	190.64	20.51	237.69
4	168.62	7214.67	172.48	4848.31
5	42.07	434.39	52.39	452.03
6	13.05	232.84	28.26	1014.80
7	874.97	603.49	581.46	1928.96
8	NaN	12605.20	2409.50	NaN
9	19.93	169.51	78.17	309.46
10	13.40	297.75	17.13	141.41
11	66.34	2348.93	60.28	NaN
12	146.36	3929.18	124.01	9853.74
13	4.96	86.70	11.38	33.63
14	2.15	39.36	4.97	17.24
15	18.75	392.67	35.29	451.63
16	319.90	NaN	366.57	NaN
17	13.81	81.61	23.76	174.21
18	238.37	1369.81	383.95	4559.33
19	30.17	799.55	68.35	1657.61
20	81.99	632.69	1175.49	4146.94
21	455.83	NaN	2743.00	NaN
22	223.27	8901.20	1426.83	NaN
23	65.00	1306.09	553.65	2437.76
24	3.65	91.19	7.87	55.34
25	9.64	750.86	25.49	314.50
26	76.04	5486.19	158.04	NaN
27	6.87	353.50	14.64	57.91
28	6.16	493.60	14.19	248.43
29	5.52	67.06	9.20	19.16
30	30.44	2454.93	51.14	2432.61
31	9.89	179.33	9.83	35.04
32	18.15	619.50	78.73	241.43
33	11.03	1361.16	25.27	675.14
34	30.08	2385.43	91.46	NaN
35	22.71	2463.64	35.74	2543.56
36	120.08	NaN	273.13	9350.41
37	57.23	1702.74	90.71	NaN
38	3.12	4.96	4.95	4.98
39	17.09	230.50	36.15	623.73
40	18020.00	NaN	NaN	NaN
41	3361.52	NaN	9224.50	NaN
42	392.59	NaN	455.71	NaN
43	308.46	NaN	313.39	NaN
44	NaN	NaN	NaN	NaN
45	70.47	2014.47	168.04	NaN
46	50.75	1289.76	122.84	3528.95
47	8.85	7.44	17.43	7.78
48	5.70	9.02	10.61	9.76
49	3.22	1.87	6.72	1.89
50	1.36	0.94	2.46	0.99
51	4.83	6.11	7.72	6.14
52	8.97	28.74	14.42	28.87
53	14.37	40.89	23.76	41.09
54	13.58	27.02	21.90	27.86
55	13.68	30.44	23.01	31.26
56	3.78	4.74	6.34	4.68
57	73.50	1373.68	94.70	1361.67
58	NaN	NaN	NaN	NaN
59	15.24	122.75	15.90	128.76
60	6.99	14.09	8.58	16.29

Equations

$$P_g(s)=a_{g,1}+a_{g,2}s+a_{g,3}{s}^2$$

$$P(s)=\{\begin{array}{ll}{P}_1(s)& \text{if }{s}_{\phi (1,1)}\le s\le s_{\phi (1,n_1)},\\ P_2(s)& \text{if }{s}_{\phi (2,1)}\le s\le s_{\phi (2,n_2)},\\ ⋮& \\ P_{n_g}(s)& \text{if }{s}_{\phi (n_g,1)}\le s\le s_{\phi (2,n_{n_g})}.\end{array}$$

$$\min\Bigg{[}\sum_{\begin{subarray}{c}i\in\mathcal{S\cup B}\\ m\in\mathcal{M}\\ h\in\mathcal{H}\end{subarray}}(p_{m}+y_{h})V_{i}^{+}+\sum_{\begin{subarray}{c}i\in\mathcal{S\cup W}\\ m\in\mathcal{M}\end{subarray}}q_{m}V_{i}^{-}+\sum_{\begin{subarray}{c}i\in\mathcal{S}\\ h\in\mathcal{H}\\ t\in\mathcal{T}\end{subarray}}(c_{i,i-1}^{r,h}f_{i,i-1,t}^{r,h}+c_{i,i+1}^{r,h}f_{i,i+1,t}^{r,h})\\ +\sum_{\begin{subarray}{c}j\in\mathcal{B}\\ m\in\mathcal{M}\\ h\in\mathcal{H}\\ t\in\mathcal{T}\end{subarray}}(p_{m}+y_{h}+c_{h}\tilde{d}_{j})(f_{j,\vartheta(j)-1,t}^{b,h}+f_{j,\vartheta(j)+1,t}^{b,h})\\ +\sum_{\begin{subarray}{c}k\in\mathcal{W}\\ m\in\mathcal{M}\\ h\in\mathcal{H}\\ t\in\mathcal{T}\end{subarray}}(q_{m}+c_{h}\tilde{d}_{k})(f_{\delta(k)-1,k,t}^{w,h}+f_{\delta(k)+1,k,t}^{w,h})\Bigg{]}.$$

$$f_{i-1,i,t}^{r,h}+f_{i,i+1,t}^{u,h}+\sum _{\begin{array}{c}j\in \mathcal{B}\\ \vartheta (j)=i\end{array}}{f}_{j,i+1,t}^{b,h}$$

$$f_{i+1,i,t}^{r,h}+f_{i,i-1,t}^{u,h}+\sum _{\begin{array}{c}j\in \mathcal{B}\\ \vartheta (j)=i\end{array}}{f}_{j,i-1,t}^{b,h}$$

$$\sum _{t\in \mathcal{T}}(f_{i,i+1,t}^{u,h}+f_{i,i-1,t}^{u,h})$$

$$\sum _{\begin{array}{c}t\in \mathcal{T}\end{array}}(f_{j,\vartheta (j)+1,t}^{b,h}+f_{j,\vartheta (j)-1,t}^{b,h})$$

$$\sum _{t\in \mathcal{T}}(f_{i-1,i,t}^{l,h}+f_{i+1,i,t}^{l,h})$$

$$\sum _{\begin{array}{c}t\in \mathcal{T}\end{array}}(f_{\delta (k)-1,k,t}^{w,h}+f_{\delta (k)+1,k,t}^{w,h})$$

$$\sum _{\begin{array}{c}t\in \mathcal{T}\end{array}}(f_{j,\vartheta (j)+1,t}^{b,h}+f_{j,\vartheta (j)-1,t}^{b,h})$$

$$\sum _{\begin{array}{c}t\in \mathcal{T}\end{array}}(f_{\delta (k)-1,k,t}^{w,h}+f_{\delta (k)+1,k,t}^{w,h})$$

$$-(f_{\gamma (i)-1,\gamma (i),t}^{r,h}-f_{\gamma (i)-1,\gamma (i),t}^{l,h})$$

$$f_{\gamma (i)-1,\gamma (i),t}^{r,h}-f_{\gamma (i)-1,\gamma (i),t}^{l,h}$$

$$-(f_{\gamma (i),\gamma (i)+1,t}^{r,h}-f_{\gamma (i),\gamma (i)+1,t}^{u,h})$$

$$f_{\gamma (i),\gamma (i)+1,t}^{r,h}-f_{\gamma (i),\gamma (i)+1,t}^{u,h})$$

$$-(f_{\gamma (i)+1,\gamma (i),t}^{r,h}-f_{\gamma (i)+1,\gamma (i),t}^{l,h})$$

$$f_{\gamma (i)+1,\gamma (i),t}^{r,h}-f_{\gamma (i)+1,\gamma (i),t}^{l,h}$$

$$-(f_{\gamma (i),\gamma (i)-1,t}^{r,h}-f_{\gamma (i),\gamma (i)-1,t}^{u,h})$$

$$f_{\gamma (i),\gamma (i)-1,t}^{r,h}-f_{\gamma (i),\gamma (i)-1,t}^{u,h}$$

$$f_{i,i+1,t}^{r,h}$$

$$\forall i\in \mathcal{S},(k_1,k_2)\in {\bar{\mathcal{I}}}^2,h\in \mathcal{H},t\in \mathcal{T},\gamma (k_1)\le i,i+1\le \gamma (k_2),$$

$$f_{i+1,i,t}^{r,h}$$

$$\forall i\in \mathcal{S},(k_1,k_2)\in {\bar{\mathcal{I}}}^2,h\in \mathcal{H},t\in \mathcal{T},\gamma (k_1)\le i,i+1\le \gamma (k_2),$$

$$f_{j,\vartheta (j)+1,t}^{b,h}$$

$$\forall j\in \mathcal{B},(k_1,k_2)\in {\bar{\mathcal{I}}}^2,h\in \mathcal{H},t\in \mathcal{T},\gamma (k_1)\le \vartheta (j)-1,\vartheta (j)+1\le \gamma (k_2),$$

$$f_{j,\vartheta (j)-1,t}^{b,h}$$

$$\forall j\in \mathcal{B},(k_1,k_2)\in {\bar{\mathcal{I}}}^2,h\in \mathcal{H},t\in \mathcal{T},\gamma (k_1)\le \vartheta (j)-1,\vartheta (j)+1\le \gamma (k_2),$$

$$f_{\delta (j)-1,j,t}^{w,h}$$

$$\forall j\in \mathcal{W},(k_1,k_2)\in {\bar{\mathcal{I}}}^2,h\in \mathcal{H},t\in \mathcal{T},\gamma (k_1)\le \delta (j)-1,\delta (j)+1\le \gamma (k_2),$$

$$f_{\delta (j)+1,j,t}^{w,h}$$

$$\forall j\in \mathcal{W},(k_1,k_2)\in {\bar{\mathcal{I}}}^2,h\in \mathcal{H},t\in \mathcal{T},\gamma (k_1)\le \delta (j)-1,\delta (j)+1\le \gamma (k_2).$$

$$f_{i,i+1,t}^{r,h}$$

$$f_{i+1,i,t}^{r,h}$$

$$f_{j,\vartheta (j)+1,t}^{b,h}$$

$$f_{j,\vartheta (j)-1,t}^{b,h}$$

$$f_{\delta (j)-1,j,t}^{w,h}$$

$$f_{\delta (j)+1,j,t}^{w,h}$$

$$f_{i,i+1,t}^{r,h}$$

$$f_{i+1,i,t}^{r,h}$$

$$f_{j,\vartheta (j)+1,t}^{b,h}$$

$$f_{j,\vartheta (j)-1,t}^{b,h}$$

$$f_{\delta (j)-1,j,t}^{w,h}$$

$$f_{\delta (j)+1,j,t}^{w,h}$$

$$\sum _{t=0}^u(f_{\gamma (k),\gamma (k)+1,t}^{u,h}+f_{\gamma (k),\gamma (k)-1,t}^{u,h})+M(1-y_{k,u})$$

$$\sum _{t=0}^u(f_{\gamma (k)-1,\gamma (k),t}^{l,h}+f_{\gamma (k)+1,\gamma (k),t}^{l,h})+M(1-y_{k,u})$$

$$\sum _{t=0}^u\sum _{k\in \mathcal{I}}{y}_{k,t}\ge t+1\forall u\in \setminus \{n_b\}.$$

$$y_{k,t}\ge y_{k,t-1},\forall k\in \mathcal{I},t\in \setminus \{0\}.$$

Full text

1 Introduction

Transportation infrastructures are an important sign of development and welfare in a given country Celauro et al. (2015); Herranz-Loncán (2007). While planners normally try to balance and coordinate among various transportation methods, roads and highways are considered the leading component of transportation infrastructures. Due to the fast increase of traffic volume and loadings, high costs of road construction projects, and environmental and safety impacts, road construction require continuing innovation. Road design is one of the early and significant steps in road construction projects and it influences the construction cost and other contributing constraints significantly.

Road design starts by selecting the corridor within which the road is to be constructed. From here, a so called horizontal alignment is fixed. Next vertical alignment and earthwork movement problems must be solved. In the literature, these last two steps are often combined for the sake of efficiency Hare et al. (2011); Moreb and Aljohani (2004). Some studies consider the simultaneous optimization of the horizontal alignment, vertical alignment, and earthwork movement Chew et al. (1989); Easa et al. (2002); Hassan et al. (1998); Cheng and Lee (2006). Our motivation to focus solely on the last two problems is threefold. First, the horizontal alignment problem is very complex and can currently only be solved by heuristics that do not guarantee global optimality. Moreover, the resulting roads generally still require a final detailed optimization of the vertical alignment and earthwork movement problems Lee et al. (2009). Next, when upgrading a road, the horizontal alignment is already decided and only the vertical and earthwork problems need to be solved. Finally, building a deterministic global optimization algorithm to solve the combined vertical alignment and earthwork problem is a serious step towards solving the complete problem. Specifically, such algorithms can be used in bilevel formulations that are designed to simultaneously solve horizontal, vertical and earthwork problems, as proposed in Mondal et al. (2015).

Heuristic optimization methods, such as Genetic Algorithms (GAs), have been used to solve the road alignment problem. Lee. and Cheng (2001) developed a mathematical model for vertical alignment and proposed a heuristic method to solve it. Jha et al. (2006) present a comprehensive review of the road design problem and the related challenges. They also present solution algorithms based on GAs to optimize highway alignments by adapting operators and encoding schemes. Other heuristic-based optimization methods for highway alignment optimization include Goktepe et al. (2009); Jong and Schonfeld (2003); Kim et al. (2005); Aruga (2005); Kang et al. (2012).

Another early and popular approach relied on dynamic programming. Dynamic programming was first employed Goh et al. (1988) to develop two models for vertical alignment. They use piecewise-linear segments to represent the alignment. Based on the results of the paper, the state parametrization model is more flexible to be used for solving a three-dimensional problem while the dynamic programming model has a better performance in terms of computational time and ease of formulation. In general, compared to optimizing the horizontal alignment, a dynamic programming approach has a better performance in vertical alignment optimization; however, the resulting alignment is limited to a finite set of points at each station. Hence, only a portion of the problem search space is explored Jong (1998). Other models using dynamic programming include Trietsch (1987); Nicholson et al. (1975); Li et al. (2013).

Our focus is to minimize the road construction costs (primarily excavation, embankment, and hauling costs), subject to road safety and quality constraints. While costs of land, pavement, vehicle operation, motorist time, accidents, and environmental impacts all contribute to the cost, especially when one considers highways, they have less impact when the horizontal alignment is fixed, or for logging roads. Optimizing solely road construction costs allows us to simplify our model and to take advantage of specific problem structures to speed up the computation.

The current state-of-the-art on the vertical alignment and earthwork movement problem relies on linear programming, or mixed integer linear programming. Easa (1988) proposed a two-step linear programming model to optimize vertical alignment and earthwork movement. In the first step by enumerating all feasible grades, a feasible road grade is chosen and for each grade the amount of cut and fill is calculated. In the second step the earthwork transportation cost is minimized using linear programming. This trial and error approach does not guarantee finding the global solution. Moreb (1996) proposed a linear programming model for vertical alignment which incorporates the grade selection and earthwork allocation stages. Moreb and Aljohani (2004) extended the work of Easa (1988) to handle the sharp connectivity points problem in piecewise linear models by representing the road profile as a continuous quadratic curve. Moreb (2009) proposed a model in which the road profile is represented as a one-dimensional spline. This model improves the computational efficiency and guarantee the global optimality by solving a single linear programming problem. Koch and Lucet (2010) improved the Moreb (2009) model by introducing new gap and slope constraints that reduce errors in the model.

Hare et al. (2011) developed a mixed integer linear program (MILP) model for vertical alignment in which the blocks are also taken into account. In Hare et al. (2015) road side slopes are added to the previous MILP model proposed in Hare et al. (2011) which yields a more accurate earthwork cost calculation. Hare et al. (2014) also improved the performance of the MILP model (called complete transportation graph or CTG) and studied a novel quasi network flow (QNF) model for vertical alignment optimization considering the earthwork allocations. The QNF model presented significant improvements in computational time, but at a cost of decreased accuracy in the final solution.

In this paper, we improve the QNF model of Hare et al. (2014), by presenting a new model called the multi-haul cost quasi network flow (MH-QNF) model. The MH-QNF model balances the accuracy of the CTG model in approximating the earthwork cost with the low computational cost of the QNF model. The main idea behind the new model is to differentiate between hauls of different distances in earthwork allocation. For simplicity, we consider three possible distances, short hauls, middle hauls, and long hauls; however, any number of haul distances can be accommodated by the model.

2 Model description

The road is approximated as a quadratic spline. It is split into $m$ segments indexed by $\mathcal{G}=\{1,2,...,m\}$. For all $g\in\mathcal{G}$, each segment is represented using the following equation

[TABLE]

in which $s$ is the distance from the beginning of the road.

Each segment is further split into sections so that the $g$th spline segment is made of $n_{g}$ sections indexed by the set $\mathcal{S}_{g}=\{1,2,...,n_{g}\}$. The total number of sections in a road is $n=\sum_{g\in\mathcal{G}}n_{g}$ and these sections are indexed by the set $\mathcal{S}=\{1,2,...,n\}$. The section indexes of a specific spline segment $g$ are mapped to the actual section index set $\mathcal{S}$ using the function $\varphi:(\mathcal{G}\times\mathcal{S}_{g})\rightarrow\mathcal{S}$. For example, if $\varphi(g,j)=i$ then $s_{i}=s_{\varphi(g,j)}$ for all $i\in\mathcal{S},g\in\mathcal{G},j\in\mathcal{S}_{g}$. So the spline function representing the road profile is

[TABLE]

In the vertical alignment problem, one goal is to obtain the optimal offset between the ground profile and the road profile. In addition, the required cut and fill volumes should be calculated for each offset. The offset is denoted by $u_{i},i\in\mathcal{S}$, and the cut and fill volumes of a section $i$ are denoted by $V^{+}_{i}$ and $V^{-}_{i},i\in\mathcal{S}$ respectively. (Shrinkage and swell factors for each material can be easily incorporated in the model but are not included to simplify the presentation.)

When constructing a road we need to consider borrow pits to bring in material and waste pits to dump extra materials. They are modeled as external sections and indexed by the sets $\mathcal{B}=\{1,2,...,n_{\beta}\}$ and $\mathcal{W}=\{1,2,...,n_{w}\}$ in which $n_{\beta}$ is the number of borrow pit sections and $n_{w}$ is the number of waste pit sections. A borrow pit index is mapped to the section index to which it is attached by the function $\vartheta:\mathcal{B}\rightarrow\mathcal{S}$. Similarly, the function $\delta:\mathcal{W}\rightarrow\mathcal{S}$ maps the waste pit index to the section index to which it is attached.

Access roads are required in road construction and they are used as gateways to the road being constructed. Access roads are linked to a section of the road and are indexed by the set $\mathcal{R}=\{1,2,...,n_{r}\}$ in which $n_{r}$ is the number of access roads in the construction site. The function $\varrho:\mathcal{R}\rightarrow\mathcal{S}$ maps an access road index to the section index to which it is linked. The capacity of the $i$th borrow pit (respectively waste pit) is denoted by $C^{b}_{i}$ (respectively $C^{w}_{i}$). The set $\mathcal{N}=\mathcal{S}\cup\mathcal{B}\cup\mathcal{W}=\{1,2,...,n+n_{\beta}+n_{w}\}$ represents all the indexes for sections, borrow pits, and waste pits. The dead haul distance is the distance between a borrow/waste pit and the section to which it is linked. The dead haul distance of the $i$th borrow or waste pit is denoted by $\tilde{d}_{i}$.

The side slope is defined as the steady decrease/increase of height when moving orthogonally from the road profile to the ground profile in a cut/fill section. Side slopes are represented as trapezoid shaped cross-sections and are approximated using stacked rectangles to preserve the linearity of the model, see Hare et al. (2014). The model can handle nonsymmetric section-dependent side slopes to account for local soil type and sloped terrains as well as cuts and fills at the same section.

Another important consideration in vertical alignment optimization are blocks Hare et al. (2011). Blocks are obstacles that need to be removed to access some parts of the corridor. They indicate a river or mountain, and require building a bridge or a tunnel. We define $n_{b}$ as the number of blocks in the corridor. Blocks are indexed by the set $\mathcal{I}=\{1,2,...,n_{b}\}$. The function $\gamma:\mathcal{I}\rightarrow\mathcal{S}$ maps a block index to the section index to which it is linked. To model the block removal process we use a time step $t$ which specifies the time at which a block is removed. In Hare et al. (2011), it is shown that to remove $n_{b}$ blocks we need at most $n_{b}+1$ time steps. So we define the set $\mathcal{T}=\{0,1,2,...,n_{b}\}$ to model the required time steps and we use the binary variables $y_{kt}$ for each block $k\in\mathcal{I}$ and time step $t\in\mathcal{T}$ to determine whether a block is removed or not. After a block is removed we can move material over its section, see constraints (2)–(9).

In the MH-QNF model, similarly to the QNF model Hare et al. (2014), we draw sections and pits as nodes while arcs show feasible movements between them. In the QNF model the authors use a single haul to move materials. This is equivalent to assuming a single type of earthmoving equipment. In real roadway construction sites, different equipments are required for earthmoving tasks. For example, for the short distances using a bulldozer may be more economical, while for a long distance movement a truck may be preferred. Therefore, to obtain a more realistic model, we extend the previous QNF model Hare et al. (2014) by using multiple hauling paths, which is more realistic and yields a more precise earthwork solution.

Figure 1 shows the typical flows for a section in the proposed MH-QNF model.

In this model, the material can be moved from the current section $i$ to the next section ($i+1$) or the previous one ($i-1$). Depending on the distance of a hauling task, the material movement can be performed through one of a number of hauling paths indexed by $\mathcal{H}=\{1,2,...,n_{h}\}$. (In our figures and numerical tests we shall apply $n_{h}=3$, which is based on consultation with our industry partner, whose standard software includes free-haul, over-haul, and end-haul.)

When the flow of material reaches a section, materials can be unloaded to fill the section, a cut from the section can be performed with the new materials being added to the flow, or the section remains unchanged and the material is transferred to the next section. To transfer material to the left or right we use virtual transit nodes for both left and right directions. Since we have $n_{h}$ different hauls, $n_{h}$ groups of transit nodes are denoted in the model, see Figure 1 (which represents three hauls for consistency with numerical experiments). If material is moved to the next section, the left-side transit nodes are employed; otherwise if the material is moved to the previous section the right-hand side transit nodes are used. In Figure 1, we define $f_{i,i+1,t}^{r,s}$ and $f_{i,i-1,t}^{r,s}$ (for all $i\in\mathcal{S},t\in\mathcal{T}$) as the flow of material from the left (right) short haul transit node $i$ to the left (right) short haul transit node $i+1$ (or $i-1$) at time step $t$, respectively. Similarly, we use the variables $f_{i,i+1,t}^{r,m}$ and $f_{i,i-1,t}^{r,m}$ (for all $i\in\mathcal{S},t\in\mathcal{T}$) for the middle hauling path, and the variables $f_{i,i+1,t}^{r,l}$ and $f_{i,i-1,t}^{r,l}$ (for all $i\in\mathcal{S},t\in\mathcal{T}$) for the long hauling path.

Figure 2 shows the load and unload flows that model the fill and cut volumes of materials for each section.

After cutting the materials we can transfer them to the left or right using any of the three hauling options (short, middle, or long). Similarly, a section can be filled using the materials taken from left or right transit nodes. To model these possibilities, we introduce the variables $f_{i-1,i,t}^{l,s}$ and $f_{i+1,i,t}^{l,s}$ (for all $i\in\mathcal{S},t\in\mathcal{T}$) as the load (fill) flows of materials from the left and right transit nodes for the short hauling path. The variables $f_{i-1,i,t}^{l,m}$, $f_{i+1,i,t}^{l,m}$; and $f_{i-1,i,t}^{l,l}$ and $f_{i+1,i,t}^{l,l}$ (for all $i\in\mathcal{S},t\in\mathcal{T}$) are used as the load flows of materials for the middle and long hauling paths, respectively. Similarly, the variables $f_{i,i+1,t}^{u,s}$ and $f_{i,i-1,t}^{u,s}$ (for all $i\in\mathcal{S},t\in\mathcal{T}$) are the unload (cut) flows of materials to the left and right transit nodes for the short hauling path. In the same way, the variables $f_{i,i+1,t}^{u,m}$, $f_{i,i-1,t}^{u,m}$; and $f_{i,i+1,t}^{u,l}$, $f_{i,i-1,t}^{u,l}$ (for all $i\in\mathcal{S},t\in\mathcal{T}$) are the unload flows for the middle and long hauling paths, respectively.

Each borrow pit is attached to a section. A borrow flow denotes material movement from the pit to that specific section; from there the material may move to a transit node. The variables $f_{j,\vartheta(j)+1,t}^{b,s}$ and $f_{j,\vartheta(j)-1,t}^{b,s}$ (for all $j\in\mathcal{B},t\in\mathcal{T}$) can be used to transfer materials from a borrow pit to the next or the previous section via the short haul path. Similarly, to represent the borrow flows for the middle and long hauling paths we define the variables $f_{j,\vartheta(j)+1,t}^{b,m}$, $f_{j,\vartheta(j)-1,t}^{b,m}$; and $f_{j,\vartheta(j)+1,t}^{b,l}$ and $f_{j,\vartheta(j)-1,t}^{b,l}$ (for all $j\in\mathcal{B},t\in\mathcal{T}$), respectively. Waste flows are handled similarly.

A typical real-world road construction project involves several types of materials with different excavation, embankment, and hauling costs. The proposed model is able to handle more than one material. We define the index set $\mathcal{M}=\{1,2,...,n_{m}\}$ corresponding to the indexes of four types of materials. In the present research we chose $n_{m}=4$ based on consultation with our industry partner.

The cost components of a vertical alignment problem consist of the embankment, excavation, and hauling costs of materials in the construction site. The per unit of volume excavation (cut or unload from a section) cost of material is represented by $p_{1}$, $p_{2}$, … $p_{n_{m}}$ for different types of materials, respectively. Similarly, $q_{1}$, $q_{2}$, … $q_{n_{m}}$ denote the per unit volume of embankment (fill or load to a section) cost of material for the different types of materials. Finally, for the short hauling path, the hauling cost of materials from the section $i$ to the section $i-1$ (respectively $i+1$) is defined as $c_{i,i-1}^{r,s}=c_{s}d_{i,i-1}$ (respectively, $c_{i,i+1}^{r,s}=c_{s}d_{i,i+1}$) in which $c_{s}$ is the cost of moving one unit volume of materials per unit distance for the short hauling path, and $d_{i,j}$ is the distance between sections $i$ and $j$. Similarly, for the middle hauling path we define the hauling cost components $c_{i,i-1}^{r,m}=c_{m}d_{i,i-1}$ and $c_{i,i+1}^{r,m}=c_{m}d_{i,i+1}$; and for the long hauling path we have $c_{i,i-1}^{r,l}=c_{l}d_{i,i-1}$ and $c_{i,i+1}^{r,l}=c_{l}d_{i,i+1}$. There is another cost when cutting materials from a section or a borrow pit and it is the loading cost, which is equipment dependent. Therefore, we define the loading cost components $y_{s}$, $y_{m}$, and $y_{l}$ as the cost of loading different types of construction equipments which use short, middle, and long hauling paths, respectively.

In the MH-QNF model the objective is to minimize the total cost of excavation, embankment, and hauling tasks. Therefore, the objective function is defined as follows.

[TABLE]

Since the transit nodes are only intermediate virtual points that provide transits for the flows of materials, the sum of the input flows to a transit node must be equal to the sum of the output flows from the node. Therefore, for all $i\in\mathcal{S},h\in\mathcal{H},t\in\mathcal{T}$ the flow constraints are defined as

[TABLE]

The section nodes are treated as the source and destination of the flows. Therefore, the total sum of load flows to a section or a waste pit should be equal to the total fill volume of materials from that section or waste pit. Similarly, the total sum of unload flows from a section or a borrow pit should be equal to the total cut volume of materials from that section or borrow pit. These constraints are called balance constraints and defined as

[TABLE]

The borrow and waste pit constraints ensure that the total sum of borrow/waste flows from/to a borrow/waste pit must not exceed the capacity of that borrow/waste, i.e.,

[TABLE]

Before introducing block constraints we require some additional notation. Let $\bar{\mathcal{I}}^{2}=\{(k_{1},k_{2})\in\mathcal{I}\times\mathcal{I}:k_{1}<k_{2},\text{ and }k_{1},k_{2}$ are two consecutive blocks with no access road in between them}, $\bar{\mathcal{I}}_{\leftarrow}=\{k\in\mathcal{I}:$ the set of blocks $k$ before which there is no access road}, and $\bar{\mathcal{I}}_{\rightarrow}=\{k\in\mathcal{I}:$ the set of blocks $k$ after which there is no access road}. The parameter $M_{i}$ is used to model binary constraints; in theory it can be any large enough number but for numerical reasons it is set to the largest possible cut or fill volume at a section $i$. For the sake of simplicity, we use $M$ instead of $M_{i}$.

Earth movement is not allowed between two sections when there is a block between them. This means that for each hauling path $h\in\mathcal{H}$, if a section $i$ is a block then material cannot be moved over that block or section. In other words, the whole material moved to the section $i$ from the section $i-1$ must be loaded to section $i$; i.e., for the left transit nodes $\forall h\in\mathcal{H},$ $f_{i-1,i,t}^{r,h}=f_{i-1,i,t}^{l,h}$. Moreover, if the section $i$ is a block the only material that can move from this section to the next section ($i+1$) is the material that has been excavated from section $i$, i.e., for the left transit nodes $\forall h\in\mathcal{H},$ $f_{i,i+1,t}^{r,h}=f_{i,i+1,t}^{u,h}$. In the case of the right transit nodes, $f_{i+1,i,t}^{r,h}=f_{i+1,i,t}^{l,h}$ and $f_{i,i-1,t}^{r,h}=f_{i,i-1,t}^{u,h}$ ensure that no material can move across section $i$ to the previous section. Therefore, for the left transit nodes, the constraints

[TABLE]

ensure that no material movement is allowed across a block. In the same way, for the right transit nodes the constraints

[TABLE]

ensure no material movement across a block.

Material movement is also not allowed between two blocks, before the first block, and after the last block with no access roads, until the blocks are removed. Thus, for all $h\in\mathcal{H}$, the transit, borrow, and waste flows are forbidden between two blocks with no access roads, until one of the blocks is removed, i.e.,

[TABLE]

After the last block with no access road, the transit, borrow, and waste flows are not allowed until the block is removed:

[TABLE]

In the same way, before the first block with no access road, the transit, borrow, and waste flows are not allowed until the block is removed:

[TABLE]

After excavating/embanking the required amount of earth from/to a section with a block, the block is considered removed. The block removal indicator constraints are used to satisfy this expectation as follows

[TABLE]

At the end of each time step $t\in\mathcal{T}$, at least one block should be removed so that removing all the blocks does not take more than $n_{b}+1$ time steps. The block removal enforcement constraint guarantees this,

[TABLE]

In addition, when a block is removed it should remain removed. This is done by the monotonicity constraint as follows

[TABLE]

The continuity constraints,

[TABLE]

ensure that the height and the slope of the first section of a segment is equal to the height and the slope of the last section of the previous segment.

The volume constraints

[TABLE]

where $A_{i}$ is the area of section $i$, guarantee that the total excavated/embanked volume of material from/to a section equals to the volume difference between the road profile and the ground profile for that section.

There are restrictions on the grade of the road profile to satisfy safety considerations. Slope constraints

[TABLE]

are used to constrain the spline segments within a closed interval $[G_{L},G_{U}]$, in which $G_{L}$ and $G_{U}$ are the minimum and the maximum valid grades, respectively.

Finally, We need bound constraints to restrict the domain of variables with

[TABLE]

and

[TABLE]

2.1 The CTG and QNF Models

To evaluate the performance of the MH-QNF model, we compare the efficiency, robustness, and accuracy of the proposed model with the CTG and QNF models presented in Hare et al. (2014). At this point, the QNF model of Hare et al. (2014) is quite easy to describe – the QNF model is a special case of the MH-QNF model, where the number of haul types is one ($n_{h}=1$). The CTG model, which appears in Hare et al. (2014, 2015), is slightly harder to outline. In essence, instead of hauling paths, for every pair of nodes $(i,j)$, the CTG model places an arc moving from $i$ to $j$. The flow through the arc is a continuous variable and the arc cost is computed based on $i$ and $j$. The CTG model allows a great deal of flexibility in arc-costs, but the resulting number of variables is very large. In particular, in the MH-QNF model (and hence the QNF model) the total number of variables grows linearly with the number of sections. Conversely, for the CTG model, the total number of variables grows quadratically with the number of sections. As the number of sections is typically the dominating factor in these models, the CTG model is typically much larger than the MH-QNF model.

3 Numerical results

To evaluate the performance of the MH-QNF model, we compare the efficiency, robustness and accuracy of the proposed model with both CTG and QNF models presented in Hare et al. (2014) in terms of timing and optimal cost.

The selected problem test collection consists of $60$ problems. These problems are generated by changing various parameters (i.e., number of sections, blocks, access roads, offset levels and section lengths) in 7 distinct road samples, provided by our industrial partner (see Table 1). All tests represent realistic road design problems.

We consider 5 hours timeout for all the experiments; beyond this a test is reported as unsuccessful. Linear programming feasibility tolerance and relative mixed integer programming gap tolerance parameters are set to $10^{-06}$ and 1%, respectively. The upper bound and lower bound for slope constraints are 0.1 and -0.1 ($\pm 10\%$). All other parameters are set to their default value.

After selection of the test collection, we solved the problems by different models using the same computational environment. The workstation used for experiments has an Intel(R)Xeon(R) CPU W3565 3.20GHz (4 cores and 8 threads) processor and 24.0 GB of RAM. To solve the developed mathematical models we use the academic version of the IBM ILOG CPLEX optimizer 12.51 edition http://www.cplex.com.

One of the influencing factors on our experiments is the cost components used in the MH-QNF model shown in Table 2 (specific values were selected in consultation with our industry partner). Excavation and embankment cost components depend on the different types of materials. Loading and hauling cost components depend on the type of hauling path. The Avg Value column is the average of the cost values for a cost component. All the three models are multi material and share the same excavation and embankment cost components for each material. Both MH-QNF model and the CTG model use multiple hauling paths with different hauling and loading costs. Therefore, in the case of the QNF model which has a single hauling path, for the loading and hauling cost components we need to test the model with different pairs of hauling and loading cost components.

Since the MH-QNF model is an extension of the QNF model, to get a better understanding of the importance of loading and hauling cost components, in addition to the average values of loading and hauling costs (1.067, 0.005, respectively), we test the QNF model with the other three pairs of cost values for loading and hauling cost components reported for different hauls (Table 2); i.e., we run the QNF model four times with different loading and hauling costs: average case (1.067, 0.005), short haul (0.000, 0.008), middle haul (0.600, 0.004), and long haul (2.600, 0.002). Thus, we have 4 different setups for the QNF model as follows.

• •

QNF-S: the model QNF configured with short hauling and loading cost values.

• •

QNF-M: the model QNF configured with middle hauling and loading cost values.

• •

QNF-L: the model QNF configured with long hauling and loading cost values.

• •

QNF-A: the model QNF configured with average hauling and loading cost values as shown in Table 2.

As shown in Figure 3, to evaluate the efficiency of using a variety of hauling paths over a single path in a road construction site, the loading and hauling cost components are selected such that

• •

if $\text{ hauling distance}<150m$ then short hauling path will be picked as the optimum method of transferring material;

• •

if $150m\leq\text{hauling distance}<1000m$ then middle hauling path will be picked;

• •

and finally, if $\text{hauling distance}\geq 1000m$ then the long hauling path will be selected.

(These specific values were selected after discussion with our industrial partner.)

In Hare et al. (2014) the authors propose six potential techniques to improve the speed of solution search. They extensively test all the techniques with both models and identify the setup which works the best for both models. Their results highlight two techniques as the most efficient methods. Therefore, we test all models using two different configurations. In the first configuration, we use SOS2 variables for the volume constraints and SOS1 variables for modeling the block constraints in all the models. In the second configuration, again we use SOS2 variables for the volume constraints but for the block constraints the basic technique is employed.

In general, in our experiments we use 12 model configurations. To increase the readability of the report on the results we use the following naming convention.

Models that use the basic technique for blocks:

• •

MQN-B: The MH-QNF model.

• •

CTG-B: The CTG model with multi haul cost components.

• •

QNS-B: The QNF-S model.

• •

QNM-B: The QNF-M model.

• •

QNL-B: The QNF-L model.

• •

QNA-B: The QNF-A model. 2. 2.

Models that use the SOS1 techniques for blocks:

• •

MQN-S1: The MH-QNF model.

• •

CTG-S1: The CTG model with multi haul cost components.

• •

QNS-S1: The QNF-S model.

• •

QNM-S1: The QNF-M model.

• •

QNL-S1: The QNF-L model.

• •

QNA-S1: The QNF-A model.

3.1 Reporting the results

As a first pass at examining model accuracy, we summarize the relative errors in cost values obtained by all models in Table 3. The column “opt. found” is the number of problems where the model optimum value was found before timeout occurred (regardless of relative error). The columns “min./mean/max. error” are the minimum/mean/maximum relative error over all problems solved when compared to CTG configurations (CTG-B or CTG-S1), based on only problems solved by both CTG configuration and the model examined.

Among the above 12 model configurations, none of the QNF model configurations are able to solve a single problem satisfactorily. We may therefore discard these methods as inaccurate.

Table LABEL:Table2 shows the solution times for the remaining 4 configurations applied to 60 test problems when only 1% relative error is allowed. “NaN” values indicate that the solver exceeded the 5 hours timeout.

The reliability and robustness of an optimization method is defined as the ability of the method to perform well over a wide range of optimization problems. To evaluate the reliability of the model we use two measures: success rate and computational accuracy. The success rate is defined as the number of problems in the given problem set that are successfully solved to optimality by the optimization method. We use performance profiles Dolan and Moré (2002) to evaluate the reliability and robustness of the models. While performance profiles are now prevalent in optimization benchmarking, we provide a brief review for completeness. The performance profile of a solver $s$ is defined as

[TABLE]

where $\lvert\mathcal{P}\rvert$ represents the cardinality of the problem set $\mathcal{P}$ and $r_{p,s}$ is the performance ratio defined as

[TABLE]

in which $t_{p,s}$ is the performance measure which is the computational time in our case. In this equation, for a specific problem and the best solver in terms of $t_{p,s}$ has $r_{p,s}=1$.

One of the advantages of performance profiles, is that they implicitly includes the performance ratio success rate as a reliability factor. The value of $\rho_{s}({\alpha})$ gives a sense of how promising the solutions found by an optimization algorithm are relative to the best solution found by all the optimization algorithms that are compared together.

For each solution obtained for problem number $s$ using the model $m$ the percent relative error is calculated using the following formula:

[TABLE]

in which $Obj_{s}$ is the optimal objective value obtained for the problem $s$ by the model $m$ and $Obj_{b}$ is the objective value obtained by the benchmark model $b$. In our experiments, a solution found by model $m$ is considered successful if $\lvert E_{s,m}\rvert\leq 1.0\%$.

When available, we consider the CTG model as the benchmark and we calculate the relative error for the solutions obtained using the MH-QNF model. When the CTG model does not provide a solution but the MH-QNF does, we consider the MH-QNF relative error to be $0.0\%$. Figure 4 provides the resulting performance profile.

Examining Figure 4, it is clear that MH-QNF greatly outperforms CTG. In terms of the block technique, the CTG model has a relatively better performance when it uses the basic technique for blocks by solving almost $82\%$ of the problems compared to the SOS1 version (CTG-S1) which ends up with solving $75\%$. Similarly, in the case of the MH-QNF model, the performance of the model using the basic technique (MQN-B) is considerably higher than the SOS1 technique (MQN-S1). Except in one test, MQN-B solves most of the problems faster than MQN-S1. In general, among all the configurations, MQN-B is the winner and has a better performance than the CTG configurations.

4 Conclusion

In this research, we extend a model based on a previous quasi network flow model by Hare et al. (2014) for the vertical alignment problem considering the earthwork cost. In particular, we add multiple hauling flows to the model in order to improve the accuracy of the model and provide more modeling flexibility for the users. The model considers features such as blocks and side-slopes. We compare the model with two state-of-the-art models in the field: CTG and QNF models.

The CTG model allows a great deal of flexibility in arc-costs, but the resulting number of variables is very large. The total number of variables grows quadratically with the number of sections. As the number of sections is typically the dominating factor in these models, the CTG model is typically much larger than the QNF and MH-QNF models.

In the QNF model, the total number of variables grows linearly with the number of sections. However, in QNF model there is a linear relation between the hauling cost and the distance traveled by the volume of materials (e.g., the cost of moving 1 unit of earth for 100m equals to 100 times the cost of moving 1 unit for 1m), which is not true in real world. Therefore, we proposed the MH-QNF model to integrate the low computational cost of QNF model with the accuracy of CTG in solving the vertical alignment problem. The experimental results confirm that the new model is both highly accurate and considerably faster than the CTG model.

These results also provide an important advancement to solving the complete road design problem (where earthwork, vertical alignment, and horizontal alignment are optimized simultaneously). Recent research by Mondal et al. (2015), has viewed the complete road design problem as a bi-level optimization problem, where the horizontal alignment is optimized using a non-gradient based method which calls the vertical alignment problem as a black-box. For this approach to work, it is critical that the vertical alignment is solved both quickly and accurately. The MH-QNF model provides both of these attributes.

Acknowledgments

This work was supported by Natural Sciences and Engineering Research Council of Canada (NSERC) under Collaborative Research and Development (CRD) Grant #CRDPJ 411318-10. The Grant was sponsored by Softree Technical Systems Inc. This work was supported by Discovery Grants #355571-2013 (Hare) and #298145-2013 (Lucet) from NSERC. Part of the research was performed in the Computer-Aided Convex Analysis (CA2) laboratory funded by a Leaders Opportunity Fund (LOF) from the Canada Foundation for Innovation (CFI) and by a British Columbia Knowledge Development Fund (BCKDF).

The authors acknowledge the original ideas and concepts of Dr. Donovan Hare that contributed to the quasi network flow model.