A Dynamic Model Incorporated with Delay Estimation and Variable Selection for Control of SCR De-NO_x Process

Abstract

In order to control NH₃ injection for the selective catalytic reduction of nitrogen oxide (NO_x) denitration (SCR de-NO_x) process, a model that can accurately and quickly predict outlet NO_x emissions is required. This paper presents a dynamic kernel partial least squares (KPLS) model incorporated with delay estimation and variable selection for outlet NO_x emission and investigated control strategy for NH₃ injection. First, k-nearest neighbor mutual information (KNN_MI) was used for delay estimation, and the effect of historical data lengths on KNN_MI was taken into account. Bidirectional search based on the change rate of KNN_MI (KNN_MI_CR) was used for variable selection. Delay–time difference update algorithm and feedback correction strategy were proposed. Second, the NH₃ injection compensator (NIC) and the outlet NO_x emission model constituted a correction controller. Then, its output and the output of the existing controller are added up to suitable NH₃ injection. Finally, the KNN_MI_CR method was compared with different algorithms by benchmark dataset. The field data results showed that the KNN_MI_CR method could improve model accuracy for reconstructed samples. The final model can predict outlet NO_x emissions in different operating states accurately. The control result not only meets the NO_x emissions standard (50 mg/m³) but also keeps high de-NO_x efficiency (80%). NH₃ injection and NH₃ escape are reduced by 11% and 39%.

1. Introduction

How to reduce nitrogen oxide (NO_x) emissions from coal-fired power plants has become an important task of environmental treatment. Selective catalytic reduction denitration (SCR de-NO_x) technology is one of the mature flue gas denitration technologies [1]. Many efforts have been made to optimize the SCR system, such as superior catalyst technology, enhanced de-NO_x control strategy, and byproduct limitation [2]. However, it is very difficult to control the amount of NH₃ injection accurately. First, new energy is integrated into the power grid on a large scale in recent years. The boiler load needs to be changed to stabilize the fluctuation from the output power of new energy frequently. Therefore, the wide range of operating states poses certain difficulties in NH₃ injection control. Second, the upper limitation of NO_x emission is decreased to 50 mg/m³ in China. In order to meet the national environmental standard, excessive NH₃ injection occurs easily. Hence, it is required to control NO_x emissions within the upper limitation and to avoid NH₃ escape. Third, the SCR de-NO_x system has large inertia, and it also has a large delay in NO_x analysis. When the operating state changes, the SCR system cannot adjust NH₃ injection in time based on the sampling inlet NO_x concentration, which can make the outlet NO_x emission deviate from the set value. Furthermore, the SCR de-NO_x process, which has the feature of nonlinear and multivariable coupling, is complex. When boiler load is stable, outlet NO_x emission can fluctuate widely for the reasons of flue gas temperature, flue gas flow rate, NH₃ injection and oxygen content. If boiler load changes, outlet NO_x emission changes dramatically. The proportion integration differentiation (PID) cascade control cannot achieve optimal control, which ensures high de-NO_x efficiency and low NH₃ escape. Therefore, the establishment of an accurate SCR model and the research of optimal control strategies are significant.

The current control strategy mainly uses a transfer function to control NH₃ injection. Considering the strong coupling between boiler parameter and SCR process parameter, the transfer function is difficult to reflect SCR de-NO_x process accurately. With the application of a plant-level supervisor information system, operating data can be viewed and recalled easily. Therefore, data-driven modeling and control methods are widely used. Data-driven modeling methods, such as Hammerstein–Wiener [3], radial basis function autoregressive exogenous [4], online support vector regression [5], neural network, nonlinear autoregressive exogenous and polynomial fitting [6], have been used to predict NO_x emission. Although most models are applied for diesel engines, the reaction principle is the same as the SCR de-NO_x system of a coal-fired power plant. The differences are that the boiler parameter affects outlet NO_x emissions, and the NO_x analyzer has a large lag in the coal-fired power plant. In addition, the advanced thermal management control strategies of catalysts have been developed to reduce emission concentrations [7].

Variable selection and delay estimation are investigated to simple model structure and improve model accuracy. Mutual information (MI), which denotes statistical dependence between any two variables, is a measure of variable correlation [8]. Therefore, it can express a linear and nonlinear correlation [9]. In order to improve the model accuracy, Stanišić et al. [10] used the MI method to select input variables. However, how to select an evaluation strategy and how to estimate the MI value accurately are two main difficulties in variable selection by the MI method. In their evaluation strategy, Battiti et al. proposed the mutual information feature selection (MIFS) method [11], which got a good result. After this, Peng et al. [12] and Fleuret [13] proposed a different evaluation strategy. The above methods mainly use one-dimensional MI to evaluate the variable set. When there are redundant variables, the subset of the optimal variables cannot be obtained. The literature [10] used a single evaluation function as an indicator, which often leads to an imbalance between the relevant variable and redundancy variables. Considering the large time-delay process, the effect of time-delay needs to be considered in variable selection. Ludwig et al. [14] used the MI method to estimate the time-delay and got better results in a case. However, the method ignores the effect of variable dimension on the MI value. In general, output variable y(t) not only related to input variable x_i(t) at time t, but also related to historical input x_i(t − 1), x_i(t − 2), …, x_i(t − l_i) before time t. Here, l_i is defined as the historical data length of the input variable x_i. Therefore, it is necessary to determine the historical data length of each input variable when the time-delay is estimated.

Recently, adaptive nonlinear model predictive controllers, such as model-based predictive controller [15], the generalized predictive controller [16], dynamic matrix controller [17], nonlinear autoregressive moving average with exogenous inputs (NARMAX) based predictive controller [18], the neural network-based predictive controller [19] and subspace predictive controller [20], have been used to control the NH₃ injection and got better effect compared with the conventional PID controller. However, some models (e.g., [15,16]) involved high order nonlinear optimization, which is time-consuming and gets a local minimum rather than a global minimum easily. Not all of the operating states are covered in the NARMAX prediction model. The neural network prediction model is difficult to be predicted online, quickly and accurately. The input variable of the subspace predictive model only relates to NH₃ injection. Hence, it cannot reflect other variables that affect outlet NO_x emission correctly. In addition, the control effect is also affected by the constraints, such as the NH₃ valve opening and the increment of NH₃ valve opening, which are neglected in the above literature.

The objective of this paper is to model outlet NO_x emissions under the different operating state, thereby to facilitate controller design. The kernel partial least squares (KPLS) model shows a good effect when it deals with the high dimensionality and colinearity data. In this paper, a dynamic KPLS model incorporated with variable selection and delay estimation is proposed. The k-nearest neighbor mutual information (KNN_MI) was used for delay estimation, and the effect of historical data lengths on KNN_MI was taken into account. Based on the results of delay estimation, a bidirectional search based on the change rate of KNN_MI (KNN_MI_CR) was used for variable selection. In order to improve model update accuracy, delay–time difference (DTD) update algorithm and feedback correction strategy are proposed. In order to take advantage of the existing controller, which has rapid response-ability to deal with disturbance signals, a parallel control structure that can identify the dynamic model online by using operation data is proposed. The existing conventional controller is used as the main controller. The NH₃ injection compensator (NIC) combined with the outlet NO_x emission model are used as correction controller, which can ensure the outlet NO_x emission is stable near the set value. The principle of NIC is to eliminate the adverse effects of the large inertia of the SCR system and the lag of NO_x analysis through feed-forward compensation to a certain extent.

The paper is organized as follows: In Section 2, the SCR de-NO_x process and its control strategy are described. In Section 3, the delay estimation and variable selection method are presented. The experimental results of the benchmark dataset for variable selection are presented in Section 4. The field data experiment for the SCR de-NO_x process is investigated in Section 5, and the predictive accuracy of the outlet NO_x emission model and the control effect of the proposed control strategy are analyzed. Finally, concluding remarks are given in Section 6.

2. The SCR de-NO_x Process and Its Control Strategy

2.1. The SCR de-NO_x Process

For a 1000 MW ultra-supercritical coal-fired power plant, the SCR de-NO_x system is shown in Figure 1. The temperature of the SCR de-NO_x system is controlled by the bypass gas damper of the economizer in the upstream flue. Mixed with dilution air, the ammonia gas is sprayed by the NH₃ flow rate valve through the ammonia injection grid (AIG). In the reaction zone, the selective catalytic reaction occurs on the stainless steel plates, which are supported by the V₂O₅-WO₃-TiO₂ catalyst. Finally, the NO_x in the flue gas is converted into harmless N₂ and H₂O, and thereby the flue gas is purified. The main reaction of the SCR de-NO_x process is described as follows [1]:

4NH₃ + 4NO + O₂→4 N₂ + 6H₂O

(1)

4NH₃ + 2NO₂ + O₂→3 N₂ + 6H₂O

(2)

Actually, adverse side reaction will occur as follows [1]:

2SO₂ + O₂→2SO₃

(3)

NH₃ + SO₃ + H₂O→NH₄HSO₄

(4)

2NH₃ + SO₃+H₂O→(NH₄)₂SO₄

(5)

It can be seen that the excessive NH₃ injection leads to side reaction, and de-NO_x efficiency is decreased because of the regeneration of NO_x. Excess NH₃ reacts with SO₃ in the flue gas to form NH₄HSO₄ and (NH₄)₂SO₄, which not only reduce the catalyst activity but also cause plugging and corrosion in the heated surface of the air preheater. In addition, increased NH₃ escape could lead to secondary pollution. Secondary pollution refers to new pollutants, that is, unreacted ammonia gas.

2.2. The SCR de-NO_x Process Control Strategy

Under the variable state of the power plant, the SCR de-NO_x process shows nonlinear, strong disturbance and large time-delay in NO_x analysis, so it is hard to achieve a good control effect by the existing PID control. Actually, except for controlling the temperature of the reaction zone by the bypass flue gas, it is necessary to design a better automatic control system which can control the NH₃ injection strictly. In order to ensure the most suitable NH₃/NO_x ratio, which is related to de-NO_x efficiency and ensure the NH₃ escape in a reasonable range, the proposed control system includes the main controller and the correction controller, which are in parallel operation. The main controller is the existing PID controller. The correction controller, which is shown in Figure 1, is investigated in this paper. The parallel control system is shown in Figure 2. The dynamic KPLS model is used as an outlet NO_x emission predictive model, which can be identified offline by the field data.

The correction controller includes three parts mainly: outlet NO_x emission predictive model, feedback correction and NH₃ injection compensator (NIC).

(1)

Outlet NO_x emission predictive model

In general, a dynamic data-driven model is constructed by adding the model input, such as the historical input data $x(t-1),\cdots ,x(t-w+1)$ and output data $y(t-1),\cdots ,y(t-w+1)$. Therefore, $x(t-\tau )$$x(t-\tau -1)$ and $y(t-1)$ are selected as the new input variables in this paper. As a kind of model update method, the time difference (TD) algorithm can work out the variable drift and achieve better predictive accuracy. Moreover, the data model based on the TD algorithm does not need frequent reconstruction and parameter updating [21]. The time-delay must be considered in the SCR deNO_x process. Hence, the delay-time difference (DTD) algorithm is proposed. Combining with the DTD method and KPLS model, the dynamic KPLS model, which is named DTD–KPLS, is proposed.

In general, the original input variable sample and output variable sample are used in the regression algorithm. However, $\Delta x(t)$ is calculated between the input variable x(t) at the time of t and the input variable x(t − i) at the time of t − i through the TD method first. Similarly, $\Delta y(t)$ is calculated by the TD method:

$$\Delta x(t)=x(t)-x(t-i)$$

(6)

$$\Delta y(t)=y(t)-y(t-i)$$

(7)

Then, the relationship between $\Delta x(t)$ and $\Delta y(t)$ is modeled as:

$$\Delta y(t)=\mathrm{f}\left[\Delta x(t)\right]$$

(8)

If new data $x(t^{\prime })$ is obtained, the TD of input data can be computed as:

$$\Delta x(t^{\prime })=x(t^{\prime })-x(t^{\prime }-i)$$

(9)

Thus, the TD of output data can be predicted by the training model as:

$$\Delta y(t^{\prime })=\mathrm{f}(\Delta x(t^{\prime }))$$

(10)

Finally, the predicted output value at time t′ is:

$$\widehat{y}(t^{\prime })=\Delta y(t^{\prime })+y(t^{\prime }-i)$$

(11)

When the DTD method is used, the time-delay $\tau$ needs to be considered in Equations (8) and (10). Hence, the training model becomes:

$$\Delta y(t)=\mathrm{f}\left[\Delta x(t-\tau )\right]=\mathrm{f}\left[x(t-\tau )-x(t-\tau -1)\right]$$

(12)

Hence, the actual predicted output is:

$$\widehat{y}(t^{\prime })=\mathrm{f}\left[x(t^{\prime }-\tau )-x(t^{\prime }-\tau -1)\right]+y(t^{\prime }-1)$$

(13)

In this paper, the estimation of the DTD–KPLS model from the training set is depicted as follows [22]:

• The relevant variable samples are obtained, and the preprocessing, which contains outlier elimination and data filtering, is performed;
• The input matrix X_real$\in {\mathit{R}}^{n\times m}$ and the output matrix $\mathit{Y}\in {\mathit{R}}^{n\times 1}$ are confirmed;
• The time-delay of each input variable x_i(I = 1,…,n) is estimated and the phase space reconstruction is made for X_real. Hence, the new input variable matrix X_st is obtained;
• The $\Delta \mathit{x}(t-\tau )$ and $\Delta \mathit{y}(t-1)$ of new input are calculated, and the $\Delta \mathit{y}(t)$ of output is calculated;
• The training set $\Delta {\mathit{X}}_{\mathrm{tr}}^0$ and $\Delta {\mathit{Y}}_{\mathrm{tr}}^0$ are confirmed and $\Delta {\mathit{X}}_{\mathrm{tr}}^1$ and $\Delta {\mathit{Y}}_{\mathrm{tr}}^1$ (z-score normalized) are obtained;
• The kernel matrix ${\mathit{K}}_{\mathrm{tr}}^0$ of the training set is calculated as:

  $${\mathit{K}}_{\mathrm{tr}}^0=k(\Delta x_{\mathrm{tr}},\Delta x_{\mathrm{tr}})$$

  (14)
• The kernel matrix ${\mathit{K}}_{\mathrm{tr}}^0$ of the training set is centered as:

  $${\mathit{K}}_{\mathrm{tr}}^1=\left(\mathit{I}-\frac{1}{n}{\mathbf{1}}_n{\mathbf{1}}_n^{\mathrm{T}}\right){\mathit{K}}_{\mathrm{tr}}^0\left(\mathit{I}-\frac{\mathbf{1}}{n}{\mathbf{1}}_n{\mathbf{1}}_n^{\mathrm{T}}\right)$$

  (15)
  Here, I is a unit matrix, 1_n is a matrix whose all elements are 1 and whose dimension is n.
• L is the number of latent variables, i iterates from 1 to L, and initialize the score vector ${\mathit{u}}^i$ of ${\mathit{X}}_{\mathrm{tr}}^1$ randomly.
• The score vector ${\mathit{t}}^i$ is calculated as:

  $${\mathit{t}}^i=\frac{{\mathit{K}}_{\mathrm{tr}}^1{\mathit{u}}^i}{\Vert {\mathit{K}}_{\mathrm{tr}}^1{\mathit{u}}^i\Vert }$$

  (16)
• The weight vector ${\mathit{c}}^i$ is calculated as:

  $${\mathit{c}}^i={\left({\mathit{Y}}_{\mathrm{tr}}^1\right)}^{\mathrm{T}}{\mathit{t}}^i$$

  (17)
• The score vector ${\mathit{u}}^i$ is calculated as:

  $${\mathit{u}}^i=\frac{{\mathit{Y}}_{\mathrm{tr}}^i{\mathit{c}}^i}{\Vert {\mathit{Y}}_{\mathrm{tr}}^i{\mathit{c}}^i\Vert }$$

  (18)
• Step (8) to (11) is repeated until that ${\mathit{t}}^i$ is converged;
• The matrix ${\mathit{K}}_{\mathrm{tr}}^1$ and ${\mathit{Y}}_{\mathrm{tr}}^1$ are reduced as follows: until that t and u are extracted;

  $${\mathit{K}}_{\mathrm{tr}}^{i+1}=\left[\mathit{I}-{\mathit{t}}^i{\left({\mathit{t}}^i\right)}^{\mathrm{T}}\right]{\mathit{K}}_{\mathrm{tr}}^i\left[\mathit{I}-{\mathit{t}}^i{\left({\mathit{t}}^i\right)}^{\mathrm{T}}\right]$$

  (19)

  $${\mathit{Y}}_{\mathrm{tr}}^{i+1}={\mathit{Y}}_{\mathrm{tr}}^i-{\mathit{t}}^i{\left({\mathit{t}}^i\right)}^{\mathrm{T}}{\mathit{Y}}_{\mathrm{tr}}^i$$

  (20)
• The regression coefficient B is obtained, and the regression model of the training set is obtained as follows:

  $${\widehat{\mathit{Y}}}_{\mathrm{tr}}={\mathit{K}}_{\mathrm{tr}}\mathit{B}={\mathit{K}}_{\mathrm{tr}}\mathit{U}{\left({\mathit{T}}^{\mathrm{T}}{\mathit{K}}_{\mathrm{tr}}\mathit{U}\right)}^{-1}{\mathit{T}}^{\mathrm{T}}{\mathit{Y}}_{\mathrm{tr}}$$

  (21)

Here, T and U are both matrices, which are formed by score vector t and score vector u.

The prediction for the test set is similar to that for the training set, except the computation of the test kernel matrix ${\mathit{K}}_{\mathrm{te}}^0$ and the centralization of ${\mathit{K}}_{\mathrm{te}}^0$:

$${\mathit{K}}_{\mathrm{te}}^0=k(\Delta x_{\mathrm{te}},\Delta x_{\mathrm{tr}})$$

(22)

$${\mathit{K}}_{\mathrm{te}}^1=\left({\mathit{K}}_{\mathrm{tr}}^0-\frac{1}{n}{\mathbf{1}}_{nt}{\mathbf{1}}_n^{\mathrm{T}}{\mathit{K}}_{\mathrm{te}}^0\right)\left(\mathit{I}-\frac{1}{n}{\mathbf{1}}_n{\mathbf{1}}_n^{\mathrm{T}}\right)$$

(23)

Here, nt is the number of the test set. In addition, the kernel function uses the radial basis function (RBF) kernel in this paper.

(2)

Feedback correction

In order to compensate for the prediction error, the feedback correction equation is as follows:

$$\tilde{y}(t+1)=\widehat{y}(t+1)+\Delta y(t+1)$$

(24)

$$\Delta y(t+1)=\rho \times \Delta y_0(t+1)+(1-\rho )\Delta y(t)$$

(25)

$$\Delta y_0(t+1)=y(t)-\widehat{y}(t)$$

(26)

Here, $\rho$ is 0.3, $\tilde{y}(t)$ is the corrected predicted value at time t, $\widehat{y}(t)$ is the predicted value at time t, $y(t)$ is the real value at time t.

Because the SCR de-NO_x process is very complicated, a complete mathematical model of the SCR process has still not been summarized. However, the compensation amount of NH₃ injection $\Delta y_{\mathrm{N}{\mathrm{H}}_3}$ can be calculated according to the nonlinear Equation (27):

$$\Delta y_{\mathrm{N}{\mathrm{H}}_3}(t)=\mathrm{f}(\Delta \tilde{y}(t))$$

(27)

Here,

$$\Delta y_{\mathrm{N}{\mathrm{H}}_3}(t)=y_{\mathrm{N}{\mathrm{H}}_3}^{\prime }-y_{\mathrm{N}{\mathrm{H}}_3}(t)$$

(28)

$$\Delta \tilde{y}(t)=\tilde{y}(t)-y_{set}$$

(29)

$y_{\mathrm{N}{\mathrm{H}}_3}^{\prime }$ is the demand amount of NH₃ injection at the current operating state. $y_{\mathrm{N}{\mathrm{H}}_3}(t)$ is the NH₃ injection flow rate at the time of t.

In this paper, the TD-KPLS model is used to establish the nonlinear model of the NIC, which is shown in Equation (27). The difference between the TD-KPLS model and DTD–KPLS model is that the samples of the TD-KPLS model are not reconstructed. In the model, the deviation $\Delta \tilde{y}(t)$ between the corrected DTD–KPLS model output $\tilde{y}(t)$ and the set value of outlet NO_x mission $y_{\mathrm{set}}$ is used as the input variables, and the compensation amount of NH₃ injection $\Delta y_{\mathrm{N}{\mathrm{H}}_3}(t)$ is used as the output variable. In addition, input variables also include inlet flue gas temperature and inlet NO_x concentration, which reflect the changing trend of the boiler load to a certain extent.

The adverse effects brought by the inertia characteristic of the SCR system are eliminated by feed-forward compensation. Hence, the NIC can compensate for the NH₃ injection in time to ensure that the outlet NO_x emission is stable at the set value.

In the SCR de-NO_x system, the NH₃ injection valve can be seen as a linear part, so the algorithm design and parameter setting are performed on the basis of that. However, the NH₃ injection valve has a certain dead zone and saturation zone, which shows nonlinearity. Therefore, in this paper, the middle curve in the on and off stroke is taken as the fitting curve between the NH₃ injection valve opening and the NH₃ injection flow rate.

In order to avoid valve saturation and action frequently, the following constraints of the NH₃ injection valve need to be considered:

$$U_{\mathrm{V}\min }\le \Delta u_{\mathrm{V}}(t+1|t)+u_{\mathrm{V}}(t)\le U_{\mathrm{V}\max }$$

(30)

$$u_{\mathrm{V}}{}_{\min }\le \Delta u_{\mathrm{V}}(t+1|t)\le u_{\mathrm{V}}{}_{\max }$$

(31)

Here, U_Vmax and U_Vmin are the upper and lower limits for the NH₃ valve opening. u_Vmax and u_Vmin are the upper and lower limits for the increment of the NH₃ valve opening. $\Delta u_{\mathrm{V}}(t+1|t)$ represents the increment of the NH₃ valve opening from the next time t+1 to the current time t.

3. The Delay Estimation and Variable Selection Method

3.1. Delay Estimation

Aiming at nonlinearity and large lag in SCR de-NO_x process, this paper proposes the k-nearest neighbor MI(KNN_MI) method to estimate the time-delay of the input variable and the effect of the historical data lengths on KNN_MI was taken into account. The KNN_MI method can be summarized as follows [8].

The KNN_MI of variable X and Y is:

$$\mathrm{MI}(X,Y)=\psi (k)-\langle \psi ({\mathit{n}}_x+1)+\psi ({\mathit{n}}_y+1)\rangle +\psi (n)$$

(32)

Here, n_x(i) is the sample number x_j whose distance from x_i is less than ε_i/2 strictly, ε_i/2 denotes the distance from x_i. to its k-th neighbor, and n_y(i) is calculated similarly, i ∈ [1, ⋯, n]. Ψ(x) is the digamma function, Ψ(x) = Γ(x)⁻¹dΓ(x)/dx. It meets the recursion equation Ψ(x + 1) = 1/x + Ψ(x) and Ψ(1) ≈ −0.5772156. The symbol $\langle \cdots \rangle$ denotes the mean of all variables in it.

The KNN_MI of high dimensional variables $(X_1,X_2,\cdots ,X_{\mathrm{m}},Y)$ is:

$$\mathrm{MI}(X_1,X_2,\cdots ,X_{\mathrm{m}},Y)=\psi (k)-\langle \psi ({\mathit{n}}_{x_1})+\cdots +\psi ({\mathit{n}}_{x_{\mathrm{m}}})+\psi ({\mathit{n}}_y)\rangle +m\psi (n)$$

(33)

Because of the time-delay between each input variable x_.i and the output variable y is different, the phase space of x_.i is reconstructed by embedding with time-delay ${\tau }_i\in \left[\min ({\tau }_i),\max ({\tau }_i)\right]$ ($\min ({\tau }_i)$ and $\max ({\tau }_i)$ are resolved by field experience). In addition, MI value has a relation to the length w of ${x^{\prime }}_{\cdot i}$. Therefore, the time-delay τ_i and length w are calculated at the time of t as follows:

$$\begin{array}{c}\underset{{\tau }_i={{\tau }^{\prime }}_i,w_i={w^{\prime }}_i}{\max }\mathrm{MI}\left({[x_{\cdot i}(t-{\tau }_i-w_i+1),\cdots ,x_{\cdot i}(t-{\tau }_i-1),x_{\cdot i}(t-{\tau }_i)]}^{\mathrm{T}},{[y(t-w_i+1),\cdots ,y(t-1),y(t)]}^{\mathrm{T}}\right)\\ \begin{array}{cc}\mathrm{s}.\mathrm{t}.& \min \left({\tau }_i\right)\le {\tau }_i\le \max \left({\tau }_i\right)\end{array},{\tau }_i+w_i<{\mathrm{T}}_{\max };i\in [1,m]\end{array}$$

(34)

Here, T_max is the maximum reaction time when the SCR system responds. Equation (34) involves constrained multivariable nonlinear optimization. If the number of input variables is m, 2 m parameters need to be obtained. Thus, within the above constraints, a global searcher maximizes the objective function by particle swarm optimization (PSO) algorithm. Hence, the optimal ${\tau }^{\prime }=[{{\tau }^{\prime }}_1,{{\tau }^{\prime }}_2,\cdots ,{{\tau }^{\prime }}_m]$ and $w^{\prime }=[{w^{\prime }}_1,{w^{\prime }}_2,\cdots ,{w^{\prime }}_m]$ are obtained. By estimating the time-delay ${{\tau }^{\prime }}_i$ of input variable x_.i, the reconstructed input matrix X_st is assumed as follows:

$$X_{\mathrm{st}}=\left[\begin{array}{ccc}{x}_{\cdot 1}(t-{{\tau }^{\prime }}_1-n+1)& \cdots & x_{\cdot i}(t-{{\tau }^{\prime }}_i-n+1)\\ ⋮& \ddots & ⋮\\ x_{\cdot 1}(t-{{\tau }^{\prime }}_1-1)& \cdots & x_{\cdot i}(t-{{\tau }^{\prime }}_i-1)\\ x_{\cdot 1}(t-{{\tau }^{\prime }}_1)& \cdots & x_{\cdot i}(t-{{\tau }^{\prime }}_i)\end{array}\right]$$

(35)

3.2. Variable Selection

The subset S, which has m input variables, equals to set $(X_1,X_2,\cdots ,X_m)$. The MI between S and Y is defined as:

$$\mathrm{MI}(S,Y)=\mathrm{MI}(X_1,X_2,\cdots ,X_m,Y)$$

(36)

Information entropy represents a measure of uncertainty. If the added variable d_i is irrelevant, the change of uncertainty is less. Therefore, the change rate of KNN_MI, which shown in Equation (37) is larger, when the added variable d_i is relevant:

$$R_{\Delta \mathrm{MI}}(S+d_i,Y)=\frac{\Delta \mathrm{MI}(S,d_i,Y)}{\mathrm{MI}(S,Y)}=\frac{\mathrm{MI}(X_1,X_2,\cdots ,X_m,d_i,Y)-\mathrm{MI}(X_1,X_2,\cdots ,X_m,Y)}{\mathrm{MI}(X_1,X_2,\cdots ,X_m,Y)}$$

(37)

The high-dimensional KNN_MI and the impact on original MI(S, Y) are considered in the above equation. Therefore,$R_{\Delta \mathrm{MI}}(S+d_i,Y)$ can be used to select variable subset S in forward-selection (step (2)–step (5)). In addition, the one-dimensional KNN_MI between any two variables in the variable subset S can be used to delete redundant variables in backward selection (step (6)–step (7)). Based on the change rate of KNN_MI (KNN_MI_CR), the algorithm for the variable selection is designed as follows:

(1)

One-dimensional MI(X_i,Y) between each input variable X_i and output variable Y are calculated, i = 1,…,m;

(2)

Based on one-dimensional KNN_MI values, the input variables are sorted from max to min. The most relevant variable which corresponds to the maximum KNN_MI value is selected as the initial variable set S;

(3)

Based on the most relevant variable and the sorting result of other input variables, $R_{\Delta \mathrm{MI}}(S+d_i,Y)$, which relates to the other input variable $d_i$ and the output variable Y, is calculated in sequence by Equation (37);

(4)

The calculated $R_{\Delta \mathrm{MI}}(S+d_i,Y)$ are sorted from min to max. The most relevant variable d_i, which corresponds to the minimum $R_{\Delta \mathrm{MI}}(S+d_i,Y)$ is selected. If $R_{\Delta \mathrm{MI}}(S+d_i,Y)<\alpha$ is satisfied, $S=S\cup d_i$. Here, $\alpha$ is a threshold value;

(5)

Step (3) and step (4) are repeated until all input variables were selected. The relevant variable subset S is finally obtained;

(6)

Backward selection: The one-dimensional KNN_MI between any two variables in S is calculated, and then the combination of redundant variables which has the max KNN_MI value is selected.

(7)

The redundant variables are eliminated by combining with the one-dimensional KNN_MI in step (2), and the optimal relevant variable subset S is finally obtained.

4. The Benchmark Dataset Experiment

In order to verify the proposed variable selection method, the Friedman dataset and Housing dataset [23] are used to experiment, and the predictive accuracy of the KPLS model is evaluated by the root mean square error (RMSE).

4.1. Friedman Dataset

Friedman function:

$$Y=10\sin (\pi X_1{X}_2)+20{(X_3-0.5)}^2+10X_4+5X_5+\epsilon$$

(38)

Here, the input variable $X_i(1\le i\le 10)$ is a uniform distribution between 0 and 1, ε is the Gaussian noise with a variance of 1. The variables from X₆ to X₁₀ are irrelevant variables. The variable X₁₁ = 0.5X₁ + ε is a redundant variable. The sample number is 500. The number of training sets is 338. The number of the test set is 162.

First, this paper uses the proposed variable selection method (KNN_MI_CR) to select variables. In Section 3.2, the one-dimensional MI(X_i,Y) calculated in step (1) is shown in Figure 3.

In Figure 3, if the variable selection is only based on one-dimensional KNN_MI, X₈ and X₁₀ are mistakenly selected, X₁ and X₂ are mistakenly eliminated.

In step (2), the most relevant variable X_4, which corresponds to the maximum one-dimensional KNN_MI value, is selected as the initial variable set S. Hence, the initial variable set S is X₄.

High dimensional KNN_MI and its change rate, which is calculated in step (4), are shown in Figure 4. There, +X₅ represents MI(X₄, X₅, Y), +X₃ represents MI(X₄, X₅, X₃, Y) and so on. The threshold value $\alpha$ is 0.1. Therefore, variables from X₆ to X₁₀ are eliminated, and the relevant variable subset S = (X₄, X₅, X₃, X₂, X₁₁, X₁).

In Figure 3, MI(X₁, Y), MI(X₂, Y) and MI(X₁₁, Y) are relatively close. In addition, MI(X₃, Y) and MI(X₅, Y) are also relatively close. Therefore, the relevant variable subset S could have a redundant variable. One-dimensional KNN_MI between any two variables in S is shown in Figure 5. It is obvious that MI(X₁, X₁₁) is the maximum value. In Figure 3, MI(X₁, Y) > MI(X₁₁, Y). Therefore, X₁₁ is a redundant variable. Finally, optimal relevant variable subset S = (X₄, X₅, X₃, X₂, X₁). This is consistent with Equation (38).

Second, the comparison among variable selection methods and the test results of the KPLS model are shown in

Table 1. The variable selection result and model test result for the Friedman dataset (no redundant variable).

Variable Selection Method	Variable Selection Result	RMSE
All variables	X1~X10	1.4966
MIFS [11]	X4, X5, X3, X1, X9	1.2019
mRMR [12]	X4, X5, X3, X1, X9	1.2019
CMIM [13]	X4, X5, X3, X7, X8	1.2433
BIF [24]	X4, X5, X3, X10, X9	1.2079
KNN_MI_CR	X4, X5, X3, X2, X1	1.1260

and

Table 2. The variable selection result and model test result for the Friedman dataset (contain redundant variable X₁₁).

Variable Selection Method	Variable Selection Result	RMSE
All variables	X1~X11	1.6222
MIFS [11]	X4, X5, X3, X1, X9	1.2043
mRMR [12]	X4, X5, X3, X1, X9	1.2043
CMIM [13]	X4, X5, X3, X11, X8	1.1843
BIF [24]	X4, X5, X3, X8, X10	1.3810
KNN_MI_CR	X4, X5, X3, X2, X1	1.1272

.

In Table 1 and Table 2, it can be seen that BIF, mRMR, MIFS, and CMIM have not selected all valid input variables, and the CMIM method selects redundant variable X₁₁ in the input variable subset. Compared with the other four variable selection methods, the KNN_MI_CR method selects all relevant input variables. Therefore, the KNN_MI_CR method proposed in this paper can select the input variable subset effectively.

4.2. Housing Dataset

The Housing dataset is used to predict Boston house price (defined as the output variable Y) by statistic values (defined as input variables X₁, ⋯, X₁₃), which represent 13 attributes in each house. The number of training sets is 338. The number of the test set is 168. The comparison of different variable selection methods and test results of the KPLS model are shown in

Table 3. The variable selection result and model test results for the Housing dataset.

Variable Selection Method	Variable Selection Result	RMSE
All variables	X1~X13	9.9566
MIFS [11]	X13, X6, X11, X1, X4	9.6594
mRMR [12]	X13, X10, X12, X6, X7	9.8252
CMIM [13]	X13, X10, X2, X6, X4	10.5333
BIF [24]	X13, X11, X6, X3, X5	9.7156
KNN_MI_CR	X13, X10, X4, X3, X2	9.4404

.

It can be seen from the results that the predictive accuracy of the model is highest by using the KNN_MI_CR method. Therefore, the KNN_MI_CR method contributes to establishing a more effective prediction model.

5. The Field Data Experiment for SCR de-NOx Process

5.1. Field Data and Preselected Input Variables

The inlet NO_x concentration and NH₃ injection flow rate reflect the NH₃/NO_x molar ratio, which influences the de-NO_x efficiency and the NH₃ escape directly. The NH₃ injection flow rate is controlled to match up the change of boiler load by the NH₃ valve. Inlet O₂ content influences outlet NO_x concentration and de-NO_x efficiency directly. Moreover, the change in boiler load often influences the inlet flue gas flow rate. Then, the flue gas temperature is changed through heat exchange. The inlet flue gas temperature influences the reaction rate and catalyst activity.

In a word, the outlet NO_x emission is related to many factors. In addition, the changes in boiler load, coal quality, and combustion condition could cause a large fluctuation in inlet NO_x concentration. Therefore, preselected input variables and the output variable are shown in

Table 4. The preselected input variables and output variable for SCR de-NO_x process.

Variable	Parameter	Unit
Preselected input variables	Boiler load/x1	MW
	Total coal feed rate/x2	t/h
	Total air volume /x3	t/h
	Inlet NOx concentration/x4	mg/m3
	Inlet O2 content/x5	%
	NH3 valve opening/x6	%
	NH3 injection flow rate/x7	kg/h
	Inlet flue gas temperature/x8	°C
	Inlet flue gas flow rate/x9	km3/h
Output variable	Outlet NOx emission/y	mg/m3

.

Assuming that coal quality is invariable, the selected field data cover steady state and variable state. One-dimensional linear interpolation is carried out on the measured NO_x concentration at the process of blowback. Therefore, 2000 samples were collected, and the sampling period is 10 s. Due to the length of this paper, the samples of part preselected input variables and output variables are shown in Figure 6.

5.2. Result Analysis of Delay Estimation and Variable Selection

In general, the maximum time-delay of SCR reaction which includes NO_x analysis is 120~400 s, and the maximum time-delay of boiler load which influences inlet NO_x concentration is 600 s, so the range of τ is (20,60) from x₁ to x₃, and (1,40) from x₄ to x₉. T_max in Equation (32) is 1200 s, so the range of w is (1,120). The time-delay of each input variable is different-in-different operation point. The samples from 200 point to 500 point at 800 MW operating state are taken as an example to analyze the delay estimation and variable selection methods. The size of the SCR reactor is 11.67 m × 13.95 m × 12.6 m. The optimal fitness curves of the PSO algorithm for all preselected input variables are shown in Figure 7, and the results are shown in

Table 5. The results of delay estimation.

Input Variables	KNN_MI	Optimal Result of Time-Delay	Time-Delay
x1	MI(x1, x3)	41	410 + 330 = 740 s
x2	MI(x2, x3)	40	400 + 330 = 730 s
x3	MI(x3, y)	20	200 + 330 = 530 s
x4	MI(x4, y)	33	330 s
x5	MI(x5, y)	23	230 s
x6	MI(x6, y)	24	240 s
x7	MI(x7, y)	40	400 s
x8	MI(x8, y)	40	400 s
x9	MI(x9, y)	14	140 s

.

It can be seen from Figure 7, each preselected variable reaches the optimal global solution before the 10th iteration, and the optimal parameter ${\tau }^{\prime }$ and $w^{\prime }$ are obtained. The optimal result multiplies the sampling period to obtain the time-delay of the preselected variable, which is shown in Table 5.

The effect of time-delay needs to be considered for the variable selection in SCR de-NO_x process. Therefore, based on the delay estimation result, the phase space reconstruction on each input variable is performed first, and then the KNN_MI_CR method is used to select a suitable variable.

It can be seen from

Table 6. The variable selection result and test result for SCR de-NO_x process.

Method	Samples		Reconstructed Samples
	Result	RMSE (TD-KPLS)	Result	RMSE (DTD–KPLS)
All input variable	X1~X9	4.4561 mg/m3	X1~X9	4.2619 mg/m3
MIFS [11]	X3, X4, X6	3.2071 mg/m3	X3, X4, X8	3.2089 mg/m3
mRMR [12]	X3, X5, X4	3.2101 mg/m3	X3, X5, X2	3.2168 mg/m3
CMIM [13]	X3, X7, X4	3.2061 mg/m3	X3, X7, X4	3.2003 mg/m3
BIF [24]	X3, X4, X7	3.3069 mg/m3	X3, X4, X7	3.2207 mg/m3
DISR [25]	X3, X4, X7	3.3069 mg/m3	X3, X4, X2	3.2130 mg/m3
KNN_MI_CR	X3, X5, X4	3.2101 mg/m3	X3, X7, X4	3.2003 mg/m3

that if the model samples are reconstructed by phase space, the result of variable selection is different. The DTD–KPLS model, whose variables are selected by the KNN_MI_CR method, has the highest predictive accuracy. Although the variable selection result is the same, the predictive accuracy of the DTD–KPLS model whose samples are reconstructed is improved in the BIF and CMIM method. In addition, the sequence of the selected input variable is different, and the model accuracy is different, which indicates that the consideration of high-dimensional KNN_MI is conducive to the selection of suitable input variables.

5.3. Analysis of Dynamic Model

In general, dynamic modeling is realized by adding the historical input data $x(t-1),\cdots ,x(t-w+1)$ and output data $y(t-1),\cdots ,y(t-w+1)$ as a new input of the model. The 2000 prediction samples are used to analyze model accuracy by using different input variables. The results are shown in

Table 7. The predictive accuracy of DT-KPLS and DTD–KPLS model using different input variables.

Model	Input Variables	Dimension of Input Variable	RMSE
DT-KPLS	$x(t)$	3	3.2268 mg/m3
	$y(t-1),x(t)$	4	3.2166 mg/m3
	$x(t-1),x(t)$	6	3.2054 mg/m3
	$x(t-1),y(t-1),x(t)$	7	3.2061 mg/m3
DTD–KPLS	$x(t-\tau )$	3	3.2139 mg/m3
	$y(t-1),x(t-\tau )$	4	3.2082 mg/m3
	$x(t-\tau -1),x(t-\tau )$	6	3.2035 mg/m3
	$x(t-\tau -1),y(t-1),x(t-\tau )$	7	3.2003 mg/m3

.

It can be seen from Table 7 when input variables are $x(t-\tau -1)$ , $y(t-1)$ and $x(t-\tau )$, the predictive accuracy of the model is highest. In addition, because the model input is historical input data before the current time t, the DTD–KPLS model can predict the current value of the output variable in advance accurately.

When parameter optimization is used in model control strategies, a large number of calculations need to be avoided. In this paper, the DTD–KPLS model has only one RBF kernel parameter. The predictive accuracy of the DTD–KPLS model in different kernel parameters is analyzed. The results are shown in Figure 8.

The change rate of the RMSE value was no more than 1%, as shown in Figure 8, which indicates that the predictive accuracy was not sensitive to the kernel parameter. Therefore, the DTD–KPLS model could adopt fixed parameters to avoid a large number of calculations in parameter optimization. In this paper, the single-step prediction time was 1 s, so it could meet the demand for control.

In order to compensate for the prediction error, this paper proposes a feedback correction strategy for the DTD–KPLS model. The predicted value and the corrected predicted value are shown in Figure 9.

It can be seen from Figure 9, and the DTD–KPLS model shows high predictive accuracy with a low RMSE of 3.2933 mg/m³ and a high determinate coefficient (Q²) of 0.7011 in the whole state. A large amount of NO_x was produced in a variable state. Hence the predictive accuracy was lower at the peak of the outlet NO_x emission curve. After feedback correction, the predicted values were closer to the real value. Steady state and variable state were covered in the model samples, so the DTD–KPLS model had strong self-adaptability.

5.4. Analysis of Control Effect

Two kinds of NH₃ injection control modes have mainly existed in China. One is fixed NH₃/NO_x molar ratio control, and the other one is fixed outlet NO_x value control. In this paper, the fixed outlet NO_x value control was analyzed. First, the NO_x flow rate was determined by the inlet NO_x concentration and the flue gas flow rate, and the preset molar ratio was calculated. Then, the NH₃/NO_x molar ratio was corrected based on the measured value and the set value of outlet NO_x emission. The amount of NH₃ injection was determined by the corrected NH₃/NO_x molar ratio and the calculated NO_x flow rate. Finally, the fixed outlet NO_x emission was achieved.

In order to obtain the demand amount of NH₃ injection $y_{\mathrm{N}{\mathrm{H}}_3}^{\prime }$ at the current operating state, the transfer function from the NH₃ injection to the decreased amount of outlet NO_x emission was investigated through a dynamic characteristic test as follows:

$$W(s)=-\frac{1.73e^{-21s}}{10236s^2+196s+1}$$

(39)

Then, the control structure diagram set up in MATLAB/SIMULINK is shown in Figure 10.

In the field data in Figure 6, the inlet NO_x and the flue gas flow rate fluctuated from 135~230 mg/m³ to 120~160 km³/h. During the simulation, the inlet NO_x, the flue gas flow rate and the set value of outlet NO_x were set to 160 mg/m³, 1300 km³/h and 37 mg/m³, which were consistent with the field data. PID tuning program was used to get optimal PID parameters. The results are shown in Figure 11, Figure 12 and Figure 13.

If the inlet NO_x value and the set value of outlet NO_x were fixed, as the deviation between them decreased, the PID controller reduced the actual NH₃/NO_x molar ratio; thereby, the amount of NH₃ injection was reduced gradually. After the system was stabilized, the outlet NO_x emission, the NH₃ injection flow rate and the NH₃/NO_x molar ratio were stable at 37 mg/m³(set value), 60 kg/h and 0.8857, which were close to the average value of field data. So the demand amount of NH₃ injection $y_{\mathrm{N}{\mathrm{H}}_3}^{\prime }$ corresponding to the set value of outlet NO_x emission was 59.6 kg/h, and this method could accurately simulate the fixed outlet NO_x control mode. The control process lasted for nearly 1000 s before the system became stable. Hence, the existing system could not meet real-time control because of the long adjustment time.

The predicted value of the TD-KPLS model was taken as the output of NIC in each model update. The input variables included the deviation between the corrected predicted NO_x emission value and the set value of outlet NO_x emission, inlet flue gas temperature and inlet NO_x concentration. The output variable was the amount of NH₃ injection compensation. The output of the NIC is shown in Figure 14.

Figure 14 shows that the output of the NIC was consistent with the trend of the expected value of NH₃ injection compensation. Therefore, the predictive accuracy of TD-KPLS was high. In addition, the predicted value was a mostly negative value indicating that the amount of NH₃ injection was reduced by 11% after adding the NIC. Assuming that the price of NH₃ water was about 3000 RMB per ton, the data used in this paper was taken as an example. After the application of the NIC, the cost of de-NO_x could be saved by 531.2 RMB per day.

Because the SCR de-NO_x system had large inertia and had a large delay in NO_x analysis; hence, the following analysis started from the 100th point. The control effect of outlet NO_x emission is shown in Figure 15.

Figure 15 shows the control effect of the existing controller was not ideal, and the outlet NO_x emission produced drastic fluctuation when the load changes. For example, the outlet NO_x emission exceeded the national NO_x emission standard (50 mg/m³) at the 500, 1400 and 1800 points. In contrast, after the NIC was added, the outlet NO_x emission whose fluctuation ranged in ±5 mg/m³ was stable near the set value. Therefore, the outlet NO_x emission was far lower than the national NO_x emission standard.

It can be seen from Figure 16 that the minimum value of the de-NO_x efficiency after adding the NIC (76.31%) was higher than the minimum value of the de-NO_x efficiency of the existing controller (73.79%). The average de-NO_x efficiency was 80.12% by using the existing controller, and the average de-NO_x efficiency was 80.15% after adding the NIC. This indicated that the de-NO_x efficiency was basically unchanged after adding the NIC.

In order to calculate the NH₃ escape further, the following formula was introduced.

$$p(t)=\frac{(W_m(t)-W_c(t))\times {10}^6}{17\times V_{gas}(t)}$$

(40)

Here, t represents the sampling time; W_m is the measured NH₃ injection flow rate, kg/h; W_c is the NH₃ injection consumption, kg/h; V_gas is the inlet flue gas flow rate, m³/h; p is the NH₃ escape, μg/m³.

Because the compensation amount of NH₃ injection in this paper was negative value mostly, the consumption of NH₃ injection was regarded unchanged basically. Figure 17 was the comparison between the NH₃ escape of the existing controller and the NH₃ escape after NH₃ injection compensation.

In Figure 17, after the NIC was added, although the average of de-NO_x efficiency did not change, the NH₃ escape was reduced by 39% after NH₃ injection compensation.

Actually, in order to avoid actuator saturation, the NH₃ injection valve opening was restrained by conditions (23) and (24) in the NH₃ injection control. The NH₃ valve opening of the existing controller and NH₃ valve opening after NH₃ injection compensation are shown in Figure 18.

Due to the lag of NO_x analysis, the NH₃ control valve was prevented from being over-opened or over-closed. Therefore, when the NH₃ control valve was in automatic mode, the NH₃ valve opening was limited from U_Vmin = 20% to U_Vmax = 50%. In Figure 18, the NH₃ valve opening of the existing controller reached 50% from 1300 points to 1500 points, which indicated that the actuator reached saturation in automatic control. When the NIC was used, the average of the NH₃ valve opening was reduced, thereby actuator saturation was avoided.

The upper limitation u_Vmax and the lower limitation u_Vmin of the increment of NH₃ valve opening were 2% and −2%, respectively. Figure 19 shows that the increment of NH₃ valve opening after NH₃ injection compensation was within limitation. Therefore, the NH₃ valve was protected, and the control accuracy was improved.

In all, the proposed control strategy had a good control effect. Through the NH₃ injection was adjusted by NIC continually, the controller could decrease the NH₃ injection and decrease the NH₃ escape as much as possible while ensuring that the high de-NO_x efficiency and NO_x emission meets the national NO_x emission standard. Therefore, the operation cost was reduced, and secondary pollution caused by the incomplete reaction of NH₃ was avoided.

6. Conclusions

A dynamic KPLS model incorporated with variable selection and delay estimation was proposed for the prediction of outlet NO_x emission in advance. The KNN_MI method was used to estimate the time-delay, and the effect of historical data lengths on the KNN_MI value was considered. Bidirectional search, which is based on the change rate of KNN_MI (KNN_MI_CR), was used for variable selection, which avoids using a single evaluation strategy and considers the effect of high-dimensional KNN_MI. Finally, the DTD update method and feedback correction strategy are proposed, so the fixed parameter is used in the final model. Therefore, the final model could be identified offline by using field data and have high predictive accuracy.

The proposed control system includes the main controller and the correction controller. The main controller is the existing PID controller. The combination of NIC and the outlet NO_x emission model is used as the correction controller. According to the results of the field data experiment, the proposed control strategy has several findings as follows:

(1)

In the aspect of the control effect, the outlet NO_x emission meets the national NO_x emission standard, whose fluctuation range in ±5 mg/m³ is stable near the set value, and the system maintains high de-NO_x efficiency; the NH₃ escape is reduced by 39% after NIC.

(2)

Compared with the conventional PID control, the proposed control strategy not only protects the NH₃ injection valve but also reduces the operation cost because of the reduction of NH₃ injection;

(3)

The reduction of NH₃ escape avoids operating accidents and secondary pollution.

The potential disadvantage of the proposed control strategy is that the calculation of the demand amount of NH₃ injection $y_{\mathrm{N}{\mathrm{H}}_3}^{\prime }$ at the current operating state involves a dynamic characteristic test. For different operating states, $y_{\mathrm{N}{\mathrm{H}}_3}^{\prime }$ will be different. If the operating state changes frequently, it will not be easy to determine $y_{\mathrm{N}{\mathrm{H}}_3}^{\prime }$ accurately.