Reliability-based performance analysis of mining drilling operations through Markov chain Monte Carlo and mean reverting process simulations

Abstract

In recent years, commodity prices have swiftly decreased, narrowing the profit margin for many mining operations and forcing them to find effective cost management strategies to respond to low prices. Given that equipment is one of the most significant assets of a mining company, efficient equipment utilization has strong potential to reduce costs. This paper focuses on the relationship between the number of available drilling machines based on reliability analysis and the number of holes to be created on a bench of an open pit mining operation. Since equipment availability is random in nature, a range of holes to be drilled corresponding to a specified probability level was determined. To assess the performance of the proposed approach, a case study was carried out using two stochastic modeling techniques. Evolutions of reliabilities of 10 rotary drilling machines over a specific time were simulated by Markov chain Monte Carlo and mean reverting processes, using historical data. Multiple simulations were then used for risk quantification. Results show that the proposed approach can be used as a tool to assist production scheduling and assess the associated risk.

Keywords

Open pit mining Markov chain Markov chain Monte Carlo mean reverting process reliability analysis risk quantification

1. Introduction

Even though there are alternative rock fragmentation approaches, for example, microwave assisted breakage, drilling and blasting are the most widely used technique in hard-rock open pit mining operations. Drilling is one of the primary operations in an open pit mining cycle. It is a complex operation because it is affected by several factors, such as geological structure, equipment condition, drilling parameters, and operator experience.¹ Failures in drilling equipment severely affect production schedules because they cause blasting delays and affect the subsequent production process. If limits of operations are known, production plans can be more realistic. Equipment condition is an essential factor to achieve a desirable production rate and to sustain mining operation.² Reliability analysis is the best tool to ensure the condition of the equipment. Reliability analysis helps practitioners to forecast future failure of system components and prevent unwanted stoppages.^3,4

Equipment reliability can be used as a key system performance metric during the equipment lifetime.^5,6 Field data are required to analyze system reliability and construct reliability models to facilitate decision-making. A non-repairable system is the system in which a repair is expensive or non-feasible. On the other hand, repairable systems are systems that can be restored after a failure for satisfactory operation.⁷ Time between failures and time to repair are primary data types needed to characterize system reliability.⁸ Drilling machines are repairable systems that can be restored after failure to perform desired performance.

It is impossible to predict the future availability of drilling equipment in a fleet. Therefore, there is risk (uncertainty) associated with the number of holes to be drilled. If there are delays due to lack of available equipment, subsequent processes (e.g., blasting, loading, and hauling) will also be delayed such that production targets are not attained.⁹ To address the risks arising from available drilling equipment, stochastic modeling is conducted—using past failure and repair data—to simulate the available equipment in the fleet through multiple future scenarios. The stochastic simulations methods are used to calculate a range of production quantities.¹⁰

Mining engineering involves a significant number of risks due to heterogeneities of geologic and geotechnics phenomena. Mining operations consist of consecutive activities (drilling, blasting, loading, hauling, and crushing). Given that the drilling is the first activity, delays in drilling due to equipment availability or insufficient equipment performance will also result in delays in subsequent activities. Therefore, production targets may not be attained. To address the risks arising from available drilling equipment, stochastic processes can be used to generate future probable realizations.¹⁰ Markov chain Monte Carlo (MCMC) simulation and mean reversion (MR) techniques have strong potential to quantify the risks associated with drilling operations, which can be significant.¹⁰ MCMC is a mathematical model for stochastic systems describing a series of possible events in a time series, whereby the probability of each event depends only on the condition of the preceding event.¹¹ MCMC generates probable realizations depending on the current condition. The MR theory states that variables eventually move toward and oscillate around the equilibrium level, which can be calculated from historical data.¹² Thus, MR is used to create possible scenarios depending on the historical average.

MCMC simulation fits well the nature of the problem because the change in available machines in the next time increment does not fluctuate much. In other words, the number of available pieces of equipment in the drilling operation for the next shift depends on the previous shift. Furthermore, due to degradation, equipment availability will decrease over time. This phenomenon can be modeled through MCMC. On the other hand, MR is also suited well to simulate the number of available drilling machines. It assumes that random increments are generated from a normal distribution. This is a reasonable assumption because the decrease in equipment availability is governed by the availability of initial time and the long-term mean of the process. The biggest issue in MR is to calibrate the parameters. In this research, since information from the previous benches were collected, the calibration is a relatively easy task because there is an opportunity to observe deviations of calibrated parameters from actual realizations. Two different stochastic approaches are used to see if their outcomes agree.

Although both methods simulate equipment availability for a given duration, there are differences between their applications. MCMC uses transition probabilities, and the current state depends upon the previous state. The transition probabilities are computed from previous experience. This computation is based on a formula. Therefore, the technique is quite fast. In addition, the ability to reduce multidimensional problems to a series of lower-dimensional ones is one of the most important characteristics of MCMC. However, the so-called “memoryless” character of MCMC is the biggest drawback. The assumption of exponential distribution for time to failure is also critical.^11,13 MCMC should be implemented carefully after the validity of these assumptions. On the other hand, MR is mainly used in finance to simulate long-term future prices. It requires some parameters, such as the mean reverting factor, long-term mean, and constant volatility factor. These parameters are calibrated from historical data. The values of parameters, to a certain extent, depends on the quality of the maintenance program in the mine and the heterogeneity level of rock characteristics. MR assumes that equipment availability tends to be the average availability over the time. This assumption is highly related to the size and quality of maintenance activities. Also, MR is very sensitive to outliers and data noisiness.¹²

In recent years, some researchers have focused on the performance measurement of drilling machines from different aspects. Ataei et al.¹⁴ investigated the physical and mechanical properties of rock to measure the penetration rate, unlike in this research, where we computed the penetration rate from historical data and further investigated rock and machine interaction. Basarir et al.¹⁵ developed a model to predict the performance of drilling machines using an adaptive neuro-fuzzy inference system and multiple regression. However, since the machine condition was ignored, it is difficult to have reliable outcomes. In our research, we linked production scheduling, which focuses on the quantity of material to be extracted for a given set of equipment. Furthermore, Al-Chalabi et al.¹⁶ conducted a study for an underground mine drilling rig to build a process to simulate the reliability of repairable complex systems based on historical data by ordinary Monte Carlo simulation. In their research, the previous failures do not affect future events. In other words, it is time-independent, unlike in this research, applying MCMC simulation, where time dependency was considered.

In this research, the relationship between the reliability of drilling equipment and its performance is quantified. Then, equally probable realizations of the available number of drilling machines over time series are generated through MCMC and MR, based on historical data. These realizations are used to assess the feasibility of the targeted production plans. The originality of this paper lies in proposing a risk quantification approach, which assists mine management to (1) determine production rates (mine production scheduling) based on drilling performance and (2) develop drilling equipment maintenance plans (preventive maintenance and spare part management).

2. Research methods

After collecting field data over one year, a power law model was applied for each drilling machine. The power law method, also known as Crow-AMSAA, was used to analyze repairable complex systems. The power law model parameters were calculated by Reliasoft© RGA software. To investigate the trend of the time between failure datasets, the Laplace trend test was used to show if the system behavior was improving or deteriorating. Historical data were also used to investigate the relationship between machine reliability and machine performance by JMP© statistical software. In addition, the number of available drilling machines for each shift was simulated separately by MCMC (100 simulations in ModelRisk© software) and MR (100 simulations in Microsoft Excel©) for a three-month period. Finally, the range of drillable holes was generated according to the number of available drilling machines.

2.1. Reliability analysis

Reliability analysis helps to deal with uncertainty and to make an informed decision. The general expression for the function of reliability is given as follows⁶:

R (t) = \Pr {T \geq t}

(1)

where R(t) is the reliability at time interval t, T is the time to failure of the system or item and R(t) ≥ 0, R (0) = 1. Other expressions of the reliability function are presented in the following⁷:

R (t) = 1 - F (t) = 1 - \int_{0}^{t} f (x) dx

(2)

R (t) = e^{- \int_{0}^{t} λ (t) dt}

(3)

where F(t) is the cumulative failure distribution function, f(x) is the failure probability density function, and λ(t) is the hazard rate.

The reliability expressions given above are used to determine the reliability of a system or an item for which the time to failure is characterized by statistical distributions, such as exponential, normal, or Weibull, if the system behaves “as good as new” after the repair.¹⁶ The failure process is called the renewal process. This basic model is called the Homogenous Poisson Process (HPP).¹⁶ The reliability function for three-Parameter Weibull distribution is given by Equation (4),⁵ where the three defining parameters of Weibull distribution are the shape parameter (β), also known as the Weibull slope, the scale parameter (η), and the location parameter (γ), also known as the shift parameter. These systems are known as independent and identically distributed (i.i.d.) when there is no trend at the dataset. However, for the complex systems, such as trucks, loaders, and drilling machines, the failures are dependent based on the current age of the remaining components. Therefore, most cases of complex systems are between “as good as new” and “as bad as old” conditions after the repair and the deterioration trend can be seen. This process is called the non-renewable process. To characterize the reliability of drilling machines, the Non-homogenous Poisson Process (NHPP), which is a generalization of the Poisson process, can be used instead of distributions. It has a wide applicability to model repairable systems^16,17:

R (t) = e^{- {(\frac{t - γ}{η})}^{β}}

(4)

In reliability life data analysis, failure data are needed to calculate the time between failures.⁵ A power law model was fitted to the data to analyze the life of each machine for the duration of the operation. The initial transition and steady-state matrices were created to determine the transitions between stages on shifts based on reliability analysis by the Markov chain (MC) technique. The last step was to determine the number of available drilling machines and the expected number of holes that can be drilled. The problem was solved by MCMC and MR.

Drilling machines can be modeled by the power law method. In terms of an analytical method to investigate the trend of the time between failure datasets, the Laplace trend test is used in this research. It shows whether the system behavior is improving or deteriorating.

The power law technique is used to model the system failure intensity function to manage each succeeding system failure (Equation (5)). β and λ were estimated from Equations (6) and Equation (7)¹⁸:

u (t) = λ β t^{β - 1}, t > 0

(5)

\hat{β} = \frac{\sum_{q = 1}^{K} N_{q}}{\sum_{q = 1}^{K} \sum_{i = 1}^{N_{q}} \ln (\frac{T_{q}}{X_{iq}})}

(6)

\hat{λ} = \frac{\sum_{q = 1}^{K} N_{q}}{(T_{1}^{β}) + (T_{2}^{β}) + \dots + (T_{K}^{β})}

(7)

where K is the number of drilling machines, N_q is the total number of failures for each system, and T is the observation time for each failure dataset. The estimated β value can also indicate the trend. If β = 1, there is no trend; if β >1, the system is degrading; and if β <1, the system is improving.¹⁷

The Cramer–Von Mises test is the most suitable goodness of fit test to analyze multiple repairable systems that follow the power law model.¹⁹ If a calculated result is less than the critical value from the goodness of fit test table, it fails to reject the NHPP power model. The power law model mean repair function and reliability function, defined as the probability of zero failure from time t to t+s, for NHPP, can be seen in Equations (8) and Equation (9), respectively⁶:

M (t) = \int_{0}^{t} λ (t) dt ≅ λ t^{β}

(8)

R (dt) = e^{- [M (t + dt) - M (t)]}

(9)

2.2. Markov chain

The mathematical definition of the MC can be seen in the following equation²⁰:

\begin{matrix} P {S_{t + 1} = j | S_{0} = k_{0}, S_{1} = k_{1}, . . ., S_{t - 1} = k_{t - 1}, S_{t} = i} \\ = P {S_{t + 1} = j | S_{t} = i} \end{matrix}

(10)

where S_t is a stochastic process. For reliability modeling, the probability of being in state j in time t+1(d_t) when it is in state i at time t can be formulated as follows²⁰:

\begin{matrix} P {E_{t + 1} = j | E_{t} = i} = {E 1_{j} (t + 1) | E 1_{i} (t), \\ E 2_{j} (t + 1) | E 2_{i} (t), . . ., E m_{j} (t + 1) | E m_{i} (t)} \end{matrix}

(11)

where the integers (1, 2, . . ., m) represent the number of pieces of equipment, and i and j represent the current and future states, respectively.

Equipment reliability data were used to estimate the number of available drilling machines in a certain period by the MC. The power law model tends to represent the life data of the system. It is used to estimate the parameters to make the function fit the data closely.²⁰ Equation (12) shows the probability of being in operation after one unit of time (d_t) when the equipment is working, and Equation (13) shows the probability of being under repair after one unit of time (d_t) when the equipment is working, based on time between failure datasets:

\begin{matrix} P {E_{t + dt} = j | E_{t} = i} \\ = {E m_{j} (t + dt) = 1 | E m_{i} (t) = 1} = e^{- (λ (t^{β}))} \end{matrix}

(12)

\begin{matrix} P {E_{t + dt} = j | E_{t} = i} \\ = {E m_{j} (t + dt) = 0 | E m_{i} (t) = 1} = 1 - e^{- (λ (t^{β}))} \end{matrix}

(13)

Similarly, Equation (14) is used to calculate the probability of being in operation after one unit of time (d_t) when the equipment is under the repair condition and Equation (15) is used to calculate the probability of continuing to stay under repair after one unit of time (d_t) when the equipment is under repair the condition based on time-to-repair datasets:

\begin{matrix} P {E_{t + dt} = j | E_{t} = i} \\ = {E m_{j} (t + dt) = 1 | E m_{i} (t) = 0} = 1 - e^{- (λ (t^{β}))} \end{matrix}

(14)

\begin{matrix} P {E_{t + dt} = j | E_{t} = i} \\ = {E m_{j} (t + dt) = 0 | E m_{i} (t) = 0} = e^{- (λ (t^{β}))} \end{matrix}

(15)

The one-step transition matrix for equipment was formulated from previous equations, depending on the states. Multiplication is needed to determine the probability of transitioning one state to another for more than one piece of equipment:

\begin{matrix} P {E_{t + dt} = j | E_{t} = i} \\ = {E 1_{j} (t + dt) | E 1_{i} (t)} \times \dots \times {E m_{j} (t + dt) | E m_{i} (t)} \end{matrix}

(16)

After a certain number of transitions, the system will reach a steady state that is independent of the current state and has constant probability. At this state, the transitioning probabilities in a certain state are independent of the probability distribution of the initial state.

More information about developing MC models for repairable systems and mining operations can be found in the literature.^10,11,20

2.3. Markov chain Monte Carlo

The MCMC technique was used to generate samples from complex distributions generated by the MC method. These samples were then used by the Monte Carlo method to quantify estimation and model the risk.²¹ There are many algorithms and sampling methods to implement MCMC. The Metropolis–Hastings algorithm is a frequently used way to set up the MC.²² It creates random samples from the target distribution to form a MC that conditionally depends only upon the last event.²²

For the ordinary Monte Carlo theory, samples are generated randomly and distributed as independent and identical.^22–24 However, MCMC is used to create samples that depend only on the previous sample based on the MC, which is stationary and reversible. The difference between ordinary Monte Carlo and MCMC can be seen by the formulation of variance $(σ^{2})$ in Equations (17) and (18), respectively¹¹:

σ^{2} = var {g (X)}

(17)

σ^{2} = var {g (X)} + 2 \sum_{k = 1}^{\infty} cov {g (X_{i}), g (X_{i + k})}

(18)

where g(X) is a real-valued function on the state space.

The variance of the sample at Equation (17) is i.i.d. That is, the variance of the sample in Equation (18) is the function of the variance of the previous sample.

According to the Metropolis–Hastings algorithm, the probability of a proposed move from i to j is given by the following^11,22:

r (i, j) = \frac{h (j) q (j, i)}{h (i) q (i, j)}

(19)

a (i, j) = \min (1, r (i, j))

(20)

where h is the un-normalized density, q is the conditional probability density, r is the Hastings ratio, and a is the probability of moving from i to j.

The choice of the proposal density has a significant impact on the performance of the algorithm. The convergence characteristics of the implemented MCMC will be highly related to the choice of the proposal density. To generate samples, Metropolis sampling is required. In this point, the proposal distribution q(S∣S(t–1)) and the prior distribution π (0) over the initial state in the MC need to be selected. In this research, Gaussian distribution was used for both distributions. The prior distribution is centered at zero μ = S(t–1) and (σ = 1) and the proposal distribution is centered at the previous state of the MC (μ = 0 and σ = 1).

Convergence of the MCMC to its stationary distribution is a requirement. Unfortunately, there are no universally accepted approaches to prove convergence. In this research, we used software. To our best knowledge, it utilized the technique proposed by Gelman and Rubin.¹⁹ The technique has two stages. In the first, the target distribution is estimated and using this distribution the starting point is produced. Thus, the required number of independent chains is met. In the second, the target distribution of the scalar quantity under consideration (e.g., a Student’s t distribution and the scale parameter) is re-constructed through the last k iterations.

The MCMC has a serious impact on solving in a wide range of stochastic problems. Pang et al.²⁵ used MCMC to estimate wind speed distribution, and Malhotra²⁶ used it for applications in network and computer security. Similarly, Ozdemir and Kumral¹³ tested fleet efficiency and Mardia et al.²⁷ implemented MCMC to model rock fractures.

2.4. Mean reversion

MR is used to create future observations using historical data. The MR stochastic process can be formulated as follows²⁸:

d x_{t} = κ (θ - x_{t}) dt + σ \sqrt{x_{t}} d Z_{t}

(21)

where x_t is the process level at the initial time, κ is the speed of reversion, θ is the long-term mean level, σ is the volatility, and Z is the increment of standard Brownian Motion.

Parameters κ, θ, and σ can be forecast from a regression equation based on historical data as follows²⁸:

\frac{Δ x_{t}}{x_{t}} = β_{0} + β_{1} \frac{1}{x} + ϵ

(22)

The variable κ is the negative of the intercept (–β₀), θ is the negative ratio of the coefficient of 1/x_t of the intercept, and σ is the standard error of the residuals.

To simulate the generated MR model, a starting value was calculated from historical data, and independent normally distributed error values (uniformly distributed between 0 and 1) were generated using Microsoft Excel©.

MR is used to perform risk analysis and decision-making. Therefore, it has rightfully received attention from research in finance and asset management. Detailed information about using this technique can be found in the literature.^28–33

3. Case study

3.1. Reliability analysis

Historical data were collected from an open pit mine over a one-year period (Table 1). The data were obtained from 10 drilling machines, which are 320XPC rotary blasthole drills. These machines can provide up to 68 tons of bit loading and the maximum single-pass hole depth is up to 20 m. The dimension of the bit used was 12¼ inches, which is around 32 cm. The time between failures and the time under repair were recorded for each machine. In addition, under the same condition, the drilling length was calculated for each drilling machine to investigate the relationship between reliability and drilling performance.

Table 1.

Summary of historical data collected from 10 rotary drilling machines.

Machine	Total no.failures	Uptime	Downtime	Mean time between failures	Mean time under repair	Total time
		(Hours)
1	72	4428	1794	61.50	25.27	6222
2	100	4688	1751	46.89	17.69	6439
3	86	4565	1406	53.09	16.35	5971
4	85	5238	1269	61.63	14.93	6507
5	54	3455	2039	64.00	37.77	5494
6	73	4438	1536	60.81	21.34	5974
7	89	4498	1762	50.54	20.02	6260
8	111	4566	1923	41.14	17.32	6489
9	97	4902	1425	50.54	14.84	6327
10	113	4441	1737	39.31	15.51	6178

The parameters of the power law model are listed in Table 2 and the reliability plots of the machines are shown in Figure 1.

Table 2.

Parameters of the power law model for each of the 10 rotary drilling machines.

Machine	Power law parameters
	β	λ
1	1.1412	0.00036
2	1.1941	0.00019
3	1.0993	0.00024
4	1.1274	0.00015
5	1.1924	0.00031
6	1.1763	0.00019
7	1.0988	0.00037
8	1.1809	0.00031
9	1.0578	0.00028
10	1.0437	0.00049

Figure 1.

Reliability plots of the rotary drilling machines.

The reliabilities of the drilling machines vary. For example, between 0 and 1000 hours, the probability of fulfilling the intended functions of Machine 5 is only approximately 30%, whereas the probability for Machine 4 in the same time range is approximately 70%.

Once the reliability analysis and the reliability levels were obtained, the drilling length was calculated for each drilling machine based on reliability levels by regression analysis:

y_{n} = 5.2764 \times e^{0.0011 (x_{[n - (n - 20)]})}

(23)

where y_n is the time required to drill for given length (min), x is the starting point length of the drilling with the same drillbit, n is the length increment (m), and 20 is the length of the drill hole. To make a meaningful comparison, the same length for each equipment was used. As can be seen from Figure 1, all equipment reliabilities after 1400 hours are less than 60%: the lower limit of the effective working range. Therefore, the drilling length was calculated for 1400 hours of operation.

Equation (20) shows the relationship between drilling time and drilling length when the variables of the drilling machine (i.e., rotation speed (rev/min), pulldown force (MPa), and bailing air pressure (MPa)) are constant at 80, 150, and 1.6, respectively (standard drilling operation application based on the rock characteristics according to the manufacturer’s recommendation).³⁴

For simplicity, the drilling length (20 m) was converted to the number of drill holes using JMP© Statistical software to plot the relationship between reliability and drilling performance by linear regression analysis as follows (Table 3):

noh = - 934.1 + (2, 663.7 \times r)

(24)

where noh is the number of drillable holes and r is the reliability level.

Table 3.

Linear regression output to estimate the number of drill holes from the drilling length.

Term	Estimated value	p-value	R ²
Intercept	−934.1	0.0003	0.76
Reliability	2663.7	<0.0001	0.76

There is a direct association between the number of drillable holes and reliability, which is a criterion of drilling performance.

3.2. Markov chain

The number of the available drilling machines and the number of drillable holes can be generated by the MC and simulated by MCMC. There are two possible conditions for a machine: in operation (1) or under repair (0). Hence, for 10 drilling machines, there are 1024 (2¹⁰) possible states (summarized in Table 4).

Table 4.

Possible states for 10 drilling machines.

Possible states	Drilling machines
	1	2	3	4	5	6	7	8	9	10
1	0	0	0	0	0	0	0	0	0	0
2	1	0	0	0	0	0	0	0	0	0
3	0	1	0	0	0	0	0	0	0	0
4	0	0	1	0	0	0	0	0	0	0
.	.	.	.	.	.	.	.	.	.	.
.	.	.	.	.	.	.	.	.	.	.
.	.	.	.	.	.	.	.	.	.	.
1022	1	1	1	1	1	1	1	1	0	1
1023	1	1	1	1	1	1	1	1	1	0
1024	1	1	1	1	1	1	1	1	1	1

The initial probability matrix of each piece of equipment’s condition was calculated by Reliasoft© RGA software based on equipment reliability data. The probability of failure and probability of repair were obtained for 10 hours of operation, which is one shift (Table 5).

Table 5.

The probability (%) of changing conditions for all drilling machines.

Machine	1 → 0	1 → 1	0 → 1	0 → 0
	Probability of failure (%)	Probability of not failure (%)	Probability of repaired (%)	Probability of not repaired (%)
1	2.55	97.45	0.90	99.10
2	0.11	99.89	2.33	97.67
3	0.29	99.71	1.58	98.42
4	2.64	97.36	3.93	96.07
5	0.04	99.96	0.03	99.97
6	1.48	98.53	0.45	99.55
7	1.56	98.44	0.90	99.10
8	1.67	98.33	4.04	95.96
9	0.94	99.06	0.72	99.28
10	0.51	99.49	4.74	95.26

If Machine 1 is working, the probability of being under repair is 2.55%, and the probability of being in operation is 97.45% after one shift. On the other hand, if Machine 1 is under repair, the probability of being in operation is 0.90% and the probability of staying under the repair condition is 99.10%.

The initial probability matrix shows the probability of transitioning between stages on shift. Transitioning probability from one stage to another was calculated based on the condition changing probabilities (Table 5), and ModelRisk© software was used to obtain transitioning probabilities at different shifts. The one-step transition matrix (1024 × 1024) based on the initial probability matrix is summarized in Table 6. The number of available drilling machines is calculated for each state, and the states are grouped based on the number of available machines.

Table 6.

Initial transition matrix (%).

State	0	1	2	3	4	5	6	7	8	9	10
0	81.90	16.63	1.40	0.06	0.00	0.00	0.00	0.00	0.00	0.00	0.00
1	0.99	82.76	15.08	1.13	0.04	0.00	0.00	0.00	0.00	0.00	0.00
2	0.01	1.99	83.58	13.51	0.88	0.03	0.00	0.00	0.00	0.00	0.00
3	0.00	0.03	3.02	84.36	11.90	0.67	0.02	0.00	0.00	0.00	0.00
4	0.00	0.00	0.07	4.06	85.11	10.27	0.48	0.01	0.00	0.00	0.00
5	0.00	0.00	0.00	0.11	5.12	85.82	8.62	0.32	0.01	0.00	0.00
6	0.00	0.00	0.00	0.00	0.17	6.20	86.49	6.94	0.19	0.00	0.00
7	0.00	0.00	0.00	0.00	0.00	0.24	7.29	87.13	5.24	0.10	0.00
8	0.00	0.00	0.00	0.00	0.00	0.01	0.33	8.40	87.72	3.51	0.03
9	0.00	0.00	0.00	0.00	0.00	0.00	0.01	0.43	9.52	88.27	1.77
10	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.01	0.54	10.66	88.78

As an example, if any four drilling machines are in operation, the probability that any five drilling machines will be in operation after one shift (10 hours of operation) is 10.27%. On the other hand, the probability that any three drilling machines will be in operation is 4.06%, and the probability that the same number will stay in operation is 85.11%. It should be noted that the 0.00% probability in the matrices represents small but non-zero probabilities.

After a specified time, the system will reach the steady-state level where the transitioning probabilities do not change in time and the states are constant. The steady-state matrix was calculated by ModelRisk© software and is summarized in Table 7.

Table 7.

Steady-state matrix (%).

State	1	2	3	4	5	6	7	8	9	10
0	0.04	0.55	3.69	12.86	25.11	28.75	19.55	7.70	1.61	0.14
1	0.04	0.55	3.69	12.86	25.11	28.75	19.55	7.70	1.61	0.14
2	0.04	0.55	3.69	12.86	25.11	28.75	19.55	7.70	1.61	0.14
3	0.04	0.55	3.69	12.86	25.11	28.75	19.55	7.70	1.61	0.14
4	0.04	0.55	3.69	12.86	25.11	28.75	19.55	7.70	1.61	0.14
5	0.04	0.55	3.69	12.86	25.11	28.75	19.55	7.70	1.61	0.14
6	0.04	0.55	3.69	12.86	25.11	28.75	19.55	7.70	1.61	0.14
7	0.04	0.55	3.69	12.86	25.11	28.75	19.55	7.70	1.61	0.14
8	0.04	0.55	3.69	12.86	25.11	28.75	19.55	7.70	1.61	0.14
9	0.04	0.55	3.69	12.86	25.11	28.75	19.55	7.70	1.61	0.14
10	0.04	0.55	3.69	12.86	25.11	28.75	19.55	7.70	1.61	0.14

When the steady-state matrix is calculated, the number of available drilling machines can be simulated for future shifts. Table 8 shows the probabilities of having available drilling machines. The initial matrix helps generate the number of available drilling machines for the next shift. On the other hand, the steady-state matrix helps generate the number of available machines for future shifts when the system reaches the steady state.

Table 8.

The probability of having available drilling machines.

No. machines in operation	Probability (%)
0	0.00
1	0.04
2	0.55
3	3.69
4	12.86
5	25.11
6	28.75
7	19.55
8	7.70
9	1.61
10	0.14

At steady state, the probability of having five drilling machines working is 25.11% and the probability of having six drilling machines working is 28.75%.

As mentioned above, reliability data were used to create the initial matrix to implement the MC. This process shows that each drilling machine is unique. Even if all equipment is brand new at the beginning of the operation, the reliabilities and availabilities cannot be the same over time due to various factors, such as geological heterogeneity, human factors, quality of maintenance, and random events. Besides, mining operations use typically varying equipment ages. Therefore, equipment performance is different at any given time point. The probabilities of having available drilling machines will be different, as shown in Table 9.

Table 9.

Probability of having available drilling machines when machines are assumed identical.

No. machines in operation	Probability (%)
0	0.01
1	0.09
2	0.68
3	3.04
4	8.90
5	17.85
6	24.81
7	23.60
8	14.70
9	5.42
10	0.90

Tables 8 and 9 show that, particularly for long-term plans, it can be misleading to assume that all machines are identical. Therefore, a reliability analysis of each machine is necessary to obtain accurate results.

3.3. Markov chain Monte Carlo simulation

The number of available drilling machines for different shifts was generated by MCMC using ModelRisk© software based on the initial matrix and start vector, which is the initial state where eight of the 10 drilling machines are available. The estimation facilitates visualizing the variation in number of available drilling machines between shifts. Multiple scenarios were generated by replicating MCMC. Figure 2 demonstrates possible outcomes of uncertainties for 4 of 100 randomly generated scenarios of available drilling machines. Because of seasonal effects, 90 consecutive days (180 shifts) were simulated. Over the time, similar results were obtained by different simulations.

Figure 2.

Number of available drilling machines simulated by MCMC.

The general trend falls between five and eight available drilling machines. The degradation trend is evident, as expected.

Once 100 scenarios were generated and the number of drillable holes was calculated for 180 shifts. The parameters of the drilling machine, rotation speed, pulldown force, and bailing air pressure, were set at 80 rpm, 150 MPa, and 1.6 MPa, respectively. It was assumed that drill bits are changed for every 1400 m (70 holes), which is the level of effective drilling. Results are illustrated in Figure 3.

Figure 3.

Number of drillable holes simulated by MCMC for 180 shifts.

3.4. Mean reversion simulation

In the same manner, the number of available drilling machines was simulated by MR using Microsoft Excel©. The parameters of the simulation calculated from historical data are presented in Table 10.

Table 10.

Parameters of mean reversion simulation (see Equation (23)).

Number of available machines at initial time (X₀)	8
Speed (κ)	108
Long-term mean (θ)	6.6
Volatility (σ)	6

Multiple scenarios were generated by the MR process and 100 randomly generated simulations were obtained for 180 shifts. The four simulations that exemplify the uncertainties (Figure 4) show that the available number of drilling machines fluctuates between four and eight for most of the shifts, similar to the MCMC results.

Figure 4.

Number of available drilling machines simulated by MR for 180 shifts.

The number of available holes was also calculated for the same drilling circumstances (Figure 5).

Figure 5.

The number of drillable holes simulated by MR for 180 shifts.

The comparison between MCMC and MR is presented in Table 11. The results show that their means are almost same but MCMC has lower standard deviation than MR.

Table 11.

Comparison between Markov chain Monte Carlo (MCMC) and mean reversion (MR).

Number of drillable holes
Parameter	MCMC	MR
Mean	59,325	59,235
Standard deviation	778	966

According to the MCMC simulation, the probability that the number of drillable holes is likely to be between 59,000 and 61,000 is 94% (Figure 4). MR results show 87% probability that the number of drillable holes is between the same ranges (Figure 5). There is 7% difference between two methods to create a production schedule. There will be around 3% deviation according to given interval (2000 out of around 60,000 holes). Given the nature of a mining operation, this is quite acceptable to install production capacity.

To validate the approach, the simulation results and actual realizations are compared. In actual application, the number of drilled holes was 60,391. As can be seen, this result is within the range obtained by simulations. Therefore, as long as there will not be a significant change in rock characteristics, the rest of the operation can be designed with respect to the simulation results.

4. Conclusion

This paper presents an approach based on a combination of reliability analysis, MC theory, and MR to simulate the number of available drilling machines and the number of drillable holes for a given probability level. Using simulation results, a range of production rates can be generated for a pre-specified probability. Firstly, the power law model was applied to time between failure data to determine β and λ parameters. Then, reliability analysis was conducted to characterize the behavior of drilling machines. Finally, the number of available machines was simulated by the MCMC technique based on reliability analysis and the MR process using historical data. The results of the MCMC simulation indicated a more than 80% probability that five to eight drilling machines will be available for every 10-hour period (one shift). Using these processes will facilitate more accurate decision-making for production scheduling and risk management. It can be seen from the results of the number of drillable holes that there will be approximately only 3% deviation to the planned drilling schedule by using both MCMC and MR processes.

Also, the association between drilling machine reliability level and performance was quantified for 10 drilling machines. The direct relationship between reliability and performance was demonstrated by regression analysis. The comparison of the reliability of similar drilling machines was illustrated by statistical analysis.

The paper also discusses the assumption that all machines have identical reliability. The MC results showed that this assumption could be misleading if long-term plans are considered. In this way, the necessity of implementing the reliability analysis for the decision-making mechanism and risk management were shown.

Future research will use these research results to develop preventive maintenance and short-term mine planning.

Footnotes

Authors’ Note

Omer Faruk Ugurlu is also affiliated with Department of Mining and Materials Engineering, McGill University, Canada.

Funding

This work was supported by the Natural Sciences and Engineering Research Council of Canada (NSERC) (ID: 461514).

ORCID iDs

Omer Faruk Ugurlu

Mustafa Kumral

Author biographies

Omer Faruk Ugurlu has a PhD in mining engineering from McGill University. He currently works as a university lecturer at Istanbul University-Cerrahpasa. His research focuses on statistics, paste fill system design, surface mining, equipment reliability, mine risk analysis, and occupational health and safety in mines.

Mustafa Kumral has a PhD in mining engineering from the University of Leeds. He currently works as an associate professor at McGill University. He has over 25 years’ experience in mining research and education. His research areas are mineral economics, mine reliability, and asset management. He has authored many journal and conference papers on these topics.

References

Eren

Real-time-optimization of drilling parameters during drilling operations. Doctoral Dissertation. Middle East Technical University, Turkey, 2010, pp.1–14.

Hall

Daneshmend

LK.

Reliability modelling of surface mining equipment: data gathering and analysis methodologies. Int J Surf Mining Reclamat Environ 2003; 17: 139–155.

Rahimdel

Ataei

Kakaei

, et al. Reliability analysis of drilling operation in open pit mines. Arch Mining Sci 2013; 58: 569-578.

Peralta

Sasmito

Kumral

Reliability effect on energy consumption and greenhouse gas emissions of mining hauling fleet towards sustainable mining. J Sustain Mining 2016; 15: 85–94.

Ebeling

CE.

An introduction to reliability and maintainability engineering. Long Grove, Illinois: Tata McGraw-Hill Education, 2004, pp.2-30.

Elsayed

EA.

Reliability engineering. Hoboken, New Jersey, USA: John Wiley & Sons, 2012, pp.1-73.

Wang

Coit

DW.

Repairable systems reliability trend tests and evaluation. In: Butler

(ed) reliability and maintainability symposium, IEEE proceedings, Alexandria, VA, USA, 24 January 2005, Paper no. 11D4, pp.416-421. Danvers, Massachusetts: IEEE.

Rahimdel

Ataei

Khalokakaei

, et al. Reliability-based maintenance scheduling of hydraulic system of rotary drilling machines. Int J Mining Sci Technol 2013; 23: 771-775.

Ozdemir

Kumral

A system-wide approach to minimize the operational cost of bench production in open-cast mining operations. Int J Coal Sci Technol 2019; 6: 84–94.

10.

Yarmuch

Epstein

Cancino

, et al. Evaluating crusher system location in an open pit mine using Markov chains. Int J Mining Reclamat Environ 2017; 31: 24-37.

11.

Geyer

Introduction to Markov Chain Monte Carlo. Handbook Markov Chain Monte Carlo 2011; 10: 3-48.

12.

Pilipovic

Energy risk: valuing and managing energy derivatives. New York: McGraw Hill Professional, 2007, pp.60-63.

13.

Ozdemir

Kumral

Stochastic assessment of the material haulage efficiency in the earthmoving industry. J Constr Eng Manag 2017; 143: 05017013.

14.

Ataei

Kakaie

Ghavidel

, et al. Drilling rate prediction of an open pit mine using rock mass drillability index. Int J Rock Mech Mining Sci 2014; 73: 130–138.

15.

Basarir

Tutluoglu

Karpuz

Penetration rate prediction for diamond bit drilling by adaptive neuro-fuzzy inference system and multiple regressions. Eng Geol 2014; 15: 1–9.

16.

Al-Chalabi

Hoiseinie

Lundberg

. Monte Carlo reliability simulation of underground mining drilling rig. In: Current trends in reliability, availability, maintainability and safety. Switzerland: Springer, 2016, pp.633-643.

17.

Crow

LH.

Poisson process. In: F Proschan RJ Serfling (eds) Reliability analysis for complex, repairable systems in reliability and biometry. Philadelphia, PA: SIAM, 1975, pp.379-410.

18.

Birolini

Reliability engineering: theory and practice. Berlin, Heidelberg: Springer, 2010, pp.458-488.

19.

Gelman

Rubin

DB.

Inference from iterative simulation using multiple sequences. Stat Sci 1992; 7: 457-472.

20.

Hillier

Lieberman

GJ.

Markov chains. In: Kane

(ed) Advance praise for introduction to operation research. 7th ed. New York: McGraw-Hill Higher Education, 2001, pp.802-833.

21.

Jonannes

Nicholas

. MCMC Methods for Continuous-Time Financial Econometrics. Handbook of Financial Econometrics: Applications 2010; 2: 1–71.

22.

Andrieu

De Freitas

Doucet

, et al. An introduction to MCMC for machine learning. Mach Learn 2003; 50: 5-43.

23.

Kumral

Sari

YA.

Simulation-based mine extraction sequencing with chance constrained risk tolerance. Simulation 2017; 93: 527–539.

24.

Ugurlu

Kumral

. Optimization of drill bit replacement time in open-cast coal mines. Int J Coal Sci Technol 2019; 6: 399–407.

25.

Pang

Forster

Troutt

MD.

Estimation of wind speed distribution using Markov chain Monte Carlo techniques. J Appl Meteorol 2001; 40: 1476–1484.

26.

Malhotra

Markov Chain Monte Carlo Models, Gibbs Sampling, & Metropolis algorithm for high-dimensionality complex stochastic problems: applications in network and computer security. computational quantitative finance-risk management. Ithaca, NY: Global Risk Management Network, LLC, 2014.

27.

Mardia

Nyirongo

Walder

, et al. Markov chain Monte Carlo implementation of rock fracture modelling. Math Geol 2007; 39: 355-381.

28.

Hong

Sarkar

Commodity betas with mean reverting output prices. J Bank Fin 2008; 32: 1286–1296.

29.

Carcano

Falbo

Stefani

Speculative trading in mean reverting markets. Eur J Oper Res 2005; 163: 132–144.

30.

Exley

Mehta

Smith

Mean reversion. In: finance and investment conference, Brussels, 28 June 2004, pp.1-31.

31.

Swishchuk

AN.

Explicit option pricing formula for a mean-reverting asset in energy market. J Numer Appl Math 2008; 1: 216-233.

32.

Triantafyllopoulos

Montana

Dynamic modeling of mean-reverting spreads for statistical arbitrage. Computat Manag Sci 2011; 8: 23-49.

33.

Zhang

Zhuang

Yin

Stochastic algorithms and numerics for mean-reverting asset trading. Stochast Anal Stochast Syst Appl Fin 2011; 10: 245-261.

34.

Ugurlu

Kumral

. Cost optimization of drilling operations in open-pit mines through parameter tuning. Qual Technol Quantit Manag 2020; 17: 173–185.