Moving Horizon Estimation of Sulfur Concentrate Grade Based on Kinetic Models Under Multiple Working Conditions

This article proposed a moving horizon approach to predict sulfur concentrate grade in a sulfur flotation process. In this approach, the grade prediction problem is formulated as a moving horizon estimation problem. To build the nominal model for the estimation, the kinetic of sulfur flotation process are firstly studied. The unknown variable in the kinetic model, i.e., the flux of overflow, is obtained using machine vision technology. Moreover, to account for the multi-modal characteristic of sulfur flotation process, the kinetic model parameters are identified for different working conditions. To eliminate the effect of, a robust parameter identification approach is adopted. Finally, the kinetic model is embedded in the moving horizon estimation framework, which reconstructs the concentrate grade using from online measured process variables. Experimental results demonstrate the feasibility and efficient of the proposed method.


I. INTRODUCTION
Sulphur concentrate grade is a key index for metallurgical enterprises to evaluate the flotation performance [1]. Due to the lack of an economical online detection device, the sulfur concentrate grade is conventionally obtained by discrete-delayed chemical test. The test process is complicated and time-consuming. The human operators cannot obtain the sulfur concentrate grade online. Therefore, one can only adjust the manipulated variables based on the last testing result and operational experience [2], [3]. In addition, the complexity of process dynamics and the time-varying characteristics of feeding minerals increased the difficulty of stable and optimal operation, resulting in low concentrate grade and severe fluctuations, which is a waste of valuable resources and decrease of production efficiency. As the foundation of operational optimization, an online estimation approach that provides accurate estimate of sulfur concentrate grade and supports operational decision is essential in increasing the production efficiency [4].
The associate editor coordinating the review of this manuscript and approving it for publication was F. K. Wang . The concentrate grade estimation problem has been a hot topic of mineral flotation processes. On one hand, the flotation kinetics is the foundation for concentrate grade estimation. It has been extensively studied since the middle of last century. Various types of models concerning both the dynamics of flotation process and the sub-processes of flotation have been proposed, e.g., flotation rate model [5], mass balance model [6], attachment time model [7], etc. A static simulation of characteristic parameters (suspended minerals, mineralized froth) based on the mass balance equation was proposed to predict the concentrate grade and flotation recovery rate [8]. Based on the hydrodynamic principles in the flotation cell, a mechanism model of concentrate grade and froth recovery rate was established [9], [10]. In [11] and [12], the flotation rate equation and the principle of mass balance were combined to establish the mechanism model of concentrate grade and flotation recovery rate. In [13], based on a laboratory study, the relationship model between flotation recovery rate and concentrate grade was obtained to realize the measurement of recovery rate and concentrate grade. However, most of these studies are based on certain assumptions that some key process variables can be obtained online. For example, the pulp level is difficult to obtain due to lack of reliable and easy-to-maintain measurement devices. In order to obtain the pulp level, some experimental studies are conducted in a transparent flotation cell. However, in practice, the flotation cell is not transparent. On the other hand, as froth image features can reflect flotation performance, a lot of studies on the concentrate grade estimation based on machine learning have been carried out [14], [15]. M. Massinaei et al. proposed a machine-learning based approach for analysis of performance parameters of coal flotation process [16]. Liu et al. proposed a method for extracting froth color and texture based on multi-resolution and multivariate, and achieved the measure of concentrate grade by partial least squares [17]. Hargrave et al. analyzed the relationship between froth characteristics and flotation performance, and established a relationship model between concentrate grade and froth characteristics [18]. Marais and Aldrich established a multi-layer neural network model based on froth texture characteristics to predict platinum concentrate grade and recovery rate [19]. Kaartinen et al. used principal component analysis to analyze the relationship between zinc flotation image feature and concentrate grade, and a model for estimating the concentrate grade based on partial least squares was established to solve the problem that the concentrate grade cannot be measured continuously online by X-ray fluorescence analysis [20]. Bergh et al. analyzed the characteristics of froth images according to multivariate statistical methods, diagnosed the fault conditions of the flotation process online, and established a grade prediction model based on partial least squares, based on which a predictive control method is proposed to control the flotation process [21]. The above results provides experiences in constructing kinetic model and froth image feature based concentrate estimation model. However, some practical issues in mineral flotation processes are not considered in the above studies, e.g, online measure of pulp level [22], multiple working conditions. Therefore, a hypothesis is proposed: can the flotation kinetics and froth image features be combined in a framework to provide accurate estimation of concentrate grade?
For practical sulfur flotation process, due to the high temperature, acid mist and strong corrosiveness of the sulfur flotation site, it is difficult to use the commonly physical liquid level gauges in this harsh environment for a long time. Therefore, the overflow concentrate slurry flow and froth depth of the sulfur flotation cannot be measured online. In addition, the sulfur flotation process has multiple working conditions, under which the kinetic model parameters take different values. Therefore, in this article, through the analysis of the characteristics of different working conditions, kinetic models of the sulfur flotation process under multiple working conditions are established first. The features of froth image are then used to predict the froth depth and flow rate of overflow concentrate slurry. Finally, the established kinetic model is used as the nominal model of a moving horizon estimation framework, which is designed to achieve online measurement of the sulfur concentrate grade.
Actually, the concentrate grade estimation problem also exist in other continuous industrial processes. Therefore, the resulted moving horizon estimation framework provides an alternative solution for the estimation of key performance index of complex industrial processes with multiple working conditions and unmeasurable key process variables.
The rest of the paper is organized as follows. Section II presents the problem analysis. Section III gives a detailed description of kinetic modeling of sulfur flotation process. Section IV presents the kinetic model parameter identification under multiple working conditions. Section V describes the online estimation method of sulfur concentrate grade. Experimental results are discussed in section VI. Section VII concludes with future areas for research.

II. PROBLEM ANALYSIS
The sulfur flotation is a physical process, which mainly uses the natural hydrophobicity of elemental sulfur, blowing air in the flotation tank, and generating bubbles through the stirring of the flotation machine. The sulfur element with strong hydrophobicity will be attached to the surface of the bubbles and float up to form a froth layer. Finally, sulfur concentrate products are obtained through overflow. Since no flotation reagent is added during the entire sulfur flotation process, the flotation conditions are adjusted mainly by adjusting the height of the ore slurry surface while controlling the air flow to maintain the uniform and proper froth depth. The level of sulfur concentrate is an important indicator of sulfur flotation performance, and also determines the sales price of sulfur concentrate products. In order to maximize economic benefits, enterprises must maximize the concentrate grade of sulfur. Therefore, the concentrate grade of sulfur is an important indicator of the effect for sulfur flotation.
At present, the level of automation in the domestic flotation production process is generally not high, and there is no effective on-line detection equipment for the sulfur concentrate grade. The sulfur concentrate grade is mainly obtained by the on-site test of workers on site. However, the sulfur test is different from the general metal grade test, which requires complicated processes such as filtration, drying and crushing. The entire test process is complicated and takes a long time. Therefore, the current value of concentrate grade cannot be obtained in real time at the sulfur flotation site. This severely lagging detection method of sulfur concentrate grade makes the operator unable to adjust the working conditions in time, which makes it difficult to achieve the optimal value of the sulfur concentrate grade. Therefore, how to measure the concentrate grade of sulfur online under the existing conditions is a difficult problem to be solved in the process of sulfur flotation.
The concentrate grade of sulfur is determined by the content of sulfur and impurities in the sulfur concentrate. According to the principle of sulfur flotation process kinetics [23], the kinetic models of sulfur flotation process based on the mass balance include the composition of sulfur and impurities in the slurry, the flotation rate and other information in flotation cell. Therefore, the kinetic models can be used to obtain the information required to calculate the concentrate grade of sulfur. However, due to the rich dynamic characteristics of the sulfur flotation process, a single kinetic model cannot fully describe the sulfur flotation process. Therefore, it is necessary to classify the operating conditions of the sulfur flotation and obtain the kinetic models under different operating conditions. In addition, some parameters in the kinetic models cannot be obtained in real time, such as the flow rate of overflow concentrate slurry and the froth depth, leading to the difficulties for obtaining the kinetic models of the sulfur flotation process. Due to the time-varying characteristics of the sulfur flotation process, the parameters of the kinetic model are not invariable under various operating conditions. Therefore, it is necessary to select a suitable online measurement method to overcome the measurement deviation caused by the time-varying characteristics of the sulfur flotation process, so as to obtain the accurate measurement of sulfur concentrate grade.
The visual characteristics of the flotation surface froth contain a wealth of information describing the flotation state. The characteristics of each image are strongly correlated with the sulfur flotation process conditions and the froth depth, so that digital image processing technology can be introduced into the flotation condition recognition, the measurement of froth depth and overflow concentrate slurry flow. Through continuous online acquisition of froth characteristics, froth features based working condition recognition model and online measurement model of froth depth are established, seen in literature [22], [24] for details, which provides real-time reliable process parameter and working condition information for the online measurement of concentrate grade in sulfur flotation process.

A. KINETIC MODELING OF SULFUR FLOTATION PROCESS
The establishment of the kinetic models of the sulfur flotation process is the basic condition and the first step to realize the online measurement of sulfur concentrate grade. In this section, the kinetic models of sulfur flotation process are firstly established from the perspective of mass balance. For the problem that the overflow flow in the kinetic model cannot be directly obtained, the overflow flow information is obtained through machine vision. The relationship between the flotation rate and the froth depth was obtained through analysis. Finally, the kinetic models can be obtained which can be used to calculate the grade of sulfur concentrate.

1) KINETIC MODELS OF FLOTATION PROCESS BASED ON MASS BALANCE
Consider a single flotation cell (Figure 1), if the pulp is divided into two types of minerals, i.e., sulfur and impurities, it can be known from mass balance principle that the dynamic characteristics of sulfur flotation process can be described by the kinetic equation (1), where M 1 is the mass of sulfur in the flotation cell, and M 2 is the mass of impurities in the flotation cell, and satisfies: where M is the total mass of the pulp in the flotation cell. K 1 and K 2 are the flotation rates of sulfur and impurities, respectively. Q T is tailing flow rate, Q C is the flow rate of the overflow concentrate slurry, Q F is the flow rate of the ore, ε g is the air occupation rate, G s1 is the mass of sulfur in the ore per unit time, G s2 is the mass of impurities in the ore per unit time, S c is the cross-sectional area of the flotation cell, h is the slurry level (h = H -h F , H is the height of the flotation cell, h F is the froth depth). The concentrate grade of the sulfur flotation cell can be calculated by the following equation: where G SC1 = K 1 M 1 is the mass of sulfur in the produced concentrate per unit time, G SC2 = K 2 M 2 is the mass of impurities in the produced concentrate per unit time. If consider the dynamics of M , then where If define where 0 < α < 1 is the proportion of sulfur in the flotation cell, then, Compared with Eqs (1) and (3), equation (7) has fewer unknown variables, which is more suitable for model parameter identification. Eqs (1) and (3) directly represent sulfur and impurities in pulp, which is more suitable for qualitative analysis of sulfur flotation process. The detailed information of each parameter in the sulfur flotation kinetic models are shown in Table 1.
In equation (7), ε g and S c are constants, L C is detected every 8 hours, Q T is a manipulated variable which can be measured online, Q C and h can be obtained by machine vision (see Section III-B and III-C), then Q F can be obtained given the value of h. M can be calculated by measuring the specific gravity according to the following equation: where ρ(s) is the pulp density, s is a measure of distance from the bottom of the flotation cell. Since the density of the pulp from the entrance to the last flotation cell does not change much, and the density of the pulp at the entrance can be calculated from the continuous online measurement of the hydrometer, ρ(s) can be replaced by the density of the pulp at the entrance. Similarly, Gs can also be calculated based on Q F and specific gravity.
Most of the variables in equation (7) can be obtained, e.g, ε g , S c , Q T , Q C , h, Q F , M , G s . Therefore, the concentrate grade can be calculated if the values of K 1 , K 2 , and α are obtained.

2) FLOW RATE ESTIMATION OF THE OVERFLOW CONCENTRATE BASED ON FROTH IMAGES
According to Neethling et al. [25], most of the pulp in flotation froth is contained within the Plateau boundary, i.e., the connection network between the bubbles. The length λ of the Plateo boundary per unit volume of froth and the crosssectional area A p of the Plateo boundary of the froth surface can be calculated as: where d f is the diameter of the froth, β is the air recovery rate, η is the viscosity of the pulp, ρ is the density of the pulp, and g is the gravity constant, Q a is the air flow rate. Therefore, the volume percentage of pulp on the surface of the froth is The flow rate of the overflow concentrate froth is: where v is the overflow speed of the froth, h overflow is the thickness of the froth layer at the overflow port, l is the total length of the overflow groove. VOLUME 8, 2020 The air recovery rate can be calculated by the following equation: Therefore, the flow rate of overflow concentrate is 3

) RELATIONSHIP BETWEEN FLOTATION RATE AND FROTH DEPTH
It can be known from equation (7) that K 1 , K 2 and α directly affect the grade of sulfur concentrate: where K 1 and K 2 are the flotation rates of sulfur and impurities, respectively, α is related to the property of the feeding ore. If define, then, Therefore, increasing K 1 the flotation rate of sulfur (or reducing K 2 the flotation rate of impurities) can increase the concentrate grade, while decreasing K 1 (increasing K 2 ) will cause the concentrate grade to decrease. The relationship between the flotation rate K , the ore grade f , the froth depth h F , and the air flow rate Q a is obtained in [26]: The flotation rate of valuable minerals is larger than the flotation rate of impurities. For valuable minerals, increasing the froth depth can increase its flotation rate; however, for impurities, increasing the froth depth will reduce its flotation rate ( Figure 2 shows the relationship between the flotation rate and the froth depth under different ore conditions and air flow rates [26], seen in Table 2). This agrees with the qualitative conclusions obtained through experiments in [27] and [28]: when the froth layer is thin, it is difficult for the bubbles to merge during the upward floating process, and the pulp is easy to overflow with the froth, resulting in a low concentrate grade; when the froth layer is very thick, mergers occur between the bubbles, and the poorly hydrophobic gangue falls back to the pulp with the merger of the froth, resulting in that the concentrate grade is improved.  Therefore, for the sulfur flotation process, following equation is used to represent the relationship between the flotation rate, froth depth and the air flow rate [26]: where a i (i = 1,2,3,4) are coefficients to be identified using production data.

B. KINETIC MODELS PARAMETER IDENTIFICATION UNDER MULTIPLE WORKING CONDITIONS
Due to the influence of feeding and reaction conditions, the sulfur flotation process exhibits complex dynamic characteristics and has various working conditions. One set of model parameters is often difficult to describe the process dynamics under multiple working conditions. Therefore, in order to establish kinetic models that can cover multiple working conditions, this section analyzes the industrial data and identifies the model parameters for different working conditions.

1) RELATIONSHIP BETWEEN FLOTATION RATE AND FROTH DEPTH
By combining Eq (1), (3), and (19), following equation can be obtained: where 0 < M 1 < 9.5, 0 < M 2 < 9.5, 0 < h < 3.9, 0 < Q F < 40, 0 < Q C < 20, 0 < Q T < 20, 0 < Q a < 31, 0 < G s1 < 76, 0 < G s2 < 76, 0 < ε g < 0.08, S c = 3.46185, a 1i (i = 1,2,3,4) and a 2i (i = 1,2,3,4) are the coefficients for K 1 and K 2 , respectively. Due to the continuous changes in the feeding conditions and operational parameters, the sulfur flotation process switches between different working conditions. If the feeding ore flow rate suddenly increase or decrease, and the tailing flow is adjusted reasonably, the fluctuation of the overflow rate is not large, the froth overflow rate and froth depth are relatively stable, leading to that the fluctuation of sulfur concentrate grade is small. If the feeding ore flow rate suddenly increases and the tailing flow is not properly adjusted, the overflow rate will increase accordingly, causing the froth overflow to accelerate. Thus, the froth depth is decreased, and the sulfur concentrate grade is reduced. If the feeding ore flow rate is suddenly decreased, and the tailing flow is not adjusted, the overflow flow will be reduced accordingly, causing the froth overflow to slow down and the flotation level to decrease, resulting in no froth overflow in severe cases. Figures 3 to 18 illustrate how the sulfur concentrate grade changes with the feeding ore grade and the manipulated variables when the model parameters in (20) are in different ranges. According to the sulfur concentrate grade, the flotation process can be classified into 6 working conditions, see Table 3.

2) ROBUST IDENTIFICATION OF KINETIC MODEL PARAMETERS
Consider sulfur flotation process, substitute equation (19) into equation (7), where     where θ = (a 11 , a 12 , a 13 , a 14 , a 21 , a 22 , a 23 , a 24 , α) is the model parameter set to be identified, β = (Q a , h F ) is the air flow rate and froth depth. For each working condition, the identification of kinetic model parameters is to find an optimal set of model parameters θ, which minimizes the error between the concentrate    grade calculated using the kinetic model and the actual value: where θ is the model parameters to be identified, N is the number of data samples for identification, y i is the actual   sulfur concentrate grade, x i is the manipulated variables of the sulfur flotation process (i.e., the air flow rate and the froth depth), σ 2 i is the variance of y i , and f (x i , θ) is the sulfur concentrate grade calculated using the kinetic models.
Robust and precise identification of model parameters is the foundation of control and other model-based control applications [29]- [33]. However, the parameter identification of the sulfur flotation process kinetic models encounters problems including parameter correlation, unbalanced magnitude, outliers and multi-rate. Before using the optimization algorithm to identify the model parameters, the kinetic models and the identification samples are processed as follows, (i) Reconstruction of kinetic models: According to the analysis in Section III, M , Q T , h and G s in the differential equation (21) can be obtained online, ε g and S c are fixed values. Therefore, g(θ, β) can be calculated online. If a ij (i = 1,2; j = 1,2,3,4), α and other kinetic model parameters in equation (22) can be identified, sulfur concentrate grade can be calculated in real time. However, because a ij    (i = 1,2; j = 1,2,3,4) and α are correlated, the identification result is affected. It is necessary to reconstruct equation (22), and the principle of reconstruction is to reduce the correlation between model parameters. If design following new model parameters,â then where σ Q a , σ h F and σ Q a h F are the standard deviations of Q a , h F and Q a h F , respectively, which are used to balance the difference in magnitude between the parameters. The numerical stability of Eqs. (22) and (25) can be evaluated using condition numbers. The condition number is used to measure the degree to which the input change can influence the output, i.e., the sensitivity of the output to the input error. If the condition number is small, the error of the model parameter has little effect on the model output. If the number of conditions is large, small errors in the model parameters would cause large model output errors.
For a nonlinear function f (x) with multiple inputs, the condition number v c (x) can be calculated using the following equation, which is the supremum of the ratio between the increment of f (x) and the increment of x when x approaches infinity.
(ii) Outlier eliminating: As in (i), the accuracy of parameter identification can be improved by eliminating outliers. The removal of outliers can be achieved not only by outlier removal methods, but also by setting reasonable optimization goals. In this article, the outlier is eliminated by converting the optimization goal of the identification process from minimizing the sum of all sample errors to minimizing the median of sample errors.
(iii) Design of identification steps: By reconstruction in (i), the number of parameters is reduced. However,â 1 ,â 2 ,â 3 andâ 4 cannot be used directly to calculate sulfur concentrate grade. Therefore, it is necessary to design the identification steps to obtain the real model parameters a ij (i = 1, 2; j = 1, 2, 3, 4) and α under different working conditions. According to equation (21), the denominator of L C equals to g(θ, β), which can be obtained online. As it can be tested routinely, the real value of K 1 α can be obtained on specific time points. If denote, then According to the reconstruction principle, equation (29) can be transferred to following equation Therefore, the complete parameter identification procedure is: a. Construct identification samples based on the test data, and ensure that these time points include all types of working conditions; b. Calculate g(θ, β) of each sample using equation (21) and real-time data collected, where M is calculated from formula (8), Q T is measured by a flow meter, and Q C is obtained based on froth image and equation (14), and Gs is calculated using the feeding flow rate Q F and the specific gravity of the feeding slurry; c. Use (25) to identifyâ 1 ,â 2 ,â 3 andâ 4 under different working conditions. Set the optimization target to minimizing the median of the identification error to remove the outliers. After removing the outliers, re-identify the model parameters; d. Calculate the value of K 1 α according to the sulfur concentrate grade obtained from the test, identify α and a 1j (j = 1, 2, 3, 4) according to equation (30), calculate a 2j (j = 1, 2, 3, 4) usingâ 1 ,â 2 ,â 3 andâ 4 . In order to compare the identification performance with and without the robust identification method, 350 data samples were constructed for experiment. These data samples are from the production site and include all working conditions from C1 to C5 (As there is no concentrate overflow under working condition C6, and the concentrate grade cannot be detected, which is an abnormal working condition, so it is not considered in the experiment study), of which 300 groups are used as identification samples (60 samples of each type of working conditions), and the remaining 50 groups are used as testing samples (the 50 groups of samples include all the working conditions from C1 to C5). The basic algorithm for parameter identification is particle swarm optimization (PSO).
The results using and without using the robust identification method are shown in Figure 19 and Table 4. It can be seen that no matter the robust identification method is used or not, the kinetic models identified by the two methods can  achieve satisfying estimation of the variation trend of sulfur concentrate grade.
Compared with the method without robust identification, the accuracy of the robust identification method is not always higher on all test samples. The number of testing samples with large errors is reduced. The average relative error decreased from 8.63% to 5.79%, and the maximum relative error decreased from 23.57% to 14.54%. Therefore, by adopting the robust parameter identification method including model reconstruction, outlier eliminating, and step-by-step identification, the identification of nonlinear functions can be reduced to the identification of linear functions, which improves the robustness of identification process and ultimately improves the accuracy of the identification results.

C. ONLINE ESTIMATION OF SULFUR CONCENTRATE GRADE
In the previous section, the kinetic models of sulfur flotation process under different operating conditions were obtained through offline identification of model parameters. These models can accurately describe the process dynamics under different operating conditions. Although the sulfur flotation process switches between different working conditions, this switching is not completed instantaneously. Therefore, in the production process, the process dynamics is not changed abruptly. However, it changes continuously. The parameters in the kinetic models change slowly or rapidly. In addition, external disturbances will also affect the dynamic characteristics of the sulfur flotation process. Therefore, the kinetic models of sulfur flotation process should be time-varying models. In order to realize the on-line detection of sulfur concentrate grade, based on the kinetic models under different working conditions established in the previous section, it is necessary to study the strategy of online estimation of sulfur concentrate grades.

1) PROBLEM REFORMULATION OF REALTIME PREDICTION OF SULFUR CONCENTRATE GRADE
According to the expression of sulfur concentrate grade, the real-time measurement of sulfur concentrate grade relies on the values of process variables such as K 1 , K 2 and α. From the analysis in the previous section, it is known that in the sulfur flotation production process, the air flow rate Q a , the froth depth h F can be measured online, and Therefore, if the values of a 1j (j = 1, 2, 3, 4) and α can be obtained online, the real-time measurement of sulfur concentrate grade can be achieved. In Section IV, although the values of a ij (i = 1, 2; j = 1, 2, 3, 4) and α were obtained through identification, a 1j (j = 1, 2, 3, 4) and α change with feeding conditions and operational variables. In addition, the identification of a 1j (j = 1, 2, 3, 4) and α requires the testing result of sulfur concentrate grade which can only be obtained offline. Therefore, it is necessary to design an applicable online measurement strategy of sulfur concentrate grade. As discussed in the previous section, the denominator K 1 α + K 2 (1 − α) of the concentrate grade L c can be obtained online. If denote K 1 α + K 2 (1-α) as y, the numerator K 1 α of L c as x1, the denominator K 2 (1-α) of L c as x2, a sulfur concentrate grade estimation system can be constructed, and L C = x 1 x 1 +x 2 . If consider the impact of feeding conditions and operating variables on x 1 and x 2 ,ẋ where x = (x 1 , x 2 ) is the system state, u = (Q a , h F ) is the control input, θ is the model parameters (identified in Section IV), y is the system output,f (x, u, θ ) is a function related to feeding conditions and working conditions, i.e., f (x, u, θ ) is different under different working conditions. Discretizing equation (33), where the subscript k indicates the variable at time k. equation (34) is the nominal discrete model of equation (32).
In this article, the identified kinetic models under different working conditions are used as the nominal models of sulfur under different working conditions, If consider the parameter uncertainty, external disturbance and detection error, then equation (34) can be expressed as, where w k is the unknown model uncertainty and external disturbances, v k is the measurement error of y. Based on equation (36), the real-time prediction problem of sulfur concentrate grade can be regarded as an estimation problem, i.e., the estimation of the system state x using the system output y. Because the sulfur flotation process has various working conditions, the process dynamic have timevarying characteristics. Moving horizon estimation approach can handle constraints and perform moving horizon optimization in the time domain. Therefore, this article apply the moving horizon estimation method in the online measurement of sulfur concentrate grade.

2) ONLINE MEASUREMENT OF SULFUR CONCENTRATE GRADE BASED ON MOVING HORIZON ESTIMATION
The state estimation problem is to reconstruct the value of system state {x k } T k=0 from the system output {y k } T k=0 . In equation (36), f k (·) is known. Therefore, estimating the system state {x k } T k=0 equals to obtaining the initial system state x 0 and external disturbance sequence {w k } T −1 k=0 , which can be formulated as following problem, * where y(k; x 0 ; 0; {w j }) is the system output at time k with initial state x 0 and external disturbance w j k j=0 , L k is the stage cost function, L k : W k × V k ∈ R ≥0 for k >0, is the actual value of system state at time 0, (·) is usually expressed as (40), shown at the bottom of the next page where 0 is a symmetric positive definite matrix.
At time T , the solution is (x 0|T −1 , {ŵ k|T −1 } T −1 k=0 ). According to equation (39), since all the output information from time 0 to time T -1 is used, the estimation problem is a full information estimation problem. Its calculation complexity is at least linear with T . As T increases, the load of calculation will gradually increase. Therefore, for nonlinear systems, it is impractical to solve the estimation problem (37) online. It is required to limit the problem size. It is not difficult to find that equation (37) can be written as The estimation problem equation (37) can be transferred to following, * is the arrival cost from x 0 to z at time T-N, T −N is the reachable set at time T-N, However, for most system, especially for nonlinear systems with constraints, the arrival cost cannot be solved analytically.
One option for this problem is to approximate the arrival cost, or to transfer problem (43) tô whereẐ T −N (z) is the approximation of Z T −N (z). If denote the solution of problem (47) as, Before giving the convergence theorem, define following, Definition 1: If there exist a positive integer N 0 and a K function ϕ(·), then for any two states x 1 and x 2 , and any k > 0, which satisfies then the system is uniformly observable.
Definition 2: Consider following system then for any initial condition x 0 and time k ≥ 0, x k (k; x 0 , 0) ∈ X k , and for any ε > 0, there exist δ > 0 and a positive integer and for all x 0 ∈ X 0 , when T → ∞, then this observer is an asymptotic stable observer for system (51). Theorem 1: If following conditions holds, (i) f k (·) and h k (·) are Lipschitz continuous; (ii) the stage cost function L k (·) and the initial penalty function (·) are lower semicontinuous function; (iii) (w, v) ∈ (W k × V k ), x ∈ X 0 ,x 0 ∈ X 0 , k ≥ 0, there exist K function η(·) and γ (·), (iv) there exist an initial condition x 0|∞ , and external disturbance sequence {w k|∞ } ∞ k=0 , for all k ≥ 0, (x 0|∞ , {w k|∞ }) ∈ k (v) system (36) is uniform observable, and there exist a K functionγ (·) satisfying following, then if there exist positive number δ w and δ v , which W k ⊆ In addition, there exist {Ẑ j (·)}, which satisfies (57), as shown at the bottom of the page, then for allx 0 ∈ X 0 , the moving horizon estimator is an asymptotic stable observer.
The proof of stability and convergence follows the same lines of reasoning of [34], and is omitted here for brevity.

IV. EXPERIMENTAL STUDY
In order to verify the feasibility and performance of the method, simulation experiments were performed using actual data collected from industrial sites. The steps of this experiment are as follows.
(i) Determining the experiment period: Since the test of sulfur concentrate is performed every eight hours in the plant, the experiment period is set to 8 hours. If the current testing time point is T N , then the sulfur concentrate grade is estimated from T N to the next test time point. When the next test time point T N + 1 is reached, T N + 1 is set as the initial time and a new round of moving horizon estimation is started; (ii) Data sample construction: Firstly obtain the sulfur concentrate grade value at the initial and final time of each measurement period from the test data. Use the data collected in real time on the industrial site to calculate g(θ, β), i.e., the output y k of the system (34), then identify the current working condition, and determine the range of the state x k according to the working condition category. At the same time, determine the range of external disturbance and output detection error according to the actual working condition and experience. The interval between data calculation is 10 minutes; (iii) Moving horizon estimation: Determine the value of N (N = 10 in this experiment), construct a moving horizon estimator according to the conditions in Theorem 1 and estimate the sulfur concentrate grade.  In this experiment, a total of 50 sets of continuous testing data samples were selected. The 50 sets of continuous data included four types of working conditions, C2 to C5. The experimental results are shown in Table 5 and Figure 20. In order to illustrate the advantages of the moving horizon estimation (MHE) method, it was compared with the kinetic model under different working conditions.
From the experimental results, it can be observed that the proposed method can measure the sulfur concentrate grade more accurately. Compared with the method only based on kinetic model, by using this method, there is no large fluctuations of the estimation value during the experimental period. The estimation curve is smoother. This indicates that the proposed approach can better adapt to the change of working conditions. This is because: (1) The structure of the kinetic model is determined according to the reaction laws. The model parameters are identified for each working condition, respectively. Therefore, the kinetic model can cover the dynamics of different working conditions. (2) The proposed method uses the kinetic models as the nominal models, and also uses online data to estimate the process state and disturbance. The proposed method explicitly includes external disturbance and measurement error. Therefore, this approach can extract useful information from the online measurements to enhance the estimation performance.
However, the accuracy and stability of the proposed approach is limited by the quality of data samples and the performance of working condition classification.

V. CONCLUSION
In this article, a moving horizon approach to predict sulfur concentrate grade is proposed. The approach utilizes both flotation kinetics, froth image features and real-time process measurements. By model reformulation, outlier elimination and stepwise identification, the robust identification approach decreased the influence of uncertainties of identification samples. By identifying the model parameters under different working conditions, the kinetic models can handle the multi-modal characteristic of flotation process. By embedding the kinetic model into the moving horizon estimation framework, the sulfur concentrate grade can be reconstructed online using real-time process measurements. The experimental study illustrates the feasibility and performance of the proposed approach. The improvement of estimation performance is due to the increase of information utilization rate. Both kinetics and froth image features, historical and real-time data are utilized [35]. As more information is used in the estimation framework, the uncertainty in the estimation can be decreased. To further increase the estimation performance, future research directions include: systematically study the features related to the dynamics of sulfur concentrate grade; embed more information related to the dynamics of sulfur concentrate grade in the estimation framework, e.g., trend information [36]; increase the optimization performance of the estimation algorithm [37]; combine deep learning, reinforcement learning and control theory with the estimation framework [38]- [41]; design robust moving horizon estimation to increase the stability of estimation performance.