Research on the Prediction of Load Regulation Capacity for Supercritical Thermal Power Unit

Yong He¹, Zhuo Huang^2,3, Yaohua Tang⁴, Weimin Guo⁴, Yuanli Cai^2,3*, Haonan Jiang^2,3

¹State Grid Henan Electric Power Research Institute, Zhengzhou, China.
²Xiamen Institute of Technology, Xiamen, China.
³Faculty of Electronic and Information Engineering, Xi’an Jiaotong University, Xi’an, China.
⁴Rundian Energy Science and Technology Co., Ltd., Zhengzhou, China.

DOI: 10.4236/jpee.2020.88002                

Abstract

Both the modeling and the load regulation capacity prediction of a supercritical power plant are investigated in this paper. Firstly, an indirect identification method based on subspace identification method is proposed. The obtained identification model is verified by the actual operation data and the dynamic characteristics of the system are well reproduced. Secondly, the model is used to predict the load regulation capacity of thermal power unit. The power, main steam pressure, main steam temperature and other parameters are simulated respectively when the unit load is going up and down. Under the actual constraints, the load regulation capacity of thermal power unit can be predicted quickly.

Keywords

Supercritical Thermal Power Unit, Subspace Identification, Indirect Identification, Load Regulation

1. Introduction

Due to the different distribution of power production, the power supply and consumption in a region cannot always keep balance [1]. In order to reduce the impact on industrial production and life of human beings, dynamic dispatching of power grid is needed. Before dispatching center issues dispatching instructions, it is necessary to know the capability of the units in the whole power grid to increase power in a short time (generally 1 - 3 minutes). Therefore, one unit needs to predict how much power it can regulate in the next 1 - 3 minutes according to its current state. Normally, a load regulation capacity prediction model is set up and uploaded to the dispatching center as the basis for rapid dispatching. In order to ensure safe operation of the whole power grid system in a region and improve its resilience in case of accidents, it is necessary to accurately master the power generation load regulation capacity to implement the refined unit optimal scheduling.

Supercritical unit is a large and complex system with strong coupling, nonlinearity and multivariable characteristics. Although there are many works on the dynamic characteristics of supercritical units [2] - [9], the available models are mainly suitable for power plant simulators, which are usually nonlinear and too complex to be directly used to the control system design and load prediction. In the past, there is little related research of load regulation capacity as the reference of dispatching center. In this paper, we propose to combine the load regulation of the unit with the optimization of the controller parameters of the coordinated control system. This combination can not only ensure the stable operation of the unit, but also guide the adjustment of the parameters of the existing coordinated control system of the power plant unit. Thus the quality of power generation can be improved and the maximum economic benefits can be produced.

This paper deals with the modelling of supercritical thermal power system and the regulation of unit load. The rest of this paper is arranged as follows. Based on the principle of subspace identification, an indirect model identification method for supercritical power plants is proposed in Section 2. With the designed model predictive controller, Section 3 investigates the load regulation of a supercritical thermal power unit. Some concluding remarks are drawn in Section 4.

2. Model Identification of Supercritical Thermal Power System

In order to predict the load capacity of supercritical thermal power unit, the mathematical model of thermal power unit should be established firstly. Due to the complex physical and chemical processes involved in thermal power units, it is very difficult to establish a complete mechanism model. Subspace identification is an advanced system identification method, which has been successfully applied in many fields. Compared with the traditional system identification methods, subspace identification method has many advantages, including no parameterization, no iterative optimization, and only QR decomposition and other simple linear algebra tools are involved. In addition, the state space model of the system is directly obtained with the subspace identification method, which is convenient for the state estimation and control algorithm design of multivariable systems. Based on the actual input and output data of a thermal power unit, this paper uses subspace identification method to explore the system modeling of thermal power unit.

2.1. Basic Principle of Subspace Identification Method

There are three basic algorithms for subspace identification, namely MOESP (Multivariable out-put Error State space), N4SID (Numerical algorithm for Subspace State space System Identification) and CVA (Canonical Variate Analysis), and they have different features. In this paper, N4SID [10] [11] open-loop subspace identification method is used to identify the supercritical thermal power unit system model. Subspace identification method adopts the form of state space to describe the system to be identified. Hankle matrices are constructed and calculated according to the input and output data. The core idea is to decompose the system into a deterministic part and a stochastic part by using the nature of oblique projection so as to eliminate the noise [12].

Consider the following linear time-invariant system:

$x\left(k+1\right)=Ax\left(k\right)+Bu\left(k\right)+w\left(k\right)$ (1)

$y\left(k\right)=Cx\left(k\right)+Du\left(k\right)+v\left(k\right)$ (2)

where $x\left(k\right)\in R^n$ is the system state at time k, $u\left(k\right)\in R^p$, $y\left(k\right)\in R^m$ are the input and output signals at time k respectively. $w\left(k\right)$ and $v\left(k\right)$ indicate the process and measurement noises respectively. Without the noise terms, a deterministic system is described.

From (1) and (2) while considering the deterministic case, one can easily get

$\left[\begin{array}{c}y\left(k\right)\\ y\left(k+1\right)\\ ⋮\\ y\left(k+i\right)\end{array}\right]=\left[\begin{array}{c}C\\ CA\\ ⋮\\ C{A}^i\end{array}\right]x\left(k\right)+\left[\begin{array}{cccc}D& 0& 0& 0\\ CB& D& 0& 0\\ ⋮& ⋮& \ddots & 0\\ C{A}^{i-1}B& C{A}^{i-2}B& \cdots & D\end{array}\right]\left[\begin{array}{c}u\left(k\right)\\ u\left(k+1\right)\\ ⋮\\ u\left(k+i\right)\end{array}\right]$ (3)

According to the standard expression of subspace identification, Equation (3) can be decomposed into past and future parts, namely

$X_f=A^i{X}_p+{\Delta }_i{U}_p+{\Delta }_i^s{E}_p$ (4)

$Y_p={\Gamma }_i{X}_p+H_i^d{U}_p+H_i^s{E}_p$ (5)

$Y_f={\Gamma }_i{X}_f+H_i^d{U}_f+H_i^s{E}_f$ (6)

where

$U_p=\left(\begin{array}{cccc}u\left(0\right)& u\left(1\right)& \cdots & u\left(j-1\right)\\ u\left(1\right)& u\left(2\right)& \cdots & u\left(j\right)\\ ⋮& ⋮& \ddots & ⋮\\ u\left(i-1\right)& u\left(i\right)& \cdots & u\left(i+j-2\right)\end{array}\right)\in R^{pi\times j}$

$U_f=\left(\begin{array}{cccc}u\left(i\right)& u\left(i+1\right)& \cdots & u\left(i+j-1\right)\\ u\left(i+1\right)& u\left(i+2\right)& \cdots & u\left(i+j\right)\\ ⋮& ⋮& \ddots & ⋮\\ u\left(2i-1\right)& u\left(2i\right)& \cdots & u\left(2i+j-2\right)\end{array}\right)\in R^{pi\times j}$

$Y_p$ and $Y_f$ are constructed on the same way as $U_P$ and $Y_f$.

The observability matrix ${\Gamma }_i$, the lower triangle Toeplitz matrix $H_i^d$ and other matrices are defined as follows:

${\Gamma }_i=\left[\begin{array}{c}C\\ CA\\ C{A}^2\\ ⋮\\ C{A}^i\end{array}\right],H_i^d=\left[\begin{array}{cccc}D& & & \\ CB& D& & \\ ⋮& \ddots & \ddots & \\ C{A}^{i-2}B& C{A}^{i-3}B& \cdots & D\end{array}\right],H_i^s=\left[\begin{array}{cccc}0& & & \\ C& 0& & \\ ⋮& \ddots & \ddots & \\ C{A}^{i-2}& C{A}^{i-3}& \cdots & 0\end{array}\right]$

${\Delta }_i=\left[\begin{array}{cccc}{A}^{i-1}B& A^{i-2}B& \cdots & B\end{array}\right],{\Delta }_i^s=\left[\begin{array}{cccc}{A}^{i-1}& A^{i-2}& \cdots & I\end{array}\right]$

First, we project $Y_f$ onto the orthogonal complement subspace $U_f^{\perp }$ of $U_f$, and obtain

$Y_f/U_f^{\perp }={\Gamma }_i{X}_f/U_f^{\perp }+H_i^d{U}_f/U_f^{\perp }+H_i^s{E}_f/U_f^{\perp }$ (7)

Multiply the weighting matrices $W_1$ and $W_2$ on both sides of Equation (7), and through the correlation and orthogonality, the above formula can be simplified as

$W_1{Y}_f/U_f^{\perp }{W}_2=W_1{\Gamma }_i{X}_f/U_f^{\perp }{W}_2+W_1{H}_i^s{E}_f{W}_2$ (8)

The weighting matrices $W_1$ and $W_2$ should meet the following conditions:

1) $rank\left(W_1{\Gamma }_i\right)=rank\left({\Gamma }_i\right)$ ;

2) $rank\left(X_f/U_f^{\perp }{W}_2\right)=rank\left(X_f\right)$ ;

3) $W_1{H}_i^s{E}_f{W}_2=0$.

The third condition requires that $W_2$ is independent of $E_f$, so Equation (8) may be further simplified as

$\Phi =W_1{Y}_f/U_f^{\perp }{W}_2=W_1{\Gamma }_i{X}_f/U_f^{\perp }{W}_2$ (9)

Through SVD decomposition, Equation (9) can be denoted as

$\Phi =\left(\begin{array}{cc}{U}_1& U_2\end{array}\right)\left(\begin{array}{cc}{S}_1& 0\\ 0& 0\end{array}\right)\left(\begin{array}{c}{V}_1^{\text{T}}\\ V_2^{\text{T}}\end{array}\right)$ (10)

The order of the system is determined according to the order of matrix $S_1$.

Using MATLAB notation, the matrix $C$ can be directly derived from ${\Gamma }_i$ as follows:

$C={\Gamma }_i\left(1:m,:\right)$ (11)

According to the following property

$\underset{\_}{{\Gamma }_i}\cdot A=\overline{{\Gamma }_i}$ (12)

where $\underset{\_}{{\Gamma }_i}$ is the matrix ${\Gamma }_i$ deleting the last i rows, and $\overline{{\Gamma }_i}$ is the matrix ${\Gamma }_i$ deleting the first i rows. By left multiplying ${\underset{\_}{{\Gamma }_i}}^{†}$ of Equation (12), we can obtain matrix $A$ as follows:

$A={\underset{\_}{{\Gamma }_i}}^{†}\overline{{\Gamma }_i}$ (13)

If both sides of Equation (6) multiply ${\Gamma }_i^{†}$ by left and $U_f^{†}$ by right, we can get

${\Gamma }_i^{†}{Y}_f{U}_f^{†}={\Gamma }_i^{†}{H}_i^d$ (14)

Let the left side of the above equation be denoted as $P$ and ${\Gamma }_i^{†}$ be denoted as $L$, then

$\left(\begin{array}{cccc}{P}_1& P_2& \cdots & P_i\end{array}\right)=\left({\begin{array}{cccc}{L}_1& L_2& \cdots & L\end{array}}_i\right)\times \left[\begin{array}{cccc}D& 0& 0& 0\\ CB& D& 0& 0\\ ⋮& ⋮& \ddots & 0\\ C{A}^{i-2}B& C{A}^{i-3}B& \cdots & D\end{array}\right]$ (15)

Equation (15) constitutes a set of overdetermined linear equations, which can be solved by the least square method to obtain matrices $B$ and $D$.

2.2. Indirect Identification Method

The main controlled variables of a supercritical thermal power unit are unit load, main steam pressure and main steam temperature. A large number of simulation experiments show that based on the subspace identification method discussed in the previous subsection, the system load and main steam pressure can be identified satisfactorily. However, if the input and output data are used to directly identify the main steam temperature, the effect is poor. A typical set of direct identification results is shown in Figure 1. The data are normalized and the sampling period is 5 s, which applies to all the figures in this paper.

To overcome the above problem in direct identification, we have made a detailed analysis on the process and propose the following indirect identification method. Firstly, the mathematical model between the middle point temperature and the main control inputs is identified by using the fuel feed, the water supply and the opening degree of the steam turbine valve. Then the main steam temperature of the system is obtained by the steam flow at the middle point passing through the multi-stage superheater.

The main steam temperature model obtained by the indirect identification method can meet the requirements of system control. Figure 2 shows the main

Figure 1. Direct identification results of main steam temperature.

steam temperature identifying results obtained by the proposed indirect identification method. All-day operation data of a power plant are used in the system identification.

Taking a thermal power plant in Henan Province as an example, the basic process and results based on the indirect identification method are given below. The model structure is shown in Figure 3.

Firstly, according to the actual operation data of the system, the subsystem from the turbine valve opening, total fuel quantity, and total water supply to the system load and main steam pressure are identified. It is a module with three inputs and two outputs, as shown in Figure 3. The related system matrices are as follows:

$\begin{array}{l}{A}_{\text{1}}{}_{\text{2}}=\left[\begin{array}{cc}0.7773& -0.1209\\ -0.1516& 0.8871\end{array}\right]\\ B_{\text{1}}{}_{\text{2}}=\left[\begin{array}{c}\begin{array}{ccc}-0.001713& -0.00007509& -0.001384\end{array}\\ \begin{array}{ccc}-0.001267& -0.0005914& -0.001289\end{array}\end{array}\right]\\ C_{\text{1}}{}_{\text{2}}=\left[\begin{array}{cc}29.96& -72.14\\ 40.14& -77.42\end{array}\right]\\ D_{\text{1}}{}_{\text{2}}=\left[\begin{array}{c}\begin{array}{ccc}0& 0& 0\end{array}\\ \begin{array}{ccc}0& 0& 0\end{array}\end{array}\right]\end{array}$

Figure 2. Indirect identification results of main steam temperature.

Figure 3. Block diagram of indirect identification model.

Secondly, the subsystem from the turbine valve opening, total fuel feed, total water supply to the intermediate point temperature is identified. It is a module with three inputs and one output, as shown in Figure 3. The obtained system matrices are as follows:

$\begin{array}{l}{A}_{\text{03}}=\left[\begin{array}{cc}0.9619& -0.02362\\ -0.06153& 0.8526\end{array}\right]\\ B_{\text{03}}=\left[\begin{array}{ccc}-0.0009411& -0.0007184& -0.0004354\\ -0.004808& -0.002631& -0.002633\end{array}\right]\\ C_{\text{03}}=\left[\begin{array}{cc}16.21& -12.66\end{array}\right]\\ D_{\text{03}}=\left[\begin{array}{ccc}0& 0& 0\end{array}\right]\end{array}$

Thirdly, the subsystem from the intermediate point temperature, air-coal ratio, desuperheating water flow to the main steam temperature is identified. It is a module with three inputs and one output in Figure 3. The identified system matrices are as follows:

$\begin{array}{l}{A}_{\text{3}}=\left[\begin{array}{cc}0.8611& -0.02323\\ 0.1741& 0.9902\end{array}\right]\\ B_{\text{3}}=\left[\begin{array}{c}\begin{array}{ccc}-0.0003134& -0.0006446& -0.0004444\end{array}\\ \begin{array}{ccc}-0.0005305& -0.000952& -0.0006654\end{array}\end{array}\right]\\ C_{\text{3}}=\left[\begin{array}{cc}16.41& 20.84\end{array}\right]\\ D_{\text{3}}=\left[\begin{array}{ccc}0& 0& 0\end{array}\right]\end{array}$

The above separated models are repeatedly verified by using different operation data, which include almost all the characteristics of the normal operation. Based on the above results, the total system model of the thermal power unit is obtained as follows:

$A=\left[\begin{array}{cccccc}0.9862& -0.0064& 0& 0& 0& 0\\ -0.01069& 0.9911& 0& 0& 0& 0\\ 0& 0& 0.9996& -0.045& 0& 0\\ 0& 0& 0.02414& 0.919& 0& 0\\ 0& 0& -\text{0}\text{.000602}& \text{1}\text{.02}e-\text{05}& 0.9955& -0.05659\\ 0& 0& \text{0}\text{.00169}& -\text{2}\text{.87}e-\text{05}& 0.05911& 0.8846\end{array}\right]$

$B=\left[\begin{array}{ccccc}0& 0.0001118& 5.615e-05& 4.623e-06& 0\\ 0& 0.0001065& 4.608e-05& \text{2}\text{.23}e-\text{06}& 0\\ 0& -0.0001031& -0.0001032& -\text{0}\text{.000118}& 0\\ 0& -\text{0}\text{.000166}& -0.0004942& \text{8}\text{.7}e-\text{05}& 0\\ 0.0001053& 0& 0& 0& -\text{9}\text{.55}e-\text{06}\\ 0.0002139& 0& 0& 0& \text{4}\text{.85}e-\text{05}\end{array}\right]$

$C=\left[\begin{array}{cccccc}106.5& 3.506& 0& 0& 0& 0\\ 70.69& 68.54& 0& 0& 0& 0\\ 0& 0& 0& 0& 31.04& -0.6756\end{array}\right]$

$D=\left[\begin{array}{ccccc}0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\end{array}\right]$

It is seen that the overall identification model of the unit has five inputs and three outputs. Based on the above model, a lot of experiments are carried out. The actual operation data of the unit in different time periods are used to verify the model and very satisfactory results are obtained. Figure 4 shows a set of actual output and the model output.

It is shown that the load deviation is about −10 - 10 MW, the main steam pressure deviation is about −0.2 - 0.2 MPa, and the main steam temperature deviation is about −5˚C - 5˚C. The accuracy completely meets the requirements of the model prediction control method. The effectiveness of the proposed indirect identification is well verified.

3. Prediction of Load Regulation Capacity

The load regulation capacity of a thermal power unit is directly related to the controller used in the system. The thermal power plant unit in our research adopts the model predictive control, which is also realized in the prediction of the actual load regulation capacity of the system. Model predictive control is the most promising advanced control method in industrial process control. According to the specific performance index, the actual control problem is transformed into the solution of a constrained optimization problem [13]. The elementary idea is to solve the optimization problem in real time according to the current sampling results, and update the control effort at the current time.

3.1. Load Regulation Strategy

In the actual operation process of the system, there are two control strategies to increase or decrease the generation load. One is the “turbine following boiler” strategy. When the generating load changes, the boiler side of the unit will adjust the corresponding control inputs according to the generating demand. By increasing and decreasing the boiler energy storage, the purpose of matching the boiler capacity with the demand of the generating unit can be achieved. The other one is the “boiler following turbine” strategy. The generator unit directly adjust the opening of the steam turbine valve according to the actual power generation demand, which ensures that the generating capacity of the unit will increase or decrease rapidly. The control goal of the boiler side is to maintain the stability of the main steam temperature and main steam pressure in the boiler.

We have adopted the “turbine following boiler” strategy. By adjusting the actual energy storage at the boiler side to adjust the future load change of the system, i.e. using predictive control method, can quickly adjust the input of relevant fuel and auxiliary raw materials. Also, the opening of the turbine valve can be gradually adjusted. Thus the system energy output changes rapidly, and at the same time, the stability of the main steam temperature and the slow change of the main steam pressure in the system can be maintained.

3.2. Load Rise Forecast

In the case that the thermal power unit needs to increase the power rapidly, set

(a)(b)(c)

Figure 4. Identified model output vs actual plant output. (a) System load; (b) Main steam pressure; (c) Main steam temperature.

the load setting value of the power generation as large as possible, so as to show the ultimate capacity of the system to increase the load without violating the relevant constraints. Because the upper and lower bounds of each variable and the upper and lower bounds of the rate of change are limited, especially the constraints on the main steam temperature (±2˚C) and the main steam pressure (±10 MPa), the controller ensures that the system can maintain stable operation state when the load is increased rapidly.

In the following simulation, the unit controller adopts the constrained model predictive control method based on the identification model. The main control signals during load raising of the unit are shown in Figure 5, including valve opening (%), coal consumption (t/h), total water supply (t/h), air-coal ratio and desuperheating water (t/h), respectively.

Figure 5 shows that in the controlled input constraint range, in order to make full use of the energy storage and accelerate the initial load response speed of the unit, the valve opening gradually increases. The valve opening is the main means of rapid regulation of the system. Its steady-state value is maintained between 85% and 90%, and the valve opening returns to the original state after achieving the purpose of rapid load regulation. Coal and water are the main control inputs of the generating unit. When the unit is loaded, the coal feed and water supply increase rapidly. The coal feed and water supply determine the internal variable state when the load of the unit is finally stable. The air-coal ratio is an important parameter affecting the unit load. In order to ensure full combustion and less heat dissipation of coal, the air-coal ratio is set as an empirical constant. Desuperheating water is the last process of regulating the main steam temperature, and it is an auxiliary means. Considering its regulating effect on the main steam temperature, there is no strict restriction here.

Figure 5. System inputs at load rise.

The system outputs at load rise period are shown in Figure 6, including unit load (MW), main steam pressure (MPa) and main steam temperature (˚C).

It can be seen from Figure 6 that when the system load needs to increase rapidly, the main steam temperature and main steam pressure will change accordingly. It is not difficult to find that coal and water have a certain lag effect on the unit load. The variation range of main steam temperature and main steam pressure meet the given constraints. Under various constraint conditions, the maximum load change is a gradual rise processing, which reflects the load raising capacity of the unit.

3.3. Load Reduction Forecast

In case of rapid power reduction, the setting value of the generation load is set to a downward step as large as possible. In order to facilitate the comparison with the results of load rise, the same constraints and processing methods are adopted for each variable during load decrease. The simulation results of the control inputs are shown in Figure 7. It can be seen that the control inputs vary within the constraint range. In order to reduce the unit load rapidly, the valve opening decreases rapidly and then gradually returns to the original state. The coal feed and water supply decrease rapidly, the air-coal ratio remains unchanged, and the desuperheating water has no strict restriction.

The controlled variables in load reduction are shown in Figure 8. When the system load drops rapidly, the change of main steam temperature and main steam pressure is kept within certain range. On the premise of satisfying the constraints of each input, the maximum load reduction is a gradual decline curve, which reflects the load reduction capacity of the unit.

Figure 6. System outputs at load rise.

Figure 7. System input at load reduction.

Figure 8. System output at load reduction.

4. Conclusion

Based on the principle of subspace identification, an indirect identification method for supercritical units is proposed in this paper. Compared with the actual operation data, the accuracy of the obtained identification model is good enough and the dynamic characteristics of the system are well reproduced. Based on the model obtained by subspace identification and the actual model predictive control algorithm, the short-term prediction of the system load regulation capacity is realized while keeping the key parameters of the unit meet the relevant constraints. A large number of experimental data show that the established identification model and load regulation algorithm are correct and effective.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	China Electricity Power Committee (2015) China’s Electric Power Industry Annual Development Report. China Market Press, Beijing.
[2]	Zhang, C.X. (1987) Boiler Dynamic Characteristics and Its Mathematical Model. China Water & Power Press, Beijing.
[3]	Ponce, I.U., Bentsman, J., Orlov, Y.and Aguilar, L.T. (2015) Generic Nonsmooth Output Synthesis: Application to a Coal-Fired Boiler/Turbine Unit with Actuator Dead Zone. IEEE Transactions on Control Systems Technology, 23, 2117-2128. https://doi.org/10.1109/TCST.2015.2399672
[4]	Ghabraei, S., Moradi, H. and Vossoughi, G. (2015) Multivariable Robust Adaptive Sliding Mode Control of an Industrial Boiler-Turbine in the Presence of Modeling Imprecisions and External Disturbances: A Comparison with Type-I Servo Controller. ISA Transactions, 58, 398-408. https://doi.org/10.1016/j.isatra.2015.04.010
[5]	Wu, X., Shen, J., Li, Y. and Li, K.Y. (2014) Hierarchical Optimization of Boiler-Turbine Unit Using Fuzzy Stable Model Predictive Control. Control Engineering Practice, 30, 112-123. https://doi.org/10.1016/j.conengprac.2014.03.004
[6]	Kong, X., Liu, X. and Lee, K.Y. (2015) Nonlinear Multivariable Hierarchical Model Predictive Control for Boiler-Turbine System. Energy, 93, 309-322. https://doi.org/10.1016/j.energy.2015.09.030
[7]	Ma, L., Lee, K.Y. and Wang, Z. (2016) Intelligent Coordinated Controller Design for a 600MW Supercritical Boiler Unit Based on Expanded-Structure Neural Network Inverse Models. Control Engineering Practice, 53, 194-201. https://doi.org/10.1016/j.conengprac.2015.09.002
[8]	Liu, X.J., Kong, X.B., Hou, G.L. and Wang, J.H. (2013) Modeling of a 1000MW Power Plant Ultra-Supercritical Boiler System Using Fuzzy-Neural Network Methods. Energy Conversion and Management, 65, 518-527. https://doi.org/10.1016/j.enconman.2012.07.028
[9]	Wang, G., Yan, W., Chen, S., Zhang, X. and Shao, H. (2014) Multi-Model Predictive Control of Ultra-Supercritical Coal-Fired Power Unit. Chinese Journal of Chemical Engineering, 22, 782-787. https://doi.org/10.1016/j.cjche.2014.05.005
[10]	Li, H. and Cai, Y. (2016) N4SID Based Subspace Identification Method for Normal Temperature Continuous Transonic Wind Tunnel System. Proceedings of the 28th Chinese Control and Decision Conference, Yinchuan, 28-30 May 2016, 1394-1398. https://doi.org/10.1109/CCDC.2016.7531202
[11]	Van Overschee, P. and De Moor, B. (1994) N4SID: Subspace algorithms for the Identification of Combined Deterministic-Stochastic Systems. Automatica, 30, 75-93. https://doi.org/10.1016/0005-1098(94)90230-5
[12]	Verdult, V. and Verhaegen, M. (2002) Subspace Identification of Multivariable Linear Parameter-Varying Systems. Automatica, 38, 805-814. https://doi.org/10.1016/S0005-1098(01)00268-0
[13]	Qian, J., Zhao, J. and Xu, Z. (2007) Predictive Control. Chemical Industry Press, Beijing.