# Long-Term Room and Pillar Mine Production Planning Based on Fuzzy 0-1 Linear Programing and Multicriteria Clustering Algorithm with Uncertainty

Authors: Miloš V. Gligorić; Zoran M. Gligorić; Čedomir R. Beljić; Suzana M. Lutovac; Vesna M. Damnjanović
Journal: Mathematical Problems in Engineering
 (2019)

---

## 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.


---

## Body

## 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 ofTMC, and standard deviation of ore grade in theTMC. 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 ofTMCs having the value of attributes closely to the values of technological requirements. In the purpose of creation suchTMCs (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 ofTMCs. 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 eachTMC 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 ofTMCs.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 ofTMCs. 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μA of a crisp set A⊆X assigns a value either 0 or 1 to each member inX. This function can be generalized to a function μA~ such that value assigned to the element of the universal setX falls within a specified range, i.e., μA~:X→0,1. The assigned value indicates the membership grade of the element in the setA. The function μA~ is called the membership function and the set A~=x,μA~x:x∈X defined by μA~x for each x∈X is called a fuzzy set [18, 19].A fuzzy numberA~=a,b,c is said to be a triangular fuzzy number if its membership function is given by(1)μA~x=x-ab-a,a≤x≤bx-cb-c,b≤x≤c0,otherwiseFor more details of arithmetic operations on triangular fuzzy numbers, see [16, 18].The absolute value of the triangular fuzzy numberA~=a,b,c is denoted by A~ and defined as follows [16]:(2)A~x=max⁡A~x,-A~x,ifx≥00,ifx<0A ranking function is a function R:FR→R, whereF(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 A~=a,b,c be a triangular fuzzy number then [18](3)RA~=a+2b+c4Linear 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:(4)max∑j=1naj,bj,cjxjsubject to(5)∑j=1nqijxj≤pi,1≤i≤m;1≤j≤n(6)xj∈0,1Since variablesxj and coefficientsqij 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:(7)maxR∑j=1naj,bj,cjxjExample 1.
Let us consider the following fuzzy objective function and convert it by the proposed method:(8)max1,6,9x1+2,3,8x2+4,6,7x3The fuzzy objective function may be written as follows:(9)maxR1x1+2x2+4x3,6x1+3x2+6x3,9x1+8x2+7x3Using arithmetic operations on triangular fuzzy numbers and (3), the fuzzy objective function is converted into the following crisp objective function:(10)max1422x1+16x2+23x3Suppose the solution of this objective function, with respect to a given set of constraints, isx1=x3=1 andx2=0. Fuzzy optimal value of our objective function is obtained by puttingx1,x2, andx3 in (8). The value of the given objective function is(11)1,6,9∙1+2,3,8∙0+4,6,7∙1=5,12,16An 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]:(12)x-0A~=∫acxμA~xdx∫acμA~xdxwhere x-0A~ is the defuzzified value. The defuzzification formula of triangular fuzzy number A~=a,b,c is(13)x-0A~=a+b+c3This formula will be used in this paper. Defuzzified (crisp) value of our objective function is x-05,12,16=11.


### 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:(14)F~V~,X=∑i=1N∑t=1Tv~itdxit=∑i=1N∑t=1Tv~it1+d~txit→maxsubject to(15)∑t=1Txit=1,∀i∈N(16)∑i=1Nxit=1,∀t∈T(17)xit+∑s=1Sxs,t+1≤1,∀i∈N,∀t∈T,s∈S(18)xit∈0,1,∀i∈N,∀t∈T(19)F~-I~≥0where v~itd is fuzzy present value of the technological mining cut, discounted value,v~it is fuzzy economic value of the technological mining cut in timet,xit is binary variable, which equals 1 if and only if the technological mining cuti is mined in timet,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 fromxit, i.e., this set is composed of nonneighbouring mining cuts,d~ is fuzzy discount rate,T is the planning time,I~ 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 fromt tot+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 Figure1. Minerals contained in pillars are nonrecoverable and therefore are not included in the ore reserves of the mine [21].Figure 1
Room and pillar mining method [22].Technological mining cut (TMC) is a part of mineral deposit characterized with respect to the given set of technological requirements (criteria). It means thatTMC 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 ofTMCs having the value of attributes closely to the values of technological requirements. Hence, the first step concerns partition of the deposit in adequate number ofTMCs. 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 suchTMCs (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,B=bhh=1,2,…,H, and each block is characterized by the block attribute vector BAVh=ajhj=1,2,..,A;∀h∈H, whereH is the total number of blocks andA total number of attributes, such as block tonnage, ore grade, etc.A setTMCi=blil<H;∀i∈N is defined as a subset ofB and called technological mining cut.EachTMC is characterized by the mining cut attribute vector MCAVi=aicc=1,2,…,K;∀i∈N, whereK 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); C=Cjj=1,2,…,K.At last, creation of technological mining cuts can be mathematically formulated as a multiobjective partition problem, wherein theTMCs must meet technological performance criteria, subject to the given criteria constraints:(20)max︸a1c∈TMC1fa1c=f1a11,f2a12,…,fKa1Kmax︸a2c∈TMC2fa2c=f1a21,f2a22,…,fKa2K⋮max︸aic∈TMCifaic=f1ai1,f2ai2,…,fKaNKsubject to(21)acl≤ac≤acu,c=1,2,…,Kwherefaic is the ultimate relative closeness (URC) of theith technological mining cut to the positive ideal technological solution, taking into account all technological criteria,acl is the lower bound of value of thecth technological criterion,acu is the upper bound of value of thecth 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:(22)B=⋃i=1NTMCi=TMC1∪TMC2∪…∪TMCiTMCα∩TMCβ=∅,α≠βCreation of the set TMCi=blil<H;∀i∈N 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 vectorMCAVi=aicc=1,2,…,K;∀i∈N and Cteh=cjtehj=1,2,…,K, where Cteh 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 whichTMC is to grow can be concisely expressed by the following decision making matrix:(23)D~=x~ijN×K=TMCCC1C2⋯CjTMC1x~11x~12⋯x~1jTMC2x~21x~22⋯x~2j⋮⋮⋮⋱⋮TMCix~i1x~i2⋯x~ijwhere x~ij is the estimated value of technological mining cutTMCi with respect to the technological criteriaCj. Note that there is a difference between required technological vector and vector of technological criteria, Cteh≠C, 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 ofx~ij with weights (w~j) of technological criteria.(24)P~=p~ijN×K=r~ij∙w~jN×K=x~ij∑i=1Nx~ij∙w~jN×K,∀j∈1,K;∑j=1Kw~j=1To 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:(25)e~j=-1ln⁡N∑i=1Nr~ij∙ln⁡r~ij,∀j∈1,K;0≤e~j≤1The objective weight for each criterion is given by the following equation:(26)w~j=1-e~j∑j=1K1-e~j,∀j∈1,K;0≤w~j≤1The fuzzy positive ideal solution A~+ and negative one A~- is defined as (27)A~+=max︸i=1,2,…,Np~ij∣j∈J,min︸i=1,2,…,Np~ij∣j∈J′=p~1+,p~2+,…,p~K+(28)A~-=min︸i=1,2,…,Np~ij∣j∈J,max︸i=1,2,…,Np~ij∣j∈J′=p~1-,p~2-,…,p~K-(29)max︸i=1,2,…,Np~ij=max︸i=1,2,…,Naij,bij,cij=max︸i=1,2,…,Nbij,∀j∈1,K(30)min︸i=1,2,…,Np~ij=min︸i=1,2,…,Naij,bij,cij=min︸i=1,2,…,Nbij,∀j∈1,KwhereJ = {j = 1,2,...,K∣j associated with criteria that should be maximized};J′ = {j = 1,2,...,K∣j associated with criteria that should be minimized}.The distance from eachTMC to A~+ and A~- is calculated according to the following equations:(31)di+=∑j=1Kdp~ij,p~j+,∀i∈1,N(32)di-=∑j=1Kdp~ij,p~j-,∀i∈1,Nwhere dp~ij,p~j is the distance measurement between two fuzzy triangular numbers, calculated as follows: (33)dp~ij,p~j=13aij-aj2+bij-bj2+cij-cj2The relative closeness coefficient of eachTMC is calculated as(34)RCCi=di-di-+di+,∀i∈1,NDecision on whichTMC is to grow is making according to the following selection rule:(35)TMCgrow=maxRCC1,RCC2,…,RCCNThe second stage problem that considers which block should be added to the selected cluster (TMCgrow) is defined by the following form:(36)badd=maxRCC1add+O1add,RCC2add+O2add,+⋯+,RCCMadd+OMaddwhereRCCMadd is the new relative closeness coefficient of the newTMC obtained after adding themth neighbouring block to theTMCgrow and neighbouring block is the block that has at least one common edge withTMCgrow,OMadd 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:(37)D~block=φ~mjM×K=TMCgrow∪bCC1C2⋯CjTMCgrow∪b1φ~11φ~12⋯φ~1jTMCgrow∪b2φ~21φ~22⋯φ~2j⋮⋮⋮⋱⋮TMCgrow∪bmφ~m1φ~m2⋯φ~mjwhereφ~mj is the new estimated value of the TMCgrow∪bm; m=1,2,..,M with respect to the technological criteriaCj.The valueφ~mj is estimated after making the union of the attribute vector of the TMCgrow and the attribute vector of the neighbouring block:(38)φ~mj=MCAVgrow∪BAVm=acgrowc=1,2,…,K∪ajmj=1,2,…,A;∀m∈M;K=AWith 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]:(39)OMadd=Bbef-BmaftBbef=Bmcom+Bmngh+Bmwaste-3Bbef,∀m∈1,MwhereBbef is number of mineable blocks surrounding theTMCgrow before adding the newmth block,Bmaft is number of mineable blocks surrounding theTMCgrow after adding the newmth block.Bmcom is number of common edges betweenTMCgrow and block to be added,Bmngh is number of common edges between block to be added and remainingTMCs,Bmwaste 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 setBbef and Bmaft it is necessary to meet the following two spatial constraints:(i)
only blocks having at least one common edge withTMCgrow can be added to theTMCgrow,(ii)
TMCi,∀i∈1,N, 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 anyTMC cannot be element of Bbef and Bmaft, respectively. Suppose theTMC1 is selected to grow; see Figure 2.Figure 2
Spatial plan ofTMCs.According to spatial constraints only blocks 20, 32, and 33 can be elements of the setBbef because block 26 violates spatial homogeneity of theTMC2. If we add block 33 to theTMC1 than only blocks 20, 32, and 34 can be elements of the set Bmaft because block 46 violates spatial homogeneity of theTMC3. Hence, for TMC1∪b33 value of the penalty function is equal to 0. Note, mutual takeover of blocks byTMCs is allowed, but spatial homogeneity ofTMCs must be always preserved.In our model, the block attribute vectorBAVh=ajhj=1,2,..,A;∀h∈H is composed of the following components:(40)BAVh=a~1h,a2,1h,a2,2h,…,a2,γhj=1+γ;∀h∈Hwherea~1h is ore tonnage in blockh (t), expressed as fuzzy triangular number,a2,γh is grade of thehth block with respect to the γth metal (%),γ is total number of metal concentrates beneficiated from the ore. For polymetallic ore, γ>1.The mining cut attribute vectorMCAVi=aicc=1,2,…,K;∀i∈N is composed of the following components:(41)MCAVi=aicc=1,2,…,K;∀i∈N=a~i1,ai,2,ai,3,1,ai,3,2,…,ai,3,γc=2+γ;∀i∈Nwherea~i1 is ore tonnage in theTMC (t), expressed as fuzzy triangular number,ai2 is compactness of theTMC,ai,3,γ is standard deviation of the grade in theTMC, with respect to the γth metal (%).Ore tonnage in theTMC represents a total sum of ore tonnage in blocks contained withinTMC:(42)a~i1=∑h=1lia~1h,∀i∈1,N,l1+l2+⋯+lN=HCompactness of theTMC is defined as the following ratio: of the square of the perimeter and the area of theTMC,(43)ai2=Pi2Ai=Pi2li∙be2,∀i∈1,NwherePi is perimeter of theith TMC,Ai is area of theith TMC,li is total number of blocks in theith TMC,be is length of the block edge (m).Standard deviation of the grade in theTMC, with respect to the γth metal is calculated as follows:(44)ai,3,γ=1li∑h=1lia-2,γ-a2,γh2,li<H,∀i∈1,N,∀γ∈1,YThe required technological vector Cteh=cjtehj=1,2,…,K includes the following components: (45)Cteh=cjtehj=1,2,…,K=c~1teh,c2teh,c3,1teh,c3,2teh,…,c3,γtehj=2+γwherec~1teh is annual capacity of production (t/year), expressed as fuzzy triangular number,c2teh is desired or target value of compactness ofTMC, and it is set up to 16,c3,γteh 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):(46)c~1teh=1N∑h=1Ha~1hTarget value of compactness of theTMC is expressed by the Schwartzberg's index of the simple square geometric shape [33]: (47)c2teh=4e2e2=16wheree 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 theTMC, with respect to the γth metal. It is calculated as follows and there is no target value for this component:(48)c3,γteh=ai,3,γThe main aim of the vector of technological criteria is to enable creation ofTMCs so that the ultimate similarity between vector MCAVi=aicc=1,2,…,K;∀i∈N and Cteh=cjtehj=1,2,…,K is maximized. It means that eachTMC must meet technological requirements as maximum as possible. Vector of technological criteria C=Cjj=1,2,…,K is composed of the following components:(49)C=Cjj=1,2,…,K=C~1,C2,C3,1,C3,2,…,C3,γj=2+γwhere C~1 is absolute distance between annual capacity of production and ore tonnage in theTMC (t),C2 is absolute distance between target value of compactness and compactness ofTMC,C3,γ is standard deviation of the grade in theTMC with respect to the γth metal (%).Value of technological mining cutTMCi with respect to the technological criteria C~1 is calculated as follows, and it should be minimized:(50)x~i1=a-0a~i1c-0c~1teh×100%,∀i∈1,NValue of technological mining cutTMCi with respect to the technological criteria C2 is calculated as follows, and it should be minimized:(51)xi2=16-ai2,∀i∈1,NValue of technological mining cutTMCi with respect to the technological criteria C3,γ is as follows, and it should be minimized:(52)xi,3,γ=ai,3,γ,∀i∈1,N,∀γ∈1,YThe same calculations are applied when we define values of the decision making matrix D~block=φ~mjM×K.Set of the technological criteria constraintsacl≤ac≤acu,c=1,2,…,K is composed of only constraint related to the criterion C~1: (53)-10%≤a1≤+10%-10%≤1-a-0a~i1c-0c~1teh×100%≤+10%;∀i∈1,Nwhere a-0a~i1 is defuzzified value of the ore tonnage in theTMC and c-0c~1teh 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 graphMD=(B,E), where B=bhh=1,2,…,H is a set of mineable blocks, and E=euv∣u,v=1,2,..,H;u≠v is a set of edges witheuv representing the common edge between blocku and blockv [34, 35]. The problem of creation ofTMCs is solved by the multicriteria partitioning of graphMD. This approach attempts to address situation in which the creation ofTMCs should simultaneously maximize the ultimate relative closeness (URC) of the eachTMC to the technological requirements, with respect to technological criteria and constraints. Given a set B=bhh=1,2,…,H, and total number ofTMCs, N≥2, then the model of creation ofTMCs can be formulated as follows:(54)Z=max︸∀i∈1,NURCi=max︸∀i∈1,Ndi-⋃h=1li<Hbihdi-⋃h=1li<Hbih+di+⋃h=1li<Hbihsubject to(55)∑i=1Nbih=1,∀h∈1,H(56)a~1l≤1-defuzzified∑h=1Ha~1ihbihdefuzzifiedc~1teh×100%≤a~1u,∀i∈1,N(57)bih∈0,1,∀h∈H,∀i∈Nwhere ⋃h=1li<Hbih represents the union of mineable blocks contained in the ultimateithTMC. Constraint (55) forces that a mineable block can belong to one and only oneTMC. It means no intersection betweenTMCs, TMCα∩TMCβ=∅,α≠β. Constraint (56) is already explained, and it does not allow anyTMC to be terminated, TMCi≠∅,∀i∈1,N. Constraint (57) is related to binary variable bih, which equals 1 if and only if the mineable blockh belongs to theithTMC.Here we present the algorithm of partitioning the graphMD, and it is based on the constrained polygonal spatial clustering algorithm [24–26]. Code of the algorithm is represented by the Algorithm 1.Algorithm 1:Algorithm of the TMC creation.
1. Set iteration toξ=1; 2. Create a set;TMCi=TMC1,TMC2,…,TMCN 3. Create a set;C=Cjj=1,2,…,K=C~1,C2,C3,1,C3,2,…,C3,γj=2+γ 4. Select the best technological mining cut (BTMC) to grow; B T M C = T M C g r o w = m a x R C C 1 , R C C 2 , … , R C C N 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 BMB∉NMB if: for ξ+1, BMB∈TMCαgrow for ξ+2, BMB∈TMCβgrow for ξ+3, bh∈TMCαgrow This is to avoid algorithm gets stuck. 6. Select the best mineable blocks (BMB) from NMB to add to BTMC; B M B = m a x R C C 1 a d d + O 1 a d d , R C C 2 a d d + O 2 a d d , + ⋯ + , R C C M a d d + O M a d d 7. AddBMB to BTMC and create new technological mining cut; N T M C = B T M C ∪ B M B 8. Update the state of the mining cut attribute vector of theNTMC 9. If: 9.1 no free mineable blocks 9.2 -10%≤1-a-0a~i1c-0c~1teh×100%≤+10%;∀i∈1,N 9.3 ai,3,γ,∀i∈1,N,∀γ∈1,Y → converge than Stop 10. Else go to Step1 and set ξ+1Our 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 ofN 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 initialTMC). Obviously, the initialTMC is composed of one and only one mineable block. From the set of peripheral blocks we select the firstN blocks according to the decreasing order of values of relative closeness coefficient. Once the initialTMCs are selected, they begin to be alive and the process ofTMCs growing can run. EachTMC is grown by adding neighbouring blocks to him one by one until the desired state ofTMC is achieved. Growing ofTMCs is an iterative process in nature. At the beginning of the each iteration we first select whichTMC is to be grown. Selection is based on the measuring the multidimensional distance between current and desired state of the eachTMC. It is expressed by the relative closeness coefficient andTMC having the largest value ofRCC 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 bestTMC but with adding the penalty function. Mutual takeover of the block between two neighbouringTMCc is allowed but homogeneity of eachTMC must always be preserved. During the process of growing, the infinity mutual takeover can arise. If the block has mutually been acquired between twoTMCs 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 MCAVi=aicc=1,2,…,K;∀i∈N 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 theTMC is based on the three following main components: metal price, costs, and discount rate.One of the most influencing variables on the economic value ofTMC 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 variableX. If the value of that variable is governed by the laws of probability, then variableX can be treated as a stochastic variable. Suppose that we monitored values ofX at equal interval, ∆τi=τi-τi-1=const.; i=2,3,…,T. 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, Xτ=x1,x2,…,xT, 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 dividingN objects into a specified number ofM disjoint groups that are also called classes or clusters. Fuzzy C-mean algorithm is based on minimization of the following least-squared errors function: (58)F=U,X,C=∑τ=1T∑m=1Muτmωxτ-cm2subject to(59)∑m=1Muτm=1,τ=1,2,…,T(60)0≤uτm≤1,m=1,2,…,M,τ=1,2,…,T(61)0<∑τ=1Tuim<m,m=1,2,…,MwhereX is the vector of monitored metal prices; X=x1,x2,…,xT⫅RT,C is the vector of class centers; C=c1,c2,…,cM⫅RM,U is the fuzzy partition matrix; U=uτmT×M,ω is the coefficient of fuzzification and we take value equal to 2.The objective functionF is iteratively minimized. The iteration process stops until Fj+1-Fj<δ, wherej is the number of iteration and δ represents the minimum amount of improvement. For more details see [36–39]. Suppose that partitioningF(X,C) has been done and the sequence of obtained centers is sorted in an ascending order; c1<c2<⋯<cM.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 tripletA~=a,b,c, where parametersa,b, andc, 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: (1) a≤b≤c; (2) membership function is increasing in the interval a,b and decreasing in the interval b,c. Range boundaries and value of code of linguistic variables are shown in Table 1.Table 1
The transformation of linguistic variables to triangular fuzzy numbers and range code.
Linguistic Triangular fuzzy number Membership Range code value variable function y m , m = 1,2 , … , 7 VVL a < c 1 , b = c 1 , c 1 < c ≤ c 1 + c 2 2 0 ; 1 ; 0 y 1 = 1  VL c 1 + c 2 2 < a , b = c 2 , c 2 < c ≤ c 2 + c 3 2 0 ; 1 ; 0 y 2 = 2  L c 2 + c 3 2 < a , b = c 3 , c 3 < c ≤ c 3 + c 4 2 0 ; 1 ; 0 y 3 = 3  M c 3 + c 4 2 < a , b = c 4 , c 4 < c ≤ c 4 + c 5 2 0 ; 1 ; 0 y 4 = 4  H c 4 + c 5 2 < a , b = c 5 , c 5 < c ≤ c 5 + c 6 2 0 ; 1 ; 0 y 5 = 5  VH c 5 + c 6 2 < a , b = c 6 , c 6 < c ≤ c 6 + c 7 2 0 ; 1 ; 0 y 6 = 6  VVH c 6 + c 7 2 < a , b = c 7 , c 7 < c 0 ; 1 ; 0 y 7 = 7Suppose a monthly monitored metal price time series is givenP=pττ=1,2,…,T and a range code vector Y=ymm=1,2,…,M is defined on it. Transformation of monitored series to the fuzzy state series (range code series) is performed according to the following tagging mechanism:(62)PY∀τ∈1,T=1,ifminpττ=1,2,…,T≤pτ≤c1+c22m,ifcm-1+cm2<pτ≤cm+cm+12M,ifcm-1+cm2<pτ<maxpττ=1,2,…,T,m=2,3,…,M-1At last, the monitored series is expressed as a range code time series: (63)P=pττ=1,2,…,T→Zτ,m∀τ∈1,T=zmτ∀τ∈1,T,m=1,2,…,MFor example, Z=2,4,3,1,….,5,7,6. 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]: (64)dzm=αln⁡z-m-ln⁡zmzmdτ+σzmdWwherez-m 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,dW 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 zmτ is performed by using exact discrete-time expression as follows:(65)zmτ=exp⁡ln⁡zmτ-1e-α∆τ+ln⁡z-m-σ22α1-e-α∆τ+N0,1σ1-e-2α∆τ2αwhere∆τ is constant time interval from point τ to τ+1, one month,N0,1 is normally distributed random variable.Calibration of the mean reverting process is based on the following regression:(66)zmτ+1=β0+β1zmτ+εThe speed of mean reverting (α) equals the negative of the slope, long-term equilibrium (z-m) is the intercept divided by the speed, and the rate of range code volatility (σ) is standard deviation from the regression.Letzm=zmτ,τ=0,1,…,T denote an evolution path of range code values with spot values zmτ, where zmτ is calculated by (65). Simulation of (65), for every τ=1,2,…,T-1, is composed of the following two steps: (a) do calculation of (65) and (b) for obtained result determine the range code value. Suppose the zmsimτ is the value obtained by (65) at point τ, then corresponding range code value is assigned as follows:(67)zmτ∀τ∈1,T=1,if0≤zmsimτ≤1+22m,ifym-1+ym2<zmsimτ≤ym+ym+12M,ifzmsimτ>ym-1+ym2,m=2,3,…,M-1Note 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:(68)ZS×Tsim=zmsτS×T→Eq.67YS×Tsim=ymsτS×T→Table1P~S×Tsim=p~msτS×TwhereS 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 20th month, uponS simulations.Figure 3
One simulated evolution path of the range code values on a monthly span.Figure 4
Distribution of range code values obtained by 500 simulations (the20th month).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; I~=a,b,c,a<b<c,b-a<c-b.Operating costs are also very significant source of uncertainties related to the economic value ofTMC. 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]:(69)dcτ=μcτdτ+σdWτwhereμ andσ are the drift and volatility and Wτ is a normalized Brownian motion. Solution of this equation is as follows and it describes behaviour of operating costs:(70)cτ=cτ-1exp⁡μ-σ22∆τ+N0,1σ∆τ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:(71)CS×Tsim=csτS×TBearing 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; d~=a,b,c,a<b<c.Upon discussion on all main components separately, we have to calculate the economic value ofTMC 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:(72)V~conS×Tsim=mcon∙mmr∙P~S×Tsim=mcon∙mmr∙p~msτS×Twherep~msτ is unit selling price of the metal in period τ and simulations 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]:(73)V~hS×Tsim=a~1h∑γ=1Ya2,γh∙MγmγconV~conS×Tsim-CS×Tsim,∀h∈1,HwhereY is the total number of metal concentrates beneficiated from the ore,Mγ is mill recovery rate of the γth metal (%).Space of economic value of mineable block can be represented by the following fuzzy matrix:(74)V~hS×Tsim=v~sτhS×T=a,b,c11ha,b,c12h…a,b,c1τha,b,c21ha,b,c22h…a,b,c2τh⋮⋮⋱⋮a,b,cs1ha,b,cs2h…a,b,cSTh,∀h∈1,HExpected economic values of mineable block trough the planning time are obtained by averaging each column of the matrix V~hS×Tsim. (75)EV~h∀τ∈1,T=1S∑s=1Sa,b,csτh∀τ∈1,T,∀h∈1,HIf 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. (76)EEV~τhtt=1,T=112∑τ=112EV~τhtt=1,T,∀h∈1,HLet V~th,∀h∈1,H,t=1,T be the expected economic value of mineable block through the planning time, then its discounted value for every year is defined as follows: (77)V~thd=V~th1+d~1V~th1+d~2⋯V~th1+d~t,∀h∈1,H,t=1,TFinally, the coefficient v~itd which represents the fuzzy present value of the technological mining cut (discounted value) is calculated by the following equation: (78)v~itd=∑h=1liV~ithd,∀i∈1,N,t=1,T,l1+l2+⋯+lN=HIt means that discounted economic value of theTMC is a sum of discounted economic values of the mineable blocks contained inTMC.


## 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μA of a crisp set A⊆X assigns a value either 0 or 1 to each member inX. This function can be generalized to a function μA~ such that value assigned to the element of the universal setX falls within a specified range, i.e., μA~:X→0,1. The assigned value indicates the membership grade of the element in the setA. The function μA~ is called the membership function and the set A~=x,μA~x:x∈X defined by μA~x for each x∈X is called a fuzzy set [18, 19].A fuzzy numberA~=a,b,c is said to be a triangular fuzzy number if its membership function is given by(1)μA~x=x-ab-a,a≤x≤bx-cb-c,b≤x≤c0,otherwiseFor more details of arithmetic operations on triangular fuzzy numbers, see [16, 18].The absolute value of the triangular fuzzy numberA~=a,b,c is denoted by A~ and defined as follows [16]:(2)A~x=max⁡A~x,-A~x,ifx≥00,ifx<0A ranking function is a function R:FR→R, whereF(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 A~=a,b,c be a triangular fuzzy number then [18](3)RA~=a+2b+c4Linear 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:(4)max∑j=1naj,bj,cjxjsubject to(5)∑j=1nqijxj≤pi,1≤i≤m;1≤j≤n(6)xj∈0,1Since variablesxj and coefficientsqij 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:(7)maxR∑j=1naj,bj,cjxjExample 1.
Let us consider the following fuzzy objective function and convert it by the proposed method:(8)max1,6,9x1+2,3,8x2+4,6,7x3The fuzzy objective function may be written as follows:(9)maxR1x1+2x2+4x3,6x1+3x2+6x3,9x1+8x2+7x3Using arithmetic operations on triangular fuzzy numbers and (3), the fuzzy objective function is converted into the following crisp objective function:(10)max1422x1+16x2+23x3Suppose the solution of this objective function, with respect to a given set of constraints, isx1=x3=1 andx2=0. Fuzzy optimal value of our objective function is obtained by puttingx1,x2, andx3 in (8). The value of the given objective function is(11)1,6,9∙1+2,3,8∙0+4,6,7∙1=5,12,16An 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]:(12)x-0A~=∫acxμA~xdx∫acμA~xdxwhere x-0A~ is the defuzzified value. The defuzzification formula of triangular fuzzy number A~=a,b,c is(13)x-0A~=a+b+c3This formula will be used in this paper. Defuzzified (crisp) value of our objective function is x-05,12,16=11.


## 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:(14)F~V~,X=∑i=1N∑t=1Tv~itdxit=∑i=1N∑t=1Tv~it1+d~txit→maxsubject to(15)∑t=1Txit=1,∀i∈N(16)∑i=1Nxit=1,∀t∈T(17)xit+∑s=1Sxs,t+1≤1,∀i∈N,∀t∈T,s∈S(18)xit∈0,1,∀i∈N,∀t∈T(19)F~-I~≥0where v~itd is fuzzy present value of the technological mining cut, discounted value,v~it is fuzzy economic value of the technological mining cut in timet,xit is binary variable, which equals 1 if and only if the technological mining cuti is mined in timet,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 fromxit, i.e., this set is composed of nonneighbouring mining cuts,d~ is fuzzy discount rate,T is the planning time,I~ 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 fromt tot+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 Figure1. Minerals contained in pillars are nonrecoverable and therefore are not included in the ore reserves of the mine [21].Figure 1
Room and pillar mining method [22].Technological mining cut (TMC) is a part of mineral deposit characterized with respect to the given set of technological requirements (criteria). It means thatTMC 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 ofTMCs having the value of attributes closely to the values of technological requirements. Hence, the first step concerns partition of the deposit in adequate number ofTMCs. 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 suchTMCs (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,B=bhh=1,2,…,H, and each block is characterized by the block attribute vector BAVh=ajhj=1,2,..,A;∀h∈H, whereH is the total number of blocks andA total number of attributes, such as block tonnage, ore grade, etc.A setTMCi=blil<H;∀i∈N is defined as a subset ofB and called technological mining cut.EachTMC is characterized by the mining cut attribute vector MCAVi=aicc=1,2,…,K;∀i∈N, whereK 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); C=Cjj=1,2,…,K.At last, creation of technological mining cuts can be mathematically formulated as a multiobjective partition problem, wherein theTMCs must meet technological performance criteria, subject to the given criteria constraints:(20)max︸a1c∈TMC1fa1c=f1a11,f2a12,…,fKa1Kmax︸a2c∈TMC2fa2c=f1a21,f2a22,…,fKa2K⋮max︸aic∈TMCifaic=f1ai1,f2ai2,…,fKaNKsubject to(21)acl≤ac≤acu,c=1,2,…,Kwherefaic is the ultimate relative closeness (URC) of theith technological mining cut to the positive ideal technological solution, taking into account all technological criteria,acl is the lower bound of value of thecth technological criterion,acu is the upper bound of value of thecth 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:(22)B=⋃i=1NTMCi=TMC1∪TMC2∪…∪TMCiTMCα∩TMCβ=∅,α≠βCreation of the set TMCi=blil<H;∀i∈N 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 vectorMCAVi=aicc=1,2,…,K;∀i∈N and Cteh=cjtehj=1,2,…,K, where Cteh 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 whichTMC is to grow can be concisely expressed by the following decision making matrix:(23)D~=x~ijN×K=TMCCC1C2⋯CjTMC1x~11x~12⋯x~1jTMC2x~21x~22⋯x~2j⋮⋮⋮⋱⋮TMCix~i1x~i2⋯x~ijwhere x~ij is the estimated value of technological mining cutTMCi with respect to the technological criteriaCj. Note that there is a difference between required technological vector and vector of technological criteria, Cteh≠C, 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 ofx~ij with weights (w~j) of technological criteria.(24)P~=p~ijN×K=r~ij∙w~jN×K=x~ij∑i=1Nx~ij∙w~jN×K,∀j∈1,K;∑j=1Kw~j=1To 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:(25)e~j=-1ln⁡N∑i=1Nr~ij∙ln⁡r~ij,∀j∈1,K;0≤e~j≤1The objective weight for each criterion is given by the following equation:(26)w~j=1-e~j∑j=1K1-e~j,∀j∈1,K;0≤w~j≤1The fuzzy positive ideal solution A~+ and negative one A~- is defined as (27)A~+=max︸i=1,2,…,Np~ij∣j∈J,min︸i=1,2,…,Np~ij∣j∈J′=p~1+,p~2+,…,p~K+(28)A~-=min︸i=1,2,…,Np~ij∣j∈J,max︸i=1,2,…,Np~ij∣j∈J′=p~1-,p~2-,…,p~K-(29)max︸i=1,2,…,Np~ij=max︸i=1,2,…,Naij,bij,cij=max︸i=1,2,…,Nbij,∀j∈1,K(30)min︸i=1,2,…,Np~ij=min︸i=1,2,…,Naij,bij,cij=min︸i=1,2,…,Nbij,∀j∈1,KwhereJ = {j = 1,2,...,K∣j associated with criteria that should be maximized};J′ = {j = 1,2,...,K∣j associated with criteria that should be minimized}.The distance from eachTMC to A~+ and A~- is calculated according to the following equations:(31)di+=∑j=1Kdp~ij,p~j+,∀i∈1,N(32)di-=∑j=1Kdp~ij,p~j-,∀i∈1,Nwhere dp~ij,p~j is the distance measurement between two fuzzy triangular numbers, calculated as follows: (33)dp~ij,p~j=13aij-aj2+bij-bj2+cij-cj2The relative closeness coefficient of eachTMC is calculated as(34)RCCi=di-di-+di+,∀i∈1,NDecision on whichTMC is to grow is making according to the following selection rule:(35)TMCgrow=maxRCC1,RCC2,…,RCCNThe second stage problem that considers which block should be added to the selected cluster (TMCgrow) is defined by the following form:(36)badd=maxRCC1add+O1add,RCC2add+O2add,+⋯+,RCCMadd+OMaddwhereRCCMadd is the new relative closeness coefficient of the newTMC obtained after adding themth neighbouring block to theTMCgrow and neighbouring block is the block that has at least one common edge withTMCgrow,OMadd 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:(37)D~block=φ~mjM×K=TMCgrow∪bCC1C2⋯CjTMCgrow∪b1φ~11φ~12⋯φ~1jTMCgrow∪b2φ~21φ~22⋯φ~2j⋮⋮⋮⋱⋮TMCgrow∪bmφ~m1φ~m2⋯φ~mjwhereφ~mj is the new estimated value of the TMCgrow∪bm; m=1,2,..,M with respect to the technological criteriaCj.The valueφ~mj is estimated after making the union of the attribute vector of the TMCgrow and the attribute vector of the neighbouring block:(38)φ~mj=MCAVgrow∪BAVm=acgrowc=1,2,…,K∪ajmj=1,2,…,A;∀m∈M;K=AWith 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]:(39)OMadd=Bbef-BmaftBbef=Bmcom+Bmngh+Bmwaste-3Bbef,∀m∈1,MwhereBbef is number of mineable blocks surrounding theTMCgrow before adding the newmth block,Bmaft is number of mineable blocks surrounding theTMCgrow after adding the newmth block.Bmcom is number of common edges betweenTMCgrow and block to be added,Bmngh is number of common edges between block to be added and remainingTMCs,Bmwaste 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 setBbef and Bmaft it is necessary to meet the following two spatial constraints:(i)
only blocks having at least one common edge withTMCgrow can be added to theTMCgrow,(ii)
TMCi,∀i∈1,N, 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 anyTMC cannot be element of Bbef and Bmaft, respectively. Suppose theTMC1 is selected to grow; see Figure 2.Figure 2
Spatial plan ofTMCs.According to spatial constraints only blocks 20, 32, and 33 can be elements of the setBbef because block 26 violates spatial homogeneity of theTMC2. If we add block 33 to theTMC1 than only blocks 20, 32, and 34 can be elements of the set Bmaft because block 46 violates spatial homogeneity of theTMC3. Hence, for TMC1∪b33 value of the penalty function is equal to 0. Note, mutual takeover of blocks byTMCs is allowed, but spatial homogeneity ofTMCs must be always preserved.In our model, the block attribute vectorBAVh=ajhj=1,2,..,A;∀h∈H is composed of the following components:(40)BAVh=a~1h,a2,1h,a2,2h,…,a2,γhj=1+γ;∀h∈Hwherea~1h is ore tonnage in blockh (t), expressed as fuzzy triangular number,a2,γh is grade of thehth block with respect to the γth metal (%),γ is total number of metal concentrates beneficiated from the ore. For polymetallic ore, γ>1.The mining cut attribute vectorMCAVi=aicc=1,2,…,K;∀i∈N is composed of the following components:(41)MCAVi=aicc=1,2,…,K;∀i∈N=a~i1,ai,2,ai,3,1,ai,3,2,…,ai,3,γc=2+γ;∀i∈Nwherea~i1 is ore tonnage in theTMC (t), expressed as fuzzy triangular number,ai2 is compactness of theTMC,ai,3,γ is standard deviation of the grade in theTMC, with respect to the γth metal (%).Ore tonnage in theTMC represents a total sum of ore tonnage in blocks contained withinTMC:(42)a~i1=∑h=1lia~1h,∀i∈1,N,l1+l2+⋯+lN=HCompactness of theTMC is defined as the following ratio: of the square of the perimeter and the area of theTMC,(43)ai2=Pi2Ai=Pi2li∙be2,∀i∈1,NwherePi is perimeter of theith TMC,Ai is area of theith TMC,li is total number of blocks in theith TMC,be is length of the block edge (m).Standard deviation of the grade in theTMC, with respect to the γth metal is calculated as follows:(44)ai,3,γ=1li∑h=1lia-2,γ-a2,γh2,li<H,∀i∈1,N,∀γ∈1,YThe required technological vector Cteh=cjtehj=1,2,…,K includes the following components: (45)Cteh=cjtehj=1,2,…,K=c~1teh,c2teh,c3,1teh,c3,2teh,…,c3,γtehj=2+γwherec~1teh is annual capacity of production (t/year), expressed as fuzzy triangular number,c2teh is desired or target value of compactness ofTMC, and it is set up to 16,c3,γteh 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):(46)c~1teh=1N∑h=1Ha~1hTarget value of compactness of theTMC is expressed by the Schwartzberg's index of the simple square geometric shape [33]: (47)c2teh=4e2e2=16wheree 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 theTMC, with respect to the γth metal. It is calculated as follows and there is no target value for this component:(48)c3,γteh=ai,3,γThe main aim of the vector of technological criteria is to enable creation ofTMCs so that the ultimate similarity between vector MCAVi=aicc=1,2,…,K;∀i∈N and Cteh=cjtehj=1,2,…,K is maximized. It means that eachTMC must meet technological requirements as maximum as possible. Vector of technological criteria C=Cjj=1,2,…,K is composed of the following components:(49)C=Cjj=1,2,…,K=C~1,C2,C3,1,C3,2,…,C3,γj=2+γwhere C~1 is absolute distance between annual capacity of production and ore tonnage in theTMC (t),C2 is absolute distance between target value of compactness and compactness ofTMC,C3,γ is standard deviation of the grade in theTMC with respect to the γth metal (%).Value of technological mining cutTMCi with respect to the technological criteria C~1 is calculated as follows, and it should be minimized:(50)x~i1=a-0a~i1c-0c~1teh×100%,∀i∈1,NValue of technological mining cutTMCi with respect to the technological criteria C2 is calculated as follows, and it should be minimized:(51)xi2=16-ai2,∀i∈1,NValue of technological mining cutTMCi with respect to the technological criteria C3,γ is as follows, and it should be minimized:(52)xi,3,γ=ai,3,γ,∀i∈1,N,∀γ∈1,YThe same calculations are applied when we define values of the decision making matrix D~block=φ~mjM×K.Set of the technological criteria constraintsacl≤ac≤acu,c=1,2,…,K is composed of only constraint related to the criterion C~1: (53)-10%≤a1≤+10%-10%≤1-a-0a~i1c-0c~1teh×100%≤+10%;∀i∈1,Nwhere a-0a~i1 is defuzzified value of the ore tonnage in theTMC and c-0c~1teh 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 graphMD=(B,E), where B=bhh=1,2,…,H is a set of mineable blocks, and E=euv∣u,v=1,2,..,H;u≠v is a set of edges witheuv representing the common edge between blocku and blockv [34, 35]. The problem of creation ofTMCs is solved by the multicriteria partitioning of graphMD. This approach attempts to address situation in which the creation ofTMCs should simultaneously maximize the ultimate relative closeness (URC) of the eachTMC to the technological requirements, with respect to technological criteria and constraints. Given a set B=bhh=1,2,…,H, and total number ofTMCs, N≥2, then the model of creation ofTMCs can be formulated as follows:(54)Z=max︸∀i∈1,NURCi=max︸∀i∈1,Ndi-⋃h=1li<Hbihdi-⋃h=1li<Hbih+di+⋃h=1li<Hbihsubject to(55)∑i=1Nbih=1,∀h∈1,H(56)a~1l≤1-defuzzified∑h=1Ha~1ihbihdefuzzifiedc~1teh×100%≤a~1u,∀i∈1,N(57)bih∈0,1,∀h∈H,∀i∈Nwhere ⋃h=1li<Hbih represents the union of mineable blocks contained in the ultimateithTMC. Constraint (55) forces that a mineable block can belong to one and only oneTMC. It means no intersection betweenTMCs, TMCα∩TMCβ=∅,α≠β. Constraint (56) is already explained, and it does not allow anyTMC to be terminated, TMCi≠∅,∀i∈1,N. Constraint (57) is related to binary variable bih, which equals 1 if and only if the mineable blockh belongs to theithTMC.Here we present the algorithm of partitioning the graphMD, and it is based on the constrained polygonal spatial clustering algorithm [24–26]. Code of the algorithm is represented by the Algorithm 1.Algorithm 1:Algorithm of the TMC creation.
1. Set iteration toξ=1; 2. Create a set;TMCi=TMC1,TMC2,…,TMCN 3. Create a set;C=Cjj=1,2,…,K=C~1,C2,C3,1,C3,2,…,C3,γj=2+γ 4. Select the best technological mining cut (BTMC) to grow; B T M C = T M C g r o w = m a x R C C 1 , R C C 2 , … , R C C N 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 BMB∉NMB if: for ξ+1, BMB∈TMCαgrow for ξ+2, BMB∈TMCβgrow for ξ+3, bh∈TMCαgrow This is to avoid algorithm gets stuck. 6. Select the best mineable blocks (BMB) from NMB to add to BTMC; B M B = m a x R C C 1 a d d + O 1 a d d , R C C 2 a d d + O 2 a d d , + ⋯ + , R C C M a d d + O M a d d 7. AddBMB to BTMC and create new technological mining cut; N T M C = B T M C ∪ B M B 8. Update the state of the mining cut attribute vector of theNTMC 9. If: 9.1 no free mineable blocks 9.2 -10%≤1-a-0a~i1c-0c~1teh×100%≤+10%;∀i∈1,N 9.3 ai,3,γ,∀i∈1,N,∀γ∈1,Y → converge than Stop 10. Else go to Step1 and set ξ+1Our 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 ofN 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 initialTMC). Obviously, the initialTMC is composed of one and only one mineable block. From the set of peripheral blocks we select the firstN blocks according to the decreasing order of values of relative closeness coefficient. Once the initialTMCs are selected, they begin to be alive and the process ofTMCs growing can run. EachTMC is grown by adding neighbouring blocks to him one by one until the desired state ofTMC is achieved. Growing ofTMCs is an iterative process in nature. At the beginning of the each iteration we first select whichTMC is to be grown. Selection is based on the measuring the multidimensional distance between current and desired state of the eachTMC. It is expressed by the relative closeness coefficient andTMC having the largest value ofRCC 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 bestTMC but with adding the penalty function. Mutual takeover of the block between two neighbouringTMCc is allowed but homogeneity of eachTMC must always be preserved. During the process of growing, the infinity mutual takeover can arise. If the block has mutually been acquired between twoTMCs 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 MCAVi=aicc=1,2,…,K;∀i∈N 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 theTMC is based on the three following main components: metal price, costs, and discount rate.One of the most influencing variables on the economic value ofTMC 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 variableX. If the value of that variable is governed by the laws of probability, then variableX can be treated as a stochastic variable. Suppose that we monitored values ofX at equal interval, ∆τi=τi-τi-1=const.; i=2,3,…,T. 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, Xτ=x1,x2,…,xT, 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 dividingN objects into a specified number ofM disjoint groups that are also called classes or clusters. Fuzzy C-mean algorithm is based on minimization of the following least-squared errors function: (58)F=U,X,C=∑τ=1T∑m=1Muτmωxτ-cm2subject to(59)∑m=1Muτm=1,τ=1,2,…,T(60)0≤uτm≤1,m=1,2,…,M,τ=1,2,…,T(61)0<∑τ=1Tuim<m,m=1,2,…,MwhereX is the vector of monitored metal prices; X=x1,x2,…,xT⫅RT,C is the vector of class centers; C=c1,c2,…,cM⫅RM,U is the fuzzy partition matrix; U=uτmT×M,ω is the coefficient of fuzzification and we take value equal to 2.The objective functionF is iteratively minimized. The iteration process stops until Fj+1-Fj<δ, wherej is the number of iteration and δ represents the minimum amount of improvement. For more details see [36–39]. Suppose that partitioningF(X,C) has been done and the sequence of obtained centers is sorted in an ascending order; c1<c2<⋯<cM.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 tripletA~=a,b,c, where parametersa,b, andc, 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: (1) a≤b≤c; (2) membership function is increasing in the interval a,b and decreasing in the interval b,c. Range boundaries and value of code of linguistic variables are shown in Table 1.Table 1
The transformation of linguistic variables to triangular fuzzy numbers and range code.
Linguistic Triangular fuzzy number Membership Range code value variable function y m , m = 1,2 , … , 7 VVL a < c 1 , b = c 1 , c 1 < c ≤ c 1 + c 2 2 0 ; 1 ; 0 y 1 = 1  VL c 1 + c 2 2 < a , b = c 2 , c 2 < c ≤ c 2 + c 3 2 0 ; 1 ; 0 y 2 = 2  L c 2 + c 3 2 < a , b = c 3 , c 3 < c ≤ c 3 + c 4 2 0 ; 1 ; 0 y 3 = 3  M c 3 + c 4 2 < a , b = c 4 , c 4 < c ≤ c 4 + c 5 2 0 ; 1 ; 0 y 4 = 4  H c 4 + c 5 2 < a , b = c 5 , c 5 < c ≤ c 5 + c 6 2 0 ; 1 ; 0 y 5 = 5  VH c 5 + c 6 2 < a , b = c 6 , c 6 < c ≤ c 6 + c 7 2 0 ; 1 ; 0 y 6 = 6  VVH c 6 + c 7 2 < a , b = c 7 , c 7 < c 0 ; 1 ; 0 y 7 = 7Suppose a monthly monitored metal price time series is givenP=pττ=1,2,…,T and a range code vector Y=ymm=1,2,…,M is defined on it. Transformation of monitored series to the fuzzy state series (range code series) is performed according to the following tagging mechanism:(62)PY∀τ∈1,T=1,ifminpττ=1,2,…,T≤pτ≤c1+c22m,ifcm-1+cm2<pτ≤cm+cm+12M,ifcm-1+cm2<pτ<maxpττ=1,2,…,T,m=2,3,…,M-1At last, the monitored series is expressed as a range code time series: (63)P=pττ=1,2,…,T→Zτ,m∀τ∈1,T=zmτ∀τ∈1,T,m=1,2,…,MFor example, Z=2,4,3,1,….,5,7,6. 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]: (64)dzm=αln⁡z-m-ln⁡zmzmdτ+σzmdWwherez-m 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,dW 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 zmτ is performed by using exact discrete-time expression as follows:(65)zmτ=exp⁡ln⁡zmτ-1e-α∆τ+ln⁡z-m-σ22α1-e-α∆τ+N0,1σ1-e-2α∆τ2αwhere∆τ is constant time interval from point τ to τ+1, one month,N0,1 is normally distributed random variable.Calibration of the mean reverting process is based on the following regression:(66)zmτ+1=β0+β1zmτ+εThe speed of mean reverting (α) equals the negative of the slope, long-term equilibrium (z-m) is the intercept divided by the speed, and the rate of range code volatility (σ) is standard deviation from the regression.Letzm=zmτ,τ=0,1,…,T denote an evolution path of range code values with spot values zmτ, where zmτ is calculated by (65). Simulation of (65), for every τ=1,2,…,T-1, is composed of the following two steps: (a) do calculation of (65) and (b) for obtained result determine the range code value. Suppose the zmsimτ is the value obtained by (65) at point τ, then corresponding range code value is assigned as follows:(67)zmτ∀τ∈1,T=1,if0≤zmsimτ≤1+22m,ifym-1+ym2<zmsimτ≤ym+ym+12M,ifzmsimτ>ym-1+ym2,m=2,3,…,M-1Note 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:(68)ZS×Tsim=zmsτS×T→Eq.67YS×Tsim=ymsτS×T→Table1P~S×Tsim=p~msτS×TwhereS 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 20th month, uponS simulations.Figure 3
One simulated evolution path of the range code values on a monthly span.Figure 4
Distribution of range code values obtained by 500 simulations (the20th month).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; I~=a,b,c,a<b<c,b-a<c-b.Operating costs are also very significant source of uncertainties related to the economic value ofTMC. 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]:(69)dcτ=μcτdτ+σdWτwhereμ andσ are the drift and volatility and Wτ is a normalized Brownian motion. Solution of this equation is as follows and it describes behaviour of operating costs:(70)cτ=cτ-1exp⁡μ-σ22∆τ+N0,1σ∆τ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:(71)CS×Tsim=csτS×TBearing 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; d~=a,b,c,a<b<c.Upon discussion on all main components separately, we have to calculate the economic value ofTMC 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:(72)V~conS×Tsim=mcon∙mmr∙P~S×Tsim=mcon∙mmr∙p~msτS×Twherep~msτ is unit selling price of the metal in period τ and simulations 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]:(73)V~hS×Tsim=a~1h∑γ=1Ya2,γh∙MγmγconV~conS×Tsim-CS×Tsim,∀h∈1,HwhereY is the total number of metal concentrates beneficiated from the ore,Mγ is mill recovery rate of the γth metal (%).Space of economic value of mineable block can be represented by the following fuzzy matrix:(74)V~hS×Tsim=v~sτhS×T=a,b,c11ha,b,c12h…a,b,c1τha,b,c21ha,b,c22h…a,b,c2τh⋮⋮⋱⋮a,b,cs1ha,b,cs2h…a,b,cSTh,∀h∈1,HExpected economic values of mineable block trough the planning time are obtained by averaging each column of the matrix V~hS×Tsim. (75)EV~h∀τ∈1,T=1S∑s=1Sa,b,csτh∀τ∈1,T,∀h∈1,HIf 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. (76)EEV~τhtt=1,T=112∑τ=112EV~τhtt=1,T,∀h∈1,HLet V~th,∀h∈1,H,t=1,T be the expected economic value of mineable block through the planning time, then its discounted value for every year is defined as follows: (77)V~thd=V~th1+d~1V~th1+d~2⋯V~th1+d~t,∀h∈1,H,t=1,TFinally, the coefficient v~itd which represents the fuzzy present value of the technological mining cut (discounted value) is calculated by the following equation: (78)v~itd=∑h=1liV~ithd,∀i∈1,N,t=1,T,l1+l2+⋯+lN=HIt means that discounted economic value of theTMC is a sum of discounted economic values of the mineable blocks contained inTMC.


## 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 Figures5 and 6 and Tables 2 and 3.Table 2
Input parameters.
Parameter Value Number of blocks 115  Block dimension 25×25×5 m  Pillar dimension 8×8×5 m  Planning time 5 years; 60 months  Number of technological mining cuts 5  Metal content of concentrate  Lead 70%  Zinc 50%  Metal recovery rate  Lead 95 %  Zinc m c o n % - 8 ∙ 100 m c o n % ≤ 85 % ; m c o n ∙ m r → m c o n % - 8 m c o n % - 8 ∙ 100 m c o n % > 85 % ; m c o n ∙ m r → 85 % 42 %  Mill recovery rate  Lead 92 %  Zinc 80 %  Capital costs ($) 1 2  0 0 0  0 00 1 3  0 0 0  0 00 1 5  0 0 0  0 00  Capacity of production (t/year)-triangular fuzzy number 15 5  2 04 17 2  4 49 18 9  6 94  Capacity constraint (t/year) ± 10%  Fuzzy C-mean clustering-coefficient of fuzzification 2  Metal price linguistic variablesM=7 very very low (VVL), very low (VL), low (L), medium (M), high (H), very high (VH), very very high (VVH)  Discount rate-triangular fuzzy number (%) 1 0 11 1 2Table 3
Reserves of mineable blocks.
Block Reserves (t) Block Reserves (t) triangular fuzzy number triangular fuzzy number 1 5301 5891 6480 59 6458 7176 7893 2 6844 7604 8365 60 6954 7727 8499 3 5081 5646 6210 61 6293 6992 7691 4 7064 7849 8634 62 7119 7910 8701 5 7615 8461 9307 63 6128 6809 7489 6 6954 7727 8499 64 7284 8094 8903 7 7725 8583 9442 65 5962 6625 7287 8 5715 6350 6984 66 8977 9974 10972 9 6747 7497 8247 67 6128 6809 7489 10 5508 6120 6732 68 7119 7910 8701 11 7367 8186 9004 69 6293 6992 7691 12 6871 7635 8398 70 6954 7727 8499 13 7491 8323 9156 71 6458 7176 7893 14 6747 7497 8247 72 6789 7543 8297 15 7615 8461 9307 73 6623 7359 8095 16 6623 7359 8095 74 6623 7359 8095 17 7615 8461 9307 75 6789 7543 8297 18 5907 6564 7220 76 6458 7176 7893 19 6871 7635 8398 77 6954 7727 8499 20 5715 6350 6984 78 6293 6992 7691 21 7064 7849 8634 79 7119 7910 8701 22 5522 6135 6749 80 6128 6809 7489 23 6816 7574 8331 81 7284 8094 8903 24 6678 7421 8163 82 9130 10144 11158 25 7367 8186 9004 83 6334 7038 7742 26 6541 7268 7994 84 7243 8048 8853 27 7505 8339 9172 85 6486 7206 7927 28 6458 7176 7893 86 7092 7880 8667 29 6637 7375 8112 87 6637 7375 8112 30 6279 6977 7674 88 6940 7711 8482 31 6816 7574 8331 89 6789 7543 8297 32 6100 6778 7456 90 6789 7543 8297 33 6995 7772 8550 91 6940 7711 8482 34 5921 6579 7237 92 6637 7375 8112 35 7174 7971 8768 93 7092 7880 8667 36 6486 7206 7927 94 6486 7206 7927 37 7092 7880 8667 95 7243 8048 8853 38 6637 7375 8112 96 8895 9883 10872 39 6940 7711 8482 97 5921 6579 7237 40 6789 7543 8297 98 6995 7772 8550 41 6789 7543 8297 99 6100 6778 7456 42 6940 7711 8482 100 6816 7574 8331 43 6637 7375 8112 101 6279 6977 7674 44 7092 7880 8667 102 6637 7375 8112 45 6486 7206 7927 103 6458 7176 7893 46 7243 8048 8853 104 6458 7176 7893 47 6334 7038 7742 105 6637 7375 8112 48 7394 8216 9038 106 6279 6977 7674 49 6183 6870 7557 107 6816 7574 8331 50 6128 6809 7489 108 6100 6778 7456 51 7119 7910 8701 109 6678 7421 8163 52 6293 6992 7691 110 7229 8033 8836 53 6954 7727 8499 111 6816 7574 8331 54 6458 7176 7893 112 6954 7727 8499 55 6789 7543 8297 113 6954 7727 8499 56 6623 7359 8095 114 7092 7880 8667 57 6623 7359 8095 115 6816 7574 8331 58 6789 7543 8297 ∗Bold italic font indicates peripheral blocks.Figure 5
Ore body with regular pattern of square pillars.Figure 6
Pb and Zn grade distribution (%).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, initialTMCs (seeds) were defined; see Figure 7.Figure 7
InitialTMCs (seeds).Having selected initialTMCs algorithm proceeded to do iterations and characteristically phases of transition are presented by Figure 8.Figure 8
Results of the algorithm after: (a) 50 iterations. (b) 100 iterations. (c) 150 iterations. (d) 200 iterations.
(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 theTMC (see Figure 9).Figure 9
Final partition of the ore body.Stopping condition 9.2 related to the standard deviation of lead and zinc in theTMC is represented by Figure 10.Figure 10
Convergence of standard deviation of lead and zinc metal in theTMC.Realized technological mining cut attribute values that consider ore tonnage and compactness of theTMC are shown in Table 4.Table 4
Ore tonnage and compactness ofTMC.
TMC Ore tonnage (t) Compactness Production capacity error (%)  1 152640 39.13 -1.65 169601 186561  2 153217 30.72 -1.28 170241 187266  3 154940 29.39 -0.17 172156 189371  4 155573 62.78 +0.24 172859 190145  5 159649 48.16 +2.86 177388 195127 ∗(a,b,c) triangular fuzzy number is expressed as column.Mean absolute percentage error of the production capacity is(79)1.65+1.28+0.17+0.24+2.865=1.24%and 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 ofTMCs achieved by algorithm are suitable. Algorithm succeeded to create threeTMCs 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.Table 5
Ore grade characteristics ofTMC.
TMC Mean value of the grade (%) Standard deviation of the grade (%) Coefficient of the grade variation (%) Lead Zinc Lead Zinc Lead Zinc 1 2.263 3.364 0.325 1.198 14.392 35.607 2 1.801 3.832 0.292 0.418 16.222 10.931 3 2.170 3.937 0.381 0.820 17.565 20.823 4 1.950 2.573 0.169 0.928 8.710 36.079 5 2.136 4.636 0.267 1.262 12.531 27.226Grade 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 inTMC1 andTMC4 coefficient of zinc grade variation is greater than 32.7% but not so significantly. Realized distribution of lead and zinc grade inTMCs enables mineral processing planners to create much more flexible plans because they are informed in advance about expected grade distribution.Since the all attributes ofTMCs 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.Table 6
Class centers of metal price.
Metal price Class center (m=1,2,…,7) c 1 c 2 c 3 c 4 c 5 c 6 c 7 Lead ($/t) 1176 1689 1863 2064 2205 2381 2665 Zinc ($/t) 1236 1816 2014 2176 2344 2695 3244Figure 11
Monitored price on monthly span: (a) zinc; (b) lead.
(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 Table7 and Figure 12.Table 7
The transformation of linguistic variables.
Linguistic variable Triangular fuzzy number of metal price ($/t) Membership function Range code valueym,m=1,2,…,7 Lead Zinc VVL (963 1176 1433) (1101 1236 1526) 0 ; 1 ; 0 y 1 = 1 VL (1433 1689 1776) (1526 1816 1915) 0 ; 1 ; 0 y 2 = 2 L (1776 1863 1963) (1915 2014 2095) 0 ; 1 ; 0 y 3 = 3 M (1963 2064 2134) (2095 2176 2260) 0 ; 1 ; 0 y 4 = 4 H (2134 2205 2293) (2260 2344 2519) 0 ; 1 ; 0 y 5 = 5 VH (2293 2381 2523) (2519 2695 2970) 0 ; 1 ; 0 y 6 = 6 VVH (2523 2665 3080) (2970 3244 3533) 0 ; 1 ; 0 y 7 = 7Figure 12
Linguistic variables expressed by triangular fuzzy numbers: (a) lead; (b) zinc.
(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.Figure 13
Monitored series expressed by fuzzy states and range codes: (a) lead price as a fuzzy state, (b) lead price as a range code, (c) zinc price as a fuzzy state, and (d) zinc price as a range code.
(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 Table8.Table 8
Parameters of simulation.
Process Value Lead metal pricemean reverting process-monthly span spot value (τ=0)2361 $/t→range code 6→(2293 2381 2523) $/tequilibriumrange code 4→(1963 2064 2134) $/tspeed of mean reversion 0.11251volatility 0.18491  Zinc metal pricemean reverting process-monthly span spot value (τ=0)3060 $/t→range code 7→(2970 3244 3533)equilibriumrange code 4→(2095 2176 2260)speed of mean reversion 0.06133volatility 0.16416  Operating costsItô-Doob stochastic process - monthly span spot value (τ=0) 50 $/tdrift 0,0025volatility 0,0125  Time span of simulation (τ=1,2,…,60), T=60 months  Number of simulations 500One simulation of lead and zinc metal price, unit operating cost and resulting economic value of mineable block 1, is represented in Figure14.Figure 14
One possible artificial scenario of economic value of block 1.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 Figure15.Figure 15
Expected economic values of block 1.Transformation of forecasted economic values of block 1, from monthly to annual level, is represented by Figure16.Figure 16
Forecasted economic values of block 1-annual span.Discounted or present economic values of mineable block 1 are represented by Figure17.Figure 17
Present economic values of block 1.Applying the process of simulation on all blocks simultaneously, we obtain their present economic values and we can calculate the present economic values ofTMCs. According to the results of mineral deposit partition, following present economic values of every TMC are shown in Table 9.Table 9
Present economic values ofTMC.
TMC Value ($) Year 1 Year 2 Year 3 Year 4 Year 5 1 5104734 3944737 3044364 2443462 1924433 6591226 5356427 4305812 3576349 2926394 8228167 6904516 5704581 4871250 4096627  2 5583555 4336349 3366286 2691067 2117861 7190095 5863841 4730052 3916424 3201541 8951887 7536746 6246798 5320340 4472001  3 6213929 4891545 3857917 3139633 2524390 7923489 6519182 5317593 4456306 3694805 9796030 8298845 6935676 5958908 5056129  4 2976598 2108873 1450337 1038819 686021 4063918 3133980 2351322 1834430 1374755 5275030 4267887 3358689 2754948 2200817  5 6198802 4899506 3882193 3166488 2555475 7882521 6503582 5321728 4466255 3712085 9723915 8256167 6917828 5949400 5056656 ∗(a,b,c) triangular fuzzy number is expressed as column.By this way, every coefficientv~itd,i=1,2,..,5;t=1,2,…,5, of the objective function (see (14)) is defined. Solution of the fuzzy objective function isx14=1;x22=1;x31=1;x45=1;x53=1 while the remains are equal to zero. It means that optimal plan of production isTMC3 (year 1) →TMC2 (year 2) →TMC5 (year 3) →TMC1 (year 4) →TMC4 (year 5) (see Figure 18).Figure 18
Optimal plan of production.Value of the objective function isF~=175619552406016531322673$. Value of (19) is(80)F~-I~=175619552406016531322673-120000001300000015000000=25619551106016619322673$This 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 theTMCs (clusters) of the ore body. Two-stage fuzzy multicriteria clustering algorithm as iterative process is applied for creation ofTMCs where each cluster meets the following technological requirements: annual capacity of production, compactness of the shape ofTMC, and standard deviation of ore grade in theTMC.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 eachTMCs 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 ofTMCs, 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.



---

*Source: 3078234-2019-06-27.xml*