# Evaluation of Underground Zinc Mine Investment Based on Fuzzy-Interval Grey System Theory and Geometric Brownian Motion

Authors: Zoran Gligoric; Lazar Kricak; Cedomir Beljic; Suzana Lutovac; Jelena Milojevic
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

Underground mine projects are often associated with diverse sources of uncertainties. Having the ability to plan for these uncertainties plays a key role in the process of project evaluation and is increasingly recognized as critical to mining project success. To make the best decision, based on the information available, it is necessary to develop an adequate model incorporating the uncertainty of the input parameters. The model is developed on the basis of full discounted cash flow analysis of an underground zinc mine project. The relationships between input variables and economic outcomes are complex and often nonlinear. Fuzzy-interval grey system theory is used to forecast zinc metal prices while geometric Brownian motion is used to forecast operating costs over the time frame of the project. To quantify the uncertainty in the parameters within a project, such as capital investment, ore grade, mill recovery, metal content of concentrate, and discount rate, we have applied the concept of interval numbers. The final decision related to project acceptance is based on the net present value of the cash flows generated by the simulation over the time project horizon.


---

## Body

## 1. Introduction

If we take into consideration that underground mining projects are planned and constructed in an uncertain physical and economic environment, then evaluation of such projects is truly interdisciplinary in nature.Mine investments provide a good example of irreversible investment under uncertainty. Irreversible investment requires more careful analysis because, once the investment takes place, it cannot be recouped without a significant loss of value. Engineering economics is a widely used economic technique for the evaluation of engineering projects. Within it, different methods can be used to make the best decision, that is, whether to accept a project or not.There is a considerable literature dedicated to the problem of mining project evaluation. Samis et al. use Real Options Monte Carlo Simulation to examine the valuation of a multiphase copper-gold project in the presence of a windfall profits tax [1]. Dimitrakopoulos applies the Monte Carlo technique (conditional simulation) to quantify geological uncertainty such as ore grade and tonnage [2]. Topal uses different techniques to estimate the value of the mining project. The major challenge of project evaluation is how to deal with the uncertainty involved in capital investment. Discounted cash flow (DCF) methods, decision trees (DT), Monte Carlo simulation (MCS), and real options (RO) are commonly used for evaluating mining projects [3]. Dessureault et al. use real options pricing as a method for the flexible valuation of a mining project. This paper presents the methods that can be used for the calculation of process and project volatility in operations and provides practical applications from mining operations in USA and Canada [4]. Elkington et al. noted that uncertainty is intrinsic to all mining projects and should be planned for by providing operating and strategic flexibility [5]. Trigeorgis presents a decision-tree model for a mining project in which the present value of the remaining cash flows is uncertain [6]. Samis and Poulin provide a related decision-tree model where mineral price is the underlying source of uncertainty [7, 8].Prior to initialization, a mining project is often evaluated by calculating its net present value (NPV). The NPV is defined as the discounted difference between the expected value of project revenues and costs over the life of the project. The NPV is the preferred criterion of project profitability since it reflects the net contribution to the owner’s equity considering his cost of capital. We propose a simulation approach to incorporate uncertainty in the NPV calculations. Simulation of future zinc metal prices is performed by fuzzy-interval grey system theory. The dynamic nature of operating costs is described by the stochastic process called geometric Brownian motion. In this way, we obtain the probability distribution of operating costs for every year of the project and after that we transform them into adequate interval numbers. The remaining risk factors such as capital investment, ore grade, mill recovery, metal content of concentrate, and discount rate are also quantified by interval numbers using expert knowledge (estimation). When these interval numbers are incorporated in the NPV calculation, we obtain the interval-valued NPV, that is, the project value at risk.The main purpose of this study is to provide an efficient and easy way of strategic decision making, particularly in small underground mining companies. We were motivated by the fact that, in our country, as one of the many developing countries, there are mainly small underground mining companies employing just two or three mining engineers who are responsible for both the production maintenance and strategic planning. In such an environment, mining engineers do not have time to create adequate procedures for decision making, particularly for the decisions influenced by highly volatile parameters such as metal price. There are many stochastic methods for treating the uncertainty of metal prices (e.g., Mean Reversion Process), but if we want to apply them, it is necessary to collect a lot of historical data and run complex regression analysis in order to define the parameters of the simulation process. Interval grey theory can handle problems with unclear information very precisely. Its concept is intuitive and simple to understand for mining engineers. In order to build the forecasting model, only a few data are needed.The proposed model is tested on a hypothetical example, which is similar to many real case studies, and the experiment results verify the rationality and effectiveness of the method.


## 2. Preliminaries of Interval-Valued Differential Equations

### 2.1. Basic Concepts of Fuzzy Set Theory

Various theories exist for describing uncertainty in the modelling of real phenomena and the most popular one is fuzzy set theory [9]. In this paper we applied the concept of the interval-valued possibilistic mean of fuzzy number [10–12].In the classical set theory, an element either belongs or does not belong to a given set. By contrast, in fuzzy set theory, a fuzzy subsetA
~ defined on a universe of discourse X is characterized by a membership function μ
A
~
(
x
)
, which maps each element in A
~ with a real number in the unit interval. Generally, this can be expressed as μ
A
~
(
x
)
:
X
→
[
0,1
], where the value μ
A
~
(
x
) is called the degree of membership of the element x in the fuzzy set A
~. If the universal set X is fixed, a membership function fully determines a fuzzy set. In fuzzy set theory, classical sets are usually called crisp sets.Definition 1.
Leta
1, a
2, and a
3 be real numbers such that a
1
<
a
2
<
a
3. A set A
~ with membership function
(1)
μ
A
~
(
x
)
=
{
x
-
a
1
a
2
-
a
1
,
a
1
≤
x
≤
a
2
x
-
a
3
a
3
-
a
2
,
a
2
≤
x
≤
a
3
0
,
otherwise
is called a fuzzy triangular number and is denoted as A
~
=
(
a
1
,
a
2
,
a
3
). In geometric interpretations, the graph of A
~
(
x
) is a triangle with its base on the interval [
a
1
,
a
3
] and vertex at x
=
a
2.Fuzzy sets can also be represented via theirγ-levels.Definition 2.
Aγ-level set of a fuzzy set A
~ is defined by [
A
~
]
γ
=
{
x
∈
X
∣
μ
A
~
(
x
)
≥
γ
} if γ
>
0 and [
A
~
]
γ
=
cl
⁡
{
x
∈
X
∣
μ
A
~
(
x
)
>
0
} (the closure of the support of A
~) if γ
=
0.In particular, a fuzzy setA
~ is a fuzzy number if and only if the  γ-levels are nested nonempty compact intervals [
A
*
γ
,
A
*
γ
]. This property is the basis for the lower-upper representation of values of the  γ-levels [13]. A fuzzy number A
~ is completely defined by a pair of functions A
*
,
A
*
:
[
0,1
]
→
X, defining the end-points of  γ-levels (
A
*
γ
=
A
*
(
γ
)
;
A
*
γ
=
A
*
(
γ
)
) and satisfying the following conditions:(1)
A
* is a bounded nondecreasing left-continuous function on [
0,1
],
(2)
A
* is a bounded nonincreasing left-continuous function on [
0,1
],
(3)
A
*
(
γ
)
≤
A
*
(
γ
) for all 0
≤
γ
≤
1.According to Dubois and Prade [14], the interval-valued probabilistic mean of a fuzzy number A
~, with  γ levels A
~
γ
=
[
a
γ
,
b
γ
]
,
γ
∈
[
0,1
] (see Figure 1) is the interval E
(
A
~
)
=
[
E
*
(
A
~
)
,
E
*
(
A
~
)
], where
(2)
E
*
(
A
~
)
=
∫
0
1
a
γ
d
γ
,
E
*
(
A
~
)
=
∫
0
1
b
γ
d
γ
.Carlsson and Fuller [11] introduced the interval-valued possibilistic mean of a fuzzy number A
~ as the interval M
(
A
~
)
=
[
M
*
(
A
~
)
,
M
*
(
A
~
)
]. First, we note that from the equality
(3)
M
-
(
A
~
)
∶
=
∫
0
1
γ
(
a
γ
+
b
γ
)
d
γ
=
∫
0
1
γ
(
(
a
γ
+
b
γ
)
/
2
)
d
γ
∫
0
1
γ
d
γ
it follows that M
-
(
A
~
) is nothing else but the level-weighted average of the arithmetic means of all  γ-sets; that is, the weight of the arithmetic mean of a
γ and b
γ is just  γ. Second, we can rewrite M
-
(
A
~
) as
(4)
M
-
(
A
~
)
∶
=
∫
0
1
γ
(
a
γ
+
b
γ
)
d
γ
=
2
∫
0
1
γ
a
γ
d
γ
+
2
∫
0
1
γ
b
γ
d
γ
2
=
1
2
(
∫
0
1
γ
a
γ
d
γ
1
/
2
+
∫
0
1
γ
b
γ
d
γ
1
/
2
)
=
1
2
(
∫
0
1
γ
a
γ
d
γ
∫
0
1
γ
d
γ
+
∫
0
1
γ
b
γ
d
γ
∫
0
1
γ
d
γ
)
.
Third, let us take a closer look at the right-hand side of the equation for M
-
(
A
~
). The first quantity, denoted by M
*
(
A
~
), can be reformulated as
(5)
M
*
(
A
~
)
=
2
∫
0
1
γ
a
γ
d
γ
=
∫
0
1
γ
a
γ
d
γ
∫
0
1
γ
d
γ
=
∫
0
1
Pos
[
A
≤
a
γ
]
a
γ
d
γ
∫
0
1
Pos
[
A
≤
a
γ
]
d
γ
=
∫
0
1
Pos
[
A
≤
a
γ
]
×
min
⁡
[
A
]
γ
d
γ
∫
0
1
Pos
[
A
≤
a
γ
]
d
γ
,
where Pos denotes possibility; that is,
(6)
Pos
[
A
~
≤
a
γ
]
=
∏
(
(
-
∞
,
a
γ
]
)
=
sup
⁡
︸
u
≤
a
γ
A
~
(
u
)
=
γ
.
So M
*
(
A
~
) is nothing else but the lower possibility-weighted average of the minima of the  γ-sets, and this is why we call it the lower possibilistic mean value of A
~.Figure 1
Fuzzy number andγ-level cut.In a similar manner, we introduceM
*
(
A
~
), the upper possibilistic mean value of A
~:
(7)
M
*
(
A
~
)
=
2
∫
0
1
γ
b
γ
d
γ
=
∫
0
1
γ
b
γ
d
γ
∫
0
1
γ
d
γ
=
∫
0
1
Pos
[
A
≥
b
γ
]
b
γ
d
γ
∫
0
1
Pos
[
A
≥
b
γ
]
d
γ
=
∫
0
1
Pos
[
A
≥
b
γ
]
×
max
⁡
[
A
]
γ
d
γ
∫
0
1
Pos
[
A
≥
b
γ
]
d
γ
,
where we have used equality
(8)
Pos
[
A
~
≥
b
γ
]
=
∏
(
[
b
γ

,

∞
)
)
=
sup
⁡
︸
u
≥
b
γ
A
~
(
u
)
=
γ
.
The lower possibilistic mean M
*
(
A
~
) is the weighted average of the minima of the  γ-levels of A
~. Similarly, the upper possibilistic mean M
*
(
A
~
) is the weighted average of the maxima of the  γ-levels of A
~. According to Carlsson, the fuzzy number can now be expressed as follows:
(9)
M
*
(
A
~
)
=
2
∫
0
1
γ
a
γ
d
γ
,
M
*
(
A
~
)
=
2
∫
0
1
γ
b
γ
d
γ
.Definition 3.
LetA
~
=
(
a
2
,
α
,
β
) be a triangular fuzzy number with centre a
2, left-width α
>
0, and right-width β
>
0 (see Figure 1); then a
γ and b
γ are computed as follows:
(10)
1
:
γ
=
α
:
(
a
γ
-
(
a
2
-
α
)
)
⟶
a
γ
=
a
2
-
(
1
-
γ
)
α
,
1
:
γ
=
β
:
(
(
a
2
)
-
b
γ
)
⟶
b
γ
=
a
2
+
(
1
-
γ
)
β
.
Now, theγ-level of A
~ is computed by
(11)
[
A
~
]
γ
=
[
a
2
-
(
1
-
γ
)
α
,
a
2
+
(
1
-
γ
)
β
]
,
∀
γ
∈
[
0,1
]
.
That is,
(12)
M
*
(
A
~
)
=
2
∫
0
1
γ
[
a
2
-
(
1
-
γ
)
]
α
d
γ
=
a
2
-
α
3
,
M
*
(
A
~
)
=
2
∫
0
1
γ
[
a
2
+
(
1
-
γ
)
]
β
d
γ
=
a
2
+
β
3
,
and therefore
(13)
M
(
A
~
)
=
[
M
*

(
A
~
)
,
M
*

(
A
~
)
]
=
[
a
2
-
α
3
,
a
2
+
β
3
]
.
Obviously, M
(
A
~
) is a closed interval bounded by the lower and upper possibilistic mean values of A
~. The crisp possibilistic mean value of A
~ is defined as the arithmetic mean of its lower possibilistic and upper possibilistic mean values; that is,
(14)
M
-
(
A
~
)
=
M
*

(
A
~
)
+
M
*

(
A
~
)
2
.
According to the above way of transformation (see (13)), we obtain an interval number having both a lower bound a
l and upper bound a
u, A
=
[
a
l
,
a
u
], where a
l
=
a
2
-
α
/
3, and a
u
=
a
2
+
β
/
3.From interval arithmetic, the following operations of interval numbers are defined as follows.Definition 4.
For∀
k
>
0 and a given interval number A
=
[
a
l
,
a
u
]
,
a
l
,
a
u
∈
R,
(15)
k
×
[
a
l
,
a
u
]
=
[
k
×
a
l
,
k
×
a
u
]
and for ∀
k
<
0,
(16)
-
k
×
[
a
l
,
a
u
]
=
[
-
k
×
a
u
,
-
k
×
a
l
]
.Definition 5.
For any two interval numbers[
a
l
,
a
u
] and [
b
l
,
b
u
],
(17)
[
a
l
,
a
u
]
+
[
b
l
,
b
u
]
=
[
a
l
+
b
l
,
a
u
+
b
u
]
,
[
a
l
,
a
u
]
-
[
b
l
,
b
u
]
=
[
a
l
-
b
u
,
a
u
-
b
l
]
.Definition 6.
For any two interval numbers[
a
l
,
a
u
] and [
b
l
,
b
u
], the multiplication is defined as follows:
(18)
[
a
l
,
a
u
]
×
[
b
l
,
b
u
]
=
[
min
⁡
(
a
l
b
l
,
a
l
b
u
,
a
u
b
l
,
a
u
b
u
)
,
max
⁡
⁡
(
a
l
b
l
,
a
l
b
u
,
a
u
b
l
,
a
u
b
u
)
]
.Definition 7.
For any two interval numbers[
a
l
,
a
u
] and [
b
l
,
b
u
], the division is defined as follows:
(19)
[
a
l
,
a
u
]
÷
[
b
l
,
b
u
]
=
[
min
⁡
(
a
l
b
l
,
a
l
b
u
,
a
u
b
l
,
a
u
b
u
)
,
max
⁡
(
a
l
b
l
,
a
l
b
u
,
a
u
b
l
,
a
u
b
u
)
]
,
0
∉
[
b
l
,
b
u
]
.


### 2.2. Interval-Valued Differential Equations

In this section we consider an interval-valued differential equation of the following form:(20)
y
˙
=
f
(
t
,
y
(
t
)
)
,
y
(
t
0
)
=
y
0
,
where f
:
[
a
,
b
]
×
E
→
E with f
(
t
,
y
(
t
)
)
=
[
f
l
(
t
,
y
(
t
)
)
,
f
u
(
x
,
y
(
t
)
)
] for y
(
t
)
∈
E
,
y
(
t
) =
[
y
l
(
t
)
,
y
u
(
t
)
]
,
y
0
=
[
y
0
l
,
y
0
u
]. Note that we consider only Hukuhara differentiable solutions; that is, there exists δ
>
0 such that there are no switching points in [
t
0
,
t
0
+
δ
] [15, 16].Definition 8.
Letf
:
[
a
,
b
]
→
E be Hukuhara differentiable at t
0
∈
]
a
,
b
[. We say that f is (i)-Hukuhara differentiable at t
0 if
(21)
f
˙
(
t
0
)
=
[
f
˙
l
(
t
0
)
,
f
˙
u
(
t
0
)
]
and that f is (ii)-Hukuhara differentiable at t
0 if
(22)
f
˙
(
t
0
)
=
[
f
˙
u
(
t
0
)
,
f
˙
l
(
t
0
)
]
.
The solution of the differential equation (20) depends on the choice of the Hukuhara derivative ((i) or (ii)). To solve the interval-valued differential equation it is necessary to reduce the interval-valued differential equation to a system of ordinary differential equations [17–20].Let[
y
(
t
)
]
=
[
y
l
(
t
)
,
y
u
(
t
)
]. If y
(
t
) is (i)-Hukuhara differentiable, then D
1
y
(
t
)
=
[
y
˙
l
(
t
)
,
y
˙
u
(
t
)
] transforms (20) into the following system of ordinary differential equations:
(23)
y
˙
l
(
t
)
=
f
l
(
t
,
y
(
t
)
)
,
y
l
(
t
0
)
=
y
0
l
,
y
˙
u
(
t
)
=
f
u
(
t
,
y
(
t
)
)
,
y
u
(
t
0
)
=
y
0
u
.
Also, if y
(
t
) is (ii)-Hukuhara differentiable, then D
2
y
(
t
)
=
[
y
˙
u
(
t
)
,
y
˙
l
(
t
)
] transforms (20) into the following system of ordinary differential equations:
(24)
y
˙
l
(
t
)
=
f
u
(
t
,
y
(
t
)
)
,
y
l
(
t
0
)
=
y
0
l
,
y
˙
u
(
t
)
=
f
l
(
t
,
y
(
t
)
)
,
y
u
(
t
0
)
=
y
0
u
,
where f
(
t
,
y
(
t
)
)
=
[
f
l
(
t
,
y
(
t
)
)
,
f
u
(
t
,
y
(
t
)
)
].


## 2.1. Basic Concepts of Fuzzy Set Theory

Various theories exist for describing uncertainty in the modelling of real phenomena and the most popular one is fuzzy set theory [9]. In this paper we applied the concept of the interval-valued possibilistic mean of fuzzy number [10–12].In the classical set theory, an element either belongs or does not belong to a given set. By contrast, in fuzzy set theory, a fuzzy subsetA
~ defined on a universe of discourse X is characterized by a membership function μ
A
~
(
x
)
, which maps each element in A
~ with a real number in the unit interval. Generally, this can be expressed as μ
A
~
(
x
)
:
X
→
[
0,1
], where the value μ
A
~
(
x
) is called the degree of membership of the element x in the fuzzy set A
~. If the universal set X is fixed, a membership function fully determines a fuzzy set. In fuzzy set theory, classical sets are usually called crisp sets.Definition 1.
Leta
1, a
2, and a
3 be real numbers such that a
1
<
a
2
<
a
3. A set A
~ with membership function
(1)
μ
A
~
(
x
)
=
{
x
-
a
1
a
2
-
a
1
,
a
1
≤
x
≤
a
2
x
-
a
3
a
3
-
a
2
,
a
2
≤
x
≤
a
3
0
,
otherwise
is called a fuzzy triangular number and is denoted as A
~
=
(
a
1
,
a
2
,
a
3
). In geometric interpretations, the graph of A
~
(
x
) is a triangle with its base on the interval [
a
1
,
a
3
] and vertex at x
=
a
2.Fuzzy sets can also be represented via theirγ-levels.Definition 2.
Aγ-level set of a fuzzy set A
~ is defined by [
A
~
]
γ
=
{
x
∈
X
∣
μ
A
~
(
x
)
≥
γ
} if γ
>
0 and [
A
~
]
γ
=
cl
⁡
{
x
∈
X
∣
μ
A
~
(
x
)
>
0
} (the closure of the support of A
~) if γ
=
0.In particular, a fuzzy setA
~ is a fuzzy number if and only if the  γ-levels are nested nonempty compact intervals [
A
*
γ
,
A
*
γ
]. This property is the basis for the lower-upper representation of values of the  γ-levels [13]. A fuzzy number A
~ is completely defined by a pair of functions A
*
,
A
*
:
[
0,1
]
→
X, defining the end-points of  γ-levels (
A
*
γ
=
A
*
(
γ
)
;
A
*
γ
=
A
*
(
γ
)
) and satisfying the following conditions:(1)
A
* is a bounded nondecreasing left-continuous function on [
0,1
],
(2)
A
* is a bounded nonincreasing left-continuous function on [
0,1
],
(3)
A
*
(
γ
)
≤
A
*
(
γ
) for all 0
≤
γ
≤
1.According to Dubois and Prade [14], the interval-valued probabilistic mean of a fuzzy number A
~, with  γ levels A
~
γ
=
[
a
γ
,
b
γ
]
,
γ
∈
[
0,1
] (see Figure 1) is the interval E
(
A
~
)
=
[
E
*
(
A
~
)
,
E
*
(
A
~
)
], where
(2)
E
*
(
A
~
)
=
∫
0
1
a
γ
d
γ
,
E
*
(
A
~
)
=
∫
0
1
b
γ
d
γ
.Carlsson and Fuller [11] introduced the interval-valued possibilistic mean of a fuzzy number A
~ as the interval M
(
A
~
)
=
[
M
*
(
A
~
)
,
M
*
(
A
~
)
]. First, we note that from the equality
(3)
M
-
(
A
~
)
∶
=
∫
0
1
γ
(
a
γ
+
b
γ
)
d
γ
=
∫
0
1
γ
(
(
a
γ
+
b
γ
)
/
2
)
d
γ
∫
0
1
γ
d
γ
it follows that M
-
(
A
~
) is nothing else but the level-weighted average of the arithmetic means of all  γ-sets; that is, the weight of the arithmetic mean of a
γ and b
γ is just  γ. Second, we can rewrite M
-
(
A
~
) as
(4)
M
-
(
A
~
)
∶
=
∫
0
1
γ
(
a
γ
+
b
γ
)
d
γ
=
2
∫
0
1
γ
a
γ
d
γ
+
2
∫
0
1
γ
b
γ
d
γ
2
=
1
2
(
∫
0
1
γ
a
γ
d
γ
1
/
2
+
∫
0
1
γ
b
γ
d
γ
1
/
2
)
=
1
2
(
∫
0
1
γ
a
γ
d
γ
∫
0
1
γ
d
γ
+
∫
0
1
γ
b
γ
d
γ
∫
0
1
γ
d
γ
)
.
Third, let us take a closer look at the right-hand side of the equation for M
-
(
A
~
). The first quantity, denoted by M
*
(
A
~
), can be reformulated as
(5)
M
*
(
A
~
)
=
2
∫
0
1
γ
a
γ
d
γ
=
∫
0
1
γ
a
γ
d
γ
∫
0
1
γ
d
γ
=
∫
0
1
Pos
[
A
≤
a
γ
]
a
γ
d
γ
∫
0
1
Pos
[
A
≤
a
γ
]
d
γ
=
∫
0
1
Pos
[
A
≤
a
γ
]
×
min
⁡
[
A
]
γ
d
γ
∫
0
1
Pos
[
A
≤
a
γ
]
d
γ
,
where Pos denotes possibility; that is,
(6)
Pos
[
A
~
≤
a
γ
]
=
∏
(
(
-
∞
,
a
γ
]
)
=
sup
⁡
︸
u
≤
a
γ
A
~
(
u
)
=
γ
.
So M
*
(
A
~
) is nothing else but the lower possibility-weighted average of the minima of the  γ-sets, and this is why we call it the lower possibilistic mean value of A
~.Figure 1
Fuzzy number andγ-level cut.In a similar manner, we introduceM
*
(
A
~
), the upper possibilistic mean value of A
~:
(7)
M
*
(
A
~
)
=
2
∫
0
1
γ
b
γ
d
γ
=
∫
0
1
γ
b
γ
d
γ
∫
0
1
γ
d
γ
=
∫
0
1
Pos
[
A
≥
b
γ
]
b
γ
d
γ
∫
0
1
Pos
[
A
≥
b
γ
]
d
γ
=
∫
0
1
Pos
[
A
≥
b
γ
]
×
max
⁡
[
A
]
γ
d
γ
∫
0
1
Pos
[
A
≥
b
γ
]
d
γ
,
where we have used equality
(8)
Pos
[
A
~
≥
b
γ
]
=
∏
(
[
b
γ

,

∞
)
)
=
sup
⁡
︸
u
≥
b
γ
A
~
(
u
)
=
γ
.
The lower possibilistic mean M
*
(
A
~
) is the weighted average of the minima of the  γ-levels of A
~. Similarly, the upper possibilistic mean M
*
(
A
~
) is the weighted average of the maxima of the  γ-levels of A
~. According to Carlsson, the fuzzy number can now be expressed as follows:
(9)
M
*
(
A
~
)
=
2
∫
0
1
γ
a
γ
d
γ
,
M
*
(
A
~
)
=
2
∫
0
1
γ
b
γ
d
γ
.Definition 3.
LetA
~
=
(
a
2
,
α
,
β
) be a triangular fuzzy number with centre a
2, left-width α
>
0, and right-width β
>
0 (see Figure 1); then a
γ and b
γ are computed as follows:
(10)
1
:
γ
=
α
:
(
a
γ
-
(
a
2
-
α
)
)
⟶
a
γ
=
a
2
-
(
1
-
γ
)
α
,
1
:
γ
=
β
:
(
(
a
2
)
-
b
γ
)
⟶
b
γ
=
a
2
+
(
1
-
γ
)
β
.
Now, theγ-level of A
~ is computed by
(11)
[
A
~
]
γ
=
[
a
2
-
(
1
-
γ
)
α
,
a
2
+
(
1
-
γ
)
β
]
,
∀
γ
∈
[
0,1
]
.
That is,
(12)
M
*
(
A
~
)
=
2
∫
0
1
γ
[
a
2
-
(
1
-
γ
)
]
α
d
γ
=
a
2
-
α
3
,
M
*
(
A
~
)
=
2
∫
0
1
γ
[
a
2
+
(
1
-
γ
)
]
β
d
γ
=
a
2
+
β
3
,
and therefore
(13)
M
(
A
~
)
=
[
M
*

(
A
~
)
,
M
*

(
A
~
)
]
=
[
a
2
-
α
3
,
a
2
+
β
3
]
.
Obviously, M
(
A
~
) is a closed interval bounded by the lower and upper possibilistic mean values of A
~. The crisp possibilistic mean value of A
~ is defined as the arithmetic mean of its lower possibilistic and upper possibilistic mean values; that is,
(14)
M
-
(
A
~
)
=
M
*

(
A
~
)
+
M
*

(
A
~
)
2
.
According to the above way of transformation (see (13)), we obtain an interval number having both a lower bound a
l and upper bound a
u, A
=
[
a
l
,
a
u
], where a
l
=
a
2
-
α
/
3, and a
u
=
a
2
+
β
/
3.From interval arithmetic, the following operations of interval numbers are defined as follows.Definition 4.
For∀
k
>
0 and a given interval number A
=
[
a
l
,
a
u
]
,
a
l
,
a
u
∈
R,
(15)
k
×
[
a
l
,
a
u
]
=
[
k
×
a
l
,
k
×
a
u
]
and for ∀
k
<
0,
(16)
-
k
×
[
a
l
,
a
u
]
=
[
-
k
×
a
u
,
-
k
×
a
l
]
.Definition 5.
For any two interval numbers[
a
l
,
a
u
] and [
b
l
,
b
u
],
(17)
[
a
l
,
a
u
]
+
[
b
l
,
b
u
]
=
[
a
l
+
b
l
,
a
u
+
b
u
]
,
[
a
l
,
a
u
]
-
[
b
l
,
b
u
]
=
[
a
l
-
b
u
,
a
u
-
b
l
]
.Definition 6.
For any two interval numbers[
a
l
,
a
u
] and [
b
l
,
b
u
], the multiplication is defined as follows:
(18)
[
a
l
,
a
u
]
×
[
b
l
,
b
u
]
=
[
min
⁡
(
a
l
b
l
,
a
l
b
u
,
a
u
b
l
,
a
u
b
u
)
,
max
⁡
⁡
(
a
l
b
l
,
a
l
b
u
,
a
u
b
l
,
a
u
b
u
)
]
.Definition 7.
For any two interval numbers[
a
l
,
a
u
] and [
b
l
,
b
u
], the division is defined as follows:
(19)
[
a
l
,
a
u
]
÷
[
b
l
,
b
u
]
=
[
min
⁡
(
a
l
b
l
,
a
l
b
u
,
a
u
b
l
,
a
u
b
u
)
,
max
⁡
(
a
l
b
l
,
a
l
b
u
,
a
u
b
l
,
a
u
b
u
)
]
,
0
∉
[
b
l
,
b
u
]
.


## 2.2. Interval-Valued Differential Equations

In this section we consider an interval-valued differential equation of the following form:(20)
y
˙
=
f
(
t
,
y
(
t
)
)
,
y
(
t
0
)
=
y
0
,
where f
:
[
a
,
b
]
×
E
→
E with f
(
t
,
y
(
t
)
)
=
[
f
l
(
t
,
y
(
t
)
)
,
f
u
(
x
,
y
(
t
)
)
] for y
(
t
)
∈
E
,
y
(
t
) =
[
y
l
(
t
)
,
y
u
(
t
)
]
,
y
0
=
[
y
0
l
,
y
0
u
]. Note that we consider only Hukuhara differentiable solutions; that is, there exists δ
>
0 such that there are no switching points in [
t
0
,
t
0
+
δ
] [15, 16].Definition 8.
Letf
:
[
a
,
b
]
→
E be Hukuhara differentiable at t
0
∈
]
a
,
b
[. We say that f is (i)-Hukuhara differentiable at t
0 if
(21)
f
˙
(
t
0
)
=
[
f
˙
l
(
t
0
)
,
f
˙
u
(
t
0
)
]
and that f is (ii)-Hukuhara differentiable at t
0 if
(22)
f
˙
(
t
0
)
=
[
f
˙
u
(
t
0
)
,
f
˙
l
(
t
0
)
]
.
The solution of the differential equation (20) depends on the choice of the Hukuhara derivative ((i) or (ii)). To solve the interval-valued differential equation it is necessary to reduce the interval-valued differential equation to a system of ordinary differential equations [17–20].Let[
y
(
t
)
]
=
[
y
l
(
t
)
,
y
u
(
t
)
]. If y
(
t
) is (i)-Hukuhara differentiable, then D
1
y
(
t
)
=
[
y
˙
l
(
t
)
,
y
˙
u
(
t
)
] transforms (20) into the following system of ordinary differential equations:
(23)
y
˙
l
(
t
)
=
f
l
(
t
,
y
(
t
)
)
,
y
l
(
t
0
)
=
y
0
l
,
y
˙
u
(
t
)
=
f
u
(
t
,
y
(
t
)
)
,
y
u
(
t
0
)
=
y
0
u
.
Also, if y
(
t
) is (ii)-Hukuhara differentiable, then D
2
y
(
t
)
=
[
y
˙
u
(
t
)
,
y
˙
l
(
t
)
] transforms (20) into the following system of ordinary differential equations:
(24)
y
˙
l
(
t
)
=
f
u
(
t
,
y
(
t
)
)
,
y
l
(
t
0
)
=
y
0
l
,
y
˙
u
(
t
)
=
f
l
(
t
,
y
(
t
)
)
,
y
u
(
t
0
)
=
y
0
u
,
where f
(
t
,
y
(
t
)
)
=
[
f
l
(
t
,
y
(
t
)
)
,
f
u
(
t
,
y
(
t
)
)
].


## 3. Model of the Evaluation

### 3.1. The Concept of Evaluation

Economic evaluation of a mine project requires estimation of the revenues and costs throughout the defined lifetime of the mine. Such evaluation can be treated as strategic decision making under multiple sources of uncertainties. Therefore, to make the best decision, based on the information available, it is necessary to develop an adequate model incorporating the uncertainty of the input parameters. The model should be able to involve a common time horizon, taking the characteristics of the input variables that directly affect the value of the proposed project.The model is developed on the basis of full discounted cash flow analysis of an underground zinc mine project. The operating discounted cash flows are usually estimated on an annual basis. Net present value of investment is used as a key criterion in the process of mine project estimation. The expected net present value of the project is a function of the variables as(25)
E
(
NPV
∣
Q
,
P
,
G
,
M
,
C
,
I
,
r
,
t
)
≥
0
,
where Q denotes the production rate (capacity); P denotes the zinc metal price; G is the grade; M is mill recovery; C is operating costs; I is capital investment; r is discount rate, and t is the lifetime of the project; that is, the period in which the cash flow is generated.In this paper, we treat in detail only the variability of metal prices and operating costs, without intending to decrease the significance of the remaining parameters. These parameters are taken into account on the basis of expert knowledge (estimation).


### 3.2. Forecasting the Revenue of the Mine

Most mining companies realize their revenues by selling metal concentrates as a final product. Estimating mineral project revenue 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 value of the metal concentrate can be expressed as follows:(26)
V
con
=
P
·
(
m
con
-
m
mr
)
,
where P is metal price ($
/
t), m
con is metal content of concentrate (%), m
mr is metal recovery ratio (%), and
(27)
m
con
-
m
mr
=
{
(
m
con
%
-
8
)
·
100
m
con
%
≤
85
%
;
m
con
%
-
8
(
m
con
%
-
8
)
·
100
m
con
%
>
85
%
;
85
%
.
Annual mine revenue is calculated according to the following equation:
(28)
R
year
=
Q
year
·
V
con
·
G
·
M
m
con
,
where Q
year is annual ore production (t
/
year), Gis grade of the ore mined (%), M is mill recovery (%).Annual ore production is derived from the mining project schedule and is defined as crisp value. The concept of grade(
G
) is defined as the ratio of useful mass of metal to the total mass of ore and its critical value fluctuates over deposit space and can be estimated by experts and defined as interval number G
=
[
G
l
,
G
u
]. Mill recovery (
M
) is related to the flotation as the most widely used method for the concentration of fine grained minerals. It can also be defined as interval number M
=
[
M
l
,
M
u
]. Metal content represents the quality of concentrate and we also apply the concept of an interval number to define it:
m
con
=
[
m
con
l
,
m
con
u
].The major external source of risk affecting mine revenue is related to the uncertainty about market behaviour of metal prices. Forecasting the precise future state of the metal price is a very difficult task for mine planners. To predict future metal prices, we apply the concept based on the transformation of historical metal prices into adequate fuzzy-interval numbers and grey system theory.The forecasting model of metal prices is composed of the following steps.Step 1.
Create the setPDF
=
{
pdf
j
}
,
j
=
1,2
,
…
,
N, where pdf
j is the probability density function of metal prices for every historical year. The minimum number of elements of the set is four: j
min
⁡
=
1,2
,
3,4.Step 2.
Transform the set PDF into the setTFN
=
{
(
a
j
,
b
j
,
c
j
)
}, where (
a
j
,
b
j
,
c
j
) is an adequate fuzzy triangular number.Step 3.
Transform the set TFN into the setINT
=
{
[
a
j
l
,
a
j
u
]
}, where [
a
j
l
,
a
j
u
] is an adequate interval number.Step 4.
Using grey system prediction theory, create a grey differential equation of typeG
M
(
1,1
), that is, the first-order variable grey derivative.Step 5.
Testing ofG
M
(
1,1
) by residual error testing and the posterior error detection method.


#### 3.2.1. Analysis of Historical Metal Prices

For every historical year it is necessary to define a probability density function with the following characteristics: shape by histogram, mean value(
µ
j
), and standard deviation (
σ
j
). In this way, we obtain the sequence of probability density functions of P; P
j
~
(
pdf
j
,
µ
j
,
σ
j
)
,
j
=
1,2
,
…
,
N, where N is the total number of historical years.


#### 3.2.2. Fuzzification of Metal Prices

The sequence of obtainedpdf
j of P
j can be transformed into a sequence of triangular fuzzy numbers of P
j; P
j
~
TFN
j, j
=
1,2
,
…
,
N; that is, P
1
~
pdf
1
→
P
1
~
TFN
1; P
2
~
pdf
2
→
P
2
~
TFN
2
;
…; P
N
~
pdf
N
→
P
N
~
TFN
N. The method of transformation is based on the following facts: the support of the membership function and the pdf are the same, and the point with the highest probability (likelihood) has the highest possibility. For more details, see Swishchuk et al. [21]. The uncertainty in the P parameter is modelled by a triangular fuzzy number with the membership function which has the support of µ
j
-
2
σ
j
<
P
j
<
µ
j
+
2
σ
j
,
j
=
1,2
,
…
,
N, set up for around 95% confidence interval of distribution function. If we take into consideration that the triangular fuzzy number is defined as a triplet (
a
1
j
,
a
2
j
,
a
3
j
), then a
1
j and a
3
j are the lower bound and upper bound obtained from the lower and upper bound of 5% of the distribution, and the most promising value a
2
j is equal to the mean value of the distribution. For more details, see Do et al. [22].


#### 3.2.3. Metal Prices as Interval Numbers

The sequence of obtainedTFN
j of P
j is transformed into a sequence of interval numbers of P
j
;
P
j
~
IN
T
j
,
j
=
1,2
,
…
,
N; that is, P
1
~
TF
N
1
→
P
1
~
IN
T
1
;
P
2
~
TF
N
2
→
P
2
~
IN
T
2
;
…;P
N
~
TF
N
N
→
P
N
~
IN
T
N. The method of transformation is described in Section 2.1 (basic concepts of fuzzy set theory). According to this method of transformation, we obtain set
(29)
P
j
~
IN
T
j
;
P
j
=
{
[
P
j
l
,
P
j
u
]
}
=
{
[
a
2
j
-
α
j
3
,
a
2
j
+
β
j
3
]
}
,
j
=
1,2
,
…
,
N
.


#### 3.2.4. Forecasting Model of Metal Prices

The grey model is a powerful tool for forecasting the behaviour of the system in the future. It has been successfully applied to various fields since it was proposed by Deng [23–31]. In this paper we use a one-variable first-order differential grey equation, G
M
(
1,1
). The essence of G
M
(
1,1
) is to accumulate the original data (historical metal prices) in order to obtain regular data. By setting up the grey differential equation, we obtain the fitted curve in order to predict the future states of the system.Definition 9.
Assume that(30)
P
(
0
)
(
t
j
)
=
{
[
P
(
0
)
l
(
t
j
)
,
P
(
0
)
u
(
t
j
)
]
}
=
{
[
P
(
0
)
l
(
t
1
)
,
P
(
0
)
u
(
t
1
)
]
,
[
P
(
0
)
l
(
t
2
)
,
P
(
0
)
u
(
t
2
)
]
,
…
,
[
P
(
0
)
l
(
t
N
)
,
P
(
0
)
u
(
t
N
)
]
}
is the original series of interval metal prices obtained by transformation. Sampling interval isΔ
t
j
=
t
j
-
t
j
-
1
=
1
year.Definition 10.
LetP
(
1
)
(
t
j
)
=
{
[
P
(
1
)
l
(
t
j
)
,
P
(
1
)
u
(
t
j
)
]
} be a new sequence generated by the accumulated generating operation (AGO), where
(31)
P
(
1
)
l
(
t
j
)
=
∑
j
=
1
N
P
(
0
)
l
(
t
j
)
,
P
(
1
)
u
(
t
j
)
=
∑
j
=
1
N
P
(
0
)
u
(
t
j
)
.
In the process of forecasting metal prices, P
(
1
)
l
(
t
j
) and P
(
1
)
u
(
t
j
) are the solutions of the following grey differential equation:
(32)
d
[
P
(
1
)
l
(
t
)
,
P
(
1
)
u
(
t
)
]
d
t
+
[
q
l
,
q
u
]
·
[
P
(
1
)
l
(
t
)
,
P
(
1
)
u
(
t
)
]
=
[
w
l
,
w
u
]
.
Obviously, (32) is an interval-valued differential equation (see (20)).
To get the values of parametersq
=
[
q
l
,
q
u
] and w
=
[
w
l
,
w
u
], the least square method is used as follows:
(33)
⊗
[
q
w
]
=
⊗
(
[
B
T
B
]
-
1
B
T
Y
(
0
)
)
,
where
(34)
⊗
B
=
[
-
1
2
(
P
(
0
)
(
1
)
+
P
(
0
)
(
2
)
)
1
-
1
2
(
P
(
0
)
(
2
)
+
P
(
0
)
(
3
)
)
⋮
1
⋮
-
1
2
(
P
(
0
)
(
N
-
1
)
+
P
(
0
)
(
N
)
)
1
]
⊗
Y
(
0
)
=
[
P
(
0
)
(
2
)
P
(
0
)
(
3
)
⋯
P
(
0
)
(
N
)





]
T
.
Note that the sign ⊗ indicates the interval number.
IfP
˙
(
1
)
(
t
) is considered as (i)-Hukuhara differentiable, then the following system of ordinary differential equations is as follows:
(35)
P
˙
(
1
)
l
(
t
)
=
w
l
-
q
u
·
P
(
1
)
u
(
t
)
,
P
l
(
t
0
)
=
P
(
0
)
l
(
1
)
,
P
˙
(
1
)
u
(
t
)
=
w
u
-
q
l
·
P
(
1
)
l
(
t
)
,
P
u
(
t
0
)
=
P
(
0
)
u
(
1
)
.
The solution of this system (forecasted equation) is as follows:
(36)
P
^
(
1
)
l
(
t
j
+
1
)
=
1
2
·
e
q
l
·
q
u
·
t
·
Z
1
+
1
2
·
e
-
q
l
·
q
u
·
t
·
Z
2
+
w
u
q
l
,
j
=
1,2
,
…
,
N
,
P
^
(
1
)
u
(
t
j
+
1
)
=
1
2
·
e
q
l
·
q
u
·
t
·
Z
3
+
1
2
·
e
-
q
l
·
q
u
·
t
·
Z
4
+
w
l
q
u
,
j
=
1,2
,
…
,
N
,
where
(37)
Z
1
=
P
(
0
)
l
(
1
)
-
q
u
·
P
(
0
)
u
(
1
)
q
l
·
q
u
+
w
l
q
l
·
q
u
-
w
u
q
l
,
Z
2
=
P
(
0
)
l
(
1
)
+
q
u
·
P
(
0
)
u
(
1
)
q
l
·
q
u
-
w
l
q
l
·
q
u
-
w
u
q
l
,
Z
3
=
P
(
0
)
u
(
1
)
-
P
(
0
)
l
(
1
)
·
q
l
·
q
u
q
u
+
w
u
q
l
·
q
u
-
w
l
q
u
,
Z
4
=
P
(
0
)
u
(
1
)
+
P
(
0
)
l
(
1
)
·
q
l
·
q
u
q
u
-
w
u
q
l
·
q
u
-
w
l
q
u
.
To obtain the forecasted value of the primitive (original) metal price data at time (
t
j
+
1
), the inverse accumulated generating operation (IAGO) is used as follows:
(38)
⊗
P
^
(
0
)
(
t
1
)
=
⊗
P
(
0
)
(
t
1
)
⊗
P
^
(
0
)
(
t
j
+
1
)
=
⊗
P
^
(
1
)
(
t
j
+
1
)
-
⊗
P
^
(
1
)
(
t
j
)
,
j
=
1,2
,
…
,
N
.


#### 3.2.5. The Model Accuracy

Definition 11.
For a given interval grey number⊗
(
[
a
l
,
a
u
]
), it is common to take a whitening value ⊗
~
(
[
a
l
,
a
u
]
)
=
ω
·
a
l
+
(
1
-
ω
)
·
a
u for ∀
ω
∈
[
0,1
]. Furthermore, if ω
=
0.5, it is called equal weight average whitenization [23].
Residual error testing is composed of the calculation of relative error and absolute error between⊗
~
P
(
0
)
(
t
j
) and ⊗
~
P
^
(
0
)
(
t
j
) based on the following formulas:
(39)
Δ
ε
(
t
j
)
=
⊗
~
P
(
0
)
(
t
j
)
-
⊗
~
P
^
(
0
)
(
t
j
)
,
|
Δ
ε
(
t
j
)
|
=
|
⊗
~
P
(
0
)
(
t
j
)
-
⊗
~
P
^
(
0
)
(
t
j
)
|
.
The posterior error detection method means calculation of the standard deviation of original metal price series (S
1) and standard deviation of absolute error (S
2):
(40)
S
1
=
∑
j
=
1
N
(
⊗
~
P
(
0
)
(
t
j
)
-
⊗
~
P
-
(
0
)
(
t
j
)
)
2
N
-
1
,
S
2
=
∑
j
=
1
N
(
Δ
ε
(
t
j
)
-
Δ
ε
-
(
t
j
)
)
2
N
-
1
.
The variance ratio is equal to
(41)
C
vr
=
S
2
S
1
.
The standard of judgment is represented as follows [24]:
(42)
C
vr
=
{
<

0.35
;
Excellence
<

0.60
;
Pass
<

0.65
;
Reluctant pass
≥

0.65
;
No pass.


### 3.3. Volatility of Costs

Capital development in an underground mine consists of shafts, ramps, raises, and lateral transport drifts required to access ore deposits with expected utility greater than one year. This is the development required to start up the ore production and to haul the ore to the surface. Experience with investments in capital development might show that such expenditures can run considerably higher than the estimates, but it is quite unlikely that actual costs will be lower than estimated. Thus, the interval number might represent capital investment for the project:I
=
[
I
l
,
I
u
].Operating costs 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. These costs are generally more difficult to estimate than capital costs for most mining ventures. If we take into consideration that production will be carried out for many years, then it is very important to predict the future states of operating costs. Although there is some intention to create a correlation between metal price and operating cost, it is very hard to define it, since price and cost vary continuously and are different over time. At the project level, there will not be a perfect correlation between price and cost because of adjustments to variables such as labour, energy, explosives, and fuel, as well as other material expenditures that are supplied by industries that are not directly linked to metal price fluctuations. In order to protect themselves, suppliers are offering short-term contracts to mines that are the opposite of traditional long-term contracts. Some components of the operating cost such as inputs used for mineral processing are usually purchased at market prices that fluctuate monthly, annually, or even over shorter periods.The uncertainties related to the future states of operating costs are modelled with a special stochastic process, 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 and Sambandham [32]).In this model, we apply a continuous time process using the Itô-Doob type stochastic differential equation to describe the movement of operating costs. A general stochastic differential equation takes the following form:(43)
d
C
t
=
μ
(
C
t
,
t
)
d
t
+
σ
(
C
t
,
t
)
d
W
t
,
C
t
0
=
C
0
.
Here, t
≥
t
0, W
t is a Brownian motion, and C
t
>
0; this is the cost process. C
t is called the geometric Brownian motion (GBM), which is the solution of the following linear Itô-Doob type stochastic differential equation:
(44)
d
C
t
=
μ
C
t
d
t
+
σ
C
t
d
W
t
,
where μ and σ are some constants (μ is called the drift and σ is called the volatility) and W
t is a normalized Brownian motion. Using the Itô-Doob formula applied to f
(
C
t
)
=
ln
⁡
C
t, we can solve (44). The solution of (44) is given by the exact discrete-time equation for C
t:
(45)
C
t
=
C
t
-
1
·
e
{
(
μ
-
σ
2
/
2
)
Δ
t
+
N
(
0,1
)
σ
Δ
t
}
,
where N
(
0,1
) is the normally distributed random variable and Δ
t
=
1 (year). Equation (45) describes an operating cost scenario with spot costs C
t. By simulating C
t, we obtain operating costs for every year. Simulated values of the costs are obtained by performing the following calculations:
(46)
C
=
C
t
s
=
[
C
1
1
C
2
1
…
C
T
1
C
1
2
C
2
2
…
C
T
2
⋮
⋮
⋮
⋮
C
1
S
C
2
S
…
C
T
S
]
,
HHHHH
s
=
1,2
,
…
,
S
;
t
=
1,2
,
…
,
T
,
where S denotes the number of simulations and T the number of project years. In the space C
t
s each row represents one simulated path of costs over the project time, while each column represents simulated values of costs for every year. The main objective of using simulation is to determine the distribution of the C for every year of the project. In this way we obtain the sequence of probability density functions of C; C
t
~
(
pd
f
t
,
µ
t
,
σ
t
)
,
t
=
1,2
,
…
,
T.Applying the same concept of metal prices transformation, we obtain the future sequence of operating costs expressed by interval numbers;C
(
t
i
)
=
[
C
l
(
t
i
)
,
C
u
(
t
i
)
]
,
i
=
1,2
,
…
,
T, where T is the total project time.


### 3.4. Criterion of the Evaluation

The net present value (NPV) of the mine project is an integral evaluation criterion that recognizes the time effect of money over the life-of-mine. It is calculated as a difference between the sum of discounted values of estimated future cash flows and the initial investment and can be defined as follows:(47)
⊗
NPV
=
∑
t
=
1
T
(
(
Q
·
{
⊗
P
(
t
p
+
t
i
)
·
(
m
con
-
m
mr
)
·
(
G
·
M
/
m
con
)
-
⊗
C
(
t
p
+
t
i
)
}
)
×
(
(
1
+
⊗
r
)
t
i
)
-
1
)
-
⊗
I
,
where Q is annual ore production (t/year), r is discount rate, T is the number of periods for the life of the investment, I is initial (capital) investment, and t
p is preproduction time, time needed to prepare deposit to be mined (construction time)Finally, the last parameter that can be expressed by an interval number is the discount rate. Discounted cash flow methods are widely used in capital budgeting; however, determining the discount rate as a crisp value can lead to erroneous results in most mine project applications. A discount rate range can be established in a way which is either just acceptable (maximum value) or reasonable (minimum value);r
=
[
r
l
,
r
u
].A positive NPV will lead to the acceptance of the project and a negative NPV rejects it; that is,⊗
~
NPV
≥
0. Consider
(48)
⊗
~
(
[
NPV
l
,
NPV
u
]
)
=
ω
·
NPV
l
+
(
1
-
ω
)
·
NPV
u
,
HHHHHHHHHHHHHHHHHHHHH
ω
∈
[
0,1
]
.
Weight whitenization of interval NPV is obtained by ω
NPV
=
0.5.


## 3.1. The Concept of Evaluation

Economic evaluation of a mine project requires estimation of the revenues and costs throughout the defined lifetime of the mine. Such evaluation can be treated as strategic decision making under multiple sources of uncertainties. Therefore, to make the best decision, based on the information available, it is necessary to develop an adequate model incorporating the uncertainty of the input parameters. The model should be able to involve a common time horizon, taking the characteristics of the input variables that directly affect the value of the proposed project.The model is developed on the basis of full discounted cash flow analysis of an underground zinc mine project. The operating discounted cash flows are usually estimated on an annual basis. Net present value of investment is used as a key criterion in the process of mine project estimation. The expected net present value of the project is a function of the variables as(25)
E
(
NPV
∣
Q
,
P
,
G
,
M
,
C
,
I
,
r
,
t
)
≥
0
,
where Q denotes the production rate (capacity); P denotes the zinc metal price; G is the grade; M is mill recovery; C is operating costs; I is capital investment; r is discount rate, and t is the lifetime of the project; that is, the period in which the cash flow is generated.In this paper, we treat in detail only the variability of metal prices and operating costs, without intending to decrease the significance of the remaining parameters. These parameters are taken into account on the basis of expert knowledge (estimation).


## 3.2. Forecasting the Revenue of the Mine

Most mining companies realize their revenues by selling metal concentrates as a final product. Estimating mineral project revenue 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 value of the metal concentrate can be expressed as follows:(26)
V
con
=
P
·
(
m
con
-
m
mr
)
,
where P is metal price ($
/
t), m
con is metal content of concentrate (%), m
mr is metal recovery ratio (%), and
(27)
m
con
-
m
mr
=
{
(
m
con
%
-
8
)
·
100
m
con
%
≤
85
%
;
m
con
%
-
8
(
m
con
%
-
8
)
·
100
m
con
%
>
85
%
;
85
%
.
Annual mine revenue is calculated according to the following equation:
(28)
R
year
=
Q
year
·
V
con
·
G
·
M
m
con
,
where Q
year is annual ore production (t
/
year), Gis grade of the ore mined (%), M is mill recovery (%).Annual ore production is derived from the mining project schedule and is defined as crisp value. The concept of grade(
G
) is defined as the ratio of useful mass of metal to the total mass of ore and its critical value fluctuates over deposit space and can be estimated by experts and defined as interval number G
=
[
G
l
,
G
u
]. Mill recovery (
M
) is related to the flotation as the most widely used method for the concentration of fine grained minerals. It can also be defined as interval number M
=
[
M
l
,
M
u
]. Metal content represents the quality of concentrate and we also apply the concept of an interval number to define it:
m
con
=
[
m
con
l
,
m
con
u
].The major external source of risk affecting mine revenue is related to the uncertainty about market behaviour of metal prices. Forecasting the precise future state of the metal price is a very difficult task for mine planners. To predict future metal prices, we apply the concept based on the transformation of historical metal prices into adequate fuzzy-interval numbers and grey system theory.The forecasting model of metal prices is composed of the following steps.Step 1.
Create the setPDF
=
{
pdf
j
}
,
j
=
1,2
,
…
,
N, where pdf
j is the probability density function of metal prices for every historical year. The minimum number of elements of the set is four: j
min
⁡
=
1,2
,
3,4.Step 2.
Transform the set PDF into the setTFN
=
{
(
a
j
,
b
j
,
c
j
)
}, where (
a
j
,
b
j
,
c
j
) is an adequate fuzzy triangular number.Step 3.
Transform the set TFN into the setINT
=
{
[
a
j
l
,
a
j
u
]
}, where [
a
j
l
,
a
j
u
] is an adequate interval number.Step 4.
Using grey system prediction theory, create a grey differential equation of typeG
M
(
1,1
), that is, the first-order variable grey derivative.Step 5.
Testing ofG
M
(
1,1
) by residual error testing and the posterior error detection method.


### 3.2.1. Analysis of Historical Metal Prices

For every historical year it is necessary to define a probability density function with the following characteristics: shape by histogram, mean value(
µ
j
), and standard deviation (
σ
j
). In this way, we obtain the sequence of probability density functions of P; P
j
~
(
pdf
j
,
µ
j
,
σ
j
)
,
j
=
1,2
,
…
,
N, where N is the total number of historical years.


### 3.2.2. Fuzzification of Metal Prices

The sequence of obtainedpdf
j of P
j can be transformed into a sequence of triangular fuzzy numbers of P
j; P
j
~
TFN
j, j
=
1,2
,
…
,
N; that is, P
1
~
pdf
1
→
P
1
~
TFN
1; P
2
~
pdf
2
→
P
2
~
TFN
2
;
…; P
N
~
pdf
N
→
P
N
~
TFN
N. The method of transformation is based on the following facts: the support of the membership function and the pdf are the same, and the point with the highest probability (likelihood) has the highest possibility. For more details, see Swishchuk et al. [21]. The uncertainty in the P parameter is modelled by a triangular fuzzy number with the membership function which has the support of µ
j
-
2
σ
j
<
P
j
<
µ
j
+
2
σ
j
,
j
=
1,2
,
…
,
N, set up for around 95% confidence interval of distribution function. If we take into consideration that the triangular fuzzy number is defined as a triplet (
a
1
j
,
a
2
j
,
a
3
j
), then a
1
j and a
3
j are the lower bound and upper bound obtained from the lower and upper bound of 5% of the distribution, and the most promising value a
2
j is equal to the mean value of the distribution. For more details, see Do et al. [22].


### 3.2.3. Metal Prices as Interval Numbers

The sequence of obtainedTFN
j of P
j is transformed into a sequence of interval numbers of P
j
;
P
j
~
IN
T
j
,
j
=
1,2
,
…
,
N; that is, P
1
~
TF
N
1
→
P
1
~
IN
T
1
;
P
2
~
TF
N
2
→
P
2
~
IN
T
2
;
…;P
N
~
TF
N
N
→
P
N
~
IN
T
N. The method of transformation is described in Section 2.1 (basic concepts of fuzzy set theory). According to this method of transformation, we obtain set
(29)
P
j
~
IN
T
j
;
P
j
=
{
[
P
j
l
,
P
j
u
]
}
=
{
[
a
2
j
-
α
j
3
,
a
2
j
+
β
j
3
]
}
,
j
=
1,2
,
…
,
N
.


### 3.2.4. Forecasting Model of Metal Prices

The grey model is a powerful tool for forecasting the behaviour of the system in the future. It has been successfully applied to various fields since it was proposed by Deng [23–31]. In this paper we use a one-variable first-order differential grey equation, G
M
(
1,1
). The essence of G
M
(
1,1
) is to accumulate the original data (historical metal prices) in order to obtain regular data. By setting up the grey differential equation, we obtain the fitted curve in order to predict the future states of the system.Definition 9.
Assume that(30)
P
(
0
)
(
t
j
)
=
{
[
P
(
0
)
l
(
t
j
)
,
P
(
0
)
u
(
t
j
)
]
}
=
{
[
P
(
0
)
l
(
t
1
)
,
P
(
0
)
u
(
t
1
)
]
,
[
P
(
0
)
l
(
t
2
)
,
P
(
0
)
u
(
t
2
)
]
,
…
,
[
P
(
0
)
l
(
t
N
)
,
P
(
0
)
u
(
t
N
)
]
}
is the original series of interval metal prices obtained by transformation. Sampling interval isΔ
t
j
=
t
j
-
t
j
-
1
=
1
year.Definition 10.
LetP
(
1
)
(
t
j
)
=
{
[
P
(
1
)
l
(
t
j
)
,
P
(
1
)
u
(
t
j
)
]
} be a new sequence generated by the accumulated generating operation (AGO), where
(31)
P
(
1
)
l
(
t
j
)
=
∑
j
=
1
N
P
(
0
)
l
(
t
j
)
,
P
(
1
)
u
(
t
j
)
=
∑
j
=
1
N
P
(
0
)
u
(
t
j
)
.
In the process of forecasting metal prices, P
(
1
)
l
(
t
j
) and P
(
1
)
u
(
t
j
) are the solutions of the following grey differential equation:
(32)
d
[
P
(
1
)
l
(
t
)
,
P
(
1
)
u
(
t
)
]
d
t
+
[
q
l
,
q
u
]
·
[
P
(
1
)
l
(
t
)
,
P
(
1
)
u
(
t
)
]
=
[
w
l
,
w
u
]
.
Obviously, (32) is an interval-valued differential equation (see (20)).
To get the values of parametersq
=
[
q
l
,
q
u
] and w
=
[
w
l
,
w
u
], the least square method is used as follows:
(33)
⊗
[
q
w
]
=
⊗
(
[
B
T
B
]
-
1
B
T
Y
(
0
)
)
,
where
(34)
⊗
B
=
[
-
1
2
(
P
(
0
)
(
1
)
+
P
(
0
)
(
2
)
)
1
-
1
2
(
P
(
0
)
(
2
)
+
P
(
0
)
(
3
)
)
⋮
1
⋮
-
1
2
(
P
(
0
)
(
N
-
1
)
+
P
(
0
)
(
N
)
)
1
]
⊗
Y
(
0
)
=
[
P
(
0
)
(
2
)
P
(
0
)
(
3
)
⋯
P
(
0
)
(
N
)





]
T
.
Note that the sign ⊗ indicates the interval number.
IfP
˙
(
1
)
(
t
) is considered as (i)-Hukuhara differentiable, then the following system of ordinary differential equations is as follows:
(35)
P
˙
(
1
)
l
(
t
)
=
w
l
-
q
u
·
P
(
1
)
u
(
t
)
,
P
l
(
t
0
)
=
P
(
0
)
l
(
1
)
,
P
˙
(
1
)
u
(
t
)
=
w
u
-
q
l
·
P
(
1
)
l
(
t
)
,
P
u
(
t
0
)
=
P
(
0
)
u
(
1
)
.
The solution of this system (forecasted equation) is as follows:
(36)
P
^
(
1
)
l
(
t
j
+
1
)
=
1
2
·
e
q
l
·
q
u
·
t
·
Z
1
+
1
2
·
e
-
q
l
·
q
u
·
t
·
Z
2
+
w
u
q
l
,
j
=
1,2
,
…
,
N
,
P
^
(
1
)
u
(
t
j
+
1
)
=
1
2
·
e
q
l
·
q
u
·
t
·
Z
3
+
1
2
·
e
-
q
l
·
q
u
·
t
·
Z
4
+
w
l
q
u
,
j
=
1,2
,
…
,
N
,
where
(37)
Z
1
=
P
(
0
)
l
(
1
)
-
q
u
·
P
(
0
)
u
(
1
)
q
l
·
q
u
+
w
l
q
l
·
q
u
-
w
u
q
l
,
Z
2
=
P
(
0
)
l
(
1
)
+
q
u
·
P
(
0
)
u
(
1
)
q
l
·
q
u
-
w
l
q
l
·
q
u
-
w
u
q
l
,
Z
3
=
P
(
0
)
u
(
1
)
-
P
(
0
)
l
(
1
)
·
q
l
·
q
u
q
u
+
w
u
q
l
·
q
u
-
w
l
q
u
,
Z
4
=
P
(
0
)
u
(
1
)
+
P
(
0
)
l
(
1
)
·
q
l
·
q
u
q
u
-
w
u
q
l
·
q
u
-
w
l
q
u
.
To obtain the forecasted value of the primitive (original) metal price data at time (
t
j
+
1
), the inverse accumulated generating operation (IAGO) is used as follows:
(38)
⊗
P
^
(
0
)
(
t
1
)
=
⊗
P
(
0
)
(
t
1
)
⊗
P
^
(
0
)
(
t
j
+
1
)
=
⊗
P
^
(
1
)
(
t
j
+
1
)
-
⊗
P
^
(
1
)
(
t
j
)
,
j
=
1,2
,
…
,
N
.


### 3.2.5. The Model Accuracy

Definition 11.
For a given interval grey number⊗
(
[
a
l
,
a
u
]
), it is common to take a whitening value ⊗
~
(
[
a
l
,
a
u
]
)
=
ω
·
a
l
+
(
1
-
ω
)
·
a
u for ∀
ω
∈
[
0,1
]. Furthermore, if ω
=
0.5, it is called equal weight average whitenization [23].
Residual error testing is composed of the calculation of relative error and absolute error between⊗
~
P
(
0
)
(
t
j
) and ⊗
~
P
^
(
0
)
(
t
j
) based on the following formulas:
(39)
Δ
ε
(
t
j
)
=
⊗
~
P
(
0
)
(
t
j
)
-
⊗
~
P
^
(
0
)
(
t
j
)
,
|
Δ
ε
(
t
j
)
|
=
|
⊗
~
P
(
0
)
(
t
j
)
-
⊗
~
P
^
(
0
)
(
t
j
)
|
.
The posterior error detection method means calculation of the standard deviation of original metal price series (S
1) and standard deviation of absolute error (S
2):
(40)
S
1
=
∑
j
=
1
N
(
⊗
~
P
(
0
)
(
t
j
)
-
⊗
~
P
-
(
0
)
(
t
j
)
)
2
N
-
1
,
S
2
=
∑
j
=
1
N
(
Δ
ε
(
t
j
)
-
Δ
ε
-
(
t
j
)
)
2
N
-
1
.
The variance ratio is equal to
(41)
C
vr
=
S
2
S
1
.
The standard of judgment is represented as follows [24]:
(42)
C
vr
=
{
<

0.35
;
Excellence
<

0.60
;
Pass
<

0.65
;
Reluctant pass
≥

0.65
;
No pass.


## 3.2.1. Analysis of Historical Metal Prices

For every historical year it is necessary to define a probability density function with the following characteristics: shape by histogram, mean value(
µ
j
), and standard deviation (
σ
j
). In this way, we obtain the sequence of probability density functions of P; P
j
~
(
pdf
j
,
µ
j
,
σ
j
)
,
j
=
1,2
,
…
,
N, where N is the total number of historical years.


## 3.2.2. Fuzzification of Metal Prices

The sequence of obtainedpdf
j of P
j can be transformed into a sequence of triangular fuzzy numbers of P
j; P
j
~
TFN
j, j
=
1,2
,
…
,
N; that is, P
1
~
pdf
1
→
P
1
~
TFN
1; P
2
~
pdf
2
→
P
2
~
TFN
2
;
…; P
N
~
pdf
N
→
P
N
~
TFN
N. The method of transformation is based on the following facts: the support of the membership function and the pdf are the same, and the point with the highest probability (likelihood) has the highest possibility. For more details, see Swishchuk et al. [21]. The uncertainty in the P parameter is modelled by a triangular fuzzy number with the membership function which has the support of µ
j
-
2
σ
j
<
P
j
<
µ
j
+
2
σ
j
,
j
=
1,2
,
…
,
N, set up for around 95% confidence interval of distribution function. If we take into consideration that the triangular fuzzy number is defined as a triplet (
a
1
j
,
a
2
j
,
a
3
j
), then a
1
j and a
3
j are the lower bound and upper bound obtained from the lower and upper bound of 5% of the distribution, and the most promising value a
2
j is equal to the mean value of the distribution. For more details, see Do et al. [22].


## 3.2.3. Metal Prices as Interval Numbers

The sequence of obtainedTFN
j of P
j is transformed into a sequence of interval numbers of P
j
;
P
j
~
IN
T
j
,
j
=
1,2
,
…
,
N; that is, P
1
~
TF
N
1
→
P
1
~
IN
T
1
;
P
2
~
TF
N
2
→
P
2
~
IN
T
2
;
…;P
N
~
TF
N
N
→
P
N
~
IN
T
N. The method of transformation is described in Section 2.1 (basic concepts of fuzzy set theory). According to this method of transformation, we obtain set
(29)
P
j
~
IN
T
j
;
P
j
=
{
[
P
j
l
,
P
j
u
]
}
=
{
[
a
2
j
-
α
j
3
,
a
2
j
+
β
j
3
]
}
,
j
=
1,2
,
…
,
N
.


## 3.2.4. Forecasting Model of Metal Prices

The grey model is a powerful tool for forecasting the behaviour of the system in the future. It has been successfully applied to various fields since it was proposed by Deng [23–31]. In this paper we use a one-variable first-order differential grey equation, G
M
(
1,1
). The essence of G
M
(
1,1
) is to accumulate the original data (historical metal prices) in order to obtain regular data. By setting up the grey differential equation, we obtain the fitted curve in order to predict the future states of the system.Definition 9.
Assume that(30)
P
(
0
)
(
t
j
)
=
{
[
P
(
0
)
l
(
t
j
)
,
P
(
0
)
u
(
t
j
)
]
}
=
{
[
P
(
0
)
l
(
t
1
)
,
P
(
0
)
u
(
t
1
)
]
,
[
P
(
0
)
l
(
t
2
)
,
P
(
0
)
u
(
t
2
)
]
,
…
,
[
P
(
0
)
l
(
t
N
)
,
P
(
0
)
u
(
t
N
)
]
}
is the original series of interval metal prices obtained by transformation. Sampling interval isΔ
t
j
=
t
j
-
t
j
-
1
=
1
year.Definition 10.
LetP
(
1
)
(
t
j
)
=
{
[
P
(
1
)
l
(
t
j
)
,
P
(
1
)
u
(
t
j
)
]
} be a new sequence generated by the accumulated generating operation (AGO), where
(31)
P
(
1
)
l
(
t
j
)
=
∑
j
=
1
N
P
(
0
)
l
(
t
j
)
,
P
(
1
)
u
(
t
j
)
=
∑
j
=
1
N
P
(
0
)
u
(
t
j
)
.
In the process of forecasting metal prices, P
(
1
)
l
(
t
j
) and P
(
1
)
u
(
t
j
) are the solutions of the following grey differential equation:
(32)
d
[
P
(
1
)
l
(
t
)
,
P
(
1
)
u
(
t
)
]
d
t
+
[
q
l
,
q
u
]
·
[
P
(
1
)
l
(
t
)
,
P
(
1
)
u
(
t
)
]
=
[
w
l
,
w
u
]
.
Obviously, (32) is an interval-valued differential equation (see (20)).
To get the values of parametersq
=
[
q
l
,
q
u
] and w
=
[
w
l
,
w
u
], the least square method is used as follows:
(33)
⊗
[
q
w
]
=
⊗
(
[
B
T
B
]
-
1
B
T
Y
(
0
)
)
,
where
(34)
⊗
B
=
[
-
1
2
(
P
(
0
)
(
1
)
+
P
(
0
)
(
2
)
)
1
-
1
2
(
P
(
0
)
(
2
)
+
P
(
0
)
(
3
)
)
⋮
1
⋮
-
1
2
(
P
(
0
)
(
N
-
1
)
+
P
(
0
)
(
N
)
)
1
]
⊗
Y
(
0
)
=
[
P
(
0
)
(
2
)
P
(
0
)
(
3
)
⋯
P
(
0
)
(
N
)





]
T
.
Note that the sign ⊗ indicates the interval number.
IfP
˙
(
1
)
(
t
) is considered as (i)-Hukuhara differentiable, then the following system of ordinary differential equations is as follows:
(35)
P
˙
(
1
)
l
(
t
)
=
w
l
-
q
u
·
P
(
1
)
u
(
t
)
,
P
l
(
t
0
)
=
P
(
0
)
l
(
1
)
,
P
˙
(
1
)
u
(
t
)
=
w
u
-
q
l
·
P
(
1
)
l
(
t
)
,
P
u
(
t
0
)
=
P
(
0
)
u
(
1
)
.
The solution of this system (forecasted equation) is as follows:
(36)
P
^
(
1
)
l
(
t
j
+
1
)
=
1
2
·
e
q
l
·
q
u
·
t
·
Z
1
+
1
2
·
e
-
q
l
·
q
u
·
t
·
Z
2
+
w
u
q
l
,
j
=
1,2
,
…
,
N
,
P
^
(
1
)
u
(
t
j
+
1
)
=
1
2
·
e
q
l
·
q
u
·
t
·
Z
3
+
1
2
·
e
-
q
l
·
q
u
·
t
·
Z
4
+
w
l
q
u
,
j
=
1,2
,
…
,
N
,
where
(37)
Z
1
=
P
(
0
)
l
(
1
)
-
q
u
·
P
(
0
)
u
(
1
)
q
l
·
q
u
+
w
l
q
l
·
q
u
-
w
u
q
l
,
Z
2
=
P
(
0
)
l
(
1
)
+
q
u
·
P
(
0
)
u
(
1
)
q
l
·
q
u
-
w
l
q
l
·
q
u
-
w
u
q
l
,
Z
3
=
P
(
0
)
u
(
1
)
-
P
(
0
)
l
(
1
)
·
q
l
·
q
u
q
u
+
w
u
q
l
·
q
u
-
w
l
q
u
,
Z
4
=
P
(
0
)
u
(
1
)
+
P
(
0
)
l
(
1
)
·
q
l
·
q
u
q
u
-
w
u
q
l
·
q
u
-
w
l
q
u
.
To obtain the forecasted value of the primitive (original) metal price data at time (
t
j
+
1
), the inverse accumulated generating operation (IAGO) is used as follows:
(38)
⊗
P
^
(
0
)
(
t
1
)
=
⊗
P
(
0
)
(
t
1
)
⊗
P
^
(
0
)
(
t
j
+
1
)
=
⊗
P
^
(
1
)
(
t
j
+
1
)
-
⊗
P
^
(
1
)
(
t
j
)
,
j
=
1,2
,
…
,
N
.


## 3.2.5. The Model Accuracy

Definition 11.
For a given interval grey number⊗
(
[
a
l
,
a
u
]
), it is common to take a whitening value ⊗
~
(
[
a
l
,
a
u
]
)
=
ω
·
a
l
+
(
1
-
ω
)
·
a
u for ∀
ω
∈
[
0,1
]. Furthermore, if ω
=
0.5, it is called equal weight average whitenization [23].
Residual error testing is composed of the calculation of relative error and absolute error between⊗
~
P
(
0
)
(
t
j
) and ⊗
~
P
^
(
0
)
(
t
j
) based on the following formulas:
(39)
Δ
ε
(
t
j
)
=
⊗
~
P
(
0
)
(
t
j
)
-
⊗
~
P
^
(
0
)
(
t
j
)
,
|
Δ
ε
(
t
j
)
|
=
|
⊗
~
P
(
0
)
(
t
j
)
-
⊗
~
P
^
(
0
)
(
t
j
)
|
.
The posterior error detection method means calculation of the standard deviation of original metal price series (S
1) and standard deviation of absolute error (S
2):
(40)
S
1
=
∑
j
=
1
N
(
⊗
~
P
(
0
)
(
t
j
)
-
⊗
~
P
-
(
0
)
(
t
j
)
)
2
N
-
1
,
S
2
=
∑
j
=
1
N
(
Δ
ε
(
t
j
)
-
Δ
ε
-
(
t
j
)
)
2
N
-
1
.
The variance ratio is equal to
(41)
C
vr
=
S
2
S
1
.
The standard of judgment is represented as follows [24]:
(42)
C
vr
=
{
<

0.35
;
Excellence
<

0.60
;
Pass
<

0.65
;
Reluctant pass
≥

0.65
;
No pass.


## 3.3. Volatility of Costs

Capital development in an underground mine consists of shafts, ramps, raises, and lateral transport drifts required to access ore deposits with expected utility greater than one year. This is the development required to start up the ore production and to haul the ore to the surface. Experience with investments in capital development might show that such expenditures can run considerably higher than the estimates, but it is quite unlikely that actual costs will be lower than estimated. Thus, the interval number might represent capital investment for the project:I
=
[
I
l
,
I
u
].Operating costs 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. These costs are generally more difficult to estimate than capital costs for most mining ventures. If we take into consideration that production will be carried out for many years, then it is very important to predict the future states of operating costs. Although there is some intention to create a correlation between metal price and operating cost, it is very hard to define it, since price and cost vary continuously and are different over time. At the project level, there will not be a perfect correlation between price and cost because of adjustments to variables such as labour, energy, explosives, and fuel, as well as other material expenditures that are supplied by industries that are not directly linked to metal price fluctuations. In order to protect themselves, suppliers are offering short-term contracts to mines that are the opposite of traditional long-term contracts. Some components of the operating cost such as inputs used for mineral processing are usually purchased at market prices that fluctuate monthly, annually, or even over shorter periods.The uncertainties related to the future states of operating costs are modelled with a special stochastic process, 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 and Sambandham [32]).In this model, we apply a continuous time process using the Itô-Doob type stochastic differential equation to describe the movement of operating costs. A general stochastic differential equation takes the following form:(43)
d
C
t
=
μ
(
C
t
,
t
)
d
t
+
σ
(
C
t
,
t
)
d
W
t
,
C
t
0
=
C
0
.
Here, t
≥
t
0, W
t is a Brownian motion, and C
t
>
0; this is the cost process. C
t is called the geometric Brownian motion (GBM), which is the solution of the following linear Itô-Doob type stochastic differential equation:
(44)
d
C
t
=
μ
C
t
d
t
+
σ
C
t
d
W
t
,
where μ and σ are some constants (μ is called the drift and σ is called the volatility) and W
t is a normalized Brownian motion. Using the Itô-Doob formula applied to f
(
C
t
)
=
ln
⁡
C
t, we can solve (44). The solution of (44) is given by the exact discrete-time equation for C
t:
(45)
C
t
=
C
t
-
1
·
e
{
(
μ
-
σ
2
/
2
)
Δ
t
+
N
(
0,1
)
σ
Δ
t
}
,
where N
(
0,1
) is the normally distributed random variable and Δ
t
=
1 (year). Equation (45) describes an operating cost scenario with spot costs C
t. By simulating C
t, we obtain operating costs for every year. Simulated values of the costs are obtained by performing the following calculations:
(46)
C
=
C
t
s
=
[
C
1
1
C
2
1
…
C
T
1
C
1
2
C
2
2
…
C
T
2
⋮
⋮
⋮
⋮
C
1
S
C
2
S
…
C
T
S
]
,
HHHHH
s
=
1,2
,
…
,
S
;
t
=
1,2
,
…
,
T
,
where S denotes the number of simulations and T the number of project years. In the space C
t
s each row represents one simulated path of costs over the project time, while each column represents simulated values of costs for every year. The main objective of using simulation is to determine the distribution of the C for every year of the project. In this way we obtain the sequence of probability density functions of C; C
t
~
(
pd
f
t
,
µ
t
,
σ
t
)
,
t
=
1,2
,
…
,
T.Applying the same concept of metal prices transformation, we obtain the future sequence of operating costs expressed by interval numbers;C
(
t
i
)
=
[
C
l
(
t
i
)
,
C
u
(
t
i
)
]
,
i
=
1,2
,
…
,
T, where T is the total project time.


## 3.4. Criterion of the Evaluation

The net present value (NPV) of the mine project is an integral evaluation criterion that recognizes the time effect of money over the life-of-mine. It is calculated as a difference between the sum of discounted values of estimated future cash flows and the initial investment and can be defined as follows:(47)
⊗
NPV
=
∑
t
=
1
T
(
(
Q
·
{
⊗
P
(
t
p
+
t
i
)
·
(
m
con
-
m
mr
)
·
(
G
·
M
/
m
con
)
-
⊗
C
(
t
p
+
t
i
)
}
)
×
(
(
1
+
⊗
r
)
t
i
)
-
1
)
-
⊗
I
,
where Q is annual ore production (t/year), r is discount rate, T is the number of periods for the life of the investment, I is initial (capital) investment, and t
p is preproduction time, time needed to prepare deposit to be mined (construction time)Finally, the last parameter that can be expressed by an interval number is the discount rate. Discounted cash flow methods are widely used in capital budgeting; however, determining the discount rate as a crisp value can lead to erroneous results in most mine project applications. A discount rate range can be established in a way which is either just acceptable (maximum value) or reasonable (minimum value);r
=
[
r
l
,
r
u
].A positive NPV will lead to the acceptance of the project and a negative NPV rejects it; that is,⊗
~
NPV
≥
0. Consider
(48)
⊗
~
(
[
NPV
l
,
NPV
u
]
)
=
ω
·
NPV
l
+
(
1
-
ω
)
·
NPV
u
,
HHHHHHHHHHHHHHHHHHHHH
ω
∈
[
0,1
]
.
Weight whitenization of interval NPV is obtained by ω
NPV
=
0.5.


## 4. Numerical Example

The management of a small mining company is evaluating the opening of a new zinc deposit. The recommendations from the prefeasibility study suggest the following:(i)
the underground mine development system connecting the ore body to the surface is based on the combination of ramp and horizontal drives. This system is used for the purpose of ore haulage by dump trucks and conduct intake fresh air. Contaminated air is conducted to the surface by horizontal drives and declines,
(ii)
they suggest purchasing new mining equipment.The completion of this project will cost the company about 3500 000 USD over three years. At the beginning of the fourth year, when construction is completed, the new mine will produce 100,000t/year during 5 years of production.The input parameters required for the project evaluation are given in Table1. Note that the situation is hypothetical and the numbers used are to permit calculation.Table 1
Input parameters.
Mine ore production (t/year)
100 000
Zinc grade (%)
[4.2, 5.0] = [0.042, 0.050]
Metal content of concentrate (%)
m
con
-
m
mr
=
{
(
m
con
%
-
8
)
·
100
m
con
%
≤
85
%
;
m
con
%
-
8
(
m
con
%
-
8
·
100
)
m
con
%
>
85
%
;
85
%
[47, 52] = [0.47, 0.52]
[39, 42] = [0.39, 0.42]
Mill recovery (%)
[75, 80] = [0.75, 0.80]
Metal price ($/t)
2009
2010
2011
2012
2013
January
1203
2415
2376
1989
2031
February
1118
2159
2473
2058
2129
March
1223
2277
2341
2036
1929
April
1388
2368
2371
2003
1856
May
1492
1970
2160
1928
1831
June
1555
1747
2234
1856
1839
July
1583
1847
2398
1848
1838
August
1818
2047
2199
1816
1896
September
1879
2151
2075
2010
1847
October
2071
2374
1871
1904
1885
November
2197
2283
1935
1912
1866
December
2374
2287
1911
2040
1975
Initial (capital) investment ($)
[3000000, 4000000]
Operating costs, geometric Brownian motion,yearly time resolution ($/t), equation (39)Number of simulations
Spot value 32; drift 0.020;cost volatility rate 0.10;
S
=
500
Construction period (year)
32014
;
2015
;
2016
Mine life (year)
52017
;
2018
;
…
;
2021
Discount rate (%)
[7, 10] = [0.07, 0.10]The interval values of historical metal prices are calculated by Steps1, 2, and 3 of the metal prices forecasting model. The values obtained are as represented in Table 2 and Figure 2.Table 2
Transformationpdf
j
→
TFN
j
→
INT
j.
2009
2010
2011
2012
2013
pdf
j
μ
j
1658
2160
2195
1950
1910
σ
j
410
216
207
83
92
TFN
j
a
1
j
=
μ
j
-
2
σ
j
838
1728
1781
1784
1726
a
2
j
=
μ
j
1658
2160
2195
1950
1910
a
3
j
=
μ
j
+
2
σ
j
2478
2592
2609
2116
2094
α
j
=
β
j
820
432
414
166
184
INT
j
P
j
l
1385
2016
2057
1895
1849
P
j
u
1931
2304
2333
2005
1971Figure 2
Historical zinc metal prices expressed as triangular fuzzy numbers.Based on data in Table2 and (36), the fuzzy-interval AGOGM
(
1,1
) model of zinc metal prices is set up as follows:
(49)
P
^
(
1
)
l
(
t
j
+
1
)
=
-
14908.63
·
e
0.0454
·
t
-
59327.37
·
e
-
0.0454
·
t
+
75620.82
,
P
^
(
1
)
u
(
t
j
+
1
)
=
11055.73
·
e
0.0454
t
-
43995.13
·
e
-
0.0454
·
t
+
34871.36
.
To test the precision of the model, relative error and absolute error of the model are calculated and the results are represented in Table 3.Table 3
Relative and absolute error of the model.
Year
Original values ($/t)
Simulated values ($/t)
Whitenizationω = 0.5
Relative error
Absolute error
Relative error (%)
Original
Simulated
2009
[1385, 1931]
[1385, 1931]
1658
1658
0
0
0
2010
[2016, 2304]
[1941, 2466]
2160
2204
−44
44
−2.03
2011
[2057, 2333]
[1792, 2404]
2195
2098
97
97
+4.42
2012
[1895, 2005]
[1647, 2346]
1950
1997
−47
47
−2.41
2013
[1849, 1971]
[1505, 2293]
1910
1899
11
11
+0.57The standard deviation of the original metal price series (S
1) and standard deviation of absolute error (S
2) are 216.77 and 37.93, respectively. The variance ratio of the model is C
vr
=
0.175. These show that the obtained model has good forecasting precision to predict the zinc metal prices.According to (28), annual mine revenues are represented in Table 4.Table 4
Mine revenues over production period 2017–2021.
2017
2018
2019
2020
2021
R
year (mill’s USD)
[2.279, 7.466]
[1.974, 7.360]
[1.672, 7.269]
[1.374, 7.192]
[1.079, 7.131]Uncertainty related to the unit operating costs is quantified according to (45), that is, by geometric Brownian motion. Figure 3 represents seven possible paths (scenarios) of the unit operating costs over the project time.Figure 3
Seven simulated operating cost paths on a yearly time resolution.The results of the simulations and transformations are represented in Table5.Table 5
Simulation of the unit operating costs.
2013
2014
2015
2016
2017
2018
2019
2020
2021
Simulation 1
32
35.24
36.47
40.67
46.07
55.00
54.14
55.69
48.59
⋮
⋮
⋮
⋮
⋮
⋮
⋮
⋮
⋮
⋮
Simulation 500
32
31.17
35.17
37.91
37.66
36.27
37.15
34.34
36.08
pdf
i
μ
i
32
32.41
33.12
33.51
34.05
34.56
35.16
35.80
36.35
σ
i
0.0
3.16
4.55
5.70
6.60
7.49
8.45
9.42
10.00
TFN
i
a
1
i
=
μ
i
-
2
σ
i
32
26.09
24.02
22.11
20.85
19.58
18.26
16.96
16.35
a
2
i
=
μ
i
32
32.41
33.12
33.51
34.05
34.56
35.16
35.80
36.35
a
3
i
=
μ
i
+
2
σ
i
32
38.73
42.22
44.91
47.25
49.54
52.06
54.64
56.35
α
i
=
β
i
0.0
6.32
9.10
11.40
13.20
14.98
16.90
18.84
20.00
INT
i
C
i
l ($/t)
32
30.30
30.08
29.71
29.65
29.56
29.52
29.52
29.68
C
i
u ($/t)
32
34.52
36.16
37.32
38.45
39.56
42.09
42.08
43.01Annual operating costs are represented in Table6.Table 6
Operating costs over production period 2017–2021.
2017
2018
2019
2020
2021
C
year (mill’s USD)
[2.965, 3.845]
[2.956, 3.956]
[2.956, 4.029]
[2.952, 4.208]
[2.968, 4.301]The discounted cash flow of the project is represented in Table7.Table 7
Discounted cash flow over production period 2017–2021.
2017
2018
2019
2020
2021
R
year (mill’s USD)
[2.279, 7.466]
[1.974, 7.360]
[1.672, 7.269]
[1.374, 7.192]
[1.079, 7.131]
C
year (mill’s USD)
[2.965, 3.845]
[2.956, 3.956]
[2.956, 4.029]
[2.952, 4.208]
[2.968, 4.301]
R
year
-
C
year
[−1.566, 4.501]
[−1.982, 4.404]
[−2.357, 4.313]
[−2.834, 4.240]
[−3.222, 4.163]
Discount factor
[1.07, 1.10]
[1.14, 1.21]
[1.23, 1.33]
[1.31, 1.46]
[1.40, 1.61]
Present value
[−1.423, 4.206]
[−1.638, 3.863]
[−1.772, 3.506]
[−1.941, 3.236]
[−2.001, 2.973]According to (47) and data in Table 7, the net present value of the project is as follows:
(50)
⊗
NPV
=
[
-
8.775
,
17.784
]
-
[
3.000
,
4.000
]
=
[
-
11.775
,
13.784
]
mill’sUSD.
Weight whitenization of interval NPV is obtained by ω
NPV
=
0.5 and the white value is NPV = 1.004 million USD. This means the project is accepted.


## 5. Conclusion

The combined effect of market volatility and uncertainty about future commodity prices is posing higher risks to mining businesses across the globe. In such times, knowing how to unlock value by maximizing the value of resources and reserves through strategic mine planning is essential. In our country, small mining companies are faced with many problems but the primary problem is related to the shortage of capital for investment. In such an environment every mining venture must be treated as a strategic decision supported by adequate analysis. The developed economic model is a mathematical representation of project evaluation reality and allows management to see the impact of key parameters on the project value. The interaction between production, costs, and capital is highly complex and changes over time, but needs to be accurately modelled so as to provide insights around capital configurations of that business.The evaluation of a mining venture is made very difficult by uncertainty on the input variables in the project. Metal prices, costs, grades, discount rates, and countless other variables create a high risk environment to operate in. The incorporation of risk into modelling will provide management with better means to deal with uncertainty and the identification and quantification of those factors that most contribute to risk, which will then allow mitigation strategies to be tested. The model brings forth an issue that has the dynamic nature of the assessment of investment profitability. With the fuzzy-interval model, the future forecast can be done from the beginning of the process until the end.From the results obtained by numerical example, it is shown that fuzzy-interval grey system theory can be incorporated into mine project evaluation. The variance ratioC
vr
=
0.175 shows that the metal prices forecasting model is credible to predict the future values of the most important external parameter. The operating costs prediction model, based on geometric Brownian motion, gives the same result that we get if we use scenarios; however, it does not require us to simplify the future to the limited number of alternative scenarios.With interval numbers, the end result will be interval NPV, which is the payoff interval for the project. Using the weight whitenization of the interval NPV, we obtain the payoff crisp value for the project. This value is the value at risk, helping the management of the company to make the right decision.



---

*Source: 914643-2014-07-23.xml*