Closed-Loop Dynamic Blending Optimization Based on Variational Bayesian and its Application in Industry

There exist uncertainties in raw material components and distribution parameters in blending process. Considering the dynamic uncertainty of distribution parameters and the dynamic change of raw material inventory, a general closed-loop dynamic blending optimization (CDBO) is proposed. The proposed method consists of both data-driven system and model-based system. In the data-driven system, the feedback optimization problem is reconstructed into a linear regression problem, and the variational Bayesian is proposed to obtain the mean and variance of the distribution. Then, according to the incoming raw materials, the mean value of each distribution parameter is adjusted by expert rules. Last, in the model-based system, blending ratio is obtained by chance-constrained programming. To verify the effectiveness of CDBO, a detailed derivation process of variable Bayesian, expert rule and chance-constrained programming are established for zinc smelting process. Numerical studies and an industrial application for zinc smelting process are presented to demonstrate the advantages of the proposed method. Compared with the manual blending, the volatility index of chemical tests data is greatly reduced and the compliance rates of zinc and lead is increased by 6.7% and 3.3%, respectively.

The main problem of blending optimization is uncertainty, which includes two meanings, one is the uncertainty in the component of raw materials, and the other is uncertainty of distribution parameters. Due to the wide range of feed sources and the variety of raw materials, there are The associate editor coordinating the review of this manuscript and approving it for publication was Emanuele Crisostomi .
varying degrees of uncertainty in the component of raw materials in most processing industries [1], [2]. In the practical production, the uncertainty of raw material components has a great impact on the price and availability of production [5]. And chance-constrained programming is widely used to solve the problem of blending optimization under uncertainty in industry, including ironmaking [2], zinc smelting [1], oil refining [4], [8], [10] and so on. However, these studies are based on the open-loop method without considering the uncertainty of distribution parameters. The distribution parameters during practical production are not necessarily accurate, and the inventory of raw materials is dynamic. Therefore, the distribution parameters of each component should be adjusted dynamically according to the chemical tests data and the incoming raw materials.
The problem of closed-loop blending optimization is of great importance in academia and industry. For the deterministic optimization model, Chèbre et al. proposed a control algorithm for on-line blending systems, while the uncertainties are treated in closed-loop by an estimator of the components properties [25]. On the basis of Chèbre's study, Chen and Yang proposed a double loop optimization strategy to reduce the re-blending times of final product [10]. However, many industrial blending systems do not have on-line data, and many of them are off-line tests with large lag. Therefore, a closed-loop dynamic optimization method is proposed to solve the uncertainty of distributed parameters in chance constrained programming. Since the number of measurements is smaller than the number of optimization parameters in industrial process, the distributed parameter cannot be obtained directly from the measurements for blending optimization. However, the problem can be formulated in a Bayesian framework, which provides certain distinct advantages, including less interactions, and full utilization of prior knowledge than other optimization methods. The most important thing is that Bayesian framework can provide a posterior distribution rather than a point estimate and, therefore, providing an estimate of the uncertainty, which, can be used as a feedback mechanism for adapting the data acquisition process [16]. Furthermore, all required model parameters can be estimated, including mean and variance, not just mean [17]. There are two kinds of reasoning algorithms in Bayesian method [15], [16], [17] dealing with complex posterior distribution. One is Markov chain Monte Carlo (MCMC) sampling which is time-consuming when the number of unknown variables is large [15]. The other is based on approximation techniques [15], [16], [17], [18], [19], such as variational Bayesian or expectation propagation which work by directly utilizing another easily estimated distribution to approximate the desired posterior distribution.
Aiming at the uncertainty of raw material composition, the dynamic change of distribution parameters, and the small amount of chemical test data, closed-loop dynamic blending optimization (CDBO) method is proposed by putting the variational Bayesian into the feedback mechanism. Firstly, the chemical tests data and the blending ratio are fed back, and the variational Bayesian is used to adjust the distribution parameters of each component in each mine. Then, expert rules are established to adjust the average value of each component in each bin according to new raw material. Finally, the blending ratio is obtained by solving the chance constrained programming with the distribution parameters. To demonstrate the efficacy of the CDBO, an application and detailed variational Bayesian in Zinc smelting process are presented. The main contributions of this study are as follows: 1) A general closed-loop dynamic blending optimization under uncertainty is proposed by using the variational Bayesian as the feedback mechanism to solve the uncertainty of raw material composition and distribution parameters.
2) The feedback optimization and detailed variational Bayesian of zinc smelting are established to solve the dynamic change of distribution parameters in zinc smelting process.
3) The verification idea of zinc smelting closed-loop dynamic blending optimization is proposed, which includes offline data-driven system verification and overall verification of blending field trial operation.
4) This method has been successfully applied to solve the problem of blending optimization in zinc smelting, which reduces the volatility index of chemical tests data and improves the compliance rates of zinc and lead.
The rest of this paper is organized as follows. The general closed-loop dynamic blending optimization (CDBO) is introduced simply in Section 2. The application and detailed inference process of CDBO in zinc smelting are presented in Section 3. The application of CDBO in zinc smelting is verified in Section 4. Section 5 concludes the study.

II. GENERAL CLOSED-LOOP DYNAMIC BLENDING OPTIMIZATION
In this section, we first provide a brief introduction to general closed-loop dynamic blending optimization as shown in Figure 1, which includes two parts: data-driven system (including variational Bayesian and expert rules) and modelbased system.  [20]. And the feedback optimization model of distribution parameters is established by feedback based on the blending ratio and chemical tests data. We employ variational Bayesian to approximate the posterior by a simpler distribution and then finding the member of that family which minimizes the Kullback-Leibler (KL) divergence to the true posterior [15]. For conjugate variables, it is easy to obtain closed-form updates [18]. For nonconjugate variables, Laplace variational inference [18] which uses a Talor approximation around the maximum a posterior (MAP) point to construct a Gaussian proxy for the posterior is proposed.

2) EXPERT RULES
Based on the feedback optimization mechanism, the distributed parameters are updated. However, the new raw materials will change the composition parameters of each bin. VOLUME 11, 2023 Field experts have rich experience in this aspect, so expert rules are used for secondary adjustment. The key to the success of an expert system lies in the knowledge level of system domain experts and their understanding of domain knowledge. For the new incoming raw materials, experts in the field of blending adjust the mean value of distribution parameters of each bin according to the quantity and quality of incoming materials and inventory. Therefore, the expert rules of different blending optimization can be established according to different fields. The general expression of the rule is: IF Condition THEN Conclusion. Generally speaking, conditions and conclusions can be composed of multiple statements connected by AND and OR.

B. CHANCE-CONSTRAINED PROGRAMMING
Through the data-driven system, the distribution parameters of each component can be obtained. And chance-constrained programming (CCP) is a powerful paradigm for optimization under uncertainty, whose goal is to meet the constraint conditions with a certain probability. The optimal proportioning can be obtained by putting the distribution parameters into the CCP. The general form of the blending problem is as follows min h(X 1, X 2, X 3,w,λ,ξ ) where h is the objective function; g = (g 1 ,g 2 , · · · , g m ) represents a constraint mapping and parameter 1−α is the set risk level; X 1,X 2 and X 3 represent decision variables.w,λ andξ represent random variables with different distributions.

III. APPLICATION OF CDBO IN ZINC SMELTING
Based on the CDBO, a closed-loop dynamic blending optimization model in zinc smelting process is established.
In the data-driven system, the variable Bayesian model of zinc smelting process with its detailed derivation and expert rules are established. In the model-based system, the chance constrained programming model and hybrid intelligent optimization method for zinc smelting process are given. According to reference [1], the distribution of each component in each bin of zinc smelting is shown in Table 1. The initial values of each distribution are obtained by statistical method. The blending process is to mix different raw materials in proportion. The feedback optimization model of distribution parameters and the matrix form are Eq. (2) and Eq. (3), respectively.
where the vector z is the chemical tests data, M is the number of chemical tests data at a time; X 2 is the corresponding ratio of normal distribution; The vector w is the composition of ore bin, which is normal distribution, and the dimension is N 1; X 2 is the corresponding ratio of lognormal distribution; λ is the composition of ore bin, which is lognormal distribution; And the dimension is N2; It can be seen from Table 1 that N 1 and N 2 of zinc, lead and silica are different;The noise is assumed to be independent Gaussian noise The conditional distribution of z is given by where I M is an identity matrix of order M . For each parameter w j and λ j , Gaussian prior and Lognormal prior are set as follows: where θ −1 j is the variance. We impose conjugate priors on the parameters θ j and τ , namely, θ j ∼ Gamma(θ j |c 0 , d 0 ), τ ∼ Gamma(τ |a 0 , b 0 ), where a 0 , b 0 , c 0 and d 0 are hyperparameters. The Bayesian graph and the probability density function of relevant random variables are shown in Figure 2.
Subsequently, the joint distribution of all variables is as follows: 496 VOLUME 11, 2023 Our goal is to obtain the posterior distribution of w, λ, θ and τ . According to the Bayesian theorem, we have However, the margin distribution (9) is computationally infeasible, because integral is very difficult to calculate. The basic idea of variational Bayesian is to employ variational distribution q(w, λ, τ, θ) to approximate the posterior one. Then Kullback-Leibler (KL) divergence is used to measure the difference between the two distributions [15]. Based on variational Bayesian, the original problem is transformed into the following optimization problem KL (q(w, λ, τ, θ) p(w, λ, τ, θ|z)) (10) According to [1], variational distributions for each item are independent. The variational distribution can be expressed as Then, the optimal variational posterior distribution can be obtained by the Theorem 1. And the main steps of the solution process are listed in Table 2.
Theorem 1. The optimal variational posterior distributions of our model are as follows.
For conjugate part: where where the expression of E[λλ T ] is Eq. (S1) of Appendix A. For nonconjugate part: where where λ j can be obtained by MAP.
Proof. Please see Appendix A for details.

2) EXPERT RULES FOR ZINC SMELTING
In zinc smelting, the distribution parameters of each ore bin are changed due to the incoming zinc concentrate. The experts have rich field experience on how to deal with the mean value of each component distribution. Therefore, the industry experience of experts in this paper is input into the system rule base. If a new zinc concentrate is put into storage, the average value of each component is updated. The input is inventory with the average value of each component, the weight and composition of new warehousing raw materials, and the output is the updated mean value of each component. The expression is as follows.
where y1 is the average content of each component in each inventory. y2 is the declared value of each component provided by the supplier. Different rules are made for different ore bin and different composition. Take zinc in highquality ore as an example, the rule table is shown in Table 3. The inventory is M and the weight of incoming material is L. Pr where J is the introduced intermediate variable; Constraint (1) denotes minimize the pessimistic value of the Original objective function, which is to minimize the zinc content per ton; Constraint (2) indicated that the minimum percentage of zinc content in the mixed zinc concentrates should be larger than the technical requirement. Constraint (3) indicated that the maximum percentage of each impurity in the blend should be met. Constraint (4) limited the quantity of the mixed zinc concentrates and constraint (5) limited the quantity of each zinc concentrate. And a hybrid intelligent optimization method (HIOM) [1] composed of Monte Carlo, neural network and genetic algorithm is used to solve MiniMin CCP. The Flowchart of the hybrid intelligent optimization method is shown in Figure 3. Finally, the optimized ratio can be obtained for blending.

IV. SIMULATION AND APPLICATION VERIFICATION IN ZINC SMELTING
At present, according to the experience of workers, trial and error method is used to obtain the proportion of ingredients in zinc smelters. And there are three main problems. First, they could only obtain a feasible solution, not an optimized one; Second, they don't directly consider uncertainty of each component of zinc concentrates; Last, they don't take full use of the feedback from the chemical tests data. Our method solves these problems to some extent.  To verify the effectiveness of the closed-loop dynamic blending optimization model in zinc smelting process, the verification is carried out in two steps: Firstly, the effectiveness of the data-driven system is verified. Based on the historical data of manual blending, the calculated value Z c of blending results is obtained by using the optimized average value u of each component in each bin and the actual manual blending ratio. Secondly, the validity of the datadriven system was verified by comparing with the test data; Then, the overall verification of the system is carried out. The zinc smelting blending optimization system had been put into operation for a period of time, and the results of manual blending and actual optimization are statistically analyzed and compared.

A. VERIFICATION OF DATA-DRIVEN SYSTEM
The historical data are collected of manual blending to verify the data-driven system. The manual blending of a smelter in China is divided into two shifts in the morning and the evening, and the test results are delayed by one day. The block diagram is shown in Figure 4. First, the artificial ingredients and the corresponding chemical tests data are fed back. Second, the mean and variance of distribution parameters of each component in each bin are adjusted by using the variational Bayesian method. Then, according to the new raw    According to the data of August 2020, we set the parameters of expert rules. The chemical tests data and corresponding calculated values Z c of zinc, lead and silica in one month are shown in Figure 5. The accuracy within each error range is shown in Table 5. According to the statistics, the accuracy between the chemical tests data and the calculated value of zinc is 85.5% within the range of plus or minus 0.4, the accuracy of silica is 82.3% within the range of plus or minus 0.25 and the accuracy of lead is 83.9% within the range of plus or minus 0.2. Thus, it can be seen the effect of datadriven system is good.
To verify the adaptability of the data-driven system, the data are tested in September. The chemical tests data and corresponding calculated values Z c are shown in Figure 6. The accuracy within each error range is shown in Table 6. The accuracy between the chemical tests data and the calculated value of zinc, silicon and lead are 80%, 86.7% and 81.7%, respectively. Therefore, the data-driven system has strong adaptability.
However, there are some points with large errors, as shown in the green and turquoise dotted box in Figure 5 and Figure 6. The main reasons are as follows: 1) For the green dotted box, the chemical analysis method is delayed for about one week, so the composition of the ore bin is adjusted according to the declared value of supplier. However, there are large errors between the declared values (DV) and the chemical tests data (TD) of raw materials, as shown in Table 7. 2) For VOLUME 11, 2023  turquoise dotted box, there are two situations: One situation is when new raw materials are put into the warehouse during the blending operation the proportion cannot be adjusted; Another situation is that the new raw materials have been put into the warehouse, but the batcher have not placed the data into the system in time. All of these will lead to the components of ore bin not adjusted in time, leading to the jump of test results and calculation results.

B. VERIFICATION OF CLOSED-LOOP DYNAMIC BLENDING OPTIMIZATION SYSTEM
To verify the closed-loop dynamic blending optimization, we have compiled the zinc smelting blending optimization software for field trial operation verification. The blending optimization software is shown in Figure 7. The blending software framework is shown in Figure 7a. First, data can be obtained from the enterprise MES system and stored in the local SQL server. Then, we can distribute zinc concentrate to each warehouse (Figure 7c), get the optimized ratio ( Figure 7b) and statistical analysis of historical results (Figure 7d) through the human-computer interface. Last, the variational Bayesian model, blending optimization model and hybrid intelligent optimization algorithm are run on the big data platform. This project is led by a university with cooperation from an IT company. We are mainly responsible for human-computer interface, SQL database and optimization algorithm development based on big data platform. The human-computer interface is developed in C# language. The database is managed by SQL Server 2012. Big data platform supports different development languages such as Python, R, Pyspark and so on to implement custom According to the industrial demand, it is hoped that the more times each component can meet the requirements and the fluctuation of blending results will be small. Due to the uncertainty of each component, we define two statistical indicators to reflect the quality of the blending results, which are compliance rates and Volatility index. The compliance rate (CR) represents the proportion of zinc, lead and silica reaching the standard each month in the total number of tests. At the same time, we define the volatility index (VI)  to describe the volatility of the results. It is defined as where i is zinc, lead and silica; n is the total amount of data; data i j is the jth data of i element; data i is the average value of i element.
To compare proposed method (Pro) with the actual artificial blending (Art), a total of 60 blending operations were carried out, including using artificial blending in the first half of the month, and blending optimization system in the second half of the month using the blending optimization system. Through field operation, the chemical tests data of zinc, lead and silica are shown in Fig. 8. The corresponding compliance rates (CR) and volatility index (VI) are shown in Table 8. It can be seen that the fluctuation of zinc, lead and silica is small, and the compliance rate of zinc and lead is increased by 6.7% and 3.3% respectively. This can improve product quality and stabilize the production process. At the same time, it can reduce the number of maintenance, indirectly reduce the production cost and improve the production efficiency. In summary, the experiments reveal the superiority of the method recommended in this paper over artificial ingredients in terms of all the metrics.

V. CONCLUSION
In the blending process, especially for the process industry, There exist uncertainties in raw material components and distribution parameters. And the chemical tests data have the characteristics of large lag and small sample size. In addition, the inventory is dynamic. Thus, a closed-loop dynamic blending optimization was proposed based on variational Bayesian. The recommended method consists of two parts: data-driven system and model-based system. Hence, the distribution parameters of each component and the blending ratio could be obtained by data-driven system and modelbased system, respectively. To verify the effectiveness of the CDBO, a detailed variable Bayesian model derivation, expert rule and chance constrained programming model were established for zinc smelting process. After the numerical analysis of field historical data, it could be included that the data-drive system could meet the needs of the field. Compared with the manual blending, the volatility index of chemical tests data is greatly reduced and the compliance rates of zinc and lead is increased by 6.7% and 3.3%, respectively.

APPENDICES A PROOF OF THEOREM 1
In this part, a detailed variational derivation is provided for Theorem 1.
At first, the variational Bayesian method is introduced to derive the optimal solution of the general model. The optimal variational distribution is expressed as where x is all the variables to be inferred and x j is the jth variable; z is the observed data. Based on the KL divergence, we have The evidence lower bound (ELBO) is defined as where E q(x) denotes the expectation operator on the (w.r.t.) variational distribution q(x). So, Eq. (21) can be equivalent to The log p(z) is a constant that does't affect the q(x), so minimizing KL divergence is equivalent to maximizing ELBO. Since the variational distribution is independent, ELBO can be expressed as VOLUME 11, 2023 where E −q(x j ) denotes the expectation operator w.r.t. all variational distributions except q(x j ). The const. represents all items that do not depend on x j . Therefore, the variational distribution for x j should satisfy Therefore, the variational distribution q(z) can be attained by coordinately updating q(x j ). For conjugate variables, it is easy to obtain closed-form updates [18].
where r = z − X 2 * λ. Obviously, q(w) is also a normal distribution. We denote it by N (w|v, ψ) with Obviously, q(τ ) is still a Gamma distribution. Where Obviously, θ j is still a Gamma distribution. Where Laplace variational inference [18] which uses a Talor approximation around the maximum a posterior (MAP) point to construct a Gaussian proxy for the posterior is proposed. For example, for a nonconjugate variable λ, define the function g(λ) A Taylor expansion around λ gives where λ is the MAP of g(λ); H ( λ) is the Hessian of g(λ) evaluated at λ, Since λ is the maximum of the g(λ), its gradient ∇g(λ)| λ= λ is zero. So, Eq. (35) can be approximated by Therefore, g(λ) can be approximated by shown at the bottom of the next page, where Log and ln are base-e logarithm and y = z − X 1 * w. Define the function g(λ j ) to contain the terms inside the update And we define the approximate distribution of q(λ j ) as N (u j , j ). λ j can be obtained by MAP. Based on Eq. (37) and Eq. (39), so u = λ j and Thus, the approximate update for q(θ) is to set it to The mean E[·] and variance V [·] of lognormal distribution can be expressed as From Eq. (41) to (43), it can be concluded that: .5f i +0.5f j ), therefore, the expression of E[λλ T ] is Eq. (S1), shown at the bottom of the page. (1, ψ), where 1 represents the degree of freedom. Then, there is  exp(e 1 +e 2 +0.5f 1 +0.5f 2 ) · · · exp(e 1 +e N 2 +0.5f 1 +0.5f N 2 ) exp(e 2 +e 1 +0.5f 2 +0.5f 1 ) exp(2e 2 +2f 2 ) · · · exp(e 2 +e N 2 +0.5f 2 +0.5f N 2 ) . . . . . . . . . . . .

5) According to 3) and 4), we have
exp(e 1 +e N 2 +0.5f 1 +0.5f N 2 ) exp(e 2 +e N 2 +0.5f 2 +0.5f N 2 ) · · · exp(2e N 2 +2f N 2 )      (Sl) VOLUME 11, 2023 6) According to 4), we have E w T X 1 T X 1 * w = E trace w T X YUDONG LI received the B.S. and M.S. degrees in automation and electrical engineering from Central South University. He is currently pursuing the Ph.D. degree in control science and engineering with Northeastern University. His research interests include machine learning and industrial big data. VOLUME 11, 2023