A mixed-integer linear programming approach for ...

10 downloads 87065 Views 2MB Size Report
Aug 25, 2016 - Mixed-Integer Linear Programming model, integrated with a cutting plane method, ... Automation in Construction 71 (2016) 314–324.
Automation in Construction 71 (2016) 314–324

Contents lists available at ScienceDirect

Automation in Construction journal homepage: www.elsevier.com/locate/autcon

A mixed-integer linear programming approach for temporary haul road design in rough-grading projects Chaojue Yi, Ming Lu ⁎ Construction Engineering and Management, Department of Civil & Environmental Engineering, University of Alberta, Canada

a r t i c l e

i n f o

Article history: Received 6 December 2015 Received in revised form 11 June 2016 Accepted 15 August 2016 Available online 25 August 2016 Keywords: Haul road layout design Mixed-integer linear programming Cutting plane method Optimization

a b s t r a c t Temporary haul road plays a pivotal part in achieving cost-efficiency and successful project delivery in heavy civil and industrial construction. Temporary haul road layout design has been empirically performed by field managers, superintendents or even truck drivers largely based on experience instead of science. Previous research endeavors in earthmoving research devised analytical algorithms to minimize the total earthmoving cost and developed simulation models to optimize earthmoving resources and processes. In the same domain, this research introduces an optimization methodology for temporary haul road layout design in order to facilitate mass earthmoving operations and improve construction performances on a typical rough-grading site. A Mixed-Integer Linear Programming model, integrated with a cutting plane method, is established for the identified problem. In particular, the cutting plane method is introduced to refine the optimization formulation by maintaining field accessibility and haul road continuity, thus ensuring practical feasibility of the analytical solution. Performances of the newly devised algorithms were benchmarked against genetic algorithms in solving ten test cases. Further, feasibility of the proposed methodology was evaluated based on a real-world case study. In conclusion, the proposed methodology is capable of tackling the temporary haul road layout design problem with high computing efficiency and delivering optimal results that are ready for field implementation. © 2016 Elsevier B.V. All rights reserved.

1. Introduction To shape the topography of the field to a specified grading design, rough grading operations are performed as the preliminary construction work on industrial and infrastructure projects, consisting of excavation, loading, haulage, unloading, grading, and compaction of earth materials [1]. In particular, site grading operations employ a large fleet of heavy equipment and represent about 25% of the total construction cost on large industrial and infrastructural projects [2]; in a separate study concerning mine haul road maintenance, the truck haulage cost is found to account for up to 50% of the total direct cost of earthworks [3]. Given multiple sources and destinations (i.e. cut areas vs. fill areas) within the site, meticulous planning of site grading operations is crucial to materialize productivity and cost efficiency. To improve hauling efficiency and reduce hauling cost during earthmoving operations in the field, temporary haul road is generally built with lower quality of material (such as pit run or gravel) or directly laid on rough ground [4]. The majority of temporary haul road on a

⁎ Corresponding author. E-mail addresses: [email protected] (C. Yi), [email protected] (M. Lu).

http://dx.doi.org/10.1016/j.autcon.2016.08.022 0926-5805/© 2016 Elsevier B.V. All rights reserved.

rough grading site will be removed once the grading work is completed. In the current research, temporary haul road is simply classified into two categories: high grade vs. low grade. High grade haul road is gravel surfaced while low grade haul road is simply built on rough ground which may require frequent maintenance. According to the research conducted in [5], the construction of temporary haul road of higher grade can result in 50% increase in truck hauling speed while reducing road maintenance cost by 65%; thus, an effective temporary haul road layout design including a limited length of higher grade haul road potentially leads to total direct cost saving by as much as 15%. In regions where high open-pit mining activity takes place, regulations and design code for permanent haul road are well established with regards to alignment, surface and materials in order to ensure safe and efficient mining operations [6]. For instance, the road width needs be 3 to 4 times the width of the widest heavy hauler the road carries. To the contrary, design code for temporary haul road on rough-grading projects is generally non-existent. With limited time available, field managers need to decide whether to build higher grade temporary haul road and how to design temporary haul road layout largely based on experience instead of science; the overarching constraint is to ensure that the benefits resulting from shorter truck cycle time and lower equipment operating cost must exceed the expenses of building higher-grade temporary haul road [5].

C. Yi, M. Lu / Automation in Construction 71 (2016) 314–324

2. Literature review Optimization of the operational cost in mass earthworks has been focused on: (1) fleet balancing for achieving the maximum earthmoving productivity [7–11]; (2) investigating hauling strategies under geological constraints [12–15]; (3) determining earthmoving volumes from cut areas to fill areas on site [11,16–18]; and (4) identifying the proper sequence to execute earthmoving jobs [2,14,19,20]. It is worth mentioning that temporary haul road alignment design has not been explicitly defined as an optimization problem in previous construction engineering research. To a certain degree, the proposed temporary haul road layout design problem (THRLDP) is analogous to the classical transit route network design problem (TRNDP) which is formalized in traffic engineering. The focus of TRNDP is mainly on the optimization of multiple objectives characterizing the efficiency of transportation networks under operational and resource constraints [21]. Similarly, both TRNDP and THRLDP are intended to model entity flows (Passenger Flow vs. Earth Flow) by applying heuristics or analytical methods given a specific scope of the system [22]. In general, the traffic demand in THRLDP, which is defined in terms of entity flows of certain type and volume, can be predetermined once the haul road layout is designed. Based on (1) a balanced site grading design and (2) a temporary haul road layout design, earthmoving volumes and truck cycles from cut areas to fill areas can be specified by applying linear programming techniques. This differs from TRNDP which entails the prediction of traffic demand based on historical data [2,5]. Table 1 contrasts THRLDP against TRNDP in detail. 3. Problem statement At present, despite the fact that significant endeavors have been devoted to optimization of earthmoving operations, research is scant in regard to optimizing the layout design of temporary haul road in mass earthworks so as to minimize total earthmoving costs, while also catering to practical requirements in terms of connecting haul road on site with access road or permanent road offsite. Based on consultation with experienced industry professionals in Alberta who specialize in earthworks design and construction, it is identified that designing the temporary haul road layout in a timely and effective fashion presents a distinctive challenge to field engineers and managers. Conflicting constraints need to be balanced: on one hand, building the higher-grade haul road can increase the hauling efficiency, thus shortening project duration and reducing duration-dependent costs; on the other hand, the haul road construction cost would cause the contractor to considerably raise the field overhead. Meanwhile, the practical needs to maintain haul road continuity and accessibility make formulation of analytical solutions difficult. Temporary haul road design needs to avoid road sections of higher grade being sporadically scattered on site. Further, temporary haul road on site needs to be connected to the access road offsite via a single or multiple site entrances [5]. The present

315

research is intended to overcome the above identified challenges and address field application needs by developing an analytical optimization approach to design temporary haul road networks on rough grading projects, which minimizes total project cost while satisfying practical constraints of continuity and accessibility in temporary haul road design. In this research, a novel quantitative methodology is proposed for the identified problem by seamlessly integrating a cutting plane method into Mixed-Integer Linear Programming (MILP). In the following sections, mathematical formulation and optimization procedures are illustrated step by step with examples. First, relevant factors and assumptions considered in the problem definition are summarized as follows: 1. The site grading design, expressed as a cut and fill matrix, is already completed and remains unchanged. It is common practice that civil engineering professionals apply civil information modeling software (such as Autodesk Civil 3D) to achieve a cut-fill balanced grading design with zero or near zero net volume of earthworks. Unbalanced grading design scenarios having significant extra dump volume or borrow volume will be out of the scope of the current paper. 2. Choosing the truck hauling path on a particular job is dictated by the optimal cut-to-fill volume allocation derived from planning [5]. Except for trucks, equipment transit time from one cell to another in the field (including excavators, graders, and compactors) is ignored in the optimization formulation, because it is relatively trivial when compared to the truck hauling time. 3. To streamline the mathematical formulation, the average speed of a truck hauling empty on a particular haul road section would be the same as that of a truck hauling a full load. In other words, in the present research, the truck hauling speed only depends on the type of haul road being of higher grade (gravel surfaced) or of lower grade (rough ground.), irrespective of truck's load state. 4. “Truck loading and unloading plus waiting” time is generally considered in earthmoving cost estimating and fleet balancing studies; nonetheless, it is not considered in the formulation of the present problem, as the effect of such time is not as relevant as the truck haul time in the context of temporary haul road layout design. 4. Methodology overview The proposed methodology is shown in a flowchart in Fig. 1. Based on the site grading design, the field is divided into cells for determining the cut or fill volume of each particular cell in a unified volume measure (e.g. bank cubic meter), resulting in the earthwork volume matrix. Next, empirical parameters and practical constraints are evaluated, with basic cost data estimated as model inputs. Then, the optimum haul road layout design can be achieved by minimizing the cost function defined in the Mixed-Integer Programming formulation. A cutting plane method is integrated into MILP (C-MILP in short) in order to add “connectivity” constraints to haul road layout design

Table 1 Contrasting main features of two problems: TRNDP vs. THRLDP. Category

TRNDP

THRLDP

Flow type Objective

Passenger flow Multiple objectives (e.g. Minimizing operation cost, maximizing service efficiency, maximizing utilization rate and etc.) Route layouts, frequencies, fleet availability, demand characteristics [23] Operating environment (e.g. Network structure), operational strategies and rules, and available resources [23] Multi-source-multi-destination Estimated from historical statistics Heuristic methods or analytical methods

Earth flow Single objective (i.e. Minimizing operation cost)

Variables Constraints Flow type Flow quantity Optimization methods

Haul road layouts, traffic demand governed by grading design and earthworks volume allocation Operating environment (e.g. Network structure), haul road continuity & accessibility, operational strategies and rules, natural blocks, and available resources Multi-source-multi-destination Determined by applying mathematical methods [24] Heuristic methods or analytical methods

316

C. Yi, M. Lu / Automation in Construction 71 (2016) 314–324

optimization. As a result, the proposed optimization model is capable of generating the solution in the form of a feasible temporary haul road design in a two dimensional site layout. Validating the proposed methodology entails two methods: (1) based on computer-simulated test-bed cases, the optimum solution obtained is benchmarked against related methods published in literature in terms of computing efficiency and optimization performances; (2) based on real-world project information, the optimum solution obtained is confirmed by the project manager who was involved in the actual design.

5. Input data The input data as needed for implementing the proposed methodology include cut and fill volume data of each cell in the site area, average haul speeds of trucks on haul road of different grades, and project-specific empirical parameters or field data, as listed below: • Construction and removal cost of the gravel-surfaced road ($/km); • Maintenance cost ($/km) for gravel-surfaced road and for rough ground haul road, respectively; • Mean truck-haul speed on “gravel-surfaced” road; • Mean truck-haul speed on “rough ground” road; • Truck volume capacity and truck quantity; • Hourly rate of equipment; • Hourly rate of labor; • Hourly rate of project overhead; • Working efficiency factor;

• Nearby permanent road and access road with specified connection locations to the temporary haul road in the field. 6. Grid network structure Grid-based systems have been commonly applied to symbolize the road network structure [23]. The site is generally divided into square cells with cell width being 150 m or 200 m for the purpose of temporary haul road design on a rough grading site [24,25]. In the current research, the grid network structure for haul road design consists of 150 m × 150 m square cells. The conceptual model of a possible haul road layout design is demonstrated in Fig. 2. For each cell, the centroid is simplified to be the cell's geometric center; thus, the potential road layout design can be denoted by road links, each connecting the centroids of two adjacent cells with straight-line sections. The road type of each link can be distinguished by a dot line for “rough ground” by default and by a solid line for “gravel surfaced”, as shown in Fig. 2. It should be noted that the “road links” are segments between the centroids of any two adjacent cells, instead of between any two cell centroids; two cells are deemed adjacent only if they are “immediately adjacent” or “diagonally adjacent”. Based on Graph Theory, an undirected graph is written as G= (V,E), meaning that G consists of node-set V and edge-set E. We use v to represent a node (centroid of a cell) and (i,j) to represent an edge (road link) in the haul road network. Note ∀v ∈ V, and all ∀ (i,j)∈E where i and j are end nodes of edge (i,j). The type of haul road is denoted by a Boolean parameter x(i, j), which equals to “0” given “rough ground” haul road; equals to “1” given “gravel surfaced” haul road (e.g. in Fig. 2. x(1,2) =0 ; x(6,7) = 1).

Fig. 1. Flowchart to overview the proposed optimization methodology.

C. Yi, M. Lu / Automation in Construction 71 (2016) 314–324

317

The second constraint is to restrict f(i,j) to a real number, as shown in Eq. (4). f ði; jÞ ∈R

ð4Þ

All f(i,j) is non-negative according to its definition, as constrained by Eq. (5). f ði; jÞ ≥0

Given the case of rough-grading site design, optimal cut-to-fill assignments can be obtained by the optimization methodology and the results are shown in Fig. 3 (“+” means cut while “−” means fill). In short, optimization of traffic volume on each road link provides the basis for further optimizing temporary haul road design.

Fig. 2. Site partitioning strategy and road links to represent haul road design.

7. Cut to fill volume assignment Determining the traffic volume on each road link by a scientific method is the prerequisite to haul road design. Note traffic volume assignment is generally defined as a cut–fill cell pairing, which represents a specified volume of material is excavated from a specified cut location and hauled to a specified fill location. The optimum traffic volume on each road link can be determined by applying linear programming techniques [5,17]. In order to assign the optimal traffic volume on each road link, the truck haul time is considered as the key decision variable, instead of the haul distance, because truck hauling speed and road maintenance cost differ markedly on haul road of different grades. The haul time is calculated simply as t(i, j) = s(i, j)/v(i, j), where t(i, j) is the truck haul time per truck load on edge (i, j); s(i, j) is the distance of edge (i, j); and v(i,j) is haul speed of a truck on edge (i, j). Then, the average haul time in the field is defined in Eq. (1): 2∙ ∑ t avg ¼

h

ði; jÞ∈A

 i f ði; jÞ ∙ t ði; jÞ −xði; jÞ ∙Δt ði; jÞ   ∑ f ði; jÞ

ð5Þ

ð1Þ

8. Cost function With more high-grade road links, the average unit haul time is reduced at the expense of potentially increasing the cost of building and maintaining the temporary haul road network on site. Hence, the optimization of temporary haul road design becomes a tradeoff problem: how to achieve a balance between project-duration-dependent cost and road construction plus maintenance cost. The objective is to minimize the total cost by combining the two costs and identifying the optimum equilibrium between them. The total cost is given as Eq. (6). C total ¼ C duration þ C road

ð6Þ

The project duration dependent cost is estimated by Eq. (7) r   C duration ¼ ∑ α p ∙T

ð7Þ

p¼1

ði; jÞ∈A

where tavg is the averaged round-trip haul time per truck load, assuming that truck returns on the same path in each cycle, f(i,j) is the allocated earthwork quantity on an edge (i, j), and ∑ ðf ði; jÞ Þ equals to the ði; jÞ∈A

total earthwork quantity in cubic meters. Δt(i,j) is the time saving per truck load on edge (i, j) if a truck travels on “gravel-surfaced” road in contrast with “rough-ground” road. x(i,j) is a Boolean input parameter. Note as tavg is the averaged round-trip haul time per truck load; assuming truck returns on the same path in each cycle, ‘2’ is applied to account for the round trip effect. Next, the optimal traffic volume on each road link can be determined by solving for the decision variable f(i,j) in the linear programming based mathematical model: 2∙ ∑ Min



h

ði; jÞ∈A

 i f ði; jÞ ∙ t ði; jÞ −xði; jÞ ∙Δt ði; jÞ   : ∑ f ði; jÞ

ð2Þ

ði; jÞ∈A

Z = objective function, which equals to the average haul time as defined in Eq. (1). The optimization is to obtain the allocated earthwork quantity on all edges resulting in the minimum averaged haul time [5]. f+(v) is the earth volume flowing into node v; f−(v) is the earth volume flowing out of node v; and Dv is the earth demand (cut or fill) as per cut and fill design of node. Thus, Eq. (3) is given to constrain the volume flow conforming to the cut and fill design. þ



f ðvÞ−f ðvÞ ¼ Dv

ð3Þ

where αp signifies duration related rates such as crew hourly or daily rate, equipment hourly or daily rate and hourly or daily rates of miscellaneous overhead; p is a number series from 1 to r for denoting all the relevant duration-dependent rates. T is the total project duration estimated by Eq. (8): ∑





ði; jÞ∈A

f ði; jÞ

 ∙t avg

f ∙n∙c

ð8Þ

where ∑ ðf ði; jÞ Þ is the total earthwork quantity in cubic meters; n is ði; jÞ∈A

the number of trucks (fleet size), c is the volume capacity of one truck in cubic meters, and f is the operations efficiency factor (commonly taken as 45-min h or 0.75 in construction estimating). Eq. (8) is given to determine the total duration, which allows for multiple trucks to work concurrently. Combined with Eqs. (1), (8) can be written as Eq. (9) 



 i f ði; jÞ ∙ t ði; jÞ −xði; jÞ ∙Δt ði; jÞ ði; jÞ∈A ði; jÞ∈A   ∙ T¼ f ∙n∙c ∑ f ði; jÞ ði; jÞ∈A h  i 2∙ ∑ f ði; jÞ ∙ t ði; jÞ −xði; jÞ ∙Δt ði; jÞ ∑

¼

f ði; jÞ

2∙ ∑

ði; jÞ∈A

f ∙n∙c

h

ð9Þ

Then, factorized by Eq. (9) and combined with Eq. (7), the project duration dependent cost can be further simplified to Eq. (10) where α

318

C. Yi, M. Lu / Automation in Construction 71 (2016) 314–324

(1) Cut and fill design

(2) Optimized earthwork volume allocation

Fig. 3. Illustration of grading design model and earthwork volume allocation. (1) Cut and fill design. (2) Optimized earthwork volume allocation.

is the combined duration-dependent cost parameter. C duration ¼ α∙ ∑

h

ði; jÞ∈A

 i f ði; jÞ ∙ t ði; jÞ −xði; jÞ ∙Δt ði; jÞ

ð10Þ

Where α¼

  r 2∙ ∑p¼1 α p

ð11Þ

f ∙n∙c

The haul road construction plus maintenance cost is defined by Eq. (12), C road ¼ C construction þ C maintenance:

ð12Þ

The haul road construction cost is proportional to the length of highgrade road as given in Eq. (13), where β is the unit construction plus removal cost of the gravel-surfaced road. Note the construction and removal cost of the rough-ground road is insignificant hence ignored.   C construction ¼ β∙ ∑ sði; jÞ ∙xði; jÞ

ð13Þ

ði; jÞ∈A

Despite much lower construction cost, rough-ground road is much more expensive to maintain due to safety-related risks and high wear and tear on tires and trucks. Note in a previous study the maintenance cost was simply defined as a function of the proportion of the length of gravel-surfaced haul road over the total length of haul road in the field [5]. As the maintenance cost varies with truck traffic volumes, the haul road to carry lesser traffic volumes generally incurs lower maintenance cost. As the traffic volume on each road link can be fixed by cutfill earthwork volume allocation, the maintenance cost of haul road is given as in Eq. (14) factoring in the traffic volume and per km maintenance cost for rough-ground road and gravel-surfaced road, respectively. C maintenance ¼ ∑

ði; jÞ∈A

    ð1Þ ð2 Þ σ 1 ∙sði; jÞ ∙f ði; jÞ þ ∑ σ 2 ∙sði; jÞ ∙f ði; jÞ ði; jÞ∈A

ð14Þ

Where σ1 ¼

σ2 ¼

C avg−rough Q1

ð15Þ

C avg−temp Q2

ð16Þ

of rough ground haul road in one day. Similarly, σ2 is the per km haul road maintenance cost if the truck hauls on gravel-surfaced haul road, which depends on Cavg−temp ( the maintenance cost per km gravel-surfaced haul road in one day considering average traffic volume and typical haulers) and Q2 ( the average traffic volume in truckloads passing the 1 km long section of gravel-surfaced haul road in one day). (1) (2) (1) f(i,j) and f(i,j) result from cut to fill assignment optimization: f(i,j) is the allocated earthwork quantity in connection with rough ground haul (2) road, while f(i,j) is the allocated earthwork quantity associated with sections of gravel-surfaced haul road. Note the existence of two extreme cases: (1) if there is no temporary (2) gravel-surfaced haul road to be built, then f(i,j) = 0; (2) if trucks haul on (1) gravel-surfaced haul road across the entire site, then f(i,j) = 0. Eq. (10) can be further transformed by denoting “traffic scenario” on rough-ground haul road and gravel-surfaced haul road separately, as in Eq. (17): C duration ¼ α 

X ð1Þ ð1Þ ð2Þ ð2Þ ði; jÞ∈Að f ði; jÞ  t ði; jÞ þ f ði; jÞ  t ði; jÞ

ð17Þ

Combining Eqs. (13), (14), and (17), the cost function is finally given as Eq. (18),    ð1Þ ð1Þ ð2Þ ð2Þ f ði; jÞ ∙t ði; jÞ þ f ði; jÞ ∙t ði; jÞ þ β∙ ∑ sði; jÞ ∙xði; jÞ ði; jÞ∈A ði; jÞ∈A     ð1Þ ð2Þ þ ∑ σ 1 ∙sði; jÞ ∙f ði; jÞ þ ∑ σ 2 ∙sði; jÞ ∙ f ði; jÞ

C total ¼ α∙ ∑



ði; jÞ∈A

ði; jÞ∈A

ð18Þ

9. Mixed-Integer Linear Programming Model With the development of mathematical programming theory and the advance of off-the-shelf computing software, complicated problems can be solved analytically based on mathematical model formulations using integer variables, nonlinear constraints, asymmetry structured networks etc. [26–28]. Advanced analytical methods have also been applied in transit route network design problems [28,29]. The temporary haul road layout design problem is defined in this research and formulated as a Mixed-Integer Linear Programming model by following the similar strategy applied in [30].

9.1. Decision variables

σ1 is the per km haul road maintenance cost per truckload if the truck hauls on rough ground. σ1 depends on Cavg − rough, which is the maintenance cost per km rough ground road in one day considering average traffic volume by typical haulers (e.g. CAT 777), and Q1, which is the average traffic volume in truckloads passing the 1 km long section

To extend the optimal cut-to-fill assignment problem, a new combinatorial optimization problem is formulated factoring the decision variable as of whether to construct temporary haul road of high grade, as given in Eq. (19).  xði; jÞ ¼

0 1

“rough ground”link “gravel surfaced”link

ð19Þ

C. Yi, M. Lu / Automation in Construction 71 (2016) 314–324 (1) (2) Two additional variables are declared as f(i,j) and f(i,j) . The subscript “(i,j)” indicates the edges (i.e. road links) on the site layout graph. The superscripts “(1)” and “(2)” indicate traffic volume being allocated on either rough ground haul road or gravel-surfaced haul road, respectively.

The proposed model in this study is intended to directly minimize the total cost based on the proposed cost function in Eq. (18). Then the objective function can be given in Eq. (20).    ð1Þ ð1Þ ð2Þ ð2Þ f ði; jÞ ∙t ði; jÞ þ f ði; jÞ ∙t ði; jÞ þ β∙ ∑ sði; jÞ ∙xði; jÞ ði; jÞ∈A ði; jÞ∈A     ð1Þ ð2Þ þ ∑ σ 1 ∙ sði; jÞ ∙f ði; jÞ þ ∑ σ 2 ∙sði; jÞ ∙f ði; jÞ

Min Z ¼ α∙ ∑



ði; jÞ∈A

ði; jÞ∈A

ð1Þ

ð1Þ

ð2Þ

ð20Þ

ð2Þ

Z = objective function; α∙ ∑ ðf ði; jÞ ∙t ði; jÞ þ f ði; jÞ ∙t ði; jÞ Þ is to summaði; jÞ∈A

rize the project-time-dependent cost including labor, equipment and miscellaneous overhead expenses; note, α is related with haul efficiency, unit rates, number of trucks in the fleet, and truck capacity [referring to Eq. (7)]. Similarly, β∙ ∑ ðsði; jÞ ∙xði; jÞ Þ is to represent the high-grade ði; jÞ∈A

haul road construction cost, depending on the per km construction cost β and the distance of haul road of high grade to be built (gravel-surð1Þ

ð2Þ

faced). In addition, ∑ ðσ 1 ∙sði; jÞ ∙f ði; jÞ Þ þ ∑ ðσ 2 ∙sði; jÞ ∙ f ði; jÞ Þ is to denote ði; jÞ∈A

ði; jÞ∈A

the haul road maintenance cost, which depends on per km maintenance cost parameters (σ1 or σ2) the distance of road sections of particular type and the earth volume flow passing each road section. 9.3. Constraints As explained in the mathematical model for optimal traffic volume assignment, constraints are imposed in order to restrict the earth volume flows on the haul road network model conforming to the cut and fill design. þ

negative. ð1Þ

ð24Þ

ð2Þ

ð25Þ

f ði; jÞ ≥0 f ði; jÞ ≥0

9.2. Objective function



f ðvÞ−f ðvÞ ¼ Dv

ð21Þ

The traffic volumes being allocated on rough ground haul road or on gravel-surfaced haul road are real numbers, as demonstrated in Eqs. (22) and (23). ð1Þ

ð22Þ

ð2Þ

ð23Þ

f ði; jÞ ∈R f ði; jÞ ∈R

Eqs. (24) and (25) maintain that the traffic volume being allocated on rough ground haul road or on gravel-surfaced haul road is non

(1) Original label

319

9.4. Cutting plane method The cutting plane method is an integer programming technique which is intended for iteratively refining the feasible solution set by introducing linear inequalities (termed cuts) [24]. In this temporary haul road layout design problem, the solution obtained from the above-formulated MILP-based mathematical model does not address constraints of road continuity and field accessibility. In the current research, the cutting plane method is integrated into MILP so as to impose appropriate cutting-planes to the MILP model. Those cutting planes are defined as valid inequalities derived based on the concept of connected components. In graph theory, a connected component is a subgraph belonging to the main graph in which any two nodes are connected by edges; note a solo node with no connected edges can also be characterized as a special case of connected component, which is shown in Fig. 4(1). In the case of temporary haul road design problem, a node corresponds to the center point of each cell while an edge represents a potential road link connecting two cells. A connected component can be defined as node set, in which any node is connected to at least another node in the same set to form an edge. If the node set has only one node, the solo node does not connect with any other node in the system. In this study, the connected component is defined as per Eq. (26), where all v in Vm are initially assigned with a unique label Lm , (v) which belongs to a label set Lm; m is the label set number (e.g. the initially assigned label is expressed as L1,(v)) ;and X(E) is the edge set containing all the high-grade road links; Then, connected components are further defined as sub-sets of nodes where nodes in each set has the same label, all (i, j) ∈ X(E). Lm;ðiÞ ¼ Lm;ð jÞ

ð26Þ

By applying Eq. (26), Fig. 4(1) is further updated by changing node labels, shown in Fig. 4(2). Then, the connected component(s) with the largest number of nodes [i.e. the longest connected component(s)] is chosen and the particular node set containing these nodes is denoted as the set S. Once S is selected, valid inequalities, as demonstrated in Eq. (27), are added to the mathematical model as part of the constraints, where S ∈ x(V) . x(V) is the node set of high-grade road links; the edge cut of a node set include those edges that do not belong to the current node set but can be immediately linked with a node in the current set. For example, the edge cut

(2) Updated label

Fig. 4. Connected components. (1) Original label. (2) Updated label.

320

C. Yi, M. Lu / Automation in Construction 71 (2016) 314–324

9.5. Demonstration case As the main purpose is to demonstrate how the cutting plane method is applied, parameters α and β in the objective function are arbitrarily set in the case study. Each node is initially assigned a label from 1 to 12, defined as the label set L1 ={1, 2, …12}. Solving the proposed MILP-based mathematical model (identified as Model No.1), the haul road layout design is generated (Fig. 6(1)). Nevertheless, the high-grade haul road layout is neither accessible (at entrance cell 1) nor continuous (road link 3–7 is not connected to the other road links). Following Eq. (26), connected components can be identified. The process of defining connected components is shown in Eq. (28) L1;ð9Þ ¼ L1;ð5Þ ; L1;ð10Þ ¼ L1;ð5Þ ; L1;ð11Þ ¼ L1;ð5Þ ; L1;ð7Þ ¼ L1;ð3Þ

ð28Þ

Then the node set with L1 , (v) = L1 , (5) has the largest number of nodes and is chosen to be the connected component. Along with the accessibility requirement (trucks access the site from permanent road via entrance at Cell 1), two valid inequalities, given in Eqs. (29) and (30), are added to the original mathematical model as constraints. The updated mathematical model is identified as Model No.2. x1 ðδð1ÞÞ ¼ x1 ð1−2; 1−5; 1−6Þ≥1

ð29Þ

x1 ðδðf5; 9; 10; 11gÞÞ ¼ x1 ð5−1; 5−2; 5−6; 9−6; 10−6; 11−6; 11−7; 11−8; 11−12Þ≥1

Fig. 5. Flowchart of integrating a cutting plane method into MILP.

ð30Þ of connected component 2 in Fig. 4(1) includes eleven edges, namely: (3, 2), (3, 4), (3, 6), (3, 8), (7, 2), (7, 4) (7, 6) (7, 8) (7,10) (7, 11) (7, 12). δ(S) is used to denote the edge cut set of S; xm(δ(S)) is to add up the road type variable (0/1) associated with each edge in δ(S) . Thus , xm(δ(S)) ≥ 1 means that at least one edge in the edge cut of S is included in S in the next generation of the solution. With valid inequalities added, the updated haul road layout design can be obtained by solving the mathematical model. xm ðδðSÞÞ≥1

ð27Þ

It is important to mention that site accessibility can be guaranteed by adding a valid inequality in connection with a solo node set, defined as the entrance set. The node itself is a connected component, which is also referred to as the entrance component. By repeating the above process, the final solution can be obtained, which is illustrated in Fig. 5. In the following section an illustration case is presented based on the example given in Fig. 3 in order to elaborate the complete formulation.

1) First run

Where x1(δ(1)) ≥ 1 means the objective function is subject to the constraint that at least one edge in the edge cut of node set {1} should be included in the next generation of solution forming part of the high-grade haul road network. Similarly, x1(δ({5, 9, 10, 11})) ≥ 1 means the objective function is subject to the constraint that at least one edge in the edge cut of node set {5, 9, 10, 11} should be present in the next generation of the solution. Then, by solving Model No.2, the updated high-grade haul road network is generated [Fig. 6(2)]. The solution is accessible but still not fully connected (road link 1-2 is detached from the others). Each node is then assigned a second label from 1 to 12, defined as a set L2 = {1, 2, …12}. Following Eq. (26), new connected components can be identified as given in Eq. (31) L2;ð9Þ ¼ L2;ð5Þ ; L2;ð10Þ ¼ L2;ð5Þ ; L2;ð7Þ ¼ L2;ð5Þ ; L2;ð2Þ ¼ L2;ð1Þ

ð31Þ

Then node set with L2,(v) = L2 ,(5) has the largest number of nodes and hence is chosen to be the connected component. One extra valid inequality, as given in Eq. (32), is added to Model No.2 as constraint according to Eq. (26). Then the mathematical model is further updated

2) Second run

3) Third run

Fig. 6. C-MILP optimization process and solution illustration. 1) First run 2) Second run 3) Third run.

C. Yi, M. Lu / Automation in Construction 71 (2016) 314–324

321

type: 64-bit Operating System. Python Version 3.4 [31] with a Python module integrated with CPLEX Version 12.61 [32] was utilized to code the algorithms in the current case study. MILP MCGA A total of ten test sets were randomly generated for assessing perforTotal CPU Total Site layout Total earthwork CPU mances of the two algorithms and contrasting them. A summary of test time (s) cost ($) time (s) cost ($) Test no. (m× n) volumes (m3) settings and results is given in Table 2. The test sets each represent a 1 4×8 1,261,800 0.61 237,720 28.87 264,879 rough grading project of particular size (i.e. the project size increases 2 4×9 1,362,800 1.12 299,429 35.46 325,264 with a larger number of cells), with randomly-simulated cut and fill vol3 4 × 10 1,946,900 1.58 264,234 50.87 285,673 4 5×9 2,160,500 6.75 304,785 103.85 352,687 umes. For each test case, the basic cell in the grid has a constant size, 5 5 × 10 2,188,900 7.57 350,798 285.36 387,452 namely: 150 m by 150 m; while, the m × n site layout means that the 6 5 × 11 2,787,900 8.33 394,662 361.25 409,856 site has m rows and n columns of cells (e.g. the site in Fig. 3 is a 3 × 4 lay7 6 × 10 2,697,300 11.22 426,944 365.87 449,562 out); the total earthwork volume represents the sum of cut volumes, 8 7 × 11 2,672,900 19.11 419,685 821.63 438,561 9 8 × 12 3,178,800 345.04 435,341 – – which is equal to the sum of fill volumes in each test case. 10 9 × 15 3,672,000 977.22 525,558 – – Performance assessment was conducted by independently applying the two algorithms to solve each “test set” problem. Note, in test set 9 and 10, MCGA did not converge after excessive computation time, thus no results were reported. Based on the test results from the first to Model No.3. 8 sets, MCGA spent much more computing time than the proposed x2 ðδðf5; 7; 9; 10gÞÞ method. In addition, the total cost resulting from MCGA is on average ¼ x2 ð5−1; 5−2; 5−6; 9−6; 10−6; 7−2; 7−3; 7−4; 7−6; 7−8; 7−11Þ≥1 11% higher than that from MILP being proposed. ð32Þ Table 2 Test problem details and result comparison.

The node set with L2,(v) = L2 ,(5) is chosen to be the connected component, which means the updated solution must have at least one edge in the edge cut of the node set with L2 ,(v) = L2,(5) . By solving Model No.3, a new high-grade haul road layout is generated, as shown in Fig. 6(3). The layout is both accessible and interconnected. The optimization process halts and the latest solution is taken as the final optimum design. 10. Validation 10.1. Numerical simulation In previous research, Liu [4] utilized an advanced version of GA named Multi-Generation Compete Genetic Algorithm (MCGA) to design the temporary haul road layout as part of site grading design. The MCGA algorithm allows earlier generations of solutions to participate in the current iteration of competition. It is emphasized that the enhanced version of GA still takes a non-deterministic amount of time in converging to solutions while still having good chances of converging to local minima. In order to prove a new algorithm is superior to a proven one, common practice is to apply numeric investigation based on the same problem [19,20]. The results from the newly proposed method and the GA method [4] based on different cases are contrasted in terms of (1) the computing time required to solve haul road design problems of various sizes, and (2) the derived optimum solution itself. The computer platform used in the cross validation had the following configurations: Intel Core i7-5500 U CPU, up to 3.00 GHz; Memory: 8.0 GB; System

10.2. Practical case A practical rough grading project was utilized for further validating the proposed method. The rough grading project was the preliminary work package of a camp site construction in Fort McMurray, AB. The field was around 120 ha, divided into 48 cells each being 150 m by 150 m. The project had a total amount of 335,600 m3 of earth to be handled from cut and fill. The cut or fill volume of each cell along with the cell identification number is shown in Fig. 7. The rough ground (0–2 m depth top soil) in the field was mainly composed of sand and glacial till. The “rough ground” haul road in the current case would lead to slower hauling speed, higher maintenance cost and potential safety hazard. To ensure operations safety and control project cost, the contract stipulated that temporary haul road be constructed at the cost of the contractor. The parameters relevant to the MILP model formulation were evaluated based on field data and empirical data with the help of the project manager of the contractor company. The construction plus removal cost of gravel-surfaced temporary haul road was estimated at $17,500/km; given average traffic volume with typical haulers for rough ground haul road, the maintenance cost in one day was estimated at $600/km given average traffic volume of 150 truckloads; the maintenance cost in one day was estimated at $150/km with average traffic volume of 180 truckloads per day. The mean truck-haul speed on gravel-surfaced haul road was 36 km/h; the mean haul speed on rough ground haul road was 24 km/h; the truck volume capacity was 30 m3 (CAT 777); the hourly rate for all the equipment was $3000/h; the hourly rate for all the personnel was $2000/h; the hourly rate for field overhead was

Fig. 7. Cut and fill volumes of each cell on practical case study.

322

C. Yi, M. Lu / Automation in Construction 71 (2016) 314–324

(1) Empirical design

(2) Optimal design

Fig. 8. Practical design vs. optimum design of temporary haul roads for case study (1) Empirical design. (2) Optimal design.

$200/h; the operation efficiency factor was 0.75; and 8 trucks of the same type made up the fleet. The temporary haul road layout which was designed by the field engineer is shown in Fig. 8(1). Based on the empirical design, the total earthmoving cost was budgeted at approximately $ 400,000 CAD including haul road construction plus maintenance. Two shortfalls were observed in the empirical design, namely: (1) no gravel-surfaced haul road section connects to the permanent road at the entrance; and (2) scientific proof of the cost-effectiveness of this haul road design was absent. The C-MILP model was coded in Python Version 3.4 following the proposed optimization procedure (Fig. 1). By executing functions in CPLEX Version 12.61, an optimum solution for the case study was produced, as shown in Fig. 8(2). Throughout the optimization process, a total of 63 solutions (62 intermediate results plus 1 the final result) were generated automatically by computer. The program terminated when both accessibility and continuity of temporary haul road design were achieved. Design layouts for all 63 solutions are shown in Table B.1 in Appendix B. The resulting temporary haul road design is smoothly connected to the entrance; on the other hand, the as-designed alignment of temporary haul road in the field consists of two parallel straightline sections, interconnected by a road link in the middle. The total cost was determined to be $ 361,854 CAD. The optimum haul road layout would potentially reduce the total budget $ 400,000 CAD by nearly 10%. 11. Conclusion This research has introduced an analytical method for optimizing temporary haul road layout design in planning rough-grading earthworks projects, aimed to minimize the earthmoving cost while satisfying site accessibility and haul road continuity requirements. The problem formulation starts by applying Mixed-Integer Linear Programming (MILP), which is supplemented by a cutting plane method which imposes additional constraints allowing for field accessibility and haul road continuity. Performances of the analytical method being proposed were cross-checked by comparing against a GA-based methodology developed in previous research. Results from numerical simulations showed that the proposed method outperformed GA in terms of computing efficiency and reliability of finding optimum solutions. Then, based on a real-world rough grading project, the effectiveness of the proposed methodology was further validated by improving the high grade temporary haul road layout design. The optimum design resulting from C-MILP can save approximately the total budget by 10% while maintaining accessibility and continuity constraints in the temporary haul road layout design. By applying the innovative cutting plane method based MILP in this research, the empirical design problem of temporary haul road design for rough grading projects in construction engineering has been, to

some extent, defined as an analytical problem for optimization. Computer programming automation has also been materialized to generate the optimal temporary haul road layout design within a deterministic time frame. With regards to current practice improvement, this research caters for the immediate needs of earthworks contractors for enhancing the temporary haul road design on rough grading projects, potentially leading to field productivity improvement and better budget control. Furthermore, aided by analytical methodology and computer automation, engineers who may have limited field experience in temporary road network design would be more confident and competent to cope with the real world challenges. In order to deliver field-ready solutions to an empirical temporary design problem in a fully analytical and automated fashion, further improvements of the present research will be required in the future. The current problem definition can be further extended in order to cope with the more complicated scenario of unbalanced cut and fill design in the real world. Full integration and automation of the proposed methodology would also be worthy of pursuit by linking the site grading and temporary haul road design system with databases which track construction cost performances; as such all the empirical cost data used as input parameters in the current model will be replaced with historical data that are more relevant and accurate. In addition, the temporary haul road design in this study is performed on a relatively flat grading site without considering road blocks such as waterways, deep valleys or high hills, which keep heavy trucks from passing. Nonetheless, such road blocks in a real-world construction site add further constraints and would compound the optimization formulation of the proposed problem. Acknowledgements Authors are grateful to Sam Johnson, Project Manager, Graham Construction & Engineering Inc., for sharing field experiences in temporary haul road design and providing valuable insight, information, and data in this research. Appendix A. Notation x(i,j) G V E v (i, j) t(i,j) s(i,j) v(i,j) f(i,j) tavg

Boolean variable which is determined by whether to build higher-grade temporary haul road on edge (i, j); undirected graph; node set; edge set; node; directed edge; truck haul time per truck load on edge (i, j); distance of edge (i,j); haul speed of a truck on edge (i, j); allocated earthwork quantity on an edge (i, j) averaged round-trip haul time per truck load

C. Yi, M. Lu / Automation in Construction 71 (2016) 314–324

323

Table B.1 Interim layouts in optimization process. Run number

Intermediate layout

Run number

Intermediate layout

Run number

1

22

43

2

23

44

3

24

45

4

25

46

5

26

47

6

27

48

7

28

49

8

29

50

9

30

51

10

31

52

11

32

53

12

33

54

13

34

55

14

35

56

15

36

57

16

37

58

17

38

59

18

39

60

19

40

61

20

41

62

21

42

63

Intermediate layout

324

C. Yi, M. Lu / Automation in Construction 71 (2016) 314–324

Δt(i,j)

difference in haul time per truck load on edge (i,j) if a truck travels on “gravel-surfaced” road, instead of “rough-ground” road; f+(v) earth volume flowing into node v; f−(v) earth volume flowing out of node v; Dv earth demand (cut or fill) as per cut and fill design of node v; Ctotal total costs; Cduration project duration related costs; Croad haul road construction and maintenance costs; Cconstruction construction and removal costs of the gravel-surfaced road; Cmaintenance haul road maintenance costs; p a number series from 1 to r for denoting all the relevant duration-dependent rates; αp different duration related rates; T total haul duration; n truck number; c volume capacity of one truck in cubic meters; f operations efficiency factor; α combined duration-dependent cost parameter; β unit construction and removal cost of the gravel-surfaced road; σ1 unit maintenance cost if trucks haul on rough ground; σ2 unit maintenance cost if trucks haul on gravel-surfaced haul road; Q1 average traffic volume measured in truckloads of typical haulers using rough ground in one day; Q2 average traffic volume measured in truckloads using gravelsurfaced haul road in one day; Cavg−temp maintenance cost per km given average traffic volume with typical haulers using gravel-surfaced haul road in one day; Cavg−rough maintenance cost per km given the average traffic volume with typical haulers using rough ground in one day; (1) f(i,j) allocated earthwork quantity on rough ground haul road; (2) f(i,j) allocated earthwork quantity on gravel-surfaced haul road; m label set number; L label set; x(E) edge set of higher-grade (gravel-surfaced) road; x(V) node set of higher-grade (gravel-surfaced) road links; S set(s) of connected component(s) with the largest number of same label (i.e. longest connected component(s)); δ(S) edge cut set of S Appendix B. Practical case study optimization process and specifications

References [1] B.D.A. Prata, A stochastic colored petri net model to allocate equipments for earthmoving operations, J. Inf. Technol. Constr. 13 (2008) 476–490 http://itcon.org/ data/works/att/2008_28.content.05477.pdf. [2] W.L. Hare, V.R. Koch, Y. Lucet, Models and algorithms to improve earthwork operations in road design using mixed integer linear programming, Eur. J. Oper. Res. 215 (2011) 470–480, http://dx.doi.org/10.1016/j.ejor.2011.06.011. [3] R. Thompson, A. Visser, Mine haul road maintenance management systems, J. South. Afr. Inst. Min. Metall. (2003) 303–312. [4] C. Liu, Optimization of temporary haul road design and earthmoving job planning based on site rough-grading designMaster's Thesis University of Alberta, Canada, 2014http://dx.doi.org/10.1007/s13398-014-0173-7.2. [5] C. Liu, M. Lu, Optimizing earthmoving job planning based on evaluation of temporary haul road networks design for mass earthworks projects, J. Constr. Eng. Manag. (2015) 1–14, http://dx.doi.org/10.1061/(ASCE)CO.1943-7862.0000940. [6] D. Tannant, B. Regensburg, Guidelines for Mine Haul Road Design, 2001 1–10 https://circle.ubc.ca/bitstream/id/90332/haul_road_design_.

[7] F. Cheng, Y. Wang, X. Ling, Multi-Objective Dynamic Simulation-Optimization for Equipment Allocation of Earthmoving Operations, Construction Research Congress, 2010 328–338, http://dx.doi.org/10.1061/41109(373)33. [8] F.Y.Y. Ling, Y. Ke, M.M. Kumaraswamy, M. Asce, S. Wang, Key relational contracting practices affecting performance of public construction projects in China, J. Constr. Eng. Manag. (2013) 1–12, http://dx.doi.org/10.1061/(ASCE)CO.1943-7862.0000781. [9] M. Marzouk, O. Moselhi, Multiobjective optimization of earthmoving operations, J. Constr. Eng. Manag. 130 (2004) 105–113, http://dx.doi.org/10.1061/(ASCE)CO. 1943-7862.0000781. [10] O. Moselhi, A. Alshibani, Crew optimization in planning and control of earthmoving operations using spatial technologies, J. Inf. Technol. Constr. 12 (2007) 121–137 http://www.itcon.org/2007/7. [11] O. Moselhi, A. Alshibani, Optimization of earthmoving operations in heavy civil engineering projects, J. Constr. Eng. Manag. 135 (2009) 948–954, http://dx.doi.org/10. 1061/(ASCE)0733-9364(2009)135:10(948). [12] H. Gwak, C. Yi, D. Lee, Stochastic optimal haul-route searching method based on genetic algorithm and simulation, J. Comput. Civ. Eng. 04015029 (2015) 1–12, http:// dx.doi.org/10.1061/(ASCE)CP.1943-5487.0000496. [13] S. Kang, J. Seo, GIS method for haul road layout planning in large earthmoving projects: framework and case study, J. Constr. Eng. Manag. 139 (2013) 424, http://dx. doi.org/10.1061/(ASCE)CO.1943-7862.0000561. [14] G. L'Heureux, M. Gamache, F. Soumis, Mixed integer programming model for short term planning in open-pit mines, Min. Technol. 122 (2013) 101–109, http://dx.doi. org/10.1179/1743286313Y.0000000037. [15] J. Son, K. Mattila, D. Myers, Determination of haul distance and direction in mass excavation, J. Constr. Eng. Manag. 131 (2005) 302–309, http://dx.doi.org/10.1061/ (ASCE)0733-9364(2005)131:3(302. [16] R.X. de Lima, E.F.N. Júnior, B. de Athayde Prata, J. Weissmann, Distribution of materials in road earthmoving and paving: a mathematical programming approach, J. Constr. Eng. Manag. 139 (2012) 1046–1054, http://dx.doi.org/10.1061/(ASCE)CO. 1943-7862.0000666. [17] D. Li, C. Liu, M. Lu, Optimizing earthwork hauling plan with minimum cost flow network, International Construction Specialty Conference University of British Columbia, Vancouver, Canada. June 7-10, 2015, pp. 1–9, http://dx.doi.org/10.14288/ 1.0076334. [18] Y. Ji, F. Seipp, A. Borrmann, S. Ruzika, E. Rank, Mathematical modeling of earthwork optimization problems, Proceeding of the International Conference on Computing in Civil and Building Engineering (ICCCBE), 2010 http://www.cms.bgu.tum.de/publications/paper_Ji_ICCCBE2010.pdf. [19] R. Burdett, E. Kozan, R. Kenley, Block models for improved earthwork allocation planning in linear infrastructure construction, Eng. Optim. (2014) 1–23, http://dx. doi.org/10.1080/0305215X.2014.892594. [20] R.L. Burdett, E. Kozan, An integrated approach for earthwork allocation, sequencing and routing, Eur. J. Oper. Res. 238 (2015) 741–759, http://dx.doi.org/10.1016/j.ejor. 2014.04.036. [21] P. Chakroborty, Genetic algorithms for optimal urban transit network design, Comput. Aided Civ. Infrastruct. Eng. 18 (2003) 184–200, http://dx.doi.org/10. 1111/1467-8667.00309. [22] C.K. Wong, I.W.H. Fung, C.M. Tam, Comparison of using mixed-integer programming and genetic algorithms for construction site facility layout planning, J. Constr. Eng. Manag. 136 (2010) 1116–1128, http://dx.doi.org/10.1061/(ASCE)CO.1943-7862. 0000214. [23] K. Kepaptsoglou, M. Karlaftis, Transit route network design problem: review, J. Transp. Eng. 135 (2009) 491–505, http://dx.doi.org/10.1061/(ASCE)0733947X(2009)135:8(491). [24] R. Vaughan, Optimum polar networks for an urban bus system with a many-tomany travel demand, Transp. Res. B Methodol. 20 (1986) 215–224. [25] Alberta Infrastructure and Transportation, Chapter I access management guideline, (2005). https://www.transportation.alberta.ca/951.htm. [26] C. Huang, C.K. Wong, C.M. Tam, Optimization of tower crane and material supply locations in a high-rise building site by mixed-integer linear programming, Autom. Constr. 20 (2011) 571–580, http://dx.doi.org/10.1016/j.autcon.2010.11.023. [27] C. Huang, C.K. Wong, Optimisation of site layout planning for multiple construction stages with safety considerations and requirements, Autom. Constr. 53 (2015) 58–68, http://dx.doi.org/10.1016/j.autcon.2015.03.005. [28] H. Cancela, A. Mauttone, M.E. Urquhart, Mathematical programming formulations for transit network design, Transp. Res. B Methodol. 77 (2015) 17–37, http://dx. doi.org/10.1016/j.trb.2015.03.006. [29] P. Luathep, A. Sumalee, W.H.K. Lam, Z. Li, H.K. Lo, Global optimization method for mixed transportation network design problem: a mixed-integer linear programming approach, Transp. Res. B Methodol. 45 (2011) 808–827, http://dx.doi.org/10. 1016/j.trb.2011.02.002. [30] A.A.B. Pritsker, L.J. Waiters, P.M. Wolfe, Multiproject scheduling with limited resources: a zero-one programming approach, Manag. Sci. 16 (1969) 93–108. [31] PYTHON. [Computer Software]. PYTHON, Version 3.4, Python Software Foundation, 9450 SW Gemini Dr., ECM# 90772, Beaverton, OR 97008, USA, 2015. [32] CPLEX. [Computer Software]. CPLEX, Version 12.6, IBM Corporation, 1 New Orchard Road, Armonk, New York 10504-1722, United States, 2015.