Flue Gas Oxygen Content Model Based on Bayesian Optimiza-Tion Main-Compensation Ensemble Algorithm in Municipal Solid Waste Incineration Process

1. Introduction

The total amount of municipal solid waste (MSW) in developing countries such as China has grown rapidly with the advancement of urbanization [1], resulting in the phenomenon of “garbage siege” in many cities [2]. To promote the goal of “carbon peak carbon neutralization” [3], the treatment of MSW is urgent. At present, MSW incineration (MSWI) technology has been widely used in the world because of its advantages of harmlessness, reduction and resource utilization [4,5,6]. However, it also emits pollutants. The use of artificial intelligence (AI) technology for pollution reduction research is one of the current hot spots.

In the MSWI process, flue gas oxygen content is one of the key controlled variables that can characterize the incineration efficiency and pollutant emissions [7]. When flue gas oxygen content is too low, the combustion efficiency of MSW will be greatly reduced. And many toxic and harmful gases such as dioxin (DXN) and sulfur dioxide (SO₂) will be produced [8]. When flue gas oxygen content is too high, it will increase the emission of combustion NOx gas [8]. Therefore, it is necessary to accurately control flue gas oxygen content. How to construct an effective flue gas oxygen content model is the basis for its accurate control, which is of great significance for the smooth, efficient and environmentally friendly operation of MSWI process.

At present, the modeling of flue gas oxygen content for MSWI process mainly includes mechanism modeling and data-driven modeling [9,10]. Due to the complex mechanism reaction and nonlinear dynamic characteristics, it is difficult to establish an accurate mathematical model of flue gas oxygen content [11]. In the context of AI and information construction, a large amount of industrial data containing mechanism knowledge can be obtained, which makes data-driven modeling methods receive more and more attention [12]. Data-driven modeling is mainly divided into two categories: neural network (NN)-based algorithm [13,14] and tree-based algorithm [15]. The former is weaker than the latter in terms of interpretability. For complex industrial process modeling, there are many reports based on NN modeling. Ma et al. [16] proposed a modeling method of flue gas oxygen content in power plants based on stacked target enhanced autoencoder and long short-term memory (LSTM). Qiao et al. [6] proposed an improved LSTM (ILSTM) modeling method. In the references [17,18,19], a Takagi-Sugeno (TS) fuzzy neural network (TSFNN) controlled object modeling method based on multiple-input multiple output (MIMO) was proposed. Since NN models are composed of many neurons, the connection between neurons is different. With the increase of NN depth and width, the number of parameters increases rapidly, resulting in large memory consumption and resource consumption. Compared with NN, tree models have better generalization and interpretability in terms of small sample modeling. In researches of tree-based modeling, Tang et al. [20] proposed a deep forest regression modeling method based on cross-layer full connection for soft measurement of key process parameters such as product quality and environmental pollution indicators in complex industrial processes. Xia et al. [21] proposed the furnace temperature modeling of MSWI process based on linear regression decision tree (LRDT). However, the MSWI process is complex, and the data has certain uncertainty. And small changes in the input data may lead to completely different tree structures, which makes the model difficult to repeat and explain in some cases, and will reduce the modeling accuracy. Therefore, Xia et al. [22] proposed an ensemble T-S fuzzy regression tree (EnTSFRT) algorithm, introduced fuzzy reasoning [23] to deal with uncertainty, and illustrated the effectiveness of the EnTSFRT model through experiments.

In order to further improve the performance of the model, Yang et al. [24] proposed a dynamic model of a hybrid structure of selective catalytic reduction (SCR) denitrification system. Luo et al. [25] proposed an integrated modeling method that combines the mechanism model with the data-driven model to realize the roller kiln temperature modeling. Dong et al. [26] proposed a hybrid model based on mechanism modeling and data-driven modeling for the construction of leaching system models. The research showed that the compensation model of secondary modeling of the model error can effectively extract the effective information in the error [27]. This main-compensation model mechanism can adopt the same model or different models. For example, gradient boosting decision tree (GBDT) used the same decision tree to construct a multi-level error compensation model [28], and reference [29] adopted NN compensated linear master model. Recenlty, long short term memory network (LSTM) has beed used to model pollutant emission concentration such as NOx and CO in MSWI process with satisfactory performance [30,31]. At present, the main-compensation ensemble controlled object model for flue gas oxygen content has not been reported.

In addition, hyperparameter adjustment based on manual approach is complex and time-consuming in the modeling process, and it is easy to fall into local optimum. To solve this problem, Ma et al. [32] proposed the particle swarm optimization algorithm to optimize the multi-task NN model. Zhou et al. [33] proposed Bayesian optimization (BO)-based online robust parameter selection for sequential support vector regression. Deng et al. [34] used the improved sparrow search algorithm to optimize the squeeze casting process parameters intelligently.

In summary, the article proposes the BO-EnTSFRT-LSTM algorithm. The algorithm takes primary air volume, secondary air volume, feeder speed average, drying grate speed average and ammonia injection amount as the manipulated variables (MVs). At the same time, considering that there are many hyperparameters in the main-compensation ensemble model, BO is proposed to optimize the hyperparameters. Firstly, the main model is constructed based on EnTSFRT. Secondly, the error of the EnTSFRT model is taken as the true value to construct the LSTM-based model. Finally, BO is used to optimize the hyperparameters.

The innovation of the article is reflected in: (1) The main modeling based on tree algorithm is proposed to improve the interpretability of the controlled object model. (2) A method of hyperparameters selection of the main-complement ensemble model based on BO algorithm is proposed. It reduces the complexity and time-consuming of manual hyperparameter adjustment approach and enhances the model performance. (3) A main-complement ensemble model based on BO optimization for flue gas oxygen content is constructed for the first time. It provides support for related intelligent control research.

The rest of the organizational structure of this article is arranged as follows. Section Ⅱ describes the MSWI process for flue gas oxygen content , the overall modeling strategy and the implementation of the algorithm. In section Ⅲ, the experimental research is carried out, and the experimental results are analyzed to verify the feasibility and superiority of the proposed method. Section Ⅳ summarizes and proposes future research directions.

2. Materials and Methods

2.1. MSWI Process Description for Key Controlled Variables

This study is based on an MSWI plant of typical grate furnace in Beijing. Its process flow is shown in Figure 1.

Figure 1 shows that this MSWI procss mainly includes six process stages: solid waste fermentation, solid waste combustion, waste heat exchange, flue gas cleaning, flue gas emission, and steam power generation [35,36]. The detailed description is as follows.

(1)Solid waste fermentation stage. The original MSW contains a large amount of water, which is not conducive to combustion. Therefore, it needs to be bio-fermented in a solid waste reservoir. This process can increase the calorific value of the solid portion by about 30% with about 5-7 days. Then, the fermented MSW transfers it to the hopper and pushes into the incinerator by the feeder and enters the solid waste combustion stage.

(2)Solid waste combustion stage. This stage is to convert MSW into high temperature under the coupling interaction of multi-phase such as solid-gas-liquid and multi-field such as heat flow force. It can be divided into three phases: drying, combustion and burnout. The dring phase is realized based on the furnace’s heat radiation phenomenon and the preheated primary air drying function. During combustion phase, MSW is ignited and produces flue gas. After about 1.5 hours of high-temperature combustion, the combustible components are fully burned and the enters the burnout phase. The non-combustible ash is pushed out of the furnace by the burnout grate.

(3)Waste heat exchange stage. Firstly, the high temperature flue gas is initially cooled by the water wall. Secondly, heat energy is transferred to the boiler by radiation and convection using equipment such as superheaters, evaporators and economizers. Then, the water in the boiler is converted into high-pressure superheated steam and enters the steam generation stage. Finally, the outlet flue gas temperature of the boiler is reduced to 200 °C. In this stage, the cooling rate needs to be strictly controlled to prevent the re-generation of pollutants.

(4)Steam power generation stage. Using the high temperature steam generated by the waste heat boiler. Converting mechanical energy into electrical energy to achieve self-sufficiency in plant-level electricity consumption and external supply of surplus electricity. At the same time, it also can realize resource utilization and obtain economic benefits.

(5)Flue gas cleaning stage. Firstly, the denitrification system is used to remove NOx. Secondly, the acid gas is neutralized by semi-dry deacidification process. Then, activated carbon is used to adsorb DXN and heavy metals in flue gas. Finally, the particles, neutralizing reactants and activated carbon adsorbates in the flue gas are removed by the bag filter.

(6)Flue gas emission stage. The flue gas that meets the national emission standards is discharged into the atmosphere through the chimney through the induced draft fan.

The key controlled variables involved in the solid waste combustion stage are furnace temperature, flue gas oxygen content, steam flow and combustion line. This stage includes three processes: drying, combustion and burnout. In the combustion process, a strong oxidation reaction occurs, indicating that the combustible component and oxygen undergo a complete combustion reaction. Obviously, the oxygen content (air supply volume) and MSW feed rate (grate speed) in the combustion process are crucial to the combustion process. In the burnout process, under the action of high temperature and primary air, the coke will oxidize with oxygen. Therefore, maintaining the stability of flue gas oxygen content is of great significance for the stable operation of combustion process.

2.2. Flue Gas Oxygen Content Control Description

In the MSWI process, flue gas oxygen content is a characterization of the excess oxygen air coefficient. It can characterize the combustion state to a certain extent [1]. When flue gas oxygen content is too low and the excess air ratio is too low, a large amount of toxic and harmful gases will be produced in addition to the increase of heat loss and the decrease of combustion efficiency [7]. When flue gas oxygen content is too large and the excess air coefficient is too large, it shows that the amount of air is too large. The excess air can take away a lot of heat and dust, i.e., high flue gas oxygen content, also increase the emission of NOx pollutants [7].

Normally, the flue gas oxygen content at the outlet of the waste heat boiler is taken as one of the key controlled variables to reflect the MSWI process. Only when flue gas oxygen content is controlled within a certain range, it is possible to make the MSW and combustible flue gas in the furnace fully burned. Thus, the emission concentration of secondary pollutants can be ensure to meet the environmental standards. Therefore, it is of great significance to establish an accurate controlled object model of flue gas oxygen content to support the research of related intelligent control algorithms to realize the smooth, efficient and environmentally friendly operation of MSWI process [8].

The MSWI process is complex, and there are many variables that affect the O₂ at the outlet of the waste heat boiler, and some of them are coupled. Therefore, first of all, the variables are preliminarily selected according to the expert experience [37]. The primary air volume, secondary air volume, feeder speed average, drying grate speed average and ammonia injection amount are selected as MVs. In addition, because the mechanism model is difficult to construct, the construction of data-driven model is the current research hotspot. Without affecting the normal operation of the industrial site, collecting process data for modeling at the same time is the first problem to be solved in building the model.

2.3. Data Acquisition Devices Description and Experimental Data Analysis

Usually, in order to ensure the safe and reliable operation of the MSWI process, DCS systems with closed characteristics in industrial sites are strictly integrated with external devices for data collection and algorithm testing. To address this issue, we collected data using an edge verification platform of MSWI power plant equipped with a safety isolation acquisition device as shown in Figure 2.

The device in Figure 1 facilitates one-way data transmission, mitigating interference issues during the data acquisition process. The industrial data for the experiment are obtained from an actual MSWI power plant in Beijing, China, with a collection frequency of 1 second. 16 hours of continuous process data from 8:00 to 24:00 on a certain day are used. The data are processed through moving averages to eliminate outliers, 857 samples were generated and used to build the FT controlled object model. Further, we divide these samples into training set and testing set in a ratio of 3:1.

There are five key MVs (primary air volume, secondary air volume, feeder speed average, drying grate speed average and ammonia injection amount) and one controlled variable (CV) (flue gas oxygen content). Pearson correlation coefficient ( PCC ) is used to evaluate the correlation between the key MVs and the CV, as shown in Figure 3.

Figure 3 shows that there is a strong linear relationship between ammonia injection amount and flux gas oxygen content for the modeling data used in this article. Obviously, the flux gas oxygen content of the actual MSWI process should have a stronger mapping relationship with other MVs. Since PCC can only measure the linear relationship between MVs and CVs, Figure 3 also indicates that there is a complex nonlinear or even uncertain mapping relationship between flux gas oxygen content and MVs. Therefore, it is necessary to adopt a modeling strategy that can characterize strong uncertainty and complex mapping relationships to construct a controlled object model for flux gas oxygen content.

2.4. Modeling Strategy

Based on the previous analysis, it can be concluded that constructing a flue gas oxygen content model for the MSWI process is the foundation for implementing stable control algorithms. In this article, an interpretable flue gas oxygen content based on BO-EnTSFRT-LSTM is proposed. In which, the main model is constructed by using EnTSFRT to model the uncertainty, the LSTM is used to construct the compensation model by using the main model error as the true value, and the hyperparameters of the constructed main-complement ensemble model are optimized based on the BO algorithm. The overall modeling strategy is shown in Figure 4.

In Figure 3, where $u$ denotes MVs, $\mathit{y}$ represents the true value of flue gas oxygen content, ${\widehat{\mathit{y}}}_{\mathrm{main}}$ represents the output value of EnTSFRT main model, ${\widehat{e}}_{\mathrm{main}}$ represents the error value of EnTSFRT main model, ${\widehat{e}}_{\mathrm{comp}}$ denotes the output value of LSTM compensation model, ${\widehat{e}}_{\mathrm{mix}}$ denotes the error value of LSTM compensation model, $\widehat{\mathit{y}}$ represents the output value of BO-EnTSFRT-LSTM. The relationships among above the variables are as follow:

(1)

(2)

In addition, $\mathrm{lr}$ , $\mathrm{bs}$ , $\lambda$ , ${\mathrm{LSM}}_{\mathsf{\alpha }}$ , $\mathrm{MinSamples}1$ , $\mathrm{Tn}$ , $\beta$ , $\upsilon$ and $S$ in Figure 4 represent the hyperparameters of the main-complement ensemble model that the BO algorithm needs to optimize. They are learning rate, number of batches, the intermediate variable in the recursive calculation process, regularization coefficient, the minimum number of samples and the number of decision trees of the main model and the learning rate, regularization coefficient and maximum number of iterations of the compensation model, respectively.

The functions of different modules are as follows:

1) Main model construction module based on EnTSFRT: It is obtained by constructing multiple TSFRT sub-models with parallel fusion mode;

2) Compensation model construction module based on LSTM: It is used to construct the compensation model with the error true value of the EnTSFRT model;

3) Hyperparameter optimization module based on BO algorithm: The hyperparameters of the constructed main-complement ensemble model are optimized based on the BO algorithm.

2.5. Algorithm Implementation

2.5.1. Main Model Construction Module Based on EnTSFRT

The EnTSFRT model requires parallel fusion of multiple TSFRT sub-models. The TSFRT sub-model mainly includes a screening layer and a fuzzy inference layer. Among them, the screening layer is used for feature screening, and the fuzzy inference layer uses TS fuzzy inference.

where $\mathit{D}={\left\{{\mathit{u}}_n,{\mathit{y}}_n\right\}}_{n=1}^N\in {ℝ}^{N\times M}$ denotes the input of EnTSFRT, ${\mathit{u}}_n$ denotes MVs, $n$ denotes the $n\mathrm{th}$ input sample, $N$ represents the total number of data set samples, $M$ represents dimension of the input feature, ${\mathit{U}}^j$ represents the input variable matrix of the $j\mathrm{th}$ subset sample, ${\mathit{D}}^j$ denotes the $j\mathrm{th}$ training subset obtained by random feature sampling of $\mathit{D}$ , ${\mu }_{\mathrm{CS}}^{\mathit{v}+1}(\cdot )$ and ${\mu }_{\mathrm{CS}}^{\mathit{v}+2}(\cdot )$ denote the membership degree on the non-leaf nodes of $\mathit{v}+1$ and $\mathit{v}+2$ , $f_{\mathrm{TSFRT}}$ denotes the TSFRT sub-model.

2.5.1.1. Training Subsets Partition Submodule

The training subset ${\mathit{D}}^j$ of $M$ features is obtained by feature random sampling of $\mathit{D}$ . It can be expressed as:

(3)

where ${\mathit{U}}^j={\left\{{\mathit{u}}_n^j\right\}}_{\mathrm{n}=1}^N$ denotes the input variable matrix, ${\mathit{y}}^j={\left\{{\mathit{y}}_n^j\right\}}_{n=1}^N$ represents the true value vector of flue gas oxygen content.

Correspondingly, from the perspective of input features, ${\mathit{u}}^j$ with $M$ input features can be expressed as:

(4)

Therefore, a subset denoted as $\left\{{\mathit{D}}^1,\cdots ,{\mathit{D}}^j,\cdots {\mathit{D}}^{\mathrm{Tn}}\right\}$ can be obtained.

2.5.1.2. Sub-Model Construction Submodule

The construction of the $jth$ TSFRT sub-model is described as an example.

The data set ${\mathit{D}}^j$ is used as the input of the screening layer. And the best segmentation position is determined. Each recursive coordinate point $\left(j,m\right)$ in ${\mathit{D}}^j$ is obtained. And its mean square error (MSE) value is calculated to minimize MSE to determine the segmentation node. The estimation process of the loss function value is as follows:

(5)

where $\Phi$ represents the loss value of data set ${\mathit{D}}^j$ , $f_{\mathrm{MSE}}^j\left(\cdot \right)$ denotes the MSE value of the left subset ${\mathit{D}}_{\mathrm{L}}^j$ (or right subset ${\mathit{D}}_{\mathrm{R}}^j$ ), ${\mathit{y}}_{\mathrm{L}}^j$ denotes the true vector of the left subset ${\mathit{D}}_{\mathrm{L}}^j$ , ${\mathit{y}}_{\mathrm{R}}^j$ denotes the true vector of the right subset ${\mathit{D}}_{\mathrm{R}}^j$ , ${\mu }_{cs}^{\nu }\left({\mathit{u}}^j\right)$ denotes the clear membership function of the input ${\mathit{u}}^j$ on the $\nu \mathrm{th}$ non-leaf node, ${\vartheta }_{{\nu }_{\mathrm{L}}}^j$ represents the average of the target values of the left subset ${\mathit{D}}_{\mathrm{L}}^j$ , ${\vartheta }_{{\nu }_{\mathrm{R}}}^j$ represents the average of the target values of the right subset ${\mathit{D}}_{\mathrm{R}}^j$ .

In this article, both ${\vartheta }_{{\nu }_{\mathrm{L}}}^j$ and ${\vartheta }_{{\nu }_{\mathrm{R}}}^j$ are represented as ${\vartheta }_t^j$ . Its calculation is as follows:

(6)

where $N_{S\mathit{u}bset}^j$ denotes the number of ${\mathit{D}}_{\mathrm{L}}^j$ (or ${\mathit{D}}_{\mathrm{R}}^j$ ).

In Equation (5), ${\mu }_{cs}^{\mathit{v}}\left({\mathit{u}}^j\right)$ is calculated as:

(7)

where ${\mu }_{cs}\left({\mathit{u}}^j\right)$ denotes the clear membership function of the input ${\mathit{u}}^j$ , ${\varphi }_{\mathit{v}}$ indicates that the segmentation threshold can be expressed as ${\varphi }_{\mathit{v}}={\mathit{u}}_{j,m}$ . It is assumed that the TSFRT sub-model in this article consists of $V$ nodes. So, the number of non-leaf nodes is ${\displaystyle \frac{V}{2}}-1$ , and the sum of membership sets is ${\left\{{\mu }_{\mathrm{CS}}^{\mathit{v}}(\cdot )\right\}}_{\mathit{v}=1}^{{\displaystyle \frac{V}{2}}-1}$ .

${\mathit{D}}^j$ is divided into left and right subsets as follows:

(8)

where ${\mu }_{cs}^1\left(\cdot \right)$ denotes the first membership degree of a clear set $C_{\mathrm{cs}}^{{\mathit{v}}_{\mathrm{leaf}}^j}$ obtained by using the minimum MSE. The first element of the clear set is defined by Equation (7). The clear set is described as:

(9)

Repeat the previous step until the minimum number of samples at the leaf node reaches the set threshold. The $t\mathrm{th}$ clear set can be expressed as:

(10)

It can be further expressed as:

(11)

The training data ${\mathit{D}}^{{\mathit{v}}_{\mathrm{leaf}}^j}$ of the node is determined according to $C_{\mathrm{cs}}^{{\mathit{v}}_{\mathrm{leaf}}^j}$ . And the T-S fuzzy inference is constructed. ${\mathit{D}}^{{\mathit{v}}_{\mathrm{leaf}}^j}$ denotes the following:

(12)

where ${\mathit{D}}^{{\mathit{v}}_{\mathrm{leaf}}^j}$ represents the training data of fuzzy reasoning, ${\mathit{v}}_{\mathrm{leaf}}^j$ denotes the ${\mathit{v}}_{\mathrm{leaf}}^j\mathrm{th}$ leaf node, and $N_{{\mathit{v}}_{\mathrm{leaf}}^j}^j$ denotes the number of samples in the ${\mathit{v}}_{\mathrm{leaf}}^j\mathrm{th}$ leaf node.

In the fuzzy inference layer, $K$ rules are defined to represent the local linear relationship between the input variable and the target. Define the $k\mathrm{th}$ rule as follows:

(13)

where ${\mathit{y}}^k$ represents the output of the consequent of the $k\mathrm{th}$ fuzzy rule, that is, the weight of the $k\mathrm{th}$ fuzzy rule. Further, we simplify Equation (13) as follows:

(14)

where ${\mathit{u}}_1^{{\mathit{v}}_{\mathrm{leaf}}^j}$ denotes a variable of $C_{\mathrm{cs}}^{{\mathit{v}}_{\mathrm{leaf}}^j}$ , $k$ denotes the membership degree of the $k\mathrm{th}$ fuzzy rule, ${\mu }_{A_1^k}^k\left(\cdot \right)$ denotes the membership function of $A_1^k$ , and describes the membership of ${\mathit{u}}_1^{{\mathit{v}}_{\mathrm{leaf}}^j}$ to $A_1^k$ .

In the article, Gaussian function is used as membership function, as follows:

(15)

where ${\mathsf{\delta }}_{m,k}$ denotes the center of the membership function, and ${\sigma }_{m,k}^2$ represents the width of membership function.

The activation intensity of the $k\mathrm{th}$ fuzzy rule of the fuzzy rule of $m$ input variables is calculated as follows:

(16)

where ${\mathsf{\alpha }}^k$ denotes the product output of the $k\mathrm{th}$ fuzzy rule.

Further normalization operation is as follows:

(17)

According to Equations (13)–(17), ${\mathit{y}}^k$ and ${\overline{\mathsf{\alpha }}}^k$ are linearly combined and accumulated. The predicted output values $\widehat{\mathit{y}}$ are obtained as follows:

(18)

where $g\left(\cdot \right)$ represents the consequent output function, and $\mathsf{\omega }$ denotes the consequent weight parameter.

Repeat the above steps until ${\displaystyle \frac{V}{2}}$ leaf nodes are obtained. By obtaining ${\displaystyle \frac{V}{2}}$ fuzzy inferences to construct a multi-input single-output model, the TSFRT simplified model of multi-input single-output can be obtained as follows:

(19)

where $f_{\mathrm{TSFRT}}^j$ denotes the $j\mathrm{th}$ TSFRT sub-model, $\mathit{U}=\left\{{\mathit{U}}^1,\cdots ,{\mathit{U}}^j,\cdots ,{\mathit{U}}^{\mathrm{Tn}}\right\}$ , ${\theta }_{\mathrm{leaf}}$ denotes the minimum number of samples, $\mathsf{\omega }=\left\{{\mathsf{\omega }}_1,{\mathsf{\omega }}_2,\cdots ,{\mathsf{\omega }}_m\right\}$ denotes the consequent weight parameter matrix, $\mathsf{\delta }$ represents the center matrix of membership function, $\mathcal{l}$ represents the width matrix of membership function.

In addition, the article uses the least squares method to update the weight of the consequent. Equation (18) can be merged as follows:

(20)

where ${\mathit{u}}^{*}=[{\overline{\mathsf{\alpha }}}_1{\mathit{u}}_1,{\overline{\mathsf{\alpha }}}_2{\mathit{u}}_2,\cdots ,{\overline{\mathsf{\alpha }}}_m{\mathit{u}}_m]\in {ℝ}^{1\times m}$ .

Given the input and output matrices, the T-S consequent weights can be estimated as follows:

(21)

where ${\mathit{U}}^{*}$ is composed of ${\mathit{u}}^{*}$ .

The least square method is used for parameter learning. The recursive calculation is as follows:

(22)

(23)

where ${\mathsf{\omega }}_0$ is the weight of the initial time, which is a random vector; $l$ denotes the $l\mathrm{th}$ recursive process; $H_l$ is the intermediate variable in the recursive process; $H_0$ is defined by $H_0=\gamma I$ , $\lambda$ denotes a large positive value, and $I$ represents the unit matrix.

2.5.1.3. Parallel Fusion Submodule

Firstly, according to the modeling process of TSFRT, total TSFRT models ${\left\{f_{\mathrm{TSFRT}}^j\right\}}_{j=1}^{\mathrm{Tn}}$ are constructed, which are defined as $F\in {ℝ}^{N\text{×}\mathrm{Tn}}$ .

Then, to estimate the weight with the minimum training error, the following optimal problem is used to calculate the pseudo-inverse:

(24)

where ${\mathrm{LSM}}_{\mathsf{\alpha }}$ denotes the regularization coefficient.

Using the ridge regression theory, that is, the Moore-Penrose inverse matrix, the weight matrix is calculated. The calculation results are as follows:

(25)

Finally, the output of EnTSFRT is as follows:

(26)

2.5.1.4. Flow Chart of Main Model Modeling

The flow chart of the construction of the main model is as follows:

Figure 6. The flow chart of the construction of the main model.

Figure 6. The flow chart of the construction of the main model.

2.5.2. Compensation Model Construction Module Based on LSTM

LSTM is a special recurrent neural network (RNN). The difference is that LSTM introduces three gates of forget gate, input gate and output gate and one cell state. The structure of LSTM is shown as follows :

Figure 7. Structure of compensation model construction module based on LSTM.

Figure 7. Structure of compensation model construction module based on LSTM.

where $t$ denotes the current moment, $h$ denotes the hidden state, $c$ denotes the cell state, $\tilde{c}$ denotes the candidate cell state, $\zeta$ represents the sigmoid function, $\psi$ represents the tanh function, $f$ denotes the value of the forget gate, $i$ denotes the value of the input gate, $o$ represents the value of the output gate, $\otimes$ represents element-by-element multiplication.

Its forward calculation and back propagation are as follows.

2.5.2.1. Forward Calculation Process

The output of forget gate, input gate, output gate, candidate cell state, cell state and hidden state can give the network memory ability of repetitive unit [7]. The forward calculation process at time $t$ are as follows.

Firstly, update the output of the forget gate as follows:

(27)

where $\mathit{W}$ represents the weight matrix, $b$ represents the bias matrix.

Secondly, the output of the input gate is updated and new candidate cell state is generated, as follows:

(28)

(29)

Further, the cell state is updated as follows:

(30)

Then, the output gate is updated and the cell state at the current moment is generated, as follows:

(31)

(32)

Finally, the output value ${\widehat{e}}_{\mathrm{comp}}\left(t\right)$ of LSTM is obtained as follows :

(33)

where ${\mathit{W}}_{o\mathit{u}t}$ is the output weight matrix.

2.5.2.2. Back Propagation Process

The loss function $J\left(t\right)$ in the training process of LSTM is defined as follows:

(34)

The weights and biases are updated based on the back propagation algorithm. Taking the update of ${\mathit{W}}_{\mathit{u}f}$ as an example, as follows:

(35)

Using the chain rule, ${\displaystyle \frac{\partial J\left(t\right)}{\partial {\mathit{W}}_{\mathit{u}f}}}$ is calculated as follows:

(36)

The calculation of ${\displaystyle \frac{\partial {\widehat{e}}_{\mathrm{comp}}\left(t\right)}{\partial {\mathit{W}}_{\mathit{u}f}\left(t\right)}}$ is as follows :

(37)

Further, ${\mathit{W}}_{\mathit{u}f}\left(t+1\right)$ can be obtained.

The update process of the remaining weights and biases is the same as the update process of ${\mathit{W}}_{\mathit{u}f}$ .

2.5.2.3. Flow Chart of Compensation Model Modeling

The flow chart of the construction of the compensation model is as follows:

Figure 8. The flow chart of the construction of the compensation model.

Figure 8. The flow chart of the construction of the compensation model.

where ${\mathit{D}}^E$ is the input of the compensation model, denoted as ${\mathit{D}}^E={\left\{\mathit{u}\left(t\right),{\widehat{e}}^{\mathrm{main}}\left(t\right)\right\}}_{t=1}^N$ ; and $s$ represents the $s\mathrm{th}$ iteration number.

2.5.3. Hyperparameter Optimization Module Based on BO Algorithm

2.5.3.1. The Principle of BO Algorithm

BO is a kind of optimization method based on machine learning. Its purpose is to minimize the objective function. In the article, the eight hyperparameters in the main complement ensemble model are selected based on the BO algorithm to achieve the purpose of minimizing the objective function. The objective function selected in this article is the negative vaue of R-square (R²) indicator.

The optimization process of BO algorithm is as follows.

Step 1: Set the boundary space

Set the value range of the hyperparameters $X=\left\{{\mathit{x}}_1;\cdots ;{\mathit{x}}_9\right\}$ of BO. It includes the upper boundary (UB) and the lower boundary (LB) of the sample space, which is set as follows:

(38)

where $\left\{{\mathit{x}}_1;\cdots ;{\mathit{x}}_8\right\}$ denotes as follows:

(39)

Step 2: Construct the initial data set

Randomly initialize $P_{\mathrm{init}}$ input samples $X=\left\{{\mathit{x}}_1;\cdots ;{\mathit{x}}_{P_{\mathrm{init}}}\right\}$ . And observe the corresponding objective function value $\mathsf{\Gamma }={\left[\mathsf{\Gamma }\left({\mathit{x}}_1\right),\cdots ,\mathsf{\Gamma }\left({\mathit{x}}_{P_{\mathrm{init}}}\right)\right]}^{\mathrm{T}}$ of the sample to form the initial data set.

Step 3: Construct / train the GPR model

The Gaussian process regression (GPR) model is constructed by using the data set. The GPR model uses the known data as a prior probability distribution that satisfies the multivariate normal distribution. It is described as follows:

(40)

where ${\mathsf{\xi }}_0=[{\mathsf{\xi }}_0({\mathit{x}}_1),\cdots ,{\zeta }_0({\mathit{x}}_p)]$ is the mean function, $p$ represents the number of data set samples at the current time, and ${\mathsf{\Sigma }}_0$ denotes the covariance matrix as follows:

(41)

where ${\sigma }^2$ represents the variance of the measurement noise, ${\mathsf{\Sigma }}_0({\mathit{x}}_a,{\mathit{x}}_b)$ is the kernel function between the $a\mathrm{th}$ and $b\mathrm{th}$ sample points. It is calculated as follows:

(42)

where ${\mathit{x}}_{ao}$ represents the $o\mathrm{th}$ dimension feature of the $a\mathrm{th}$ sample, ${\sigma }_{\mathsf{\Gamma }}$ is the kernel amplitude of the kernel function, ${\eta }_o$ denotes the length ratio of the $o\mathrm{th}$ dimension of the sample.

For a new sample point in the sample space, the posterior probability distribution of the sample point can be calculated by the prior probability distribution. According to the definition of GP, the joint Gaussian distribution $\mathsf{\Gamma }$ and the predicted value ${\widehat{\mathsf{\Gamma }}}^{*}$ of the known data also obey the Gaussian distribution, which can be expressed as:

(43)

where ${\mathsf{\Sigma }}_0^{\ast }={\left[{\mathsf{\Sigma }}_0({\mathit{x}}_1,{\mathit{x}}_1),\cdots ,{\mathsf{\Sigma }}_0({\mathit{x}}_p,{\mathit{x}}^{\ast })\right]}^{\mathrm{T}}\in {ℝ}^{p\times 1}$ represents the covariance of training testing set, and ${\mathsf{\Sigma }}_0^{\ast \ast }={\mathsf{\Sigma }}_0({\mathit{x}}^{\ast },{\mathit{x}}^{\ast })\in {ℝ}^{1\times 1}$ denotes the covariance of the testing set.

Since the training set is known, the conditional distribution for $\mathsf{\Gamma }$ given ${\widehat{\mathsf{\Gamma }}}^{*}$ can be calculated based on Equation (43) as follows:

(44)

where ${\mathsf{\xi }}_p$ and ${\sigma }_p^2$ are as follows:

(45)

where ${\mathsf{\xi }}_p$ is the posterior mean value, and ${\sigma }_p^2$ denotes the posterior variance.

Step 4: Estimate the parameters of the GPR model

We estimate the parameters in the mean function and covariance matrix. Maximum likelihood estimation (MLE) is used to estimate the above parameters, which improves the fitting and prediction performance of the model [38,39].

Step 5: Sampling

In order to find the optimal point ${\mathit{x}}_{\mathrm{op}}$ to minimize the objective function value $\mathsf{\Gamma }\left({\mathit{x}}_{\mathrm{op}}\right)$ , the expected improvement (EI) sampling function is used to sample the objective function.

Step 6: Update the data set

We calculate the objective function value $\mathsf{\Gamma }\left({\mathit{x}}^{\ast }\right)$ corresponding to the new sample point ${\mathit{x}}^{\ast }$ . Then, we update the data sets $X=\left\{{\mathit{x}}_1;\cdots ;{\mathit{x}}_p;{\mathit{x}}^{\ast }\right\}$ and $\mathsf{\Gamma }={\left[\mathsf{\Gamma }\left({\mathit{x}}_1\right),\cdots ,\mathsf{\Gamma }\left({\mathit{x}}_p\right),\mathsf{\Gamma }\left({\mathit{x}}^{\ast }\right)\right]}^{\mathrm{T}}$ .

Repeat steps 3-6 until the preset maximum number of iterations is reached. And find the optimal sample point ${\mathit{x}}_{\mathrm{op}}$ that meet the requirements from the dataset.

2.5.3.2. The Flow Chart of BO Algorithm

The flow chart of BO algorithm is as follows:

Figure 9. The flow chart of BO algorithm.

Figure 9. The flow chart of BO algorithm.

2.6. pseudo code

The pseudo code of the pseudocode algorithm used in this article is as follows:

where $\mathrm{Loop}$ represents the number of optimization iterations, $\mathsf{\Omega }$ represents the sampling function, $\mathrm{ratio}$ represents the exploration ratio, $\varsigma$ is the current minimum number of samples.

Table 1 shows that the parameters of BO process are $\mathrm{Loop}$ , $\mathsf{\Omega }$ , and $\mathrm{ratio}$ . Among them, BO establishes a probability model based on the past evaluation results of the objective function and finds the value that minimizes the objective function. The commonly used sampling functions include: expectation boosting, knowledge gradient, and entropy search and predictive entropy search. Its selection needs to consider how likely it is to obtain the optimal value at a certain point, and to evaluate whether the point can reduce the uncertainty of the Bayesian statistical model. The exploration ratio is a parameter that controls the balance between exploration and exploitation, and it needs to be adjusted according to specific optimization tasks and application scenarios.

3. Results and Discussion

3.1. Performance Metrics

The evaluation indexes used in the article are root mean squard error (RMSE), mean absolute error (MAE) and R-Square (R²), as follows:

(46)

(47)

(48)

where RMSE represents the standard deviation of the sample for the difference between the predicted and true values, reflecting the degree of dispersion of the sample. MAE represents the average absolute error between the predicted and true values. And when doing nonlinear fitting, the smaller the RMSE, the better[40]. R² is used to quantify the fitting degree of the model to the data. Generally speaking, the closer the value of R² is to 1, the better the fitting degree of the model is.

3.2. Experimental Results

The parameters of BO are set based on Thumb rule. The detaied values are follows: $\mathrm{Loop}=20$ , $\mathsf{\Omega }$ = EI function, and $\mathrm{ratio}$ = 0.5. The optimal hyperparameters are selected as $\mathrm{lr}=0.1829$ , $\lambda =23$ , $\mathrm{bs}=21$ , ${\mathrm{LSM}}_{\mathsf{\alpha }}=0.6714$ , $\mathrm{MinSamples}1=42$ , $\mathrm{Tn}=23$ , $\beta =0.0263$ , $\upsilon =0.1923$ and $S=345$ . The fitting curves of the training set and testing set of the main-compensation ensemble model are shown in Figure 10 and Figure 11.

According to Figure 10 and Figure 11, only using the EnTSFRT-based main model, the predicted value cannot fit the actual value effectively. The compensation model has a better fitting effect on the error value of the main model. At the same time, by comparing Figure 10(a)-Figure 11(a) with Figure 10(c)-Figure 11(c), it can be seen that the fitting effect of the proposed method on the training set and the testing set has a great improvement.

The weights of all trees constructed by the main model are shown in Figure 12:

The curve of the objective function in the BO optimization process is as follow:

Figure 13 shows that the value of the objective function is in the range of -0.5079 to -0.7988 in the BO optimization process.

The variation curves of the 9 hyperparameters in the BO optimization process are as follows:

Figure 14 shows the change process of hyperparameters of the main-complement ensemble model in the BO optimization process. The specific detailed analysis is as follow:

1) Figure 14(a) shows the change curves of $lr$ , which is the learning rate of the main model. When the range of $Loop$ is 0 to 50, the change of $lr$ is dramatic. When the range of $Loop$ is 50 to 280, the value of $lr$ is basically stable at 0.18. When the $Loop$ is greater than 280, it first stabilizes at about 0.3. Then changes sharply and finally tends to be stable at about 0.18.

2) Figure 14(b) shows the change curves of $\lambda$ , which is the intermediate variable in the recursive calculation process of the main model. With the increase of $Loop$ , its change trend is similar to that of $lr$ . Its value is finally stabilized at about 16.

3) Figure 14(c) shows the change curves of $bs$ , which is the the number of batches of the main model training. With the increase of $\mathrm{Loop}$ , its change trend is similar to that of $lr$ . Its value is finally stabilized at about 25.

4) Figure 14(d) shows the change curves of $LS{M}_a$ , which is the regularization coefficient of the master model. With the increase of $Loop$ , its change trend is similar to that of $lr$ . Its value is finally stabilized at about 0.48.

5) Figure 14(e) shows the change curves of $MinSamples1$ , which is the minimum number of leaf node samples of the main model. With the increase of $\mathrm{Loop}$ , the value of $MinSamples1$ first changed dramatically. Then its value is basically stable at about 41.

6) Figure 14(f) shows the change curves of $Tn$ , which is the number of decision trees of the main model. When $\mathrm{Loop}$ is 0 to 20, $Tn$ changes drastically. When $\mathrm{Loop}$ is greater than 20, the value of $Tn$ is stable at about 23.

7) Figure 14(g) shows the change curves of $\beta$ , which is the learning rate of the compensation model. With the increase of $\mathrm{Loop}$ , the change trend of its value is similar to that of $Tn$ . Its value eventually converges to about 0.026.

8) Figure 14(h) shows the change curves of $\upsilon$ , which is the regularization coefficient of the compensation model. When $\mathrm{Loop}$ is in the range of 0 to 40, its value fluctuates greatly. Then, its value converges to about 0.2. When $\mathrm{Loop}$ is greater than 280, the value is mostly around 0.01.

9) Figure 14(i) shows the change curves of $S$ , which is the iteration times of the compensation model. When $Loop$ is in the range of 0 to 50, its value changes dramatically. When $Loop$ is greater than 50, its values converge to 328, 345 and 382 respectively.

3.3. Comparison and Discussion

To further illustrate the effectiveness, the proposed method is compared with tree-based method such as random forest (RF), LRDT, EnTSFRT-LRDT and BO-EnTSFRT. The hypeparameters of each model are set as follows. In RF model, the minimum number of samples is set to 50, and the number of features is 4. In LRDT model, the minimum number of samples is set to 20, and the regularization coefficient is 0.4. In EnTSFRT-LRDT model, $lr=0.1$ , $\lambda =10$ , $bs=18$ , $LS{M}_a=0.3$ , $MinSamples1=30$ , $Tn=160$ , $MinSamples2=35$ , and $\beta =0.3$ . In EnTSFRT-LSTM model, $lr=0.1$ , $\lambda =15$ , $bs=18$ , $LS{M}_a=0.3$ , $MinSamples1=55$ , $Tn=80$ , $\beta =0.05$ , $\upsilon =0.3$ , $S=320$ . In BO-EnTSFRT model, $lr=0.8554$ , $\lambda =18$ , $bs=41$ , $LS{M}_a=0.919$ , $MinSamples1=45$ , and $Tn=200$ . In BO-EnTSFRT-LRDT model, $lr=0.024$ , $\lambda =6$ , $bs=48$ , $LS{M}_a=0.6123$ , $MinSamples1=77$ , $Tn=184$ , $MinSamples2=22$ , and $\beta =0.1346$ , where $MinSamples2$ is the minimum sample size of LRDT model.

Based on the above models, the test is performed on the testing set, and repeated 20 times. One of the fitting curves and total statistical results are as follows.

It can be seen from Figure 15 that the fitting effect of the proposed method is the best based on different methods for the testing set. At the same time, it shows that the overall performance of the proposed method is the best one.

The calculation results of performance evaluation indexes of different methods are as follows.

Table 2 show the statistical results of performance evaluation indexes of different methods in terms of training set and testing set with 20 repeated times. The proposed BO-EnTSFRT-LSTM model has the lowest average RMSE and MAE for the training set (2.6610E-02, 2.0684E-02) and the testing set (4.3991E-01, 3.1771E-01). At the same time, the proposed method has the highest R² on the training set (9.9896E-01) and the testing set (7.5649E-01±2.4804E-04). By comparing the method proposed in this article with EnTSFRT-LSTM and BO-ENTSFRT, the important role of BO algorithm and LSTM-based compensation model is illustrated respectively. In addition, the comparison between the proposed method and the BO-EnTSFRT-LRDT further illustrates the superiority of the LSTM algorithm. The proposed BO-EnTSFRT-LRDT method utilizes EnTSFRT to handle uncertainty in industrial modeling data, fits error data using LSTM, and avoids hyperparameters that are prone to overfitting using BO, resulting in optimal modeling performance. By comparing with other methods, the superiority of the proposed method is proved.

Through the above comparative experiments and causal analysis, it indicates that the proposed method has the best performance compared to BO-EnTSFRT-LSTM-based modeling methods used in this study.

3.4. Hyperparameter Discussion

In the article, the single factor sensitivity analysis of the hyperparameters optimized by BO algorithm is carried out. Since the objective function of the BO algorithm is the negative value of the R² of the ensemble model, the hyperparameters are changed according to a certain step size in the range of hyperparameters and negative value of R² of the testing set are recorded for analysis. In the article, a total of 10 sets of experiments are set for each hyperparameter. The interval setting is as shown in Table 3:

The selection of step size should consider the sensitivity of hyperparameters and the limitation of computing resources. In this article, the experiment is carried out in a fixed step size according to the setting interval of the hyperparameters based on author’s debug experience and PC hardware performance. For $lr$ of the main model, the initial value is set to 0.1, and the step size is 0.1; the initial value of $\lambda$ is 3 and the step size is 3; the initial value of $\mathrm{bs}$ is 5 and the step size is 5; the initial value of $LS{M}_a$ is set to 0.1, and the step size is 0.1; the initial value of $MinSamples1$ is 10 and the step size is 10; the initial value of $Tn$ is 20, and the step size is 20. For the $\beta$ of the compensation model, the initial value is 0.01 and the step size is 0.01; the initial value of $\upsilon$ is 0.05, and the step size is 0.1; the initial value of $S$ is 220, and the step size is 400. In the experiment, the single factor analysis strategy is used. By performing 20 repeated experiments, the mean and variance change curves of RMSE on the testing set are shown in Figure 16.

Figure 16 shows a single factor-based hyperparameter analysis for mean and variance of R² in terms of average value of 20 times, and the results indicate that different hyperparameters of the constructed BO-EnTSFRT-LSTM model have different effects. The specific detailed analysis is as follow:

1) Figure 16(a) shows the effect of $lr$ on the performance of the BO-EnTSFRT-LSTM model. When $lr$ is less than 0.8, the mean value of R² fluctuates in a small range. When $lr$ is greater than 0.8, R² decreases rapidly. It shows that it is easier to meet the requirements when $lr$ is small. If the learning rate is too large, it may lead to skipping the optimal solution, which makes the model unable to achieve the best performance.

2) Figure 16(b) shows the effect of $\lambda$ on the performance of the BO-EnTSFRT-LSTM model. When $\lambda$ is in the range of 0 to 15, the mean value of R² generally shows a rapid growth. And when its value is in the range of 15 to 27, only a small fluctuation is produced. At the same time, the variance in this range also shows an overall downward trend. It shows that it is easier to meet the conditions when the value of $\lambda$ is larger.

3) Figure 16(c) shows the effect of $bs$ on the performance of the BO-EnTSFRT-LSTM model. With the increase of $bs$ , the mean and variance of R² have relatively large fluctuations. It shows that the influence of $bs$ on the performance of the model is more complicated and needs to be determined in combination with other hyperparameters.

4) Figure 16(d) shows the effect of $LS{M}_a$ on the performance of the BO-EnTSFRT-LSTM model. With the increase of $LS{M}_a$ , the mean value of R² increases rapidly. When $LS{M}_a$ is in the range of 0.3 to 0.7, its value shows a decreasing trend as a whole. When $LS{M}_a$ is 0.7, the mean value of R² is the smallest (0.7029). And then the mean value of R² increases rapidly. At the same time, the variance of R² also shows a downward trend.

5) Figure 16(e) shows the effect of $MinSamples1$ on the performance of the BO-EnTSFRT-LSTM model. When the value of $MinSamples1$ is greater than 30, the mean and variance of R² fluctuate in a small range. And when $MinSamples1=40$ , the mean value of R² is the largest (0.7077) and the variance of R² is small (0.0008209). It shows that the larger $MinSamples1$ is easier to obtain higher performance model.

6) Figure 16(f) shows the effect of $\mathrm{Tn}$ on the performance of the BO-EnTSFRT-LSTM model. When $Tn$ is in the range of 20 to 80 and 100 to 160, the mean value of R² decreases rapidly. And when tn is 20, the mean value of R² is the largest (0.7096) and the variance of R² is small (0.0008537). It shows that it is easier to obtain a high performance model when $Tn$ is small.

7) Figure 16(g) shows the effect of $\beta$ on the performance of the BO-EnTSFRT-LSTM model. When the value of $\beta$ is 0.02, the mean value of R² is largest (0.7206) and the variance value of R² is smallest (0.0006435). When the value of $\beta$ is greater than 0.07, the learning rate is too large, resulting in gradient explosion. Therefore, in order to avoid the occurrence of gradient explosion, a smaller $\beta$ should be selected.

8) Figure 16(h) shows the effect of $\upsilon$ on the performance of the BO-EnTSFRT-LSTM model. When $\upsilon$ is less than 0.8, the mean fluctuation range of R² is small. When $\upsilon$ is 0.5, the mean value of R² is the largest (0.7152) and the corresponding variance is the smallest (0.0004927). And when $\upsilon$ is greater than 0.8, the mean value of R² decreases rapidly.

9) Figure 16(i) shows the effect of $S$ on the performance of the BO-EnTSFRT-LSTM model. When $S$ is in the range of 220 to 340, the mean value of R² fluctuates in a small range. The mean value of R² increases rapidly when $S$ is in the range of 340 to 380. And when $S$ is 380, the mean value of R² is the largest (0.7306).

It should be noted that the sensitivity analysis conducted here is based on single factor hyperparameters. In fact, these numerous hyperparameters are interdependent, which is one of the reasons why BO is used for hyperparameter selection in this article.

3.5. Comprehensive Analysis

For the method proposed in this article, there are some limitations in each stage as follows.

(1) When the main model is constructed, it consists of multiple decision trees, and the training and prediction process consumes a lot of computing resources. Further research is needed on how to perform structural reduction on the constructed decision tree based main model. Moreover, how to further enhance the model’s ability to represent uncertainties in the data also needs to be strengthened.

(2) When constructing the compensation model, the problem of gradient disappearance or gradient explosion may occur. At the same time, the structure of LSTM is complex and the training process is slow. Therefore, the structure of the model can be improved by researching its lightweight algorithms in the future.

(3) In the process of hyperparameters optimization, it is easy to fall into local optimum. For high-dimensional and non-convex functions with unknown smoothness and noise, it is often difficult for BO algorithm to fit and optimize them. How to improve the BO algorithm to enhance its ability and adaptability to the above problems is a problem that needs to be solved. In addition, further consideration is needed on how to perform BO optimization on non-gaussian distributed data.

4. Conclusion

The stable control of flue gas oxygen content during MSWI process is beneficial to improve incineration efficiency and reduce pollutant emissions. Aiming at the problem of lack of high-precision and interpretable flue gas oxygen content model, the article proposes a flue gas oxygen content model based on BO-optimized main-complement ensemble algorithm. The main contributions are as follows: (1) A hybrid model integrating parallel fusion and main-complement ensemble is proposed. (2) The strategy of using BO to optimize the hyperparameters of the main-complement ensemble model is proposed. (3) An interpretable and high-precision controlled object model for flue gas oxygen content is constructed for the first time. Experiments show that the performance of the proposed method is better than other methods, which shows the feasibility and effectiveness of the proposed method. In the subsequent research, the control algorithm will be carried out on this basis to ensure the efficient and stable operation of the MSWI process.

In addition, there are some limitations in this article. Due to the use of TSFRT and LSTM, this model consumes a lot of computer resources in its training stage. At the same time, how to improve the ability and adaptability of the BO algorithm to the above problems is a problem that needs to be solved. In addition, further consideration is needed on how to perform BO optimization on non-gaussian distribution industrial data.