# Operational Efficiency Forecasting Model of an Existing Underground Mine Using Grey System Theory and Stochastic Diffusion Processes

Authors: Svetlana Strbac Savic; Jasmina Nedeljkovic Ostojic; Zoran Gligoric; Cedomir Cvijovic; Snezana Aleksandrovic
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

Forecasting the operational efficiency of an existing underground mine plays an important role in strategic planning of production. Degree of Operating Leverage (DOL) is used to express the operational efficiency of production. The forecasting model should be able to involve common time horizon, taking the characteristics of the input variables that directly affect the value of DOL. Changes in the magnitude of any input variable change the value of DOL. To establish the relationship describing the way of changing we applied multivariable grey modeling. Established time sequence multivariable response formula is also used to forecast the future values of operating leverage. Operational efficiency of production is often associated with diverse sources of uncertainties. Incorporation of these uncertainties into multivariable forecasting model enables mining company to survive in today’s competitive environment. Simulation of mean reversion process and geometric Brownian motion is used to describe the stochastic diffusion nature of metal price, as a key element of revenues, and production costs, respectively. By simulating a forecasting model, we imitate its action in order to measure its response to different inputs. The final result of simulation process is the expected value of DOL for every year of defined time horizon.


---

## Body

## 1. Introduction

Efficiency is a prerequisite for the survival of every mining company, especially in high competitive market environment such as mineral resource industry. Efficiency signifies company’s ability to meet its short or long-term goals. Operational efficiency is defined as the ratio between the input required to run production process and the output gained from the production. In the context of mining business, operational efficiency refers to the length of time until mineral assets are transformed to money. Peak operational efficiency occurs when the right combination of mineral deposit characteristics, human resources, mining technology, and mineral processing come together to optimize mining performance.Operational efficiency is related to finding the very best way of mining to produce a mineral asset. It enables management of the mining company to increase productivity, increase profitability, improve competitiveness, use freed-up capacity, and enable company to grow or increase its future market value.Mining company’s management uses different strategies to reach their defined goals. One of the most important elements of a company’s management operations is to forecast what goals are realistic and capability estimation of the company in order to achieve them. Planners try to forecast the behavior of the input variables of production environment and arrive at desirable states. They create their strategies on realistic targets drawn from these forecasts. A forecast is based on past and current production indicators, that is, business numbers. In many cases, the forecast might be burdened by some inaccuracies, so it would be a mistake to base a budget on that. It is very important to emphasize the forecasting acts serve as a basis for further planning.The forecasting process is performed in a specific environment. If we take into consideration the fact that the environment is changed over time then it is obvious that the forecasts and targets are changed as well. Management must be able to describe environment changes in order to strategically link the forecasting and planning functions, improving the performance of both.Large capital intensive projects, such as those in the mineral resource industry, are often associated with diverse sources of both endogenous and exogenous uncertainties. These uncertainties can greatly influence the operational efficiency. Having the ability to plan for these uncertainties is increasingly recognized as critical to long-term mining company success. In the mining industry in particular, the relationships between input variables that are controllable and those that are not and the physical and economic outcomes are complex and often nonlinear. Operational efficiency forecasting of mine in today’s environment is much complex than it was just a few years ago. There are typically many variables, which are directly or indirectly associated with the forecasting process.There is a considerable literature dedicated to the operational efficiency measurement. It includes many approaches, which take into account various aspects of the problem.Briciu et al. [1] applied the concept of Cost-Volume-Profit analysis in monitoring and measuring the performance of companies in the mining industry in Romania. Zhao et al. [2] evaluated operating efficiency of Chinese Coal Mining Companies by Malmqusit Productivity Index. Factor analysis was employed by Li et al. to evaluate performance of coal mine companies [3]. Four types of empirical tests have been performed (parametric mean differences, nonparametric Wilcoxon rank sum test, static regression panel estimation, and dynamic regression panel estimation) to estimate managerial and operational efficiency of privatized mining companies in Jordan [4]. A stochastic frontier analysis method was used to estimate profit efficiency in the South African mining sector [5]. The analytic hierarchy process (AHP) methodology was selected for ranking the efficiency of selected platinum mining methods [6]. Improving the efficiency of truck/shovel operations can increase the total operational efficiency of mine. To make the truck and shovel fleet more optimum on the basis of productivity improvements and associated costs reduction, it was essential to calculate the best match factors and truck/shovels assignments (or fleet size) for the system [7].This paper investigates the relationship between ore production rate, fixed costs, revenues, production costs, working days, and degree of use of production capacity, as input variables, and Degree of Operating Leverage (DOL) as indicator of operational efficiency of an existing underground mine. Multivariable grey model is used to establish the time sequence response formula of DOL. This formula enables us to predict the future values of DOL time series based only on a set of the most recent data. To decrease uncertainty of future values of the operational efficiency we use expert’s knowledge and simulation processes to find future values of input variables affecting them. Estimation of future revenues is based on the application of mean reversion process, normal and uniform distribution. Geometric Brownian motion is used to define future values of production costs. Values of ore production rate, fixed costs, working days, and degree of use of production capacity are estimated by expert’s knowledge.By simulating a forecasting system, we imitate its action in order to measure its response (output) to different inputs. The simulation allows analysts to describe the uncertainty of variables that influence the value of DOL by different time depending scenarios. The first objective of the use of simulation in the forecasting is to determine the distribution of the DOL from the variables that affect his performance, which results in the average or expected value of DOL for every year of defined time horizon.


## 2. Grey Information Systems and Multivariable Grey ModelGM(h,N)

The grey theory established by Dr. Deng includes grey relational analysis, grey modeling, prediction, and decision making of a system in which the model is unsure or the information is incomplete [8]. It provides an efficient solution to the uncertainty, multi-input, and discrete data problem.Grey models predict the future values of a time series based only on a set of the most recent data depending on the window size of the predictor. It is assumed that all data values to be used in grey models are positive, and the sampling frequency of the time series is fixed. From the simplest point of view, grey models can be viewed as curve fitting approaches [9].In grey system theory, multivariable grey model (GM(h,N)) denotes a grey model, where h is the order of the differential equation and N is the number of variables. The GM(h,N) is defined as follows [10]: (1)∑i=0haidix11kdti=∑j=2Nbjxj1k,k=1,2,…,n,where ai and bj are determined coefficients, x11k is the major sequence factor, xj1k is the influencing sequence factors, and k is the time sequence variable.In this paper we are applyingGM(1,N) [11–15]. First we separate the data into two sequences: major sequence factor (it masters the system behaviors) and influencing sequence factors (it influences the system behaviors).The major sequence factor is represented as follows:(2)X10=x101,x102,…,x10k,k=1,2,…,n.The influence sequence factors are (3)X20=x201,x202,…,x20k,X30=x301,x302,…,x30k,⋮XN0=xN01,xN02,…,xN0k.Assume the original series of data with k equal time samples is represented as X0k=x01,x02,…,x0k, k=1,2,…,n. The first-order Accumulating Generation Operator (AGO) of X(0) is defined as X1k=x11,x12,…,x1k, where x1k=∑k=1nx0k for k=1,2,…,n. This operator is used to smooth the randomness of the data and to weaken the tendency of variation. The average generation of adjacent sample sequence of X(1)(k) is as follows: Z1k=z11,z12,…,z1k, where z1k=1/2x1k+x1k-1, for k=1,2,…,n.The first-order grey differential equation of theGM(1,N) is (4)dx11kdk+ax11k=∑j=2Nbjxj1k.According to the GM(1,N), the constructed AGO sequence is (5)x10k+az11k=∑j=2Nbjzj1k.The constructed AGO sequence can be represented by the following matrix form: (6)x102x103⋮x10n=-z112x212⋯xN12-z113x213⋯xN13⋮⋮⋮-z11nx21n⋯xN1nab2⋮bN.Factors of the matrix a^=a,b2,…,bNT are obtained by using the least square method as follows: (7)a^=BTB-1BTY,where(8)Y=x102x103⋮x10n,B=-z112x212⋯xN12-z113x213⋯xN13⋮⋮⋮-z11nx21n⋯xN1n.Finally, the GM(1,N) can be expressed as follows: (9)dx11dk+ax11=b2x21+b3x31+⋯+bNxN1.The time sequence response formula, that is, the solution of (9), is(10)x^11k+1=x101-∑j=2Nbj-1a·xj1k+1e-ak+∑j=2Nbj-1a·xj1k+1,k=0,1,2,…,n-1.According to (10), we can obtain X^1k+1=x^11k+1. The fitting forecast value of xk0 can be obtained by the inverse AGO as follows:(11)x^0k=x^1k-x^1k-1,k=2,3,…,n,x^01=x01,k=0.The predicted values when k≤n can be used to check the adequacy of the model and predicted value when k>n can be used as the forecast value for the data series.The relative percentage error (RPE(j,k)) and the average relative percentage error (ARPE(j,k)) are as follows:(12)RPE=x^j0k-xj0kxj0k×100%ARPE=1n-1∑k=2nx^j0k-xj0kxj0k×100%.Grey relational (relevance) analysis is part of grey system theory that gives us information about the relevance degree of each input variable in the system. It is utilized to ascertain the primary factors that are needed to make a superiority comparison in the system [11]. The fundamental principle is to recognize the relevance degree among many factors, according to the similarity levels of the geometrical patterns of sequence curves. We can judge whether they are close or not according to the similar degree of sequence curve shape. The more similar the curve is, the higher the correlation degree between relative series is, and vice versa [11].Set the major sequence factor variable to be analyzed,X10k, k=1,2,…,n, as the dependent variable, and the influencing sequence factors, Xj0k, j=2,…,N, as the independent variables. Because various variables have different physical significances, that is, dimensions, it is difficult to obtain accurate conclusion. To eliminate dimensions from analysis, equalization method is used, and formula is (13)xj′k=xjk∑j=1Nxjk,j=1,2,…,N.The grey relational coefficient is calculated as follows [16]:(14)ξjk=min︸jkx1′k-xj′k+θ·max︸jkx1′k-xj′kΔjk+θ·max︸jkx1′k-xj′k,j=2,3,…,N,where θ refers to the discriminating coefficient and θ∈0,1; Δjk=x1′k-xj′k.The aggregated value of the grey relational coefficient is calculated as follows:(15)ρjk=1n∑k=1nξjk,j=2,3,…,N.The higher ρjk is, the bigger impact the jth variable has on the variable to be analyzed; on the contrary, the lower ρjk is, the smaller impact the jth element has on the variable to be analyzed. By this way, we can see how contribution of the input variables to output variable is changed over time.Finally, the management of the company obtains the aggregated information vectorAIVj,k=X^0j,k,RPEj,k,ARPEj,k,ξjk,ρj,k which sublimates the following important indicators:(i)
The forecasting model.(ii)
The relative percentage error of forecasting model.(iii)
The average relative percentage error of forecasting model.(iv)
The grey relational coefficient between inputs and output for thekth time point.(v)
The aggregated grey relational coefficient between inputs and output.


## 3. Forecasting Model

### 3.1. Concept of the Model

Overall operational efficiency measures have enjoyed a great deal of interest among planners analyzing the performance of an underground mine. Suppose the production of mine in periodk can be represented by mining technology set (MTS): (16)MTSk=Ik,Ok:Ik∈R+g,Ok∈R+l,Ikcan  realizeOk,where Ik=i1k,i2k,…,igk is the input mining technology vector and Ok=o1k,o2k,…,olk is the output vector which represents the indicators obtained from production. The MTS is the space of the feasible combinations of input-output vectors in period k.The input vector of the MTS is composed of many variables but they can be divided into three main subspaces: characteristics of mineral deposit, underground mining methods, and mineral processing methods. The output vector is usually composed of the following variables: realized production rate and quality of produced mineral assets.Operational efficiency is usually expressed by some kind of economic indicator. According to required type of operational efficiency, MTS can be transformed and joined to economic set (ES). This transformation means that some variables ofI(k) and O(k) are used to calculate some variables of ES.To express the operational efficiency of an existing underground mine we apply the concept of Degree of Operating Leverage (DOL). The DOL is the leverage ratio that sums up the effects of an amount of operating leverage on the company’s earnings (revenues) before interests and taxes (EBIT).The DOL during the periodk can be calculated with following formula (Weygandt et al.) [17]: (17)DOLk=REVk-PCkREVk-PCk-FCk,where REV is revenues (USD), PC is production costs (USD), and FC is fixed costs (USD).Operating leverage refers to the amount of operating fixed costs in the cost structure. DOL is the ratio of the relative change in sales. Implications are related to amplifying the changes in sales volume into larger changes in EBIT. Even in a perfectly competitive market, price plays a role in determining the DOL magnitude. Assume that a firm is currently operating with a positive profit. When market equilibrium price rises, the denominator in (17) rises reducing the DOL with no change in output, operating fixed cost, or unit variable cost. Consequently, breakeven output falls since the contribution margin is larger. DOL can vary due to changes in any of the variables appearing in (17). These variables include management determined choices (operating fixed cost and output levels); market determined parameters (price in a competitive market as time passes); and economic and engineering realities (unit variable costs, given operating fixed cost increases due to new capital acquisition) [18].The higher value of DOL indicates that the business of mining company is exposed to the greater risk.The input vectorXk=Ak∪Bk∪Ck that is used to evaluate operational efficiency is represented as the union of three following sets:A(k)
is a subset of the input mining technology variables.B(k)
is a subset of the output mining technology variables.C(k)
is a set of the external variables.According to our problem, we create the vector space of input variables as follows:(18)Xk=Xjk=x6,k=x1,1x1,2⋯x1,kx2,1x2,2⋯x2,kx3,1x3,2⋯x3,kx4,1x4,2⋯x4,kx5,1x5,2⋯x5,kx6,1x6,2⋯x6,k,where x1,k is ore production rate (t/year), x2,k is fixed costs (USD), x3,k is revenues (USD), x4,k is production costs (USD), x5,k is working days (days/year), and x6,k is degree of use of production capacity (%).Lety1,y2,…,yk be a sequence of observed DOL values on an annual basis, that is, the major sequence factor. Equation (9) will be used for expressing the forecast model of DOL of an existing underground mine as follows: (19)dy1dk+ay1=b1x11+b2x21+⋯+bNxN1.The time sequence response formula of DOL is(20)y^1k+1=y01-∑j=16bja·xj1k+1e-ak+∑j=16bja·xj1k+1,k=0,1,2,…,n-1.The fitting forecast value of yk0 can be obtained by the inverse AGO as follows:(21)y^0k=y^1k-y^1k-1,k=2,3,…,n,y^01=y01,k=0.For k>n we can make forecasts of DOL as follows:(22)y^1n+h+1=y01-∑j=16bja·xj1n+h+1e-an+h+∑j=16bja·xj1n+h+1,h=0,1,2,…,T-1,y^0n+h=y^1n+h-y^1n+h-1,h=1,2,…,T,where T is future time horizon.The forecasting model should be able to involve a common time horizon, taking the characteristics of the input variables that directly affect the value of DOL. In this section, the proposed method of forecasting of DOL under multiple uncertainty and nonuniformity of input variables is outlined. The method is based mainly on simulating multiple realizations of the uncertain variables and making forecasts using value expectations.According to the valuation function defined by (22), the expected value of DOL over future time is defined as follows: (23)Ey^0∣x1,x2,x3,x4,x5,x6,n+h.Before we start to create forecasting process, it is necessary to divide the set of input variables into two subsets; Xn+h=EXPn+h∪SIMn+h, where EXPn+h denotes the set composed of variables defined by experts, while SIMn+h denotes the set composed of variables defined by some kind of stochastic law. In our case these subsets are(24)EXPn+h=x1,x2,x5,x6,SIMn+h=x3,x4.Values of x1, x2, x5, and x6 are of deterministic nature while the others are stochastic.


### 3.2. Volatility of the Input Variables

#### 3.2.1. Volatility of Revenues

Mining companies not having smelting facilities realize their revenues by selling their metal concentrates as a final product. Estimating revenues is, indeed, a difficult and risky activity. Annual mine revenue is calculated by multiplying the number of units produced and sold during the year by the sales price per unit. The annual mine revenue (x3) is calculated according to the following equation: (25)REVyear=Qyear·∑j=1nVjcon·Gj·Mjmjcon,where Qyear is the annual production of ore (x1), Vcon is the value of the metal concentrate, G is the grade of ore mined (%), M is the mill recovery rate (%), mcon is the metal content of the concentrate (%), and n is the number of metal concentrates derived from the ore (n>1 for polymetallic deposits).The first major component of the mine revenue calculation is annual production of concentrate. One of the key variables associated with annual concentrate production is the tonnage of ore mined. Annual ore tonnage is derived from the mining project schedule and is denotedQ.The second key variable associated with determining the annual production of saleable units is the grade of the ore mined. The concept of the ore grade (G) is defined as the ratio of useful mass of metal to the total mass of ore and its critical value fluctuates over space and can be estimated by normal probability density function. The use of normal distribution is based on geostatistical methods developed for evaluation of grades in mineral deposits. The deposit is sampled by means of drill-holes. Each sample gives us the information about mass of metal (grade) and we can create histogram of grades values. According to the obtained histogram adequate probability distribution function is selected. In most cases of mineral deposits normal distribution is used as representative function: (26)G~Nμ,σ2.At first sight, it seems that G does not change its value over time, but if we take into consideration the mining front advances over time through different parts of deposit, then it is also time depending. For simplicity, we adopt the values of G which belong to the same interval, defined by μ±2σ, over future time horizon.Most ores require beneficiation before saleable product can be produced. The resulting milling losses must be estimated and appropriate recovery percentages established. These recoveries are commonly estimated from a metallurgical testing program. Percentage recovery is the third basic variable which must be estimated to arrive at a final estimate of the annual production of saleable units extracted from the mine and is denotedM. A specific stochastic behavior that is used to quantify uncertainty related to mill recovery rate (M) is uniform probability density function. Perhaps it is virtually certain that mill recovery will lie in interval [a,b] but it is possible to achieve it by mineral processing control. For example, flotation is the most widely used method for the concentration of fine grained minerals. It takes advantage of the different physicochemical surface properties of minerals in particular their wettability, which can be a natural property or artificially changed by chemical reagents. By altering the hydrophobic (water repelling) or hydrophilic (water attracting) conditions of their surfaces, mineral particles suspended in water can be induced to adhere to air bubbles passing through a flotation cell or to remain in the pulp. The air bubbles pass to the upper surface of the pulp and form a froth, which, together with the attached hydrophobic minerals, can be removed. The tailings, containing the hydrophilic minerals, can be removed from the bottom of the cell. According to that, within interval [a,b] no value is more likely than any other. This gives rise to a uniform distribution: (27)M~unifa,b.It seems that M does not change its value over time, but if we take into consideration the fact that the flotation is influenced by many parameters and fact that it is very difficult to keep them constantly over time, then it is also time depending. For simplicity, we adopt the values of M which belong to the same interval [a,b] but M takes the stochastic values for every year of the time horizon.The second major component of the mine revenue calculation is unit sales price or unit metal concentrate sales price (Vcon). It directly depends on mineral asset price, metal content of the concentrate, and metal recovery rate. Estimating future mineral prices is an exercise for which a high error of estimation invariably exists. The characteristically long preproduction periods of mining projects mean that their success will be determined by mineral prices five to ten years in the future.The market risks related to mineral asset price (P) are modeled with a special stochastic process, a mean reversion process. The mean reversion process has economic logic; for example, although the commodity prices have sensible short-term oscillations, they tend to revert back to a “normal” long-term equilibrium level. The past values of the changes in this uncertain factor help predict the future. We will use a model where the metal spot price is assumed to follow the stochastic process [19]:(28)dP=αln⁡P--ln⁡PPdt+σPdW.Let x=ln⁡P; applying Ito’s lemma allows the characterization of the log price by an Ornstein-Uhlenbeck stochastic mean reverting process: (29)dx=αx--xdt+σdWwith (30)x¯=ln⁡P--σ22α,where P- is the long-run equilibrium metal price, α measures the speed of mean reversion to the long-run mean log price P-, dW is an increment to a standard Brownian motion, and σ refers to the price volatility rate. The metal price adjustment mechanism is accounted for by market forces. The correct discrete-time format for the continuous-time process of mean reversion is the stationary first-order autoregressive process [20], so the sample path simulation equation for xt is performed by using exact discrete-time expression:(31)xt=xt-1e-αΔt+x-1-e-αΔt+N0,1σ1-e-2αΔt2α,where Δt is the fixed time interval from time t to t+1 and N(0,1) is the normally distributed random variable.By substituting (31) to P=ex, we have exact discrete-time equation for Pt, given by (32)Pt=eln⁡Pt-1e-αΔt+ln⁡P--σ2/2α1-e-αΔt+N0,1σ1-e-2αΔt/2α.In order to estimate the parameters of the mean reversion process, we run the following regression: (33)dxt+1=β0+β1xt+ε,where β0=αx-dt and β1=-αdt. Hence, if we regress observation dx against x, we can obtain estimates of β0 and β1. σ is the standard deviation obtained from the regression. The speed of mean reversion (α) is the negative of the slope, while the long-run equilibrium P- is the intercept estimate of that regression divided by the speed of mean reversion.LetP=Pt,t=0,1,…,T denote a price scenario with spot prices Pt, where Pt is determined by (32). Figure 1 presents a sample path of the metal price (e.g., lead) simulated using the above equation.Figure 1
One simulated lead price path on a yearly time resolution.In the purpose ofP simulation, we apply a metal price scenario for the time interval, n,n+T, with increment Δt=1.The unit value of metal concentrate sales price can be expressed as follows:(34)Vcont=fPt,mcon,mmr,where mmr is the metal recovery rate (%). The value of Vcon is calculated in different ways with respect to type of metal.


#### 3.2.2. Volatility of Production Costs

Production costs (PC) are incurred directly in the production process. These costs include the ore and waste development of individual stopes, the actual stoping activities, the mine services providing logistical support to the miners, and the milling and processing of the ore at the plant. The uncertainties related to the future states of unit production costs are modeled with a special stochastic process, the geometric Brownian motion. Certain stochastic processes are functions of a Brownian motion process and these have many applications in finance, engineering, and the sciences. Some special processes are solutions of Itô-Doob type stochastic differential equations (Ladde, Sambandham) [21].In this model, we apply a continuous time process using the Itô-Doob type stochastic differential equation to describe movement of unit production costs. A general stochastic differential equation takes the following form:(35)dCOt=ρ·COt,tdt+σ·COt,t·dWt,COt0=CO0.Here, t≥t0, Wt is a Brownian motion, and COt>0; this is the cost process.C
O
t is called the geometric Brownian motion, which is solution of the following linear Itô-Doob type stochastic differential equation: (36)dCOt=ρ·COtdt+σ·COt·dWt,where ρ is the drift, σ is the volatility, and Wt is normalized Brownian motion.Using the Itô-Doob formula applied toCOt=ln⁡COt, we can solve this equation: (37)COt=COt-1·eρ-σ2/2Δt+N0,1σΔt.Equation (37) describes an operating cost scenario with spot costs COt.LetCO=COt,t=0,1,…,T denote a cost scenario with spot costs COt, where COt is determined by (37). Figure 2 presents a sample paths (s=1,2,…,S) of the unit production cost simulated using (37) S times.Figure 2
Simulated cost paths on a yearly time resolution.In the purpose of CO simulation, we apply a cost scenario for the time interval,n,n+T, with increment Δt=1.The annual production costs are expressed as follows:(38)PCyear=Qyear·COt.


### 3.3. Variables Based on Expert’s Knowledge

To estimate an adequate value of these variables, it is necessary to get opinions of experts dealing with them. Suppose we havep experts and each of them has given his/her opinion. The final value is expressed by averaging the opinions of experts. Table 1 presents the expert estimation process.Table 1
Expert estimation of the variables.
Expert
Value
E
1
E
2
⋯
E
p
x
1,1
x
1,2
⋯
x
1
,
p
E
(
x
1
)
=
x
1,1
+
x
1,2
+
⋯
+
x
1
,
p
p
x
2,1
x
2,2
⋯
x
2
,
p
E
(
x
2
)
=
x
2,1
+
x
2,2
+
⋯
+
x
2
,
p
p
x
5,1
x
5,2
⋯
x
5
,
p
E
(
x
5
)
=
x
5,1
+
x
5,2
+
⋯
+
x
5
,
p
p
x
6,1
x
6,2
⋯
x
6
,
p
E
(
x
6
)
=
x
6,1
+
x
6,2
+
⋯
+
x
6
,
p
p


### 3.4. DOL Forecasting Model Based on Expert’s Knowledge and Simulation

The model is developed on the basis of expert’s knowledge and simulation of DOL changing over defined time horizon and takes into account the variability of input parameters. By simulating a forecasting system, we imitate its action in order to measure its response (output) to different inputs. The advantage of simulating a system is the possibility of replicating its evolution as many times as necessary in independent conditions. The simulation allows analysts to describe the uncertainty of variables that influence the value of DOL by different time depending scenarios. The values of DOL are usually forecasted on an annual basis.The first objective of the use of simulation in the forecasting is to determine the distribution of the DOL from the variables that affect his performance, which results in the average or expected value of DOL for every year of defined time horizon. The relation between uncertain variables affecting the value of DOL is described by (22).For each simulation, the input values and DOL result represent one possible state of nature. Simulated values of DOL are obtained by performing the following calculations:(39)y^1sn+h+1=y^1x1,n+h,x2,n+hs,x3,n+hs,x4,n+hs,x5,n+h,x6,n+hs,h=0,1,2,…,T-1,y^0sn+h=y^1sn+h-y^1sn+h-1,h=1,2,…,T,s=1,2,…,S,where S denotes the number of simulations.Space of simulation fors=1 and k>n, where the first six rows concern the evolution path of the input variables while the last one concerns evolution path of the DOL, can be represented as follows:(40)DOLTs=1=s=1s=1⋯s=1x1,1x1,2⋯x1,Tx2,1x2,2⋯x2,Tx3,1x3,2⋯x3,Tx4,1x4,2⋯x4,Tx5,1x5,2⋯x5,Tx6,1x6,2⋯x6,T↓↓↓↓equation22equation22equation22equation22↓↓↓↓y^11s=1y^21s=1⋯y^T1s=1↓↓↓↓y^10s=1y^20s=1⋯y^T0s=1.Once the simulations DOLTs, s=1,2,…,S, have been obtained, they can be used to estimate the distribution of the DOL for every year of the defined time horizon. The expected values of the DOL can be expressed by the following vector space: (41)EDOLTS=y^10s=1y^20s=1⋯y^T0s=1y^10s=2y^20s=2⋯y^T0s=2⋮⋮⋱⋮y^10s=Sy^20s=S⋯y^T0s=S↓↓↓↓Ey^10Ey^20⋯Ey^T0.According to (41), the management of the company obtains the forecasting information vector of DOL: (42)FIVj,T=EDOL=EDOL1EDOL2⋯EDOLT.This vector concerns the information of the possible future states of the Degree of Operating Leverage over defined time horizon.Space of simulation fors=1 and k>n, where the first seven rows concern the evolution paths of the DOL and input variables while the rest concern evolution paths of the grey relational coefficient, can be represented as follows: (43)ξj,Ts=1=s=1s=1⋯s=1y^10y^20⋯y^T0x1,1x1,2⋯x1,Tx2,1x2,2⋯x2,Tx3,1x3,2⋯x3,Tx4,1x4,2⋯x4,Tx5,1x5,2⋯x5,Tx6,1x6,2⋯x6,T↓↓↓↓equation14equation14equation14equation14↓↓↓↓ξ1,1s=1ξ1,2s=1⋯ξ1,Ts=1⋮⋮⋱⋮ξ6,1s=1ξ6,2s=1⋯ξ6,Ts=1.Once the simulations ξj,Ts, s=1,2,…,S, have been obtained, they can be used to estimate the distribution of the ξj for every year of the defined time horizon. The expected values of the ξj can be expressed by the following vector space: (44)Eξj,TS=Eξ1,1SEξ1,2S⋯Eξ1,TSEξ2,1SEξ2,2S⋯Eξ2,TSEξ3,1SEξ3,2S⋯Eξ3,TSEξ4,1SEξ4,2S⋯Eξ4,TSEξ5,1SEξ5,2S⋯Eξ5,TSEξ6,1SEξ6,2S⋯Eξ6,TS.The expected values of the aggregated grey relational coefficient are calculated as follows: (45)ρj=1h∑h=1TEξj,hS,j=1,2,…,6.According to (45), the management of the company obtains the forecasting information vector of the aggregated grey relational coefficient:(46)ρj=ρ1ρ2⋮ρ6.


## 3.1. Concept of the Model

Overall operational efficiency measures have enjoyed a great deal of interest among planners analyzing the performance of an underground mine. Suppose the production of mine in periodk can be represented by mining technology set (MTS): (16)MTSk=Ik,Ok:Ik∈R+g,Ok∈R+l,Ikcan  realizeOk,where Ik=i1k,i2k,…,igk is the input mining technology vector and Ok=o1k,o2k,…,olk is the output vector which represents the indicators obtained from production. The MTS is the space of the feasible combinations of input-output vectors in period k.The input vector of the MTS is composed of many variables but they can be divided into three main subspaces: characteristics of mineral deposit, underground mining methods, and mineral processing methods. The output vector is usually composed of the following variables: realized production rate and quality of produced mineral assets.Operational efficiency is usually expressed by some kind of economic indicator. According to required type of operational efficiency, MTS can be transformed and joined to economic set (ES). This transformation means that some variables ofI(k) and O(k) are used to calculate some variables of ES.To express the operational efficiency of an existing underground mine we apply the concept of Degree of Operating Leverage (DOL). The DOL is the leverage ratio that sums up the effects of an amount of operating leverage on the company’s earnings (revenues) before interests and taxes (EBIT).The DOL during the periodk can be calculated with following formula (Weygandt et al.) [17]: (17)DOLk=REVk-PCkREVk-PCk-FCk,where REV is revenues (USD), PC is production costs (USD), and FC is fixed costs (USD).Operating leverage refers to the amount of operating fixed costs in the cost structure. DOL is the ratio of the relative change in sales. Implications are related to amplifying the changes in sales volume into larger changes in EBIT. Even in a perfectly competitive market, price plays a role in determining the DOL magnitude. Assume that a firm is currently operating with a positive profit. When market equilibrium price rises, the denominator in (17) rises reducing the DOL with no change in output, operating fixed cost, or unit variable cost. Consequently, breakeven output falls since the contribution margin is larger. DOL can vary due to changes in any of the variables appearing in (17). These variables include management determined choices (operating fixed cost and output levels); market determined parameters (price in a competitive market as time passes); and economic and engineering realities (unit variable costs, given operating fixed cost increases due to new capital acquisition) [18].The higher value of DOL indicates that the business of mining company is exposed to the greater risk.The input vectorXk=Ak∪Bk∪Ck that is used to evaluate operational efficiency is represented as the union of three following sets:A(k)
is a subset of the input mining technology variables.B(k)
is a subset of the output mining technology variables.C(k)
is a set of the external variables.According to our problem, we create the vector space of input variables as follows:(18)Xk=Xjk=x6,k=x1,1x1,2⋯x1,kx2,1x2,2⋯x2,kx3,1x3,2⋯x3,kx4,1x4,2⋯x4,kx5,1x5,2⋯x5,kx6,1x6,2⋯x6,k,where x1,k is ore production rate (t/year), x2,k is fixed costs (USD), x3,k is revenues (USD), x4,k is production costs (USD), x5,k is working days (days/year), and x6,k is degree of use of production capacity (%).Lety1,y2,…,yk be a sequence of observed DOL values on an annual basis, that is, the major sequence factor. Equation (9) will be used for expressing the forecast model of DOL of an existing underground mine as follows: (19)dy1dk+ay1=b1x11+b2x21+⋯+bNxN1.The time sequence response formula of DOL is(20)y^1k+1=y01-∑j=16bja·xj1k+1e-ak+∑j=16bja·xj1k+1,k=0,1,2,…,n-1.The fitting forecast value of yk0 can be obtained by the inverse AGO as follows:(21)y^0k=y^1k-y^1k-1,k=2,3,…,n,y^01=y01,k=0.For k>n we can make forecasts of DOL as follows:(22)y^1n+h+1=y01-∑j=16bja·xj1n+h+1e-an+h+∑j=16bja·xj1n+h+1,h=0,1,2,…,T-1,y^0n+h=y^1n+h-y^1n+h-1,h=1,2,…,T,where T is future time horizon.The forecasting model should be able to involve a common time horizon, taking the characteristics of the input variables that directly affect the value of DOL. In this section, the proposed method of forecasting of DOL under multiple uncertainty and nonuniformity of input variables is outlined. The method is based mainly on simulating multiple realizations of the uncertain variables and making forecasts using value expectations.According to the valuation function defined by (22), the expected value of DOL over future time is defined as follows: (23)Ey^0∣x1,x2,x3,x4,x5,x6,n+h.Before we start to create forecasting process, it is necessary to divide the set of input variables into two subsets; Xn+h=EXPn+h∪SIMn+h, where EXPn+h denotes the set composed of variables defined by experts, while SIMn+h denotes the set composed of variables defined by some kind of stochastic law. In our case these subsets are(24)EXPn+h=x1,x2,x5,x6,SIMn+h=x3,x4.Values of x1, x2, x5, and x6 are of deterministic nature while the others are stochastic.


## 3.2. Volatility of the Input Variables

### 3.2.1. Volatility of Revenues

Mining companies not having smelting facilities realize their revenues by selling their metal concentrates as a final product. Estimating revenues is, indeed, a difficult and risky activity. Annual mine revenue is calculated by multiplying the number of units produced and sold during the year by the sales price per unit. The annual mine revenue (x3) is calculated according to the following equation: (25)REVyear=Qyear·∑j=1nVjcon·Gj·Mjmjcon,where Qyear is the annual production of ore (x1), Vcon is the value of the metal concentrate, G is the grade of ore mined (%), M is the mill recovery rate (%), mcon is the metal content of the concentrate (%), and n is the number of metal concentrates derived from the ore (n>1 for polymetallic deposits).The first major component of the mine revenue calculation is annual production of concentrate. One of the key variables associated with annual concentrate production is the tonnage of ore mined. Annual ore tonnage is derived from the mining project schedule and is denotedQ.The second key variable associated with determining the annual production of saleable units is the grade of the ore mined. The concept of the ore grade (G) is defined as the ratio of useful mass of metal to the total mass of ore and its critical value fluctuates over space and can be estimated by normal probability density function. The use of normal distribution is based on geostatistical methods developed for evaluation of grades in mineral deposits. The deposit is sampled by means of drill-holes. Each sample gives us the information about mass of metal (grade) and we can create histogram of grades values. According to the obtained histogram adequate probability distribution function is selected. In most cases of mineral deposits normal distribution is used as representative function: (26)G~Nμ,σ2.At first sight, it seems that G does not change its value over time, but if we take into consideration the mining front advances over time through different parts of deposit, then it is also time depending. For simplicity, we adopt the values of G which belong to the same interval, defined by μ±2σ, over future time horizon.Most ores require beneficiation before saleable product can be produced. The resulting milling losses must be estimated and appropriate recovery percentages established. These recoveries are commonly estimated from a metallurgical testing program. Percentage recovery is the third basic variable which must be estimated to arrive at a final estimate of the annual production of saleable units extracted from the mine and is denotedM. A specific stochastic behavior that is used to quantify uncertainty related to mill recovery rate (M) is uniform probability density function. Perhaps it is virtually certain that mill recovery will lie in interval [a,b] but it is possible to achieve it by mineral processing control. For example, flotation is the most widely used method for the concentration of fine grained minerals. It takes advantage of the different physicochemical surface properties of minerals in particular their wettability, which can be a natural property or artificially changed by chemical reagents. By altering the hydrophobic (water repelling) or hydrophilic (water attracting) conditions of their surfaces, mineral particles suspended in water can be induced to adhere to air bubbles passing through a flotation cell or to remain in the pulp. The air bubbles pass to the upper surface of the pulp and form a froth, which, together with the attached hydrophobic minerals, can be removed. The tailings, containing the hydrophilic minerals, can be removed from the bottom of the cell. According to that, within interval [a,b] no value is more likely than any other. This gives rise to a uniform distribution: (27)M~unifa,b.It seems that M does not change its value over time, but if we take into consideration the fact that the flotation is influenced by many parameters and fact that it is very difficult to keep them constantly over time, then it is also time depending. For simplicity, we adopt the values of M which belong to the same interval [a,b] but M takes the stochastic values for every year of the time horizon.The second major component of the mine revenue calculation is unit sales price or unit metal concentrate sales price (Vcon). It directly depends on mineral asset price, metal content of the concentrate, and metal recovery rate. Estimating future mineral prices is an exercise for which a high error of estimation invariably exists. The characteristically long preproduction periods of mining projects mean that their success will be determined by mineral prices five to ten years in the future.The market risks related to mineral asset price (P) are modeled with a special stochastic process, a mean reversion process. The mean reversion process has economic logic; for example, although the commodity prices have sensible short-term oscillations, they tend to revert back to a “normal” long-term equilibrium level. The past values of the changes in this uncertain factor help predict the future. We will use a model where the metal spot price is assumed to follow the stochastic process [19]:(28)dP=αln⁡P--ln⁡PPdt+σPdW.Let x=ln⁡P; applying Ito’s lemma allows the characterization of the log price by an Ornstein-Uhlenbeck stochastic mean reverting process: (29)dx=αx--xdt+σdWwith (30)x¯=ln⁡P--σ22α,where P- is the long-run equilibrium metal price, α measures the speed of mean reversion to the long-run mean log price P-, dW is an increment to a standard Brownian motion, and σ refers to the price volatility rate. The metal price adjustment mechanism is accounted for by market forces. The correct discrete-time format for the continuous-time process of mean reversion is the stationary first-order autoregressive process [20], so the sample path simulation equation for xt is performed by using exact discrete-time expression:(31)xt=xt-1e-αΔt+x-1-e-αΔt+N0,1σ1-e-2αΔt2α,where Δt is the fixed time interval from time t to t+1 and N(0,1) is the normally distributed random variable.By substituting (31) to P=ex, we have exact discrete-time equation for Pt, given by (32)Pt=eln⁡Pt-1e-αΔt+ln⁡P--σ2/2α1-e-αΔt+N0,1σ1-e-2αΔt/2α.In order to estimate the parameters of the mean reversion process, we run the following regression: (33)dxt+1=β0+β1xt+ε,where β0=αx-dt and β1=-αdt. Hence, if we regress observation dx against x, we can obtain estimates of β0 and β1. σ is the standard deviation obtained from the regression. The speed of mean reversion (α) is the negative of the slope, while the long-run equilibrium P- is the intercept estimate of that regression divided by the speed of mean reversion.LetP=Pt,t=0,1,…,T denote a price scenario with spot prices Pt, where Pt is determined by (32). Figure 1 presents a sample path of the metal price (e.g., lead) simulated using the above equation.Figure 1
One simulated lead price path on a yearly time resolution.In the purpose ofP simulation, we apply a metal price scenario for the time interval, n,n+T, with increment Δt=1.The unit value of metal concentrate sales price can be expressed as follows:(34)Vcont=fPt,mcon,mmr,where mmr is the metal recovery rate (%). The value of Vcon is calculated in different ways with respect to type of metal.


### 3.2.2. Volatility of Production Costs

Production costs (PC) are incurred directly in the production process. These costs include the ore and waste development of individual stopes, the actual stoping activities, the mine services providing logistical support to the miners, and the milling and processing of the ore at the plant. The uncertainties related to the future states of unit production costs are modeled with a special stochastic process, the geometric Brownian motion. Certain stochastic processes are functions of a Brownian motion process and these have many applications in finance, engineering, and the sciences. Some special processes are solutions of Itô-Doob type stochastic differential equations (Ladde, Sambandham) [21].In this model, we apply a continuous time process using the Itô-Doob type stochastic differential equation to describe movement of unit production costs. A general stochastic differential equation takes the following form:(35)dCOt=ρ·COt,tdt+σ·COt,t·dWt,COt0=CO0.Here, t≥t0, Wt is a Brownian motion, and COt>0; this is the cost process.C
O
t is called the geometric Brownian motion, which is solution of the following linear Itô-Doob type stochastic differential equation: (36)dCOt=ρ·COtdt+σ·COt·dWt,where ρ is the drift, σ is the volatility, and Wt is normalized Brownian motion.Using the Itô-Doob formula applied toCOt=ln⁡COt, we can solve this equation: (37)COt=COt-1·eρ-σ2/2Δt+N0,1σΔt.Equation (37) describes an operating cost scenario with spot costs COt.LetCO=COt,t=0,1,…,T denote a cost scenario with spot costs COt, where COt is determined by (37). Figure 2 presents a sample paths (s=1,2,…,S) of the unit production cost simulated using (37) S times.Figure 2
Simulated cost paths on a yearly time resolution.In the purpose of CO simulation, we apply a cost scenario for the time interval,n,n+T, with increment Δt=1.The annual production costs are expressed as follows:(38)PCyear=Qyear·COt.


## 3.2.1. Volatility of Revenues

Mining companies not having smelting facilities realize their revenues by selling their metal concentrates as a final product. Estimating revenues is, indeed, a difficult and risky activity. Annual mine revenue is calculated by multiplying the number of units produced and sold during the year by the sales price per unit. The annual mine revenue (x3) is calculated according to the following equation: (25)REVyear=Qyear·∑j=1nVjcon·Gj·Mjmjcon,where Qyear is the annual production of ore (x1), Vcon is the value of the metal concentrate, G is the grade of ore mined (%), M is the mill recovery rate (%), mcon is the metal content of the concentrate (%), and n is the number of metal concentrates derived from the ore (n>1 for polymetallic deposits).The first major component of the mine revenue calculation is annual production of concentrate. One of the key variables associated with annual concentrate production is the tonnage of ore mined. Annual ore tonnage is derived from the mining project schedule and is denotedQ.The second key variable associated with determining the annual production of saleable units is the grade of the ore mined. The concept of the ore grade (G) is defined as the ratio of useful mass of metal to the total mass of ore and its critical value fluctuates over space and can be estimated by normal probability density function. The use of normal distribution is based on geostatistical methods developed for evaluation of grades in mineral deposits. The deposit is sampled by means of drill-holes. Each sample gives us the information about mass of metal (grade) and we can create histogram of grades values. According to the obtained histogram adequate probability distribution function is selected. In most cases of mineral deposits normal distribution is used as representative function: (26)G~Nμ,σ2.At first sight, it seems that G does not change its value over time, but if we take into consideration the mining front advances over time through different parts of deposit, then it is also time depending. For simplicity, we adopt the values of G which belong to the same interval, defined by μ±2σ, over future time horizon.Most ores require beneficiation before saleable product can be produced. The resulting milling losses must be estimated and appropriate recovery percentages established. These recoveries are commonly estimated from a metallurgical testing program. Percentage recovery is the third basic variable which must be estimated to arrive at a final estimate of the annual production of saleable units extracted from the mine and is denotedM. A specific stochastic behavior that is used to quantify uncertainty related to mill recovery rate (M) is uniform probability density function. Perhaps it is virtually certain that mill recovery will lie in interval [a,b] but it is possible to achieve it by mineral processing control. For example, flotation is the most widely used method for the concentration of fine grained minerals. It takes advantage of the different physicochemical surface properties of minerals in particular their wettability, which can be a natural property or artificially changed by chemical reagents. By altering the hydrophobic (water repelling) or hydrophilic (water attracting) conditions of their surfaces, mineral particles suspended in water can be induced to adhere to air bubbles passing through a flotation cell or to remain in the pulp. The air bubbles pass to the upper surface of the pulp and form a froth, which, together with the attached hydrophobic minerals, can be removed. The tailings, containing the hydrophilic minerals, can be removed from the bottom of the cell. According to that, within interval [a,b] no value is more likely than any other. This gives rise to a uniform distribution: (27)M~unifa,b.It seems that M does not change its value over time, but if we take into consideration the fact that the flotation is influenced by many parameters and fact that it is very difficult to keep them constantly over time, then it is also time depending. For simplicity, we adopt the values of M which belong to the same interval [a,b] but M takes the stochastic values for every year of the time horizon.The second major component of the mine revenue calculation is unit sales price or unit metal concentrate sales price (Vcon). It directly depends on mineral asset price, metal content of the concentrate, and metal recovery rate. Estimating future mineral prices is an exercise for which a high error of estimation invariably exists. The characteristically long preproduction periods of mining projects mean that their success will be determined by mineral prices five to ten years in the future.The market risks related to mineral asset price (P) are modeled with a special stochastic process, a mean reversion process. The mean reversion process has economic logic; for example, although the commodity prices have sensible short-term oscillations, they tend to revert back to a “normal” long-term equilibrium level. The past values of the changes in this uncertain factor help predict the future. We will use a model where the metal spot price is assumed to follow the stochastic process [19]:(28)dP=αln⁡P--ln⁡PPdt+σPdW.Let x=ln⁡P; applying Ito’s lemma allows the characterization of the log price by an Ornstein-Uhlenbeck stochastic mean reverting process: (29)dx=αx--xdt+σdWwith (30)x¯=ln⁡P--σ22α,where P- is the long-run equilibrium metal price, α measures the speed of mean reversion to the long-run mean log price P-, dW is an increment to a standard Brownian motion, and σ refers to the price volatility rate. The metal price adjustment mechanism is accounted for by market forces. The correct discrete-time format for the continuous-time process of mean reversion is the stationary first-order autoregressive process [20], so the sample path simulation equation for xt is performed by using exact discrete-time expression:(31)xt=xt-1e-αΔt+x-1-e-αΔt+N0,1σ1-e-2αΔt2α,where Δt is the fixed time interval from time t to t+1 and N(0,1) is the normally distributed random variable.By substituting (31) to P=ex, we have exact discrete-time equation for Pt, given by (32)Pt=eln⁡Pt-1e-αΔt+ln⁡P--σ2/2α1-e-αΔt+N0,1σ1-e-2αΔt/2α.In order to estimate the parameters of the mean reversion process, we run the following regression: (33)dxt+1=β0+β1xt+ε,where β0=αx-dt and β1=-αdt. Hence, if we regress observation dx against x, we can obtain estimates of β0 and β1. σ is the standard deviation obtained from the regression. The speed of mean reversion (α) is the negative of the slope, while the long-run equilibrium P- is the intercept estimate of that regression divided by the speed of mean reversion.LetP=Pt,t=0,1,…,T denote a price scenario with spot prices Pt, where Pt is determined by (32). Figure 1 presents a sample path of the metal price (e.g., lead) simulated using the above equation.Figure 1
One simulated lead price path on a yearly time resolution.In the purpose ofP simulation, we apply a metal price scenario for the time interval, n,n+T, with increment Δt=1.The unit value of metal concentrate sales price can be expressed as follows:(34)Vcont=fPt,mcon,mmr,where mmr is the metal recovery rate (%). The value of Vcon is calculated in different ways with respect to type of metal.


## 3.2.2. Volatility of Production Costs

Production costs (PC) are incurred directly in the production process. These costs include the ore and waste development of individual stopes, the actual stoping activities, the mine services providing logistical support to the miners, and the milling and processing of the ore at the plant. The uncertainties related to the future states of unit production costs are modeled with a special stochastic process, the geometric Brownian motion. Certain stochastic processes are functions of a Brownian motion process and these have many applications in finance, engineering, and the sciences. Some special processes are solutions of Itô-Doob type stochastic differential equations (Ladde, Sambandham) [21].In this model, we apply a continuous time process using the Itô-Doob type stochastic differential equation to describe movement of unit production costs. A general stochastic differential equation takes the following form:(35)dCOt=ρ·COt,tdt+σ·COt,t·dWt,COt0=CO0.Here, t≥t0, Wt is a Brownian motion, and COt>0; this is the cost process.C
O
t is called the geometric Brownian motion, which is solution of the following linear Itô-Doob type stochastic differential equation: (36)dCOt=ρ·COtdt+σ·COt·dWt,where ρ is the drift, σ is the volatility, and Wt is normalized Brownian motion.Using the Itô-Doob formula applied toCOt=ln⁡COt, we can solve this equation: (37)COt=COt-1·eρ-σ2/2Δt+N0,1σΔt.Equation (37) describes an operating cost scenario with spot costs COt.LetCO=COt,t=0,1,…,T denote a cost scenario with spot costs COt, where COt is determined by (37). Figure 2 presents a sample paths (s=1,2,…,S) of the unit production cost simulated using (37) S times.Figure 2
Simulated cost paths on a yearly time resolution.In the purpose of CO simulation, we apply a cost scenario for the time interval,n,n+T, with increment Δt=1.The annual production costs are expressed as follows:(38)PCyear=Qyear·COt.


## 3.3. Variables Based on Expert’s Knowledge

To estimate an adequate value of these variables, it is necessary to get opinions of experts dealing with them. Suppose we havep experts and each of them has given his/her opinion. The final value is expressed by averaging the opinions of experts. Table 1 presents the expert estimation process.Table 1
Expert estimation of the variables.
Expert
Value
E
1
E
2
⋯
E
p
x
1,1
x
1,2
⋯
x
1
,
p
E
(
x
1
)
=
x
1,1
+
x
1,2
+
⋯
+
x
1
,
p
p
x
2,1
x
2,2
⋯
x
2
,
p
E
(
x
2
)
=
x
2,1
+
x
2,2
+
⋯
+
x
2
,
p
p
x
5,1
x
5,2
⋯
x
5
,
p
E
(
x
5
)
=
x
5,1
+
x
5,2
+
⋯
+
x
5
,
p
p
x
6,1
x
6,2
⋯
x
6
,
p
E
(
x
6
)
=
x
6,1
+
x
6,2
+
⋯
+
x
6
,
p
p


## 3.4. DOL Forecasting Model Based on Expert’s Knowledge and Simulation

The model is developed on the basis of expert’s knowledge and simulation of DOL changing over defined time horizon and takes into account the variability of input parameters. By simulating a forecasting system, we imitate its action in order to measure its response (output) to different inputs. The advantage of simulating a system is the possibility of replicating its evolution as many times as necessary in independent conditions. The simulation allows analysts to describe the uncertainty of variables that influence the value of DOL by different time depending scenarios. The values of DOL are usually forecasted on an annual basis.The first objective of the use of simulation in the forecasting is to determine the distribution of the DOL from the variables that affect his performance, which results in the average or expected value of DOL for every year of defined time horizon. The relation between uncertain variables affecting the value of DOL is described by (22).For each simulation, the input values and DOL result represent one possible state of nature. Simulated values of DOL are obtained by performing the following calculations:(39)y^1sn+h+1=y^1x1,n+h,x2,n+hs,x3,n+hs,x4,n+hs,x5,n+h,x6,n+hs,h=0,1,2,…,T-1,y^0sn+h=y^1sn+h-y^1sn+h-1,h=1,2,…,T,s=1,2,…,S,where S denotes the number of simulations.Space of simulation fors=1 and k>n, where the first six rows concern the evolution path of the input variables while the last one concerns evolution path of the DOL, can be represented as follows:(40)DOLTs=1=s=1s=1⋯s=1x1,1x1,2⋯x1,Tx2,1x2,2⋯x2,Tx3,1x3,2⋯x3,Tx4,1x4,2⋯x4,Tx5,1x5,2⋯x5,Tx6,1x6,2⋯x6,T↓↓↓↓equation22equation22equation22equation22↓↓↓↓y^11s=1y^21s=1⋯y^T1s=1↓↓↓↓y^10s=1y^20s=1⋯y^T0s=1.Once the simulations DOLTs, s=1,2,…,S, have been obtained, they can be used to estimate the distribution of the DOL for every year of the defined time horizon. The expected values of the DOL can be expressed by the following vector space: (41)EDOLTS=y^10s=1y^20s=1⋯y^T0s=1y^10s=2y^20s=2⋯y^T0s=2⋮⋮⋱⋮y^10s=Sy^20s=S⋯y^T0s=S↓↓↓↓Ey^10Ey^20⋯Ey^T0.According to (41), the management of the company obtains the forecasting information vector of DOL: (42)FIVj,T=EDOL=EDOL1EDOL2⋯EDOLT.This vector concerns the information of the possible future states of the Degree of Operating Leverage over defined time horizon.Space of simulation fors=1 and k>n, where the first seven rows concern the evolution paths of the DOL and input variables while the rest concern evolution paths of the grey relational coefficient, can be represented as follows: (43)ξj,Ts=1=s=1s=1⋯s=1y^10y^20⋯y^T0x1,1x1,2⋯x1,Tx2,1x2,2⋯x2,Tx3,1x3,2⋯x3,Tx4,1x4,2⋯x4,Tx5,1x5,2⋯x5,Tx6,1x6,2⋯x6,T↓↓↓↓equation14equation14equation14equation14↓↓↓↓ξ1,1s=1ξ1,2s=1⋯ξ1,Ts=1⋮⋮⋱⋮ξ6,1s=1ξ6,2s=1⋯ξ6,Ts=1.Once the simulations ξj,Ts, s=1,2,…,S, have been obtained, they can be used to estimate the distribution of the ξj for every year of the defined time horizon. The expected values of the ξj can be expressed by the following vector space: (44)Eξj,TS=Eξ1,1SEξ1,2S⋯Eξ1,TSEξ2,1SEξ2,2S⋯Eξ2,TSEξ3,1SEξ3,2S⋯Eξ3,TSEξ4,1SEξ4,2S⋯Eξ4,TSEξ5,1SEξ5,2S⋯Eξ5,TSEξ6,1SEξ6,2S⋯Eξ6,TS.The expected values of the aggregated grey relational coefficient are calculated as follows: (45)ρj=1h∑h=1TEξj,hS,j=1,2,…,6.According to (45), the management of the company obtains the forecasting information vector of the aggregated grey relational coefficient:(46)ρj=ρ1ρ2⋮ρ6.


## 4. Application of the Procedure

### 4.1. A Numerical Example Statement

A management of an operating underground zinc mine is trying to forecast the future states of operational efficiency for five years ahead. For this problem, the input parameters that are required for the forecasting are given in Tables2 and 3. In order to define values of production rate, fixed costs, working days, and degree of use of production capacity, three experts are included into process of estimation. Note that the situation is hypothetical and numbers used are to permit calculation. Values of the input parameters are similar to the real environment of mining.Table 2
Input parameters of observed period.
Parameter
Observed period
Year
1
2
3
4
5
DOL
1.219
1.421
1.333
1.676
1.437
Production rate (t/year)
100000
95000
100000
100000
100000
Fixed costs (USD)
1200000
1500000
1100000
1400000
1600000
Revenues (USD)
Ore grade (%)
4.11
3.96
4.04
3.82
3.87
Mill recovery rate (%)
78.52
78.01
81.05
79.02
81.56
Metal content of the concentrate (%)
53.8
53.8
53.8
53.8
53.8
Metal recovery rate (%) mmr=mcon%-8·100mcon%c∫c≤85%;mcon%-8mcon%-8·100mcon%>85%;85%  mmr=53.8-8·10053.8∫cc=85.13;85%
85
85
85
85
85
Zinc metal price (USD/t)
2160
2195
1950
1910
2160
Unit value of metal concentrate sales price (USD/t) Vcon=Pmcon-mmr=P·0.85
1836
1866
1658
1624
1836
Revenues (USD)
11065884
10773877
10137834
9145279
10823137
Production costs (USD)
Unit production costs (USD/t)
43.87
60.11
57.40
56.76
55.64
Production costs (USD)
4387000
5714000
5740000
5676000
5564000
Working days (day/year)
330
310
340
300
320
Degree of use of production capacity (%)
90
85
93
82
88Table 3
Input parameters required for simulation of DOL.
Parameter
Value
Year
6
7
8
9
10
Production rate (t/year)
100000
105000
97000
100000
110000
Fixed costs (USD)
1600000
1800000
1500000
1700000
2000000
Revenues (USD)
Ore grade (%)-normal distribution
Min. 3.45; medium 4.06; max. 4.68; volatility 0.205
Mill recovery rate (%)-uniform distribution
Min. 77; medium 78.5; max. 80; volatility 0.866
Metal content of the concentrate (%)
53.8
Metal recovery rate (%)
85
Zinc metal price (USD/t)-mean reversion process
Spot value 2113; equilibrium metal price 2277; speed ofmean reversion 0.9221; price volatility rate 0.2734
Production costs (USD)
Unit production costs (USD/t)-geometric Brownian motion
Spot value 65; drift 0.02382; cost volatility 0.09351
Working days (day/year)
330
330
340
340
340
Degree of use of production capacity (%)
86
90
92
85
93
Number of simulations
500


### 4.2. A Numerical Example Solution

Step 1 (forecasting model).
Transformations of the original data sequences by AGO are as follows:(47)Y0=1.2191.4211.3331.6761.437,Y1=1.2192.6403.9735.6497.086,X10=10000095000100000100000100000,X11=100000195000295000395000495000,X20=12000001500000110000014000001600000,X21=12000002700000380000052000006800000,X30=110658841077387710137834914527910823137,X31=1106588421839761319775954112287451946011,X40=43870005714000574000056760005564000,X41=438700010101000158410002151700027081000,X50=330310340300320,X51=33064098012801600,X60=9085938288,X61=90175268350438.Calculation of the average generation of adjacent sample sequence of Y(1)(k) is as follows: (48)Y1=1.2192.6403.9735.6497.086,Z1Y1=1.9293.3074.8126.369.Construction of the matrix B and matrix Y is as follows: (49)B=-1.21919500027000002183976110101000640175-3.30729500038000003197759515841000980268-4.812395000520000041122874215170001280350-6.369495000680000051946011270810001600438,Y=1.4211.3331.6761.437.The constructed AGO sequence of our problem is represented by the following matrix form: (50)1.4211.3331.6761.437=-1.21919500027000002183976110101000640175-3.30729500038000003197759515841000980268-4.812395000520000041122874215170001280350-6.369495000680000051946011270810001600438ab2⋮bN.According to a^=BTB-1BTY we obtain the following factors of the matrix a^=1.5379;0.0001129;-1.03377·10-6;1.304·10-7;-1.4598·10-7;-0.0042024;-0.0796303T. Values of the determined coefficients are a=1.5379; b1=0.0001129; b2=-1.03377·10-6; b3=1.304·10-7; b4=-1.4598·10-7; b5=-0.0042024; b6=-0.0796303.The AGO time sequence response formula of DOL is(51)y^1k+1=1.219-7.3462·10-5x11k+1+6.748·10-7x21k+1-8.479·10-8x31k+1+9.492·10-8x41k+1+0.0027325x51k+1+0.0517784x61k+1e-1.5379k+7.3462·10-5x11k+1-6.748·10-7x21k+1+8.479·10-8x31k+1-9.492·10-8x41k+1-0.0027325x51k+1-0.0517784x61k+1;k=0,1,2,3,4.AGO and inverse AGO values of the time sequence response formula of DOL are represented by Table 4.Table 4
Fitted values of DOL.
Error
Year
1
2
3
4
5
DOL-AGO values
2.292
3.642
5.291
6.546
DOL-inverse AGO values
1.219
1.073
1.350
1.649
1.255The adequacy of the obtained time sequence response formula is obtained by using (12) and summarized in Table 5.Table 5
Error estimation of the model.
Error
Year
1
2
3
4
5
DOL-observed values
1.219
1.421
1.333
1.676
1.437
DOL-fitted values
1.219
1.073
1.350
1.649
1.255
RPE (%)
0
24.47
1.28
1.63
12.70
ARPE (%)
10.02Step 2 (grey relevance analysis).
Applying (13) we obtain the following results of equalization (see Tables 6 and 7).Table 6
Sum of the major and influencing sequence variables.
DOL(k)
1.219
1.421
1.334
1.677
1.437
X
1(k)
100000
95000
100000
100000
100000
X
2(k)
1200000
1500000
1100000
1400000
1600000
X
3(k)
11065884
10773876
10137834
9145278
10823136
X
4(k)
4387000
5714000
5740000
5676000
5564000
X
5(k)
330
310
340
300
320
X
6(k)
90
85
93
82
88
Sum
16753305
18083273
17078268
16321662
18087545Table 7
Equalized values of the major and influencing sequence variables.
DOL′(k)
7.28E − 08
7.86E − 08
7.81E − 08
1.03E − 07
7.95E − 08
X
1
′
k
0.005969
0.005253
0.005855
0.006127
0.005529
X
2
′
k
0.071628
0.08295
0.064409
0.085776
0.088459
X
3
′
k
0.660519
0.595792
0.59361
0.560315
0.598375
X
4
′
k
0.261859
0.315983
0.3361
0.347759
0.307615
X
5
′
k
1.97E − 05
1.71E − 05
1.99E − 05
1.84E − 05
1.77E − 05
X
6
′
k
5.37E − 06
4.7E − 06
5.45E − 06
5.02E − 06
4.87E − 06Absolute values ofΔjk=DOL′k-xj′k are represented by Table 8.Table 8
Absolute values ofΔj(k).
Δ
1
k
0.005969
0.005253
0.005855
0.006127
0.005529
Δ
2
k
0.071628
0.08295
0.064409
0.085775
0.088459
Δ
3
k
0.660519
0.595792
0.59361
0.560315
0.598375
Δ
4
k
0.261859
0.315983
0.3361
0.347759
0.307615
Δ
5
k
1.96E − 05
1.71E − 05
1.98E − 05
1.83E − 05
1.76E − 05
Δ
6
k
5.30E − 06
4.62E − 06
5.37E − 06
4.92E − 06
4.79E − 06Extreme values ofΔjk are represented by Table 9.Table 9
Min. and max. ofΔj(k).
Min.
Max.
Δ
1
(
k
)
0.005253
0.006127
Δ
2
k
0.064409
0.088459
Δ
3
k
0.560315
0.660519
Δ
4
k
0.261859
0.347759
Δ
5
k
1.71E − 05
1.98E − 05
Δ
6
k
4.62E − 06
5.37E − 06Grey relational (relevance) analysis ofGM(1,6) is represented by Table 10.Table 10
Correlation matrix.
DOL(k)/Xj θ=0.5
DOL(1)
DOL(2)
DOL(3)
DOL(4)
DOL(5)
Aggregated value
X
1
0.92078
1
0.93251
0.90497
0.96797
0.94524
X
2
0.93769
0.85421
1
0.83565
0.81875
0.88926
X
3
0.89886
0.96169
0.96396
1
0.95901
0.95670
X
4
1
0.88951
0.85442
0.83532
0.90497
0.89684
X
5
0.91331
1
0.90701
0.95695
0.98009
0.95147
X
6
0.91514
1
0.90739
0.96063
0.97805
0.95224According to Table10, we obtain the rank order of relevancy of the influencing variables over observed time as follows:Year 1:X4>X2>X1>X6>X5>X3.Year 2:(X1=X5=X6)>X3>X4>X2.Year 3:X2>X3>X1>X6>X5>X4.Year 4:X3>X6>X5>X1>X2>X4.Year 5:X5>X6>X1>X3>X4>X2.According to aggregated values, we obtain the final rank order of relevancy of the influencing variables as follows:X3>X6>X5>X1>X4>X2.Step 3 (simulation of the influencing sequence variables and DOL).
Scenario of zinc price is represented by Table11.Table 11
Mean reversion process of zinc price.
Year 6
P
6
=
e
ln
⁡
2113
·
e
-
0.9221
+
ln
⁡
2277
-
0.2734
2
/
(
2
·
0.9221
)
1
-
e
-
0.9221
+
N
0,1
·
0.2734
·
1
-
e
-
2
·
0.9221
/
2
·
0.9221
Year 7
P
7
=
e
ln
⁡
P
6
·
e
-
0.9221
+
ln
⁡
2277
-
0.2734
2
/
(
2
·
0.9221
)
1
-
e
-
0.9221
+
N
0,1
·
0.2734
·
1
-
e
-
2
·
0.9221
/
2
·
0.9221
Year 8
P
8
=
e
ln
⁡
P
7
·
e
-
0.9221
+
ln
⁡
2277
-
0.2734
2
/
(
2
·
0.9221
)
1
-
e
-
0.9221
+
N
0,1
·
0.2734
·
1
-
e
-
2
·
0.9221
/
2
·
0.9221
Year 9
P
9
=
e
ln
⁡
P
8
·
e
-
0.9221
+
ln
⁡
2277
-
0.2734
2
/
(
2
·
0.9221
)
1
-
e
-
0.9221
+
N
0,1
·
0.2734
·
1
-
e
-
2
·
0.9221
/
2
·
0.9221
Year 10
P
10
=
e
ln
⁡
P
9
·
e
-
0.9221
+
ln
⁡
2277
-
0.2734
2
/
(
2
·
0.9221
)
1
-
e
-
0.9221
+
N
0,1
·
0.2734
·
1
-
e
-
2
·
0.9221
/
2
·
0.9221Scenario of unit production costs is represented by Table12.Table 12
Geometric Brownian motion of unit production costs.
Year 6
CO
6
=
65
·
e
0.02382
-
0.09351
2
/
2
+
N
0,1
·
0.09351
Year 7
CO
7
=
CO
6
·
e
0.02382
-
0.09351
2
/
2
+
N
0,1
·
0.09351
Year 8
CO
8
=
CO
7
·
e
0.02382
-
0.09351
2
/
2
+
N
0,1
·
0.09351
Year 9
CO
9
=
CO
8
·
e
0.02382
-
0.09351
2
/
2
+
N
0,1
·
0.09351
Year 10
CO
10
=
CO
9
·
e
0.02382
-
0.09351
2
/
2
+
N
0,1
·
0.09351Results of forecasting of DOL generated by the GM(1,6) (see (51)) for a single simulation are represented by Table 13.Table 13
One simulation path of DOL.
Parameter
Value
Year
6
7
8
9
10
Production rate (t/year)
100000
105000
97000
100000
110000
Fixed costs (USD)
1600000
1800000
1500000
1700000
2000000
Revenues (USD)
Ore grade (%)-normal distribution
4.07
3.78
4.31
4.01
3.86
Mill recovery rate (%)-uniform distribution
79.27
78.54
78.62
77.71
78.32
Metal content of the concentrate (%)
53.8
53.8
53.8
53.8
53.8
Metal recovery rate (%)
85
85
85
85
85
Zinc metal price (USD/t)-mean reversion process
2025
1936
2064
2409
2120
Revenues (USD)
13040642
12152559
13638024
15254913
14228955
Production costs (USD)
Unit production costs (USD/t)-geometric Brownian motion
71.56
88.33
90.77
98.01
95.87
Production costs (USD)
7156388
9274793
8805646
9801869
10546252
Working days (day/year)
330
330
340
340
340
Degree of use of production capacity (%)
86
90
92
85
93
y^1(k+1)=1.219-7.3462·10-5x11(k+1)+6.748·10-7x21(k+1)-8.479·10-8x31(k+1)+9.492·10-8x41(k+1) hhhhhhhhhh+0.0027325x51(k+1)+0.0517784x61(k+1)e-1.5379k+7.3462·10-5x11(k+1)-6.748·10-7x21(k+1) hhh+8.479·10-8x31(k+1)-9.492·10-8x41(k+1)-0.0027325x51(k+1)-0.0517784x61(k+1);                        k=5+h; h=0,1,2,3,4
DOL-AGO
7.892
8.981
9.723
10.955
12.147
DOL-inverse AGO
1.346
1.089
0.742
1.232
1.192This act is repeated five hundred times and set of five hundred possible states of nature is obtained for every year of the time horizon. The expected DOL values,E(DOL(6)),E(DOL(7)),…,E(DOL(10)), are calculated from the simulation results, that is, from the sets of possible states of nature. The forecasting of the DOL is represented by Figure 3 and summary Table 14.Table 14
Summary statistics.
Statistic parameter
Value
Year
6
7
8
9
10
Sample
500
500
500
500
500
Mean (expected) value of DOL
1.491
1.525
0.951
1.403
1.574
Median
1.469
1.494
0.945
1.397
1.571
Standard deviation
0.2512
0.2745
0.2631
0.2740
0.3163
Max.
2.470
2.543
1.766
2.471
2.940
Min.
0.9073
0.8394
0.2065
0.5978
0.5828
Range
1.562
1.703
1.559
1.872
2.357
Q (75%)
1.646
1.698
1.105
1.573
1.752
Q (25%)
1.293
1.340
0.761
1.227
1.379
Q (75%)–Q (25%)
0.352
0.358
0.344
0.346
0.373
Asymmetry (Skewness)
0.6235
0.5085
0.4129
0.3726
0.2453
Kurtosis
0.3917
0.3205
0.2378
0.7502
1.1323
Standard error
0.01123
0.01227
0.01176
0.01225
0.01414
Reliability coefficient (1-α)
0.95
0.95
0.95
0.95
0.95
Reliability factor
1.96
1.96
1.96
1.96
1.96
Reliability interval-upper limit
1.513
1.549
0.973
1.427
1.601
Reliability interval-lower limit
1.469
1.501
0.927
1.378
1.545
Mean (expected) value (%)
52.90
53.80
51.40
50.20
50.70Figure 3
Original, fitting, one path forecasting, and expected values of DOL after 500 simulations.Probability and cumulative density function of forecasted DOL for 6th year are represented by Figure4.Figure 4
Distribution of forecasted DOL for the 6th year.Step 4 (grey relevance analysis of results obtained by simulation).
Applying (13) we obtain the following results of equalization (see Tables 15 and 16).Table 15
Sum of the major and influencing sequence variables fors=1.
DOL(k)
1.346
1.089
0.742
1.232
1.192
X
1(k)
100000
105000
97000
100000
110000
X
2(k)
1600000
1800000
1500000
1700000
2000000
X
3(k)
13040642
12152559
13638024
15254913
14228955
X
4(k)
7156388
9274793
8805646
9801869
10546252
X
5(k)
330
330
340
340
340
X
6(k)
86
90
92
85
93
Sum
21897447
23332773
24041102
26857208
26885641Table 16
Equalized values of the major and influencing sequence variables fors=1.
DOL′(k)
6.15E − 08
4.67E − 08
3.09E − 08
4.59E − 08
4.43E − 08
X
1
′
k
0.004567
0.0045
0.004035
0.003723
0.004091
X
2
′
k
0.073068
0.077145
0.062393
0.063298
0.074389
X
3
′
k
0.595533
0.520836
0.567279
0.568001
0.52924
X
4
′
k
0.326814
0.397501
0.366275
0.364962
0.392263
X
5
′
k
1.51E − 05
1.41E − 05
1.41E − 05
1.27E − 05
1.26E − 05
X
6
′
k
3.93E − 06
3.86E − 06
3.83E − 06
3.16E − 06
3.46E − 06Grey relational (relevance) analysis of GM(1,6) for one path of simulation is represented by Table17.Table 17
Correlation matrix fors=1.
DOL(k)/Xj θ=0.5
DOL(1)
DOL(2)
DOL(3)
DOL(4)
DOL(5)
Aggregated value
X
1
0.876886
0.885498
0.950716
1
0.94227
0.931074
X
2
0.904383
0.87252
1
0.99112
0.893805
0.932366
X
3
0.916382
1
0.946311
0.945523
0.989839
0.959611
X
4
1
0.881448
0.930161
0.932326
0.889259
0.926639
X
5
0.893087
0.930804
0.930155
0.999411
1
0.950691
X
6
0.871197
0.879595
0.881844
1
0.944697
0.915467According to Table17, we obtain the simulated rank order of relevancy of the influencing variables over time as follows:Year 1:X4>X3>X2>X5>X1>X6.Year 2:X3>X5>X1>X4>X6>X2.Year 3:X2>X1>X3>X4>X5>X6.Year 4:(X1=X6)>X5>X2>X3>X4.Year 5:X5>X3>X6>X1>X2>X4.According to aggregated values, we obtain the final simulated rank order of relevancy of the influencing variables as follows:X3>X5>X2>X1>X4>X6.Since 500 simulations are performed, we obtain the correlation matrix as in Table18.Table 18
Expected correlation matrix obtained after 500 simulations.
DOL(k)/Xj θ=0.5
DOL(1)
DOL(2)
DOL(3)
DOL(4)
DOL(5)
Aggregated value
X
1
0.912041
0.921771
0.925088
0.925480
0.939693
0.924814
X
2
0.928070
0.896937
0.957765
0.913864
0.876159
0.914559
X
3
0.934076
0.940820
0.935899
0.950123
0.953475
0.942879
X
4
0.922037
0.922130
0.911800
0.901827
0.905551
0.912669
X
5
0.890792
0.938745
0.870543
0.895710
0.958414
0.910841
X
6
0.909324
0.925922
0.851819
0.933881
0.937900
0.911769According to Table18, we obtain the expected rank order of relevancy of the influencing variables over time and final expected rank order as follows:Year 1:X3>X2>X4>X1>X6>X5.Year 2:X3>X5>X6>X4>X1>X2.Year 3:X2>X3>X1>X4>X5>X6.Year 4:X3>X6>X1>X2>X4>X5.Year 5:X5>X3>X1>X6>X4>X2.Final:X3>X1>X2>X4>X6>X5.If we compare the final rank order of relevancy of the influencing variables obtained by observation to the final expected rank order of influencing variables obtained by simulation, it can be seen that the revenues retain the greatest relevance. Level of relevancy of the influencing variable related to working days is significantly changed from the second level in the past to the sixth in the future. Influencing variable related to the production rate is increased its relevance from the fourth level in the past to the second in the future, while the rest variables retain approximately the same level of relevancy.Key information obtained from the values of DOL can be summarized in Table19.Table 19
DOL’s indications.
Parameter
DOL value
High
Low
Variable costs (revenues and fixed costs are fixed)
High
Low
Fixed costs (revenues and variable costs are fixed)
Low
High
Revenues up (variable and fixed costs are fixed)
Greater profits
Lower profits
Revenues down (variable and fixed costs are fixed)
Greater losses
Lower lossesFrom Table19, we can see why it matters if value of DOL is high or low. When DOL is high, a change in revenues results in large change in profit or loss. On the other hand, when DOL is low, a change in revenues results in small change in profit or loss.


## 4.1. A Numerical Example Statement

A management of an operating underground zinc mine is trying to forecast the future states of operational efficiency for five years ahead. For this problem, the input parameters that are required for the forecasting are given in Tables2 and 3. In order to define values of production rate, fixed costs, working days, and degree of use of production capacity, three experts are included into process of estimation. Note that the situation is hypothetical and numbers used are to permit calculation. Values of the input parameters are similar to the real environment of mining.Table 2
Input parameters of observed period.
Parameter
Observed period
Year
1
2
3
4
5
DOL
1.219
1.421
1.333
1.676
1.437
Production rate (t/year)
100000
95000
100000
100000
100000
Fixed costs (USD)
1200000
1500000
1100000
1400000
1600000
Revenues (USD)
Ore grade (%)
4.11
3.96
4.04
3.82
3.87
Mill recovery rate (%)
78.52
78.01
81.05
79.02
81.56
Metal content of the concentrate (%)
53.8
53.8
53.8
53.8
53.8
Metal recovery rate (%) mmr=mcon%-8·100mcon%c∫c≤85%;mcon%-8mcon%-8·100mcon%>85%;85%  mmr=53.8-8·10053.8∫cc=85.13;85%
85
85
85
85
85
Zinc metal price (USD/t)
2160
2195
1950
1910
2160
Unit value of metal concentrate sales price (USD/t) Vcon=Pmcon-mmr=P·0.85
1836
1866
1658
1624
1836
Revenues (USD)
11065884
10773877
10137834
9145279
10823137
Production costs (USD)
Unit production costs (USD/t)
43.87
60.11
57.40
56.76
55.64
Production costs (USD)
4387000
5714000
5740000
5676000
5564000
Working days (day/year)
330
310
340
300
320
Degree of use of production capacity (%)
90
85
93
82
88Table 3
Input parameters required for simulation of DOL.
Parameter
Value
Year
6
7
8
9
10
Production rate (t/year)
100000
105000
97000
100000
110000
Fixed costs (USD)
1600000
1800000
1500000
1700000
2000000
Revenues (USD)
Ore grade (%)-normal distribution
Min. 3.45; medium 4.06; max. 4.68; volatility 0.205
Mill recovery rate (%)-uniform distribution
Min. 77; medium 78.5; max. 80; volatility 0.866
Metal content of the concentrate (%)
53.8
Metal recovery rate (%)
85
Zinc metal price (USD/t)-mean reversion process
Spot value 2113; equilibrium metal price 2277; speed ofmean reversion 0.9221; price volatility rate 0.2734
Production costs (USD)
Unit production costs (USD/t)-geometric Brownian motion
Spot value 65; drift 0.02382; cost volatility 0.09351
Working days (day/year)
330
330
340
340
340
Degree of use of production capacity (%)
86
90
92
85
93
Number of simulations
500


## 4.2. A Numerical Example Solution

Step 1 (forecasting model).
Transformations of the original data sequences by AGO are as follows:(47)Y0=1.2191.4211.3331.6761.437,Y1=1.2192.6403.9735.6497.086,X10=10000095000100000100000100000,X11=100000195000295000395000495000,X20=12000001500000110000014000001600000,X21=12000002700000380000052000006800000,X30=110658841077387710137834914527910823137,X31=1106588421839761319775954112287451946011,X40=43870005714000574000056760005564000,X41=438700010101000158410002151700027081000,X50=330310340300320,X51=33064098012801600,X60=9085938288,X61=90175268350438.Calculation of the average generation of adjacent sample sequence of Y(1)(k) is as follows: (48)Y1=1.2192.6403.9735.6497.086,Z1Y1=1.9293.3074.8126.369.Construction of the matrix B and matrix Y is as follows: (49)B=-1.21919500027000002183976110101000640175-3.30729500038000003197759515841000980268-4.812395000520000041122874215170001280350-6.369495000680000051946011270810001600438,Y=1.4211.3331.6761.437.The constructed AGO sequence of our problem is represented by the following matrix form: (50)1.4211.3331.6761.437=-1.21919500027000002183976110101000640175-3.30729500038000003197759515841000980268-4.812395000520000041122874215170001280350-6.369495000680000051946011270810001600438ab2⋮bN.According to a^=BTB-1BTY we obtain the following factors of the matrix a^=1.5379;0.0001129;-1.03377·10-6;1.304·10-7;-1.4598·10-7;-0.0042024;-0.0796303T. Values of the determined coefficients are a=1.5379; b1=0.0001129; b2=-1.03377·10-6; b3=1.304·10-7; b4=-1.4598·10-7; b5=-0.0042024; b6=-0.0796303.The AGO time sequence response formula of DOL is(51)y^1k+1=1.219-7.3462·10-5x11k+1+6.748·10-7x21k+1-8.479·10-8x31k+1+9.492·10-8x41k+1+0.0027325x51k+1+0.0517784x61k+1e-1.5379k+7.3462·10-5x11k+1-6.748·10-7x21k+1+8.479·10-8x31k+1-9.492·10-8x41k+1-0.0027325x51k+1-0.0517784x61k+1;k=0,1,2,3,4.AGO and inverse AGO values of the time sequence response formula of DOL are represented by Table 4.Table 4
Fitted values of DOL.
Error
Year
1
2
3
4
5
DOL-AGO values
2.292
3.642
5.291
6.546
DOL-inverse AGO values
1.219
1.073
1.350
1.649
1.255The adequacy of the obtained time sequence response formula is obtained by using (12) and summarized in Table 5.Table 5
Error estimation of the model.
Error
Year
1
2
3
4
5
DOL-observed values
1.219
1.421
1.333
1.676
1.437
DOL-fitted values
1.219
1.073
1.350
1.649
1.255
RPE (%)
0
24.47
1.28
1.63
12.70
ARPE (%)
10.02Step 2 (grey relevance analysis).
Applying (13) we obtain the following results of equalization (see Tables 6 and 7).Table 6
Sum of the major and influencing sequence variables.
DOL(k)
1.219
1.421
1.334
1.677
1.437
X
1(k)
100000
95000
100000
100000
100000
X
2(k)
1200000
1500000
1100000
1400000
1600000
X
3(k)
11065884
10773876
10137834
9145278
10823136
X
4(k)
4387000
5714000
5740000
5676000
5564000
X
5(k)
330
310
340
300
320
X
6(k)
90
85
93
82
88
Sum
16753305
18083273
17078268
16321662
18087545Table 7
Equalized values of the major and influencing sequence variables.
DOL′(k)
7.28E − 08
7.86E − 08
7.81E − 08
1.03E − 07
7.95E − 08
X
1
′
k
0.005969
0.005253
0.005855
0.006127
0.005529
X
2
′
k
0.071628
0.08295
0.064409
0.085776
0.088459
X
3
′
k
0.660519
0.595792
0.59361
0.560315
0.598375
X
4
′
k
0.261859
0.315983
0.3361
0.347759
0.307615
X
5
′
k
1.97E − 05
1.71E − 05
1.99E − 05
1.84E − 05
1.77E − 05
X
6
′
k
5.37E − 06
4.7E − 06
5.45E − 06
5.02E − 06
4.87E − 06Absolute values ofΔjk=DOL′k-xj′k are represented by Table 8.Table 8
Absolute values ofΔj(k).
Δ
1
k
0.005969
0.005253
0.005855
0.006127
0.005529
Δ
2
k
0.071628
0.08295
0.064409
0.085775
0.088459
Δ
3
k
0.660519
0.595792
0.59361
0.560315
0.598375
Δ
4
k
0.261859
0.315983
0.3361
0.347759
0.307615
Δ
5
k
1.96E − 05
1.71E − 05
1.98E − 05
1.83E − 05
1.76E − 05
Δ
6
k
5.30E − 06
4.62E − 06
5.37E − 06
4.92E − 06
4.79E − 06Extreme values ofΔjk are represented by Table 9.Table 9
Min. and max. ofΔj(k).
Min.
Max.
Δ
1
(
k
)
0.005253
0.006127
Δ
2
k
0.064409
0.088459
Δ
3
k
0.560315
0.660519
Δ
4
k
0.261859
0.347759
Δ
5
k
1.71E − 05
1.98E − 05
Δ
6
k
4.62E − 06
5.37E − 06Grey relational (relevance) analysis ofGM(1,6) is represented by Table 10.Table 10
Correlation matrix.
DOL(k)/Xj θ=0.5
DOL(1)
DOL(2)
DOL(3)
DOL(4)
DOL(5)
Aggregated value
X
1
0.92078
1
0.93251
0.90497
0.96797
0.94524
X
2
0.93769
0.85421
1
0.83565
0.81875
0.88926
X
3
0.89886
0.96169
0.96396
1
0.95901
0.95670
X
4
1
0.88951
0.85442
0.83532
0.90497
0.89684
X
5
0.91331
1
0.90701
0.95695
0.98009
0.95147
X
6
0.91514
1
0.90739
0.96063
0.97805
0.95224According to Table10, we obtain the rank order of relevancy of the influencing variables over observed time as follows:Year 1:X4>X2>X1>X6>X5>X3.Year 2:(X1=X5=X6)>X3>X4>X2.Year 3:X2>X3>X1>X6>X5>X4.Year 4:X3>X6>X5>X1>X2>X4.Year 5:X5>X6>X1>X3>X4>X2.According to aggregated values, we obtain the final rank order of relevancy of the influencing variables as follows:X3>X6>X5>X1>X4>X2.Step 3 (simulation of the influencing sequence variables and DOL).
Scenario of zinc price is represented by Table11.Table 11
Mean reversion process of zinc price.
Year 6
P
6
=
e
ln
⁡
2113
·
e
-
0.9221
+
ln
⁡
2277
-
0.2734
2
/
(
2
·
0.9221
)
1
-
e
-
0.9221
+
N
0,1
·
0.2734
·
1
-
e
-
2
·
0.9221
/
2
·
0.9221
Year 7
P
7
=
e
ln
⁡
P
6
·
e
-
0.9221
+
ln
⁡
2277
-
0.2734
2
/
(
2
·
0.9221
)
1
-
e
-
0.9221
+
N
0,1
·
0.2734
·
1
-
e
-
2
·
0.9221
/
2
·
0.9221
Year 8
P
8
=
e
ln
⁡
P
7
·
e
-
0.9221
+
ln
⁡
2277
-
0.2734
2
/
(
2
·
0.9221
)
1
-
e
-
0.9221
+
N
0,1
·
0.2734
·
1
-
e
-
2
·
0.9221
/
2
·
0.9221
Year 9
P
9
=
e
ln
⁡
P
8
·
e
-
0.9221
+
ln
⁡
2277
-
0.2734
2
/
(
2
·
0.9221
)
1
-
e
-
0.9221
+
N
0,1
·
0.2734
·
1
-
e
-
2
·
0.9221
/
2
·
0.9221
Year 10
P
10
=
e
ln
⁡
P
9
·
e
-
0.9221
+
ln
⁡
2277
-
0.2734
2
/
(
2
·
0.9221
)
1
-
e
-
0.9221
+
N
0,1
·
0.2734
·
1
-
e
-
2
·
0.9221
/
2
·
0.9221Scenario of unit production costs is represented by Table12.Table 12
Geometric Brownian motion of unit production costs.
Year 6
CO
6
=
65
·
e
0.02382
-
0.09351
2
/
2
+
N
0,1
·
0.09351
Year 7
CO
7
=
CO
6
·
e
0.02382
-
0.09351
2
/
2
+
N
0,1
·
0.09351
Year 8
CO
8
=
CO
7
·
e
0.02382
-
0.09351
2
/
2
+
N
0,1
·
0.09351
Year 9
CO
9
=
CO
8
·
e
0.02382
-
0.09351
2
/
2
+
N
0,1
·
0.09351
Year 10
CO
10
=
CO
9
·
e
0.02382
-
0.09351
2
/
2
+
N
0,1
·
0.09351Results of forecasting of DOL generated by the GM(1,6) (see (51)) for a single simulation are represented by Table 13.Table 13
One simulation path of DOL.
Parameter
Value
Year
6
7
8
9
10
Production rate (t/year)
100000
105000
97000
100000
110000
Fixed costs (USD)
1600000
1800000
1500000
1700000
2000000
Revenues (USD)
Ore grade (%)-normal distribution
4.07
3.78
4.31
4.01
3.86
Mill recovery rate (%)-uniform distribution
79.27
78.54
78.62
77.71
78.32
Metal content of the concentrate (%)
53.8
53.8
53.8
53.8
53.8
Metal recovery rate (%)
85
85
85
85
85
Zinc metal price (USD/t)-mean reversion process
2025
1936
2064
2409
2120
Revenues (USD)
13040642
12152559
13638024
15254913
14228955
Production costs (USD)
Unit production costs (USD/t)-geometric Brownian motion
71.56
88.33
90.77
98.01
95.87
Production costs (USD)
7156388
9274793
8805646
9801869
10546252
Working days (day/year)
330
330
340
340
340
Degree of use of production capacity (%)
86
90
92
85
93
y^1(k+1)=1.219-7.3462·10-5x11(k+1)+6.748·10-7x21(k+1)-8.479·10-8x31(k+1)+9.492·10-8x41(k+1) hhhhhhhhhh+0.0027325x51(k+1)+0.0517784x61(k+1)e-1.5379k+7.3462·10-5x11(k+1)-6.748·10-7x21(k+1) hhh+8.479·10-8x31(k+1)-9.492·10-8x41(k+1)-0.0027325x51(k+1)-0.0517784x61(k+1);                        k=5+h; h=0,1,2,3,4
DOL-AGO
7.892
8.981
9.723
10.955
12.147
DOL-inverse AGO
1.346
1.089
0.742
1.232
1.192This act is repeated five hundred times and set of five hundred possible states of nature is obtained for every year of the time horizon. The expected DOL values,E(DOL(6)),E(DOL(7)),…,E(DOL(10)), are calculated from the simulation results, that is, from the sets of possible states of nature. The forecasting of the DOL is represented by Figure 3 and summary Table 14.Table 14
Summary statistics.
Statistic parameter
Value
Year
6
7
8
9
10
Sample
500
500
500
500
500
Mean (expected) value of DOL
1.491
1.525
0.951
1.403
1.574
Median
1.469
1.494
0.945
1.397
1.571
Standard deviation
0.2512
0.2745
0.2631
0.2740
0.3163
Max.
2.470
2.543
1.766
2.471
2.940
Min.
0.9073
0.8394
0.2065
0.5978
0.5828
Range
1.562
1.703
1.559
1.872
2.357
Q (75%)
1.646
1.698
1.105
1.573
1.752
Q (25%)
1.293
1.340
0.761
1.227
1.379
Q (75%)–Q (25%)
0.352
0.358
0.344
0.346
0.373
Asymmetry (Skewness)
0.6235
0.5085
0.4129
0.3726
0.2453
Kurtosis
0.3917
0.3205
0.2378
0.7502
1.1323
Standard error
0.01123
0.01227
0.01176
0.01225
0.01414
Reliability coefficient (1-α)
0.95
0.95
0.95
0.95
0.95
Reliability factor
1.96
1.96
1.96
1.96
1.96
Reliability interval-upper limit
1.513
1.549
0.973
1.427
1.601
Reliability interval-lower limit
1.469
1.501
0.927
1.378
1.545
Mean (expected) value (%)
52.90
53.80
51.40
50.20
50.70Figure 3
Original, fitting, one path forecasting, and expected values of DOL after 500 simulations.Probability and cumulative density function of forecasted DOL for 6th year are represented by Figure4.Figure 4
Distribution of forecasted DOL for the 6th year.Step 4 (grey relevance analysis of results obtained by simulation).
Applying (13) we obtain the following results of equalization (see Tables 15 and 16).Table 15
Sum of the major and influencing sequence variables fors=1.
DOL(k)
1.346
1.089
0.742
1.232
1.192
X
1(k)
100000
105000
97000
100000
110000
X
2(k)
1600000
1800000
1500000
1700000
2000000
X
3(k)
13040642
12152559
13638024
15254913
14228955
X
4(k)
7156388
9274793
8805646
9801869
10546252
X
5(k)
330
330
340
340
340
X
6(k)
86
90
92
85
93
Sum
21897447
23332773
24041102
26857208
26885641Table 16
Equalized values of the major and influencing sequence variables fors=1.
DOL′(k)
6.15E − 08
4.67E − 08
3.09E − 08
4.59E − 08
4.43E − 08
X
1
′
k
0.004567
0.0045
0.004035
0.003723
0.004091
X
2
′
k
0.073068
0.077145
0.062393
0.063298
0.074389
X
3
′
k
0.595533
0.520836
0.567279
0.568001
0.52924
X
4
′
k
0.326814
0.397501
0.366275
0.364962
0.392263
X
5
′
k
1.51E − 05
1.41E − 05
1.41E − 05
1.27E − 05
1.26E − 05
X
6
′
k
3.93E − 06
3.86E − 06
3.83E − 06
3.16E − 06
3.46E − 06Grey relational (relevance) analysis of GM(1,6) for one path of simulation is represented by Table17.Table 17
Correlation matrix fors=1.
DOL(k)/Xj θ=0.5
DOL(1)
DOL(2)
DOL(3)
DOL(4)
DOL(5)
Aggregated value
X
1
0.876886
0.885498
0.950716
1
0.94227
0.931074
X
2
0.904383
0.87252
1
0.99112
0.893805
0.932366
X
3
0.916382
1
0.946311
0.945523
0.989839
0.959611
X
4
1
0.881448
0.930161
0.932326
0.889259
0.926639
X
5
0.893087
0.930804
0.930155
0.999411
1
0.950691
X
6
0.871197
0.879595
0.881844
1
0.944697
0.915467According to Table17, we obtain the simulated rank order of relevancy of the influencing variables over time as follows:Year 1:X4>X3>X2>X5>X1>X6.Year 2:X3>X5>X1>X4>X6>X2.Year 3:X2>X1>X3>X4>X5>X6.Year 4:(X1=X6)>X5>X2>X3>X4.Year 5:X5>X3>X6>X1>X2>X4.According to aggregated values, we obtain the final simulated rank order of relevancy of the influencing variables as follows:X3>X5>X2>X1>X4>X6.Since 500 simulations are performed, we obtain the correlation matrix as in Table18.Table 18
Expected correlation matrix obtained after 500 simulations.
DOL(k)/Xj θ=0.5
DOL(1)
DOL(2)
DOL(3)
DOL(4)
DOL(5)
Aggregated value
X
1
0.912041
0.921771
0.925088
0.925480
0.939693
0.924814
X
2
0.928070
0.896937
0.957765
0.913864
0.876159
0.914559
X
3
0.934076
0.940820
0.935899
0.950123
0.953475
0.942879
X
4
0.922037
0.922130
0.911800
0.901827
0.905551
0.912669
X
5
0.890792
0.938745
0.870543
0.895710
0.958414
0.910841
X
6
0.909324
0.925922
0.851819
0.933881
0.937900
0.911769According to Table18, we obtain the expected rank order of relevancy of the influencing variables over time and final expected rank order as follows:Year 1:X3>X2>X4>X1>X6>X5.Year 2:X3>X5>X6>X4>X1>X2.Year 3:X2>X3>X1>X4>X5>X6.Year 4:X3>X6>X1>X2>X4>X5.Year 5:X5>X3>X1>X6>X4>X2.Final:X3>X1>X2>X4>X6>X5.If we compare the final rank order of relevancy of the influencing variables obtained by observation to the final expected rank order of influencing variables obtained by simulation, it can be seen that the revenues retain the greatest relevance. Level of relevancy of the influencing variable related to working days is significantly changed from the second level in the past to the sixth in the future. Influencing variable related to the production rate is increased its relevance from the fourth level in the past to the second in the future, while the rest variables retain approximately the same level of relevancy.Key information obtained from the values of DOL can be summarized in Table19.Table 19
DOL’s indications.
Parameter
DOL value
High
Low
Variable costs (revenues and fixed costs are fixed)
High
Low
Fixed costs (revenues and variable costs are fixed)
Low
High
Revenues up (variable and fixed costs are fixed)
Greater profits
Lower profits
Revenues down (variable and fixed costs are fixed)
Greater losses
Lower lossesFrom Table19, we can see why it matters if value of DOL is high or low. When DOL is high, a change in revenues results in large change in profit or loss. On the other hand, when DOL is low, a change in revenues results in small change in profit or loss.


## 5. Conclusion

Multivariable grey model is used to establish the sequence time response formula which defines the relationship between ore production rate, fixed costs, revenues, production costs, working days, and degree of use of production capacity, as input variables, and Degree of Operating Leverage as indicator of operational efficiency of an existing underground mine. The management of the company obtains the aggregated information vector which sublimates the following important indicators: forecasting model, relative percentage error of forecasting model, average relative percentage error of forecasting model, the grey relational coefficient between inputs and output for defined time interval, and aggregated grey relational coefficient between inputs and output.Large capital intensive projects, such as those in the mineral resource industry, are often associated with diverse sources of both endogenous and exogenous uncertainties. These uncertainties can greatly influence the operational efficiency. Having the ability to plan for these uncertainties is increasingly recognized as critical to long-term mining company success. To decrease uncertainty of future values of the operational efficiency we use expert’s knowledge and simulation processes to find future values of input variables affecting them. Estimation of future revenues is based on the application of mean reversion process, normal and uniform distribution. Geometric Brownian motion is used to define future values of production costs. Values of ore production rate, fixed costs, working days, and degree of use of production capacity are estimated by expert’s knowledge. Simulation results in the average or expected value of Degree of Operating Leverage for every year of defined time horizon. Degree of Operating Leverage is very important indicator to the management of the mining company because it can be used as a good base for Cost-Volume-Profit analysis.The limitation of the developed model is related to the expert’s estimation of the production rate, fixed costs, working days, and degree of use of production capacity as deterministic or crisp values. In order to overcome this problem, we concern the possibility of application of interval or fuzzy numbers as a way to decrease uncertainty of these variables. It increases the complexity of the calculation but makes the model much more realistic. Another one approach of further exploration can be directed to developing submodels finding the future states of these variables and incorporation of them into main forecasting model.



---

*Source: 610307-2015-12-22.xml*