Integrating,geometallurgical,ball,mill,throughput,predictions,into,short-term,stochastic,production,scheduling,in,mining,complexes

Christian Both*,Roussos Dimitrakopoulos

COSMO-Stochastic Mine Planning Laboratory,Department of Mining and Materials Engineering,McGill University,Montreal H3A 0E8,Canada

Keywords: Geometallurgy Stochastic optimization Short-term open pit mine production scheduling Measurement while drilling Non-additivity Hardness

ABSTRACT This article presents a novel approach to integrate a throughput prediction model for the ball mill into short-term stochastic production scheduling in mining complexes.The datasets for the throughput prediction model include penetration rates from blast hole drilling(measurement while drilling),geological domains,material types,rock density,and throughput rates of the operating mill,offering an accessible and cost-effective method compared to other geometallurgical programs.First,the comminution behavior of the orebody was geostatistically simulated by building additive hardness proportions from penetration rates.A regression model was constructed to predict throughput rates as a function of blended rock properties,which are informed by a material tracking approach in the mining complex.Finally,the throughput prediction model was integrated into a stochastic optimization model for short-term production scheduling.This way,common shortfalls of existing geometallurgical throughput prediction models,that typically ignore the non-additive nature of hardness and are not designed to interact with mine production scheduling,are overcome.A case study at the Tropicana Mining Complex shows that throughput can be predicted with an error less than 30 t/h and a correlation coefficient of up to 0.8.By integrating the prediction model and new stochastic components into optimization,the production schedule achieves weekly planned production reliably because scheduled materials match with the predicted performance of the mill.Comparisons to optimization using conventional mill tonnage constraints reveal that expected production shortfalls of up to 7% per period can be mitigated this way.

Short-term mine production planning aims to make daily,weekly,or monthly operational decisions that best meet strategic production targets under existing operating conditions and constraints [1].Blom et al.[2] reviewed past advancements in shortterm planning for open pit mines,while recent developments consider the simultaneous stochastic optimization of a short-term mine production schedule and fleet management decisions in mining complexes under supply and equipment performance uncertainty [3,4].However,the incorporation of pertinent geometallurgical properties into short-term mine production scheduling has not been addressed sufficiently.

Geometallurgy describes the relationships between rock characteristics extracted from mineral deposits and its metallurgical/processing behavior further downstream in a mineral value chain.Key geometallurgical properties should thus be integrated into short-term mine production scheduling,which would lead to better-informed mine plans[5,6].Important properties to consider are the rock hardness and grindability,which largely control the throughput in the comminution circuit of the processing plant(s).The integration of a throughput prediction model is crucial for achieving short-term production targets since plant throughput directly influences metal production.Having a detailed modelling of the comminution circuit in production scheduling is also beneficial from a cost and environmental perspective since a substantial portion of a mine’s energy consumption can be attributed to the comminution of the ore in hard rock open pit mines [7].Although many geometallurgical models that predict throughput/comminution performance of the mineral reserve based on ore hardness and grindability have been developed [8-11],they are typically not designed to interact with mine production scheduling.These geometallurgical throughput models typically rely on very sparse measurements of hardness and typically ignore the non-additive nature of hardness when changing the support scale from points(measurement scale) to blocks (mining scale),as well as when materials are blended from various sources within a mining complex.It is important to note that blends of additive rock properties(e.g.,metal grades) can be correctly calculated by their (simple or weighted) averages,whereas blends of non-additive properties behave differently than their calculated average.Prediction models that rely on non-additive hardness properties can thus only be accurate in mine settings where no or very limited blending occurs.

Comminution performance of the mineral reserve is conventionally assessed by different hardness and grindability indices.Lynch et al.[12] reviewed the mostly used indices to date,including Bond’s Ball Mill Work Index (Wi),the SAG Power Index (SPI),the Drop Weight Index (DWi) and the resistance to impact breakage (A × b).To date,the limited integration of throughput/comminution performance into mine production scheduling involves a spatial model of the mentioned hardness and grindability indices[6,13-15].The spatial model of geometallurgical properties is built using geostatistical techniques and is commonly referred to as the geometallurgical model.This endeavor is challenging due to the non-additive nature of most hardness and grindability indices[16-18].Specific care has to be taken for spatial interpolation[19,20],change of support,including drill-hole compositing and upscaling to mineable units [11,21],and blending of materials in downstream processes [22].The latter directly contradicts the common assumptions of mine production scheduling models that use mixed integer linear programming.These models need to assume the additive behavior of variables to enforce blending constraints and processing capacities.The other major challenge is that information about hardness properties is typically derived from very sparse geometallurgical sampling [13,23].Part of the reason for sparse sampling is that drilling campaigns and laboratory tests on the samples’ comminution behavior are typically costly and time-consuming [24].

Several proxy models have been tested for predicting hardness variables as a function of primary rock attributes,which is also termed a primary-response framework [25].Applications include multiple linear regression(MLR)[26,27],projection pursuit regression [28],and a variety of supervised learning techniques [29].However,these proxy models are typically constructed using samples in drill-core support and thus do not guarantee the correct mapping of larger volumes (selective mining unit,SMU) [30].Furthermore,the calibration of these proxy models is typically performed on an extremely limited number of hardness samples.

In contrast to the above,the use of penetration rates stemming from measurement while drilling (MWD) data is investigated in this article to inform the mineral reserve about hardness and comminution performance.The potential of these spatially dense datasets has already been exploited in mining for enhanced blast design[31],the detection of rock fractures[32,33],rock characterization in various commodities [34-36],and other applications.Although suggested in the technical literature [24,31],MWD has not yet been utilized to create a direct link between the mineral reserve and its comminution performance in milling and grinding circuits.Such a link is established in the novel approach presented herein,where the rate of penetration (ROP) serves as one of the main spatial indicators for ball mill throughput.The utilization of ROP is motivated by the ability to indicate rock-type,strength,and alteration [34,37].

The next step in a conventional workflow of integrating throughput into mine production scheduling involves empirically derived process models,which are used to calculate the expected mill throughput rates using as input the geometallurgical block model of hardness and grindability indices.Examples of frequently used process models to date are the Bond equation [38],the SMC model[39],and the comminution economic evaluation tool(CEET)-[40].The mill throughput rates are either calculated per resource block or per geometallurgical domain [9,14].Morales et al.[15]presented a stochastic formulation for long-term single mine production scheduling accounting for uncertain throughput rates through geostatistical simulations.A linear constraint ensures that processed blocks do not exceed the available mill hours per production year.The above efforts towards integrating throughput into mine production scheduling are limited in that the nonadditive upscaling behavior of hardness is ignored,and the resulting single-mine production schedule is informed by assessing the mill throughput and economic value per mining block independently.In this way,the schedule ignores that materials are blended and processed together with other materials,which is typical in multi-pit,multi-processor mining complexes.As a result,the non-additive comminution behavior of blended materials in the short-term cannot be captured,and the components of the mining complex are not optimized simultaneously.

As an advantage over the conventional mixed integer linear programming models for production scheduling,the simultaneous optimization of mining complexes successfully takes non-linear blending processes in stockpiles and processing facilities into account to maximize net present value (NPV) [41,42].Initially developed for long-term planning[43,44],simultaneous optimization of mining complexes capitalizes on synergies between decisions that are conventionally optimized separately.Kumar and Dimitrakopoulos [45] introduced geometallurgical constraints related to hardness in the simultaneous optimization of mining complexes for long-term planning.The constraints aim to achieve a consistent throughput by sending pre-defined ratios of hard and soft rock to the processing plant(s).The mineral reserve is categorized into hard and soft rock by stochastically simulated Wi and SPI values,where the hardness ratios are chosen arbitrarily and do not guarantee optimal material selection for throughput maximization.

In the present article,an alternative approach is proposed which integrates geometallurgical hardness properties into shortterm mine production scheduling in mining complexes.More specifically,a ball mill throughput prediction model informed by blended rock attributes related to hardness is constructed,which is then integrated into the simultaneous stochastic optimization of mining complexes.The utilized datasets for throughput prediction include penetration rates from blast hole drilling and measured throughput rates of the operating ball mill,supplemented by geological domains,material types,and rock density.Furthermore,truck cycle data is used for material tracking.Given that the penetration rates informing rock hardness and the measured throughput rates of the ball mill stem from production data,the novel approach gives a potential financial advantage compared to laboratory grinding tests for throughput prediction.Additionally,the metallurgical response of the processed rock,i.e.,throughput,is measured directly in the appropriate support scale,and nonadditivity issues of hardness are alleviated using a compositional approach.Section 2 presents the newly proposed approach,including the stochastic integer programming formulation in detail.In Section 3,the proposed approach is applied at the Tropicana Gold Mining Complex in Western Australia.Conclusions follow.

Fig.1.An example of a mining complex and steps for the proposed approach for integrating a geometallurgical throughput prediction model into short-term stochastic mine production scheduling.

The proposed approach for integrating a geometallurgical throughput prediction model into short-term stochastic mine production scheduling consists of 4 consecutive steps,which are discussed in detail in this section.Fig.1 shows an overview of the required steps for the proposed approach.In the first step,the comminution behavior of the mineral reserve is geostatistically simulated by building additive hardness proportions using penetration rates from blasthole drilling,which is complemented by additional orebody information such as geological domains,material types,and density.In the second step,a material tracking approach considers all material movements from the various pits to the processing plant,including short-term,run-of-mine stockpiles,to inform the throughput prediction model about properties of blended materials that are sent to the processor.In the third step,measured throughput rates of the operating ball mill are utilized and linked with the previously obtained blended rock characteristics of the processed ore.The link is created through a multiple regression model,which predicts ball mill throughput as a function of rock attributes of blended material.Finally,the proposed throughput prediction model is integrated into short-term production scheduling.Note that the comminution circuit further includes a crusher and a high pressure grinding roll (HPGR).However,this manuscript only considers the prediction of ball mill throughput.

2.1.Spatial simulation of material hardness using penetration rates

Measurement while drilling(MWD)data is a centralized collection of mechanical performance indicators obtained from drilling machinery.The monitored performance indicators of the drilling process include,among other entries,rate of penetration,downhole pressure,rotational pressure,rotational speed,air pressure,and vibration.These measurements are routinely collected in operating mines during drilling activities such as exploration,grade control and blast hole drilling.The latter typically comprises dense drilling patterns (up to 6 m×6 m).A compositional approach of geostatistically simulating the hardness of the mineral reserve is presented next,consisting of creating point simulations of penetration rates first,and creating hardness proportions thereafter.

2.1.1.Spatial simulation of drilling rate of penetration

The rate of penetration (ROP) is used in this article for the spatial simulation of material hardness.Note that details regarding data cleaning and pre-processing of ROP entries are discussed in the case study presented in Section 3.The ROP measurements from drilling are converted to the unit seconds per meter,s/m,for downthe-hole compositing.This unit is preferred since compositing(averaging)is performed by down-hole length.The unit s/m is kept for the remainder of the article.The ROP composites form a set of n conditioning data points distributed in space{zROP(uα),α=1,···,n}.Conventional geostatistical simulation techniques[46]can be used to create a set of equiprobable spatial simulations of ROP,which is performed on a regularly spaced grid(point support).

A challenge arises for the change of support of hardness-related variables from simulated grid nodes (point support) to mineable volumes (SMU).Yan and Eaton [16] indicated that comminution behavior could be disproportionally affected by harder fractions of material.Thus,a simple averaging of ROP will create biases[11,19,30,47].An alternative approach is proposed herein,which generates additive geometallurgical,hardness-related variables.Additivity of variables becomes especially important for material tracking and assessing the properties of blended materials,which is discussed in a subsequent section.

2.1.2.Construction of additive hardness proportions

Instead of changing the support scale from grid nodes to mineable volumes,a compositional approach was proposed,which involves the creation of K proportions of softer (easier-topenetrate)and harder(harder-to-penetrate) material within a larger volume (SMU).Considering all simulated values at the point support scale within a mining block,each proportion represents a fraction of point penetration rates that fall within a specific interval of the simulated ROP.All K proportions naturally sum up to 1,or 100%,within each mining block.An illustration of hardness proportions is shown in the table on the right in Fig.2.The required intervals are delimited by a set of global ROP thresholds derived from a set of percentiles of the global cumulative histogram of ROP,as seen on the left in Fig.2.The number of required intervals(percentiles) is discussed in the case study.

2.2.Tracking blended materials through the mining complex

The main challenge for linking the simulated rock characteristics of the mineral reserve to the observed production data at a processing plant is to keep track of the materials extracted from the pits and sent to the processing facilities,as indicated in the second step of the proposed framework(Fig.1).The aim of this step is to characterize the pertinent rock attributes of blended materials entering the comminution circuit so that they can be matched later with observed responses(i.e.,throughput)of the ball mill.In addition to the created hardness proportions,other orebody attributes,such as density,geological domains,and material types of the extracted blocks,are tracked as well.The material tracking approach is schematically shown in Fig.3.The design accounts for multiple mines and includes all short-term stockpiling that occurs before materials enter the processing facilities.

Fig.2.Construction of penetration rate thresholds used to create hardness proportions per mining block.

Fig.3.Schematic representation of material tracking of hardness proportions through the mining complex,including short-term stockpiles.

The material movements in the various pits and stockpiles are replicated in detail using truck cycle data.The necessary information includes the start and end times of truck cycles,as well as spatial coordinates of truck loading and dumping locations.Three types of truck cycles are relevant for material tracking from the mines to the processing facilities: (1) ore is hauled directly from one of the pits to the crusher;(2)ore is hauled from one of the pits to one of the short-term stockpiles;(3) ore is reclaimed from a short-term stockpile and transported to the crusher.In the presented material tracking scheme (Fig.3),all truck cycles starting in the mine (Types 1 and 2) are first linked to their nearest Block ID of the simulated resource model (hardness proportions,rock type,density) using a nearest-neighbor search.Cycles that end at the crusher location (Types 1 and 3) are utilized for calculating the properties of incoming,blended materials.

While Type 1 cycles are easy to trace back to the resource model,a more elaborate tracing scheme is constructed for stockpiled material (Fig.3).Other than blending beds or homogenization facilities [48],short-term stockpiles (run-of-mine pads or fingers)may not have a fixed build/remove scheme,and their position and dimensions change over time.Thus,an initial grid of cells(10 m×10 m)is constructed,one for every stockpile,covering the area of potential material placement.When a truck dumps material on one stockpile (Type 2),the closest cell is determined by a nearest neighbor search using its dump coordinates.If this cell has been active before (previously dumped material),the associated Block ID of the truck cycle is added to a list.Otherwise,the cell is activated,and the associated Block ID of the truck cycle becomes its first entry.When material is hauled from a stockpile to the crusher (Type 3),the closest active Cell ID is found using truck loading coordinates.The respective list of Block IDs informs about all material that has been deposited there for calculation of blends(one-to-many relationship,detailed in Wambeke et al.[49]).An action of deleting all active cells in an individual stockpile is triggered after a stockpile is depleted.

Weighted proportions of all tracked properties are constructed by truck payload (weighting factor) from all Types 1 and 3 cycles that are recorded in a specific interval arriving at the crusher.The time interval for material tracking can be as small as a few minutes or can comprise several days,constrained by the minimum truck dumping frequency recorded in the fleet database.Please note that connecting truck load coordinates to the nearest block and cell IDs are a simplification for material tracking,which can be enhanced by modern blast movement and grade control techniques.

2.3.Building a throughput prediction model

A multiple linear regression model is constructed to predict ball mill throughput(response variable)as a function of rock attributes of blended materials (predictor variables or features),which have been obtained previously in the material tracking step.In the regression model,the ithobserved throughput,yiin a sample of N observations,is expressed as a linear combination of M predictor variables of the tracked rock attributes,fi,j,as shown in Eq.(1).

Note that wj,j=0,1,···,M,are the regression coefficients;and eiis a random variable representing the prediction error.The observed throughput rates,yi,are average observed ball mill throughput rates (milled ore tonnage per operating hour),collected for the same time intervals as the tracked rock attributes,fi,j.The goal of regression is to find the vector of weights,w=[w0,w1,···,wM],that minimize the sum of squared errors(SSE),

where fi=[1,fi1,···,fiM]is a vector of predictor variables for the ithobservation.The closed-form solution of the minimization of SSE with respect to the vector of weights,written in matrix form[50,51],is.

where X ∈RN×RM+1is the feature matrix;and y ∈RNthe vector of all observations.The resulting weights,,are used as input for production scheduling.A discussion of feature selection and regression results,including an analysis of correlation to ball mill throughput,are presented in Section 3.

2.4.Integration of throughput prediction model into stochastic shortterm mine production scheduling

The mathematical model which integrates the throughput prediction model into the simultaneous stochastic short-term optimization of mining complexes is presented next.The stochastic optimization model is formulated as a non-linear two-stage stochastic mixed integer program with recourse [52].Stochastic inputs are provided in the form of geostatistically simulated orebody scenarios to account for geological,or supply uncertainty in production scheduling.A general model for the simultaneous stochastic optimization of mining complexes has recently been introduced by Goodfellow and Dimitrakopoulos [43,53],which is used for the optimization of life-of-mine production schedules.However,the mathematical model needs to be extended to consider critical operational aspects of short-term planning and production scheduling.The extension shown in this article addresses the integration of the proposed throughput prediction model.A simulated annealing metaheuristic is used to solve the non-linear optimization model.Simulated annealing is chosen due to its successful utilization for open pit mine production scheduling for single open pit mines [54-56] and mining complexes [41,53].The typical runtimes documented for the simulated annealing algorithm for long-term production scheduling in mining complexes with more than 100000 blocks to be scheduled[43,44]are reduced significantly for short-term production scheduling(typically thousands of blocks scheduled),staying below 2-h run time to reach a solution.Details of the adapted simulated annealing algorithm and parameter selection for short-term planning can be seen in [4].

All new parts related to the integration of the developed ball mill throughput prediction model for short-term planning are discussed in this section,whereas detailed explanations of the general model can be retrieved from [43,53].The utilized indices and sets of the stochastic optimization model are listed in Table 1,whereas the used parameters and variables are presented in Tables 2 and 3,respectively.

Table 1 Indices and sets used in the stochastic optimization model.

Table 2 Parameters used in the stochastic optimization model.

Table 3 Variables used in the stochastic optimization model.

Objective function.The objective function shown in Eq.(4)aims to maximize the gross profit generated in the mining complex within the short-term planning horizon while simultaneously minimizing the risk of not meeting short-term production targets.

Part 1 of the objective function sums all revenues and costs generated within the mining complex.The formulation is presented in generalized form,to account for the multitude of possible revenues and costs which occur during the extraction and transformation of the rock,including but not limited to revenues from recovered metal,mining costs,processing costs and rehandling costs from stockpiled material.Note that the related quantities can be derived from linear or non-linear functions applied on other attributes in the mining complex,.For example,the quantity of recovered metal,vj,l,s,t,at a processing location,l ∈P,can be calculated as a function of total metal received,vk,l,s,t,at the processor l in a scenario s and period t,accounting for a nonlinear metal recovery function.Part 2 is a stochastic component that penalizes deviations from production targets given orebody uncertainty,such as mine production capacities,grade blending constraints,and more.The first two parts are well-established building blocks of simultaneous stochastic optimization of mining complexes,both in long-term[43]and short-term planning[4,57].Part 3 is a new component introduced in this article and penalizesdeviations from the mill tonnage that can be processed based on the throughput that is calculated using the proposed prediction model.Part 4 ensures that the resulting production schedule is mined in a connected (smooth) pattern.

Constraints.Details of the new constraints related to ball mill throughput are given below.Other included constraints are mine capacity constraints,blending constraints,reserve constraints,slope constraints,destination constraints,and material flow constraints,which are detailed in Goodfellow and Dimitrakopoulos[43,53].Note that the reserve constraints in this formulation enforce all blocks to be extracted by the end of the short-term planning horizon.

Complementing Part 3 of the objective function,new components are added to integrate the developed throughput prediction model into the simultaneous optimization of mining complexes.Consider the periods t for a short-term horizon T.The resulting short-term production schedule is informed by a set of orebody simulations,S,containing attributes that are used for throughput evaluation,j ∈HTPH,including hardness proportions and other attributes important for scheduling,such as metal grades.All throughput-related attributes of blended materials sent to the processing location,l ∈P,in period t,and orebody scenario s,are represented by vj,l,s,t,j ∈HTPH.This variable is calculated as a weighted average of block attributes,βb,j,s,extracted in t by their respective tonnage,Tonb,and sent to the processor l,shown in Eq.(5).The formulation also considers material attributes sent from locations other than mines,l'∈S,i.e.,stockpiles and external sources,to the processor l.

Note that the calculation of vj,l,s,tin Eq.(5)is non-linear,since it represents the attributes j of the blended material stream sent to the processor l from mines and stockpiles in period t under orebody scenario s,such as the grade or density of the blend,the proportions of harder and softer material sent to the processing plant,the fraction of weathered materials etc.The throughput calculation for a processing location for any given period t and orebody scenario s at location l is shown in Eq.(6).

Throughput rates are calculated utilizing the optimized weights,w*,obtained in Eq.(3).Specifically,the prediction model is queried at every iteration of the metaheuristic solution method,using the current production schedule.The throughput prediction obtained in Eq.(6) is used to assess the tonnage that can be processed by the ball mill,as shown in Eq.(7).

This constraint penalizes a deviation between the scheduled tonnage sent to the processing unit and the calculated tonnage in Eq.(7).The latter tonnage represents the ore tonnage that can actually be processed by processing location l,which is now based on pertinent rock attributes related to hardness.

In this section,the proposed approach of integrating a ball mill throughput prediction model into simultaneous stochastic shortterm production scheduling is applied at the open pit Tropicana Gold Mining Complex,which is located 330 km East-North-East of Kalgoorlie in Western Australia.The gold mining complex comprises 4 contiguous pits extending 6 km in strike length from north to south.Ore is sent to a single processing facility,consisting of a comminution circuit and a carbon in leach (CIL) facility.The comminution circuit comprises a primary crusher,a high-pressure grinding roll (HPGR),and a ball mill.

The gold mining complex is operated as a typical truck and shovel open pit operation.After fragmentation through conventional drilling and blasting,ore is hauled by truck directly to the crusher or to one of eight short-term,run-of-mine (ROM) stockpiles located in the vicinity of the primary crusher.Temporarily stockpiled material amounts to 80%-90%of processed ore;thus,it is crucial to include material of ROM stockpiles in material tracking,as shown in the case study(Section 3.2).In this case study,the datasets used for the throughput prediction model include spatially recorded penetration rates from blast hole drilling and measured throughput rates of the operating ball mill over the horizon of one year.This information is supplemented by orebody models providing geological domains,material types,and rock density.The material tracking relies on truck cycle data collected from a fleet management system.

3.1.Hardness proportions

The raw,recorded penetration rates used to create hardness proportions stem from MWD data collected by all drilling activities(blast hole,grade control,pre-split drilling)routinely performed in the mining complex.Fig.4 presents cross-sections of the recorded rate of penetration (ROP),shown in seconds per drilled meter (inverse penetration rate,Fig.4a),and drilled meters per hour(Fig.4b),which can be visually compared to the geological domains of the orebody (Fig.4c) and the weathering profile/regolith/material types of the deposit(Fig.4d).Large-scale geological structures,such as dipping geological units displayed in Fig.4c,can be clearly distinguished by consistent differences in ROP.Furthermore,the five stratigraphic regolith units or material types of the orebody (weathering profile,Fig.4d) are strongly reflected in ROP,indicating weathered,easier-to-penetrate rock towards the surface.Also,the harder-to-penetrate,fresh rock,situated in greater depth,is easily detectable by penetration rates.

It is important to note that biases in this database may exist due to several factors,including the utilization of different drill rig types,multiple rig operators,the degree of bit-wear,and varying drilling tasks.Thus,the data can contain entries that are less correlated to properties of the penetrated rock but rather relate to operational procedures.Many normalization procedures have been proposed in the technical literature [35,36,60].However,these normalization procedures need to be applied carefully.Internally created variables of the MWD system in the present dataset,such as the blastability index (BI) [31],were found to be highly biased towards individual drilling machines.Furthermore,the application of the adjusted penetration rate(APR)(Zhou et al.[35])resulted in new,machine-related biases and reduced geological correlation compared to ROP.This shows that site-specific data handling is required.The following targeted outlier removals of ROP data entries were applied on the present dataset,proving to be effective in this case study:

(1) Production drilling was utilized exclusively (blast holes),removing pre-split and grade control drilling from the dataset.The reason for this is that blast hole drilling is conducted solely with drill rigs of the same type in this mining complex.

(2) Due to pre-fragmented rock from previous blasting activities(easier to penetrate,resulting in higher ROP),the first 2 m of every hole were discarded.

(3) Data outliers due to some irregular events (enlarging the drill rod,flushing the hole,etc.) were identified and removed.

3.2.Material tracking in the mining complex

A fleet management system installed in the mining complex records each individual truck cycle in a central database and can therefore be used to replicate material movements in the whole mining complex,as described in Section 2.2.Fig.5 visualizes how material tracking in short-term stockpiles (run-of-mine fingers)is performed at the mining complex.

Fig.4.Cross sections (E-W) of penetration rates retrieved from measurement while drilling data,compared to cross sections (E-W) of geological domains and weathering profile of the orebody.

Fig.5.Top view of material tracking in short-term stockpiles(run-of-mine fingers)and a close-up of a single stockpile (RF06).

The left-hand side in Fig.5 shows a top view of all ROM Fingers(stockpiles) located in the vicinity of the primary crusher.A closeup of the sixth stockpile(RF06)shows all dumping events of stockpiled material coming from pits(red markers)and reclaimed material sent from the stockpile to the crusher (black markers).Furthermore,the created grid cells are indicated,which are used to spatially inform and distinguish material within the stockpiles depending on where it has been deposited.

3.3.Multiple regression implementation and parameter testing

The material tracking is conducted in daily time intervals for a period of one year.As part of data preparation for multiple regression,a moving average is calculated for both blended rock properties (predictor variables) and observed ball mill throughput(response variable).In this way,similar time intervals are preserved for inputs and outputs of the regression.Fig.6 shows the application of a 7-d moving average on daily average ball mill throughput observed at the mining complex.

The moving average is applied because the daily throughput rates observed at the mine show very noisy behavior.It can be seen in Fig.6 that a 7-d moving average helps emphasize trends of higher and lower throughput rates over longer periods,which are more likely connected to rock properties of the processed material.At the same time,the influence of technical disruptions and shortterm delays of ore feed is reduced.Additionally,an important reason for applying a 7-d moving average is that the prediction model is aimed to predict a ball mill throughput for production scheduling comprising weekly periods.By applying the 7-d moving average,data fits the time scale on which ball mill throughput predictions will be used for in production scheduling.

Fig.6.Observed ball mill throughput showing daily average and 7-day moving average for material tracking.

From the available 365 d,a 180-d window is retained for regression analysis,discarding about three months of production at the beginning and end of the observed interval.This ensures that (1)material deposited in the short-term stockpiles is known,and (2)material in later time intervals is retained so that it can be used for production scheduling.Table 1 shows statistical information of the dataset used for the regression analysis,which includes the tracked rock attributes of blended materials (predictor variables) and the observed ball mill throughput (response variables).Furthermore,correlations between the predictor variables and the response variable are shown in the last column of Table 4 using Pearson’s correlation coefficient shown in Eq.(9),

An important observation of Pearson’s correlation coefficients shown in Table 4 is that the strongest linear relationships exist between the constructed hardness proportions (A1-A10) and ball mill throughput.The softer hardness proportions (easier-topenetrate,A8-A10) correlate positively with ball mill throughput,indicating that the throughput increases when more of this softer material is processed in the ball mill.On the other hand,increased proportions of harder-to-penetrate materials lead to a decreased ball mill throughput,indicated by a negative correlation coefficient(A1-A6).The consistency of positive and negative correlation among hardness proportions supports the hypothesis that penetration rates recorded by drilling machines have the capacity to predict comminution performance/grindability of the rock.Other correlation coefficients in Table 4 indicate that the originating geological orebody domains (B1-B4) and average material density(D1)have less influence on ball mill throughput in this case study.However,the different degrees of weathering/regolith/material types(C1-C2)show a stronger correlation to throughput,whereas higher proportions of weathered material are positively correlated with ball mill throughput.This result is not surprising since higher degrees of weathering can cause degradation in rock competency,leading to better grindability in the comminution circuit [61].

Table 4 Pearson’s correlation coefficient between rock attributes of blended materials in the comminution circuit and observed ball mill throughput.

3.3.1.Number of hardness proportions

The construction of hardness proportions requires a choice of the number of thresholds and their distribution of percentiles,which can be chosen to be evenly or unevenly distributed.Fig.7 shows variations of splitting penetration rates into different hardness proportions and their effect on the multiple regression using only the constructed proportions as predictor variables (features).Leave-one-out cross-validation is applied [62],which consecutively builds a regression model using all data points but one,then predicts the value of the left-out data point.This process is repeated for every data point.

A first result of regression analysis in Fig.7 is that hardness proportions enable a prediction of ball mill throughput with a root mean squared error (RMSE) of less than 30 t/h.Additionally,it can be observed that splitting the cumulative histogram of penetration rates differently(even split,soft tail,hard tail)leads to different prediction behavior.The split of softer material appears to predict medium/high observed throughput rates better,whereas the split of harder material clearly predicts lower throughput rates better.Note that average-type approaches of hardness relying on single-valued hardness indicators such as Wi and SPI cannot capture these observed differences,which are driven by varying proportions of harder and softer rock.When comparing the 5hardness and 10 hardness proportions,splitting into more categories appears to have lower RMSE.

3.3.2.Generalization potential of various features

The material density,geological orebody domains,and degree of weathering are also tracked next to hardness proportions,as listed in Table 4.A comparison of their potential to predict ball mill throughput in the present case study is shown in Fig.8.Here,the various regression models are built on a data set comprising 80%of data points.The models are then tested on the remaining 20%of data points,revealing the generalization potential to unseen data of the tested feature sets.Note that testing on a consecutive time series of throughput data is more informative than a randomly picked test set since random test samples can be very similar to neighboring training sample points.Also,note that the yaxis changes the scale for some graphs in Fig.8 depending on the spread of throughput predictions.

It can be observed that the use of 5 hardness proportions shows the best results among all tested cases,shown in Fig.8a (even split),Fig.8b(soft tail),and Fig.8c(hard tail).Splitting the soft tail of the penetration rate distribution shows the highest correlation coefficient of 0.80(Fig.8b),whereas splitting the harder tail results in the lowest RMSE of 21.88 t/h (Fig.8c).

As expected from correlation analysis,adding the average rock density of blended materials (Fig.8e) as well as adding geological domains as features(Fig.8g)do not influence prediction capability.Interestingly,the addition of the weathering degrees(regolith/material types) as features leads to a higher correlation between test data and predictions(Fig.8f).This can be explained by the notable correlation between the degree of weathering and ball mill throughput (Table 4).However,the generalization capability of weathering proportions is ambiguous,leading to a large overestimation of throughput in the first 5 d of the test set and an increased RMSE of 50.38 t/h (Fig.8f).The generalization capabilities of higher splits into 10 hardness proportions(Fig.8d)and polynomial terms(Fig.8h)are poor.Overfitting is clearly visible for 10 hardness proportions,which did not yet appear for previous leaveone-out cross-validation (Fig.7) because of the small ratio of testing vs training data.The addition of second-order polynomials(Fig.8h) worsens the generalization potential further,resulting in even larger overfitting.

The approach presented in this article is general and allows the utilization of any rock attributes contributing to more accurate throughput prediction.This can also potentially include head grades of metals and other (deleterious) elements fed to the ball mill next to the geological attributes used herein.In the present case study,the most robust predictions are obtained from the utilization of 5 hardness proportions,as seen in Fig.8.This result forms the basis for the attributes used for short-term production scheduling,which is presented next.

3.4.Short-term production scheduling

For weekly short-term production scheduling of the described mining complex,three months of planned production are used as input to optimization (3874 mining blocks,12 m×12 m×12 m).Various sets of geostatistically simulated orebody scenarios serve as stochastic inputs for production scheduling to quantify geological uncertainty of metal grades and geometallurgical attributes,which is displayed for mining area A in Fig.9.Uncertainty of metal grade (Au) is accounted for by 20 equally probable orebody scenarios.The required number of orebody scenarios for mine production scheduling has been examined by previous authors,showing that 10-15 simulations provide proper stable schedules and forecasts [63,64].Note that each of the 5 hardness proportions forms an individual geometallurgical attribute for production scheduling and is represented by twenty orebody simulations each,as indicated in Fig.9.

A 12-week short-term production schedule is generated for the mining complex with the stochastic optimization model described in Section 2.4.The stochastic mixed integer programming model and the simulated annealing algorithm to solve the model are written in C++programming language and is run on a computer with an Intel Core i5-7200U CPU @ 2.50 GHz processor and 8 GB RAM,using Windows 10 Home as operating system.The algorithm is run for 200-000 iterations with a running time of 1.7 h.A total of 2 mining areas are simultaneously scheduled for production.Fig.10a shows the resulting 12-week production schedule for both mining areas and a risk analysis of the tonnage that is expected to be processed in the ball mill.More specifically,Fig.10a shows the difference between the scheduled ore tonnage (blocks scheduled for production) and the tonnage that can be processed by the ball mill using the relevant geometallurgical attributes and the integrated throughput prediction model.The new stochastic constraints used in this optimization model (Eq.(8)) are directly linked to the objective function,as explained earlier.

Fig.7.Regression analysis testing multiple splits of hardness proportions.

Fig.8.Generalization potential of various rock attributes testing on 20% unseen test data.

It can be seen in Fig.10a that the expected difference between scheduled and processed ore for the created 12-week production schedule(Fig.10a,solid line)is small in all scheduled periods,with a maximum expected deviation of 1 kt(0.6%)per period and a total deviation from production targets of less than 5 kt (0.3%) for the complete scheduling horizon (3 months).Note that the use of a set of geological orebody scenarios as input to optimization enables the evaluation of the risk of falling short of or surpassing the required ore tonnage per period.There is a remaining risk involved of over and underfeeding the ball mill indicated by P10 and P90 risk profiles (Fig.10a dashed lines).However,the spread of the profile is kept well below 5 kt per period.The result indicates that the risk of not meeting scheduled ball mill production using the proposed method is low and can be controlled.

Fig.9.Orebody simulations of mining area A serving as input for short-term production scheduling using an integrated throughput prediction model.

Fig.10.Risk analysis showing the difference between scheduled mill tonnage and realized tonnage of a production schedule using all newly developed components and using conventional mill tonnage constraints.

To evaluate the benefits of the new components for short-term production scheduling,the case study is repeated without the use of the throughput prediction model.A conventional mill capacity constraint is used instead,which limits the processed material per period based on a constant tonnage capacity.Fig.10b shows the resulting production schedule and a risk analysis of the mill tonnage that can be processed using the tonnage capacity constraint rather the throughput prediction model.Fig.10b shows that the differences between scheduled tonnage and processed tonnage can be significant by using conventional capacity constraints.The largest deviations can be recognized in periods 3 and 12.In both periods,an expected value of over 10 kt of ore cannot be processed by the ball mill,which amounts to a deficit of 7% of the planned tonnage in each of these periods.This is indicated by the P50 value of the risk profile in Fig.10b (solid line).The risk evaluation also shows that there is a small probability of underfeeding the ball mill in some periods,indicated by the P90 risk profile exceeding the balance line.However,the expected cumulative shortfall of ore over the complete scheduling horizon accumulates to an expected value of 61 kt,which is equivalent to 3% of total production that cannot be achieved.

Mill capacity constraints that limit the amount of ore tonnage sent to the processing facility per period are standard practice,both in stochastic and deterministic production scheduling.However,these tonnage constraints assume a constant,average throughput of the milling and grinding circuit for every period and thus ignore the natural heterogeneity of rock attributes of processed ore,leading to highly variable throughput in the short-term,as it can be observed from daily and weekly production data displayed in Fig.6.The newly presented and integrated prediction model,however,is capable of assessing the throughput and the associated risk profile to be expected from scheduled material sent to the processing plant.Fig.11 shows the ball mill throughput prediction for both schedules generated.

Fig.11.Risk analysis of expected ball mill throughput of a production schedule using all newly developed components and using conventional mill capacity tonnage constraints.

While variability in throughput can be observed for both schedules,the variations in expected throughput from period to period are particularly large when conventional tonnage constraints are used (Fig.11b).There is a decrease of expected throughput from 920 t/h in period 2 to 884 t/h in period 3(-4%).Considerable risks in production shortfall are the consequence,as seen in Fig.10b.These results emphasize the need for an integrated geometallurgical approach in mine production planning that creates short-term production schedules while integrating key metallurgical processes further downstream in mining complexes.

This article presents a novel 4-step approach to integrate a geometallurgical throughput prediction model of the ball mill into short-term stochastic production scheduling.The usefulness of the approach is shown in a large-scale open pit mining complex in Western Australia using real-world production data,whereas the key achievements and conclusions are as follows:

(1) Penetration rates stemming from measurement while drilling data are used for the first time for the prediction of ball mill throughput,informing the hardness and comminution behavior of the mineral reserve.

(2) The creation of hardness proportions with the use of penetration rates avoids biases typically introduced by the change of support and blending of non-additive geometallurgical properties related to hardness,whereas the hardness proportions created show the highest correlation to ball mill throughput among all tracked rock attributes.

(3) The use of 5 hardness proportions enables throughput predictions with an RMSE of less than 30 t/h and a correlation coefficient of up to 0.8.

(4) By integrating the throughput prediction model and the newly developed stochastic constraints into the stochastic optimization of mining complexes,the weekly production schedule can achieve planned production with high certainty(expected shortfall from total planned production less than 0.3%)by enabling the matching of the scheduled materials with the predicted performance of the ball mill.

(5) Future ball mill throughput rates are predicted within the optimization model as a function of blended rock properties,overcoming the shortcomings of previous geometallurgical models integrated into production scheduling that assess geometallurgical properties block by block.

(6) Risk analysis of a weekly stochastic production schedule using conventional tonnage capacity constraints reveals that an expected value of more than 7% of planned tonnage cannot be processed by the ball mill in certain periods,which amounts to an expected deficit of 3% of planned tonnage over the 3-month planning horizon.

(7) Given that the penetration rates informing rock hardness and the measured throughput rates of the ball mill stem from production data,the novel approach can give a financial advantage compared to laboratory grinding tests for throughput prediction,which are typically costly and timeconsuming to obtain.

The robustness of the proposed approach relies on several factors,which deserve further investigation in future studies.First,the developed approach relies on the ability to handle noise and biases contained within the collected measurement while drilling data.More studies linking penetration rates to comminution behaviour are welcome to test the reliability of the proposed approach,and additional machinery data such as loader cycle time and loader availability rates could be included to enhance prediction.The throughput prediction capability also depends on the accuracy of the material tracking from pits to processors.Next to the utilized truck cycle data herein,the monitoring of blast movements could be included in the future to improve tracking accuracy.Future work also aims to further improve the prediction of non-additive metallurgical responses of processing plants by utilizing supervised learning models beyond multiple linear regression,which can account for non-linear relationships between predictor and response variables.For future prediction models,measured processing parameters,such as ball mill power and particle size distributions,should be considered given their direct influence on throughput rates.In addition,production data generated from other parts of the processing plant,such as mineral separation processes,should be considered to better integrate value-driving geometallurgical properties into short-term mine production scheduling.

Acknowledgements

This work was supported by the National Sciences and Engineering Research Council of Canada (NSERC) under CDR Grant CRDPJ 500414-16,NSERC Discovery Grant 239019,and the COSMO mining industry consortium (AngloGold Ashanti,BHP,De Beers,AngloAmerican,IAMGOLD,Kinross Gold,Newmont Mining,and Vale).Special thanks are in order to AngloGold Ashanti Limited and the Tropicana Gold Mine (AngloGold Ashanti Australia Ltd 70% and manager,IGO Ltd 30%),in particular Mark Kent,Tom Wambeke,Aaron Caswell,Johan Viljoen and Louis Cloete for providing the data used in this study and long-standing collaboration.The authors would also like to thank K.G.van den Boogaart and the anonymous reviewers for their valuable comments that helped improve the manuscript.