Recurrent Neural Network-Based Hybrid Modeling Method for Digital Twin of Boiler System in Coal-Fired Power Plant

Abstract

Due to the simplified assumptions or unascertained equipment parameters, traditional mechanism models of boiler system in coal-fired power plant usually have predictive errors that cannot be ignored. In order to further improve the predictive accuracy of the model, this paper proposes a novel recurrent neural network-based hybrid modeling method for digital twin of boiler system. First, the mechanism model of boiler system is described through recurrent neural network (RNN) to facilitate training and updating parameters, while the interpretability of the model does not degenerate. Second, for the time-varying parameters in the mechanism model, the functional relationship between them and the state variables is constructed by neurons to improve the predictive accuracy. Third, the long short-term memory (LSTM) neural network model is established to describe the unascertained dynamic characteristics to compensate the predictive residual of the mechanism model. Fourth, the update architecture and training algorithm applicable to the hybrid model are established to realize the iterative optimization of model parameters. Finally, experimental results show that the hybrid modeling method proposed in this paper can improve the predictive performance of traditional models effectively.

1. Introduction

With the development of the Industry 4.0 Plan in various countries, it is inevitable to establish digital intelligent factories by using advanced methods and technologies such as big data, artificial intelligence and industrial Internet [1,2,3,4]. Digital twin, first proposed by Grieves in 2003 [5], is an effective way for interconnection between information space and the physical world, which can realize the deep integration of information technology and traditional industry [6,7,8]. Originally, NASA used digital twin to monitor and predict the flight status of spacecraft to help engineers make correct decisions [9]. Since then, digital twin has been gradually receiving attention from researchers in various fields, especially in the manufacturing industry such as aerospace [10,11,12], automotive [13,14,15] and electronic manufacturing [16,17,18]. According to the statistical results of Psarommatis [19], we have the following information. Frist, the proportion of digital twin research is approximately 24% for aerospace and automotive, about 10% for electronics and component manufacturing, and less than 1% for energy industry. This indicates that few references study digital twin for coal-fired power plants. Second, approximately 55.5% of the references focus on ways to improve and ensure the quality of the production process and products by using digital twin, and about 18% of the existing references study the applications of digital twin for process optimization and production control. These issues are often highly correlated with the accuracy of the digital twin model. Thereby, the high-fidelity model is the solid foundation and powerful support for improving system control performance, optimizing production strategies, predicting equipment health status, increasing efficiency and reducing pollutant emissions [20,21,22,23,24,25]. Although the digital twin architectures [26,27,28], hybrid modelling of thermal systems [29] and modelling of steam turbine [30] for coal-fired power plants are gradually being proposed, the further in-depth and enriched research is still needed on digital twin modeling for boiler system in coal-fired power plants.

Guided by the aforementioned obstacles and gaps, in this paper, we propose a novel hybrid modeling method to improve the accuracy of digital twin model for boiler system. The main contributions of this paper lie in the following four aspects. (1) On the basis of ensuring the interpretability of the mechanism model, in order to facilitate training and updating parameters, this paper uses the recurrent neural network to describe the mechanism dynamic characteristics of the boiler system. (2) By adding new neurons to the mechanism model, the functional relationship between time-varying parameters and input variables is properly described to decrease the errors caused by parameters uncertainty. (3) For the sake of reducing the predictive residual of the mechanism model, the LSTM-based error-compensation model is formed. (4) To achieve the iterative optimization of hybrid model parameters, the update architecture and training algorithm are established. Finally, the experimental results show that the proposed method could ensure the high fidelity of the digital twin hybrid model.

The rest of this paper is organized as follows. In section “State of the art”, the existing modeling methods and technologies is discussed in detail. The hybrid modeling method is established and described in detail in section “Recurrent neural network-based hybrid modeling method for digital twin of boiler system”. Section “Simulation results” provides the simulation results to illustrate the effectiveness of the proposed method. The conclusions and future work are discussed in section “Conclusion”.

2. State of the Art

Presently, the main modeling methods include mechanism analysis, data-driven identification, and the combination of the two. First, as the traditional modeling method, mechanism analysis can well interpret the internal operating principle of the system and fully reflect the physical characteristics of the plant [31,32]. Unfortunately, coal-fired power plant is a huge and complex system with characteristics such as large inertia, large delay, non-linearity and strong coupling due to many equipment and mutual coupling of operation processes [33,34,35]. It is inevitable to propose numerous simplified assumptions, including simplified assumptions of equipment structure, working fluid status and lumped equipment parameters, which doubtlessly leads to poor predictive accuracy of the model [36]. In addition, due to slagging or oxidation, some parameters such as heat transfer coefficient and pipeline damping coefficient may change with equipment operation, deviating from the nominal values provided by the equipment manufacturers. This may affect the accuracy of the model more or less [37]. Second, the data-driven modeling method, also called the black box modeling method, learns the dynamic characteristic information of the plant through the input and output data [38]. Theoretically, if the number of neurons is sufficient, the neural network model can approximate the non-linear function with arbitrary accuracy [39]. In recent years, based on the data-driven method, the coal-fired boiler combustion model [40,41], unburned carbon model of boiler [42] and output power prediction model [43,44] have all achieved significant improvement. However, the accuracy of the black box model depends too much on the data quality. Insufficient data, excessive data noise and incomplete coverage of datasets may lead to poor generalization performance of the model [45]. Furthermore, the parameters are not interpretable and lack the mechanism formula to describe the relationship between variables [46].

Obviously, only using the mechanism modeling or data-driven method is hard to meet the requirements of digital twin for high-fidelity models. Compared with the above two methods, the hybrid modeling scheme can not only ensure the physical explanations, but also improve the predictive accuracy of the digital twin model. Presently, there are two main approaches. One is to establish mechanism model and identify parameters through input and output data [47,48,49]. This method can compensate most of the errors caused by the unascertained parameters, but it is powerless for the model errors introduced by mechanism simplification. The other is to compensate all errors of mechanism model by using learning algorithms such as neural network [37,50,51]. Although this method could improve the predictive accuracy of the digital twin hybrid model, it ignores the improvement of the predictive performance for the mechanism model.

To deal with the above problems while ensuring the high fidelity of the digital twin model, this paper proposes a novel recurrent neural network-based hybrid modeling method which updates the mechanism parameters and data model parameters together. The main work of this paper includes the following four parts. First, the mechanism model of boiler system is established and then described through recurrent neural network (RNN) structure. In this way, the parameters in the mechanism model can be updated iteratively through the backpropagation algorithm, while the interpretability does not deteriorate. Second, the time-varying parameters in the mechanism model change with the change in input and state variables. By using neurons, the functional relationship between time-varying parameters and input variables is properly described to improve the predictive accuracy. Third, the long short-term memory (LSTM) neural network model is established to compensate the predictive residual of the mechanism model. Due to simplified assumptions, the predictive error caused by unmodeled dynamics cannot be reduced by updating the parameters in mechanism model only. Therefore, it is an effective solution to build a data-driven model to compensate the predictive residual of mechanism model. Fourth, the hybrid model-oriented update architecture and back-propagation through time (BPTT) training algorithm are established to realize the iterative optimization of parameters. Finally, the experimental results show that the RNN-based hybrid modeling method proposed in this article has better predictive performance than the mechanism model and the black box model.

3. Recurrent Neural Network-Based Hybrid Modeling Method for Digital Twin of Boiler System

In this section, the RNN-based hybrid modeling method for digital twin of boiler system is established and described in detail, which includes the design process of digital twin model (DTM) for boiler system, the construction of the overall DTM framework for the boiler system, definitions of input and output parameters for DTM and the establishment of DTM.

3.1. The Design Process of Digital Twin Model for Boiler System

To make it easier for readers to understand the design ideas of DTM for boiler system, we present the design process of DTM clearly and completely according to the unified DT design method proposed in [19], as shown in Figure 1. First, in Section 3.2.1, the purpose of the DTM is defined, and then the DTM framework with basic information is constructed in Section 3.2.2. Second, in Section 3.3, we define the input and output parameters of DTM according to the mechanism analysis of boiler system. Third, the DTM of boiler system is established through the above information and suitable technologies in Section 3.4. Fourth, based on the existing data, historical data and real-time data, we develop the DTM and optimize corresponding parameters in Section 4.1. Finally, in Section 4.2, we evaluate and analyze whether the performance of DTM is acceptable.

3.2. The Construction of DTM Framework for Boiler System

3.2.1. The Purpose of DTM

Due to the complexity of the boiler system in coal-fired power plant, it is inevitable to create simplified assumptions before mechanism modeling to ignore the high-order dynamics (usually called unmodeled dynamics) in system. Moreover, some parameters often cannot be known accurately, which may affect the predictive precision of the mechanism model. Therefore, it is necessary to establish a more accurate predictive model to ensure the high fidelity of the digital twin model. The structure of actual system can be described as

$$F(\mathit{x},\mathit{u},\mathit{p})=f_M(\mathit{x},\mathit{u},\mathit{p})+f_P(\mathit{x},\mathit{u})+f_W(\mathit{x},\mathit{u}),$$

(1)

where $\mathit{x}$, $\mathit{u}$ and $\mathit{p}$ are state vector, input vector and parameter vector, respectively.

$F(x,\mathit{u},\mathit{p})$ represents the model of actual system, $f_M(\mathit{x},\mathit{u},\mathit{p})$ is mechanism model, $f_P(\mathit{x},\mathit{u})$ represents the errors caused by inaccurate parameters, $f_W(\mathit{x},\mathit{u})$ represents the unmodeled dynamics of the system ignored due to the simplified assumptions.

The main purpose of hybrid modeling is to calculate error terms $f_P(\mathit{x},\mathit{u})$ and $f_W(\mathit{x},\mathit{u})$ by using measurement data information, and then achieve high-precision estimation of the digital twin model $F_D(\mathit{x},\mathit{u},\mathit{p})$ to the actual model $F(x,\mathit{u},\mathit{p})$. The two existing hybrid modeling methods are shown in Figure 2a,b, respectively. The function of the method in Figure 2a is to establish a mechanism model, and then obtain the estimated error term $f_P^E(\mathit{x},\mathit{u})$ through parameter identification. However, this method is powerless for the error term $f_W(\mathit{x},\mathit{u})$. The function of the method in Figure 2b is to establish a data-driven model to compensate error terms $f_P(\mathit{x},\mathit{u})$ and $f_W(\mathit{x},\mathit{u})$. Nevertheless, it neglects to improve the predictive performance of the mechanism model. This paper combines the above two methods to establish a hybrid model, as shown in Figure 2c. Most of the error term $f_P(\mathit{x},\mathit{u})$ is compensated by training mechanism model parameters, and then the data model $f_W^E(\mathit{x},\mathit{u})$ is established to approximate the error term $f_W(\mathit{x},\mathit{u})$. This method not only improves the predictive performance of the mechanism model, making it more interpretable, but also enhances the predictive accuracy of the overall digital twin model. Thereby, the digital twin model can be written as

$$F_D(\mathit{x},\mathit{u},\mathit{p})=f_M(\mathit{x},\mathit{u},\mathit{p})+f_P^E(\mathit{x},\mathit{u})+f_W^E(\mathit{x},\mathit{u}).$$

(2)

3.2.2. The DTM Framework with Basic Information

The purpose of the DTM established in this paper is to reflect the dynamic processes of the physical boiler system in real time accurately, as shown in Figure 3. The physical boiler system includes pulverizing system, furnace and drum. Initially, the pulverized coals in coal mills are sent into the furnace by the primary air and mixed with the secondary air for combustion. Then, the water in the drum is heated to saturated vapor through the radiation energy generated by fuel combustion in furnace. Finally, the saturated vapor flowing out of the drum is heated into super-heated vapor in the super heater. The inputs of each device in Figure 3 include not only the outputs of the front device, but also the outputs of other devices.

In DTM of boiler system, the dynamic model corresponding to each device is established. The usage mode of DTM is continuous because the dynamic processes of real physical system need to be reflected synchronously. To identify the suitable technologies for the DTM development, we need to propose the requirements for DTM. The first issue to address is that of the way to achieve iterative updating of DTM parameters using input and output data while ensuring the interpretability of DTM. The solution proposed in this paper is to establish the mechanism model of the boiler system and describe it through neural network structure. In this way, the parameters in the mechanism model can be updated iteratively through the backpropagation algorithm, while the interpretability will not deteriorate. Since the outputs of the previous time need to be used, it is a better way to describe the mechanism model through recurrent neural network (RNN). The second requirement for DTM is to ensure the high fidelity after adopting the simplified mechanism model. Then, the LSTM-based error-compensation model is formed to reduce the predictive residual of the mechanism model.

3.3. Definitions of Input and Output Parameters for DTM

In this subsection, we define the input and output for DTM, and then establish the discrete-time model of each block (including pulverizing system, furnace and drum) based on the mechanism analysis in [52] through the Euler method. The input and output parameters are shown in

Table 1. Input and output parameters for DTM.

Inputs	Meanings	Inputs	Meanings	Outputs	Meanings	Outputs	Meanings
$N_g$	normalized speed of coal feeder	$Q_{net,ar}$	net calorific value	$W_{cf}$	coal flow	$Q_{sl}$	heat absorption of water wall
$W_{lk}$	primary cool air flow	$T_{slvp}$	the temperature of evaporation area	$T_o$	temperature of air–coal mixture	$M_{dl}$	the mass of liquid region in drum
$W_{rk}$	primary hot air flow	$W_e$	the inlet water flow	${\rho }_b$	gas density	${\rho }_v$	saturated vapor density
$W_{cf}$	coal flow	$H_e$	the enthalpy of inlet water	$T_{gs}$	average temperature of gas in furnace	$H_w$	the enthalpy of water in drum
$W_{sa}$	secondary air flow	$Q_{sl}$	the heat absorption of water wall	$P_b$	furnace pressure	$H_r$	the enthalpy of saturated vapor–water mixture in riser
$W_{gs}$	gas flow	$P_s$	the pressure of super-heater	$O_{cp}$	gas oxygen content

. Then, the model of each block in boiler system can be derived.

• Pulverizing system

The function of the pulverizing system is to crush the coals and sent them into furnace with the primary air. The inputs include normalized speed of coal feeder $N_g$, primary cool air flow $W_{lk}$ and hot air flow $W_{rk}$. The outputs include coal flow $W_{cf}$ and temperature of air–coal mixture $T_o$. The discrete-time model can be written as

$$W_{cf}\left(k\right)=\left(1-1/K_{cf}\right)\cdot W_{cf}\left(k-1\right)+W_g/K_{cf},$$

(3)

$$T_o\left(k\right)=T_o\left(k-1\right)+\left(Q_{ai}+Q_{rc}-Q_{ao}-Q_{mo}\right)/K_T,$$

(4)

where $Q_{ai}=H_{lk}\cdot W_{lk}+H_{rk}\cdot W_{rk}$, $Q_{rc}=H_g\cdot W_g$, $Q_{ao}=C_{pa}\cdot W_{pa}\cdot T_o\left(k-1\right)$, $W_{pa}=W_{lk}+W_{rk}$, $W_g=K_g\cdot N_g$, $Q_{mo}=C_{cf}\cdot W_{cf}\left(k-1\right)\cdot \Delta T_m$ and $\Delta T_m=(T_o\left(k-1\right)-T_g)$. $K_{cf}$ and $K_T$ are inertia time constants, $W_g$ is inlet coal flow, $Q_{ai}$, $Q_{rc}$ and $Q_{out}$ are the energy of inlet air flow, inlet coal flow and outlet mixture flow, respectively. $H_{lk}$, $H_{rk}$ and $H_g$ are the enthalpy of cool air flow, hot air flow and coals, respectively. $C_{pa}$ and $C_{cf}$ are the specific heat capacity of air and coals. $T_g$ is the inlet temperature of coals. $K_g$ is the maximum inlet coal flow. $k$ is sampling time.

2.

Furnace

In the furnace, the pulverized coals mix with the secondary air to heat the water wall through combustion. The inputs include coal flow $W_{cf}$, secondary air flow $W_{sa}$, gas flow $W_{gs}$, net calorific value $Q_{net,ar}$ and the temperature of evaporation area $T_{sv}$. The outputs include gas density ${\rho }_b$, average temperature of gas in furnace $T_{gs}$, furnace pressure $P_b$, gas oxygen content $O_{cp}$ and heat absorption of water wall $Q_{sl}$. The discrete-time model can be written as

$${\rho }_b\left(k\right)={\rho }_b\left(k-1\right)+\left(W_{cf}\left(k\right)+W_{sa}\left(k\right)-W_{gs}\left(k\right)\right)/V_b,$$

(5)

$$T_{gs}\left(k\right)=T_{gs}\left(k-1\right)+\left(Q_{ci}+Q_{si}-Q_{go}-Q_{sl}\right)/K_{bq},$$

(6)

$$P_b\left(k\right)=P_b\left(k-1\right)+\left(W_{cf}\left(k\right)+W_{sa}\left(k\right)-W_{gs}\left(k\right)\right)/C_b,$$

(7)

$$O_{cp}\left(k\right)=O_{cpi}/K_o+\left(1-1/K_o\right)\cdot O_{cp}\left(k-1\right),$$

(8)

$$Q_{sl}\left(k\right)=K_{sl}\cdot \Delta T_s/K_{sq}+\left(1-1/K_{sq}\right)\cdot Q_{sl}\left(k-1\right),$$

(9)

where $Q_{ci}=Q_{net,ar}\left(k\right)\cdot W_{cf}\left(k\right)$, $Q_{si}=H_{sa}\cdot W_{sa}\left(k\right)$, $Q_{go}=C_{gs}\cdot W_{gs}\left(k\right)\cdot T_{gs}\left(k-1\right)$, $V_{sa}=W_{sa}/{\rho }_a$, $O_{cpi}=\left(V_{sa}-V_0\cdot W_{cf}\right)\cdot 21/V_{sa}$ and $\Delta T_s=T_{gs}^4\left(k-1\right)-T_{sv}^4\left(k\right)$. $V_b$ is furnace volume. $K_{bq}$, $K_o$ and $K_{sq}$ are inertia time constants. $H_{sa}$ is the enthalpy of secondary air flow, $C_{gs}$ is the specific heat capacity of gas, $C_b$ is flow capacity coefficient of furnace, $O_{cpi}$ is the inlet gas oxygen content at furnace bottom, $Q_{sc}$ is the radiant energy near fuel combustion, $V_{sa}$ is the volume flow of secondary air, $V_0$ is the theoretical air consumption, $K_{sl}$ is the radiant heat transfer coefficient, ${\rho }_a$ is the air density and 21 is the oxygen content of air.

3.

Drum

In the drum, the riser absorbs the energy from the furnace to heat the water into saturated vapor, and then the saturated vapor flows out to super-heater through the top of the drum. The inputs include the inlet water flow $W_e$ from economizer, the enthalpy of inlet water $H_e$, the heat absorption of water wall $Q_{sl}$ and the pressure of super-heater $P_s$. The outputs include the mass of liquid region in drum $M_{dl}$, saturated vapor density ${\rho }_v$, the enthalpy of water in drum $H_w$ and the enthalpy of saturated vapor-water mixture in riser $H_r$. The discrete-time model can be written as

$$M_{dl}\left(k\right)=M_{dl}\left(k-1\right)+W_e\left(k\right)-q_v\cdot W_{ro}-W_{dv},$$

(10)

$${\rho }_v\left(k\right)={\rho }_v\left(k-1\right)+V_v^{-1}\cdot \Delta W_v,$$

(11)

$$H_w\left(k\right)=H_w\left(k-1\right)+\left(Q_e+Q_{wv}-Q_d-Q_{dv}\right)/M_{dl},$$

(12)

$$H_r\left(k\right)=H_r\left(k-1\right)+\left(Q_{sl}+Q_d-Q_{ro}\right)/K_r$$

(13)

where $q_v=F_q\left(H_r,H_v,H_{wv}\right)=\left(H_r-H_{wv}\right)/\left(H_v-H_{wv}\right)$, $V_v=V_{drum}-\left(M_{dl}/{\rho }_w\right)$, $W_{vp}=q_v\cdot W_{ro}$, $W_{dv}=K_{ec}\cdot \left(P_{dr}-P_{dr}^{ex}\right)$, $H_{wv}=Tp\left({\rho }_v\right)$, $H_v=Tp\left({\rho }_v\right)$, ${\rho }_w=Tp\left({\rho }_v\right)$, $P_{dr}=Tp\left({\rho }_v\right)$, $Q_{dv}=W_{dv}\cdot H_v$, $\Delta W_v=W_{vp}+W_{dv}-W_v$, $W_v=\sqrt{\left(P_{dr}-P_s\right)/R_f}$, $Q_e=H_e\left(k\right)\cdot W_e\left(k\right)$, $Q_{wv}=H_{wv}\cdot \left(1-q_v\right)\cdot W_{ro}$, $Q_d=H_w\left(k-1\right)\cdot W_{ro}$ and $Q_{ro}=H_r\left(k-1\right)\cdot W_{ro}$. $q_v$ is the vapor content of vapor–water mixture at riser outlet, $W_{ro}$ is the vapor–water mixture flow of riser outlet, $V_v$ is the volume of vapor region in drum, $W_{dv}$ is the dynamic evaporation in drum, $W_v$ is the outlet vapor flow, $R_f$ is the pipe damping coefficient, $H_{wv}$ is the enthalpy of saturated water, $W_d$ is inlet water flow of down tube, $K_r$ is the inertia time constants of riser, $V_{drum}$ is the volume of drum, ${\rho }_w$ is water density in drum, $H_v$ is the enthalpy of saturated vapor, $K_{ec}$ is the coefficient, $P_{dr}$ is the drum pressure, $P_{dr}^{ex}$ is the drum pressure of the previous time, $Tp(\cdot )$ represents thermodynamic property function of saturated vapor and saturated water. The $P_{dr}$, $H_{wv}$, $H_v$ and ${\rho }_w$ can be derived by $Tp\left({\rho }_v\right)$. The $Q_e$, $Q_{wv}$, $Q_{dv}$, $Q_d$ and $Q_{ro}$ are the energy of inlet water flow, saturated water, dynamic vapor, inlet water flow of down tube and vapor-water mixture flow of riser outlet, respectively. Additionally, to simplify the drum model, we neglect the dynamic equation of the down tube. That is, the inlet water flow of the down tube is equal to the outlet flow of riser.

3.4. Establishment of Digital Twin Model for Boiler System

3.4.1. Recurrent Neural Network Description of Boiler System

Because some parameters in mechanism model cannot be known accurately, it is necessary to update the parameters iteratively by using learning algorithms to approximate the true values continuously. If we describe the mechanism model through neural network structure, the parameters can be easily trained through backpropagation algorithm. Before building the neural network of mechanism model, the writing rules of neuron symbols and weights in the network should be defined. Let the $w_G^m{}_{hu}^{(l)}$ represent the general form of weights, where $G$ represents the equipment block (i.e., the pulverizing block $pz$, furnace block $bo$ and drum block $vp$), $m$ represents the order number of outputs corresponding to the block, $h$, $u$ and $\left(l\right)$ represent the inputs order number of the neuron, the order number of the neuron in the corresponding layer, and the order number of layers corresponding to the neuron, respectively. Briefly speaking, the $w_G^m{}_{hu}^{(l)}$ represents the weight of the $h$th input of the $u$th neuron in the $l$th layer of the $m$th output of the block $G$. Similarly, the $b_G^m{}_u^{\left(l\right)}$ represents the bias of the $u$th neuron in the $l$th layer of the $m$th output of the block $G$. In addition, $w_G^0{}_n$ represents the weight shared by each output of block $G$, where $n$ is the order number of the weight. In the neurons, the symbols in the left represent the operation on the weighted input. The specific meanings are shown in

Table 2. Meanings of symbols in neurons.

Symbols	Operation	Symbols	Operation	Symbols	Operation
$\sum$	Summation	$\prod$	Product	$Pw$	Exponential
$Tp$	Thermodynamic property functions	$Fq$	Function of vapor content of vapor–water mixture	$Sq$	Square root

. The symbols on the right represent the output of neuron after calculation.

• Pulverizing system

According to Equations (3) and (4), the neural network structures of the pulverizing block are shown in Figure 4. In the two figures, weights mean the model parameters, including known parameters and unascertained parameters. The known weights and unascertained weights are shown in

Table 3. The known weights and unascertained weights in pulverizing block.

Known	Values	Known	Values	Unascertained	Values	Unascertained	Values
$w_{pz}^2{}_{11}^{(1)}$	$H_{lk}$	$w_{pz}^2{}_{13}^{(2)}$	$H_g$	$w_{pz}^0{}_1$	$K_g$	$w_{pz}^2{}_{11}^{(4)}$	$1/K_T$
$w_{pz}^2{}_{21}^{(1)}$	$H_{rk}$	$w_{pz}^2{}_{21}^{(3)}$	$-1$	$w_{pz}^1{}_{11}^{(2)}$	$1/K_{cf}$	$w_{pz}^2{}_{14}^{(2)}$	/
$b_{pz}^2{}_3^{(1)}$	$-T_g$	$w_{pz}^2{}_{31}^{(3)}$	$-1$	$w_{pz}^1{}_{21}^{(2)}$	$1-w_{pz}^1{}_{11}^{(2)}$	$b_{pz}^2{}_4^{(2)}$	/
$w_{pz}^2{}_{11}^{(2)}$	$C_{pa}$			$w_{pz}^2{}_{12}^{(2)}$	$C_{cf}$

. In Equation (4), the heat dissipation of pulverizing system is ignored for simplification. Therefore, neuron $Q_{bu}^{pz}$ is added to represent the heat dissipation and is assumed to be a function of the inlet coal flow $W_g$, that is, $Q_{bu}^{pz}=w_{pz}^2{}_{14}^{(2)}\cdot W_g+b_{pz}^2{}_4^{(2)}$. The value or initial value of the parameters is provided in the “Simulation results” section. In addition, it should be noted that all unmarked weights in the networks are equal to 1.

2.

Furnace

Based on Equations (5)–(9), the neural network structures of the furnace block are shown in Figure 5, Figure 6 and Figure 7. The known weights and unascertained weights are shown in

Table 4. The known weights and unascertained weights in furnace block.

Known	Values	Known	Values	Unascertained	Values	Unascertained	Values
$w_{bo}^0{}_1$	$-1$	$w_{bo}^2{}_{13}^{\left(1\right)}$	$H_{sa}$	$w_{bo}^1{}_{11}^{\left(2\right)}$	$1/V_b$	$w_{bo}^2{}_{14}^{\left(1\right)}$	$C_{gs}$
$w_{bo}^2{}_{32}^{\left(2\right)}$	$-1$	$w_{bo}^2{}_{42}^{\left(2\right)}$	$-1$	$w_{bo}^3{}_{11}^{\left(2\right)}$	$1/C_b$	$w_{bo}^4{}_{11}^{\left(2\right)}$	$V_0$
$w_{bo}^4{}_{21}^{\left(2\right)}$	$-21\cdot {\rho }_{sa}$	$b_{bo}^4{}_1^{\left(2\right)}$	$21$	$w_{bo}^4{}_{11}^{\left(3\right)}$	$1/K_o$	$w_{bo}^4{}_{21}^{\left(3\right)}$	$1-w_{bo}^4{}_{11}^{\left(3\right)}$
$w_{bo}^5{}_{21}^{\left(2\right)}$	$-1$			$w_{bo}^5{}_{11}^{\left(4\right)}$	$1/K_{sq}$	$w_{bo}^5{}_{21}^{\left(4\right)}$	$1-w_{bo}^5{}_{11}^{\left(4\right)}$
				$w_{bo}^2{}_{11}^{\left(2\right)}$	/	$w_{bo}^5{}_{11}^{\left(1\right)}$	/
				$b_{bo}^5{}_1^{\left(1\right)}$	/

. In furnace block, generally speaking, $K_Q=1/K_{bq}$ and $K_{sl}$ are time-varying parameters, and their specific dynamic characteristics and functional relationship cannot be known accurately. Thus, by adding neurons to the network, these two parameters are expressed as $K_Q=w_{bo}^2{}_{11}^{\left(2\right)}\cdot {\rho }_b^{-1}$ and $K_{sl}=w_{bo}^5{}_{11}^{\left(1\right)}\cdot W_{gs}+b_{bo}^5{}_1^{\left(1\right)}$, respectively.

3.

Drum

Based on Equations (10)–(13), the neural network structures of the drum block are shown in Figure 8, Figure 9 and Figure 10. The known weights and unascertained weights are shown in

Table 5. The known weights and unascertained weights in drum block.

Known	Values	Known	Values	Known	Values	Unascertained	Values
$w_{vp}^0{}_1$	$-1$	$w_{vp}^0{}_2$	$K_{ec}$	$w_{vp}^1{}_{21}^{\left(4\right)}$	$-1$	$w_{vp}^0{}_3$	$W_{ro}$
$w_{vp}^1{}_{31}^{\left(4\right)}$	$-1$	$w_{vp}^2{}_{11}^{\left(3\right)}$	$-1$	$b_{vp}^2{}_1^{\left(3\right)}$	$V_{drum}$	$w_{vp}^2{}_{14}^{\left(3\right)}$	$1/R_f$
$w_{vp}^2{}_{24}^{\left(2\right)}$	$-1$	$w_{vp}^2{}_{32}^{\left(4\right)}$	$-1$	$w_{vp}^3{}_{12}^{\left(3\right)}$	$-1$	$w_{vp}^4{}_{11}^{\left(3\right)}$	$1/K_r$
$b_{vp}^3{}_2^{\left(3\right)}$	$1$	$w_{vp}^3{}_{21}^{\left(5\right)}$	$-1$	$w_{vp}^3{}_{41}^{\left(5\right)}$	$-1$
$w_{vp}^4{}_{11}^{\left(2\right)}$	$-1$

. The specific thermodynamic functions $Tp(\cdot )$ are provided in the “Simulation results” section.

3.4.2. LSTM-Based Error Compensation Model for Digital Twin

As previously mentioned, the predictive error caused by unmodeled dynamics cannot be reduced by updating the parameters in mechanism model only. Moreover, in practice, most of the predictive errors caused by unascertained parameters can be reduced through the mechanism parameter identification, but it is difficult to eliminate the errors completely. For this reason, it is necessary to establish a data-driven model to compensate the predictive residual of mechanism model. Long short-term memory (LSTM) network is a special kind of RNN which can solve the problem of gradient disappearance and explosion during long sequence training [53]. In other words, compared with common RNN, LSTM has better performance in long time series prediction [54,55,56].

Therefore, for each block, the error-compensation model is established based on LSTM. The inputs of the LSTM are the same as that of the corresponding block in mechanism network, and the outputs are the compensation for the predictive error of the mechanism model.

3.4.3. Parameter Updating Structure for Digital Twin Model of Boiler System

The parameter updating structure for digital twin hybrid model is shown in Figure 11. According to the input ${\mathit{x}}_G$, the actual system, mechanism model and error compensation model generate the real output ${\mathit{y}}_G$, predictive output of mechanism model ${\tilde{\mathit{y}}}_G^{mc}$ and that of error compensation model ${\tilde{\mathit{y}}}_G^{ec}$, respectively. The subscript $G$ represents the equipment block. The parameters in mechanism model are updated based on the difference between ${\tilde{\mathit{y}}}_G^{mc}$ and ${\mathit{y}}_G$, that is, ${\mathit{E}}_G^{mc}={\tilde{\mathit{y}}}_G^{mc}-{\mathit{y}}_G$. This is because we hope that the predictive output of the mechanism model can approach the real value as much as possible. For the error compensation model, we hope that it can reduce the predictive residual of the mechanism model as much as possible. Therefore, the parameters in compensation model are updated based on the difference between ${\tilde{\mathit{y}}}_G^{hy}$ and ${\mathit{y}}_G$, which can be expressed as ${\mathit{E}}_G^{ec}={\tilde{\mathit{y}}}_G^{hy}-{\mathit{y}}_G$. The advantages of this structure can be summed up in the following three aspects. First, the model parameters are updated to improve the predictive performance of mechanism model. Second, due to the improvement in predictive accuracy, the mechanism model has better interpretability. Third, the error compensation model is added and updated to ensure the high fidelity of the digital twin model.

3.4.4. Training Algorithm for Digital Twin Model of Boiler System

Because the mechanism model and error compensation model mentioned above are established based on RNN, Backpropagation Through Time (BPTT) method [57], the most common optimization algorithm in RNN [58], is adopted in this paper for iterative training and updating of network parameters in these two models. The time sequence propagation structure of the network for mechanism model is shown in Figure 12. Assume that the total length of the time sequence is $T$. The directions of the black arrow represent the forward propagation through the network and time. On the contrary, that of the red arrow represent the backpropagation of errors through the network and time. The general form of the network equation at a certain moment $t$ can be written as

$${\tilde{\mathit{y}}}_G^{mc}\left(t\right)=\mathit{f}\left({\mathit{x}}_G\left(t\right),{\tilde{\mathit{y}}}_G^{mc}\left(t-1\right),{\mathit{W}}_G,{\mathit{b}}_G\right),$$

(14)

where ${\mathit{W}}_G$ and ${\mathit{b}}_G$ are the weight and bias matrices of block $G$. Since the loss function at time $t$ is related to all previous moments, using the chain rule, the gradient of $L_G^{mc}$ with respect to ${\mathit{W}}_G$ is written as

$$\frac{\partial L_G^{mc}}{\partial {\mathit{W}}_G}={\displaystyle \sum _{t=0}^T{\displaystyle \sum _{k=0}^t\frac{\partial L_G^{mc}\left(t\right)}{\partial {\tilde{\mathit{y}}}_G^{mc}\left(t\right)}\left({\displaystyle \prod _{j=k+1}^t\frac{\partial {\tilde{\mathit{y}}}_G^{mc}\left(j\right)}{\partial {\tilde{\mathit{y}}}_G^{mc}\left(j-1\right)}}\right)\frac{\partial {\tilde{\mathit{y}}}_G^{mc}\left(k\right)}{\partial {\mathit{W}}_G}}}.$$

(15)

Similarly, we can also obtain

$$\frac{\partial L_G^{mc}}{\partial {\mathit{b}}_G}={\displaystyle \sum _{t=0}^T{\displaystyle \sum _{k=0}^t\frac{\partial L_G^{mc}\left(t\right)}{\partial {\tilde{\mathit{y}}}_G^{mc}\left(t\right)}\left({\displaystyle \prod _{j=k+1}^t\frac{\partial {\tilde{\mathit{y}}}_G^{mc}\left(j\right)}{\partial {\tilde{\mathit{y}}}_G^{mc}\left(j-1\right)}}\right)\frac{\partial {\tilde{\mathit{y}}}_G^{mc}\left(k\right)}{\partial {\mathit{b}}_G}}}.$$

(16)

Finally, the parameters can be iteratively updated through the gradient descent method, which can be expressed as ${\mathit{W}}_G={\mathit{W}}_G-\eta \cdot \partial L_G^{mc}/\partial {\mathit{W}}_G$ and ${\mathit{b}}_G={\mathit{b}}_G-\eta \cdot \partial L_G^{mc}/\partial {\mathit{b}}_G$. $\eta$ is the learning rate.

4. Simulation Results

In this section, the parameter training results of the DTM are analyzed and the predictive performance of different models is compared. Before model training, the preprocessing of data and parameters is required. It should be noted that almost all parameters in gradient terms have physical meaning, so it is necessary to determine the reasonable range for time-varying parameters that need to be updated. Additionally, due to the large difference in the order of magnitude of each neuron (i.e., state variable), it is easy to cause gradient explosion or gradient disappearance in backpropagation processes. Therefore, the range of corresponding parameters or variables in gradient terms should be determined and then their normalized values need to be derived.

4.1. Preprocessing of Data and Parameters

4.1.1. Mechanism Model

The values of parameters in three blocks are shown in

Table 6. The value of parameters in mechanism models.

Parameter	Value	Parameter	Value	Parameter	Value	Parameter	Value	Parameter	Value
$H_{lk}$	41	$H_{rk}$	276	$T_g$	30	$C_{pa}$	1.022	$C_{cf}$	1.3
$H_g$	56.4	$K_g$	100	$K_{cf}$	2.5	$K_{ec}$	−600	$H_{sa}$	200
$K_T$	1667	$w_{pz}^2{}_{14}^{(2)}$	0.001	$b_{pz}^2{}_4^{(2)}$	0.001	$K_r$	$2.97\times {10}^4$	$C_{gs}$	1.5
$V_0$	5.25	${\rho }_{sa}$	1.416	$V_b$	16,855	$b_{bo}^5{}_1^{\left(1\right)}$	$5\times {10}^{-7}$	$w_{bo}^2{}_{11}^{\left(2\right)}$	$5\times {10}^{-5}$
$C_b$	14.286	$K_o$	20	$K_{sq}$	1.25	$w_{bo}^5{}_{11}^{\left(1\right)}$	$-2.5\times {10}^{-10}$	$V_{drum}$	118.53
$W_{ro}$	1290	$R_f$	$2.27\times {10}^{-6}$

. For unascertained parameters, the initial value is considered. According to the thermodynamic properties of saturated water and saturated vapor, functions $Tp(·)$ can be written as

$$P_{dr}=A_{dr}\cdot {\rho }_v^3+B_{dr}\cdot {\rho }_v^2+C_{dr}\cdot {\rho }_v+D_{dr},$$

(17)

$$H_v=A_v\cdot {\rho }_v^3+B_v\cdot {\rho }_v^2+C_v\cdot {\rho }_v+D_v,$$

(18)

$$H_{wv}=A_{wv}\cdot {\rho }_v^4+B_{wv}\cdot {\rho }_v^3+C_{wv}\cdot {\rho }_v^2+D_{wv}\cdot {\rho }_v+E_{wv},$$

(19)

$${\rho }_w=A_w\cdot H_w^3+B_w\cdot H_w^2+C_w\cdot H_w+D_w,$$

(20)

where $A_{dr}=1.104\times {10}^{-6}$, $B_{dr}=-8.206\times {10}^{-4}$, $C_{dr}=0.2265$, $D_{dr}=-0.2319$, $A_v=2.331\times {10}^{-5}$, $B_v=-6.728\times {10}^{-3}$, $C_v=-2.151$, $D_v=2860$, $A_{wv}=-7.772\times {10}^{-7}$, $B_{wv}=4.381\times {10}^{-4}$, $C_{wv}=-0.102$, $D_{wv}=14.12$, $E_{wv}=870.3$, $A_w=-1.458\times {10}^{-7}$, $B_w=4.729\times {10}^{-4}$, $C_w=-0.8501$ and $D_w=1355$. Then, the derivatives $\partial P_{dr}/\partial {\rho }_v$, $\partial H_v/\partial {\rho }_v$, $\partial H_{wv}/\partial {\rho }_v$ and $\partial {\rho }_w/\partial H_w$ in gradient terms can be easily derived. The range of time-varying parameters, inputs and outputs are shown in

Table 7. The range of time-varying parameters, inputs and outputs.

Parameter	Range	Std. Range	Parameter	Range	Std. Range	Parameter	Range	Std. Range
$N_g$	$\left[0.54,0.75\right]$	$\left[0,1\right]$	$W_{lk}$	$\left[19,32\right]$	$\left[0,1\right]$	$W_{rk}$	$\left[78,92\right]$	$\left[0,1\right]$
$W_{cf}$	$\left[54,80\right]$	$\left[0,1\right]$	$T_o$	$\left[120,145\right]$	$\left[0,1\right]$	$w_{pz}^1{}_{11}^{\left(1\right)}$	$\left[90,105\right]$	$\left[0,1\right]$
$w_{pz}^1{}_{11}^{\left(2\right)}$	$\left[0.1,0.5\right]$	$\left[0,1\right]$	$w_{pz}^2{}_{12}^{\left(2\right)}$	$\left[1.2,2.2\right]$	$\left[0,1\right]$	$w_{pz}^2{}_{11}^{\left(4\right)}$	${\left[4,9\right]}^{{\rm E}-3}$	$\left[0,1\right]$
$w_{pz}^2{}_{14}^{\left(2\right)}$	$\left[-1,1\right]$	$\left[-1,1\right]$	$b_{pz}^2{}_4^{\left(2\right)}$	$\left[-2,2\right]$	$\left[-1,1\right]$	$W_{sa}$	$\left[500,700\right]$	$\left[0,1\right]$
$W_{gs}$	$\left[590,770\right]$	$\left[0,1\right]$	$Q_{net,ar}$	$\left[19,21\right]$	$\left[0,1\right]$	$T_{slvp}$	$\left[329,332\right]$	$\left[0,1\right]$
${\rho }_b$	$\left[1.3,1.4\right]$	$\left[0,1\right]$	$T_{gs}$	${\left[1,1.1\right]}^{{\rm E}+3}$	$\left[0,1\right]$	$P_b$	${\left[1.01,1.02\right]}^{{\rm E}+5}$	$\left[0,1\right]$
$O_{cp}$	$\left[2.6,3.8\right]$	$\left[0,1\right]$	$Q_{sl}$	${\left[5.2,6.2\right]}^{{\rm E}+5}$	$\left[0,1\right]$	$w_{bo}^1{}_{11}^{\left(2\right)}$	${\left[5,8.5\right]}^{{\rm E}-5}$	$\left[0,1\right]$
$w_{bo}^2{}_{14}^{\left(1\right)}$	$\left[1,1.6\right]$	$\left[0,1\right]$	$w_{bo}^2{}_{11}^{\left(2\right)}$	${\left[2,8\right]}^{{\rm E}-5}$	$\left[0,1\right]$	$w_{bo}^3{}_{11}^{\left(2\right)}$	$\left[0.03,0.08\right]$	$\left[0,1\right]$
$w_{bo}^4{}_{11}^{\left(2\right)}$	$\left[4.5,6\right]$	$\left[0,1\right]$	$w_{bo}^4{}_{11}^{\left(3\right)}$	$\left[0.01,0.09\right]$	$\left[0,1\right]$	$w_{bo}^5{}_{11}^{\left(1\right)}$	${\left[-5,5\right]}^{{\rm E}-10}$	$\left[-1,1\right]$
$w_{bo}^5{}_{11}^{\left(4\right)}$	$\left[0.7,1\right]$	$\left[0,1\right]$	$b_{bo}^5{}_1^{\left(1\right)}$	${\left[-6,6\right]}^{{\rm E}-7}$	$\left[-1,1\right]$	$W_e$	$\left[440,562\right]$	$\left[0,1\right]$
$H_e$	$\left[1280,1350\right]$	$\left[0,1\right]$	$P_s$	$\left[17.8,18.3\right]$	$\left[0,1\right]$	$M_{dl}$	${\left[3.78,3.87\right]}^{{\rm E}+4}$	$\left[0,1\right]$
${\rho }_v$	$\left[137,149\right]$	$\left[0,1\right]$	$H_w$	$\left[1581,1596\right]$	$\left[0,1\right]$	$H_r$	$\left[2010,2074\right]$	$\left[0,1\right]$
$w_{vp}^0{}_3$	$\left[1280,1295\right]$	$\left[0,1\right]$	$w_{vp}^2{}_{14}^{\left(3\right)}$	${\left[4.2,4.5\right]}^{{\rm E}+5}$	$\left[0,1\right]$	$w_{vp}^4{}_{11}^{\left(3\right)}$	${\left[3.2,4\right]}^{{\rm E}-5}$	$\left[0,1\right]$

. It should be noted that the range such as $\left[a\times {10}^{\pm c},b\times {10}^{\pm c}\right]$ is written as ${\left[a,b\right]}^{{\rm E}\pm c}$ for simplicity. In the training processes, the positive and negative characteristics of parameters affect the direction of gradient terms. Therefore, for parameters with different ranges, different normalization ranges and methods need to be adopted to retain their positive and negative signs. In general, the normalized interval of parameters can be divided into three categories, including $\left[0,1\right]$, $\left[-1,1\right]$ and $\left[-1,0\right]$. First, for parameter with interval $\left[a,b\right]$, where $a\ge 0$ and $b>0$, the normalized interval is set as $\left[0,1\right]$. The normalization method can be written as $\overline{x}=\left(x-a\right)/\left(b-a\right)$, where $\overline{x}$ is normalized value of $x$. Second, for parameter with interval $\left[a,b\right]$, where $a<0$, $b>0$ and $a=-b$, the normalized interval is set as $\left[-1,1\right]$. The normalization method can be written as $\overline{x}=x/b$. Third, for parameter with interval $\left[a,b\right]$, where $a<0$ and $b<0$, the normalized interval is set as $\left[-1,0\right]$. The normalization method can be written as $\overline{x}=-\left(\left|x\right|-\left|b\right|\right)/\left(\left|a\right|-\left|b\right|\right)$. In Table 7, the normalized (std.) range of parameters such as $w_{pz}^2{}_{14}^{\left(2\right)}$, $b_{pz}^2{}_4^{\left(2\right)}$, $w_{bo}^5{}_{11}^{\left(1\right)}$ and $b_{bo}^5{}_1^{\left(1\right)}$ is $\left[-1,1\right]$, and that of other parameters is $\left[0,1\right]$. Moreover, the learning rate $\eta$ is set to 0.001 and the threshold of gradient is set to one.

4.1.2. Error Compensation Model

The range of outputs in error compensation model is shown in

Table 8. The range of outputs in error compensation model.

Parameter	Range	Std. Range	Parameter	Range	Std. Range	Parameter	Range	Std. Range
$\Delta {\rho }_b^{ec}$	${\left[-5,1\right]}^{{\rm E}-3}$	$\left[0,1\right]$	$\Delta T_{gs}^{ec}$	$\left[-3,4\right]$	$\left[0,1\right]$	$\Delta P_b^{ec}$	$\left[-4.5,0.5\right]$	$\left[0,1\right]$
$\Delta O_{cp}^{ec}$	$\left[-0.04,0.02\right]$	$\left[0,1\right]$	$\Delta Q_{sl}^{ec}$	${\left[-6,8\right]}^{{\rm E}+3}$	$\left[0,1\right]$	$\Delta M_{dl}^{ec}$	$\left[-150,150\right]$	$\left[0,1\right]$
$\Delta {\rho }_v^{ec}$	$\left[-1,4\right]$	$\left[0,1\right]$	$\Delta H_w^{ec}$	$\left[-8,15\right]$	$\left[0,1\right]$	$\Delta H_r^{ec}$	$\left[-8,15\right]$	$\left[0,1\right]$

. The range of inputs is the same as that of the mechanism model, which is shown in Table 7. Since the error compensation model is established based on the general LSTM network, the normalized range can be all selected as $\left[0,1\right]$, without considering the sign of outputs. The LSTM network structure is set to 1-1-1, that is, an input layer, one hidden layer with 90 neurons and an output layer. The learning rate is set to 0.001, and the threshold of gradient is set to one. The loss function is defined as

$$L_G^{ec}={\displaystyle \sum _{t=0}^T{L}_G^{ec}\left(t\right)}={\displaystyle \sum _{t=0}^T\frac{1}{2}{\left({\tilde{\mathit{y}}}_G^{hy}\left(t\right)-{\mathit{y}}_G\left(t\right)\right)}^T\left({\tilde{\mathit{y}}}_G^{hy}\left(t\right)-{\mathit{y}}_G\left(t\right)\right)},$$

(21)

where ${\tilde{\mathit{y}}}_G^{hy}\left(t\right)={\tilde{\mathit{y}}}_G^{mc}\left(t\right)+{\tilde{\mathit{y}}}_G^{ec}\left(t\right)$. The ${\tilde{\mathit{y}}}_G^{ec}$ and ${\tilde{\mathit{y}}}_G^{hy}$ represent the predictive outputs of error compensation model and that of digital twin hybrid model of block $G$, respectively.

4.2. Model Training and Analysis of Predictive Performance

4.2.1. Model Training

The training and testing datasets used in this paper are collected from the 600 MW coal-fired power plant simulator. The training process of parameters in three blocks is shown in Figure 13. To display the parameters of each block in the same picture, the iterative curve of normalized parameters is provided. The training process of some parameters shows a trend of first falling and then rising, which may be related to the initial value of block parameters. Additionally, the convergence speed of some parameters is slower than that of other parameters, which may be related to the setting of threshold and the normalized range of derivatives in gradient terms. The convergence (con.) value, normalized (std.) value and actual value of all parameters are listed in

Table 9. The convergence value, normalized value and actual value of all parameters.

Parameter	Std. Value	Con. Value	Actual Value	Parameter	Std. Value	Con. Value	Actual Value
$w_{pz}^1{}_{11}^{\left(1\right)}$	0.4666	96.9988	97	$w_{pz}^1{}_{11}^{\left(2\right)}$	0.3694	0.2477	0.245
$w_{pz}^2{}_{12}^{\left(2\right)}$	0.8617	2.0617	1.988	$w_{pz}^2{}_{11}^{\left(4\right)}$	0.8076	$8\times {10}^{-4}$	$7.93\times {10}^{-4}$
$w_{pz}^2{}_{14}^{\left(2\right)}$	0.1774	0.1774	/	$b_{pz}^2{}_4^{\left(2\right)}$	0.1836	0.3672	/
$w_{bo}^1{}_{11}^{\left(2\right)}$	0.3447	$6.21\times {10}^{-5}$	$6.13\times {10}^{-5}$	$w_{bo}^2{}_{14}^{\left(1\right)}$	0.5169	1.3101	1.31
$w_{bo}^2{}_{11}^{\left(2\right)}$	0.1650	$2.99\times {10}^{-5}$	$3.38\times {10}^{-5}$	$w_{bo}^3{}_{11}^{\left(2\right)}$	0.4461	0.0523	0.05
$w_{bo}^4{}_{11}^{\left(2\right)}$	0.6089	5.4134	5.4135	$w_{bo}^4{}_{11}^{\left(3\right)}$	0.2476	0.0298	0.0287
$w_{bo}^5{}_{11}^{\left(1\right)}$	−0.3812	$-1.91\times {10}^{-10}$	/	$w_{bo}^5{}_{11}^{\left(4\right)}$	0.9328	0.9798	0.99
$b_{bo}^5{}_1^{\left(1\right)}$	0.9738	$5.84\times {10}^{-7}$	/	$w_{vp}^0{}_3$	0.3334	1285.01	1283.85
$w_{vp}^2{}_{14}^{\left(3\right)}$	0.3048	429144.41	429150	$w_{vp}^4{}_{11}^{\left(3\right)}$	0.8271	$3.862\times {10}^{-5}$	$3.85\times {10}^{-5}$

. Due to the simplification of mechanism model established in this paper, some corresponding parameters, especially dynamic parameters such as $K_{cf}$, $K_T$, $K_o$, $K_{sq}$ and $K_r$, cannot be detected in the simulator directly. Therefore, for these parameters, we obtain their approximate true values through the derivation and calculation of other relevant parameters in the simulator. The convergence value of most parameters is very close to the real value, except for parameters $w_{bo}^2{}_{11}^{\left(2\right)}$ and $w_{vp}^0{}_3$. The reasons for the large deviation of these two parameters may be the complex gradient terms, the inappropriate normalized range of parameters and derivatives, and the insufficient number of additional neurons. The additional neurons such as $w_{pz}^2{}_{14}^{\left(2\right)}$, $b_{pz}^2{}_4^{\left(2\right)}$, $w_{bo}^5{}_{11}^{\left(1\right)}$ and $b_{bo}^5{}_1^{\left(1\right)}$ are used for mapping the relationship between time-varying parameters and state variables. Because they have no physical meaning, there are no corresponding actual values in the simulator. With regard to the training of error compensation models, we use the deep learning toolbox of MATLAB 2022a. The training processes are shown in Figure 13.

4.2.2. Analysis of Predictive Performance

The predictive results $W_{cf}$ and $T_o$ of pulverizing system are shown in Figure 14, including the curves of actual value, mechanism model with nominal parameters (MMNP), LSTM-based empirical model (EM-LSTM), GRU-based empirical model (EM-GRU) and hybrid model with parameter identification only (HMPIO). The upper subgraphs in Figure 14 depict the overall trend of outputs. As a supplement, the lower subgraphs show the details of the corresponding output as well as the difference between models. Because the time-varying parameters in MMNP are the initial nominal values and have not been updated, there are differences between them and the true values. Undoubtedly, this leads to unsatisfactory predictive accuracy. The EM-LSTM and EM-GRU are data-driven models, which learn the dynamic characteristic information of the system through the input and output data. The HMPIO can be said to be a mechanism model in the form of neural networks, which not only has interpretability, but also makes parameters easy to train through BPTT algorithm. Compared to MMNP, HMPIO achieves higher predictive accuracy by training the parameters to approximate the true values gradually. The iterative training results of the parameters are shown in Figure 13. According to the evaluation indexes proposed by Psarommatis [59], the performance of digital twin model can be analyzed comprehensively. The results are listed in

Table 10. The evaluation indexes of digital twin model for pulverizing system.

Indexes	Output	MMNP	EM-LSTM	EM-GRU	HMPIO
AOP (%)	$W_{cf}$	85.51	2.66	1.41	0.06
	$T_o$	75.01	0.2	0.5	0.17
ACVAR	$W_{cf}$	0.0438	0.0391	0.0394	0.0363
	$T_o$	0.6771	0.0147	0.0149	0.0148
GDTA (%)	global	80.26	1.43	0.95	0.12
GAVAR	global	0.3605	0.0269	0.0271	0.0256
ART (s)	global	$2.84\times {10}^{-5}$	0.0044	0.0041	$7.43\times {10}^{-6}$
RTVAR	global	$2.20\times {10}^{-7}$	$2.27\times {10}^{-5}$	$9.29\times {10}^{-6}$	$5.23\times {10}^{-8}$

. The indexes include the accuracy of each output parameter (AOP), global DT accuracy (GDTA), the accuracy variance of each output parameter (ACVAR), the global accuracy variance (GAVAR), the response time (ART) and the response time variance (RTVAR). From the results in Table 10, for pulverizing system, the HMPIO has the highest predictive accuracy and fastest response time. In addition, for this block, the predictive accuracy of HMPIO can meet the requirements, so no error compensation model (ECM) is added. The predictive results of furnace, including ${\rho }_b$, $T_{gs}$, $P_b$, $O_{cp}$ and $Q_{sl}$, are shown in Figure 15 and Figure 16. The predictive error of MMNP is much larger than that of other methods, especially for $O_{cp}$ and $Q_{sl}$. As for the results of HMPIO for variables ${\rho }_b$ and $P_b$, there are predictive errors in both steady and dynamic processes. Part of the reason is that the parameter optimization of the mechanism model can make it close to the true value, which significantly reduces the predictive error, but still cannot be consistent with the true value completely. Another part of the reason may be that there are errors which cannot be eliminated due to the simplification of mechanism model. Because of the differences in the optimization results of parameter $w_{bo}^2{}_{11}^{\left(2\right)}$ that cannot be ignored, there are large predictive errors in dynamic processes of HMPIO for $T_{gs}$, which leads to the dynamic error of $Q_{sl}$ indirectly. Therefore, the ECM is added to compensate the error of HMPIO. The green line in the figures represents the hybrid model with ECM. The hybrid model (HM) proposed in this paper reduces the predictive residual of HMPIO effectively. The results of GDTA in

Table 11. The evaluation indexes of digital twin model for furnace.

Indexes	Output	MMNP	EM-LSTM	EM-GRU	HMPIO	HM
AOP (%)	${\rho }_b$	5.23	0.80	0.63	3.99	0.01
	$T_{gs}$	9.64	0.07	0.07	0.28	0.02
	$P_b$	3.68	0.53	0.38	2.69	0.01
	$O_{cp}$	61.55	0.4	0.36	0.24	0.04
	$Q_{sl}$	164.94	0.13	0.07	0.56	0.05
ACVAR	${\rho }_b$	0.0018	$7.09\times {10}^{-4}$	$8.39\times {10}^{-4}$	0.0011	$7.63\times {10}^{-4}$
	$T_{gs}$	0.0074	0.0012	0.0012	0.0012	0.0013
	$P_b$	0.0014	$4.99\times {10}^{-4}$	$5.94\times {10}^{-4}$	$8.05\times {10}^{-4}$	$5.41\times {10}^{-4}$
	$O_{cp}$	0.1985	$4.22\times {10}^{-5}$	$4.69\times {10}^{-5}$	$7.22\times {10}^{-5}$	$4.51\times {10}^{-5}$
	$Q_{sl}$	0.6221	0.0025	0.0025	0.0028	0.0025
GDTA (%)	global	49.01	0.39	0.30	1.55	0.03
GAVAR	global	0.1662	$9.90\times {10}^{-4}$	0.0010	0.0012	0.0010
ART (s)	global	$3.25\times {10}^{-6}$	0.0048	0.0045	$2.69\times {10}^{-6}$	0.0024
RTVAR	global	$4.51\times {10}^{-9}$	$1.42\times {10}^{-6}$	$1.10\times {10}^{-6}$	$2.01\times {10}^{-9}$	$3.15\times {10}^{-7}$

demonstrate that the HM has the highest accuracy, only one tenth of EM. Due to the ECM, the response time of HM is slower than that of MMNP and HMPIO, but only half of that of EM. Although the improvement of HM accuracy comes at the cost of sacrificing some response time, it can still meet the requirements for the boiler system. Similarly, the predictive results of drum, including $M_{dl}$, ${\rho }_v$, $H_w$ and $H_r$, are shown in Figure 17 and Figure 18. The reason why the curve of MMNP for output $M_{dl}$ has been rising is that the initial value of parameter $w_{vp}^2{}_{14}^{\left(3\right)}$ is inaccurate, resulting in the mass imbalance of liquid in drum. Due to the mutual conversion and calculation between vapor and liquid, the drum is a complex non-linear block, which leads to poor predictive performance of the simplified mechanism model MMNP, even for HMPIO. For outputs $M_{dl}$ and $H_w$, even though HMPIO has better predictive effect than MMNP, it still has obvious predictive errors, especially in the dynamic processes. Meanwhile, as to output $M_{dl}$, the predictive result of EM is also not very satisfactory. Comparing the results of these models, we can see that the hybrid model has relatively excellent predictive performance. From the results of GDTA in

Table 12. The evaluation indexes of digital twin model for drum.

Indexes	Output	MMNP	EM-LSTM	EM-GRU	HMPIO	HM
AOP (%)	$M_{dl}$	35.72	3.07	4.80	11.92	0.48
	${\rho }_v$	3.20	0.27	0.19	0.22	0.09
	$H_w$	9.14	0.42	0.52	2.48	0.62
	$H_r$	2.24	0.24	0.22	0.85	0.18
ACVAR	$M_{dl}$	0.0350	0.0022	0.0029	0.0035	0.0021
	${\rho }_v$	0.0022	0.0016	0.0016	0.0017	0.0017
	$H_w$	0.0080	0.0010	0.0010	0.0038	0.0014
	$H_r$	0.0041	0.0039	0.0038	0.0039	0.0038
GDTA (%)	global	12.58	1.00	1.43	3.87	0.34
GAVAR	global	0.0123	0.0022	0.0023	0.0032	0.0023
ART (s)	global	$4.75\times {10}^{-6}$	0.0044	0.0046	$7.24\times {10}^{-6}$	0.0022
RTVAR	global	$1.45\times {10}^{-9}$	$6.41\times {10}^{-7}$	$8.70\times {10}^{-6}$	$5.19\times {10}^{-9}$	$1.04\times {10}^{-7}$

, we can derive that the accuracy of HM has been improved by 66% and 76% compared to EM-LSTM and EM-GRU, respectively. Finally, we provide the flexibility comparison in

Table 13. The flexibility index for digital twin model.

Question	MMNP	HMPIO	EM	HM
Q1: Can the DTM be used in other use cases without any changes? (YES (Y) or NO (N))	N (−0.8)	N (−0.8)	Y (+0.8)	N (−0.8)
Q2: Does the DTM require specific experiments to develop the DT (S) or can it be developed using real-time or historical data (NS)?	NS (+0.5)	NS (+0.5)	NS (+0.5)	NS (+0.5)
Q3: Can the number of input parameters be increased?	N (−0.6)	Y (+0.6)	Y (+0.6)	Y (+0.6)
Q3.1: Are there any limitations to increasing the number of input parameters? (if Q3 is YES)	/	N (+0.2)	N (+0.2)	N (+0.2)
Q3.3: Does increasing the number of input parameters require modification of the DTM? (if Q3 is YES)	/	Y (−0.2)	N (+0.2)	Y (−0.2)
Q3.4: Does increasing the number of input parameters affect the output parameters? (if Q3 is YES)	/	Y (−0.2)	Y (−0.2)	Y (−0.2)
Q4: Can the number of input parameters be decreased?	N (−0.4)	N (−0.4)	Y (+0.4)	N (−0.4)
Q4.1: Does decreasing the number of input parameters affect the output parameters? (if Q4 is YES)	/	/	Y (−0.2)	/
Q4.2: Does decreasing the number of input parameters require modification of the DTM? (if Q4 is YES)	/	/	N (+0.2)	/
Q5: Can the DTM handle more than one output parameter at the same time? (No need for aggregation of multiple parameters into a single value.)	Y (+0.6)	Y (+0.6)	Y (+0.6)	Y (+0.6)
Q6: Can the number of output parameters be increased?	N (−0.5)	N (−0.5)	Y (+0.5)	N (−0.5)
Q6.1: Are there any limitations on increasing the number of output parameters in the DTM? (if Q6 is YES)	/	/	Y (−0.2)	/
Q7: Can the DTM adapt to new operational conditions?	Y (+0.8)	Y (+0.8)	Y (+0.8)	Y (+0.8)
Q7.1: Is the adaptation of DTM performed automatically(A)/manually(M)? (if Q7 is YES)	A (+0.3)	A (+0.3)	A (+0.3)	A (+0.3)
Q7.2: Does the adaptation of the DTM require additional data? (if Q7 is YES)	Y (−0.3)	Y (−0.3)	Y (−0.3)	Y (−0.3)
Q8: Can the DTM be executed using a conventional hardware computer?	Y (+0.3)	Y (+0.3)	Y (+0.3)	Y (+0.3)
Q9: Does the use of the DTM require special training?	N (+0.3)	N (+0.3)	N (+0.3)	N (+0.3)
Total score	5	6	9.6	6

according to reference [59]. From the total score, EM has the highest flexibility, followed by HMPIO and HM, and MMNP has the lowest. HM, MMNP and HMPIO are not as flexible as EM due to their mechanism model limiting the number and relationship of input and output parameters. Although the interpretability of HM comes at the cost of sacrificing some flexibility, compared to EM, HM improves the accuracy and response time of digital twin models effectively. In conclusion, the hybrid model established in this paper has relatively excellent predictive performance, which ensures the high fidelity of digital twin model for boiler system. However, the datasets used for training and verifying the model in this paper is limited and may not cover all working conditions. Therefore, the generalization performance of hybrid model, especially ECM, needs to be further studied in future work.

5. Conclusions and Future Work

In this paper, an RNN-based hybrid modeling method for digital twin of boiler system is established to improve the predictive accuracy of mechanism model. The mechanism model of main blocks in boiler system is described through RNN to facilitate training and updating parameters. By adding neurons to mechanism model, the mapping between time-varying parameters and state variables is constructed. In order to compensate the predictive residual of the mechanism model, the LSTM-based error compensation model is established. Moreover, the update architecture and training algorithm for hybrid model are built to realize the iterative optimization of model parameters. Finally, the experimental results show that the hybrid model has better predictive performance.

However, there are some shortcomings that need to be further studied. First, for a few parameters, the optimization results are not satisfactory. In future work, we will explore corresponding solutions. Second, in the training algorithm, the calculation of some gradient terms is complex, which affects the optimization results of parameters. In the future, we will try to develop an effective approach to optimize the calculation of gradient. Finally, the datasets used for training and verifying the model are collected from a 600 MW power plant simulator. In the future, we hope to further validate the hybrid model based on the actual power plant datasets.