# Optimizing Open-Pit Truck Route Based on Minimization of Time-Varying Transport Energy Consumption

Authors: Guangwei Liu; Senlin Chai
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

This paper addresses a special truck routing optimization problem in open-pit mines based on the minimization of time-varying transport energy consumption. A mixed-integer programming model is formulated to clearly describe the engineering problem, and a series of constraints are deduced to strengthen the model. To ensure that the model has time-varying characteristics, a method to estimate time-varying parameters is proposed by using pattern recognition and trend surface estimation. This time-varying resistance coefficient is mainly used to describe the process of road damage caused by frequent rolling of heavy trucks on the road surface. At the same time, in order to make the truck routing converge to the optimal energy consumption solution quickly, some definitions and properties are then provided based on stochastic theory, and a strategy to improve the computational efficiency of the model is proposed using these properties. Finally, an improved genetic algorithm is designed to solve the model. The results of experiments show that the proposed algorithm is effective and efficient.


---

## Body

## 1. Introduction

Truck’s routing along open-pit mines is an important link connecting all process flows and external transportation tasks for mineral and waste, which is one of the key factors that will directly influence the cost of transportation systems [1, 2]. Therefore, it has long been a focus of research in the fields of mining system engineering, logistics planning, and truck dispatching [3, 4]. According to the statistics, it is generally believed that transportation costs account for more than 50% of total operation costs in open-pit mines that use trucks as the major means of transportation [5–7]. Thus, most scholars, including the authors of this paper, think it economically useful to optimize the routes taken by trucks using a specific optimization algorithm to reduce overall transportation costs of open-pit mines.The transport routing optimization problem in open-pit mines, which ultimately aims to optimize global transportation cost under specific physical and economic conditions, is a kind of multiobjective combinatorial optimization problem. It has attracted the interest of scholars in the fields of logistics planning and truck dispatching for open-pit mines owing to its effect on reducing transportation costs. At present, methods used to optimize routes for open-pit mines can be divided into two categories. One of them is to try to use the static equivalent distance as a weight or selection criterion to search for the optimal route. For example, White and Olson [8] used the minimization of the sum of the transport distance converted to the plane to determine the optimal route for all locations. Chang et al. [9] proposed a mixed-integer programming model by considering different transport revenues, which is different from White’s decision-making method based on the shortest distance. Adenso-Diaz et al. [10] noted that there are frequent changes to transport routes in coal mining and proposed a method to automatically update the transport network. Li et al. [11] proposed a method to determine the optimal route for water trucks by minimizing the cost of operational and time delays. Hu [12] used the genetic algorithm and the harmony searching algorithm to solve routing problems of high computational complexity. Chen and Han [13], Sun and Liu [14], and Li et al. [15] used the swarm optimization algorithm to solve the routing problem for open-pit mines based on a static road transportation network. Methods of the second kind used to optimize trucks transportation route for open-pit mines involve optimizing routes by solving a mathematical programming problem under multiple constraints. Choi et al. [16] proposed a method to optimize routes for haulage based on the multicriterion evaluation and least cost path analysis and noted that few studies have been conducted on the optimal haul routes of trucks because, in previous studies, road surface conditions had been assumed to remain constant throughout its lifecycle. Choi and Nieto [17] proposed a modified least cost path analysis algorithm that is applicable to areas with all types of roads and can compare the influence of terrain and curve variation on transportation efficiency. Furthermore, a series of models have been developed in logistics [18–20], transportation [21–24], and Internet routing [25, 26], where each has particular focuses, advantages, and disadvantages. However, these methods are challenging to directly apply to the problem of route optimization for open-pit mine, because, in practical scenarios, transport conditions are constrained by a variety of comprehensive factors, and the transport cost model, which is stochastic and dynamic, is not equivalent to the shortest physical distance between nodes or the single static constants associated with transport cost. In particular, frequent rolling of heavy trucks for open-pit mines leads to damage to the road surfaces and causes changes in its characteristics, which is not considered in most routing optimization problems. Therefore, in a sense, the algorithm used should contain the following features: first, it must guarantee that the solution is an optimal route meeting the basic goal of minimizing transportation cost. Second, the transportation cost model should completely consider multiple, comprehensive factors, such as pavement quality, vehicle type, speed, load, rode slope, road maintenance, and driving technology. Finally, it should rapidly converge to the global optimal solution.In light of the above considerations, the problem of finding the shortest path in terms of time-varying transport energy consumption under several constraints is examined in this paper. First, we provide a brief overview of the characteristics of the road transportation network in open-pit mines and establish a mixed-integer programming model to describe the state of global time-varying transport. Second, to ensure that the model has time-varying characteristics, we propose a method to calculate resistance in different road segments. Finally, to improve the computational efficiency of the time-varying model, we propose several optimization strategies based on stochastic theory and improve the genetic algorithm to solve the mixed-integer programming model.The remainder of this paper is organized as follows: In Section2, the core problem of routing optimization in open-pit mines is introduced, the mixed-integer programming model is proposed, and a novel technique to estimate time-varying parameters is presented. In Section 3, we introduce an optimization strategy based on stochastic theory and provide sufficient evidence of the correctness of the improvement strategy. In Section 4, we improve the genetic algorithm and detailed the design process of the model. In Section 5, a simulation of an open-pit mine is conducted using raw data to show the application of the proposed method. Finally, the conclusions of this study are offered in Section 6.


## 2. Problem Description and Mathematical Model

### 2.1. Problem Description

In an open-pit mine with trucks as the main transport equipment, there are usually several routing nodes in the road transport system as shown in Figure1. In general, the routing optimization algorithm first builds transport network diagrams based on the nodes, and static network analysis is carried out using a directed graph that employs the minimum physical distance or transport energy consumption. It is worth noting that although these algorithms can achieve the given goal of route optimization to some extent, they cannot optimize the energy consumption of the selected routes. This is mainly because the energy consumption of trucks is directly influenced by the pavement quality, vehicle type, load capacity, road slope, driving technology, and the driver’s view, which lead to transportation energy comsumption having certain time characteristics. Thus, developing the time-varying expression of energy consumption is a crucial task for creating an optimization model with the target of dynamic energy consumption.Figure 1
3D model of open-pit mine (the main body of the figure is a three-dimensional (3D) model of an open-pit mine, and the left side shows a partial stripping transportation system of the north slope; the right side is a satellite image mainly used to illustrate the spatial relationship between the road and the bench).In the problem scenario that we consider, it is generally believed that the energy consumption of trucks is dominated by the work needed to overcome driving resistance. Thus, if optimization is based on the minimum global time-varying energy consumed for transport, we must first quantitatively analyze the fluctuation in driving resistance in different states. In general, for mining trucks, this resistance fluctuation effect can be divided into three parts [27]: rolling resistance, slope resistance, and air resistance, and the immediate cause of this fluctuation can be attributed to the following two factors: the randomness of the selected route and the time-varying fluctuation of the resistance coefficient owing to frequent damage to the road surface. The former is easy to understand, as the difference between the resistance coefficient in different road segments and the type of truck on the selected route inevitably leads to a variation in time-varying resistance. The time-varying fluctuation of the resistance coefficient on the same segment is mainly caused by the frequent rolling damage owing to heavy trucks in the maintenance period, where this kind of fluctuation exhibits extremely periodicity. To ensure that the model has higher logical rigor than the traditional route optimization model, the basis of decision-making of the proposed model is considered from the above two aspects.


### 2.2. Mixed-Integer Programming Model Based on Time-Varying Transport

To define and formulate this problem explicitly, the following parameters are introduced first:V=1,2,3,…,n: set of nonempty routing nodesE: sets of directed edges between adjacent nodesei,j=i,j∈E: edge between adjacent nodes i and j, i,j∈V,Fei,j: sets of resistance functions between adjacent nodes (units: kN)ΔH: sets of elevation differences between adjacent nodes (units: m)GV,E,F,ΔH: the directed graph of the road transport system networkNk: number of moving trucks, k∈K, K represents the type of vehicle∑k=1KNei,jk: number of moving trucks between adjacent routing nodes, k∈KMPS,PEO: the amount of waste transported to dumps along the route; PS: loading point; PE: unloading point (units: mm3)Mik: the load capacity of the ith type k truck (units: m3)mik: the weight of the ith type k truck (units: t)tPS,PEO: scheduled time for the dumping workloadv¯ei,j: average driving speed of trucks between adjacent routing nodes (km/h)ΔLei,j: the length of route between adjacent nodes (units: m)Sei,jE: the limiting space of safe driving between adjacent routing nodes (units: m)Kei,j: traffic density between adjacent nodes (units: vehicle/km)lk: length of k type vehicle, k∈K (units: m)Kei,jE: traffic density limit between adjacent nodes (units: vehicle/km)KE: traffic density limit in the selected route, i∈V (units: vehicle/km)N: road capacity of the selected route (units: mm3)p: number of lanes, 0.5 for dual lanes and one for single laneskb: imbalance coefficient of vehicle drivingq: car volume (units: m3)kr: truck working time utilization coefficientAik: windshield area of ith type k truck (units: m2)r: slope resistance coefficientC: air resistance coefficientρ: atmospheric density (units: kg/m2)α: road gradientft: function of the coefficient of rolling resistanceTo represent routing node decisions, we introduced the following decision variables:(1)xij=1,if node was used;0,otherwise,for i∈V.Based on the above notations, the open-pit transport routing optimization problem can be formulated as the following mixed-integer programming model:(2)Min∑i=1NFeij⋅xij⋅ΔLei+∑i=1,k∈KNΔheij⋅xij⋅Mik+mik,for j=i+1;i,j∈V.According to resistance characteristics described in the literature [28], the constraints of the mixed-integer programming model can be formulated as follows:(3)s.tFeij≈Ffeijt+Freij+Fweij=feijt⋅cosα⋅Mik+mik+r⋅Mik+mik+12⋅C⋅Aik⋅ρ⋅v¯ei,j2,(4)Kei,jE=Lei,j+lk⋅1000Sei,jE⋅Lei,j,(5)KE=∑i=1n−1Kei,jE⋅Lei,j∑i=1n−1Lei,j,(6)KE≥MPS,PEOtPS,PEO⋅∑i=1n∑kKMik/n⋅K⋅∑i=1n−1Lei,j⋅xi,(7)Kei,jE≥MPS,PEO⋅Lei,jtPS,PEO⋅∑i=1n∑kKMik/n⋅K⋅∑i=1n−1Lei,j⋅xi2,(8)N=1000⋅∑i=1n−1v¯ei,j⋅xi⋅p⋅kb⋅Sei,jE⋅lk⋅Mik⋅kr⋅tPS,PEO∑i=1n−1xi,(9)MPS,PEO≤N,(10)∑i=1n−1xi=B,i∈V;B⊂V,(11)2≤B≤n−2.The objective function (2) of the integer programming model minimizes the total time-varying transportation work, i.e., the sum of the time-varying transportation work done to overcome resistance along the arc of each route. Constraint (3) defines the formula for time-varying resistance along the transportation route, and external factors that cause the difference in time variation are converted into corresponding parameter indices to be quantified. Constraint (4) defines the formula for the limit of traffic along the road segments based on an adequate driving distance, whereas Constraint (5) can be used to assign weights by route length to derive the overall traffic density of the selected route. Constraint (6) provides a lower bound for traffic density through a proper road capacity determined by the amount of scheduled waste. Constraint (7) provides the limits on traffic density along the road segments to prevent a shortage of road capacity along the selected segment. Constraint (8) defines the road capacity, and Constraint (9) provides conditions on it determined by the capacity of the selected route in the given planning period. Constraints (10) and (11) limit the number of selected nodes to prevent loops in the selected route.


### 2.3. A Method to Estimate Parameters of Time-Varying Coefficient of Rolling Resistance

Through formula (2), it is easy to see that the most difficult parameter to calculate in the model is the time-varying coefficient of rolling resistance feijt, mainly because of the temporal fluctuations of this parameter on different road segments and maintenance cycles. Therefore, to better simulate the time-varying characteristics of the coefficient of rolling resistance along a road used by heavy trucks, multilabel pattern classification and trend surface technology were used to construct the parameter estimation algorithm. The basic idea underlying the overall design of the algorithm is that we first use multilabel pattern classification to build a pattern recognition model that could identify the type of road surface based on attribute data and then use the average driving speed, cumulative load, and continuous intervals of operation to identify the type of pavement. Finally, the coefficient of rolling resistance is predicted by using the time-sensitive fluctuation of the attribute index. The process of implementation of the algorithm is as below.


#### 2.3.1. Enumerating Types of Road Surfaces and Collecting Parameters of Pavement Indices

Tests of the coefficient of rolling resistance [27, 28] have shown that it is mainly related to the type of road surface, driving speed, tire structure, tire material, and tire pressure. Although these factors have a certain effect on resistance coefficient tests, the impact is weaker than that of road damage caused by frequent use by heavy trucks in open-pit mines. Therefore, in our problem scenario, we analyze the problem of fluctuations in the coefficient of resistance caused by frequent damage to the pavement. For this basic assumption, three factors are summarized: the average driving speed of road segment (the driving speed limits of different segments are shown in Table 1), the amount of accumulated waste (during the period of pavement maintenance), and the interval of the previous maintenance cycle.Table 1
Maximum driving speeds for various road types.
Road types
Maximum driving speed (km/h)
Loading truck
Unloading truck
+8%
0
−8%
+8%
0
−8%
Transport artery
35
40
40
45
50
50
Transport semi-main
30
35
35
40
45
45
Temporary track
25
30
30
35
40
40
Turn back road
15
20
15
20
25
20In the process of establishing the type of surface pattern recognition algorithm, the most critical operation involves classifying pavement types according to the differences in the coefficient of rolling resistance to establish a variety of labels of pavement type. Therefore, we first need to consider how to effectively create tags for the road surface.For this reason, we first classify periodic damage to the road, which is divided into eight types of road surfaces in Table2 (it is necessary to consider types of surface according to the specific condition of the mine, and the index in Table 2 represents only the results of classification of the Zhahanaoer Open-Pit Mine).Table 2
Rolling resistance for various types of road surfaces (Zhahanaoer Open-Pit Mine).
Various types of road surfaces
Rolling resistance
(1) Very hard, smooth roadway or dirt surface, no penetration or flexing
0.010∼0.018
(2) Hard, smooth, stabilized surfaced roadway without penetration under load, watered, maintained
0.018∼0.020
(3) Firm, smooth, rolling roadway with dirt or light surfacing, flexing surfacing, flexing slightly, maintained slightly, maintained fairly regularly, watered
0.020∼0.030
(4) Dirt roadway, rutted or flexing under load, little maintenance, no water, 25 mm tire penetration or flexing
0.030∼0.045
(5) Dirt roadway, rutted or flexing under load, little maintenance, no water, 50 mm tire penetration or flexing
0.045∼0.060
(6) Rutted dirt roadway, soft under travel, no maintenance, no stabilization, 100 mm tire penetration or flexing
0.060∼0.075
(7) Rutted dirt roadway, soft under travel, no maintenance, no stabilization, 200 mm tire penetration and flexing
0.075∼0.1
(8) Very soft, muddy, rutted roadway, 300 mm tire penetration, no flexing
0.1∼0.14The multilabel classification algorithm is supported by training data, and thus another task of algorithm design is to collect them. In this section, we used the Zhahanaoer Open-Pit Mine as an example to show the format of the training data in Figure2.Figure 2
The partial training data.


#### 2.3.2. Multilabel Pattern Classification Model

The core step of the parameter estimation algorithm is to classify the training data using pattern tags formed by the type of surface. As this kind of algorithm [29–31] is relatively mature, this part provides only a brief introduction to the classification principle under the current scenario.For the multilabel classification scenario considered in this paper, the first hypothesis is that the training sample isXi,yi,i=1,2.l, xi1,xi2,xi3∈Xi, xi∈R,yi∈+1,−1, and l denotes the total number of training samples. In this given multilabel case, such a nonlinear classification problem can be mapped to a high-dimensional space by the kernel function of an SVM (support vector machine) to achieve a linear partition in high-dimensional space. According to the Mercer condition, the optimal decision function of this classification model can be expressed as follows:(12)fx=sgn∑i=1lyiaiKx⋅xi+b.The above SVM technology is only a binary-class method that could find the type of goal in a variety of types of road surfaces. Therefore, to ensure that the algorithm has the capacity for multilabel classification based on the binary class, a decision tree SVM was introduced to carry out multilabel classification. The basic idea of the method is to set up an SVM model at each level of the decision tree and select one of the subclasses at each level. The above operation is repeated until all subclasses have been chosen. Figure3 shows a schematic diagram of four class principles based on the decision tree and SVM technology.Figure 3
Four-classification logic.


#### 2.3.3. Trend Surface Model and Estimating the Coefficient of Rolling Resistance

The goal of the multilabel pattern classification algorithm is to find the type of road surface estimated by using the index parameters at any given time. The main task here is to use grouping data to estimate the coefficient of rolling resistance on the basis of the type of road surface. The method of estimation used here is not unique, but we find that data for the Zha mine are suitable for the quadratic trend surface analysis method in this paper through the detection of goodness of fit and significance. We thus only introduce the basic principle of quadratic trend surface analysis in the next section. The trend surface model can be expressed as follows:(13)y^i=b0+b1xi2+b2xi3+b3xi22+b4xi2xi3+b5xi32,where y^i is the estimated coefficient of rolling resistance and xi1, which denotes the average speed of the segment, has little effect on the periodic damage analysis, because of which we remove it to reduce the number of dimensions of the problem.We now need only use the least squares’ model to fit the parameters of the model to construct an available trend surface. The least squares’ model can be written as below, and the trend surface parameter estimation process is shown in Figure4:(14)minyi−y^i2.Figure 4
The trend surface estimation of the quadratic polynomial.


## 2.1. Problem Description

In an open-pit mine with trucks as the main transport equipment, there are usually several routing nodes in the road transport system as shown in Figure1. In general, the routing optimization algorithm first builds transport network diagrams based on the nodes, and static network analysis is carried out using a directed graph that employs the minimum physical distance or transport energy consumption. It is worth noting that although these algorithms can achieve the given goal of route optimization to some extent, they cannot optimize the energy consumption of the selected routes. This is mainly because the energy consumption of trucks is directly influenced by the pavement quality, vehicle type, load capacity, road slope, driving technology, and the driver’s view, which lead to transportation energy comsumption having certain time characteristics. Thus, developing the time-varying expression of energy consumption is a crucial task for creating an optimization model with the target of dynamic energy consumption.Figure 1
3D model of open-pit mine (the main body of the figure is a three-dimensional (3D) model of an open-pit mine, and the left side shows a partial stripping transportation system of the north slope; the right side is a satellite image mainly used to illustrate the spatial relationship between the road and the bench).In the problem scenario that we consider, it is generally believed that the energy consumption of trucks is dominated by the work needed to overcome driving resistance. Thus, if optimization is based on the minimum global time-varying energy consumed for transport, we must first quantitatively analyze the fluctuation in driving resistance in different states. In general, for mining trucks, this resistance fluctuation effect can be divided into three parts [27]: rolling resistance, slope resistance, and air resistance, and the immediate cause of this fluctuation can be attributed to the following two factors: the randomness of the selected route and the time-varying fluctuation of the resistance coefficient owing to frequent damage to the road surface. The former is easy to understand, as the difference between the resistance coefficient in different road segments and the type of truck on the selected route inevitably leads to a variation in time-varying resistance. The time-varying fluctuation of the resistance coefficient on the same segment is mainly caused by the frequent rolling damage owing to heavy trucks in the maintenance period, where this kind of fluctuation exhibits extremely periodicity. To ensure that the model has higher logical rigor than the traditional route optimization model, the basis of decision-making of the proposed model is considered from the above two aspects.


## 2.2. Mixed-Integer Programming Model Based on Time-Varying Transport

To define and formulate this problem explicitly, the following parameters are introduced first:V=1,2,3,…,n: set of nonempty routing nodesE: sets of directed edges between adjacent nodesei,j=i,j∈E: edge between adjacent nodes i and j, i,j∈V,Fei,j: sets of resistance functions between adjacent nodes (units: kN)ΔH: sets of elevation differences between adjacent nodes (units: m)GV,E,F,ΔH: the directed graph of the road transport system networkNk: number of moving trucks, k∈K, K represents the type of vehicle∑k=1KNei,jk: number of moving trucks between adjacent routing nodes, k∈KMPS,PEO: the amount of waste transported to dumps along the route; PS: loading point; PE: unloading point (units: mm3)Mik: the load capacity of the ith type k truck (units: m3)mik: the weight of the ith type k truck (units: t)tPS,PEO: scheduled time for the dumping workloadv¯ei,j: average driving speed of trucks between adjacent routing nodes (km/h)ΔLei,j: the length of route between adjacent nodes (units: m)Sei,jE: the limiting space of safe driving between adjacent routing nodes (units: m)Kei,j: traffic density between adjacent nodes (units: vehicle/km)lk: length of k type vehicle, k∈K (units: m)Kei,jE: traffic density limit between adjacent nodes (units: vehicle/km)KE: traffic density limit in the selected route, i∈V (units: vehicle/km)N: road capacity of the selected route (units: mm3)p: number of lanes, 0.5 for dual lanes and one for single laneskb: imbalance coefficient of vehicle drivingq: car volume (units: m3)kr: truck working time utilization coefficientAik: windshield area of ith type k truck (units: m2)r: slope resistance coefficientC: air resistance coefficientρ: atmospheric density (units: kg/m2)α: road gradientft: function of the coefficient of rolling resistanceTo represent routing node decisions, we introduced the following decision variables:(1)xij=1,if node was used;0,otherwise,for i∈V.Based on the above notations, the open-pit transport routing optimization problem can be formulated as the following mixed-integer programming model:(2)Min∑i=1NFeij⋅xij⋅ΔLei+∑i=1,k∈KNΔheij⋅xij⋅Mik+mik,for j=i+1;i,j∈V.According to resistance characteristics described in the literature [28], the constraints of the mixed-integer programming model can be formulated as follows:(3)s.tFeij≈Ffeijt+Freij+Fweij=feijt⋅cosα⋅Mik+mik+r⋅Mik+mik+12⋅C⋅Aik⋅ρ⋅v¯ei,j2,(4)Kei,jE=Lei,j+lk⋅1000Sei,jE⋅Lei,j,(5)KE=∑i=1n−1Kei,jE⋅Lei,j∑i=1n−1Lei,j,(6)KE≥MPS,PEOtPS,PEO⋅∑i=1n∑kKMik/n⋅K⋅∑i=1n−1Lei,j⋅xi,(7)Kei,jE≥MPS,PEO⋅Lei,jtPS,PEO⋅∑i=1n∑kKMik/n⋅K⋅∑i=1n−1Lei,j⋅xi2,(8)N=1000⋅∑i=1n−1v¯ei,j⋅xi⋅p⋅kb⋅Sei,jE⋅lk⋅Mik⋅kr⋅tPS,PEO∑i=1n−1xi,(9)MPS,PEO≤N,(10)∑i=1n−1xi=B,i∈V;B⊂V,(11)2≤B≤n−2.The objective function (2) of the integer programming model minimizes the total time-varying transportation work, i.e., the sum of the time-varying transportation work done to overcome resistance along the arc of each route. Constraint (3) defines the formula for time-varying resistance along the transportation route, and external factors that cause the difference in time variation are converted into corresponding parameter indices to be quantified. Constraint (4) defines the formula for the limit of traffic along the road segments based on an adequate driving distance, whereas Constraint (5) can be used to assign weights by route length to derive the overall traffic density of the selected route. Constraint (6) provides a lower bound for traffic density through a proper road capacity determined by the amount of scheduled waste. Constraint (7) provides the limits on traffic density along the road segments to prevent a shortage of road capacity along the selected segment. Constraint (8) defines the road capacity, and Constraint (9) provides conditions on it determined by the capacity of the selected route in the given planning period. Constraints (10) and (11) limit the number of selected nodes to prevent loops in the selected route.


## 2.3. A Method to Estimate Parameters of Time-Varying Coefficient of Rolling Resistance

Through formula (2), it is easy to see that the most difficult parameter to calculate in the model is the time-varying coefficient of rolling resistance feijt, mainly because of the temporal fluctuations of this parameter on different road segments and maintenance cycles. Therefore, to better simulate the time-varying characteristics of the coefficient of rolling resistance along a road used by heavy trucks, multilabel pattern classification and trend surface technology were used to construct the parameter estimation algorithm. The basic idea underlying the overall design of the algorithm is that we first use multilabel pattern classification to build a pattern recognition model that could identify the type of road surface based on attribute data and then use the average driving speed, cumulative load, and continuous intervals of operation to identify the type of pavement. Finally, the coefficient of rolling resistance is predicted by using the time-sensitive fluctuation of the attribute index. The process of implementation of the algorithm is as below.


### 2.3.1. Enumerating Types of Road Surfaces and Collecting Parameters of Pavement Indices

Tests of the coefficient of rolling resistance [27, 28] have shown that it is mainly related to the type of road surface, driving speed, tire structure, tire material, and tire pressure. Although these factors have a certain effect on resistance coefficient tests, the impact is weaker than that of road damage caused by frequent use by heavy trucks in open-pit mines. Therefore, in our problem scenario, we analyze the problem of fluctuations in the coefficient of resistance caused by frequent damage to the pavement. For this basic assumption, three factors are summarized: the average driving speed of road segment (the driving speed limits of different segments are shown in Table 1), the amount of accumulated waste (during the period of pavement maintenance), and the interval of the previous maintenance cycle.Table 1
Maximum driving speeds for various road types.
Road types
Maximum driving speed (km/h)
Loading truck
Unloading truck
+8%
0
−8%
+8%
0
−8%
Transport artery
35
40
40
45
50
50
Transport semi-main
30
35
35
40
45
45
Temporary track
25
30
30
35
40
40
Turn back road
15
20
15
20
25
20In the process of establishing the type of surface pattern recognition algorithm, the most critical operation involves classifying pavement types according to the differences in the coefficient of rolling resistance to establish a variety of labels of pavement type. Therefore, we first need to consider how to effectively create tags for the road surface.For this reason, we first classify periodic damage to the road, which is divided into eight types of road surfaces in Table2 (it is necessary to consider types of surface according to the specific condition of the mine, and the index in Table 2 represents only the results of classification of the Zhahanaoer Open-Pit Mine).Table 2
Rolling resistance for various types of road surfaces (Zhahanaoer Open-Pit Mine).
Various types of road surfaces
Rolling resistance
(1) Very hard, smooth roadway or dirt surface, no penetration or flexing
0.010∼0.018
(2) Hard, smooth, stabilized surfaced roadway without penetration under load, watered, maintained
0.018∼0.020
(3) Firm, smooth, rolling roadway with dirt or light surfacing, flexing surfacing, flexing slightly, maintained slightly, maintained fairly regularly, watered
0.020∼0.030
(4) Dirt roadway, rutted or flexing under load, little maintenance, no water, 25 mm tire penetration or flexing
0.030∼0.045
(5) Dirt roadway, rutted or flexing under load, little maintenance, no water, 50 mm tire penetration or flexing
0.045∼0.060
(6) Rutted dirt roadway, soft under travel, no maintenance, no stabilization, 100 mm tire penetration or flexing
0.060∼0.075
(7) Rutted dirt roadway, soft under travel, no maintenance, no stabilization, 200 mm tire penetration and flexing
0.075∼0.1
(8) Very soft, muddy, rutted roadway, 300 mm tire penetration, no flexing
0.1∼0.14The multilabel classification algorithm is supported by training data, and thus another task of algorithm design is to collect them. In this section, we used the Zhahanaoer Open-Pit Mine as an example to show the format of the training data in Figure2.Figure 2
The partial training data.


### 2.3.2. Multilabel Pattern Classification Model

The core step of the parameter estimation algorithm is to classify the training data using pattern tags formed by the type of surface. As this kind of algorithm [29–31] is relatively mature, this part provides only a brief introduction to the classification principle under the current scenario.For the multilabel classification scenario considered in this paper, the first hypothesis is that the training sample isXi,yi,i=1,2.l, xi1,xi2,xi3∈Xi, xi∈R,yi∈+1,−1, and l denotes the total number of training samples. In this given multilabel case, such a nonlinear classification problem can be mapped to a high-dimensional space by the kernel function of an SVM (support vector machine) to achieve a linear partition in high-dimensional space. According to the Mercer condition, the optimal decision function of this classification model can be expressed as follows:(12)fx=sgn∑i=1lyiaiKx⋅xi+b.The above SVM technology is only a binary-class method that could find the type of goal in a variety of types of road surfaces. Therefore, to ensure that the algorithm has the capacity for multilabel classification based on the binary class, a decision tree SVM was introduced to carry out multilabel classification. The basic idea of the method is to set up an SVM model at each level of the decision tree and select one of the subclasses at each level. The above operation is repeated until all subclasses have been chosen. Figure3 shows a schematic diagram of four class principles based on the decision tree and SVM technology.Figure 3
Four-classification logic.


### 2.3.3. Trend Surface Model and Estimating the Coefficient of Rolling Resistance

The goal of the multilabel pattern classification algorithm is to find the type of road surface estimated by using the index parameters at any given time. The main task here is to use grouping data to estimate the coefficient of rolling resistance on the basis of the type of road surface. The method of estimation used here is not unique, but we find that data for the Zha mine are suitable for the quadratic trend surface analysis method in this paper through the detection of goodness of fit and significance. We thus only introduce the basic principle of quadratic trend surface analysis in the next section. The trend surface model can be expressed as follows:(13)y^i=b0+b1xi2+b2xi3+b3xi22+b4xi2xi3+b5xi32,where y^i is the estimated coefficient of rolling resistance and xi1, which denotes the average speed of the segment, has little effect on the periodic damage analysis, because of which we remove it to reduce the number of dimensions of the problem.We now need only use the least squares’ model to fit the parameters of the model to construct an available trend surface. The least squares’ model can be written as below, and the trend surface parameter estimation process is shown in Figure4:(14)minyi−y^i2.Figure 4
The trend surface estimation of the quadratic polynomial.


## 2.3.1. Enumerating Types of Road Surfaces and Collecting Parameters of Pavement Indices

Tests of the coefficient of rolling resistance [27, 28] have shown that it is mainly related to the type of road surface, driving speed, tire structure, tire material, and tire pressure. Although these factors have a certain effect on resistance coefficient tests, the impact is weaker than that of road damage caused by frequent use by heavy trucks in open-pit mines. Therefore, in our problem scenario, we analyze the problem of fluctuations in the coefficient of resistance caused by frequent damage to the pavement. For this basic assumption, three factors are summarized: the average driving speed of road segment (the driving speed limits of different segments are shown in Table 1), the amount of accumulated waste (during the period of pavement maintenance), and the interval of the previous maintenance cycle.Table 1
Maximum driving speeds for various road types.
Road types
Maximum driving speed (km/h)
Loading truck
Unloading truck
+8%
0
−8%
+8%
0
−8%
Transport artery
35
40
40
45
50
50
Transport semi-main
30
35
35
40
45
45
Temporary track
25
30
30
35
40
40
Turn back road
15
20
15
20
25
20In the process of establishing the type of surface pattern recognition algorithm, the most critical operation involves classifying pavement types according to the differences in the coefficient of rolling resistance to establish a variety of labels of pavement type. Therefore, we first need to consider how to effectively create tags for the road surface.For this reason, we first classify periodic damage to the road, which is divided into eight types of road surfaces in Table2 (it is necessary to consider types of surface according to the specific condition of the mine, and the index in Table 2 represents only the results of classification of the Zhahanaoer Open-Pit Mine).Table 2
Rolling resistance for various types of road surfaces (Zhahanaoer Open-Pit Mine).
Various types of road surfaces
Rolling resistance
(1) Very hard, smooth roadway or dirt surface, no penetration or flexing
0.010∼0.018
(2) Hard, smooth, stabilized surfaced roadway without penetration under load, watered, maintained
0.018∼0.020
(3) Firm, smooth, rolling roadway with dirt or light surfacing, flexing surfacing, flexing slightly, maintained slightly, maintained fairly regularly, watered
0.020∼0.030
(4) Dirt roadway, rutted or flexing under load, little maintenance, no water, 25 mm tire penetration or flexing
0.030∼0.045
(5) Dirt roadway, rutted or flexing under load, little maintenance, no water, 50 mm tire penetration or flexing
0.045∼0.060
(6) Rutted dirt roadway, soft under travel, no maintenance, no stabilization, 100 mm tire penetration or flexing
0.060∼0.075
(7) Rutted dirt roadway, soft under travel, no maintenance, no stabilization, 200 mm tire penetration and flexing
0.075∼0.1
(8) Very soft, muddy, rutted roadway, 300 mm tire penetration, no flexing
0.1∼0.14The multilabel classification algorithm is supported by training data, and thus another task of algorithm design is to collect them. In this section, we used the Zhahanaoer Open-Pit Mine as an example to show the format of the training data in Figure2.Figure 2
The partial training data.


## 2.3.2. Multilabel Pattern Classification Model

The core step of the parameter estimation algorithm is to classify the training data using pattern tags formed by the type of surface. As this kind of algorithm [29–31] is relatively mature, this part provides only a brief introduction to the classification principle under the current scenario.For the multilabel classification scenario considered in this paper, the first hypothesis is that the training sample isXi,yi,i=1,2.l, xi1,xi2,xi3∈Xi, xi∈R,yi∈+1,−1, and l denotes the total number of training samples. In this given multilabel case, such a nonlinear classification problem can be mapped to a high-dimensional space by the kernel function of an SVM (support vector machine) to achieve a linear partition in high-dimensional space. According to the Mercer condition, the optimal decision function of this classification model can be expressed as follows:(12)fx=sgn∑i=1lyiaiKx⋅xi+b.The above SVM technology is only a binary-class method that could find the type of goal in a variety of types of road surfaces. Therefore, to ensure that the algorithm has the capacity for multilabel classification based on the binary class, a decision tree SVM was introduced to carry out multilabel classification. The basic idea of the method is to set up an SVM model at each level of the decision tree and select one of the subclasses at each level. The above operation is repeated until all subclasses have been chosen. Figure3 shows a schematic diagram of four class principles based on the decision tree and SVM technology.Figure 3
Four-classification logic.


## 2.3.3. Trend Surface Model and Estimating the Coefficient of Rolling Resistance

The goal of the multilabel pattern classification algorithm is to find the type of road surface estimated by using the index parameters at any given time. The main task here is to use grouping data to estimate the coefficient of rolling resistance on the basis of the type of road surface. The method of estimation used here is not unique, but we find that data for the Zha mine are suitable for the quadratic trend surface analysis method in this paper through the detection of goodness of fit and significance. We thus only introduce the basic principle of quadratic trend surface analysis in the next section. The trend surface model can be expressed as follows:(13)y^i=b0+b1xi2+b2xi3+b3xi22+b4xi2xi3+b5xi32,where y^i is the estimated coefficient of rolling resistance and xi1, which denotes the average speed of the segment, has little effect on the periodic damage analysis, because of which we remove it to reduce the number of dimensions of the problem.We now need only use the least squares’ model to fit the parameters of the model to construct an available trend surface. The least squares’ model can be written as below, and the trend surface parameter estimation process is shown in Figure4:(14)minyi−y^i2.Figure 4
The trend surface estimation of the quadratic polynomial.


## 3. Property of and Strategies for the Optimal Solution

### 3.1. Definition and Property

In Section2, we define the open-pit routing optimization model based on the minimization of time-varying transport energy consumption. By observing the model, it is easy to see that the optimal solution designed in Section 2 improves the applicability of the traditional model due to the following two aspects: on the one hand, the model quantifies the fluctuation of transport energy consumption in a dynamic way. On the other hand, a method of estimation based on a statistical learning strategy is introduced to estimate time-varying parameters in the model. Although these improvements can perfectly simulate engineering practice, at the same time, these meticulous improvements incur significant computational performance problems. Thus, to solve these conflicts effectively, it is necessary to further introduce some concepts and properties of the optimization in advance.Definition 1.
In the network graphG=V,E,WT, V denotes a nonempty set of nodes, E denotes a set of arcs, and Wt denotes the set of time-varying transport energy consumption of the route, for t∈0,tm. In the probability space Ω,F,P, the function wt is a real-valued function depending on time and is defined on Ω, i.e., there is a real value wt for any t defined on Ω. For any real x, the subset of Ωt:wt≤x that makes wt≤x is an event, that is, F is measurable. Thus, for the stochastic process of resistance doing work, the time-varying transport energy consumption is a discrete stochastic variable dependent on time and G is a dynamic stochastic network.
It is noteworthy that the once the time-varying transport energy consumption is considered a stochastic variable, we can construct its probability density function and redefine the concept of the global time-varying shortest route as follows:Definition 2.
In the road network graphG, a set of stochastic variables can be expressed as V=w1t,w2t,…,wnt. The joint probability density function constructed by this set of variables can be expressed as fw1t,w2t,…,wnt. Thus, the probability density function for any route can be expressed by fr=∫…∫fw1t,…,wntdw1t,…,dwnt. Naturally, the global shortest route can be defined as the maximum probability of routes satisfying the event (li≤l1, li≤l2, … , li≤lk). li denotes the only route from the source node to the destination node, and i≠k.
It is noteworthy that even though the concept of the global shortest route is redefined in Definition 2, probability density estimation for each route remained difficult. Therefore, to simplify the calculation further, we summarize the following two optimization properties by using statistics:Property 1.
For any two stochastic variablesX and Y, EX is smaller than EY, which is equivalent to that PX≤Y is higher than 0.5.Proof.
Property 1 demonstrates that when two stochastic variables are compared, the probability of the one with a smaller mathematical expectation is greater than that of the other. To demonstrate the practical correctness of this property, the proof of is given as shown in formulae (15)–(27).Suppose that random variablesX and Y have the following expectation relationship shown in formula (15):(15)EE≤EY.Its equivalent form is as follows:(16)EX−EY≤0⟺EX−Y≤0.We introduce another random variableZ with the following equivalent relation with X and Y:(17)Z=X−Y.Therefore, the proposition can be rewritten as(18)EZ≤0.From the definition of a random variable, the expectation ofZ can be express as(19)EZ=∫−∞+∞Z⋅PZdZ=∫−∞0Z⋅PZdZ+∫0+∞Z⋅PZdZ.According to the scope of the domain, the integral is rewritten as follows:(20)∫−∞+∞Z⋅PZdZ=∫−∞0Z⋅PZ≤0dZ︸A+∫0+∞Z⋅PZ≥0dZ︸B.According to the known conditions of Property 1, we give the sign criteria as follows:(21)∵Z<0,PZ>0,Z>0,PZ>0.And, the initial conditions are given by formula (23):(22)∵EZ≤0.Therefore, the following absolute value relationship can be obtained:(23)∴∫−∞0Z⋅PZ≤0dZ≥∫0+∞Z⋅PZ≥0dZ.Therefore, the probability eventsZ≤0 and Z≥0 must have the following relation:(24)∴PZ≤0≥PZ≥0.Because the full probability event ofZ is one,(25)∵PZ≤0+PZ≥0=1.The following formula can be derived from formulas (29) and (30):(26)PZ≤0≥1−PZ≤0⟺2∗PZ≤0≥1.We thus get the following:(27)PZ≤0≥0.5⟺PX≤Y≥0.5.Property 2.
provided only a prior condition for the comparison between stochastic variables. To further illustrate the practical effectiveness of this property in the optimization of the algorithm, by combining with a realistic optimization scenario, the second property of the optimization problem is provided.Property 3.
When searching for the optimal route along routesL and K, the routes correspond to events L≤K and L>K. Using the conclusion of Property 1, we see that these two events have the following properties: the smaller the expected value of a route, the higher the likelihood that it becomes the global optimal route.


### 3.2. Optimization Strategies

Based on the theory of optimization given in Section3.1, this section describes the optimization strategies for realistic engineering scenarios.Assuming there werel feasible routes in the road network, the expression for the time-varying transport energy consumption for these routes constituted a set of l-dimensional stochastic variables Ω=W1t,W2t,…,Wlt. The joint probability density function of these variables could be expressed as follows:(28)fr=fw1t,…,wlt=∫…∫fw1t,…,wlt⋅dw1t,…,dwlt.Thus, the traditional open-pit truck routing optimization problem could be transformed into the event probability problem shown as follows:(29)Pγwγt≤w1t,…,wγt≤wlt=∫Gwγt≤w1t,…,wγt≤wlt…∫fw1t,…,wlt⋅dw1t,…,dwlt.Finally, considering the computational complexity of the mass transport network, the expectation of these stochastic variables is used in an algorithm of adaptive values according to Properties 1 and 2, and an improved genetic algorithm is used to improve the efficiency of route searching.


## 3.1. Definition and Property

In Section2, we define the open-pit routing optimization model based on the minimization of time-varying transport energy consumption. By observing the model, it is easy to see that the optimal solution designed in Section 2 improves the applicability of the traditional model due to the following two aspects: on the one hand, the model quantifies the fluctuation of transport energy consumption in a dynamic way. On the other hand, a method of estimation based on a statistical learning strategy is introduced to estimate time-varying parameters in the model. Although these improvements can perfectly simulate engineering practice, at the same time, these meticulous improvements incur significant computational performance problems. Thus, to solve these conflicts effectively, it is necessary to further introduce some concepts and properties of the optimization in advance.Definition 1.
In the network graphG=V,E,WT, V denotes a nonempty set of nodes, E denotes a set of arcs, and Wt denotes the set of time-varying transport energy consumption of the route, for t∈0,tm. In the probability space Ω,F,P, the function wt is a real-valued function depending on time and is defined on Ω, i.e., there is a real value wt for any t defined on Ω. For any real x, the subset of Ωt:wt≤x that makes wt≤x is an event, that is, F is measurable. Thus, for the stochastic process of resistance doing work, the time-varying transport energy consumption is a discrete stochastic variable dependent on time and G is a dynamic stochastic network.
It is noteworthy that the once the time-varying transport energy consumption is considered a stochastic variable, we can construct its probability density function and redefine the concept of the global time-varying shortest route as follows:Definition 2.
In the road network graphG, a set of stochastic variables can be expressed as V=w1t,w2t,…,wnt. The joint probability density function constructed by this set of variables can be expressed as fw1t,w2t,…,wnt. Thus, the probability density function for any route can be expressed by fr=∫…∫fw1t,…,wntdw1t,…,dwnt. Naturally, the global shortest route can be defined as the maximum probability of routes satisfying the event (li≤l1, li≤l2, … , li≤lk). li denotes the only route from the source node to the destination node, and i≠k.
It is noteworthy that even though the concept of the global shortest route is redefined in Definition 2, probability density estimation for each route remained difficult. Therefore, to simplify the calculation further, we summarize the following two optimization properties by using statistics:Property 1.
For any two stochastic variablesX and Y, EX is smaller than EY, which is equivalent to that PX≤Y is higher than 0.5.Proof.
Property 1 demonstrates that when two stochastic variables are compared, the probability of the one with a smaller mathematical expectation is greater than that of the other. To demonstrate the practical correctness of this property, the proof of is given as shown in formulae (15)–(27).Suppose that random variablesX and Y have the following expectation relationship shown in formula (15):(15)EE≤EY.Its equivalent form is as follows:(16)EX−EY≤0⟺EX−Y≤0.We introduce another random variableZ with the following equivalent relation with X and Y:(17)Z=X−Y.Therefore, the proposition can be rewritten as(18)EZ≤0.From the definition of a random variable, the expectation ofZ can be express as(19)EZ=∫−∞+∞Z⋅PZdZ=∫−∞0Z⋅PZdZ+∫0+∞Z⋅PZdZ.According to the scope of the domain, the integral is rewritten as follows:(20)∫−∞+∞Z⋅PZdZ=∫−∞0Z⋅PZ≤0dZ︸A+∫0+∞Z⋅PZ≥0dZ︸B.According to the known conditions of Property 1, we give the sign criteria as follows:(21)∵Z<0,PZ>0,Z>0,PZ>0.And, the initial conditions are given by formula (23):(22)∵EZ≤0.Therefore, the following absolute value relationship can be obtained:(23)∴∫−∞0Z⋅PZ≤0dZ≥∫0+∞Z⋅PZ≥0dZ.Therefore, the probability eventsZ≤0 and Z≥0 must have the following relation:(24)∴PZ≤0≥PZ≥0.Because the full probability event ofZ is one,(25)∵PZ≤0+PZ≥0=1.The following formula can be derived from formulas (29) and (30):(26)PZ≤0≥1−PZ≤0⟺2∗PZ≤0≥1.We thus get the following:(27)PZ≤0≥0.5⟺PX≤Y≥0.5.Property 2.
provided only a prior condition for the comparison between stochastic variables. To further illustrate the practical effectiveness of this property in the optimization of the algorithm, by combining with a realistic optimization scenario, the second property of the optimization problem is provided.Property 3.
When searching for the optimal route along routesL and K, the routes correspond to events L≤K and L>K. Using the conclusion of Property 1, we see that these two events have the following properties: the smaller the expected value of a route, the higher the likelihood that it becomes the global optimal route.


## 3.2. Optimization Strategies

Based on the theory of optimization given in Section3.1, this section describes the optimization strategies for realistic engineering scenarios.Assuming there werel feasible routes in the road network, the expression for the time-varying transport energy consumption for these routes constituted a set of l-dimensional stochastic variables Ω=W1t,W2t,…,Wlt. The joint probability density function of these variables could be expressed as follows:(28)fr=fw1t,…,wlt=∫…∫fw1t,…,wlt⋅dw1t,…,dwlt.Thus, the traditional open-pit truck routing optimization problem could be transformed into the event probability problem shown as follows:(29)Pγwγt≤w1t,…,wγt≤wlt=∫Gwγt≤w1t,…,wγt≤wlt…∫fw1t,…,wlt⋅dw1t,…,dwlt.Finally, considering the computational complexity of the mass transport network, the expectation of these stochastic variables is used in an algorithm of adaptive values according to Properties 1 and 2, and an improved genetic algorithm is used to improve the efficiency of route searching.


## 4. Heuristic Solution to the Model

### 4.1. Design of the Genetic Algorithm Solution

The open-pit truck route optimization problem is a kind of TSP (traveling salesman problem), which is a typical problem of combinatorial optimization. The outstanding feature of this kind of combinatorial optimization problem is that the search space increased sharply with an increase in the scale of the problem. Thus, in the practical open-pit route optimization process, the traditional method of enumeration struggles to obtain the optimal solution. To ensure that the algorithm has better global optimization ability and adaptability, a genetic algorithm [24, 32, 33] is introduced to solve the problem. It consists of the following steps:Step A1: Gene Coding. The encoding process first constructs a global unique real tag (χ) for each node, and the feasible route sets in the real scene are mapped by a sequence of these tags on genes at different locations. For example, in a certain route, there are seven nodes between the source and the destination, where the number of source nodes is 14 and that of destination nodes is 27. Therefore, this routing gene sequence can be encoded as follows:(30)Χγ=χ1,…,χ7=14,16,17,19,23,25,27.Step A2: Generating the Initial Population. In the transportation network, the shortest route in terms of Euclidean distance is not necessarily the optimal route in terms of transport energy consumption, but it has the potential to become the optimal route compared with other routes. Therefore, to better control genes in terms of the direction of diversity and quality, this paper uses the shortest Euclidean distance as the prior condition for population initialization and proposes the following heuristic population initialization strategy:Step A2-1: calculate the Euclidean distance between nodes segmentdi,j, and assign it as the distance weight to routing arcωi,jd=di,j.Step A2-2: calculate the number of degrees of the output of each nodeoid and maintain a vector v⇀=Vi,Vi+1 that represents its direction.Step A2-3: add the source node to the position of the first gene and search for the adjacent node with the minimum valueωi,jd.Step A2-4: judge the out-degree of the node. Ifoid>1 and −v⇀∉V⇀, add the searched node to the position of the second gene. Add v⇀ to V⇀ and make oid−=1.Step A2-5: loop Step 2 until the desired population size is reached.Step A3: Calculating the Fitness of Individuals in a Population. For the open-pit routing optimization problem, we hope to obtain the minimum value of the objective function. For this reason, the expectation of the time-varying transport energy consumption along a route is used as an adaptive value, and the specific formal expression is as follows: fitnessi=D/Ewit, where D is a large constant intended to prevent fitness from getting close to zero as Ewit was too large, which leads to local optima.Step A4: Gene Selection. In nature, the role of natural selection is to retain the best adaptive environment and eliminate the suboptimal. In the optimization of the genetic algorithm, gene selection acts as the law of nature. The algorithm uses fitness as an index to evaluate the quality of genes in the process of breeding populations, and the genetic mechanism is used to control the genetic quality of the entire population to improve the optimization capability of the algorithm. Therefore, to preserve the best member of the parent generation, the population gene is selected using the elite selection strategy and roulette, which consists of the following steps:Step A4-1: selectε% of excellent individuals in the population as elite and directly copy them into the population of the next generation.Step A4-2: calculate the fitness of the individual in the remaining population1−ε% and its adaptive weight according to the formula λi=Fξi/∑j=1nFξj, which is considered the probability that the individual may be chosen.Step A4-3: using the roulette method, divide the probability internally. Create random numbersγε0,1 and ensure that the total number of random numbers is identical to the remaining population and simulate roulette behavior.Step A4-4: according to the formula∑k=1r−1pk≤γ<∑k=1rpk, the frequency of each interval is calculated, the high-frequency route is retained, and the low-frequency route is eliminated.Step A5: Heuristic Crossover Operator with Correction Strategy. The crossover method proposed in this paper uses single-point intersection. For a gene misallocated by single-point crossing, population initialization is used for its heuristic correction, and the correction method is shown in Figure 5.Step A6. Mutation Operator. The mutation operation used in this paper consists of two main parts. In the first, the algorithm first simulates two random numbers. One of them is used to specify the position of the gene for mutation, and the other is used to refer to the genetic code following mutation. In this step, the selected genetic code needs to be checked to ensure that the individual gene does not appear to be in the wrong position. A diagram of the mutation operation is shown in Figure 6.Step A7. Termination Conditions. If the termination condition is met, the iteration is stopped, or is moved to Step 2.Figure 5
Cross-operation.Figure 6
Mutation operation.


### 4.2. Parameter Optimization Strategies

To avoid the phenomenon of premature convergence of the algorithm in the initial iteration, the algorithm redesigns the crossover and mutation parameters to ensure that they could be adaptively adjusted according to the fluctuation of fitness.Strategy 1.
Adaptive adjustment of the crossover probability. The crossover operator is important to ensure that the algorithm produced new offspring, and the probability of the crossover operationPc directly restricts the rate of population update. The main effect is that when the rate of update of the population is large, it is unfavorable for the population to retain superior individuals. Conversely, when the rate of update of the population is small, the population could not evolve rapidly. Therefore, it is challenging to avoid the above problems if we do not consider the iterative process but use the same probability to handle all crossover operations. To ensure that the algorithm could obtain better opportunities for updates in the early iterations and better retain excellent characteristics of the population in later iterations, we design the following iterative relationship:(31)Pc=α−k/G3∗Pcm+k/G3∗Pc′2,fi−f¯>δ,minPc,fi−f¯≤δ,for α∈1,1.5, Pcm=0.8. k denotes the number of iterations, G denotes the total number of iterations, and Pc′ denotes the results of iterations of the crossover probability of the previous generation. In the first iteration, the initial value of Pc′ can be specified as 0.75. fi denotes the average fitness of the population in the ith iteration and f¯ the average fitness of the entire population in the iteration process. δ denotes the limit of the population’s fitness and average fitness.Strategy 2.
Adaptive adjustment of mutation probability. In the entire genetic algorithm, gene mutation is a key operation to control the diversity of the population. Its iterative nature is often the opposite of the crossover operation, and its characteristic is that if the mutation rate is too high in early iterations, it will lead to the algorithm getting caught during blindly searching and further lose the genetic characteristics of excellent individuals. Conversely, if the mutation rate is too low in later iterations, it is easy for the algorithm to fall into a local optimal solution, which makes it difficult to be found the global optimal solution. Therefore, to avoid the above problems, we design the following iterative relationship concerning mutation probability:(32)Pm=Pm′+Pm0−Pm′∗D∗fi−f¯fmax−f¯+kG,for Pm0=0.05, Pm0 denotes the initial value of the mutation probability, Pm′ denotes the results of iteration of the mutation probability of the previous generation, constant D denotes the equilibrium coefficient, fmax denotes the maximum fitness of the population of the previous generation, and the other parameters are the same as above.


## 4.1. Design of the Genetic Algorithm Solution

The open-pit truck route optimization problem is a kind of TSP (traveling salesman problem), which is a typical problem of combinatorial optimization. The outstanding feature of this kind of combinatorial optimization problem is that the search space increased sharply with an increase in the scale of the problem. Thus, in the practical open-pit route optimization process, the traditional method of enumeration struggles to obtain the optimal solution. To ensure that the algorithm has better global optimization ability and adaptability, a genetic algorithm [24, 32, 33] is introduced to solve the problem. It consists of the following steps:Step A1: Gene Coding. The encoding process first constructs a global unique real tag (χ) for each node, and the feasible route sets in the real scene are mapped by a sequence of these tags on genes at different locations. For example, in a certain route, there are seven nodes between the source and the destination, where the number of source nodes is 14 and that of destination nodes is 27. Therefore, this routing gene sequence can be encoded as follows:(30)Χγ=χ1,…,χ7=14,16,17,19,23,25,27.Step A2: Generating the Initial Population. In the transportation network, the shortest route in terms of Euclidean distance is not necessarily the optimal route in terms of transport energy consumption, but it has the potential to become the optimal route compared with other routes. Therefore, to better control genes in terms of the direction of diversity and quality, this paper uses the shortest Euclidean distance as the prior condition for population initialization and proposes the following heuristic population initialization strategy:Step A2-1: calculate the Euclidean distance between nodes segmentdi,j, and assign it as the distance weight to routing arcωi,jd=di,j.Step A2-2: calculate the number of degrees of the output of each nodeoid and maintain a vector v⇀=Vi,Vi+1 that represents its direction.Step A2-3: add the source node to the position of the first gene and search for the adjacent node with the minimum valueωi,jd.Step A2-4: judge the out-degree of the node. Ifoid>1 and −v⇀∉V⇀, add the searched node to the position of the second gene. Add v⇀ to V⇀ and make oid−=1.Step A2-5: loop Step 2 until the desired population size is reached.Step A3: Calculating the Fitness of Individuals in a Population. For the open-pit routing optimization problem, we hope to obtain the minimum value of the objective function. For this reason, the expectation of the time-varying transport energy consumption along a route is used as an adaptive value, and the specific formal expression is as follows: fitnessi=D/Ewit, where D is a large constant intended to prevent fitness from getting close to zero as Ewit was too large, which leads to local optima.Step A4: Gene Selection. In nature, the role of natural selection is to retain the best adaptive environment and eliminate the suboptimal. In the optimization of the genetic algorithm, gene selection acts as the law of nature. The algorithm uses fitness as an index to evaluate the quality of genes in the process of breeding populations, and the genetic mechanism is used to control the genetic quality of the entire population to improve the optimization capability of the algorithm. Therefore, to preserve the best member of the parent generation, the population gene is selected using the elite selection strategy and roulette, which consists of the following steps:Step A4-1: selectε% of excellent individuals in the population as elite and directly copy them into the population of the next generation.Step A4-2: calculate the fitness of the individual in the remaining population1−ε% and its adaptive weight according to the formula λi=Fξi/∑j=1nFξj, which is considered the probability that the individual may be chosen.Step A4-3: using the roulette method, divide the probability internally. Create random numbersγε0,1 and ensure that the total number of random numbers is identical to the remaining population and simulate roulette behavior.Step A4-4: according to the formula∑k=1r−1pk≤γ<∑k=1rpk, the frequency of each interval is calculated, the high-frequency route is retained, and the low-frequency route is eliminated.Step A5: Heuristic Crossover Operator with Correction Strategy. The crossover method proposed in this paper uses single-point intersection. For a gene misallocated by single-point crossing, population initialization is used for its heuristic correction, and the correction method is shown in Figure 5.Step A6. Mutation Operator. The mutation operation used in this paper consists of two main parts. In the first, the algorithm first simulates two random numbers. One of them is used to specify the position of the gene for mutation, and the other is used to refer to the genetic code following mutation. In this step, the selected genetic code needs to be checked to ensure that the individual gene does not appear to be in the wrong position. A diagram of the mutation operation is shown in Figure 6.Step A7. Termination Conditions. If the termination condition is met, the iteration is stopped, or is moved to Step 2.Figure 5
Cross-operation.Figure 6
Mutation operation.


## 4.2. Parameter Optimization Strategies

To avoid the phenomenon of premature convergence of the algorithm in the initial iteration, the algorithm redesigns the crossover and mutation parameters to ensure that they could be adaptively adjusted according to the fluctuation of fitness.Strategy 1.
Adaptive adjustment of the crossover probability. The crossover operator is important to ensure that the algorithm produced new offspring, and the probability of the crossover operationPc directly restricts the rate of population update. The main effect is that when the rate of update of the population is large, it is unfavorable for the population to retain superior individuals. Conversely, when the rate of update of the population is small, the population could not evolve rapidly. Therefore, it is challenging to avoid the above problems if we do not consider the iterative process but use the same probability to handle all crossover operations. To ensure that the algorithm could obtain better opportunities for updates in the early iterations and better retain excellent characteristics of the population in later iterations, we design the following iterative relationship:(31)Pc=α−k/G3∗Pcm+k/G3∗Pc′2,fi−f¯>δ,minPc,fi−f¯≤δ,for α∈1,1.5, Pcm=0.8. k denotes the number of iterations, G denotes the total number of iterations, and Pc′ denotes the results of iterations of the crossover probability of the previous generation. In the first iteration, the initial value of Pc′ can be specified as 0.75. fi denotes the average fitness of the population in the ith iteration and f¯ the average fitness of the entire population in the iteration process. δ denotes the limit of the population’s fitness and average fitness.Strategy 2.
Adaptive adjustment of mutation probability. In the entire genetic algorithm, gene mutation is a key operation to control the diversity of the population. Its iterative nature is often the opposite of the crossover operation, and its characteristic is that if the mutation rate is too high in early iterations, it will lead to the algorithm getting caught during blindly searching and further lose the genetic characteristics of excellent individuals. Conversely, if the mutation rate is too low in later iterations, it is easy for the algorithm to fall into a local optimal solution, which makes it difficult to be found the global optimal solution. Therefore, to avoid the above problems, we design the following iterative relationship concerning mutation probability:(32)Pm=Pm′+Pm0−Pm′∗D∗fi−f¯fmax−f¯+kG,for Pm0=0.05, Pm0 denotes the initial value of the mutation probability, Pm′ denotes the results of iteration of the mutation probability of the previous generation, constant D denotes the equilibrium coefficient, fmax denotes the maximum fitness of the population of the previous generation, and the other parameters are the same as above.


## 5. Computational Experiments

### 5.1. Application Examples and Algorithm Analysis

To test the performance of the mathematic model and the proposed algorithm with improvement strategies, we wrote all relevant programs under in C# .NET to solve for route optimization in the Zhahanaoer Open-pit Mine. The experimental platform and parameters are as follows:(1)
Experimental platform: all experiments were performed on a computer with a 2.9 GHz Intel Core I5 CPU, 8 GB of RAM, and Win 7 X64 operating system.(2)
GA parameters: population size,Popsize=150, the maximum number of iterations, Gmax=150, initial value of cross-probability, Pcm=0.8, and the initial value of the probability of mutation, Pm0=0.05.(3)
Other parameters of the calculation: the Zha mine uses two types of trucks: 1200 HP and 2000 HP. The capacity of the trucks was 60 m3 and 90 m3, respectively. The slope is r=0.08, the area of the wind shield is A=52.48m2 and 61.05m2, respectively, and the waste bulk density is 1.68t/m3.


#### 5.1.1. Application Examples

To intuitively express the effect of the optimization of the algorithm, we select 5 sets of test cases and calculate the energy consumption by generations. These five groups of examples are the transportation routings from mining and stripping sites to unloading sites of Zhahanaoer Mine in 2018. The energy consumption relation in the process of searching the optimal solution for each group of routings is as shown in Figure7.Figure 7
Energy consumption by generations.From Figure7, we see that the solution quickly converges to the minimum value of the transport energy consumption in the iterative process. This shows that the algorithm could quickly find the optimal solution for global energy consumption, and that is effective for solving such problems.


#### 5.1.2. Analysis of the Effect of Parameter Optimization

To further demonstrate the effectiveness and advantages of parameters optimization, we use all test cases to build simulation experiments, and the results of a comparison are shown in Figure8.Figure 8
Energy consumption by generations (parameters not optimized).Figure8 shows that the algorithm for the parameters that is not optimized is more likely to fall into a local optimum. In cases 1 and 4, it could not effectively find the global optimal solution until the end of the iteration. However, the other sets are also slow to converge. Based on the above analysis, we see that the improvement in parameters of the genetic algorithm plays an active role in accelerating its optimization.


### 5.2. Analysis of the Superiority of the Adopted Algorithm

To further demonstrate the superiority of the algorithm in solving the problem of route optimization, this section reports a comparative analysis of the optimization results obtained by the improved algorithm and the traditional algorithm. The data are shown in Table3.Table 3
Data for comparison of the effect of algorithmic optimization.
Case
Load status
Contrast algorithm
Node
Average speed (Km/h)
Haul distance (Km)
Transport energy (KJ)
Execution time (s)
N1
1
IGA
29
37
1.39
746932.9
1.59
Dijkstra/PSO
30
1.27
776826.4
7.07/4.72
N2
1
IGA
31
34
1.17
677493.7
2.02
Dijkstra/PSO
27
1.08
704826.5
8.77/3.04
N3
1
IGA
43
36
1.68
898642.1
3.97
Dijkstra/PSO
32
1.59
914732.3
11.87/3.73
N4
1
IGA
59
36
1.93
959358.6
4.01
Dijkstra/PSO
30
1.87
967439.5
15.02/5.35
N5
1
IGA
74
27
1.94
915632.7
5.16
Dijkstra/PSO
22
1.67
979261.8
21.17/7.79From Table3, we see that these algorithms find the optimal route in practical settings. However, there are obvious differences in terms of energy consumption and execution efficiency of these algorithms. We see that Dijkstra’s algorithm and the PSO focus on the optimal problem of a static network and obtain the global equivalent minimum distance quickly. However, this minimum equivalent distance is not equivalent to the minimum energy consumed. Therefore, the algorithm used here is more reliable than traditional optimization algorithms when resistance is constantly changing owing to frequent damage to the road surface.


## 5.1. Application Examples and Algorithm Analysis

To test the performance of the mathematic model and the proposed algorithm with improvement strategies, we wrote all relevant programs under in C# .NET to solve for route optimization in the Zhahanaoer Open-pit Mine. The experimental platform and parameters are as follows:(1)
Experimental platform: all experiments were performed on a computer with a 2.9 GHz Intel Core I5 CPU, 8 GB of RAM, and Win 7 X64 operating system.(2)
GA parameters: population size,Popsize=150, the maximum number of iterations, Gmax=150, initial value of cross-probability, Pcm=0.8, and the initial value of the probability of mutation, Pm0=0.05.(3)
Other parameters of the calculation: the Zha mine uses two types of trucks: 1200 HP and 2000 HP. The capacity of the trucks was 60 m3 and 90 m3, respectively. The slope is r=0.08, the area of the wind shield is A=52.48m2 and 61.05m2, respectively, and the waste bulk density is 1.68t/m3.


### 5.1.1. Application Examples

To intuitively express the effect of the optimization of the algorithm, we select 5 sets of test cases and calculate the energy consumption by generations. These five groups of examples are the transportation routings from mining and stripping sites to unloading sites of Zhahanaoer Mine in 2018. The energy consumption relation in the process of searching the optimal solution for each group of routings is as shown in Figure7.Figure 7
Energy consumption by generations.From Figure7, we see that the solution quickly converges to the minimum value of the transport energy consumption in the iterative process. This shows that the algorithm could quickly find the optimal solution for global energy consumption, and that is effective for solving such problems.


### 5.1.2. Analysis of the Effect of Parameter Optimization

To further demonstrate the effectiveness and advantages of parameters optimization, we use all test cases to build simulation experiments, and the results of a comparison are shown in Figure8.Figure 8
Energy consumption by generations (parameters not optimized).Figure8 shows that the algorithm for the parameters that is not optimized is more likely to fall into a local optimum. In cases 1 and 4, it could not effectively find the global optimal solution until the end of the iteration. However, the other sets are also slow to converge. Based on the above analysis, we see that the improvement in parameters of the genetic algorithm plays an active role in accelerating its optimization.


## 5.1.1. Application Examples

To intuitively express the effect of the optimization of the algorithm, we select 5 sets of test cases and calculate the energy consumption by generations. These five groups of examples are the transportation routings from mining and stripping sites to unloading sites of Zhahanaoer Mine in 2018. The energy consumption relation in the process of searching the optimal solution for each group of routings is as shown in Figure7.Figure 7
Energy consumption by generations.From Figure7, we see that the solution quickly converges to the minimum value of the transport energy consumption in the iterative process. This shows that the algorithm could quickly find the optimal solution for global energy consumption, and that is effective for solving such problems.


## 5.1.2. Analysis of the Effect of Parameter Optimization

To further demonstrate the effectiveness and advantages of parameters optimization, we use all test cases to build simulation experiments, and the results of a comparison are shown in Figure8.Figure 8
Energy consumption by generations (parameters not optimized).Figure8 shows that the algorithm for the parameters that is not optimized is more likely to fall into a local optimum. In cases 1 and 4, it could not effectively find the global optimal solution until the end of the iteration. However, the other sets are also slow to converge. Based on the above analysis, we see that the improvement in parameters of the genetic algorithm plays an active role in accelerating its optimization.


## 5.2. Analysis of the Superiority of the Adopted Algorithm

To further demonstrate the superiority of the algorithm in solving the problem of route optimization, this section reports a comparative analysis of the optimization results obtained by the improved algorithm and the traditional algorithm. The data are shown in Table3.Table 3
Data for comparison of the effect of algorithmic optimization.
Case
Load status
Contrast algorithm
Node
Average speed (Km/h)
Haul distance (Km)
Transport energy (KJ)
Execution time (s)
N1
1
IGA
29
37
1.39
746932.9
1.59
Dijkstra/PSO
30
1.27
776826.4
7.07/4.72
N2
1
IGA
31
34
1.17
677493.7
2.02
Dijkstra/PSO
27
1.08
704826.5
8.77/3.04
N3
1
IGA
43
36
1.68
898642.1
3.97
Dijkstra/PSO
32
1.59
914732.3
11.87/3.73
N4
1
IGA
59
36
1.93
959358.6
4.01
Dijkstra/PSO
30
1.87
967439.5
15.02/5.35
N5
1
IGA
74
27
1.94
915632.7
5.16
Dijkstra/PSO
22
1.67
979261.8
21.17/7.79From Table3, we see that these algorithms find the optimal route in practical settings. However, there are obvious differences in terms of energy consumption and execution efficiency of these algorithms. We see that Dijkstra’s algorithm and the PSO focus on the optimal problem of a static network and obtain the global equivalent minimum distance quickly. However, this minimum equivalent distance is not equivalent to the minimum energy consumed. Therefore, the algorithm used here is more reliable than traditional optimization algorithms when resistance is constantly changing owing to frequent damage to the road surface.


## 6. Conclusions

This study solved the problem of truck route optimization in open-pit mines under the influence of a fluctuation in resistance. A mixed-integer programming model was formulated, and time-varying parameters were estimated and used to simulate the characteristics of the time-varying rolling resistance in different segments. To solve the problem of computational complexity caused by the addition of time-varying characteristics, some effective improvement strategies based on stochastic theory were proposed, and an improved genetic algorithm was used to solve the model. Finally, the validity of the algorithm was verified through numerical experiments.



---

*Source: 6987108-2019-08-28.xml*