# Artificial Bee Colony Optimization of NOx Emission and Reheat Steam Temperature in a 1000 MW Boiler

Authors: Xian-hua Gao; Zhi-gang Su
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

This paper puts forward a new viewpoint on optimization of boiler combustion, namely, reducing NOx emission while maintaining higher reheat steam temperature rather than reducing NOx emission while improving boiler efficiency like traditional practices. Firstly, a set ofmultioutputs nonlinear partial least squares (MO-NPLS) models are established as predictors to predict these two indicators. To guarantee better predictive performance, repeated double cross-validation (rdCV) strategy is proposed to identify the structure as well as parameters of the predictors. Afterward, some controllable process variables, taken as inputs of the predictors, are then optimized by minimizing NOx emission and maximizing reheat steam temperature via multiobjective artificial bee colony (MO-ABC). Results show that our rdCV-MO-NPLS model with MO-ABC optimization methods can reduce NOx emission synchronously and improve reheat steam temperature effectively compared with nondominated sorting genetic algorithm II (NSGA-II) and combustion adjustment experimental data on a real 1000 MW boiler.


---

## Body

## 1. Introduction

Boiler combustion optimization of reducing nitrogen oxide (NOx, for short) has become a hot topic during the past three decades. There exist several currently typical ways in coping with the boiler combustion optimization problem. The first popular method is numerical simulation that is based on computational fluid dynamics (CFD) [1–3]. However, the optimal set points explored by numerical simulation cannot guarantee optimal value in practice due to inconsistencies between designed conditions in CFD and practical conditions. To consider this, it is interesting to take into account the abundant information originating from practical running data. Under such circumstance, the data-driven method is motivated.The development of data-driven methods can be divided into two stages. At the first stage, researchers only focus on NOx emission reduction by using evolutional algorithms (EAs) to search the optimal set point of controllable process variables. The typical EAs include genetic algorithm (GA), particle swarm optimization (PSO), and ant colony optimization (ACO). Among them, GA is frequently used on boiler combustion optimization. For example, in [4, 5], artificial neural network (ANN) combined with GA is proposed to reduce NOx emission. And in [6], support vector machine (SVM) and conventional GA (CGA) are used on NOx emission reduction. Researchers can notice that only NOx emission is considered as the optimization objective in the above researches, but ignoring that the combustion system is a hybrid system which is coupled with other subsystems. The reduction of NOx usually plays great influence on other important parameters such as reheat steam temperature (Tr, °C). That is, the optimization of boiler combustion system is essentially a multiobjective optimization problem.This issue with multiobjective optimization in boiler combustion systems has gained great attention at the second stage. Traditionally, boiler combustion optimization is implemented to reduce NOx emission and increase boiler efficiency by EAs (for example, see works [7–10]). ANNs with multiobjective GA (MOGA) are employed on boiler combustion optimization. Although the above methods have made significant contributions on reducing NOx emission while improving boiler efficiency, there still exist some unconsidered issues. Firstly, reheat steam temperature is always lower than its set value due to the requirement of meeting electric peak regulation. This phenomenon not only imposes a bad influence on the whole plant efficiency but also poses a threat to safe control of the boiler combustion system. Secondly, the precise calculation of boiler efficiency is hard to complete actually due to some practical issues. The third problem is the lack of effective modelling data, as historical running data are heavily coupled and incomplete. Therefore, it is urgent to improve reheat steam temperature.Up to date, to the best of our knowledge, there is few literature on boiler combustion optimization to reduce NOx emission and improve reheat steam temperature simultaneously. Motivated by the above statements, this paper aims to reduce NOx emission while achieving higher reheat steam temperature by proposing a jointed optimization method. More precisely, a Multiple Output version of Nonlinear Partial Least Squares (MO-NPLS) regression model is first proposed to predict behaviors of NOx emission and reheat steam temperature synchronously. To achieve better performance, the structure and parameters of MO-NPLS are identified by repeated double cross validation (rdCV) rather than traditional k-fold cross-validation strategy. Then, a multiobjective artificial bee colony (MO-ABC) algorithm is adopted to search optimal inputs of predictors that can achieve lower NOx emission and higher reheat steam temperature from Pareto front. The modelling data are from a combustion adjustment experiment conducted on a real 1000 MW coal-fired boiler. It will be seen that our proposed method outperforms combustion adjustment experimental data and commonly used nondominated sorting genetic algorithm II (NSGA-II). Therefore, the main contributions of this paper can be summarized in three-fold:(1)
Both NOx emission and reheat steam temperature are firstly proposed in boiler combustion optimization. This strategy can reduce NOx emission while maintaining higher reheat steam temperature to guarantee the safety in boiler combustion system control as well as the improvement of whole plant efficiency.(2)
A MO-NPLS method with rdCV strategy is raised to establish the predictive model of both NOx emission and reheat steam temperature. This algorithm is powerful due to its strong generalization ability especially with small number of samples.(3)
A MO-ABC-based multiobjective optimization method is proposed on boiler combustion optimization. On the basis of the rdCV-MO-NPLS modelling method, our proposed MO-ABC algorithm can reduce NOx emission synchronously and improve reheat steam temperature effectively compared with NSGA-II and the experimental data.The rest of the paper is organized as follows. In Section2, the 1000 MW boiler combustion system and data sources are described. Section 3 establishes predictors used to predict behaviors of NOx emission and reheat steam temperature by rdCV-MO-NPLS. Then, Section 4 presents the joint optimization method by MO-ABC and the comparison with NSGA-II and experimental data. Section 5 concludes this paper.


## 2. Description of Combustion System

In this section, the boiler combustion system is introduced, together with the dynamic characteristics of combustion system. And then, the data source used for modelling establishment and validation is presented.


### 2.1. Boiler Combustion System

Figure1 illustrates a 1000 MW tower-type ultrasupercritical boiler adopting the reheat system and air-staged combustion technique. As can be seen from Figure 1, a selective catalytic reduction (SCR) system is installed at the exit of convection flue gas and two air preheaters are arranged beneath the SCR system. Along vertical direction, the entire furnace is divided into three zones including main combustion zone, reduction zone, and burning-out zone. Furthermore, in the main combustion zone, twelve-level pulverized coal burners are equipped and correspond to six mills indexed by letters from A to F that are with burners A1 and A2 at the bottom to F1 and F2 at the top, respectively.Figure 1
Schematic diagram of a 1000 MW boiler combustion system: (1) boiler; (2) burners.Generally, different coal allocations decided by the coal feeding rate of each mill will lead to different temperature fields and flame centers in the furnace. For instance, if more coal is supplied by bottom mills, the flame center will move downward. As a result, NOx production will reduce due to the extension of reductive zone; however, the reheat steam temperature (Tr) will decrease too because of less quality of heat. In the present work, although the coal feeding rate is a directly manipulated parameter, thermal load (Qi) is monitored instead by considering the diversity of coal types as follows:(1)Qi=BiB0⋅Q0,i∈A,B,…,F,where B0 is the total coal fed by mills and Q0 indicates the total thermal load.In contrast to coal allocations, primary air is usually not adjusted because its main task is to transport pulverized power from mills, i.e., to match the pulverized coal. Secondary air, including that in the main combustion zone and separated overfire air (SOFA) in the burning-out zone, can be adjusted to improve combustion conditions in the boiler. To enhance reducibility in the main combustion zone so as to reduce NOx production, one can decrease secondary air in the main combustion zone with increase of SOFA ratio (αSOFA) as the total air volume remains constant. Evidently, this parameter represents the ratio of secondary air to total air. More interestingly, higher reheat steam temperature can be achieved with higher αSOFA, but excessive αSOFA will affect the combustion performance in the main combustion zone. However, the volume of secondary air in the main combustion zone is rather difficult to monitor due to tens of dampers in power plant. Instead, αSOFA is used and defined as(2)αSOFA=SOFAtotal air,where total air = primary air + secondary air. Considering that there is about 20% running oxygen (O2) in total air. Therefore, the higher the running O2 is, the greater the total air volume will be.In addition, notice that the higher the runningO2 value is, the better combustion effect will be. Then, both NOx emission and reheat steam temperature will increase subsequently and vice versa. This means there exists an optimal α value, namely, optimal running O2.In summary, to reduce NOx emission and increase reheat steam temperature synchronously heavily relies on the distribution of thermal loads (Qi), running O2, and αSOFA in the furnace.Remark 1.
The burners at vertical direction are tilted within a range of ±30° to tune the reheat steam temperature (Tr,°C) [11]. At low load condition, even if burners have been already fully titled up, the reheat steam temperature will be still significantly lower than its desirable setpoint. This will impose a bad influence on whole power plant efficiency (the whole plant efficiency ηcp=ηb⋅ηp⋅ηi⋅ηm⋅ηg where ηi is plant thermal efficiency which depends on reheat steam temperature).


### 2.2. Data Source

The value of reheat steam temperature from actual operational conditions is not accurate as it is controlled by a closed-loop temperature system. Therefore, we adopted the data collected from a combustion adjustment experiment shown in Supplementary file of the 1000 MW power plant described in Section2.1. Consider the problem of excessive NOx emission and lower reheat steam temperature is more severe under lower loads. As a result, only 650 MW, 800 MW, and 850 MW conditions are discussed in the following. Corresponding data are listed in Appendix A (Tables A1–A3) in the Supplementary file.Remark 2.
Reheat steam temperature in the experiment was not controlled by a closed-loop way, and hence, its values are real. Therefore, all collected experimental data are more useful than data collected in practical running way when validating the performance of our method that will be presented later.


## 2.1. Boiler Combustion System

Figure1 illustrates a 1000 MW tower-type ultrasupercritical boiler adopting the reheat system and air-staged combustion technique. As can be seen from Figure 1, a selective catalytic reduction (SCR) system is installed at the exit of convection flue gas and two air preheaters are arranged beneath the SCR system. Along vertical direction, the entire furnace is divided into three zones including main combustion zone, reduction zone, and burning-out zone. Furthermore, in the main combustion zone, twelve-level pulverized coal burners are equipped and correspond to six mills indexed by letters from A to F that are with burners A1 and A2 at the bottom to F1 and F2 at the top, respectively.Figure 1
Schematic diagram of a 1000 MW boiler combustion system: (1) boiler; (2) burners.Generally, different coal allocations decided by the coal feeding rate of each mill will lead to different temperature fields and flame centers in the furnace. For instance, if more coal is supplied by bottom mills, the flame center will move downward. As a result, NOx production will reduce due to the extension of reductive zone; however, the reheat steam temperature (Tr) will decrease too because of less quality of heat. In the present work, although the coal feeding rate is a directly manipulated parameter, thermal load (Qi) is monitored instead by considering the diversity of coal types as follows:(1)Qi=BiB0⋅Q0,i∈A,B,…,F,where B0 is the total coal fed by mills and Q0 indicates the total thermal load.In contrast to coal allocations, primary air is usually not adjusted because its main task is to transport pulverized power from mills, i.e., to match the pulverized coal. Secondary air, including that in the main combustion zone and separated overfire air (SOFA) in the burning-out zone, can be adjusted to improve combustion conditions in the boiler. To enhance reducibility in the main combustion zone so as to reduce NOx production, one can decrease secondary air in the main combustion zone with increase of SOFA ratio (αSOFA) as the total air volume remains constant. Evidently, this parameter represents the ratio of secondary air to total air. More interestingly, higher reheat steam temperature can be achieved with higher αSOFA, but excessive αSOFA will affect the combustion performance in the main combustion zone. However, the volume of secondary air in the main combustion zone is rather difficult to monitor due to tens of dampers in power plant. Instead, αSOFA is used and defined as(2)αSOFA=SOFAtotal air,where total air = primary air + secondary air. Considering that there is about 20% running oxygen (O2) in total air. Therefore, the higher the running O2 is, the greater the total air volume will be.In addition, notice that the higher the runningO2 value is, the better combustion effect will be. Then, both NOx emission and reheat steam temperature will increase subsequently and vice versa. This means there exists an optimal α value, namely, optimal running O2.In summary, to reduce NOx emission and increase reheat steam temperature synchronously heavily relies on the distribution of thermal loads (Qi), running O2, and αSOFA in the furnace.Remark 1.
The burners at vertical direction are tilted within a range of ±30° to tune the reheat steam temperature (Tr,°C) [11]. At low load condition, even if burners have been already fully titled up, the reheat steam temperature will be still significantly lower than its desirable setpoint. This will impose a bad influence on whole power plant efficiency (the whole plant efficiency ηcp=ηb⋅ηp⋅ηi⋅ηm⋅ηg where ηi is plant thermal efficiency which depends on reheat steam temperature).


## 2.2. Data Source

The value of reheat steam temperature from actual operational conditions is not accurate as it is controlled by a closed-loop temperature system. Therefore, we adopted the data collected from a combustion adjustment experiment shown in Supplementary file of the 1000 MW power plant described in Section2.1. Consider the problem of excessive NOx emission and lower reheat steam temperature is more severe under lower loads. As a result, only 650 MW, 800 MW, and 850 MW conditions are discussed in the following. Corresponding data are listed in Appendix A (Tables A1–A3) in the Supplementary file.Remark 2.
Reheat steam temperature in the experiment was not controlled by a closed-loop way, and hence, its values are real. Therefore, all collected experimental data are more useful than data collected in practical running way when validating the performance of our method that will be presented later.


## 3. Estimations of NOx and Reheat Steam Temperature Based on MO-NPLS with rdCV

PLS/NPLS algorithms have been applied in engineering (for example, see [12–14]). However, due to the limited number of samples, the prediction accuracy of PLS/NPLS based on traditional k-fold cross validation may decrease. Motivated by this and literature [15], we propose a multiple output version of NPLS (MO-NPLS) model with repeated double cross-validation (rdCV) strategy in Section 3.1. The predictive model is necessarily used to predict the behaviors of NOx emission and reheat steam temperature and then identify and validate it in different cases of loads in Section 3.2.


### 3.1. MO-NPLS Model with rdCV

Give inputsxi=xi1,…,ximT∈Rm and outputs yi=yi1,yi2T∈R2 with i=1,…,n. Then, let X=x1,…,xnT∈Rn×m and Y=y1,…,yn∈Rn×2 be the number of n observations for inputs and outputs, respectively. Evidently, Xij indicates the i‐th observation of the j‐th input, and this is held on for output Y. The superscript “T” represents transpose. To make the results independent from units, all attributes of data X,Y are normalized into range −1,1 according to the z-score rule and are denoted by E0,F0 with E0∈Rn×m and F0∈Rn×2.Definetj∈Rn and uj∈Rnj=1,…,m are the principal components of X and Y, respectively. Similarly as traditional nonlinear PLS [16, 17], rdCV-MO-NPLS is in nature a recursive modelling strategy, aiming to establish relationship between tj and ujj=1,…,m firstly and then implement the regression relationship of X and Y.Letu1 be any column of F0 at the beginning. Then, one can estimate of the first components (t1, u1) as follows.Firstly,t1 is calculated by(3)t1=E0w1w1Tw1,where w1∈Rm×1 is the weight of matrix E0 and obtained by w1T=u1TE0/u1Tu1.Then, one can fit the polynomial relationship betweent1 and u1 by the least squares method. Specifically, assume the relationship between t1 and u1 is quadratic, namely,(4)u1=ft1+e1=c1,0I+c1,1t1+c1,2t12+e1,where e1∈Rn is the residual vector which is independent and identically distributed. I∈Rn is a vector of ones. Take Z=I,t1,t12∈Rn×3 and c1=c1,0,c1,1,c1,2T∈R3, then (4) can be rewritten by a compact form:(5)u1=Zc1+e1.So, one can get the polynomial coefficientc1 by the least squares method:(6)c1=ZTZ−1ZTu1,where (ZTZ) represents the inverse matrix of ZTZ.After coefficientsc1 are fitted by formula (6), then r1, i.e., the estimate of u1, can be calculated by(7)r1=u^1=ft1=c1,0I+c1,1t1+c1,2t12.Then, updateu1 (mark u1,new) by(8)u1,new=F0q1q1Tq1,where q1∈R2×1 are weights of matrix F0 and determined by q1T=r1TF0/r1Tr1. Next, w1 needs to be updated according to [18].At this point, we can recalculate the vectort1 (mark t1,new) by formula (3) with the updated w1.Check the convergence oft1: if t1,new−t1/t1≥10−6, refit the polynomial coefficient c1 between u1,new and t1,new. Otherwise, calculate the regression coefficient p1 by(9)p1T=t1,newTE0t1,newTt1,new.Finally, the residual ofE0 and F0 (mark E1 and F1) that will be used in the next iteration can be induced by(10)E1=E0−t1,newp1T,F1=F0−r1q1T,with regression coefficient pj.ReplaceE0 and F0 with E1 and F1 and formulas (3)–(10) are repeated until all m principal components (tj,uj) with j=1,…,m are extracted.After allm principal components have been extracted, optimum number afinal (afinal≤m) can be selected. Usually, k-fold strategy is applied, whereas rdCV strategy is adopted here for achieving more robust predictive ability, which has been proven of PLS in [15]. This is the reason why we call our model rdCV-MO-NPLS. Once afinal is determined, the prediction of outputs Y can thus be obtained according to(11)Y^=Y¯+sY,∑k=1afinalrkqkT,k=1,…,afinal,where <⋅,⋅> means point-product operation of two vectors with equal dimensionality; Y¯=Y¯1,Y¯2 with Y¯j=1/n∑i=1nYij,j=1,2; and sY=sY1,sY2 with sYj=1/n∑i=1nYij−Y¯j2.From the above interpretations, the rdCV-MO-NPLS is summarized in Algorithm1.Algorithm 1: rdCV-MO-NPLS.
Input: dataX,Y and parameters Nrep,SEGtest,SEGcalib;
Output: optimal number of principal componentafinal and rdCV-MO-NPLS model.
for rep=1 to Nrep do
Split datasetS equally into SEGtest groups, i.e., S=S1,S2,…,SSEGtest;
for i=1 to SEGtest do
SelectSi as Test set, and construct Calibration sets S−i = S−Si;
Split Calibration set equally intoSEGcalib groups, S−i=V1,V2,…,VSEGcalib;
for j=1 to SEGcalib do
SelectVj as Validation set and V−j=S−i−Vj as Training set;
Fit MO-NPLS based onV−j with, respectively, a=1,2,…,m of principal components,
Apply the MO-NPLS models toVj and get predictive Y^ja for a=1,2,…,m;
Calculate mean square errorMSEja=1/nMSE∑k=1nMSEYjk−Y^jak2, with the number of objects in the used validation set    nMSE and the output objects of validation set Yj;
Estimate optimum principal componentsaopt1i, …, aoptki according to Y^a based on standard error method [15] (Here, more than   one aopt may be selected (i.e., k≥1) because of different confidence interval);
Make MO-NPLS models based on Calibration setS−i with a=aopt1i,…,aoptki;
Test fitted models onSi and obtain a group of k predictions Y^ka as well as k biases;
Find the smallest bias and determine the optimal principal componentaopti,rep.
One can getaopt with the number of SEGtest after completing the outer loop.
Totally, after a complete rdCV run, we can getaopt with the number of (Nrep×SEGtest). The final optimum of principal component afinal is the one with highest frequency in aopti,rep.
Identify rdCV-MO-NPLS model with databaseS and afinal to get model parameters [w, p, c, q] based on (3)–(10).


### 3.2. Model Identification and Validation

#### 3.2.1. Model Identification

In this section, rdCV-MO-NPLS is identified, respectively, at 650 MW, 800 MW, and 850 MW. In each case, identification as well as rdCV data are shown in TableA1 (from lines 1 to 45), Table A2 (lines 1 to 30), and Table A3 (lines 1 to 90). In cases of 650 MW and 850 MW, six inputs were considered, including C-F level thermal loads (X1 to X4, MW), running O2 (X5), and αSOFA (X6, %), whereas seven inputs were considered in 800 MW, which are B-F level thermal loads (X1 to X5, MW), running O2 (X6), and αSOFA (X7, %). In all cases, the outputs are(12)Y=NOxX,TrX.Figure2 displays the results of applying rdCV. It can be seen that the highest percentage in three cases are, respectively, 0.32 (aopt=3), 0.38 (aopt=2), and 0.256 (aopt = 2). Hence, the final optimum principal component afinal in each case is 3, 2, and 2, respectively. With the help of afinal, the rdCV-MO-NPLS model can be identified as Case 1 to Case 3.Figure 2
Percentage of principal componentsaopt: (a) 650 MW, (b) 800 MW, and (c) 850 MW.
(a)
(b)
(c)For better comparison, 5-fold cross validation (We call 5-fold-MO-NPLS model) is conducted by the same identification data as the rdCV method. After 5-fold cross validation process, the optimum principal componentafinal of 650 MW, 800 MW, and 850 MW cases is 4, 1, and 4, respectively. Then, the validation results of 5-fold-MO-NPLS model with their afinal will be seen in the next section.Case 1.
(rdCV-MO-NPLS of 650 MW load condition). Sinceafinal = 3, the first three components (t1, u1), (t2, u2), and (t3, u3) are selected. Normalize the training data presented in Table A1 (from lines 1 to 45) into E0=e1,…,e6 and F0 according to z‐score rule. Corresponding parameters of (w,p,c,q) shown in Table 1 as well as (t1, u1), (t2, u2), and (t3, u3) are identified based on (3)–(11). The rdCV-MO-NPLS predictive model is in form of(13)t1=E0w1=E00.1084,0.7642,−0.2669,−0.5667,−0.0943,0.0555T,t2=E1w2=E0I−w1p1Tw2=E0−0.2094,−0.2329,−0.6923,−0.7341,−0.1408,0.1400T,t3=E2w3=E0I−w1p1TI−w2p2Tw3=E00.8325,−0.4500,0.1520,1.2745,0.0707,−0.0796T,f=r1q1T+r2q2T+r3q3T=−0.1112+1.9943t1+0.6504t12+0.626t2+0.3445t22+0.2332t3−0.4783t32,−0.0403+1.9161t1+0.6248t12−0.9299t2−0.5116t22+0.1726t3−0.3540t32≔f1,f2,Y^=Y^1,Y^2=Y¯1,Y¯1+sY1,sY2,f1,f2=509.5+55.8445f1,581.2+6.1878f2.Table 1
Parameters of rdCV-MO-NPLS model at 650 MW.
h
1
w
1
p
1
c
1
q
1
0.1084
−1.0821
−0.2037
0.7211
0.7642
1.7506
2.7657
0.6928
−0.2669
−2.0513
0.9019
−0.5667
1.7544
−0.0943
−1.3794
0.0555
1.7356
h
2
w
2
p
2
c
2
q
2
−0.2521
1.3719
−0.0598
0.5585
−0.5344
−1.1655
1.1211
−0.8295
−0.5870
0.1064
0.6167
−0.5105
−1.2466
0.1779
0.0342
0.1181
1.2109
h
3
w
3
p
3
c
3
q
3
0.7658
1.7521
0.0860
0.8038
−0.0421
−0.7898
0.2901
0.5949
−0.5367
0.0805
−0.5950
0.3318
−0.6007
0.0995
−1.1659
0.0608
−0.2728Case 2.
(rdCV-MO-NPLS of 800 MW load condition). The first two components (t1, u1) and (t2, u2) are chosen. Based on data as shown in Table A2 (lines 1 to 30), and identified parameters of rdCV-MO-NPLS shown in Table 2, the rdCV-MO-NPLS model can be established as(14)t1=E0w1=E0−0.4416,−0.3083,0.3598,−0.1139,−0.1758,−0.5304,−0.5053T,t2=E1w2=E0I−w1p1Tw2=E00.4468,0.0730,−0.3284,−0.0332,0.3795,0.4261,0.6668T,f=r1q1T+r2q2T=−0.1626+1.104t1+0.2081t12+0.2959t2+0.0180t22,−0.1275+0.9698t12+0.1828t1−0.3007t2−0.0183t22≔f1,f2,Y^=Y^1,Y^2=Y¯1,Y¯2+sY1,sY2,f1,f2=451.9+27.3115f1,585.7+4.7006f2.Table 2
Parameters of rdCV-MO-NPLS model at 800 MW.
h
1
w
1
p
1
c
1
q
1
−0.4416
−0.0237
−0.2057
0.7513
−0.3083
−0.2094
1.4695
0.6600
0.3598
0.0880
0.2770
−0.1139
0.2170
−0.1758
−0.0283
−0.5304
−1.0054
0.5053
0.7517
h
2
w
2
p
2
c
2
q
2
0.3187
0.9644
−0.0116
0.7014
−0.0164
0.1064
0.4219
−0.7128
−0.2241
−0.2031
0.0257
−0.0662
−0.4456
0.3286
−0.9193
0.2724
0.3618
0.8133
1.0018Case 3.
(rdCV-MO-NPLS of 850 MW load condition). Similar as Case 2, the first two components (t1, u1) and (t2, u2) were selected on the basis of Table A3 (from lines 1–90). Additionally, identified parameters of rdCV-MO-NPLS are shown in Table 3. The rdCV-MO-NPLS model is(15)t1=E0w1=E00.2595,0.4561,−0.3307,−0.5786,−0.0490,0.5273T,t2=E0I−w1p1Tw2=E00.4214,0.4823,0.2326,−1.2139,0.0331,0.0308T,f=r1q1T+r2q2T=0.0143+1.3261t1+0.0230t12+0.5570t2−0.0635t22,0.0049+0.7791t1−0.0135t12−0.9594t2+0.1095t22≔f1,f2,Y^=Y^1,Y^2=Y¯1,Y¯2+sY1,sY2,f1,f2=464.2+35.8762f1,594.8+2.9182f2.Table 3
Parameters of rdCV-MO-NPLS model at 850 MW.
h
1
w
1
p
1
c
1
q
1
0.2595
0.3142
0.0113
0.8622
0.4561
0.0651
1.5380
0.5066
−0.3307
−0.6749
−0.0266
−0.5786
0.0298
−0.0490
−0.2947
0.05273
1.2676
h
2
w
2
p
2
c
2
q
2
0.1698
−1.6923
0.0092
0.5021
0.0401
−2.2468
1.1094
−0.8648
0.5532
−0.5951
−0.1266
−0.6529
−2.4330
0.0806
−0.2785
−0.4804
−0.2926


#### 3.2.2. Model Validation

In general, generalization ability of a prediction model is more important than its fitting ability. Therefore, it is necessary to validate aforementioned models to get their generalization ability. Validation samples in three cases are seen in TableA1 (from lines 46 to 53), Table A2 (from lines 31 to 36), and Table A3 (from lines 91 to 96) in Supplementary file, respectively.To state the superiority of the rdCV-MO-NPLS model, another three well-known methods: 5-fold-MO-NPLS, support vector regression (SVR), and artificial neural network (ANN) are employed to establish models of NOx emission and reheated steam temperature. All methods are on the same identification and validation data. Figures3–5 show the predictive results for different methods. The prediction error sum (PRESS=∑i=1nPRESSY^i−Yi/Yi with prediction sample number nPRESS) representing generalization ability and simulation error sum (SS=∑i=1nSSY^i−Yi/Yi with simulation sample number nPRESS) reflecting approximate ability is presented in Table 4.Figure 3
Comparisons among different predictors at 650 MW.Figure 4
Comparisons among different predictors at 800 MW.Figure 5
Comparisons among different predictors at 850 MW.Table 4
PRESS and SS values in different methods at 650 MW, 800 MW, and 850 MW1,2.
Loads
Methods
PRESS
NOx
PRESS
Tr
PRESS
SS
NOx
SS
Tr
SS
650 MW
rdCV-MO-NPLS
0.1346
0.0193
0.1539
0.8458
0.1138
0.9596
5-fold-MO-NPLS
0.1354
0.0192
0.1546
0.8458
0.1138
0.9596
SVR
0.1407
0.0294
0.1701
0.8705
0.1347
1.0052
ANN
0.1360
0.0336
0.1696
0.7806
0.1342
0.9148
800 MW
rdCV-MO-NPLS
0.0985
0.0140
0.1125
0.5368
0.1139
0.6507
5-fold-MO-NPLS
0.1397
0.0098
0.1495
0.5368
0.1139
0.6507
SVR
0.1318
0.0295
0.1613
1.3079
0.1580
1.4659
ANN
0.1633
0.0214
0.1847
1.1438
0.2416
1.3854
850 MW
rdCV-MO-NPLS
0.0931
0.0141
0.1072
2.4469
0.2517
2.6986
5-fold-MO-NPLS
0.1068
0.0136
0.1204
2.4469
0.2517
2.6986
SVR
0.1073
0.0165
0.1238
2.5070
0.1665
2.6735
ANN
0.1416
0.0139
0.1555
2.9540
0.2784
3.2324
1The values in bold are the optimal in their columns. 2PRESS=PRESSNOx+PRESSTr; SS=SSNOx+SSTr.It can be seen from Figures3–5 that the prediction of rdCV-MO-NPLS approaches the real values the most while the deviation among 5-fold-MO-NPLS, SVR, or ANN and real value is larger in the three working load conditions. Furthermore, from Table 4, it can be concluded that the PRESS value of rdCV-MO-NPLS are lower than the other three methods in all the three cases, which means rdCV-MO-NPLS has stronger predictive capability. This verified that firstly, rdCV can achieve more robust ability than the k-fold strategy; secondly, SVR can only establish the SISO model, as the boiler combustion system is strongly coupled, so the accuracy of predictive model will be reduced due to decoupling. Finally, as for ANN, it can deal with the MIMO system but is prone to overfitting easily.Remark 3.
Compared with k-fold cross validation, rdCV can achieve more robust ability. This also shows that rdCV can be successfully used in NPLS. Besides, rdCV-MO-NPLS is more powerful than SVR and ANN in generalization ability on the modelling of NOx emission and reheat steam temperature.


## 3.1. MO-NPLS Model with rdCV

Give inputsxi=xi1,…,ximT∈Rm and outputs yi=yi1,yi2T∈R2 with i=1,…,n. Then, let X=x1,…,xnT∈Rn×m and Y=y1,…,yn∈Rn×2 be the number of n observations for inputs and outputs, respectively. Evidently, Xij indicates the i‐th observation of the j‐th input, and this is held on for output Y. The superscript “T” represents transpose. To make the results independent from units, all attributes of data X,Y are normalized into range −1,1 according to the z-score rule and are denoted by E0,F0 with E0∈Rn×m and F0∈Rn×2.Definetj∈Rn and uj∈Rnj=1,…,m are the principal components of X and Y, respectively. Similarly as traditional nonlinear PLS [16, 17], rdCV-MO-NPLS is in nature a recursive modelling strategy, aiming to establish relationship between tj and ujj=1,…,m firstly and then implement the regression relationship of X and Y.Letu1 be any column of F0 at the beginning. Then, one can estimate of the first components (t1, u1) as follows.Firstly,t1 is calculated by(3)t1=E0w1w1Tw1,where w1∈Rm×1 is the weight of matrix E0 and obtained by w1T=u1TE0/u1Tu1.Then, one can fit the polynomial relationship betweent1 and u1 by the least squares method. Specifically, assume the relationship between t1 and u1 is quadratic, namely,(4)u1=ft1+e1=c1,0I+c1,1t1+c1,2t12+e1,where e1∈Rn is the residual vector which is independent and identically distributed. I∈Rn is a vector of ones. Take Z=I,t1,t12∈Rn×3 and c1=c1,0,c1,1,c1,2T∈R3, then (4) can be rewritten by a compact form:(5)u1=Zc1+e1.So, one can get the polynomial coefficientc1 by the least squares method:(6)c1=ZTZ−1ZTu1,where (ZTZ) represents the inverse matrix of ZTZ.After coefficientsc1 are fitted by formula (6), then r1, i.e., the estimate of u1, can be calculated by(7)r1=u^1=ft1=c1,0I+c1,1t1+c1,2t12.Then, updateu1 (mark u1,new) by(8)u1,new=F0q1q1Tq1,where q1∈R2×1 are weights of matrix F0 and determined by q1T=r1TF0/r1Tr1. Next, w1 needs to be updated according to [18].At this point, we can recalculate the vectort1 (mark t1,new) by formula (3) with the updated w1.Check the convergence oft1: if t1,new−t1/t1≥10−6, refit the polynomial coefficient c1 between u1,new and t1,new. Otherwise, calculate the regression coefficient p1 by(9)p1T=t1,newTE0t1,newTt1,new.Finally, the residual ofE0 and F0 (mark E1 and F1) that will be used in the next iteration can be induced by(10)E1=E0−t1,newp1T,F1=F0−r1q1T,with regression coefficient pj.ReplaceE0 and F0 with E1 and F1 and formulas (3)–(10) are repeated until all m principal components (tj,uj) with j=1,…,m are extracted.After allm principal components have been extracted, optimum number afinal (afinal≤m) can be selected. Usually, k-fold strategy is applied, whereas rdCV strategy is adopted here for achieving more robust predictive ability, which has been proven of PLS in [15]. This is the reason why we call our model rdCV-MO-NPLS. Once afinal is determined, the prediction of outputs Y can thus be obtained according to(11)Y^=Y¯+sY,∑k=1afinalrkqkT,k=1,…,afinal,where <⋅,⋅> means point-product operation of two vectors with equal dimensionality; Y¯=Y¯1,Y¯2 with Y¯j=1/n∑i=1nYij,j=1,2; and sY=sY1,sY2 with sYj=1/n∑i=1nYij−Y¯j2.From the above interpretations, the rdCV-MO-NPLS is summarized in Algorithm1.Algorithm 1: rdCV-MO-NPLS.
Input: dataX,Y and parameters Nrep,SEGtest,SEGcalib;
Output: optimal number of principal componentafinal and rdCV-MO-NPLS model.
for rep=1 to Nrep do
Split datasetS equally into SEGtest groups, i.e., S=S1,S2,…,SSEGtest;
for i=1 to SEGtest do
SelectSi as Test set, and construct Calibration sets S−i = S−Si;
Split Calibration set equally intoSEGcalib groups, S−i=V1,V2,…,VSEGcalib;
for j=1 to SEGcalib do
SelectVj as Validation set and V−j=S−i−Vj as Training set;
Fit MO-NPLS based onV−j with, respectively, a=1,2,…,m of principal components,
Apply the MO-NPLS models toVj and get predictive Y^ja for a=1,2,…,m;
Calculate mean square errorMSEja=1/nMSE∑k=1nMSEYjk−Y^jak2, with the number of objects in the used validation set    nMSE and the output objects of validation set Yj;
Estimate optimum principal componentsaopt1i, …, aoptki according to Y^a based on standard error method [15] (Here, more than   one aopt may be selected (i.e., k≥1) because of different confidence interval);
Make MO-NPLS models based on Calibration setS−i with a=aopt1i,…,aoptki;
Test fitted models onSi and obtain a group of k predictions Y^ka as well as k biases;
Find the smallest bias and determine the optimal principal componentaopti,rep.
One can getaopt with the number of SEGtest after completing the outer loop.
Totally, after a complete rdCV run, we can getaopt with the number of (Nrep×SEGtest). The final optimum of principal component afinal is the one with highest frequency in aopti,rep.
Identify rdCV-MO-NPLS model with databaseS and afinal to get model parameters [w, p, c, q] based on (3)–(10).


## 3.2. Model Identification and Validation

### 3.2.1. Model Identification

In this section, rdCV-MO-NPLS is identified, respectively, at 650 MW, 800 MW, and 850 MW. In each case, identification as well as rdCV data are shown in TableA1 (from lines 1 to 45), Table A2 (lines 1 to 30), and Table A3 (lines 1 to 90). In cases of 650 MW and 850 MW, six inputs were considered, including C-F level thermal loads (X1 to X4, MW), running O2 (X5), and αSOFA (X6, %), whereas seven inputs were considered in 800 MW, which are B-F level thermal loads (X1 to X5, MW), running O2 (X6), and αSOFA (X7, %). In all cases, the outputs are(12)Y=NOxX,TrX.Figure2 displays the results of applying rdCV. It can be seen that the highest percentage in three cases are, respectively, 0.32 (aopt=3), 0.38 (aopt=2), and 0.256 (aopt = 2). Hence, the final optimum principal component afinal in each case is 3, 2, and 2, respectively. With the help of afinal, the rdCV-MO-NPLS model can be identified as Case 1 to Case 3.Figure 2
Percentage of principal componentsaopt: (a) 650 MW, (b) 800 MW, and (c) 850 MW.
(a)
(b)
(c)For better comparison, 5-fold cross validation (We call 5-fold-MO-NPLS model) is conducted by the same identification data as the rdCV method. After 5-fold cross validation process, the optimum principal componentafinal of 650 MW, 800 MW, and 850 MW cases is 4, 1, and 4, respectively. Then, the validation results of 5-fold-MO-NPLS model with their afinal will be seen in the next section.Case 1.
(rdCV-MO-NPLS of 650 MW load condition). Sinceafinal = 3, the first three components (t1, u1), (t2, u2), and (t3, u3) are selected. Normalize the training data presented in Table A1 (from lines 1 to 45) into E0=e1,…,e6 and F0 according to z‐score rule. Corresponding parameters of (w,p,c,q) shown in Table 1 as well as (t1, u1), (t2, u2), and (t3, u3) are identified based on (3)–(11). The rdCV-MO-NPLS predictive model is in form of(13)t1=E0w1=E00.1084,0.7642,−0.2669,−0.5667,−0.0943,0.0555T,t2=E1w2=E0I−w1p1Tw2=E0−0.2094,−0.2329,−0.6923,−0.7341,−0.1408,0.1400T,t3=E2w3=E0I−w1p1TI−w2p2Tw3=E00.8325,−0.4500,0.1520,1.2745,0.0707,−0.0796T,f=r1q1T+r2q2T+r3q3T=−0.1112+1.9943t1+0.6504t12+0.626t2+0.3445t22+0.2332t3−0.4783t32,−0.0403+1.9161t1+0.6248t12−0.9299t2−0.5116t22+0.1726t3−0.3540t32≔f1,f2,Y^=Y^1,Y^2=Y¯1,Y¯1+sY1,sY2,f1,f2=509.5+55.8445f1,581.2+6.1878f2.Table 1
Parameters of rdCV-MO-NPLS model at 650 MW.
h
1
w
1
p
1
c
1
q
1
0.1084
−1.0821
−0.2037
0.7211
0.7642
1.7506
2.7657
0.6928
−0.2669
−2.0513
0.9019
−0.5667
1.7544
−0.0943
−1.3794
0.0555
1.7356
h
2
w
2
p
2
c
2
q
2
−0.2521
1.3719
−0.0598
0.5585
−0.5344
−1.1655
1.1211
−0.8295
−0.5870
0.1064
0.6167
−0.5105
−1.2466
0.1779
0.0342
0.1181
1.2109
h
3
w
3
p
3
c
3
q
3
0.7658
1.7521
0.0860
0.8038
−0.0421
−0.7898
0.2901
0.5949
−0.5367
0.0805
−0.5950
0.3318
−0.6007
0.0995
−1.1659
0.0608
−0.2728Case 2.
(rdCV-MO-NPLS of 800 MW load condition). The first two components (t1, u1) and (t2, u2) are chosen. Based on data as shown in Table A2 (lines 1 to 30), and identified parameters of rdCV-MO-NPLS shown in Table 2, the rdCV-MO-NPLS model can be established as(14)t1=E0w1=E0−0.4416,−0.3083,0.3598,−0.1139,−0.1758,−0.5304,−0.5053T,t2=E1w2=E0I−w1p1Tw2=E00.4468,0.0730,−0.3284,−0.0332,0.3795,0.4261,0.6668T,f=r1q1T+r2q2T=−0.1626+1.104t1+0.2081t12+0.2959t2+0.0180t22,−0.1275+0.9698t12+0.1828t1−0.3007t2−0.0183t22≔f1,f2,Y^=Y^1,Y^2=Y¯1,Y¯2+sY1,sY2,f1,f2=451.9+27.3115f1,585.7+4.7006f2.Table 2
Parameters of rdCV-MO-NPLS model at 800 MW.
h
1
w
1
p
1
c
1
q
1
−0.4416
−0.0237
−0.2057
0.7513
−0.3083
−0.2094
1.4695
0.6600
0.3598
0.0880
0.2770
−0.1139
0.2170
−0.1758
−0.0283
−0.5304
−1.0054
0.5053
0.7517
h
2
w
2
p
2
c
2
q
2
0.3187
0.9644
−0.0116
0.7014
−0.0164
0.1064
0.4219
−0.7128
−0.2241
−0.2031
0.0257
−0.0662
−0.4456
0.3286
−0.9193
0.2724
0.3618
0.8133
1.0018Case 3.
(rdCV-MO-NPLS of 850 MW load condition). Similar as Case 2, the first two components (t1, u1) and (t2, u2) were selected on the basis of Table A3 (from lines 1–90). Additionally, identified parameters of rdCV-MO-NPLS are shown in Table 3. The rdCV-MO-NPLS model is(15)t1=E0w1=E00.2595,0.4561,−0.3307,−0.5786,−0.0490,0.5273T,t2=E0I−w1p1Tw2=E00.4214,0.4823,0.2326,−1.2139,0.0331,0.0308T,f=r1q1T+r2q2T=0.0143+1.3261t1+0.0230t12+0.5570t2−0.0635t22,0.0049+0.7791t1−0.0135t12−0.9594t2+0.1095t22≔f1,f2,Y^=Y^1,Y^2=Y¯1,Y¯2+sY1,sY2,f1,f2=464.2+35.8762f1,594.8+2.9182f2.Table 3
Parameters of rdCV-MO-NPLS model at 850 MW.
h
1
w
1
p
1
c
1
q
1
0.2595
0.3142
0.0113
0.8622
0.4561
0.0651
1.5380
0.5066
−0.3307
−0.6749
−0.0266
−0.5786
0.0298
−0.0490
−0.2947
0.05273
1.2676
h
2
w
2
p
2
c
2
q
2
0.1698
−1.6923
0.0092
0.5021
0.0401
−2.2468
1.1094
−0.8648
0.5532
−0.5951
−0.1266
−0.6529
−2.4330
0.0806
−0.2785
−0.4804
−0.2926


### 3.2.2. Model Validation

In general, generalization ability of a prediction model is more important than its fitting ability. Therefore, it is necessary to validate aforementioned models to get their generalization ability. Validation samples in three cases are seen in TableA1 (from lines 46 to 53), Table A2 (from lines 31 to 36), and Table A3 (from lines 91 to 96) in Supplementary file, respectively.To state the superiority of the rdCV-MO-NPLS model, another three well-known methods: 5-fold-MO-NPLS, support vector regression (SVR), and artificial neural network (ANN) are employed to establish models of NOx emission and reheated steam temperature. All methods are on the same identification and validation data. Figures3–5 show the predictive results for different methods. The prediction error sum (PRESS=∑i=1nPRESSY^i−Yi/Yi with prediction sample number nPRESS) representing generalization ability and simulation error sum (SS=∑i=1nSSY^i−Yi/Yi with simulation sample number nPRESS) reflecting approximate ability is presented in Table 4.Figure 3
Comparisons among different predictors at 650 MW.Figure 4
Comparisons among different predictors at 800 MW.Figure 5
Comparisons among different predictors at 850 MW.Table 4
PRESS and SS values in different methods at 650 MW, 800 MW, and 850 MW1,2.
Loads
Methods
PRESS
NOx
PRESS
Tr
PRESS
SS
NOx
SS
Tr
SS
650 MW
rdCV-MO-NPLS
0.1346
0.0193
0.1539
0.8458
0.1138
0.9596
5-fold-MO-NPLS
0.1354
0.0192
0.1546
0.8458
0.1138
0.9596
SVR
0.1407
0.0294
0.1701
0.8705
0.1347
1.0052
ANN
0.1360
0.0336
0.1696
0.7806
0.1342
0.9148
800 MW
rdCV-MO-NPLS
0.0985
0.0140
0.1125
0.5368
0.1139
0.6507
5-fold-MO-NPLS
0.1397
0.0098
0.1495
0.5368
0.1139
0.6507
SVR
0.1318
0.0295
0.1613
1.3079
0.1580
1.4659
ANN
0.1633
0.0214
0.1847
1.1438
0.2416
1.3854
850 MW
rdCV-MO-NPLS
0.0931
0.0141
0.1072
2.4469
0.2517
2.6986
5-fold-MO-NPLS
0.1068
0.0136
0.1204
2.4469
0.2517
2.6986
SVR
0.1073
0.0165
0.1238
2.5070
0.1665
2.6735
ANN
0.1416
0.0139
0.1555
2.9540
0.2784
3.2324
1The values in bold are the optimal in their columns. 2PRESS=PRESSNOx+PRESSTr; SS=SSNOx+SSTr.It can be seen from Figures3–5 that the prediction of rdCV-MO-NPLS approaches the real values the most while the deviation among 5-fold-MO-NPLS, SVR, or ANN and real value is larger in the three working load conditions. Furthermore, from Table 4, it can be concluded that the PRESS value of rdCV-MO-NPLS are lower than the other three methods in all the three cases, which means rdCV-MO-NPLS has stronger predictive capability. This verified that firstly, rdCV can achieve more robust ability than the k-fold strategy; secondly, SVR can only establish the SISO model, as the boiler combustion system is strongly coupled, so the accuracy of predictive model will be reduced due to decoupling. Finally, as for ANN, it can deal with the MIMO system but is prone to overfitting easily.Remark 3.
Compared with k-fold cross validation, rdCV can achieve more robust ability. This also shows that rdCV can be successfully used in NPLS. Besides, rdCV-MO-NPLS is more powerful than SVR and ANN in generalization ability on the modelling of NOx emission and reheat steam temperature.


## 3.2.1. Model Identification

In this section, rdCV-MO-NPLS is identified, respectively, at 650 MW, 800 MW, and 850 MW. In each case, identification as well as rdCV data are shown in TableA1 (from lines 1 to 45), Table A2 (lines 1 to 30), and Table A3 (lines 1 to 90). In cases of 650 MW and 850 MW, six inputs were considered, including C-F level thermal loads (X1 to X4, MW), running O2 (X5), and αSOFA (X6, %), whereas seven inputs were considered in 800 MW, which are B-F level thermal loads (X1 to X5, MW), running O2 (X6), and αSOFA (X7, %). In all cases, the outputs are(12)Y=NOxX,TrX.Figure2 displays the results of applying rdCV. It can be seen that the highest percentage in three cases are, respectively, 0.32 (aopt=3), 0.38 (aopt=2), and 0.256 (aopt = 2). Hence, the final optimum principal component afinal in each case is 3, 2, and 2, respectively. With the help of afinal, the rdCV-MO-NPLS model can be identified as Case 1 to Case 3.Figure 2
Percentage of principal componentsaopt: (a) 650 MW, (b) 800 MW, and (c) 850 MW.
(a)
(b)
(c)For better comparison, 5-fold cross validation (We call 5-fold-MO-NPLS model) is conducted by the same identification data as the rdCV method. After 5-fold cross validation process, the optimum principal componentafinal of 650 MW, 800 MW, and 850 MW cases is 4, 1, and 4, respectively. Then, the validation results of 5-fold-MO-NPLS model with their afinal will be seen in the next section.Case 1.
(rdCV-MO-NPLS of 650 MW load condition). Sinceafinal = 3, the first three components (t1, u1), (t2, u2), and (t3, u3) are selected. Normalize the training data presented in Table A1 (from lines 1 to 45) into E0=e1,…,e6 and F0 according to z‐score rule. Corresponding parameters of (w,p,c,q) shown in Table 1 as well as (t1, u1), (t2, u2), and (t3, u3) are identified based on (3)–(11). The rdCV-MO-NPLS predictive model is in form of(13)t1=E0w1=E00.1084,0.7642,−0.2669,−0.5667,−0.0943,0.0555T,t2=E1w2=E0I−w1p1Tw2=E0−0.2094,−0.2329,−0.6923,−0.7341,−0.1408,0.1400T,t3=E2w3=E0I−w1p1TI−w2p2Tw3=E00.8325,−0.4500,0.1520,1.2745,0.0707,−0.0796T,f=r1q1T+r2q2T+r3q3T=−0.1112+1.9943t1+0.6504t12+0.626t2+0.3445t22+0.2332t3−0.4783t32,−0.0403+1.9161t1+0.6248t12−0.9299t2−0.5116t22+0.1726t3−0.3540t32≔f1,f2,Y^=Y^1,Y^2=Y¯1,Y¯1+sY1,sY2,f1,f2=509.5+55.8445f1,581.2+6.1878f2.Table 1
Parameters of rdCV-MO-NPLS model at 650 MW.
h
1
w
1
p
1
c
1
q
1
0.1084
−1.0821
−0.2037
0.7211
0.7642
1.7506
2.7657
0.6928
−0.2669
−2.0513
0.9019
−0.5667
1.7544
−0.0943
−1.3794
0.0555
1.7356
h
2
w
2
p
2
c
2
q
2
−0.2521
1.3719
−0.0598
0.5585
−0.5344
−1.1655
1.1211
−0.8295
−0.5870
0.1064
0.6167
−0.5105
−1.2466
0.1779
0.0342
0.1181
1.2109
h
3
w
3
p
3
c
3
q
3
0.7658
1.7521
0.0860
0.8038
−0.0421
−0.7898
0.2901
0.5949
−0.5367
0.0805
−0.5950
0.3318
−0.6007
0.0995
−1.1659
0.0608
−0.2728Case 2.
(rdCV-MO-NPLS of 800 MW load condition). The first two components (t1, u1) and (t2, u2) are chosen. Based on data as shown in Table A2 (lines 1 to 30), and identified parameters of rdCV-MO-NPLS shown in Table 2, the rdCV-MO-NPLS model can be established as(14)t1=E0w1=E0−0.4416,−0.3083,0.3598,−0.1139,−0.1758,−0.5304,−0.5053T,t2=E1w2=E0I−w1p1Tw2=E00.4468,0.0730,−0.3284,−0.0332,0.3795,0.4261,0.6668T,f=r1q1T+r2q2T=−0.1626+1.104t1+0.2081t12+0.2959t2+0.0180t22,−0.1275+0.9698t12+0.1828t1−0.3007t2−0.0183t22≔f1,f2,Y^=Y^1,Y^2=Y¯1,Y¯2+sY1,sY2,f1,f2=451.9+27.3115f1,585.7+4.7006f2.Table 2
Parameters of rdCV-MO-NPLS model at 800 MW.
h
1
w
1
p
1
c
1
q
1
−0.4416
−0.0237
−0.2057
0.7513
−0.3083
−0.2094
1.4695
0.6600
0.3598
0.0880
0.2770
−0.1139
0.2170
−0.1758
−0.0283
−0.5304
−1.0054
0.5053
0.7517
h
2
w
2
p
2
c
2
q
2
0.3187
0.9644
−0.0116
0.7014
−0.0164
0.1064
0.4219
−0.7128
−0.2241
−0.2031
0.0257
−0.0662
−0.4456
0.3286
−0.9193
0.2724
0.3618
0.8133
1.0018Case 3.
(rdCV-MO-NPLS of 850 MW load condition). Similar as Case 2, the first two components (t1, u1) and (t2, u2) were selected on the basis of Table A3 (from lines 1–90). Additionally, identified parameters of rdCV-MO-NPLS are shown in Table 3. The rdCV-MO-NPLS model is(15)t1=E0w1=E00.2595,0.4561,−0.3307,−0.5786,−0.0490,0.5273T,t2=E0I−w1p1Tw2=E00.4214,0.4823,0.2326,−1.2139,0.0331,0.0308T,f=r1q1T+r2q2T=0.0143+1.3261t1+0.0230t12+0.5570t2−0.0635t22,0.0049+0.7791t1−0.0135t12−0.9594t2+0.1095t22≔f1,f2,Y^=Y^1,Y^2=Y¯1,Y¯2+sY1,sY2,f1,f2=464.2+35.8762f1,594.8+2.9182f2.Table 3
Parameters of rdCV-MO-NPLS model at 850 MW.
h
1
w
1
p
1
c
1
q
1
0.2595
0.3142
0.0113
0.8622
0.4561
0.0651
1.5380
0.5066
−0.3307
−0.6749
−0.0266
−0.5786
0.0298
−0.0490
−0.2947
0.05273
1.2676
h
2
w
2
p
2
c
2
q
2
0.1698
−1.6923
0.0092
0.5021
0.0401
−2.2468
1.1094
−0.8648
0.5532
−0.5951
−0.1266
−0.6529
−2.4330
0.0806
−0.2785
−0.4804
−0.2926


## 3.2.2. Model Validation

In general, generalization ability of a prediction model is more important than its fitting ability. Therefore, it is necessary to validate aforementioned models to get their generalization ability. Validation samples in three cases are seen in TableA1 (from lines 46 to 53), Table A2 (from lines 31 to 36), and Table A3 (from lines 91 to 96) in Supplementary file, respectively.To state the superiority of the rdCV-MO-NPLS model, another three well-known methods: 5-fold-MO-NPLS, support vector regression (SVR), and artificial neural network (ANN) are employed to establish models of NOx emission and reheated steam temperature. All methods are on the same identification and validation data. Figures3–5 show the predictive results for different methods. The prediction error sum (PRESS=∑i=1nPRESSY^i−Yi/Yi with prediction sample number nPRESS) representing generalization ability and simulation error sum (SS=∑i=1nSSY^i−Yi/Yi with simulation sample number nPRESS) reflecting approximate ability is presented in Table 4.Figure 3
Comparisons among different predictors at 650 MW.Figure 4
Comparisons among different predictors at 800 MW.Figure 5
Comparisons among different predictors at 850 MW.Table 4
PRESS and SS values in different methods at 650 MW, 800 MW, and 850 MW1,2.
Loads
Methods
PRESS
NOx
PRESS
Tr
PRESS
SS
NOx
SS
Tr
SS
650 MW
rdCV-MO-NPLS
0.1346
0.0193
0.1539
0.8458
0.1138
0.9596
5-fold-MO-NPLS
0.1354
0.0192
0.1546
0.8458
0.1138
0.9596
SVR
0.1407
0.0294
0.1701
0.8705
0.1347
1.0052
ANN
0.1360
0.0336
0.1696
0.7806
0.1342
0.9148
800 MW
rdCV-MO-NPLS
0.0985
0.0140
0.1125
0.5368
0.1139
0.6507
5-fold-MO-NPLS
0.1397
0.0098
0.1495
0.5368
0.1139
0.6507
SVR
0.1318
0.0295
0.1613
1.3079
0.1580
1.4659
ANN
0.1633
0.0214
0.1847
1.1438
0.2416
1.3854
850 MW
rdCV-MO-NPLS
0.0931
0.0141
0.1072
2.4469
0.2517
2.6986
5-fold-MO-NPLS
0.1068
0.0136
0.1204
2.4469
0.2517
2.6986
SVR
0.1073
0.0165
0.1238
2.5070
0.1665
2.6735
ANN
0.1416
0.0139
0.1555
2.9540
0.2784
3.2324
1The values in bold are the optimal in their columns. 2PRESS=PRESSNOx+PRESSTr; SS=SSNOx+SSTr.It can be seen from Figures3–5 that the prediction of rdCV-MO-NPLS approaches the real values the most while the deviation among 5-fold-MO-NPLS, SVR, or ANN and real value is larger in the three working load conditions. Furthermore, from Table 4, it can be concluded that the PRESS value of rdCV-MO-NPLS are lower than the other three methods in all the three cases, which means rdCV-MO-NPLS has stronger predictive capability. This verified that firstly, rdCV can achieve more robust ability than the k-fold strategy; secondly, SVR can only establish the SISO model, as the boiler combustion system is strongly coupled, so the accuracy of predictive model will be reduced due to decoupling. Finally, as for ANN, it can deal with the MIMO system but is prone to overfitting easily.Remark 3.
Compared with k-fold cross validation, rdCV can achieve more robust ability. This also shows that rdCV can be successfully used in NPLS. Besides, rdCV-MO-NPLS is more powerful than SVR and ANN in generalization ability on the modelling of NOx emission and reheat steam temperature.


## 4. Combustion Optimization

Artificial bee colony (ABC) is a well-known intelligent optimization algorithm proposed in 2005 [19]. It has been applied more and more on engineering optimization in recent years [20–23] because of its simplicity to implement as it uses fewer control parameters [24, 25] and robust ability to get out of a local minimum [26]. In this section, a multiobjective version of ABC (MO-ABC) [25] is used to optimize NOx emission and reheat steam temperature.


### 4.1. Joint Optimization via Multiobjective Artificial Bee Colony Algorithm

The objective of combustion optimization is to minimize NOx emission and maximize reheat steam temperature (Tr) through searching optimal inputs (x) which is m-dimensional, namely, thermal loads (Qi), running O2, and αSOFA based on rdCV-MO-NPLS models established in Section 3. Considering that these two objectives are contradictory, we can only find a set of nondominated solutions. Objective function can be described as follows:(16)minxY=Y1x,Y2x,s.t.Xmin≤x≤Xmax,i=1,…,m,where Y1x is the NOx emission. Y2x is the inverse of Tr by converting the problem of maximizing Tr to minimize its inverse. The values of Y1x and the inverse of Y2x of 650 MW, 800 MW, and 850 MW can be calculated from formulas (13)–(15), respectively. Xmin=Xmin,1,…,Xmin,mT∈Rm and Xmax=Xmax,1,…,Xmax,mT∈Rm are the lower and upper bounds of inputs in the rdCV-MO-NPLS model.To solve problem 4.1, a MO-ABC is applied step by step as follows:Step 1.
(solutions initialization). New variables namedNp and FoodNumber will be defined as the population size and the number of food sources, respectively. Since each employed bee represents a food source and the number of employed bees and onlooker bees is equal, the FoodNumber can be set to half of the population (i.e., Np/2) [25]. Each food source will be initialized through a function randomly generated in input space according to(17)xid=Xmin,d+rand0,1Xmax,d−Xmin,d,where i is the index of the food source with the range of i=1,2,…, FoodNumber; d is the dimension index of each food source with the range of d=1,2,…,m; rand [0,1] is a random number distributed uniformly over the interval [0,1]; and Xmin,d and Xmax,d are the lower and upper bounds of the dth dimension, respectively.
In addition, assign each food source with a trial indextriali, i=1,2,…, FoodNumber; triali is initialized to 0 for each food source.Step 2.
(solutions evolve in employed bee phase). The initial solutions in (17) first evolved through the help of employed bees according to(18)vi,d=xi,d+ϕi,d⋅xi,d−xk,d,where k≠i and k=1,2,…,FoodNumber and d=1,2,…,m. ϕi,d is a random number distributed uniformly between [−1, 1]. If vi dominates xi, i.e., NOxvi≤NOxxi and 1/Trvi≤1/Trxi (denoted as vi≺xi), update xi to vi. If not, the update was unsuccessful, the trial for the old food source, xi, is incremented by one.Step 3.
(solutions evolve in onlooker bee phase). The quality of solutionxi, i.e., fitness of solution xi, is evaluated by(19)fiti=domxiFoodNumber,where domxi is the number of solutions that are dominated by food source xi. The probability for each food source xi advertised by the corresponding employed bee will be calculated in a roulette way according to(20)Pi=fiti∑k=1FoodNumberfitk.
Once a solution is selected, this solution will evolve to a new position according to formula (18) if this new position dominates its previous one. Otherwise, its trial value is increased by one.Step 4.
(solutions evolve in scout bee phase). There is at most one scout bee in the colony. This means if the maximaltrial value in the colony reaches the limit, then the corresponding food source will turn to a scout bee and after doing a random search according to formula (17). At the same time, its trial value will be reset to zero. The limit is determined by(21)limit=Np2⋅m.Step 5.
(solution archive updates). A fixed-size archive (e.g.,FoodNumber) is conserved to hold the best solutions. Specifically, after the tth cycle, we can get the nondominated solutions. Combine them with solutions which are held in archive at t−1th cycle. Then, select the top FoodNumber solutions by nondomination strategy to update archive.
All the above steps are repeated until the maximum iteration,Cmax, is met. The main process of MO-ABC is summarized in Figure 6.Figure 6
Schematic diagram of the MO-ABC algorithm.Remark 4.
The specific values ofXmin and Xmax are shown in Appendix B (Tables B1–B3) in Supplementary part, and it is determined by history running and experimental data as well as security conditions as the usual practice of data-driven methods (for example, see work [4]).


### 4.2. Optimizing Results and Analysis

To implement MO-ABC as shown above, some parameters should be preset includingNp and limit. Here, Np was set to 100 by a trivial way. limit was set to 300 (at 650 MW and 850 MW) and 350 (at 800 MW), respectively, based on formula (21).Furthermore, another popular multiobjective EA (MOEA), NSGA-II, is also conducted for comparison. Values ofNp and [Xmin,Xmax] in NSGA-II are set the same as MO-ABC.In the present work, we use the number of function evaluations (FEs) as a measure criterion to determine Cmax. Set that both algorithms are terminated after 100,000 FEs. Define CmaxNSGA−II and CmaxMO−ABC as the maximum iteration of NSGA-II and MO-ABC algorithms, respectively. Since FEsNSGA−II=Np∗CmaxNSGA−II=100∗CmaxNSGA−II=100,000, then we can get CmaxNSGA−II = 1000. In addition, the upper bound of the fitness evaluations in basic ABC is defined and given as FEs≤SN+SN+1∗CmaxMO−ABC+SN. In each iteration, there are SN employed bees, SN onlooker bees, and at most 1 scout bee [27]. In this paper, substitute “SN=Np/2=50,FEs≤100,000” into above equation, then we can get CmaxMO−ABC=990.The optimization results of MO-ABC and NAGA-II are presented in Figures7–9.Figure 7
Pareto frontier by MOEAs at 650 MW.Figure 8
Pareto frontier by MOEAs at 800 MW.Figure 9
Pareto frontier by MOEAs at 850 MW.To assess the solution sets shown in Figures7–9, a comprehensive quality indicator (QI), i.e., hypervolume (HV) [28], is applied here. The HV indicator represents the volume of the area enclosed by solution set and a reference point, which can combine the quality of convergence, spread, cardinality, and uniformity of solution set and can be defined as follows.Given a solution setA and a reference point r, HV can be calculated as(22)HVA=λ∪a∈Axa≺x≺r,where λ is the Lebesgue measure. Put it simply, the HV value of a set can be seen as the volume of the union of the hypercubes determined by each of its solutions and the reference point (as the left-bottom vertex and the right-top vertex, respectively).Specifically, for a two-objective optimization problem, HV can be calculated by the following steps. Firstly, the nondominated points in setA are sorted by descending order on their first dimensional objective values. And then, HV is calculated by(23)HVA=∑i=1nobj2pi−obj2r⋅obj1pi−obj1pi−1,where objipj represents the value of the jth solution in A on their ith dimensional objective. obj1p0 is set to obj1r. n is the number of solutions in set A.Ther and HV are illustrated in the third and fourth column of Table 5.Table 5
HV indicator and time consumption by MOEAs.
Loads
MOEAs
r
HV
Time consuming (s)
650 MW
MO-ABC
(1200, 0.002)
0.3574
36.503
NSGA-II
(1200, 0.002)
0.2892
37.021
800 MW
MO-ABC
(500, 0.002)
0.1481
38.105
NSGA-II
(500, 0.002)
0.1062
37.929
850 MW
MO-ABC
(500, 0.002)
0.1239
30.856
NSGA-II
(500, 0.002)
0.0590
29.499Table5 shows that the HV indicators of MO-ABC are higher than NSGA-II in all three load conditions. That is, the solution sets of MO-ABC all outperform NSGA-II (illustrated in Figures 7–9) in the comprehensive performance with comparable time consumption (shown in the last column of Table 5). This means MO-ABC is more applicable to boiler combustion optimization than NSGA-II.Furthermore, we can present our findings in our current study whether the nondominated solutions by MO-ABC shown in Figures7–9 are better than real values (i.e., experimental values). Among them, six recommended nondominated solutions with “+” (the principle of selecting these recommended nondominated solutions is to keep NOx as low as possible while maintaining the reheat steam temperature above its designed value) are shown in Tables 6–8 (No. 1–6) as well as the experiment values of NOx emission and reheat steam temperature (No. 7).Table 6
Recommended setpoints for inputs at 650 MW.
No
.
NOx
Tr
Thermal load (C)
Thermal load (D)
Thermal load (E)
Thermal load (F)
α
SOFA
RunningO2
P
NOx
P
Tr
1
416.2
578.6
146.6
167.5
154.7
205.8
0.225
4.805
6.26
0.28
2
436.6
580.4
143.8
171.0
150.0
209.9
0.209
4.742
1.67
0.59
3
403.2
577.1
146.6
167.6
155.9
205.7
0.229
4.672
9.19
0.02
4
427.8
579.2
146.6
171.3
147.1
209.4
0.190
4.782
3.65
0.38
5
441.3
581.7
140.1
169.5
147.1
210.0
0.203
4.791
0.61
0.81
6
439.5
582.2
142.1
170.4
147.1
209.2
0.211
4.423
1.01
0.90
7
444
577Table 7
Recommended setpoints for inputs at 800 MW.
No
.
NOx
Tr
Thermal load (B)
Thermal load (C)
Thermal load (D)
Thermal load (E)
Thermal load (F)
α
SOFA
RunningO2
P
NOx
P
Tr
1
50.2
588.7
123.1
156.0
166.6
172.9
203.4
0.253
2.211
87.9
0.46
2
57.3
589.6
123.1
155.9
166.6
173.2
203.4
0.253
2.237
86.2
0.61
3
217.6
598.1
123.1
155.4
166.5
173.0
203.4
0.210
2.237
47.8
2.06
4
250.7
599.2
122.4
155.9
166.6
173.5
203.4
0.208
2.236
39.9
2.25
5
440.3
603.1
123.1
157.3
166.2
173.1
203.4
0.196
2.224
56.1
2.92
6
479.4
603.2
123.1
160.0
166.2
173.1
203.4
0.196
2.236
15.0
2.94
7
416.9
586Table 8
Recommended setpoints for inputs at 850 MW.
No
.
NOx
Tr
Thermal load (C)
Thermal load (D)
Thermal load (E)
Thermal load (F)
α
SOFA
RunningO2
P
NOx
P
Tr
1
138.0
592.1
151.1
164.4
206.4
251.3
0.243
2.171
67.1
0.87
2
148.0
593.4
152.6
164.4
205.7
251.3
0.197
2.171
64.7
1.09
3
181.9
599.8
151.3
164.4
186.3
251.3
0.184
2.229
56.6
2.18
4
224.9
602.1
153.1
164.4
186.3
251.3
0.184
2.948
46.4
2.57
5
362.3
606.7
151.1
164.4
186.3
251.3
0.201
4.479
13.6
3.36
6
365.1
607.1
151.1
164.4
186.3
251.0
0.185
4.478
13.0
3.42
7
419.5
587In order to quantify the degree of reduction of NOx and increase of Tr relative to the experimental value, percentage improvement (PNOx and PTr) are defined as follows:(24)PNOx=NOxM−NOxeNOxe×100%,PTr=TrM−TreTre×100%,where NOxM and NOxe are the NOx values obtained by MO-ABC and experiment, respectively. TrM and Tre are Tr values obtained by MO-ABC and experiment, respectively. Results are shown in the last two columns of Tables 6–8.By comparing results in Figures7–9 and Tables 6–8, we can conclude that optimal NOx emission and reheat steam temperature achieved by MO-ABC is better than that explored directly from experimental ones (i.e., real data) based on the percentage improvement. This verified that our new viewpoint on combustion optimization is meaningful. To this end, joint optimization of NOx emission and reheat steam temperature can bring much more benefits to the combustion system with the points of synthetic optimization view.


## 4.1. Joint Optimization via Multiobjective Artificial Bee Colony Algorithm

The objective of combustion optimization is to minimize NOx emission and maximize reheat steam temperature (Tr) through searching optimal inputs (x) which is m-dimensional, namely, thermal loads (Qi), running O2, and αSOFA based on rdCV-MO-NPLS models established in Section 3. Considering that these two objectives are contradictory, we can only find a set of nondominated solutions. Objective function can be described as follows:(16)minxY=Y1x,Y2x,s.t.Xmin≤x≤Xmax,i=1,…,m,where Y1x is the NOx emission. Y2x is the inverse of Tr by converting the problem of maximizing Tr to minimize its inverse. The values of Y1x and the inverse of Y2x of 650 MW, 800 MW, and 850 MW can be calculated from formulas (13)–(15), respectively. Xmin=Xmin,1,…,Xmin,mT∈Rm and Xmax=Xmax,1,…,Xmax,mT∈Rm are the lower and upper bounds of inputs in the rdCV-MO-NPLS model.To solve problem 4.1, a MO-ABC is applied step by step as follows:Step 1.
(solutions initialization). New variables namedNp and FoodNumber will be defined as the population size and the number of food sources, respectively. Since each employed bee represents a food source and the number of employed bees and onlooker bees is equal, the FoodNumber can be set to half of the population (i.e., Np/2) [25]. Each food source will be initialized through a function randomly generated in input space according to(17)xid=Xmin,d+rand0,1Xmax,d−Xmin,d,where i is the index of the food source with the range of i=1,2,…, FoodNumber; d is the dimension index of each food source with the range of d=1,2,…,m; rand [0,1] is a random number distributed uniformly over the interval [0,1]; and Xmin,d and Xmax,d are the lower and upper bounds of the dth dimension, respectively.
In addition, assign each food source with a trial indextriali, i=1,2,…, FoodNumber; triali is initialized to 0 for each food source.Step 2.
(solutions evolve in employed bee phase). The initial solutions in (17) first evolved through the help of employed bees according to(18)vi,d=xi,d+ϕi,d⋅xi,d−xk,d,where k≠i and k=1,2,…,FoodNumber and d=1,2,…,m. ϕi,d is a random number distributed uniformly between [−1, 1]. If vi dominates xi, i.e., NOxvi≤NOxxi and 1/Trvi≤1/Trxi (denoted as vi≺xi), update xi to vi. If not, the update was unsuccessful, the trial for the old food source, xi, is incremented by one.Step 3.
(solutions evolve in onlooker bee phase). The quality of solutionxi, i.e., fitness of solution xi, is evaluated by(19)fiti=domxiFoodNumber,where domxi is the number of solutions that are dominated by food source xi. The probability for each food source xi advertised by the corresponding employed bee will be calculated in a roulette way according to(20)Pi=fiti∑k=1FoodNumberfitk.
Once a solution is selected, this solution will evolve to a new position according to formula (18) if this new position dominates its previous one. Otherwise, its trial value is increased by one.Step 4.
(solutions evolve in scout bee phase). There is at most one scout bee in the colony. This means if the maximaltrial value in the colony reaches the limit, then the corresponding food source will turn to a scout bee and after doing a random search according to formula (17). At the same time, its trial value will be reset to zero. The limit is determined by(21)limit=Np2⋅m.Step 5.
(solution archive updates). A fixed-size archive (e.g.,FoodNumber) is conserved to hold the best solutions. Specifically, after the tth cycle, we can get the nondominated solutions. Combine them with solutions which are held in archive at t−1th cycle. Then, select the top FoodNumber solutions by nondomination strategy to update archive.
All the above steps are repeated until the maximum iteration,Cmax, is met. The main process of MO-ABC is summarized in Figure 6.Figure 6
Schematic diagram of the MO-ABC algorithm.Remark 4.
The specific values ofXmin and Xmax are shown in Appendix B (Tables B1–B3) in Supplementary part, and it is determined by history running and experimental data as well as security conditions as the usual practice of data-driven methods (for example, see work [4]).


## 4.2. Optimizing Results and Analysis

To implement MO-ABC as shown above, some parameters should be preset includingNp and limit. Here, Np was set to 100 by a trivial way. limit was set to 300 (at 650 MW and 850 MW) and 350 (at 800 MW), respectively, based on formula (21).Furthermore, another popular multiobjective EA (MOEA), NSGA-II, is also conducted for comparison. Values ofNp and [Xmin,Xmax] in NSGA-II are set the same as MO-ABC.In the present work, we use the number of function evaluations (FEs) as a measure criterion to determine Cmax. Set that both algorithms are terminated after 100,000 FEs. Define CmaxNSGA−II and CmaxMO−ABC as the maximum iteration of NSGA-II and MO-ABC algorithms, respectively. Since FEsNSGA−II=Np∗CmaxNSGA−II=100∗CmaxNSGA−II=100,000, then we can get CmaxNSGA−II = 1000. In addition, the upper bound of the fitness evaluations in basic ABC is defined and given as FEs≤SN+SN+1∗CmaxMO−ABC+SN. In each iteration, there are SN employed bees, SN onlooker bees, and at most 1 scout bee [27]. In this paper, substitute “SN=Np/2=50,FEs≤100,000” into above equation, then we can get CmaxMO−ABC=990.The optimization results of MO-ABC and NAGA-II are presented in Figures7–9.Figure 7
Pareto frontier by MOEAs at 650 MW.Figure 8
Pareto frontier by MOEAs at 800 MW.Figure 9
Pareto frontier by MOEAs at 850 MW.To assess the solution sets shown in Figures7–9, a comprehensive quality indicator (QI), i.e., hypervolume (HV) [28], is applied here. The HV indicator represents the volume of the area enclosed by solution set and a reference point, which can combine the quality of convergence, spread, cardinality, and uniformity of solution set and can be defined as follows.Given a solution setA and a reference point r, HV can be calculated as(22)HVA=λ∪a∈Axa≺x≺r,where λ is the Lebesgue measure. Put it simply, the HV value of a set can be seen as the volume of the union of the hypercubes determined by each of its solutions and the reference point (as the left-bottom vertex and the right-top vertex, respectively).Specifically, for a two-objective optimization problem, HV can be calculated by the following steps. Firstly, the nondominated points in setA are sorted by descending order on their first dimensional objective values. And then, HV is calculated by(23)HVA=∑i=1nobj2pi−obj2r⋅obj1pi−obj1pi−1,where objipj represents the value of the jth solution in A on their ith dimensional objective. obj1p0 is set to obj1r. n is the number of solutions in set A.Ther and HV are illustrated in the third and fourth column of Table 5.Table 5
HV indicator and time consumption by MOEAs.
Loads
MOEAs
r
HV
Time consuming (s)
650 MW
MO-ABC
(1200, 0.002)
0.3574
36.503
NSGA-II
(1200, 0.002)
0.2892
37.021
800 MW
MO-ABC
(500, 0.002)
0.1481
38.105
NSGA-II
(500, 0.002)
0.1062
37.929
850 MW
MO-ABC
(500, 0.002)
0.1239
30.856
NSGA-II
(500, 0.002)
0.0590
29.499Table5 shows that the HV indicators of MO-ABC are higher than NSGA-II in all three load conditions. That is, the solution sets of MO-ABC all outperform NSGA-II (illustrated in Figures 7–9) in the comprehensive performance with comparable time consumption (shown in the last column of Table 5). This means MO-ABC is more applicable to boiler combustion optimization than NSGA-II.Furthermore, we can present our findings in our current study whether the nondominated solutions by MO-ABC shown in Figures7–9 are better than real values (i.e., experimental values). Among them, six recommended nondominated solutions with “+” (the principle of selecting these recommended nondominated solutions is to keep NOx as low as possible while maintaining the reheat steam temperature above its designed value) are shown in Tables 6–8 (No. 1–6) as well as the experiment values of NOx emission and reheat steam temperature (No. 7).Table 6
Recommended setpoints for inputs at 650 MW.
No
.
NOx
Tr
Thermal load (C)
Thermal load (D)
Thermal load (E)
Thermal load (F)
α
SOFA
RunningO2
P
NOx
P
Tr
1
416.2
578.6
146.6
167.5
154.7
205.8
0.225
4.805
6.26
0.28
2
436.6
580.4
143.8
171.0
150.0
209.9
0.209
4.742
1.67
0.59
3
403.2
577.1
146.6
167.6
155.9
205.7
0.229
4.672
9.19
0.02
4
427.8
579.2
146.6
171.3
147.1
209.4
0.190
4.782
3.65
0.38
5
441.3
581.7
140.1
169.5
147.1
210.0
0.203
4.791
0.61
0.81
6
439.5
582.2
142.1
170.4
147.1
209.2
0.211
4.423
1.01
0.90
7
444
577Table 7
Recommended setpoints for inputs at 800 MW.
No
.
NOx
Tr
Thermal load (B)
Thermal load (C)
Thermal load (D)
Thermal load (E)
Thermal load (F)
α
SOFA
RunningO2
P
NOx
P
Tr
1
50.2
588.7
123.1
156.0
166.6
172.9
203.4
0.253
2.211
87.9
0.46
2
57.3
589.6
123.1
155.9
166.6
173.2
203.4
0.253
2.237
86.2
0.61
3
217.6
598.1
123.1
155.4
166.5
173.0
203.4
0.210
2.237
47.8
2.06
4
250.7
599.2
122.4
155.9
166.6
173.5
203.4
0.208
2.236
39.9
2.25
5
440.3
603.1
123.1
157.3
166.2
173.1
203.4
0.196
2.224
56.1
2.92
6
479.4
603.2
123.1
160.0
166.2
173.1
203.4
0.196
2.236
15.0
2.94
7
416.9
586Table 8
Recommended setpoints for inputs at 850 MW.
No
.
NOx
Tr
Thermal load (C)
Thermal load (D)
Thermal load (E)
Thermal load (F)
α
SOFA
RunningO2
P
NOx
P
Tr
1
138.0
592.1
151.1
164.4
206.4
251.3
0.243
2.171
67.1
0.87
2
148.0
593.4
152.6
164.4
205.7
251.3
0.197
2.171
64.7
1.09
3
181.9
599.8
151.3
164.4
186.3
251.3
0.184
2.229
56.6
2.18
4
224.9
602.1
153.1
164.4
186.3
251.3
0.184
2.948
46.4
2.57
5
362.3
606.7
151.1
164.4
186.3
251.3
0.201
4.479
13.6
3.36
6
365.1
607.1
151.1
164.4
186.3
251.0
0.185
4.478
13.0
3.42
7
419.5
587In order to quantify the degree of reduction of NOx and increase of Tr relative to the experimental value, percentage improvement (PNOx and PTr) are defined as follows:(24)PNOx=NOxM−NOxeNOxe×100%,PTr=TrM−TreTre×100%,where NOxM and NOxe are the NOx values obtained by MO-ABC and experiment, respectively. TrM and Tre are Tr values obtained by MO-ABC and experiment, respectively. Results are shown in the last two columns of Tables 6–8.By comparing results in Figures7–9 and Tables 6–8, we can conclude that optimal NOx emission and reheat steam temperature achieved by MO-ABC is better than that explored directly from experimental ones (i.e., real data) based on the percentage improvement. This verified that our new viewpoint on combustion optimization is meaningful. To this end, joint optimization of NOx emission and reheat steam temperature can bring much more benefits to the combustion system with the points of synthetic optimization view.


## 5. Conclusions

This paper puts a new point on optimizing NOx emission and reheat steam temperature simultaneously using a joint optimization method, in which a set of rdCV-MO-NPLS strategies were proposed as predictors in three cases of load such as 650 MW, 800 MW, and 850 MW. Compared with another three well-known methods (k-fold-MO-NPLS, SVR, and ANN), this method can achieve higher predictive accuracies. Then, MO-ABC was applied to search the optimal set point of controllable process variables to reduce NOx emission and improve reheat steam temperature. Results showed that our joint optimization of boiler combustion with MO-ABC provided a set of tradeoff solutions and outperformed that obtained by NSGA-II and experimental data. This implies that our proposed method on boiler combustion optimization can guarantee higher economy as well as safety control of combustion systems.



---

*Source: 6287291-2019-11-11.xml*