Modeling and Finite-Horizon MPC for a Boiler-Turbine System Using Minimal Realization State-Space Model

Abstract

This paper aims to address a finite-horizon model predictive control (MPC) for non-linear drum-type boiler-turbine system using a system-identification method. Considering that the strong state coupling of a non-linear mechanism model, the subspace identification method is first utilized to obtain a linear state-space model, and transformed into an input–output model. By taking the inputs and outputs of the input–output model as system states, an augmented non-minimal state-space (NMSS) model of state measurable is constructed. In order to reduce the computation burden, the augmented NMSS model is further transformed into a canonical formulation by adopting a Kalman decomposition. Based on the minimal realization state-space model, the MPC controller is parameterized as a finite-horizon optimization problem. Finally, simulations are performed and evaluated the performance of the proposed method, and the simulation results show that: the linear model approximate the non-linear system accurately; the proposed MPC method can achieve a satisfactory stable control performance; and the computation time 18.388 s for the overall optimization problem also illustrates the real-time performance effectively.

1. Introduction

The thermal power unit, which mainly converts the chemical fuel (e.g., coal, oil, or gas) into electrical energy, is one of the most important energy conversion devices. For a thermal power unit, the fundamental requirement is to guarantee the pressure and water in the drum within the allowable range, while meeting the load demand of electric power. However, the complex dynamic characteristics of the boiler-turbine result in many difficulties for system modeling and control strategy design [1,2,3].

In practice, the operation of thermal power unit combats many challenges, such as severe non-linearity and multivariable coupling, which are the potential causes for the deteriorating performance of stability of thermal power unit. Over the past several decades, in order to eliminate or decrease the effects caused by these problems, plenty of studies have been conducted in industrial and academia; e.g., [4,5,6]. In [4], the authors utilized the fuzzy method to study the modeling problem of a drum-type non-linear boiler-turbine system. The author in [7] investigated the coordinated control problem for boiler-turbine system of coal-fired power plant based on the fuzzy control method. In [8], a multivariable non-linear exponential ARX model was utilized to characterize the non-linear dynamic features of thermal power plants over the whole operating range. The works of [9,10] consider the modeling problem of a thermal power unit which is characterized as a multi-input multi-output radial basis function-based auto-regressive model with exogenous inputs. In fact, these literature show the drawbacks of controlling boiler-turbine system: linearization is made based on the assumption that part of non-linearity can be ignored, complex controller design is not easy to implement, and the simple control methods cannot achieve accurate control.

Model predictive control (MPC), as an advanced control technique proposed by industry circle, has been recognized for an effective and promising method being applied in various fields, ranging from petrochemicals, refining, textiles, and autonomous vehicles (e.g., [11,12,13,14,15,16,17]). In [18], the author designed a robust MPC controller for an automatic voltage regulator against uncertainties. Ref. [19] proposes a heuristic algorithm of MPC to tune the parameters for the robotic manipulator. MPC optimizes an optimal control sequence based on the prediction of future outputs within prediction horizon, and only implements the first control input. The other strength of MPC is that it can directly deal with various hard or soft constraints. Motivated by this outstanding advantage, many researchers have utilized MPC to handle the control problems of thermal power unit (e.g., [6,20,21]). In [6], the authors presented two different step-response models to design dynamic matrix control strategy for a boiler-turbine system. The results showed that, for controller design of the boiler-turbine system, the step-response model based on test data is more effective than a linearized model. In [20], an MPC method merged with the genetic algorithm was introduced to tackle the boiler-turbine control problem. Recently, a hierarchical control framework was proposed in [21] to design an optimal control strategy for a boiler-turbine system, described by the Takagi–Sugeno fuzzy model. The work of [22] considers a fuzzy MPC method for the boiler-turbine system subject to disturbances and uncertainties. Apart from focusing on the engineering problems, some studies also concern the economic and environmental issues of thermal power unit operation. In [2], the authors investigated a hierarchical MPC of a boiler-turbine system which considers the plant-wide economic process optimization and regulatory process control simultaneously. Ref. [23] studies an optimization control method that utilizes MPC for the set-point optimization to improve the safety and economic performance for boiler-turbine system. In [24], from the perspective of scheduling and optimization, the authors presented a zone tracking MPC method, which considers system economic performance during the transient state while always prioritizing unit load demand tracking. Ref. [25] considers the set-point tracking problem for a boiler-turbine unit represented by a non-linear model. In order to obtain an efficient control scheme, the state-space model is on-line linearized at the current operating point, and utilized for prediction and control policy optimization. Unfortunately, the linearized model approach usually assumes that the non-linear part can be ignored, which may cause model mismatch, and results in inaccurate control and performance deterioration. More recently, considering the complexity on the modeling for boiler-turbine system, since artificial intelligence techniques show potential solutions for several complex problems, the academic circle attempt to introduce these techniques, especially in non-linear system modeling [10,26,27] and controller design [28,29,30,31].

In the control literature, there is a abundance of model formulations being utilized in various research directions [32,33,34,35,36], where the state-space model is still the most widely studied models for optimal control, robust control and MPC for its advantage at simplicity implementation of linear system. However, the shortcoming of the state-space model is that it requires states that are measurable when applying the state feedback control, or designs a state estimator when applying the output feedback control. In contrast, the non-minimal state-space (NMSS) model can effectively realize state feedback control by augmenting measured outputs, inputs, and their past values, while avoiding the design of state estimator. Few papers related on NMSS model can be seen in the literature (e.g., [37,38,39,40]). In [39], an improved NMSS-based MPC for multivariable system, which employs a non-zero-pole decoupling approach, was proposed. Based on an augmented NMSS model, the authors in [38] designed a fractional-order MPC strategy for the temperature control model in the industrial heating furnace. For the operation of a boiler-turbine system, an extensive of measurement data can be easily stored and obtained, thus the subspace identification modeling method has been applied for the modeling, multivariable controller design and optimization. In [41], the authors developed a data-driven modeling strategy and designed an MPC controller for boiler-turbine system by dividing the system into several local regions according to the operation range based on the subspace identification and multi-model methods. The subspace identification method is able to estimate the state-space model without non-linear iterative calculation and model parameterization. Despite the identified state-space model widely utilized, their states of these models are meaningless. In order to handle this problem, the idea of augmenting the measured outputs, inputs and their past values of the identified state-space model to construct a NNSS model is a potential approach so that each state variable is measurable or exactly known. Usually, when designing a MPC controller based on an NMSS model directly, the computation burden of MPC optimization problem is large. In fact, utilizing the concept of minimal realization formulation is an alternative approach, but it is rarely reported in the literature. On the other hand, many MPC implementations of NMSS model use heuristic method without stability consideration, this may leads to sub-optimal performance. It is, therefore, necessary to consider optimality to guarantee the control performance, this formulated the main standpoint of this paper.

Motivated by the above discussion, considering the inherent non-linear dynamics of boiler-turbine system and its effect for designing an effective controller, the work hopes to present a clear and intuitive path on how the identification method and MPC can be applied. Consider the influence of non-linear dynamics, we first linearize the non-linear system using an identification method. By adopting the idea of NMSS model and its minimal realization formulation, the state information will be always known and observers are not needed. The main contributions are in the identified model utilized for infinite-horizon MPC (including model minimal realization, terminal control law design, and nominal stability taken into account) and in the unique control structure to achieve good tracking ability, which can be highlighted as follows:

• The numerical algorithm for subspace state-space identification (N4SID) is utilized for a drum-type boiler-turbine system to obtain a linearized state-space model. By taking the inputs and outputs of the state-space model from the subspace identification method as system states, an augmented NMSS model with state measurable is constructed to avoid a state observer;
• The augmented NMSS model is transformed into a canonical formulation by adopting a Kalman decomposition in order to reduce the computation burden of controller parameter optimization;
• Based on the minimal realization state-space model, an MPC controller is transformed into solving a finite-horizon optimization problem, where the cost function is composed by a finite horizon cost and terminal cost. The nominal stability of finite-horizon MPC are also guaranteed for the resulting model.

In the presented control structure, we show a development and combination of MPC and minimal realization model based on previous works [23,42] for a boiler-turbine system. Compared to the literature in [2,23,25,42], the primary novelty and difference of this paper reside in constructing a state measurable, minimal-realization state-space model for the non-linear boiler-turbine system by using identification method to tackle the process non-linearity that the maintains the simplicity of the framework for linear MPC using state-space model, so that the designed MPC controller also has good tracking ability and implementation efficiency.

This rest of this paper is organized as follows. Section 2 shows the facility description of a thermal power unit and its physical process. Section 3 details the identification process of NMSS model. Section 4 presents the finite-horizon MPC method based on NMSS model. Section 5 verifies the effectiveness of the proposed method. Finally, Section 6 summarizes this paper.

2. The Overview for Drum-Type Boiler-Turbine

2.1. Description of the Process

Figure 1 shows a simplified schematic of a drum-type boiler-turbine system that usually includes combustion system, steam-water system, control system, electrical system, condensing system, and so on. The combustion system and steam-water system generate the steam to a certain extent high temperature and pressure. The electric system completes the energy conversion, i.e., realizing the transformation from thermal energy, mechanical energy to electric energy. The control system ensures the safe and economic operation of each subsystem.

The operation process of boiler-turbine system is briefly introduced. Initially, the treated pulverized coal from a coal mill is conveyed to the boiler by belt conveying technology. Pulverized coal is burned in the furnace to heat the boiler in order to convert the water into steam. Steam enters the drum and continues to be heated, and part of steam the high pressure cylinder after primary heating. The drum has the energy store function. When the load changes, it has a certain buffer effect on the imbalance between evaporation and water supply, and rate of change on steam pressure. Steam is reheated in order to improve the thermal efficiency and enters the medium pressure cylinder. Turbine generators are driven by steam from medium pressure cylinders.

2.2. Boiler-Turbine Non-Linear Dynamic Model

In this paper, the classical non-linear model of a 160 MW boiler-turbine generator unit firstly investigated in [43] is considered

$$\left\{\begin{array}{c}{\dot{x}}_1=-0.0018u_2{x}_1^{9/8}+0.9u_1-0.15u_3\hfill \\ {\dot{x}}_2=\left(0.073u_2-0.016\right)x_1^{9/8}-0.1x_2\hfill \\ {\dot{x}}_3=\left(141u_3-\left(1.1u_2-0.19\right)x_1\right)\hfill \\ y_1=x_1\hfill \\ y_2=x_2\hfill \\ y_3=0.05\left(0.13073x_3+100a_{cs}+q_e/9-67.975\right)\hfill \end{array}\right.$$

(1)

where $x_1$, $x_2$, and $x_3$, respectively, represent drum steam pressure, electric power, and fluid density. $u_1,u_2$, and $u_3$ denote the valve positions for fuel flow, steam control, and feed-water flow, respectively. $y_3$ is the drum water level determined by

$$\begin{array}{c}{q}_e=\left(0.854u_2-0.147\right)x_1+45.59u_1-2.514u_3-2.096\hfill \\ {\alpha }_{cs}=\frac{\left(1-0.001538x_3\right)(0.8x_1-25.6)}{x_3(1.0394-0.0012304x_1)}\hfill \end{array}$$

where ${\alpha }_{cs}$ and $q_e$ are the steam quality and evaporation rate, respectively.

In practice, considering the actuators’ manipulate limitation, the inputs are usually required to satisfy the following magnitude constraints:

$$\begin{array}{c}0\le u_q\le 1,q=1,2,3\hfill \end{array}$$

(2)

and have to satisfy the following rate constraints:

$$\left\{\begin{array}{c}-0.007\le {\dot{u}}_1\le 0.007\hfill \\ -2\le {\dot{u}}_2\le 2\hfill \\ -0.05\le {\dot{u}}_3\le 0.05\hfill \end{array}\right.$$

(3)

3. Model Identification Design

The non-linear model (1) properly describes the dynamic characterizes of the boiler-turbine system, however, it encounters great difficulties in controller design due to state coupling interaction and non-linearity. In the sequel, in order to obtain a linear model, the system identification will be utilized to characterize the system dynamic behavior, as shown in the following section.

3.1. Identification Test Signal

We choose the classical generalized binary noise (GBN) proposed by [44] as the test signal. The value of GBN takes $-a$ or a. At each predetermined conversion time t, the GBN signal is converted based on the following rules:

$$\left\{\begin{array}{c}\mathrm{Prob}\{u\left(t\right)=-u(t-1)\}=p_{sw}\hfill \\ \mathrm{Prob}\{u\left(t\right)=u(t-1)\}=1-p_{sw}\hfill \end{array}\right.,$$

(4)

where $p_{sw}$ denotes the conversion probability. Since the distribution is an independent staggered distribution with probability of $p_{sw}$ at each conversion time t, the mean-value of the signal is 0. The average conversion time is

$$E\left[T_{sw}\right]=\frac{T_{min}}{p_{sw}},$$

(5)

where $T_{min}$ denotes the minimum conversion time. It means that the time of GBN signal remains unchanged during the sampling process. The power spectrum of GBN is

$${\mathsf{\Phi }}_u\left(\omega \right)=\frac{\left(1-q^2\right)T_{min}}{1-2qcos{T}_{min}\omega +q^2},q=1-2p_{sw}.$$

(6)

Remark 1.

Since the plant exhibits severe non-linear dynamics, it is difficult to obtain an accurate linear model using linearization. In order to fully obtain the prior knowledge of process behavior, the plant to be identified should be effectively excited. The excitability conditions of the input signal for the boiler-turbine system (1) are: (a) the signal length should be large enough; (b) the signal changes quickly with a large amplitude, which can generate sufficient excitation to the system. From the generation method of GBN, it can seen that GBN have the merits in setting signal length arbitrarily and having a minimum amplitude factor so that it is available to apply for multivariable model identification.

3.2. Identification and Modeling Performance

Static or steady-state model is utilized to calculate the equilibrium points, which determines operation points and generates set values. For the non-linear model (1), in a steady-state working point, the state variables are constant, therefore,

$$\left\{\begin{array}{c}{\dot{x}}_1={\dot{x}}_2={\dot{x}}_3=0\hfill \\ y_3=0\hfill \end{array}\right.$$

(7)

One can calculate several nominal working points based on (1) and (7). In this paper, we consider the case that model (1) represents a 160 MW unit, and the steady-state working points is $x_{\mathrm{eq}}={[115,85,402.759]}^{\mathrm{T}}$, $u_{\mathrm{eq}}={[0.4147,0.7787,0.5436]}^{\mathrm{T}}$, and $y_{\mathrm{eq}}={[115,85,0]}^{\mathrm{T}}$, as shown in [6]. It represents $50\%$ operation condition. For this working point, the corresponding step-response the non-linear model (1) are plotted in Figure 2.

Based on the non-linear model (1), an open-loop identification experiment is performed to produce the outputs, where Figure 3 shows the state equation in non-linear model (1) using Maltalb Simulink. For the open-loop experiment, the test signal (GBN signal) is set as the process control input. Considering that the rules for test signal stated in Remark 2, choosing the GBN signal follows that: the amplitude should not be too large, and its disturbance to the steady-state operation should be small so that the process variables do not exceed the operation limit. In the identification experiment, since there have three control inputs, we set the sampling time $T_s=T_{min}=1$ and $p_{sw}=1/15$ with different amplitude to obtain the open-loop inputs. By these parameter settings, it is able to obtain a good low-pass GBN signal that can fully excite model information. Figure 4 shows a GBN signal utilized in the identification experiment.

Operate the simulation system by superposing the GBN signals in Figure 4 to the respective control valve as control input, it is able to measure and record $n=2500$ sets of output data. For the resulting input and output data, it is split into two parts: the first 2000 sets will be utilized for identification, and the rest 500 sets will be utilized for prediction. In this procedure, it takes about 5.889 s for the identification at Matlab R2021a platform with a Intel i7-10700 CPU 2.90 GHz processor.

Based on the first 2000 sets input and output data, the N4SID is applied to estimate the system matrix parameters, mainly including the following two steps: (a) constructing a Hankel matrix based on the input and output data, and performing the singular value decomposition (SVD) on the projection of the Hankel matrix to obtain the estimated state sequence; (b) applying the least square method for the estimated state sequence to obtain system matrices $(A,B,C,D)$. In the identification, the row number of Hankel matrix is set as $i=10$ to obtain the distribution of the singular value histogram, as illustrated in Figure 5. It can be seen that the system information mainly distributes in the first three singular values, which indicates the optimal order of the model as $n=3$. Therefore, the resulting identified model is

$$\left\{\begin{array}{c}x(k+1)=Ax\left(k\right)+Bu\left(k\right)\hfill \\ y^{\prime }\left(k\right)=Cx\left(k\right)+Du\left(k\right)\hfill \end{array}\right.$$

(8)

where the system matrices are

$$\begin{array}{c}A=\left[\begin{array}{ccc}0.9086& 0.0266& -0.0797\\ 0.1084& 0.9915& 0.0376\\ 0.0077& 0.0217& 0.8732\end{array}\right],\hfill \\ B=\left[\begin{array}{ccc}-0.0174& -2.5918& -0.0190\\ -0.5269& -0.0137& 0.0768\\ -0.1285& 0.0858& 0.0616\end{array}\right],\hfill \\ C=\left[\begin{array}{ccc}-0.2205& -1.5809& -0.4900\\ 8.7136& -0.8193& 2.4134\\ -0.0047& 0.0059& -0.0235\end{array}\right],\hfill \\ D=\left[\begin{array}{ccc}0.1078& 0.0177& -0.0175\\ -0.0191& -0.8039& -0.0264\\ 0.2541& -0.0958& -0.0212\end{array}\right].\hfill \end{array}$$

In Figure 6, the comparison between of prediction values of controlled variables (CVs) and their real values are plotted. It is shown that the linearized model can accurate behave the non-linear model (1).

3.3. Minimal Realization NMSS Model

Letting $y\left(k\right)=y^{\prime }\left(k\right)-Du\left(k\right)$, model (8) can be naturally transformed into

$$\left\{\begin{array}{c}x(k+1)=Ax\left(k\right)+Bu\left(k\right)\hfill \\ y\left(k\right)=Cx\left(k\right)\hfill \end{array}\right.$$

(9)

Model (9) can be rewritten as the following input–output model:

$$\mathcal{A}\left(z^{-1}\right)y\left(k\right)=\mathcal{B}\left(z^{-1}\right)u\left(k\right)$$

(10)

where

$$\begin{array}{c}\mathcal{A}\left(z^{-1}\right)=I+a_1{z}^{-1}+\cdots +a_p{z}^{-p}\hfill \\ \mathcal{B}\left(z^{-1}\right)=b_1{z}^{-1}+\cdots +b_q{z}^{-q}\hfill \end{array}$$

where $z^{-1}$ is the backward shift operator, i.e., $z^{-1}{y}^{\prime }\left(k\right)=y(k-1)$. $a\left(z^{-1}\right)$ and $b\left(z^{-1}\right)$ are the matrix polynomials, $\{a_1,a_2,\cdots ,a_p\}$ and $\{b_1,b_2,\cdots ,b_q\}$ are the coefficient matrices with proper dimension. Based on model (10), it is obvious that the order of input and output are $n_y=n_u=3$. Hence, the coefficient matrices of input–output model (10) are

$$\begin{array}{c}{a}_1=-2.7214I_3,a_2=2.4734I_3,a_3=-5.623I_3\hfill \\ b_1=\left[\begin{array}{ccc}0.9673& 0.0496& -0.7013\\ -0.5216& 19.9122& -2.0765\\ -0.1258& -0.1155& 0.0397\end{array}\right]\hfill \\ b_2=\left[\begin{array}{ccc}0.0397& -0.1399& 0.6438\\ 1.0945& -37.8389& 1.8942\\ 0.1236& -0.0548& -0.0385\end{array}\right]\hfill \\ b_3=\left[\begin{array}{ccc}0.3278& 0.0855& -0.1970\\ -0.5278& 17.5940& -0.5774\\ -0.0403& 0.0413& 0.0124\end{array}\right]\hfill \end{array}$$

Defining the intermediate state vector $x_n\left(k\right)=[y\left(k\right),y(k-1),\cdots ,y(k-p+1),$$u(k-1),$${u(k-2),\cdots ,u(k-q+1)]}^{\mathrm{T}}$, model (10) can be rewritten as the following augmented NMSS model:

$$\left\{\begin{array}{c}{x}_n(k+1)=A_n{x}_n\left(k\right)+B_{n}u\left(k\right)\hfill \\ y(k+1)=C_n{x}_n(k+1)\hfill \end{array}\right.$$

(11)

where

$$\begin{array}{c}{A}_n=\left[\begin{array}{ccccccccc}-a_1& -a_2& \cdots & -a_{p-1}& -a_p& b_2& \cdots & b_{q-1}& b_q\\ I& 0& \cdots & 0& 0& 0& \cdots & 0& 0\\ 0& I& \cdots & 0& 0& 0& \cdots & 0& 0\\ ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& ⋮& I& 0& 0& ⋮& 0& 0\\ 0& 0& \cdots & 0& 0& 0& \cdots & 0& 0\\ 0& 0& \cdots & 0& 0& I& \cdots & 0& 0\\ ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& \cdots & 0& 0& 0& \cdots & I& 0\end{array}\right],\hfill \\ B_n={\left[\begin{array}{cccccccc}{b}_1& 0& \cdots & 0& I& 0& \cdots & 0\end{array}\right]}^{\mathrm{T}},\hfill \\ C_n=\left[\begin{array}{cccccccc}I& 0& 0& \cdots & 0& 0& 0& 0\end{array}\right].\hfill \end{array}$$

Defining a new state vector $\widehat{x}\left(k\right)={[x_n\left(k\right),y\left(k\right)]}^{\mathrm{T}}$, the augmented NMSS model (11) is reconstructed as

$$\left\{\begin{array}{c}\widehat{x}(k+1)=\widehat{A}\widehat{x}\left(k\right)+\widehat{B}u\left(k\right)\hfill \\ y\left(k\right)=\widehat{C}\widehat{x}\left(k\right)\hfill \end{array}\right.$$

(12)

where

$$\widehat{A}=\left[\begin{array}{cc}{A}_n& 0\\ C_n{A}_n& I\end{array}\right],\widehat{B}=\left[\begin{array}{c}{B}_n\\ C_n{B}_n\end{array}\right],\widehat{C}=\left[0I\right]$$

Next, the augmented NMSS model (12) can be rearranged into a canonical form by utilizing a Kalman decomposition that splits the states into following four parts:

• States which are both controllable and observable;
• States which are controllable and unobservable;
• States which are observable and uncontrollable;
• States which are uncontrollable and unobservable.

By applying controllability and observability transformations on the augmented NMSS model (12), it is shown that the augmented NMSS model (12) is transformed into the following complete Kalman decomposition form:

$$\left\{\begin{array}{c}\overline{x}(k+1)=\left[\begin{array}{cccc}{\overline{A}}_{\mathrm{co}}& 0& {\overline{A}}_{13}& 0\\ {\overline{A}}_{21}& {\overline{A}}_{\mathrm{c}\overline{\mathrm{o}}}& {\overline{A}}_{23}& {\overline{A}}_{24}\\ 0& 0& {\overline{A}}_{\overline{\mathrm{c}}\mathrm{o}}& 0\\ 0& 0& {\overline{A}}_{43}& {\overline{A}}_{\overline{\mathrm{c}\mathrm{o}}}\end{array}\right]\overline{x}\left(k\right)+\left[\begin{array}{c}{\overline{B}}_{\mathrm{co}}\\ {\overline{B}}_{\mathrm{c}\overline{\mathrm{o}}}\\ 0\\ 0\end{array}\right]u\left(k\right)\hfill \\ y\left(k\right)=\left[\begin{array}{cccc}{\overline{C}}_{\mathrm{co}}& 0& {\overline{C}}_{\overline{\mathrm{c}}\mathrm{o}}& 0\end{array}\right]\overline{x}\left(k\right)\hfill \end{array}\right.$$

(13)

where $\overline{x}\left(k\right)={[{\overline{x}}_{\mathrm{co}}\left(k\right),{\overline{x}}_{\mathrm{c}\overline{\mathrm{o}}}\left(k\right),{\overline{x}}_{\overline{\mathrm{c}}\mathrm{o}}\left(k\right),{\overline{x}}_{\overline{\mathrm{c}}\overline{\mathrm{o}}}\left(k\right)]}^{\mathrm{T}}$. Therefore, taking the controllable and observable part from the Kalman decomposition yields the following minimal realization state-space model:

$$\left\{\begin{array}{c}{\overline{x}}_{\mathrm{co}}(k+1)={\overline{A}}_{\mathrm{co}}{\overline{x}}_{\mathrm{co}}\left(k\right)+{\overline{B}}_{\mathrm{co}}u\left(k\right)\hfill \\ y\left(k\right)={\overline{C}}_{\mathrm{co}}{\overline{x}}_{\mathrm{co}}\left(k\right)\hfill \end{array}\right.$$

(14)

Remark 2.

Identification just determines the model order, but the resulting state vector of identified model (8) are of dimensionless (without physical significance). Augmenting the historic input and output of model (8) as the system state ensures the state variable can be measured and calculated. Transforming the NMSS model into a minimal realization formulation ensures its controllability and observability, which is also helpful to reduce the computation burden and avoid the curse of dimensionality.

3.4. Control Objectives

For the minimal realization model (14), the control objectives are to design an effective controller such that the closed-loop system is asymptotically stable, and the constraints are always satisfied over the control horizon. The two objectives are to be achieved by manipulating the control input u to steer the state to the equilibrium point, while the control inputs are required to maintain a desired range. Since the original non-linear system has been transformed into a linearized model, it ensures that the resulting system state of the linearized NMSS model is measurable. Subsequently, it is easy to design a state feedback-based MPC controller so that the drum pressure $y_1$, electric power $y_2$, and the drum level $y_3$ converge to the equilibrium point, which will be illustrated in the sequel.

4. Finite-Horizon MPC Formulation

In order to perform MPC, an optimization problem of the system to be controlled is needed. It includes a prediction model of the physical system and its constraint, as well as a performance function. The following subsections will give the details of designing a finite-horizon MPC for the minimal realization state-space model (14). The proposed control structure is shown in Figure 7. MPC is an advanced computer control technology that on-line optimizes both the current values of control moves and the whole future control sequence over the control horizon. The MPC controller main has three ingredients: prediction model, constraints, and cost function.

Consider the following input and state constraints:

$$\begin{array}{c}-\underline{u}\le u(k+i)\le \overline{u},\hfill \\ -\underline{\psi }\le \mathsf{\Psi }{\overline{x}}_{\mathrm{co}}(k+i+1)\le \overline{\psi },i\ge 0,\hfill \end{array}$$

(15)

where $\underline{u}:={[{\underline{u}}_1,{\underline{u}}_2,\cdots ,{\underline{u}}_{n_u}]}^{\mathrm{T}}$, $\overline{u}:={[{\overline{u}}_1,{\overline{u}}_2,\cdots ,{\overline{u}}_{n_u}]}^{\mathrm{T}}$, ${\underline{u}}_j>0$, ${\overline{u}}_j>0$, $j\in \{1,\cdots ,n_u\}$; $\underline{\psi }:={[{\underline{\psi }}_1,{\underline{\psi }}_2,\cdots ,{\underline{\psi }}_g]}^{\mathrm{T}}$, $\overline{\psi }:={[{\overline{\psi }}_1,{\overline{\psi }}_2,\cdots ,{\overline{\psi }}_g]}^{\mathrm{T}}$, ${\psi }_s>0$, ${\overline{\psi }}_s>0$, $s\in \{1,\cdots ,g\}$; $\mathsf{\Psi }\in {\Re }^{g\times n_{\mathrm{co}}}$.

4.1. Performance Function and Terminal Control Law Design

Consider the following finite-horizon cost function

$$\begin{array}{cc}\hfill J\left({\overline{x}}_{\mathrm{co}}\left(k\right)\right)=& \sum _{i=0}^{N-1}\left[\parallel {\overline{x}}_{\mathrm{co}}{(k+i|k)\parallel }_W^2+{\parallel u(k+i|k)\parallel }_R^2\right]+\hfill \\ \hfill & \parallel {\overline{x}}_{\mathrm{co}}{(k+N|k)\parallel }_P^2,\hfill \end{array}$$

(16)

where W, R, and P are the positive-definite weighting matrices; N represents the control horizon. The terminal weighting matrix P is assumed to satisfy [45]

$$P\ge {(A_{\mathrm{co}}+B_{\mathrm{co}}F)}^{\mathrm{T}}P(A_{\mathrm{co}}+B_{\mathrm{co}}F)+W+F^{\mathrm{T}}RF,$$

(17)

where F is the feedback gain matrix.

Letting $X=P^{-1}$, $F=Y{X}^{-1}$ and utilizing Schur complement on (17) obtain

$$\left[\begin{array}{cccc}X& ★& ★& ★\\ A_{\mathrm{co}}X+B_{\mathrm{co}}Y& X& ★& ★\\ W^{1/2}X& 0& I& ★\\ R^{1/2}Y& 0& 0& I\end{array}\right]\ge 0$$

(18)

Lemma 1

([46]). Suppose that there exist symmetric matrices $X=P^{-1}$, $\{Z,\mathsf{\Gamma }\}$ and matrix Y satisfying (17) and

$$\begin{array}{c}\left[\begin{array}{cc}Z& ★\\ Y^{\mathrm{T}}& X\end{array}\right]\ge 0,Z_{jj}\le {\overline{u}}_{j,inf}^2,j\in \{1,\cdots ,n_u\}\hfill \end{array}$$

(19)

$$\begin{array}{c}\left[\begin{array}{cc}X& ★\\ \mathsf{\Psi }(A_{\mathrm{co}}X+B_{\mathrm{co}}Y)& \mathsf{\Gamma }\end{array}\right]\ge 0,{\mathsf{\Gamma }}_{ss}\le {\overline{\psi }}_{s,inf}^2,s\in \{1,\cdots ,g\}\hfill \end{array}$$

(20)

where $u_{j,inf}=min\{{\underline{u}}_j,{\overline{u}}_j\}$, ${\psi }_{s,inf}=min\{{\underline{\psi }}_s,{\overline{\psi }}_s\}$; $Z_{jj}\left({\mathsf{\Gamma }}_{ss}\right)$ is the $j\left(s\right)$-element of $Z\left(\mathsf{\Gamma }\right)$. When ${\overline{x}}_{\mathrm{co}}(k+N)\in {\epsilon }_P=\{z\in {\Re }^n|z^{\mathrm{T}}Pz\le 1\}$ is satisfied and $u(k+i+N)=Y{X}^{-1}{\overline{x}}_{\mathrm{co}}(k+i+N)$ is applied, $\forall i\ge 0$, then the closed-loop system is said to be asymptotically stable, $x(k+i+N)$, $\forall i\ge 0$ always lies inside of ellipsoid ${\epsilon }_P$, and constraint (15) holds for all $i\ge N$.

Remark 3.

Since model (14) has a minimal realization formulation, it implies that one can directly state feedback control law, without designing an observer. In Lemma 1, we show the details for designing a terminal control law using model (14) and handling the constraints based on the invariant set. In deriving the inequality from (17) to linear matrix inequality (LMI) (18), it is a congruence transformation that does not introduce any conservatism.

4.2. Overall Optimization Problem and Stability Analysis

The optimal control sequence $\left\{u\right(k\left|k\right),\cdots ,u(k+N-1|k\left)\right\}$ is obtained by solving the following finite-horizon MPC optimization problem at each time k:

$$\underset{u\left(k\right|k),\cdots ,u(k+N-1\left|k\right)}{min}J\left({\overline{x}}_{\mathrm{co}}\left(k\right)\right),$$

(21)

$$\begin{array}{cc}\hfill \mathrm{s}.\mathrm{t}.& -\underline{u}\le u(k+i|k)\le \overline{u},i\in \{0,1,\cdots ,N-1\}\hfill \end{array}$$

(22)

$$-\underline{\psi }\le \mathsf{\Psi }{\overline{x}}_{\mathrm{co}}(k+i+1|k)\le \overline{\psi },i\in \{0,1,\cdots ,N-1\}$$

(23)

$$\parallel {\overline{x}}_{\mathrm{co}}{(k+N|k)\parallel }_P^2\le 1.$$

(24)

The feasible initial state set is a set of initial state points ${\overline{x}}_{\mathrm{co}}\left(0\right)$ under the open-loop optimal control sequence $\left\{u\right(k\left|k\right),\cdots ,u(k+N-1|k\left)\right\}$ yielded by the optimization problem (21)–(24), such that the system state ${\overline{x}}_{\mathrm{co}}\left(k\right)$ lies in the terminal region ${\epsilon }_P$ at instant $k+N$, which is defined as

$$\begin{array}{cc}\hfill & \mathcal{F}(P,N)=\{{\overline{x}}_{\mathrm{co}}\left(0\right)\in {\Re }^n|\exists u\left(i\right)\in \mathcal{U},i\in \{0,\cdots ,N-1\}\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& {\overline{x}}_{\mathrm{co}}(i+1)\in \mathcal{X},{\overline{x}}_{\mathrm{co}}\left(N\right)\in {\epsilon }_P\}.\hfill \end{array}$$

(25)

Following from the above discussions, the following invariant property result is naturally obtained for optimization problem (21)–(24).

Lemma 2.

Assume that ${\overline{x}}_{\mathrm{co}}\left(k\right)\in \mathcal{F}(P,N)$. There exist $\kappa >0$ and $u(k+i|k)\in \mathcal{U},i\in \{0,1,\cdots ,N-1\}$, such that ${\left|u(k+i|k)\right|}^2\le \kappa {\left|{\overline{x}}_{\mathrm{co}}\left(k\right)\right|}^2$, ${\overline{x}}_{\mathrm{co}}(k+i+1|k)\in \mathcal{X},i\in \{0,1,\cdots ,N-1\}$, and ${\overline{x}}_{\mathrm{co}}(k+N|k)\in {\epsilon }_P$.

For the finite MPC, the terminal cost function, the terminal-state region and the terminal control law, shown in the above sections, are the keys to guarantee the close-loop stability. Note that the non-linear system is described by a linear state-space model in the minimal realization formulation, so the nominal stability should be discussed. Based on the presented results, in the following, we will present the stability results of the closed-loop system, composed by model (14) and control moves $\left\{u\right(k\left|k\right),\cdots ,u(k+N-1|k\left)\right\}$.

Theorem 1

(Nominal stability). Suppose that (18)–(20) and ${\overline{x}}_{\mathrm{co}}\left(0\right)\in \mathcal{F}(P,N)$ are satisfied, then (22)–(24) are always feasible for any $k\ge 0$. The closed-loop system is said to be exponential stability by implementing the optimal control input $u^{*}\left(k\right|k)$ in a receding horizon way.

Proof. 

Suppose that, at each k, the optimal solution for problem is denoted as $u^{*}(k+i|k)$. Based on Lemma 2, a feasible solution at $k+1$ can be constructed as

$$\left\{\begin{array}{c}u(k+i|k+1)=u^{*}(k+i|k),i\in \{1,\cdots ,N-1\},\hfill \\ u(k+N|k+1)=Fx(k+N|k+1)\hfill \end{array}\right.$$

(26)

Considering the invariant property, it is clear that (21)–(24) are feasible for any $k\ge 0$. For the purpose of establishing the exponential stability, we should illustrate that there exist $a,b,c(0\le a,b,c\le \infty )$, such that

$${a\parallel x\left(k\right)\parallel }_2^2\le J^{*}\left(x\left(k\right)\right)⩽{b\parallel x\left(k\right)\parallel }_2^2,\Delta J^{*}\left(x(k+1)\right)<-c{\parallel x\left(k\right)\parallel }_2^2$$

(27)

where $\Delta J^{*}\left(x(k+1)\right)=J^{*}\left(x(k+1)\right)-J^{*}\left(x\left(k\right)\right)$. Once condition (27) holds, then $J^{*}\left(x\left(k\right)\right)$ can serve as the Lyapunov function of exponential stability. It is clear that $J^{*}\left(x\left(k\right)\right)\ge x{\left(k\right)}^{\mathrm{T}}Wx\left(k\right)\ge {\lambda }_{min}\left(W\right){\parallel x\left(k\right)\parallel }_2^2$ holds, so it is available to choose $a={\lambda }_{min}\left(W\right)$. Following from Lemma 2, one has

$$\begin{array}{cc}\hfill J^{*}\left(x\left(k\right)\right)& \le \sum _{i=0}^{N-1}\left[{\parallel x(k+i|k)\parallel }_W^2+{\parallel u(k+i|k)\parallel }_R^2\right]+{\parallel x(k+N|k)\parallel }_P^2\hfill \\ \hfill & \le \left[(N+1){\mathcal{A}}^2(1+N\parallel B\parallel \sqrt{\kappa }{)}^2·max\{{\lambda }_{max}\left(W\right),{\lambda }_{max}\left(P\right)\}\right.\hfill \\ \hfill & \left.+N\kappa {\lambda }_{max}\left(R\right)\right]{\parallel x\left(k\right)\parallel }_2^2\hfill \end{array}$$

(28)

where $\mathcal{A}={max}_{i\in \{0,1,\cdots ,N\}}\parallel A^i\parallel$ and $\kappa >0$. Therefore, it is easy to choose $b=(N+1){\mathcal{A}}^2(1+N\parallel B\parallel \sqrt{\kappa }{)}^2·max\{{\lambda }_{max}\left(W\right),{\lambda }_{max}\left(P\right)\}+N\kappa {\lambda }_{max}\left(P\right)$. At time $k+1$, the performance cost is denoted as $\overline{J}\left(x(k+1)\right)$ under (26), then it is derived that

$$\begin{array}{cc}\hfill J^{*}\left(x\left(k\right)\right)& \ge {\parallel x\left(k\right)\parallel }_W^2+{\parallel u\left(k\right)\parallel }_R^2+\overline{J}\left(x(k+1)\right)\hfill \\ \hfill & \ge {\parallel x\left(k\right)\parallel }_W^2+{\parallel u\left(k\right)\parallel }_R^2+J^{*}\left(x(k+1)\right).\hfill \end{array}$$

(29)

Inequality (29) shows that $\Delta J^{*}\left(x(k+1)\right)\le -{\parallel x\left(k\right)\parallel }_W^2-{\parallel u\left(k\right)\parallel }_R^2\le -{\lambda }_{min}\left(W\right){\parallel x\left(k\right)\parallel }_2^2$, then it is available to choose $c={\lambda }_{min}\left(W\right)$. Hence, $J^{*}\left(x\left(k\right)\right)$ is Lyapunov function to prove the exponential stability. The proof is complete. □

Remark 4.

Considering the inherit non-linearity and the optimization mechanism of MPC, it is difficult to obtain the numerical solution of decision sequence $\left\{u\right(k\left|k\right),u(k+1|k,\cdots ,u(k+N-1|k\left)\right\}$ even in the case of applying a quadratic performance index. Moreover, the computation burden is always the core issue to be solved. In order to avoid a complex non-linear optimization problem and lead to a large computation burden, we identify the non-linear system to obtain a linear model and transform it into a minimal realization formulation to decrease the system dimension, while guaranteeing fully state measurability and observability.

Remark 5.

For the calculation of identification and controller parameter optimization, error mainly comes from model-mismatch and numerical arithmetic. In the linearization, the non-linear system is assumed to operate a relatively ideal working environment by ignoring isentropic efficiency or isothermal waveforms during the expansion of the working medium in the turbogenerator, this may lead to modeling error during identification process. It is the main reason for affecting the performance of controller. However, the designed MPC applies the feedback control law and optimizes the control sequence in a receding horizon to cope with this error, and maintain good robustness. For the numerical arithmetic, there exists numerical error because of calculation precision.

5. Simulation Results and Analysis

This section aims to demonstrate the effectiveness of the finite-horizon MPC based on the minimal realization model (14) and discuss the simulations results to quantify the advantages and disadvantages.

5.1. Simulation Settings

Model described by (14) is utilized for the simulated process. The platform of Matlab 2021a installed on a computer with Intel i7-10700 CPU 2.90 GHz, 16 G RAM is utilized for simulation with a sampling time $T_{s,MPC}=1$.

In the discrete-time domain and MPC implementation, the manganite constraints (2) and the rate constraints (3) on inputs are rewritten as

$$\left\{\begin{array}{c}0\le u_1(k+i|k)\le 1,i=0,\dots ,N\hfill \\ 0\le u_2(k+i|k)\le 1,i=0,\dots ,N\hfill \\ 0\le u_3(k+i|k)\le 1,i=0,\dots ,N\hfill \end{array}\right.$$

(30)

and

$$\left\{\begin{array}{c}-0.007\le \Delta u_1(k+i|k)\le 0.007,i=0,\dots ,N\hfill \\ -1\le \Delta u_2(k+i|k)\le 0.02,i=0,\dots ,N\hfill \\ -0.05\le \Delta u_3(k+i|k)\le 0.05,i=0,\dots ,N\hfill \end{array}\right.$$

(31)

where $\Delta u_n(k+i|k)=u_n(k+i|k)-u_n(k+i-1|k)$ for all $n=1,2,3$.

Consider the control objective described in Section 3, the following two cases will be presented, so that the results, efficiency and merits of the designed control method are well understood. Case 1 is the scenario for applying classical linear quadratic regulator (LQR) method [47], and Case 2 represent the scenario for applying the proposed control method. In the two cases, both cases take the same weighing matrices $W=I,R=I$ to make a fair comparison.

Case 1.

Since the LQR can not deal with constraints, Equations (30) and (31) will not be imposed.

Case 2.

Based on the resulting minimal realization model$\{{\overline{A}}_{\mathrm{co}},{\overline{B}}_{\mathrm{co}},{\overline{C}}_{\mathrm{co}}\}$, the following optimization problem is off-line solved

$$\underset{X,Y}{max}trace\left(X\right)s.t.\left(18\right)–\left(20\right)$$

(32)

Problem (32) is a standard convex optimization problem that can be directly solved by Yalmip toolbox to obtain X. One has the following terminal weighting matrix

$$P=\left[\begin{array}{ccccc}1.1852& 0.0051& 0.0125& -0.1211& 0.1697\\ 0.0051& 2.3474& -0.0809& 0.0757& -0.1275\\ 0.0125& -0.0809& 2.7050& -0.5774& -0.7972\\ -0.1211& 0.0757& -0.5774& 5.2791& -0.1514\\ 0.1697& -0.1275& -0.7972& -0.1514& 5.4558\end{array}\right].$$

In Case 2, besides the modeling and constraints, the control parameters should also be set. An important parameter is control horizon N since it affects control quality and computational complexity. Increasing N will improve the optimality, but increases computation burden. For making a balance between optimality and computation, $N=5$ is properly determined under several comparison experiments. Based on the above parameter settings, problem (21)–(24) is online solved at each time k to optimize a control sequence $\left\{u\right(k\left|k\right),u(k+1|k,\cdots ,u(k+N-1|k\left)\right\}$, and perform the first control input $u\left(k\right|k)$.

5.2. Simulation Results

Closed-loop simulation profiles generated by the two cases are shown in Figure 8, Figure 9, Figure 10, Figure 11, Figure 12 and Figure 13. Following from the responses in Figure 8 and Figure 9, it is obvious that the LQR method can not stabilize model (14). However, the control input signals and output responses of the designed MPC controller plotted in Figure 10, Figure 11, Figure 12 and Figure 13 converge to the steady-state working point. In Figure 10, the constraint of each control action and its rate constraint are always satisfied, this verifies the advantages of MPC at handling physical constraints. Figure 11, Figure 12 and Figure 13 show the outputs (electric power, drum steam pressure, and drum water level) that can eventually maintain stable control by applying the optimal control inputs, respectively. In practice, amplitude constraints have to be considered to guarantee a physical safe operation for the plant, and rate constraints avoids extensive controller parameter tuning. From the comparison, it can be seen that the MPC controller based on the minimal realization state-space model steer the inputs and outputs of the closed-loop system to the equilibrium point, which verify the effectiveness of the proposed MPC. In conclusion, the simulations results illustrate that control performance of MPC controller is superior to a LQR controller.

Further important investigation is to evaluate the performances of the developed models in real-time applications, so we consider the computation complexity of MPC optimization problem. There mainly involve two problems: off-line optimization (32), and online optimization (21)–(24). They are convex optimization problem composed by linear matrix inequalities (LMIs). Usually, the interior point method is applied to solve the feasible solution of LMIs. It is a polynomial complexity algorithm. The complexity is proportional to ${\mathfrak{K}}^3\mathfrak{L}$ for a given $\alpha$, where $\mathfrak{K}$ and $\mathfrak{L}$ represent the number of scalar LMI variables and the number of scalar rows in LMIs, respectively. It is obvious that the system dimensionality and control horizon affect the complexity of an optimization problem. As stated in Remark 4, the linear model utilized for optimization is transformed into a minimal realization formulation, so we ensure the eventually MPC problem of good efficiency from modeling. In the simulation, the total computation time for carrying out off-line and online stages is 18.388 s. Following the computation time, the implementation efficiency of a linear finite-horizon MPC based on the developed minimal realization model in the real-time performance has been effectively examined.

6. Conclusions

The focus of this paper is to investigate the control problem of the boiler-turbine system based on a system identification method and the finite-horizon MPC. The subspace identification method is firstly employed to obtain a linear state-space model that can properly reflects the dynamic characteristics. By transforming the state-space model into the input–output model and extending its inputs and outputs as a new augmented state, a NMSS model of state measurable is established, and transformed into a minimal realization formulation in order to facilitate state feedback controller and also reduce the computation complexity. A finite-horizon MPC controller, which utilizes a cost function incorporating with a finite-horizon cost and a terminal cost, is designed to guarantee the nominal stability. The physical constraints are handled by an invariant ellipsoid constraint with the terminal control law. Through theoretical proof and comparison simulation, the approach is proven to give an effective control strategy for non-linear system to maintain satisfactory performance as well as a well implementation efficiency. For future works, an interesting plan is to study other issues and make improvements, such as considering the robustness against parametric uncertainty and economic optimization problem.