Comparison of Alternative Strategies Estimating the Kinetic Reaction Rate of the Gold Cyanidation Leaching Process.

It is very important to establish an accurate process model for implementing further control and optimization of the gold cyanidation leaching process. Unfortunately, the important kinetic reaction rates affecting process operation are unmeasurable, and moreover, their estimation is an ill-posed inverse problem due to the fact that the noise in concentration measurements is easy to be amplified and propagated into the rate estimates by derivative operation. In this paper, the alternative strategies (finite difference, polynomial fitting, Savitzky-Golay filter, wavelet decomposition, and Tikhonov regularization) estimating the kinetic reaction rate of the gold cyanidation leaching process are investigated in detail. The simulation results show that the direct finite difference leads to poor estimating results for the noisy case and the other strategies are capable of avoiding the noise amplification and improving the estimating results to some extent. In all of the investigated strategies, the Tikhonov regularization leads to the satisfactory and acceptable estimating results in both the noiseless and noisy cases, which will lay an important foundation for the subsequent model identification, production index prediction, and operation optimization.

Introduction

As one of the two extractive metallurgy technologies, hydrometallurgy is a process of using a solvent to separate the metal from raw mineral materials by chemical reactions, which has the advantages of high metal recovery, good flexibility, simple equipment, good environmental conditions, low investment, quick effect, and good recovery of associated components. Moreover, hydrometallurgy processes are easier to fulfill continuous and automated production and meet the requirements of the sustainable development of mineral resources.[1−3] In a typical hydrometallurgical plant, these often consist of the following fundamental operation units: ore comminution, size classification, gravity concentration, pulp dewatering, leaching, recovery, elution, electrolytic extraction, melting, and casting.[4,5] As the most important procedure of hydrometallurgy processes, the quality of leaching solution directly affects the final economic benefits strongly. Although leaching technology in hydrometallurgy has reached a relatively advanced level, the level of automation in the real leaching process is still not high and still remains in the way of offline test and analysis in laboratory, manual adjustment by experience and manual control, which leads to low production efficiency, high consumption of materials, high production cost, and low profit.[6−8] Therefore, under the condition that the quality of the leaching products is ensured, reasonable distribution of raw materials and then reduction of the total consumption, production cost, and improvement of economic efficiency have become serious problems that need to be solved urgently in hydrometallurgy industrial processes, which are solved based on the online prediction of key variables and process optimization.[9−13]

To implement the advanced optimization and control on the leaching process successfully, the key production index (gold recovery or gold grade in cyanide residue) is necessarily measured and used to guide the further production operation. However, the traditional concentration sensors are not reliable when running in the complex industrial conditions for a long time. Moreover, the subsequent maintenance cost is too high for the production enterprise to afford. The soft measurement technique is an appropriate alternative to solve the above prediction problem of key production index instead of the hard sensors.[14−16] The mechanism model of the leaching process is composed of the mass conservation equations of key components, namely, gold in the solid, gold in the liquid, and cyanide ion in the liquid, as well as the corresponding kinetic reaction rate models.[17−21] Where the mass conservation equations of gold and cyanide ion are deterministic, both of the kinetic reaction rate models are unknown and only an approximate empirical model can be obtained. To identify the empirical model parameters accurately, the model outputs (kinetic reaction rates) should be obtained in advance. Unfortunately, in real industrial process, the gold dissolution rate and the cyanide-ion consumption rate are unmeasurable or measuring them is costly in the real plant, which brings inevitable troubles in modeling the mechanism model. Therefore, both of the above rates need to be estimated in advance accurately according to the measurable process data. Whereas, the rate estimation is an ill-posed inverse problem and the noise in concentration measurements is easy to be propagated into the rate estimates by derivative operation.[22,23] The most commonly used method to obtain the discrete data derivative in process engineering is the known finite difference method,[24−26] in which the current time data as well as one or several nearby data are used together to approximate the current time derivative. When the time interval is short enough and the data noise is absent, the direct finite difference method is a simple and effective tool to estimate the rate by the concentration measurements. To avoid the noise amplification by direct numerical differentiation, several effective smoothing techniques have been developed to obtain accurate derivative estimates of noisy data.[23,27−32] One of the most widely used procedures for differentiating noisy experimental data is the polynomial fitting (PF) method, in which a smooth polynomial function is fitted to the experimental data first, and then the fitted polynomial function is differentiated analytically to obtain the derivatives directly.[33,34] As a particular case of polynomial fitting methods, the popular Savitzky–Golay (SG) method is generally used to smooth and preprocess the experimental measurement data in industrial practice.[27] By the local least-squares regression method, an nth-order polynomial is fitted to each interior data point and m of its uniformly spaced neighboring points on either side. Recently, a new method of processing the time–concentration data of reaction kinetics based on the Tikhonov (TIK) regularization method has been proposed by Yeow et al.,[23,32,35] in which the kinetic reaction rate models are not necessary and the noise amplification and propagation can be kept under control effectively as much as possible. In this paper, the commonly used methods for estimating the derivatives in industrial practice are investigated detailedly and a systematic performance comparison of various alternative strategies estimating the kinetic reaction rates of gold cyanidation leaching process is provided.

The rest of this paper is organized as follows. Section introduces the gold cyanidation leaching process under investigation and gives the mechanism model in detail, which is used to simulate the real gold cyanidation leaching process and prepare the simulated data used for evaluating the performance of various alternative estimating strategies. Then, a detailed description of various alternative approaches estimating derivatives is given in Section , which includes the traditional finite difference method, the polynomial fitting method, the Savitzky–Golay filter method, the wavelet decomposition (WD) method, and the Tikhonov regularization method. In Section , first, the mechanism model of Section is used to simulate the real leaching plant and generate some necessary data sets used by the subsequent simulation tests. Second, the above various alternative methods are applied for estimating the kinetic reaction rates of gold cyanidation leaching plant. Finally, a detailed discussion of the above simulation test results is presented. Conclusions and some potential future research directions are presented in Section .

Process Description

Leaching Process

Due to the facts that the world’s ore grade of available gold mining is continuing to decline, the comprehensive utilization of mineral resources is becoming increasingly urgent, the requirement for environmental protection is increasingly strict, and the purity requirement of products is continuing to increase, some countries around the world have begun to attach great importance to the hydrometallurgical industry.[1,3,6,8] In 1887, the modern cyanidation gold extraction technology was first put into industrial production, which is one of the most important gold extraction methods at present. The typical cyanidation gold extraction process usually includes: gold cyanidation leaching, washing and filtration of the leached slurry, extraction of gold from the pregnant solution, and smelting of the gold slime.[1,4,17] The simplified plant flowsheet of the cyanidation gold extraction process is given in Figure . In the gold cyanidation leaching process, the insoluble solid gold is extracted from the ore into the liquid by a leaching agent (usually cyanide solution, such as NaCN because of its strong ability of dissolving gold, high stability, and low cost). The gold cyanidation leaching is essentially considered as an electrochemical corrosion process, and the solid gold in the ore is mostly oxidized and synthetized into a stable and soluble complex ion [Au(CN)2]−, the performance of which has a significant impact on the production benefit of the subsequent procedures. Gold cyanidation leaching usually occurs in continuous stirred tank reactors according to the following two chemical reaction eqs and 2.[5,18−21,36] In practice, the serial reactor network architecture is usually adopted for increasing the final gold recovery and avoiding the not-leached gold loss as much as possible. And in each leaching tank, the same chemical reaction occurs as described below

Figure 1

Simplified plant flowsheet of the cyanidation gold extraction process.

The whole reaction equation can be easily obtained by summing the left- and right-hand sides of eqs and 2, respectively

Before the gold cyanidation leaching operation, some necessary preprocessed operations need to be carried out to improve the leaching effect. In the foregoing flotation operation unit, several kinds of reagents are added to promote the flotation; however, it is detrimental to the subsequent leaching operation. Therefore, the flotation concentrate ore is first pumped into a thickener and filter press to remove the remaining harmful reagents mostly. Then, the ore pulp with an appropriate concentration is obtained by mixing the gold bearing filter mass with the recycled wastewater containing CN–, which is virulent and prohibited from discharging into the environment. Finally, the conditioned ore pulp is stored in a buffer tank and prepared for the subsequent leaching operation. The ore pulp in the buffer tank continuously overflows into the serial leaching tanks with a constant flow rate, which is controlled by the flow rate of the conditioned ore pulp pumped into the buffer tank. In each tank, the gold particle in solid reacts with the cyanide ion and dissolved oxygen in liquid and synthetizes to a stable and soluble complex ion [Au(CN)2]− to separate from other impurities, as described in eqs and 2 or eq . The reactant CN– primarily derives from the input ore pulp conditioned with the recycled wastewater including a small amount of CN– and the newly added NaCN solution with 30% concentration in each leaching tank. Another important reactant, the dissolved oxygen, derives from the pumped compressed air, which generates moderate pneumatic stirring effect to prevent the pulp from heaping up at the bottom of tanks and meanwhile enhance the reaction intensity. Then, the leached ore pulp overflows into the two-stage thickener, and after washing, the pregnant solution is stored in a pregnant solution tank and gradually pumped into the subsequent process of gold recovery by zinc.[9−11]

Mechanism Modeling

To evaluate the performance of various alternative strategies estimating the kinetic reaction rates of gold cyanidation leaching process, the real values of rates should be known. However, unfortunately, the gold dissolution rate and the cyanide-ion consumption rate are unmeasurable in real industrial process, which results in that we have to employ the existed mathematical model to simulate the real gold cyanidation leaching process. In our study, the mechanism model developed by De Andrade Lima and Hodouin is used to generate the necessary simulated data. The corresponding parameters of the gold dissolution rate and cyanide-ion consumption rate models have been identified and calibrated based on a set of real industrial data from an Australian gold processing plant. And the calibrated model has been successfully applied to the industrial gold processing plant and satisfactory results have been obtained.[17−21] Therefore, the above mechanism model is applicable and competent to the task of simulating the reality of gold cyanidation leaching process.

Without loss of generality, only the single-stage gold cyanidation leaching process is investigated in this paper. The ith-stage mechanism model of the gold cyanidation leaching process primarily consists of the following five equationswhere eqs –6 are the mass conservation equations of gold in the solid, gold in the liquid, and cyanide in the liquid, respectively. Equations and 8 are, respectively, the corresponding kinetic reaction rate models whose parameters have been identified based on the real industrial data from an Australian gold processing plant.[18,19,21] It should be noted that all of the independent variables t in the involved variables are removed in the above mechanism model for simplification. The schematic diagram of the mechanism model for gold cyanidation leaching process is illustrated in Figure . Qs and Ql are the solid and liquid flow rates, respectively, and have the relationship of Ql = Qs(1/Cw – 1) with the parameter Cw as the weight concentration of solid in the pulp. QCN is the flow rate of cyanide added into the tank. Ms, Ml are the solid and liquid holdup in the tank, respectively, which can be determined according to the formulas Ms = Qs × τ and Ml = Ql × τ with τ = V/(Qs/ρs + Ql/ρl) × 1000 as the average residence time in the tank, where V is the net active volume of the leaching tank and ρs and ρl are the corresponding solid and liquid densities. Cs, Cl, CCN, and Co are the gold grade in the ore, gold concentration in the liquid, cyanide-ion concentration in the liquid, and the dissolved oxygen concentration in the liquid, respectively. Cs represents the ideal residual gold grade in the ore after the leaching operation, which is a function of the average particle diameter d̅ of gold orewhere rAu and rCN are the kinetic reaction rates of leaching process, that is, the gold dissolution rate and the cyanide consumption rate, respectively, which are affected by Cs, CCN, Co, and d̅ and are the most key aspect for modeling the gold cyanidation leaching process.[9,10,17,19] Usually, only an approximate empirical kinetic reaction rate model can be obtained and most of model uncertainties result from the kinetic model essentially.[9,10] The above mechanism model can be easily solved using an iterative numerical method based on the Matlab software.[37,38] Finally, the gold recovery can be calculated by the following formula

Figure 2

Schematic diagram of the mechanism model for the gold cyanidation leaching process.

Kinetic Reaction Rate Estimation Based on the Concentration Measurements

As mentioned above, it is the premise of predicting the gold recovery to identify the kinetic reaction rate models and parameters accurately. However, the two kinetic reaction rates are unmeasurable in real leaching process, and therefore, both of the rates need to be estimated based on the measurable variables affecting them principally in advance. As is well known, estimating the rate based on the concentration measurements is an ill-posed inverse problem,[22,39,40] which can be easily observed by the mechanism model eqs –6. To obtain the rate rAu or rCN, we have to take the derivatives of the concentration measurements Cs or CCN with respect to time. If the improper estimating methods are adopted, the measurement noise in Cs or CCN is easily amplified and propagated into the rate estimates and therefore affect the kinetic model identification results directly. In this paper, various common alternative approaches estimating derivatives are investigated in detail, which includes the traditional finite difference method, the polynomial fitting method, the Savitzky–Golay filter method, the wavelet decomposition method, and the Tikhonov regularization method. For simplification and consistency, the problem of estimating the rate by concentration measurements is rewritten in the following general formwhere C(t) is the generalized concentration measurements at time t (obtained by transforming the mechanism model eqs –6 into the form of eq ) and r(t) is its corresponding derivative at time t, namely, the kinetic reaction rate. In the gold cyanidation leaching practice, the concentration measurements are often obtained by offline experimental analysis and hence are not sampled continuously due to the limitation of production technology, measurement sensors, and production cost. The above estimating problem is essentially to conduct a discretely numerical differentiation on the concentration measurements at discrete points in time.

Finite Difference

The simplest way to obtain the derivatives of the given data at discrete points is the finite difference method,[24,41−43] which uses the current and its nearby values to estimate the current derivative. Essentially the real unknown derivative is approximated by the difference quotient and the common derivative formulas using finite difference are shown as followswhere eqs –14 are the distinguished forward finite difference (FFD), backward finite difference (BFD), and central finite difference (CFD) formulas, respectively. t is the ith sample time of the concentration measurements between the initial time t0 and the end time t. Δt = ··· = t – t = t – t = ··· is the sample interval between any two successive sample times, namely, the step size. The above three difference formulas can be easily derived by making some simple operations on the Taylor expansions of C(t) and C(t) at the point t.[24,43] The geometrical meanings of the three approximation methods are very obvious and can be demonstrated clearly in Figure . It can be shown in Figure that the real derivative of C(t) at the point t is equal to the slope of the tangent line l, while the above approximation formulas correspond to the slopes of the secant lines AC, AB, and BC. Moreover, the approximation errors using the three finite difference methods are the higher-order infinitesimal of Δt, Δt, and (Δt)2, respectively, which means that the smaller the sample interval Δt is, the more accurate the approximation is. However, the error increased by the loss of the valid digits should be also taken into consideration at the same time. Only in terms of the accuracy, the central difference method is more preferable than the other two methods.

Figure 3

Geometrical meanings of the three finite difference approximation methods.

More generally, the approximation formula of finite difference can also be obtained using the Lagrange interpolation polynomial, which is used to approximate the function C(t) (t = t0, t1, ..., t). The Lagrange interpolation polynomials L(t) going through the different n + 1 nodes (t, C(t)) (t = t0, t1, ..., t) are shown in the following formulaAnd then the corresponding derivatives of C(t) at t = t0, t1, ..., t can be approximated by computing the derivatives of L(t) at each node.[44−46] For the special and common cases using two nodes and three nodes, the same formulas as eqs –14 can be obtained. And for the cases using three nodes and five nodes, the following two new formulas (3NFD and 5NFD) can also be derived easily[24,42,43]

Polynomial Fitting

Unlike integral, differential is used to describe the local characteristic of a function and hence is hypersensitive to the small changes of function, which are likely to result in unexpected large fluctuation in derivatives. Especially for the experimental data obtained in industrial applications, the numerical differentiation should be avoided from employing directly. In this case, it is preferable to fit the experimental data with a least-squares curve first and then make the differential operation on the fitted function. The most common method to estimate the discrete numerical derivative is the polynomial fitting.[24,34,47,48] For a set of discrete measurement data (t, C(t)) (t = t0, t1, ..., t), the task of polynomial fitting is to determine the parameters {a}(i = 0, 1, ..., m) of the mth-order polynomial P(t) = a0 + a1t + a2t2 + ··· + at, which can be realized by solving the following optimization problemwhere a = (a0,a1, ..., a)T is the decision vector, namely, the parameter vector of mth-order polynomial P(t). The optimization objective function is the residual sum of squares using the polynomial P(t) to fit the real data C(t) at the points t = t0, t1, ..., t, which can be easily solved by the well-known least-squares method.[24,48] Then, the derivatives of C(t) at the points t = t0, t1, ..., t can be approximated by the corresponding derivative values of the above-determined polynomial function P(t)

To some extent, the polynomial fitting has something common with the Lagrange interpolation polynomial. However, the interpolation problem requires that the approximated curve goes through the given n + 1 data points strictly. When the noise or even error in measurement data is large, the incorrect information behind data will be also learned by mistake (overfitting). Whereas, the polynomial fitting need not to meet the above requirement strictly, but only need to reflect the general trend of the given data as accurately as possible, in which sense the polynomial fitting method plays a role of smoothing the noisy data. In the polynomial fitting method, the polynomial order m is an important parameter to influence the fitting performance. Usually, a simple and practical approach to determine the optimal polynomial order m is to use the statistical test. If there is no statistically significant difference in the variances between two orders m + 1 and m at a significant level, then we do not use order m + 1 but only m to fit the given data.

Savitzky–Golay Filter

As a special case of the polynomial fitting, the Savitzky–Golay (SG) filter is widely used to smooth and denoise the data streams in various engineering applications, which is a digital filter method in time domain based on least-squares fitting of the local polynomials. While the noise is filtered out through the SG filter, the shape and width of the original signal still remain the same.[27,49−51] Different from the traditional polynomial fitting method that fits the overall measurement data to an mth-order polynomial function, the SG filter method conducts the similar polynomial fitter for the data points (2M + 1) in each local window, which are obtained based on a central point and several (M in each side of the central point) nearby points. The schematic diagram of the SG filter (M = 2) is demonstrated in detail in Figure . The width of the filter window is N = 2M + 1. In the Savitzky–Golay filter, the width N = 2M + 1 of the local window is an important parameter to influence the filter performance. To ensure the solvability of the local polynomial fitting problem, the width N need be larger than the local polynomial order k. Similarly to the polynomial fitting method, a simple and practical approach to determine the optimal width N is to use the statistical test. For simplification and consistency, the independent variable t in each window is transformed into z by a linear transformationwhere t̅ is the value of the central point in each window and z is the transformed independent variable, which takes the values −M, −M + 1, ..., −2, −1, 0, 1, 2, ..., M – 1, M. Now in each window, the same kth-order polynomial fitting is conducted and the fitted function S(z) = a0 + a1z + a2z2 + ··· + az is employed to approximate the value and the first-order derivative at the central point[27,51,52]

Figure 4

Schematic diagram of the SG filter.

It should be noted that there are totally 2M points not filtered in any window (M points both in the first and last windows). A simple solution is to approximate the corresponding filtered and derivative values with the fitted polynomial functions in the first and last windows or to employ some smoothing strategy. Generally, the parameters a0, a1, a2, ..., a in S(z) can also be calculated easily by referring to the tables (for various M and k) of convolution coefficients, which are obtained by a simple polynomial convolution method and reported in the reference by Savitzky and Golay.[27] The above way to calculate the parameters a0, a1, a2, ..., a by the convolution coefficient table is very convenient and time-saving.

Wavelet Decomposition

In recent years, the wavelet theory has developed very rapidly because of its better time and frequency domain characteristics compared to the Fourier transformation.[53−56] The wavelet transformation is an ideal tool to analyze and deal with the signal in time and frequency domains, which has been widely applied to various fields, such as signal analysis, data processing, image processing, speech synthesis, quantum mechanics, theoretical physics, and so on.[54−58]

The most general model for the noisy signal has the following formwhere t is the sample time, the superscript M corresponds to the measurement data, and e(t) is the Gaussian white measurement noise. The denoising objective is to suppress the noise part of the signal C(t) and to recover C(t). The aim of wavelet denoising is to distinguish between the original signal and the noisy signal, the essence of which is the function approximation problem in mathematics, that is, how to find the best approximation to the original signal in the function space spanned by the versions of scaling and translating the wavelet mother function according to the proposed criteria. The wavelet denoising method usually includes the following three basic steps:[59−62]

Wavelet decomposition—choose a wavelet and make the wavelet transforming on the noisy signal C(t) at the level N of wavelet decomposition;

Detail coefficients thresholding—perform some processing on the transformed wavelet coefficients to remove the noise, namely, for each level from 1 to N, select a threshold and apply the soft thresholding to the detail coefficients;

Reconstruction—perform the inverse wavelet transformation on the processed wavelet coefficients to obtain the denoised signal, namely, compute wavelet reconstruction based on the original approximation coefficients of level N and the modified detail coefficients of levels from 1 to N.

Finally, the above denoised data are used to estimate the kinetic reaction rates according to the derivative methods.

Tikhonov Regularization

To prevent the measurement noise from amplifying and propagating due to the derivative operation, a novel and effective estimating strategy was developed to deal with the time–concentration data with measurement noise,[23,32,35] namely, the Tikhonov regularization method, in which the rate is not obtained by the direct derivative operation on concentration but by the integration operation on the intermediate variable. Due to the smoothing effect of the integration operation, the noise in the concentration measurements is kept under control as much as possible and hence the estimated rate is less affected by measurement noise.

Assuming the initial time t0 = 0 for convenience, the following integral equation can be easily obtained according to eq

By introducing the intermediate variable f(t) = dr(t)/dt and making some proper transformation (such as integration by parts) on eq , we can obtain the following equation[23]where the superscript “C” stands for the corresponding variables by model computation. By solving the above Volterra integral eq , the unknown intermediate variable f(t), constants C(0) = Cs and r(0) can be obtained easily and then the rate r(t) can be computed by integration of f(t). In the gold cyanidation leaching practice, the grade and concentration measurements are not sampled continuously in view of the production technology, measurement sensors, and cost. Therefore, only the discrete measurements are feasible and eq is transformed into the discretized form[23,32,35]where C is the corresponding value at time t (i = 0, 1, 2, ..., n), j = 1, 2, ..., N, N is the number of uniformly spaced discretization points between 0 and t, the time interval between the adjacent discretized points is Δt = t/(N – 1), f1, f2, ..., f are the f(t) values at the N discretized points, and α (such as determined according to Simpson’s 1/3 rule) is the corresponding coefficient used to approximate the integral part in eq .[23,32,33] To reduce the influence of measurement noise, in the Tikhonov regularization method, the unknown f = [f1, f2, ..., f]T and r(0) are obtained by solving the following optimization problem[23]where the objective function J consists of two parts, namely, the square sum of the deviations between the calculated values and the measurements and the square sum of the second-order derivatives of f at all of the discretization points (except for the two end points), which is used to ensure the smoothness characteristics of f. β is the coefficient matrix used to approximate the second-order derivatives of f based on finite difference technique, which is given bywhere λ is the regularization factor used to compromise the fitting deviation and smoothness, which can be determined by the trial-and-error method, discrepancy principle, L-curve criterion method, and general cross-validation (GCV) method.[23,39,40] Then, the kinetic reaction rate can be easily computed by the integration of f = [f1, f2, ..., f]T by numerical integral method.

Results and Discussion

In this section, the proposed various strategies described in the last section are used to estimate the kinetic reaction rates of gold cyanidation leaching process, namely, the gold solution rate and the cyanide-ion consumption rate. Without loss of generality, only the gold solution rate is investigated in this paper for convenience. Due to the fact that the kinetic reaction rates are unmeasurable in the gold cyanidation leaching plant, the mechanism model presented in Section is used to simulate the actual gold cyanidation leaching process and generate the necessary simulation data. The kinetic model parameters in the above mechanism model have been identified and calibrated based on a set of actual industrial data from an Australian gold processing plant.[18−21] The process model is in good accordance with the actual plant and has better prediction performance in practice, which indicates that it is appropriate and capable of simulating the actual gold cyanidation leaching process accurately. The relevant variable and parameter values used to simulate the reality with the above mechanism model are listed in Table .[18−21] To simulate the actual industrial measuring situation, the Gaussian noise with zero mean and 5% standard deviation are also added to the output concentration measurements (Cs and CCN). The sampled interval is 0.1 h (6 min). It is obvious that if the sample interval is smaller, more process data will be collected and used to estimate the kinetic reaction rates and hence lead to more accurate estimation results, however, which is at the expense of production cost. In the practical applications, the proper sample interval should be determined according to the real characteristics of the concentration measurements. The simulated data curves of gold grade in the ore (Cs) and the corresponding generalized gold grade (C̅s) with and without measurement noise generated by the mechanism model are shown in Figure a,b, respectively. Here, the aim to introduce the generalized gold grade is to reduce the effect of measurement noise on the estimating results as much as possible because the left-hand side of the mechanism model eq still includes the noisy item Cs. The generalized gold grade can be easily obtained by transforming eq where all of the measurement noisy are introduced into the generalized gold grade C̅s. All of the following simulation tests are performed by the MATLAB software package on an Intel i5 3.1GHz PC with 4G RAM. In this paper, the estimating performances of various strategies are evaluated using the following three performance indices, namely, root-mean-square error (RMSE), mean absolute error (MAE), and mean absolute percentage error (MAPE)where rAu, is the ith real value (the nominal measurement value actually simulated and generated by the mechanism model) of the gold solution rate and r̂Au, is the ith corresponding value estimated with some strategy. Num is the data number of the estimating data set.

Table 1

Relevant Variable and Parameter Values Used To Simulate the Reality with the Mechanistic Model

Qs	QCN	Cs0	CCN0	d̅	Cw	V	ρs	ρl	Co
1.645 × 105	3 × 108	5	200	30	48.8%	412	2.8	1.0	8

Figure 5

Simulated data curves with and without measurement noise generated by the mechanism model. (a) The original data curve of the gold grade in the ore. (b) The data curve of the generalized gold grade in the ore.

First, the five finite difference methods described in eqs –14, 16, and 17 (labeled with FFD, BFD, CFD, 3NFD, and 5NFD, respectively) are used to estimate the gold solution rates for the noiseless and noisy measurement data, and the corresponding estimating results are shown in Figure a,b, respectively. Meanwhile, the estimating errors quantized with the above three performance indices are listed in Table for comparison. In can be easily observed from Figure a that when the measurement noise is not absent, the estimating results using the CFD, 3NFD, and 5NFD methods are basically in accordance with the nominal actual kinetic reaction rate values, whereas the estimating errors using the FFD and BFD methods are very obvious, which is in agreement with the theoretical error analysis and primarily due to the fact that in both of the methods only two consecutive measurement data are used to estimate the local derivative. However, for the case with measurement noise, no acceptable estimating results can be obtained using the above five finite difference methods and the oscillation is sharp, which can be obviously observed in Figure b. In this scenario, the CFD and 5NFD methods have the relatively good estimating results than others. Because of the direct derivation characteristic of finite difference methods, the noise in the measurement data is amplified and propagated into the derivative values, namely, the gold solution rates. Hence, the errors between the estimating values and the nominal actual values are extremely large, which is mainly affected by the noise level, sample interval, estimating methods, data size, and so on. If the inaccurate kinetic reaction rate values are used to identify the model parameters, no accordant model with the reality of gold cyanidation leaching process will be obtained. Therefore, the established model is not capable of predicting the key production index accurately and then guiding the production operation.

Figure 6

Corresponding estimating results of the gold solution rates for the noiseless and noisy measurement data using the five finite difference methods: (a) no noise and (b) 5% noise.

Table 2

Estimating Error Quantized with the Three Performance Indices (RMSE, MAE, MAPE)

	RMSE (no noise)	MAE (no noise)	MAPE (no noise)	RMSE (5% noise)	MAE (5% noise)	MAPE (5% noise)
FFD method	0.4242	0.2797	0.0593	1.0967	1.0310	0.3803
BFD method	0.3617	0.2323	0.0488	1.1353	0.9912	0.4048
CFD method	0.0770	0.0416	0.0078	0.4322	0.3423	0.1267
3NFD method	0.0968	0.0578	0.0110	1.8339	1.5833	0.6149
5NFD method	0.0807	0.0282	0.0036	0.5540	0.4242	0.1624
third-order PF method	0.4676	0.3644	0.1220	0.4763	0.3666	0.1247
fourth-order PF method	0.1635	0.1213	0.0399	0.3009	0.1873	0.0708
fifth-order PF method	0.1938	0.0802	0.0130	0.5470	0.2682	0.0683
sixth-order PF method	0.1725	0.0831	0.0176	0.5215	0.3298	0.1332
seventh-order PF method	0.1151	0.0597	0.0142	0.3101	0.2744	0.1055
SG filter method (k = 2, N = 3)	0.0770	0.0416	0.0078	0.4322	0.3423	0.1267
SG filter method (k = 2, N = 5)	0.1402	0.0944	0.0211	0.3112	0.2388	0.0810
SG filter method (k = 2, N = 7)	0.3430	0.2062	0.0446	0.2271	0.1741	0.0481
SG filter method (k = 2, N = 9)	0.5542	0.3432	0.0768	0.5117	0.3478	0.0871
SG filter method (k = 3, N = 5)	0.1089	0.0396	0.0061	0.6224	0.5082	0.1925
SG filter method (k = 3, N = 7)	0.4520	0.1904	0.0278	0.7736	0.5808	0.2133
SG filter method (k = 3, N = 9)	0.8427	0.3807	0.0592	1.1677	0.8358	0.2668
SG filter method (k = 3, N = 11)	1.2128	0.5927	0.1004	1.6728	1.0574	0.2587
WD and CFD method (first level)	0.0911	0.0487	0.0108	0.6068	0.4124	0.1134
WD and CFD method (second levels)	0.0921	0.0694	0.0217	0.2706	0.2161	0.0779
WD and CFD method (third levels)	0.0953	0.0704	0.0229	0.3345	0.2756	0.1111
WD and CFD method (fourth levels)	0.0958	0.0737	0.0246	0.3807	0.3091	0.1194
TIK regularization method (λ = 5 × 10–15, −9)	0.0639	0.0463	0.0125	0.2857	0.1393	0.0371
TIK regularization method (λ = 5 × 10–14, −8)	0.0653	0.0453	0.0118	0.1175	0.0797	0.0204
TIK regularization method (λ = 5 × 10–13, −7)	0.0738	0.0414	0.0090	0.3615	0.2793	0.0940
TIK regularization method (λ = 5 × 10–12, −6)	0.0978	0.0397	0.0062	0.4625	0.3561	0.1210

Second, the polynomial fitting (PF) method is employed to estimate the gold solution rate based on the generalized gold grade measurements. Instead of conducting the direct derivative operation on the noisy measurements, a polynomial function to fit the above measurements best is obtained and then the corresponding derivative values at the sample points are calculated according to the obtained polynomial function. In this scenario, the orders of the polynomial functions are selected as 3, 4, 5, 6, and 7, respectively. The polynomial fitting function curves of the generalized gold grades for the noiseless and noisy measurement data using the polynomial fitting method are shown in Figure a,b, respectively. It can be easily observed that for the noiseless case, the polynomial fitting functions with order greater than 3 have good fitting effect, and it is hard to distinguish the difference between the real values and the fitting values with the naked eye. However, for the noisy case, the variation trends of the polynomial fitting functions are identical with the generalized real noisy measurements and the fitting values at each sample time cannot be guaranteed to be equal to the real values due to the characteristic of the polynomial fitting method, which is the advantage of avoiding overfitting the noise. The fitted polynomial functions of third order to seventh order for the noisy case are easily obtained as follows (rounded off to five significant digits)Now the corresponding derivative functions of eq –37 can be easily calculated according to the polynomial derivative rules, and then, the estimated values of the gold solution rate at the sample time (t0, t1, ..., t) can be easily obtained. The corresponding estimating results of the gold solution rates for the noiseless and noisy measurement data using the polynomial fitting method are shown in Figure a,b, respectively. Meanwhile, the estimating errors quantized with the above three performance indices are also listed in Table for comparison. It can be observed from Figure and Table that for the noiseless case, the advantage of the polynomial fitting method is not very obvious and the estimating results are even slightly worse than those with the direct finite difference methods. However, for the noisy case, the superiority of the polynomial fitting method over the finite difference method is obvious that the oscillation of the estimating results because of measurement noise reduces and the three performance indices are improved a lot (such as the fourth-order PF method and the CFD method). Due to the fact that the polynomial fitting operation before the derivative operation has some smoothing effect, the derivative estimation using the preprocessed data (polynomial function values at the sample time in fact) is less affected by the measurement noise compared to the finite difference method. Moreover, the estimating curve oscillates more violently with the increase of the order in the polynomial function, which is due to the fact that the higher-order polynomial function is more complicated in structure and has enough abilities to fit the generalized gold grade data with measurement noise. At the same time, the measurement noise in the measurement data is also learned and fitted by the high-order polynomial function, which is not expected. Therefore, the fitting error and the smoothness should be considered together when determining the order in the polynomial function.

Figure 7

Polynomial fitting function curves of the generalized gold grades for the noiseless and noisy measurement data using the polynomial fitting method: (a) no noise; (b) 5% noise.

Figure 8

Corresponding estimating results of the gold solution rates for the noiseless and noisy measurement data using the polynomial fitting method: (a) no noise; (b) 5% noise.

As a special case of the polynomial fitting method, the Savitzky–Golay (SG) filter can be considered as a piecewise polynomial fitting method. In this scenario, the SG filter is used to estimate the gold solution rates based on the generalized gold grade measurements. Because the parameters such as the order k in the local polynomial function and the window width N = 2M + 1 have a significant impact on the filter effect, here in this paper, the two special cases k = 2, N = 3, 5, 7, 9 and k = 3, N = 5, 7, 9, 11 are considered in view of the characteristic of fitting the generalized measurement values locally. And here the total 2M values not filtered in the first and last windows are estimated based on its nearby data by the simple mean filter, namely, r(i) = 2r(i + 1) – r(i + 2) (i = M – 1, M – 2, ..., 0) and r(i) = 2r(i – 1) – r(i – 2) (i = n – M + 2, n – M + 3, ..., n + 1, where the data number is n + 1). The filtered curves of the generalized gold grades for the noiseless and noisy measurement data using the SG filter method are shown in Figure a–d for different parameters. And the corresponding estimating results of the gold solution rates for the noiseless and noisy measurement data using the SG filter method are shown in Figure a–d for different parameters, respectively. Meanwhile the estimating errors quantized with the above three performance indices are also listed in Table for comparison. It can be easily observed from Figure that for all of the cases, the trend of the filtered gold grades is correct, accurate, and there is less difference between the filtered and the actual values. However, for the gold solution rates, it is not the same and the slight difference between the filtered and actual values results in a large deviation on the rate estimation due to the derivative operation, which can be easily observed from Figure and Table . Here, the sample time is 0.1 h and the noise in measurement data is likely to be amplified about 10 times. In all of the noisy cases, the estimating results using the parameters k = 2 and N = 7 are the best and even better than the noiseless case, which can be attributed to the simple handling method for the 2M unfiltered values in both ends. Moreover, because the second-order polynomial function is capable of fitting the local data characteristic, the higher third-order polynomial function will fit the noisy measurements over and hence lead to bad estimating results. In general, the SG filter method is better than the above FD and PF methods if the proper parameters are selected.

Figure 9

Filtered curves of the generalized gold grades for the noiseless and noisy measurement data using the SG filter method: (a) no noise (second order); (b) 5% noise (second order); (c) no noise (third order); and (d) 5% noise (third order).

Figure 10

Corresponding estimating results of the gold solution rates for the noiseless and noisy measurement data using the SG filter method: (a) no noise (second order); (b) 5% noise (second order); (c) no noise (third order); and (d) 5% noise (third order).

Next, the wavelet decomposition (WD) method is used to denoise the generalized gold grade data and then estimate the gold solution rates with the CFD method in this scenario. Here, the Wavelet Toolbox in the soft of Matlab R2016b is used to impletment the wavelet decomposition operation for convenience.[62] In detail, the Daubechies wavelet is selected as the wavelet base, which is a family of orthogonal wavelets defining a discrete wavelet transform and characterized by a maximal number of vanishing moments for some given support. The wavelet decompositions are performed at levels N = 1, 2, 3, and 4. The soft thresholding is selected and the threshold selection rule is heursure, namely, a heuristic variant of the principle of Stein’s unbiased risk.[63,64] Moreover, the threshold is rescaled using a single estimation of level noise based on first-level coefficients. The denoised curves of the generalized gold grades for the noiseless and noisy measurement data using the WD method are shown in Figure a,b, respectively. And the corresponding estimating results of the gold solution rates for the noiseless and noisy measurement data using the CFD method are shown in Figure a,b. Meanwhile the estimating errors quantized with the above three performance indices are also listed in Table for comparison. It can be easily observed from Figures and 12 and Table that for the noiseless case, the denoised curves using the WD method are almost identical with the original data and hence the corresponding estimating results based on the above denoised data using the CFD method are expected to be as good as the above case using only the CFD method. However, for the noisy case, due to the added preprocessed operation on the noisy data, the noisy curves oscillate less obviously and turn to be smoother, which leads to better estimating results than the case estimating the derivatives directly using the CFD method. In general, the wavelet decomposition is employed as an effective tool to preprocess and smooth the noisy data.

Figure 11

Denoised curves of the generalized gold grades for the noiseless and noisy measurement data using the WD method: (a) no noise; (b) 5% noise.

Figure 12

Corresponding estimating results of the gold solution rates for the noiseless and noisy measurement data using the WD and CFD method: (a) no noise; (b) 5% noise.

Finally, the novel Tikhonov (TIK) regularization method is used to deal with the generalized gold grade data and then estimate the gold solution rates in this scenario. The corresponding estimating results of the gold solution rates for the noiseless and noisy measurement data using the TIK regularization method are shown in Figure a,b. And the backcalculated curves of the generalized gold grades for the noiseless and noisy measurement data using the TIK regularization method are shown in Figure a,b, respectively. Meanwhile, the estimating errors quantized with the above three performance indices are also listed in Table for comparison. It can be easily observed from Figure and Table that the satisfactory estimating results are obtained in both the noiseless and noisy cases, which is associated with the essential idea of the TIK regularization method that the rate is not obtained by the direct derivative operation on concentration data but by the integration operation on the intermediate variable. For the noiseless case, the objective of the estimating task is the fitting errors between the calculated values and the measurements, and hence, a relative smaller regularization factor should be selected, where λ = 5 × 10–15, 5 × 10–14, 5 × 10–13, and 5 × 10–12 are simulated in this scenario, respectively. However, for the noisy case, the objective of the estimating task emphasizes on the smooth performance of the introduced intermediate variable, and hence, a relatively greater regularization factor should be selected, where λ = 5 × 10–9, 5 × 10–8, 5 × 10–7, and 5 × 10–6 are simulated in this scenario, respectively. It is also obvious in Figure that with the continuous increase of λ, the estimation results turn out to be inaccurate increasingly. The best choice of λ simulated in the noisy case is 5 × 10–8, which leads to the good and acceptable performance indices of RMSE = 0.1175, MAE = 0.0797, and MAPE = 0.0204. The key of the TIK regularization method is to determine a proper λ that considers both the fitting errors and the smooth performance simultaneously. In practice, we can compare the backcalculated curves with the plant measurements and then select a proper λ. And other approaches include the discrepancy principle method, the L-curve criterion method, and the general cross-validation method.

Figure 13

Corresponding estimating results of the gold solution rates for the noiseless and noisy measurement data using the TIK regularization method: (a) no noise; (b) 5% noise.

Figure 14

Backcalculated curves of the generalized gold grades for the noiseless and noisy measurement data using the TIK regularization method: (a) no noise; (b) 5% noise.

In conclusion, the above strategies have their own advantages and disadvantages. In practice, the proper strategy should be selected according to the actual situations investigated, such as the level of the measurement noise, the sample time, the data number, the complicated characteristic of the process, and so on. Furthermore, two or more of the above strategies can be also used together, which is expected to give rise to more perfect estimating results and is beyond the scope of this paper.

Conclusions

It is very important to establish an accurate process model for implementing further control and optimization for gold cyanidation leaching process. Unfortunately, the important kinetic reaction rates affecting process operation are unmeasurable, and moreover, their estimation is an ill-posed inverse problem due to the fact that the noise in concentration measurements is easy to be amplified and propagated into the rate estimates by derivative operation. In this paper, the alternative strategies (finite difference, polynomial fitting, Savitzky–Golay filter, wavelet decomposition, and Tikhonov regularization) estimating the kinetic reaction rate of gold cyanidation leaching process are investigated in detail. The simulation results show that the direct finite difference leads to poor estimating results for the noisy case and the other strategies are capable of avoiding the noise amplification and improving the estimating results in some extent. In all of the investigated strategies, the Tikhonov regularization leads to the satisfactory and acceptable estimating results in both noiseless and noisy cases, which will lay an important foundation for the subsequent model identification, production index prediction, and operation optimization.

Furthermore, the above strategies have their own advantages and disadvantages. In practice, a proper strategy should be selected according to the actual situations investigated, such as the level of the measurement noise, the sample time, the complicated characteristic of the process, and so on. It should be noted that the above estimating strategies can also be applied to other industrial processes with similar characteristics. In the future, some practical implementation issues may require further investigation. Moreover, two or more of the above strategies can be also used together, which is expected to give rise to more perfect estimating results.