Abstract
Production planning in an underground mine plays a key activity in the mining company business. It is supported by the fact that mineral industry is unique and volatile environment. There are two uncertain parameters that cannot be managed by planners, metal price, and operating costs. Having ability to quantify and incorporate them in the process of planning can help companies to do their business in much easier way. We quantify these uncertainties by the simulation of mean reverting process and Itô-Doob stochastic differential equation, respectively. Mineral deposit is represented as a set of mineable blocks and room and pillar mining method is selected as a way of mining. Multicriteria clustering algorithm is used to create areas inside of mineral deposit that have technological characteristics required by the planners. We also developed a way to forecast the volatility of economic values of these areas through the planning period. Fuzzy 0-1 linear programming model is used to define the sequence of mining of these areas by maximization of the expected value of the fuzzy future cash flow. Model was tested on small hypothetical lead-zinc mineral deposit and results showed that our approach was able to solve such complex problem.
1. Introduction
Certainly, the investment environment associated with the mining industry is unique when compared with environment encountered by typical manufacturing industries. Some of characteristics which are often proclaimed as being unique are as follows. Virtually every knowledgeable observer would agree that mining ventures are extremely capital intensive. Even extremely small high grade precious metal mines employing only a handful of miners can rarely be developed for operation less than a million dollars. The amount of time required to develop a mining property for production can vary significantly. Once the occurrence of an ore deposit has been well established, it takes a number of years of intensive effort before the property is brought on stream and ore is produced on a continuous basis. In addition to the obvious risks associated with capital intensity and long lead times, there are a number of other risks associated with mining ventures. In general, these risks may be placed under the general headings of geological risks, engineering risks, economic risks, and political risks. Perhaps, the most unique aspect of the minerals industry is the fact that it deals with the extraction of a nonrenewable resource [1].
If we consider all these characteristics of the mining industry, we can conclude that production planning is essential activity which helps mining companies to do business in such environment. Planning is an optimization problem where the searching for the global optimum solution is very difficult and time consumption task. There are many approaches that try to solve this optimization problem.
O’Sullivan et al. explained very well why underground mine production planning is very difficult task and provided directions related to the application of heuristic techniques that rely on spatial and/or temporal aggregation to produce solutions [2]. Carlyle et al. used mixed-integer programming model which considers some planning constraints. Validation of this model was tested in one sector of the underground platinum mine [3].
Anani applied discrete event simulation to determine the optimal width of coal room and pillars panels under specific mining conditions. She also tested the hypothesis that heuristic preprocessing can be used to increase the computational efficiency of branch and cut solutions to the binary integer linear programming problem of room and pillar mine sequencing. The findings of her research include panel width optimization, a deterministic modelling framework that incorporates multiple mining risk in room and pillar production sequencing, and accounting for changing duty cycles in continuous miner-shuttle car matching [4]. Bakhtavar et al. used 0-1 integer programming to create model which optimizes the way of transition from open pit mine to the underground mine [5]. Nehring et al. developed a new mathematical programming model for optimization of production scheduling of a sublevel stopping operation, which significantly reduces solution times without altering results while maintaining all constraints. They represented all stope production phases by single binary variable and increased efficiency of mixed-integer programming in the process of optimization [6]. Bai et al. developed algorithm for stope design optimization at sublevel mining method. Optimization problem was treated as maximum flow over the adequate graph [7]. A general capacitated multicommodity network flow model has been used for long-term mine planning by Epstein et al. [8]. Grieco et al. applied a probabilistic mixed-integer programming approach to optimize stope in an underground mine, i.e., to define location, size, and number of active stopes with uncertainty related to ore grade and acceptable level of risk [9]. Nehring et al. integrated short and medium-term production plans by combining the short-term objective of minimizing deviation from targeted mill feed grade with the medium-term objective of maximizing net present value into a single mathematical optimization model [10]. Terblanche and Bley reduced the resolution of underground mine scheduling problem and applied mixed-integer programming to improve profitability through selective mining [11]. Kuchta et al. used mixed-integer programming to schedule Kiruna’s operations, specifically, which production blocks to mine and when to mine them to minimize deviations from monthly planned production quantities while adhering to operational restrictions [12]. Topal developed an early start and late start algorithm that defines the precedence restrictions for each mining unit in their mixed-integer linear programming model of the underground Kiruna Mine [13]. Hirschi developed a dynamic programming algorithm to supplant that trial and error practice of determining and evaluating room and pillar mining sequences. Dynamic programming has been used in mining to optimize multistage processes where production parameters are stage-specific [14]. Gligoric et al. developed a production planning model, which minimizes deviation from the acceptable rate of return, using multivariable weighted Frobenius distance function that measures the deviation from established targets. [15].
All developed mine production planning models were based on 0-1 linear programming and different methods have been used to find the extreme value of the linear objective function, for example, simplex method, simulated annealing, Branch and Bound algorithm, ant colony optimization, neural networks, etc. We applied fuzzy 0-1 linear programming to incorporate the uncertainties in the objective function and make the problem of production planning more realistic. By this way, we increase the precision of the obtained results. If we take into consideration that mine production planning belongs to the decision making field then fuzzy model really helps us to make final decision in more efficient way.
The main aim of this paper is to provide efficiency support to decision making on production planning in underground mines that use room and pillar mining method as a way of mining. Model is based on the maximization of fuzzy objective function which represents the present or discounted value of the future cash flow of production plan, with respect to the set of constraints. We consider this problem as zero-one linear programing problem in which only coefficients in the objective function are triangular fuzzy numbers. Coefficients represent the discounted economic value of the technological mining cut (TMC) which is a part of mineral deposit characterized with respect to the given set of technological requirements such as annual capacity of production, compactness of the shape of TMC, and standard deviation of ore grade in the TMC. Total number of technological mining cuts is equal to the total number of years of production.
The first step in the production planning model is related to the creation of TMCs having the value of attributes closely to the values of technological requirements. In the purpose of creation such TMCs (clusters) we developed fuzzy multicriteria clustering algorithm where uncertainties of some input data are quantified by triangular fuzzy numbers. Mining engineers uses a block model of the deposit that represents the deposit as three-dimensional array of blocks. Accordingly, clustering algorithm is applied on the set of these blocks. The second step concerns the calculation of discounted economic value of TMCs. It indicates that we are facing dynamic problem burdened with some uncertainties. These uncertainties come primary from the metal price and operating costs fluctuation through the time of planning. To estimate the future state of metal price, we developed forecasting algorithm which represents the hybrid of the fuzzy C-mean clustering algorithm and stochastic diffusion process called mean reverting process. This algorithm quantifies the future states of metal price by the fuzzy series. Operating costs are modelled by Itô-Doob stochastic differential equation. Applying concurrently simulations of these two parameters we can estimate the expected fuzzy value of each TMC for every year of the planning time. After that we discount these values by fuzzy discount rate and define the values of coefficient of objective function. Solution of the fuzzy objective function gives the order of mining of TMCs.
The proposed model is a mathematical representation of mining business reality and allows mining company management to run a dynamic optimization of the business with uncertainty. It helps mining company to survive in very risky environment.
2. Model of Production Planning
Production planning models, based on the linear programming, use block as a basic variable in the objective function. These models also use the constant values of metal price and operating costs through the planning time. It means that these models are static from the point of view of these two parameters. If we want to include fluctuation of these parameters in the objective function, then number of variables significantly increases. Suppose we have a mineral deposit contains of 1 000 blocks, and we want to mine them for 10 years with a different metal price for every year, then the number of variables is about 10 000. Our model reduces the number of variables in the objective function by creation of TMCs. It means that mentioned example would have only 100 variables obtained as years of mining to the power of two. This reduction becomes more significant when dimensions of the block are small. By decreasing the number of variables, we enable uncertainties to be included in the model. We believe that including of uncertainties is much more important than maximum value of the objective function obtained by using blocks as variables with constant values of influencing parameters.
Applying fuzzy set theory and simulation of different stochastic processes we increased flexibility of the model and made the problem more realistic. The model was tested on a small hypothetical lead-zinc mineral deposit and results showed that model can be used for solving the problem of mine production planning.
2.1. Basic Concepts of the Fuzzy Linear Programming
Fuzzy set theory, introduced by Zadeh, deals with problems in which a source of vagueness is involved and has been utilized for incorporating imprecise data into decision framework [16, 17].
The characteristic function of a crisp set assigns a value either 0 or 1 to each member in X. This function can be generalized to a function such that value assigned to the element of the universal set X falls within a specified range, i.e., . The assigned value indicates the membership grade of the element in the set A. The function is called the membership function and the set defined by for each is called a fuzzy set [18, 19].
A fuzzy number is said to be a triangular fuzzy number if its membership function is given byFor more details of arithmetic operations on triangular fuzzy numbers, see [16, 18].
The absolute value of the triangular fuzzy number is denoted by and defined as follows [16]:A ranking function is a function , where F(R) is a set of fuzzy numbers defined on set of real numbers, which maps each fuzzy number into real line, where a natural order exists. Let be a triangular fuzzy number then [18]Linear programming is one of the most frequently applied operations research techniques. In the conventional approach value of the parameters of linear programming models must be well defined and precise. However, in real world environment, this is not realistic assumption. In the real-life problems, there may exist uncertainty about the parameters. In such a situation, the parameters of linear programming problems may be represented as fuzzy numbers [18].
In this paper, we consider zero-one linear programing problem in which only coefficients in the objective function are triangular fuzzy numbers. Such problem is first converted into an adequate crisp model and after that being solved by one of the existing methods.
Suppose we have a linear programming problem with fuzzy coefficients as follows:subject toSince variables xj and coefficients qij are crisp values, it is necessary only to convert fuzzy objective function into crisp function. The process of conversion is based on the way developed by Kumar et al. [18]. Fuzzy objective function may be expressed as follows:
Example 1. Let us consider the following fuzzy objective function and convert it by the proposed method:The fuzzy objective function may be written as follows:Using arithmetic operations on triangular fuzzy numbers and (3), the fuzzy objective function is converted into the following crisp objective function:Suppose the solution of this objective function, with respect to a given set of constraints, is x1=x3=1 and x2=0. Fuzzy optimal value of our objective function is obtained by putting x1, x2, and x3 in (8). The value of the given objective function isAn important concept related to the applications of fuzzy numbers is defuzzification, which converts a fuzzy number into a crisp value. Such a transformation is not unique because different methods are possible. The most commonly used defuzzification method is the centroid defuzzification method, which is also known as center of gravity or center of area defuzzification. The centroid defuzzification method can be expressed as follows (Yager 1981) [20]:where is the defuzzified value. The defuzzification formula of triangular fuzzy number isThis formula will be used in this paper. Defuzzified (crisp) value of our objective function is .
2.2. The Model
In general terms, any economic evaluation of a mine production plan is defined by its financial outcomes. Production planning model, from the economic point of view, is defined by the following objective function and set of constraints:subject towhere is fuzzy present value of the technological mining cut, discounted value, is fuzzy economic value of the technological mining cut in time t, is binary variable, which equals 1 if and only if the technological mining cut i is mined in time t, N is number of technological mining cuts considered for planning. It is equal to the number of mining years (T), S is set of technological mining cuts whose are not accessible from xit, i.e., this set is composed of nonneighbouring mining cuts, is fuzzy discount rate, T is the planning time, is fuzzy capital investment (capital costs).
Objective function represents the present value of the future cash flow of production plan. Equations (15) and (16) ensure that each technological mining cut can be mined only once over the planning time. Equation (17) defines temporal-spatial nonconnectivity between cuts from t to t+1 and ensures the concentration of production. Equation (19) represents the decision making constraint. A positive value of (19) indicates that the obtained production plan is profitable and should be accepted.
Model is to maximize the present value of the production plan which is generated by mining the technological mining cuts over the planning time.
2.3. Creation of Technological Mining Cuts
Room and pillar mining method is designed for flat bedded deposits of limited thickness. This method is used to recover resources in open stopes. The method leaves pillars to support the hanging wall; to recover the maximum amount of ore, miners aim to leave the smallest possible pillars. Rooms and pillars are normally arranged in regular patterns. Pillars can be designed with circular or square cross sections; see Figure 1. Minerals contained in pillars are nonrecoverable and therefore are not included in the ore reserves of the mine [21].
Technological mining cut (TMC) is a part of mineral deposit characterized with respect to the given set of technological requirements (criteria). It means that TMC is a multiattribute object. Suppose the mineral deposit is divided into finite number of mineable blocks. The first step in the production planning is related to the creation of TMCs having the value of attributes closely to the values of technological requirements. Hence, the first step concerns partition of the deposit in adequate number of TMCs. To create the process of production planning more realistic, we apply the concept of fuzzy set theory for some input data. By this approach, uncertainties of input data are decreased and planning becomes much more flexible. To create such TMCs (clusters) we developed fuzzy multicriteria clustering algorithm which is based on Technique for order preference by similarity to ideal solution [23] and constrained polygonal spatial clustering algorithm [24–26].
Mining engineers uses a block model of the deposit that represents the deposit as three-dimensional array of blocks. Such model is created by applying geostatistical methods on data obtained by exploration drilling. In the process of production planning, block is defined as a basic object.
Mineral deposit can be represented as a set of mineable blocks, , and each block is characterized by the block attribute vector , where H is the total number of blocks and A total number of attributes, such as block tonnage, ore grade, etc.
A set is defined as a subset of B and called technological mining cut.
Each TMC is characterized by the mining cut attribute vector , where K is the total number of attributes and it is equal to the total number of technological criteria. Vector of technological criteria is defined by the mine production planer (decision maker); .
At last, creation of technological mining cuts can be mathematically formulated as a multiobjective partition problem, wherein the TMCs must meet technological performance criteria, subject to the given criteria constraints:subject towhere is the ultimate relative closeness (URC) of the ith technological mining cut to the positive ideal technological solution, taking into account all technological criteria, is the lower bound of value of the cth technological criterion, is the upper bound of value of the cth technological criterion.
Note, some of criteria can be excluded from the set of criteria constraints; it depends on the nature of the criterion. Solution of this problem is given as follows:Creation of the set is the two-stage fuzzy multicriteria clustering process which can be treated as two-stage multicriteria decision making process. At first stage we make decision on which cluster (TMC) is to grow, while at second which block should be added to the selected cluster. These two stages represent the one iteration, and the process of clustering is iteratively repeated until no free mineable blocks in the deposit.
Maximization of the ultimate similarity between vector and , where represents the required technological vector, is essential to clustering mineable blocks. Measure of similarity is expressed by the relative closeness coefficient. It is calculated by technique for order preference by similarity to ideal solution (TOPSIS). For detailed description of the method see [23, 27–30], and we have given brief description of its application, in the context of clustering.
Suppose that we defined number (N) of technological mining cuts. The first stage problem that considers which TMC is to grow can be concisely expressed by the following decision making matrix:where is the estimated value of technological mining cut TMCi with respect to the technological criteria Cj. Note that there is a difference between required technological vector and vector of technological criteria, , and it will be explained latter.
For simplicity of notation, we expressed all values as triangular fuzzy numbers, but some of them can be expressed as crisp value. The weighted normalized decision making matrix is computed by multiplying normalized value of with weights () of technological criteria.To avoid decision maker’s subjectivity about weights of criteria, we applied concept of the entropy method [31, 32]. Entropy value of each criterion can be calculated as follows:The objective weight for each criterion is given by the following equation:The fuzzy positive ideal solution and negative one is defined as where J = = 1,2,...,K j associated with criteria that should be ; = = 1,2,...,K j associated with criteria that should be .
The distance from each TMC to and is calculated according to the following equations:where is the distance measurement between two fuzzy triangular numbers, calculated as follows: The relative closeness coefficient of each TMC is calculated asDecision on which TMC is to grow is making according to the following selection rule:The second stage problem that considers which block should be added to the selected cluster (TMCgrow) is defined by the following form:where is the new relative closeness coefficient of the new TMC obtained after adding the mth neighbouring block to the TMCgrow and neighbouring block is the block that has at least one common edge with TMCgrow, is penalty or cost function, M is number of neighbouring mineable blocks.
The new relative closeness coefficients are calculated by applying TOPSIS on the following decision making matrix:where is the new estimated value of the ; with respect to the technological criteria Cj.
The value is estimated after making the union of the attribute vector of the and the attribute vector of the neighbouring block:With the use of the cost function our objective is to select a cluster to be grown that will preserve the maximum degree of flexibility for the other clusters to grow. In order to select a cost function that measures the reduction in flexibility on the growth of the clusters, we observe the effect of the growth of one cluster on the ability of growth of the other clusters. This cost function is as follows [24]:where is number of mineable blocks surrounding the TMCgrow before adding the new mth block, is number of mineable blocks surrounding the TMCgrow after adding the new mth block. is number of common edges between TMCgrow and block to be added, is number of common edges between block to be added and remaining TMCs, is number of common edges between block to be added and waste blocks. Waste block is block that has not grade.
When we define the set and it is necessary to meet the following two spatial constraints:(i)only blocks having at least one common edge with TMCgrow can be added to the TMCgrow,(ii), must not be divided in two or more parts; i.e., technological mining cut must be always homogeneous,
The second constraint means that any mineable block that violates the spatial homogeneity of any TMC cannot be element of and , respectively. Suppose the TMC1 is selected to grow; see Figure 2.
According to spatial constraints only blocks 20, 32, and 33 can be elements of the set because block 26 violates spatial homogeneity of the TMC2. If we add block 33 to the TMC1 than only blocks 20, 32, and 34 can be elements of the set because block 46 violates spatial homogeneity of the TMC3. Hence, for value of the penalty function is equal to 0. Note, mutual takeover of blocks by TMCs is allowed, but spatial homogeneity of TMCs must be always preserved.
In our model, the block attribute vector is composed of the following components:where is ore tonnage in block h (t), expressed as fuzzy triangular number, is grade of the hth block with respect to the γth metal (%), γ is total number of metal concentrates beneficiated from the ore. For polymetallic ore, .
The mining cut attribute vector is composed of the following components:where is ore tonnage in the TMC (t), expressed as fuzzy triangular number, is compactness of the TMC, is standard deviation of the grade in the TMC, with respect to the γth metal (%).
Ore tonnage in the TMC represents a total sum of ore tonnage in blocks contained within TMC:Compactness of the TMC is defined as the following ratio: of the square of the perimeter and the area of the TMC,where is perimeter of the ith TMC, is area of the ith TMC, is total number of blocks in the ith TMC, is length of the block edge (m).
Standard deviation of the grade in the TMC, with respect to the γth metal is calculated as follows:The required technological vector includes the following components: where is annual capacity of production (t/year), expressed as fuzzy triangular number, is desired or target value of compactness of TMC, and it is set up to 16, is standard deviation of the grade, with respect to the γth metal.
Annual capacity of production represents the quantity of ore that should be mined for one year. It is calculated as a total sum of ore tonnage in blocks divided by the total number of planning periods (number of technological mining cuts):Target value of compactness of the TMC is expressed by the Schwartzberg's index of the simple square geometric shape [33]: where e is the edge of the square or mineable block. Standard deviation of the grade, with respect to the γth metal corresponds to the standard deviation of the grade in the TMC, with respect to the γth metal. It is calculated as follows and there is no target value for this component:The main aim of the vector of technological criteria is to enable creation of TMCs so that the ultimate similarity between vector and is maximized. It means that each TMC must meet technological requirements as maximum as possible. Vector of technological criteria is composed of the following components:where is absolute distance between annual capacity of production and ore tonnage in the TMC (t), is absolute distance between target value of compactness and compactness of TMC, is standard deviation of the grade in the TMC with respect to the γth metal (%).
Value of technological mining cut TMCi with respect to the technological criteria is calculated as follows, and it should be minimized:Value of technological mining cut TMCi with respect to the technological criteria is calculated as follows, and it should be minimized:Value of technological mining cut TMCi with respect to the technological criteria is as follows, and it should be minimized:The same calculations are applied when we define values of the decision making matrix .
Set of the technological criteria constraints is composed of only constraint related to the criterion : where is defuzzified value of the ore tonnage in the TMC and is defuzzified value of the annual capacity of production (see (13)).
According to the graph theory, in this paper, the mineral deposit (MD) composed of all mineable blocks is conceived as a graph MD=(B,E), where is a set of mineable blocks, and is a set of edges with euv representing the common edge between block u and block v [34, 35]. The problem of creation of TMCs is solved by the multicriteria partitioning of graph MD. This approach attempts to address situation in which the creation of TMCs should simultaneously maximize the ultimate relative closeness (URC) of the each TMC to the technological requirements, with respect to technological criteria and constraints. Given a set , and total number of TMCs, , then the model of creation of TMCs can be formulated as follows:subject towhere represents the union of mineable blocks contained in the ultimate ith TMC. Constraint (55) forces that a mineable block can belong to one and only one TMC. It means no intersection between TMCs, . Constraint (56) is already explained, and it does not allow any TMC to be terminated, . Constraint (57) is related to binary variable , which equals 1 if and only if the mineable block h belongs to the ith TMC.
Here we present the algorithm of partitioning the graph MD, and it is based on the constrained polygonal spatial clustering algorithm [24–26]. Code of the algorithm is represented by the Algorithm 1.
| 1. Set iteration to ξ=1; | |
| 2. Create a set; | |
| 3. Create a set; | |
| 4. Select the best technological mining cut (BTMC) to grow; | |
| 5. Create the list of neighbouring mineable blocks (NMB) as candidates for the grow of | |
| BTMC such that: | |
| 5.1 block has at least one common edge with TMCgrow | |
| 5.2 none of the homogeneous of the remaining TMCs are violated | |
| 5.3 if: | |
| for ξ+1, | |
| for ξ+2, | |
| for ξ+3, | |
| This is to avoid algorithm gets stuck. | |
| 6. Select the best mineable blocks (BMB) from NMB to add to BTMC; | |
| 7. Add BMB to BTMC and create new technological mining cut; | |
| 8. Update the state of the mining cut attribute vector of the NTMC | |
| 9. If: | |
| 9.1 no free mineable blocks | |
| 9.2 | |
| 9.3 → converge | |
| than Stop | |
| 10. Else go to Step 1 and set ξ+1 |
Our two-stage fuzzy multicriteria clustering algorithm is used to partition the mineral deposit into adequate number of parts such that each part satisfies technological requirements given by the mine production planner. These parts are called technological mining cuts. Algorithm represents the iterative process which starts from the initial state with the aim of approaching the desired goal. Initialization stage means the selection of N technological mining cuts from the set of mineable blocks having at least one common edge with waste. These blocks are located along the perimeter of the mineral deposit (see Figure 2; blocks 18, 19, 24, 32, and 45 are candidates to be selected as initial TMC). Obviously, the initial TMC is composed of one and only one mineable block. From the set of peripheral blocks we select the first N blocks according to the decreasing order of values of relative closeness coefficient. Once the initial TMCs are selected, they begin to be alive and the process of TMCs growing can run. Each TMC is grown by adding neighbouring blocks to him one by one until the desired state of TMC is achieved. Growing of TMCs is an iterative process in nature. At the beginning of the each iteration we first select which TMC is to be grown. Selection is based on the measuring the multidimensional distance between current and desired state of the each TMC. It is expressed by the relative closeness coefficient and TMC having the largest value of RCC is the best to be grown. Upon the selection we proceed to pursue which block is the best to be added. From the set of neighbouring blocks we select the best one by the same approach we have done in the selection of the best TMC but with adding the penalty function. Mutual takeover of the block between two neighbouring TMCc is allowed but homogeneity of each TMC must always be preserved. During the process of growing, the infinity mutual takeover can arise. If the block has mutually been acquired between two TMCs in the three successive iterations we can say the block oscillates between them. In that case the algorithm gets stuck and the local optimum is achieved. To enable algorithm to get global optimum we exclude oscillating block from the set of neighbouring blocks and algorithm can goes on. After that, excluded block is coming back into process. Once the block was added, the mining cut attribute vector must be updated.
Stopping conditions 9.1 and 9.2 for the developed algorithm represent the point when algorithm can be executed, while conditions 9.3 allows algorithm to continue fine tuning. The process of growing goes on until no free mineable blocks and desired technological state of the each mining cut is achieved with respect to given errors.
2.4. Economic Value of the Technological Mining Cut
Generally, the economic estimation of the TMC is based on the three following main components: metal price, costs, and discount rate.
One of the most influencing variables on the economic value of TMC is the metal price. This variable belongs to the set of external variables and cannot be managed by the planners. It is primarily governed by the metal market behaviour. Ability to define the law of variable behaviour trough the time can help planners to find out much more efficient and realistic solutions. By this way we also increase the flexibility in the process of decision making. For that purpose, we developed forecasting algorithm which is based on the combination of fuzzy C-mean clustering and mean reverting process.
Consider a metal price and denote it as variable X. If the value of that variable is governed by the laws of probability, then variable X can be treated as a stochastic variable. Suppose that we monitored values of X at equal interval, ; . In this paper interval of monitoring is one month and we use symbols and T to make distinction between interval of monitoring and interval of planning (one year). Such a sequence of monitored values, , is called stochastic time series. If we assign some underlying probabilistic distribution to the time series then it becomes stochastic process.
Model of forecast does not provide the exact point estimation (crisp value) of variable, but rather the fuzzy state that the variable will be at the next point; i.e., model generates the future sequence of fuzzy states. The general concept is as follows: the monitored time series of metal price is transformed into fuzzy state series by applying the fuzzy C-mean clustering algorithm, while the future states are forecasted by stochastic diffusion process called mean reverting process. The goal of the forecasting model is to estimate the fuzzy state that the metal price will fall within one of the a priori defined states. This model is able to account for the dynamics of metal price process and distinguish increasing from a decreasing period. Therefore, efficiency of the model directly depends on the use of relevant monitored information pertaining to this process. Usage of the monitored information is primary related to the calibration of the model, i.e., to the defining of the parameters that will govern the forecasting process.
Fuzzification of monitored metal price time series is performed by fuzzy C-mean clustering algorithm. This algorithm belongs to the partitioning methods that consist of dividing N objects into a specified number of M disjoint groups that are also called classes or clusters. Fuzzy C-mean algorithm is based on minimization of the following least-squared errors function: subject towhere X is the vector of monitored metal prices; , C is the vector of class centers; , U is the fuzzy partition matrix; , ω is the coefficient of fuzzification and we take value equal to 2.
The objective function F is iteratively minimized. The iteration process stops until , where j is the number of iteration and represents the minimum amount of improvement. For more details see [36–39]. Suppose that partitioning F(X,C) has been done and the sequence of obtained centers is sorted in an ascending order; .
To describe the value (level) of metal price we use linguistic variables. A linguistic variable is variable whose values are words or sentences in natural artificial language [40]. Following seven linguistic variables are used for that purpose: very very low (VVL), very low (VL), low (L), medium (M), high (H), very high (VH) and very very high (VVH). Accordingly, number of clusters equals also seven (M=1,2,…,7). Range boundaries of each linguistic variable are defined by the transformation of linguistic variable into adequate triangular fuzzy number and corresponding range code. Range code is expressed as crisp number.
Triangular fuzzy number is defined as a triplet , where parameters a, b, and c, respectively, indicate the smallest possible value, the most promising value, and the largest possible value. This formulation is interpreted as membership function and holds the following conditions: ) ; membership function is increasing in the interval and decreasing in the interval . Range boundaries and value of code of linguistic variables are shown in Table 1.
Suppose a monthly monitored metal price time series is given and a range code vector is defined on it. Transformation of monitored series to the fuzzy state series (range code series) is performed according to the following tagging mechanism:At last, the monitored series is expressed as a range code time series: For example, . Uncertainties related to behaviour of future states of the metal price are modelled by a simulation of stochastic diffusion process, called mean reverting process. The main characteristic of this process is the tendency of the commodity prices to revert to a long-term equilibrium level after significant short-term fluctuations. Underlying process is described by the following equation [41]: where is long-term equilibrium fuzzy state expressed by the corresponding range code value, is speed of mean reverting to the long-term equilibrium range code value, is rate of range code volatility, is an increment to a standard Brownian motion.
The correct discrete-time format for the continuous-time process of mean reverting is the stationary first order autoregressive process [42] so the sample path simulation equation for is performed by using exact discrete-time expression as follows:where is constant time interval from point to , one month, is normally distributed random variable.
Calibration of the mean reverting process is based on the following regression:The speed of mean reverting () equals the negative of the slope, long-term equilibrium () is the intercept divided by the speed, and the rate of range code volatility () is standard deviation from the regression.
Let denote an evolution path of range code values with spot values , where is calculated by (65). Simulation of (65), for every , is composed of the following two steps: (a) do calculation of (65) and (b) for obtained result determine the range code value. Suppose the is the value obtained by (65) at point , then corresponding range code value is assigned as follows:Note the long-term equilibrium range code value, obtained by the regression, should also be defined by (67).
Space of simulation of future metal price and way of fuzzification is described as follows:where S represents the total number of simulations and T time span of simulation. Figure 3 represents an evolution path of the range code values while Figure 4 distribution of range code values in the month, upon S simulations.
When developing cost estimates it is important to distinguish between the types of costs being estimated. The estimator is concerned with two primary types of costs for project estimation: capital costs and operating costs. In the minerals industry, capital costs or capital investment generally means the amount of total capital dollars required to bring a mining property into production [1]. These costs include underground mine construction and mineral processing plant construction costs and equipment purchase. Underground mine construction and mineral processing plant construction are activities that can last a few years. In that period it is reasonable to expect the fluctuation of capital costs. In such environment, the right-skewed triangular fuzzy number might express uncertainties related to capital costs; .
Operating costs are also very significant source of uncertainties related to the economic value of TMC. Operating costs are incurred directly in the production process. These costs include the ore and waste development of individual stopes, the actual stopping activities, the mine services that provides logistical support to the miners, and the milling and processing of the ore at the plant [43]. Some components of operating cost such as fuel, electricity, lubricants, tyres, replacement parts, and inputs that is used for mineral processing, etc. are usually purchased at market prices that fluctuate monthly. Note when we talk about operating cost we mean unit operating cost. Uncertainties related to these costs arise as a consequence that many suppliers, that do business in market conditions as well, are offering short-term contracts to mine as a way to protect their business activities. In such environment and if we add the fact that production will be carried out for many years then it is very important to predict the future behaviour of operating costs. For that purpose, we apply a continuous-time process using the following Itô-Doob type stochastic differential equation [44]:where μ and σ are the drift and volatility and is a normalized Brownian motion. Solution of this equation is as follows and it describes behaviour of operating costs:Applying simulation on (70) we can create many artificial scenarios of operating cost behaviour trough the time. Space of simulation of operating costs is represented by the following matrix:Bearing in mind that our model is based on the maximization of the expected present value of production plan then it is obvious that discount rate has significant influence. Determining with precision the appropriate discount rate to use for financial analysis studies is a difficult, much debate topic. It is clear that the cost of acquiring investment funds must cover opportunity costs and transaction costs, must compensate for risk, and must cover anticipated inflation. Although this places some bounds on the problem, each of these items is also difficult to quantify with precision [1]. For example, if we ask experts about value of discount rate, we can expect the following answer: ”discount rate is about 10%”. This statement is burdened with vagueness rather than randomness. It can be treated as linguistic variable which can be quantified by triangular fuzzy number; .
Upon discussion on all main components separately, we have to calculate the economic value of TMC according to their joined influence.
Metal concentrate is a final product of underground mining companies that have not smelting facilities. They realize revenues by selling it on metal concentrate market. Value of the concentrate directly depends on the metal price, metal content of concentrate, and metal recovery rate. If we take into consideration that metal price uncertainty is described by the simulation and fuzzification then the unite value of concentrate is calculated as follows:where is unit selling price of the metal in period and simulation s expressed by the triangular fuzzy number ($/t), mcon is metal content of concentrate (%), mmr is metal recovery rate (%).
Economic value of mineable block is calculated as follows [1]:where Y is the total number of metal concentrates beneficiated from the ore, is mill recovery rate of the γth metal (%).
Space of economic value of mineable block can be represented by the following fuzzy matrix:Expected economic values of mineable block trough the planning time are obtained by averaging each column of the matrix . If we take into consideration that our planning model is based on the annual time span then it is necessary to transform forecasted economic values of mineable block from monthly to annual level. Let be the expected economic value of mineable block through the planning time, then its discounted value for every year is defined as follows: Finally, the coefficient which represents the fuzzy present value of the technological mining cut (discounted value) is calculated by the following equation: It means that discounted economic value of the TMC is a sum of discounted economic values of the mineable blocks contained in TMC.
3. Numerical Example
To test developed production planning model we borrowed data of small hypothetical lead-zinc ore body from Gligoric et al. [15]. Room and pillar mining method was selected as way of mining the mineral deposit. The management of mining company is looking for the best production plan with uncertainties related to the metal price and operating costs. Mineable ore reserves were estimated about 862 245 t. To start up the production it is necessary to do development of ore body and buildup mineral processing plant. All these activities are defined as follows:(i)construction of development system based on the ramp and horizontal drive that enables accessing the ore body. It will be used for purpose of ore haulage and fresh air intake,(ii)construction of horizontal drive and decline that enables exhaust of foul air after having ventilated working places of the mine,(iii)purchase of new mining equipment,(iv)buildup of the mineral processing plant,(v)purchase of mining and processing equipment.
For the purpose of planning, the ore body is represented by 115 blocks, each 25×25×5 m. The input parameters used for testing the developed model are given in Figures 5 and 6 and Tables 2 and 3.
Regular pattern of square pillars means that pillars are left unmined and distance between them is equal.
Applying the clustering algorithm on the set of peripheral blocks we selected five blocks according to the decreasing order of values of relative closeness coefficient. If we take into consideration that all peripheral blocks have the same compactness and standard deviation of the zinc and lead grade are equal to zero, then selection is based only on the absolute distance between annual capacity of production and ore tonnage in the peripheral block. By this way, initial TMCs (seeds) were defined; see Figure 7.
Having selected initial TMCs algorithm proceeded to do iterations and characteristically phases of transition are presented by Figure 8.
(a)
(b)
(c)
(d)
Finally, algorithm partitioned the original mineral deposit into five ultimate technological mining cuts with respect to the error of annual capacity of production and convergence of standard deviation of lead and zinc in the TMC (see Figure 9).
Stopping condition 9.2 related to the standard deviation of lead and zinc in the TMC is represented by Figure 10.
Realized technological mining cut attribute values that consider ore tonnage and compactness of the TMC are shown in Table 4.
Mean absolute percentage error of the production capacity isand it indicates that efficiency of the algorithm is very high. If we take into consideration the fact that mineral deposit is of very complex shape then compactnesses of TMCs achieved by algorithm are suitable. Algorithm succeeded to create three TMCs that have the similar compactness and it is 60% of the original deposit. Summary characteristics that consider the grade in the TMC are shown in Table 5. These characteristics are very important for the mineral processing planners.
Grade distribution is related to ore body genesis processes and cannot be influenced by planners. Coefficient of lead and zinc grade variation for entire mineral deposit is 17.2 and 32.7%, respectively. Algorithm also succeeded to divide ore body in parts with suitable coefficients of grade variation. We can see that only in TMC1 and TMC4 coefficient of zinc grade variation is greater than 32.7% but not so significantly. Realized distribution of lead and zinc grade in TMCs enables mineral processing planners to create much more flexible plans because they are informed in advance about expected grade distribution.
Since the all attributes of TMCs are defined we can now proceed to define their economic values. Two monitored time series represented by Figure 11 have been used to define classes of lead and zinc metal prices. Period of monitoring was 125 months. Fuzzy C-mean clustering algorithm produced the following two vectors of class centers; see Table 6.
(a)
(b)
Minimum and maximum values of the lead price are 963 $/t and 3080 $/t while for zinc they are 1101 $/t and 3533 $/t. Transformation of linguistic variables into adequate triangular fuzzy numbers and corresponding range codes according to obtained class centers is represented by Table 7 and Figure 12.
(a)
(b)
Transformation of monitored series to the fuzzy state series (range code series) is performed according to the tagging mechanism (see (62)). Monitored value was fuzzified with respect to where the maximum membership degree occurred. Results of tagging mechanism are represented by Figure 13.
(a)
(b)
(c)
(d)
To calibrate the mean reverting process, which is used to forecast the future states of lead and zinc metal price, we have run first order linear regression on range code lead price series and range code zinc price series separately. Calibration of Itô-Doob stochastic process, which is used to forecast the future states of unit operating costs, is based on the expert knowledge. Parameters needed for simulation of these diffusion processes are shown in the Table 8.
One simulation of lead and zinc metal price, unit operating cost and resulting economic value of mineable block 1, is represented in Figure 14.
After 500 simulations we obtained sample which is used to defined expected economic values of mineable block 1 through the planning time on monthly span. For every month we average the sample and expected economic values of mineable block 1 are represented by Figure 15.
Transformation of forecasted economic values of block 1, from monthly to annual level, is represented by Figure 16.
Discounted or present economic values of mineable block 1 are represented by Figure 17.
Applying the process of simulation on all blocks simultaneously, we obtain their present economic values and we can calculate the present economic values of TMCs. According to the results of mineral deposit partition, following present economic values of every TMC are shown in Table 9.
By this way, every coefficient , of the objective function (see (14)) is defined. Solution of the fuzzy objective function is x14=1; x22=1; x31=1; x45=1; x53=1 while the remains are equal to zero. It means that optimal plan of production is TMC3 (year 1) → TMC2 (year 2) → TMC5 (year 3) → TMC1 (year 4) → TMC4 (year 5) (see Figure 18).
Value of the objective function is . Value of (19) isThis difference represents the net present value of the realized production plan. Since it is a positive, production plan is profitable and should be accepted.
4. Conclusion
Having ability to optimize the production with uncertainty is recognized as critical to compete and even survive of mining company. The presence of the uncertainty makes that the model works in a dynamic environment producing the optimal results and plans. Metal price, operating costs, and annual capacity of production are identified as main sources of uncertainties and they are described by fuzzy numbers. Each mining project depends on these uncertain parameters which have a significant impact on the final financial results. By this way we incorporate risks into mine production planning.
In this paper we have presented fuzzy 0-1 linear programming model for production planning of room and pillar underground mining method. Model is based on the maximization of fuzzy objective function which represents the present or discounted value of the future cash flow of production plan, with respect to the set of constraints.
Our approach of the mine production planning contains two main parts. First part is referred on partitioning of the mineral deposit creating the TMCs (clusters) of the ore body. Two-stage fuzzy multicriteria clustering algorithm as iterative process is applied for creation of TMCs where each cluster meets the following technological requirements: annual capacity of production, compactness of the shape of TMC, and standard deviation of ore grade in the TMC.
The second part of the mine production planning model is related to calculation process of the economic value of mineable blocks using the simulation of the mean reverting process and Itô-Doob stochastic differential equation. Economic value of each mineable block is defined by metal price and operating costs with uncertainty. Accordingly, we developed concurrently simulations of these two main economical sources whereby each mineable block has the different economic value through the planning time. Based on this, we can estimate the present economic values (expected fuzzy values) of each TMCs for every year of the planning time.
Coefficients of the fuzzy linear objective function are obtained by discounting these expected values by fuzzy discount rate. Solution of the fuzzy objective function gives the order of mining of TMCs, i.e., the optimal production plan with respect to operational constraints. Bearing in mind that the value of the objective function is positive, our developed algorithm is absolutely realized and accepted for the solving more complex problems.
The model is not closed and can be extended by including the uncertainty related to the ore grade. This uncertainty is very important and significantly affects the production plan. Fuzzy numbers can also be used for this purpose but calculation of standard deviation of fuzzy numbers is a very complex task. Future exploration will be focused on the inclusion of fuzzy numbers as a way to describe ore grade uncertainty.
Data Availability
No data were used to support this study.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.