Improving Simultaneous Cooling and Power Load-Following Capability for MGT-CCP Using Coordinated Predictive Controls

Abstract

The distributed energy system is an energy supply method built around the end users, which can achieve energy sustainability and reduce emissions compared to traditional centralized energy systems. The micro gas turbine (MGT)-based combined cooling and power (CCP) system has received renewed attention as an important distributed energy system technology due to its substantial energy savings and reduced emission levels. The task of the MGT-CCP system is to quickly adapt to changes in various renewable energy sources to maintain the balance in energy supply and demand in a distributed energy system. Therefore, it is imperative to improve the load tracking capability of the MGT-CCP system with advanced control technologies toward achieving this goal. However, the difficulty of controlling the MGT-CCP system is that the MGT responds very fast while CCP responds very slowly. To this end, the dynamic characteristics and nonlinear distribution of the MGT and CCP processes are analyzed, and a coordinated predictive control strategy is proposed by utilizing the generalized predictive control for the MGT system and the Hammerstein generalized predictive control for the CCP system. The coordinated predictive control of generalized predictive control and Hammerstein generalized predictive control was implemented in an 80 kW MGT-CCP simulator to verify the effectiveness of the proposed method. The simulation results show that compared with PID and MPC, the proposed control method not only can greatly improve simultaneous cooling and power load-following capability, but also has the best control effect when accessing with renewable energy.

1. Introduction

With the continuous development of the economy, many countries are facing severe challenges in preserving the environment and reducing energy consumption so that a sustainable development of the society can be ensured. A core solution for reducing emissions and improving energy efficiency is the use of distributed energy systems (DESs) [1,2]. The DES is usually composed of a variety of modular and small-scale generation technologies, located near end-users, and can be regarded as an important complement to traditional centralized power grids. They have many advantages, such as low transmission loss, low emissions, and flexibility for utilizing a variety of energy sources, including fossil fuels and renewable energy sources (RES) [3].

As a key form of DES, micro gas turbine (MGT)-based combined cooling, heating and power systems (CCHPs) have attracted much attention worldwide because of their energy-efficiency, low-emissions and operational cost-saving characteristics [4,5,6,7]. Since the waste heat of the MGT can be continuously utilized as a heat source for a water heater or a cooling device, the average energy efficiency of the MGT-CCHP system can reach 80%, whereas the average thermal efficiency for a combined cycle gas turbine (CCGT) is approximately 50% [8]. In addition, compared with the traditional centralized power generation system, the MGT-CCHP system is smaller in size so that it can be flexibly installed in small residential or commercial districts.

At present, the research on MGT-CCHP system mainly focuses on static system configuration [9], operation optimization [10] and performance evaluation [11]. In [9], a CCHP microgrid composed of photovoltaic cell, gas turbines, gas boiler, thermal storage tank, absorption chiller, electric chiller, as well as cooling, heating, and power load is studied. In [10], a mixed-integer linear programming model was established for the optimal design of the distributed energy resource system. In [11], the performances of CCHP system in different operation modes are compared. Although the dynamic characteristics of the MGT-CCHP system have been studied in recent years [12,13], the control method used for the MGT-CCHP system is still the conventional PI/PID control. In order to improve the performance of the PI/PID controller, intelligent methods such as fuzzy logic [14,15], neural network [16], firefly algorithm [17] are used to tune the parameters of the PID controller. However, due to the complex dynamic characteristics of MGT-CCHP system, such as nonlinearity, large thermal inertia, multi-variable strong couplings, unknown disturbances and so on, the single-input-single-output (SISO) PID control is no longer applicable to meet the performance requirement, even if the parameters of the PID are well tuned.

Model predictive control (MPC) is an advanced method to deal with optimization control of multivariable, large inertia and constrained systems [18]. In order to solve the control problem of MGT-CCHP system, an observer-based model predictive control algorithm is used [19]. The simulation results show that the control performance of the MPC is obviously improved compared with the PID controller for MGT-CCHP system. A self-tuning MPC control method is proposed to improve the performance of the MPC for MGT-CCHP system [20]. In [21,22], a data-driven MPC controller is proposed to improve the operation quality of the MGT-CCHP system. In [23], a supervisory MPC controller is proposed to improve the economic efficiency. However, these MPC controllers designed for the MGT-CCHP system do not consider the nonlinearity of the system. When the MGT-CCHP system runs on a wide range of load changes, the nonlinearity of the MGT-CCHP system will lead to the performance degradation of the controller, even if the controller parameters are well tuned. Therefore, it is necessary to design a nonlinear MPC controller to improve control performance. In general, MGT-CCHP system is divided into two sets of systems in actual operation, MGT-based combined cooling and power system (MGT-CCP) and MGT-based combined heating and power system (MGT-CHP); consequently, it is more practical to study control methods for these two kinds of combined systems individually.

This paper proposes a coordinated predictive control by using a generalized predictive control (GPC) for MGT system and a Hammerstein generalized predictive control (HGPC) for the CCP system. First, the dynamic characteristics and nonlinearity distribution of the MGT and CCP processes are analyzed by step response and Vinnicombe gap (V-gap) metric, which provides a basis for the selection of advanced control strategies. Secondly, the Hammerstein identification method of the single-input and single-output (SISO) system proposed in [24] is extended to the multi-input single-output (MISO) system, which is used to model the nonlinear CCP system. Third, the HGPC method was proposed to improve the performance of the predictive control for the nonlinear CCP system. Finally, the coordinated predictive control strategy is designed with the hybrid of GPC and HGPC for the MGT-CCP system to meet the simultaneous cooling and power load-tracking requirements.

The major contribution of this study is that a coordinated predictive control strategy with GPC and HGPC is designed to overcome the nonlinear and multiple-variable control problems in the MGT-CCP system.

The minor contributions of this study are summarized below:

(1)

The nonlinearity distribution of the MGT-CCP system is first analyzed using V-gap;

(2)

The MISO Hammerstein model of the CCP system is first established to capture its nonlinearity, which is beneficial to the design of the controller;

(3)

The proposed GPC-HGPC coordinated control strategy has been validated to effectively improve the simultaneous cooling and power load tracking capability of the MGT-CCP system when accessing with renewable energy.

The proposed coordinated predictive controller is implemented on a 80 kW MGT-CCP simulator. The rest of the article is organized as follows: Section 2 analyzes the dynamics and nonlinearity of MGT-CCP system; Section 3 introduces the identification of multivariate Hammerstein ARMAX models; The Hammerstein-GPC controller and the coordinated control strategy for the MGT-CCP are proposed in Section 4. The simulation results and conclusions are given in the Section 5 and Section 6, respectively.

2. Dynamics and Nonlinearity Analysis of MGT-CCP

A MGT-CCP system, composed of an 80 kW full heat-regenerated-cycle MGT and a 75 kW lithium bromide single effect absorption chiller, is considered in this paper. The MGT-CCP system schematic diagram is shown in Figure 1. Fuel is supplied to the MGT to generate the required electric power, while the waste heat emitted from MGT is recycled to drive the refrigeration system to produce chilled water. The MGT-CCP system can generate electricity and chilled water for the power and cooling load requirements at the same time.

The study on control strategy is based on a MGT-CCP emulator, which is built on the MATLAB R2016b environment and consists of dozens of major lumped-parameter modules and other functional components. Built on three conservation laws and thermal property equations, the emulator shows quite close dynamic characteristics to the primary MGT-CCP system. The detailed process of modeling is referenced in [25] and the main design data of the MGT-CCP system is given in Appendix A.

It is worth noting that the exhaust temperature of the fully heat-regenerated MGT is 277 °C, which is much lower than the non-heat-regenerated MGT of about 615 °C. Considering the limitation of heat exchange efficiency and avoiding the corrosion of low temperature flue gas in the pipeline, the exhaust temperature of the final flue gas of MGT-CCP is designed to be 170 °C. Therefore, the thermal power generated by the chiller is about 75 kW. A ratio close to 1.0 between thermal power and electrical power means that power generation efficiency is more preferentially improved in such an all-heat regeneration MGT-CCP system.

In order to help designing the controller, the MGT-CCP system can be regarded as a double-input-and-double-output dynamic system. The manipulate variables are fuel $u_1(\mathrm{kg}/\mathrm{s})$ and refrigerant valve opening $u_2(\%)$. The output variables are power output $y_1(\mathrm{kw})$ and chilled water temperature $y_2{(}^{°}\mathrm{C})$.

Table 1. Typical operating points of MGT-CCP system.

Operating Points	${\mathit{y}}_1\left(\mathbf{kw}\right)$	${\mathit{y}}_2\left({}^{\mathbf{°}}\mathbf{C}\right)$
#1	80	8
#2	76.85	8.2
#3	73.5	8.5
#4	69.69	8.8
#5	65.18	9.1

is the typical operating points of the MGT-CCP system, where $y_1$ is power output, and $y_2$ is chilled water temperature.

Remark 1: The chilled water temperature, rather than the cooling power, is chosen as the operating point for the MGT-CCP system. In the MGT-CCP system, the cooling power is proportional to the product of the temperature difference and the chilled water flow. Since the flow rate of chilled water is assumed to be constant, from a control point of view, the chilled water temperature can represent the cooling power. Although the water flow is disturbed during operation, the chilled temperature can indicate the required cooling power generation. In addition, temperature has the advantage that it can be measured directly.

Based on the simulator, the following analysis will use two methods to expose several control issues for the MGT-CCP system. First, an open-loop step response test is performed to show the change in dynamic behavior. Second, the gap metric method is used to determine the degree of nonlinearity of the system by measuring the dynamic difference between models linearized at different operating points; thus, providing a quantitative description of the nonlinearity distribution of the MGT-CCP process.

2.1. Dynamic Analysis

Suppose the system is running at operating point #3 in Table 1. At 500 s and 3000 s, a +10% step signal of the corresponding steady-state value is added to the fuel and the refrigerant valve opening, respectively. The response of the plant is shown in Figure 2.

Figure 2 shows that the power output responds very quickly to the fuel step command because it has very fast dynamic response characteristics. In addition, Figure 2 shows that as the fuel step increases, the chilled water temperature slowly decreases to a constant value. The reason for the slow response of the chilled water is caused by exhaust gas transport and heat transfer in the chiller. It is easy to understand that the more fuel entering the MGT and the more exhaust gas entering the chiller, the greater the power of the chiller. Therefore, when the fuel is increased, the chilled water temperature is lowered. At 3000 s, the refrigerant valve opening increases, resulting in a decrease in exhaust gas entering the chiller, thus causing a decrease in chilled water temperature. As seen from Figure 1, the chiller is located behind the MGT, so the power output does not change when the refrigerant valve opening changes.

The large inertia of the chilled water temperature is due to the exhaust gas transmission and thermal transferring process in the chiller. First the lithium bromide solution is heated in the generator by the exhaust gas which is a desorption process of solution. The vapour out of the solution condenses into the refrigerant water by releasing its heat to the cooling water in the condenser. The refrigerant water flows through the throttle valve and then evaporates in the evaporator and absorbs heat from the chilled water and cools it down. All these processes in the flowchart of the chiller, shown in Figure 1, contributes to a large inertial from the fuel change to the variation of chilled water temperature. Another input the refrigerant valve is a bypass valve for the exhaust gas to the chiller, thus it influences the chilled water temperature in the same direction slowly.

To conclude, the power output is affected only by fuel, so for the safety and flexibility of MGT-CCP, a separate controller is designed for the control loop of power output. The chilled water temperature is affected by both fuel and refrigerant valve opening. There is a strong coupling between the variables. Therefore, when designing the control loop of the chilled water temperature $y_2$, it is necessary to coordinate the control with the control loop of the power output $y_1$. In addition, the chilled water temperature response has a large inertia, so it is necessary to study advanced control methods to replace the traditional PID control.

2.2. Nonlinear Analysis via Gap-Metric

Vinnicombe gap (V-gap) metric has been successfully applied to multi-model modeling and control of power plant boilers-turbine [26] and spacecraft attitude control [27], providing an appropriate partition of the operating region and the selection of local models, making it possible to use a minimum number of linear models to approximate the nonlinear behavior of the system.

The value of V-gap is a measure of the difference between two linear models. Suppose $P_1$ and $P_2$ are the transfer functions of linear models at any two operating points of the MGT-CCP system. The V-gap between $P_1$ and $P_2$ can be computed by:

$$v_g\left(P_1,P_2\right)=\max \left\{\underset{Q\in H_{\infty }}{\inf }{\Vert \left[\begin{array}{c}{M}_1\\ N_1\end{array}\right]-\left[\begin{array}{c}{M}_2\\ N_2\end{array}\right]Q\Vert }_{\infty },\underset{Q\in H_{\infty }}{\inf }{\Vert \left[\begin{array}{c}{M}_2\\ N_2\end{array}\right]-\left[\begin{array}{c}{M}_1\\ N_1\end{array}\right]Q\Vert }_{\infty }\right\}$$

(1)

where $\left(N_1,M_1\right)$ and $\left(N_2,M_2\right)$ are the elements of the normalized right coprime factorization of $P_1$ and $P_2$ as: $P_1=N_1{M}_1^{-1}$ and $P_2=N_2{M}_2^{-1}$, and $Q$ is a matrix parameter which has finite H-infinity norm. The detailed process of gap calculation is referenced in [28], which has been programmed as a $gap\left(P_1,P_2\right)$ function in the MATLAB toolbox. The value of V-gap is between 0 and 1. The greater the value of V-gap, the greater the difference between the two linear models, that is, the stronger the nonlinearity of the region.

This section proposes the use of V-gap to quantitatively analyze the nonlinearity distribution of the MGT-CCP process over the operating range (power output: 65.18 kW–80 kW; chilled water temperature: 8 °C–9.1 °C).

A local linear model is first developed for a specific operating point for gap metric calculations. At each operating point, a random identification signal with a steady-state value of ±2% for 3000 s is applied to the fuel and the refrigerant valve opening path, respectively, to generate a corresponding output signals. A linear model for each operating point in Table 1 is built using the System Identification Toolbox in MATLAB based on the resulting input-output data. The gap metric between adjacent linear models is then calculated to indicate the level of non-linearity along the corresponding region as shown in Figure 3.

Figure 3a is a nonlinearity distribution of the MGT-CCP process under varying fuel. The results show that the gaps between adjacent operating points are small and uniform, implying that one linear model is sufficient for operating points #1 through #5. Therefore, in Section 4, a linear controller is designed for the control loop of the power output.

Figure 3b shows the nonlinearity distribution of the MGT-CCP process under varying refrigerant valve opening. The gaps between the operating points #2 and #3 and #4 and #5 show extreme nonlinearity, and the gap between the operating points #3 and #4 also shows nonlinearity. The gap between the operating points #1 and #2 shows the least nonlinearity. The physical processes of the chiller include desorption and absorption of the lithium bromide solution, phase change of the liquid, and various heat-transferring processes. Its physical structure is complex, which leads to its inevitable nonlinear dynamics. Therefore, in Section 4, a nonlinear controller is designed for the control loop of the cold-water temperature.

3. Identification of Multivariable Hammerstein Model

Predictive control offers its own distinctive features for dealing with coupled variables in process control, such as large inertia, input constraints, etc. It can also handle the predictive control of nonlinear plant well with a nonlinear predictive model. Therefore, predictive control is an appropriate solution for the combined cooling and power load-following control of MGT-CCP.

In the design framework of the predictive controller, designing a prediction model takes the lead. The Hammerstein model consists of a static memoryless nonlinear element and a dynamic linear model [29,30]. It has the advantages of small computing load, strong nonlinearity-approximating ability, and favorable controller design. Since the nonlinearity of the cooling system is very strong, the Hammerstein model is taken to be the prediction model for generalized predictive controller (GPC) whose algorithm matches well with the prediction model. An effective Hammerstein model identification method for SISO system has been proposed in [24], which we extend to multi-input-single-output (MISO) systems for modeling the MGT-CCP.

Assuming that a system with $n_u$-input and single-output can be described by the following Hammerstein model:

$$A(z^{-1})y_t=B(z^{-1}){\overline{u}}_t$$

(2)

where $A(z^{-1})$ and $B(z^{-1})$ are polynomials in the shift operator $z^{-1}$ with:

$$\begin{array}{l}A(z^{-1})=1+a_1{z}^{-1}+a_2{z}^{-2}+\dots +a_n{z}^{-n},\\ B(z^{-1})=b_1{z}^{-1}+b_2{z}^{-2}+\dots +b_n{z}^{-n},b_i\in R^{1\times n_u},i=1,2,\dots ,n,\\ {\overline{u}}_t=\left[\begin{array}{c}{f}_1\left(u_t(1)\right)\\ f_2\left(u_2(2)\right)\\ ⋮\\ f_{n_u}\left(u_t(n_u)\right)\end{array}\right]=\left[\begin{array}{c}{\overline{u}}_t\left(1\right)\\ {\overline{u}}_t\left(2\right)\\ ⋮\\ {\overline{u}}_t\left(n_u\right)\end{array}\right]=\left[\begin{array}{c}{c}_{11}{u}_t(1)+c_{21}{u}_t{(1)}^2+\dots +c_{q1}{u}_t{(1)}^q\\ c_{12}{u}_t(2)+c_{22}{u}_t{(2)}^2+\dots +c_{q2}{u}_t{(2)}^q\\ ⋮\\ c_{1n_u}{u}_t(n_u)+c_{2n_u}{u}_t{(n_u)}^2+\dots +c_{q{n}_u}{u}_t{(n_u)}^q\end{array}\right].\end{array}$$

The following assumptions are made:

Assumption 1. Suppose the order of polynomial $A\left(z^{-1}\right)$ and $B\left(z^{-1}\right)$ are equal to $n$. It is assumed that $n$ and $q$ are both known, or can be determined by trial and error.

Assumption 2. The first coefficient of the function $f_j,\left(j=1,2,\cdots ,n_u\right)$ equals 1; i.e., $c_{1j}=1,\left(j=1,2,\cdots ,n_u\right)$ [24].

Assumption 3. The order of the function $f_j,\left(j=1,2,\cdots ,n_u\right)$ is $q$. Suppose $q$ is even, then there is at least one real root of (20). In this way, one is guaranteed that a value of $u_t$ can be found [30].

Then (2) can be rewritten as follows:

$$y_t=-{\displaystyle \sum _{i=1}^n{a}_i{y}_{t-i}}+{\displaystyle \sum _{i=1}^n{b}_i{\overline{u}}_{t-i}}=-{\displaystyle \sum _{i=1}^n{a}_i{y}_{t-i}}+{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^{n_u}{\displaystyle \sum _{k=1}^q{b}_{ij}{c}_{kj}}}}{u}_{t-i}{(j)}^k$$

(3)

where $b_{ij}$ represents the element of the $j$-th column of $b_i$.

Define the following vectors:

$$\begin{array}{l}{\theta }^T=\left[\begin{array}{ccccc}{m}_0& m_1& m_2& \cdots & m_n\end{array}\right]\in R^{1\times n(1+n_{u}q)},\\ m_0^T=\left[\begin{array}{cccc}{a}_1& a_2& \cdots & a_n\end{array}\right],\\ m_i^T=\left[\begin{array}{ccccccccc}{c}_{11}{b}_{i1}& c_{12}{b}_{i2}& \cdots & c_{1n_u}{b}_{i{n}_u}& \cdots & c_{q1}{b}_{i1}& c_{q2}{b}_{i2}& \cdots & c_{q{n}_u}{b}_{i{n}_u}\end{array}\right]\in R^{1\times n_{u}q},\\ {\psi }_t^T=\left[\begin{array}{ccccc}{\psi }_t(0)& {\psi }_t(1)& {\psi }_t(2)& \cdots & {\psi }_t(n)\end{array}\right],\\ {\psi }_t^T(0)=\left[\begin{array}{cccc}-y_{t-1}& -y_{t-2}& \cdots & -y_{t-n}\end{array}\right],\\ {\psi }_t^T(i)=\left[\begin{array}{ccccccccc}{u}_{t-i}(1)& u_{t-i}(2)& \cdots & u_{t-i}(n_u)& \cdots & u_{t-i}{(1)}^q& u_{t-i}{(2)}^q& \cdots & u_{t-i}{(n_u)}^q\end{array}\right],\end{array}$$

(4)

Then (3) can be rewritten into the following form:

$$y_t={\psi }_t{}^{\mathrm{T}}\theta$$

(5)

and the recursive least squares method is used to identify parameters $\theta$:

$$\begin{array}{l}{K}_m=P_{m-1}{\psi }_m{[1+{\psi }_m{}^{\mathrm{T}}{P}_{m-1}{\psi }_m]}^{-1},\\ {\widehat{\theta }}_m={\widehat{\theta }}_{m-1}+K_m[y_m-{\psi }_m{}^{\mathrm{T}}{\widehat{\theta }}_{m-1}],\\ P_m=P_{m-1}-K_m{\psi }_m{}^{\mathrm{T}}{P}_{m-1}\end{array}$$

(6)

where $m$ represents the current moment, ${\widehat{\theta }}_{m-1}$ represents the parameter estimation at the previous moment.

After obtaining the estimated value $\widehat{\theta }$ of the parameter $\theta$, the parameter $a_i(i=1,2,\cdots ,n)$ can be obtained directly from $\widehat{\theta }$, but the parameter $b_i(i=1,2,\cdots ,n)$ also needs to be separated from the parameter $c_{ij}$. Without loss of generality, let’s assume that $c_{1j}=1,(j=1,2,\cdots ,n_u)$, then $b_i(i=1,2,\cdots ,n)$ can be obtained directly from the first to $n_u$ lines of $m_i$.

After obtaining the parameter $b_i(i=1,2,\cdots ,n)$, the remaining parameter $c_{kj}(j=1,2,\cdots ,n_u,k=2,3,\cdots ,q)$ of the nonlinear function is obtained by:

$$c_{kj}=\frac{1}{n}{\displaystyle \sum _{i=1}^n\frac{m_i((k-1)n_u+j,:)}{b_{ij}}}$$

(7)

where $m_i((k-1)n_u+j,:)$ represents the $[(k-1)n_u+j]$-th row element of $m_i$.

In summary, the linear and nonlinear coefficients of the Hammerstein model (2) of the MISO system can be obtained through Equations (6) and (7), and the obtained Hammerstein model of the cooling system of the MGT-CCP will be shown and verified in Section 5.

4. Coordinated Predictive Controller Design

In this section, first we introduce the GPC algorithm, then propose the Hammerstein-GPC algorithm, and finally design the coordinated predictive control system for the MGT-CCP system.

4.1. Generalized Predictive Control

Consider the following controlled auto-regressive integrated moving average (CARIMA) model:

$$A(z^{-1})y_t=B(z^{-1})u_{t-1}+e_t/\Delta$$

(8)

where $t$ represents the discrete time, $z^{-1}$ represents the post-shift operator, $\Delta =1-z^{-1}$ represents the difference operator, $e_t$ represents noise, and $A(z^{-1})$ and $B(z^{-1})$ are polynomials in the shift operator $z^{-1}$ with $A(z^{-1})=1+a_1{z}^{-1}+a_2{z}^{-2}+\dots +a_n{z}^{-n}$, $B(z^{-1})=b_0+b_1{z}^{-1}+b_2{z}^{-2}+\dots b_{n_b}{z}^{-n_b}$.

Ignoring the noise term $e_t$, the Equation (8) can be written as follows:

$$y_t+a_1{y}_{t-1}+\cdots +a_n{y}_{t-n}=b_0{u}_{t-1}+\cdots +b_{n_b}{u}_{t-n_b-1}$$

(9)

which then can be expressed as

$$y_{t+j}=-a_1{y}_{t+j-1}-\cdots -a_n{y}_{t+j-n}+b_0{u}_{t+j-1}+\cdots +b_{n_b}{u}_{t+j-n_b-1}$$

(10)

Consider the following Diophantine equation:

$$1=E_j(z)A\Delta +z^{-1}{F}_j(z)$$

(11)

where $E_j(z)$ and $F_j(z)$ are polynomials in the shift operator $z^{-1}$ with $E_j(z^{-1})=e_{j,0}+e_{j,1}{z}^{-1}+\cdots +e_{j,j-1}{z}^{-(j-1)}$, $F_j(z^{-1})=f_{j,0}+f_{j,1}{z}^{-1}+\cdots +f_{j,n}{z}^{-n}$.

Combining the Equations (8) and (11), through some mathematical transformations, we can get the system output prediction value ${\widehat{y}}_{t+j}$ at time $t+j$:

$${\widehat{y}}_{t+j}=E_{j}B\Delta u_{t+j-1}+F_j{y}_t$$

(12)

In GPC algorithm, the optimization performance indicator at time $t$ takes the following form:

$$\min J(t)=q{\displaystyle \sum _{j=1}^P({\omega }_{t+j}-}{\widehat{y}}_{t+j}{)}^2+r{\displaystyle \sum _{j=1}^M\Delta u_{t+j-1}^2}$$

(13)

where $P$ represents prediction horizon, $M$ represents control horizon, $q$ represents the error weight coefficient, $r$ represents the control weight coefficient, and ${\omega }_{t+j}$ represents the expected reference value of the system output.

Define polynomial $G_j$:

$$G_j=E_{j}B=g_{j,0}+g_{j,1}{z}^{-1}+\cdots +g_{j,n_b+j-1}{z}^{-(n_b+j-1)}$$

(14)

Then (12) can be rewritten as follows:

$${\widehat{y}}_{t+j}=G_j\Delta u_{t+j-1}+F_j{y}_t=g_{j,0}\Delta u_{t+j-1}+g_{j,1}\Delta u_{t+j-2}+\cdots +g_{j,j-1}\Delta u_t+l_j(t)$$

(15)

where $l_j(t)=(G_j-g_{j,0}-g_{j,1}{z}^{-1}-\cdots -g_{j,j-1}{z}^{-(j-1)})\Delta u_{t+j-1}+F_j{y}_t$.

According to (15) we can easily get the following formula:

$$Y=G\Delta U+L$$

(16)

where:

$$Y=\left[\begin{array}{c}{\widehat{y}}_{t+1}\\ ⋮\\ {\widehat{y}}_{t+P}\end{array}\right],\Delta U=\left[\begin{array}{c}\Delta u_t\\ ⋮\\ \Delta u_{t+M-1}\end{array}\right],L=\left[\begin{array}{c}{l}_1(t)\\ ⋮\\ l_P(t)\end{array}\right],G=\left[\begin{array}{ccccc}{g}_{1,0}& \cdots & g_{1,0}& 0& 0\\ ⋮& & & \ddots & 0\\ g_{M,M-1}& & \cdots & & g_{1,0}\\ ⋮& & & & ⋮\\ g_{P,P-1}& & \cdots & & g_{P,P-M}\end{array}\right].$$

Because of the influence of model error and disturbance, the output prediction value of the system needs to be corrected by the actual output error on the basis of the output of the prediction model:

$$Y=G\Delta U+L+h{e}_t=G\Delta U+L+h\left(y_t-{\widehat{y}}_{t|t-1}\right)$$

(17)

where $e_t$ is the prediction error, $y_t$ is the actual output of the object, and ${\widehat{y}}_{t|t-1}$ is the predicted output of the object. According to (13) and (17), the following can be obtained:

$$\min J(t)={\left(W-G\Delta U-L-h{e}_t\right)}^{T}Q\left(W-G\Delta U-L-h{e}_t\right)+\Delta U^{T}R\Delta U$$

(18)

where $W={\left[\begin{array}{ccc}{\omega }_{t+1}& \cdots & {\omega }_{t+P}\end{array}\right]}^T$, $Q=q\cdot I_{P\times P}$ represents the error weighting matrix, $R=r\cdot I_{M\times M}$ represents the control weight matrix, and $I$ represents the unit matrix.

The output variable of the GPC algorithm is obtained according to (18):

$$u_t=u_{t-1}+\Delta u_t=u_{t-1}+\left[\begin{array}{cccc}1& 0& \cdots & 0\end{array}\right]{(G^{\mathrm{T}}QG+R)}^{-1}{G}^{\mathrm{T}}(W-L-h{e}_t)$$

(19)

4.2. Hammerstein Model-Based Generalized Predictive Control

Consider the SISO Hammerstein system as follows:

$${\overline{u}}_t=f\left(u_t\right)=c_1{u}_t+c_2{u}_t{}^2+\dots +c_q{u}_t{}^q$$

(20)

$$y_t={\displaystyle \sum _{i=1}^{n_a}{a}_i}{y}_{t-i}+{\displaystyle \sum _{j=1}^{n_b}{b}_j}{\overline{u}}_{t-j}$$

(21)

where (20) is a static nonlinear element without memory, (21) is a dynamic linear model, $n_a$ and $n_b$ are the order of the linear model, $a_i$ and $b_j$ are the parameters of the linear model, $u$ and $y$ are input and output of the system, respectively, and ${\overline{u}}_t$ is the intermediate variable.

The system model given by (20) and (21) shows that the control input $u$ affects the system output $y$ through the intermediate variable ${\overline{u}}_t$. Although the relationship between $u$ and $y$ is nonlinear, the relationship between ${\overline{u}}_t$ and $y$ is linear. Therefore, the GPC algorithm is designed with the linear model (21) first, i.e., the desired intermediate variable ${\overline{u}}_t$ is calculated through the receding-horizon optimization of the GPC algorithm. Then the control action $u_t$ is solved from nonlinear algebraic Equation (20). The structure of the Hammerstein-GPC controller is shown in Figure 4. This control strategy makes full use of the special structure of the Hammerstein model and is actually a nonlinear separation strategy, which makes the whole solution simple.

4.3. Coordinated Predictive Controller Design for the MGT-CCP System

This section presents the coordinated predictive control strategy for the MGT-CCP system, as shown in Figure 5.

Considering the physical characteristics of the MGT-CCP system, the MGT-CCP system was split into a power supply system and a cooling system. In view of the degree of nonlinearity of the two subsystems, the linear GPC is designed for the power generation system and the nonlinear Hammerstein-GPC is designed for the cooling system. The nonlinearity of the cooling system is effectively suppressed by nonlinear compensation, and the coupling between subsystems is effectively weakened by feedforward compensation. Besides, due to the inevitable modeling error, feedback correction is adopted in the proposed control strategy.

The main calculation process of the GPC-HGPC coordinated control algorithm is summarized as follows:

(1)

At time $t$, the calculation of the linear GPC controller is performed, and the optimal input increment sequence $\Delta U(1)={\left[\begin{array}{cccc}\Delta u_t(1)& \Delta u_{t+1}(1)& \cdots & \Delta u_{t+M-1}(1)\end{array}\right]}^T$ of the fuel amount is calculated, and the optimal control variable $u_t(1)=u_{t-1}(1)+\Delta u_t(1)$ of the MGT system is further obtained;

(2)

The intermediate control variable sequence $\Delta \overline{U}(1)={\left[\begin{array}{cccc}\Delta {\overline{u}}_t(1)& \Delta {\overline{u}}_{t+1}(1)& \cdots & \Delta {\overline{u}}_{t+M-1}(1)\end{array}\right]}^T$ is calculated by using the fuel quantity optimal input increment sequence $\Delta U(1)={\left[\begin{array}{cccc}\Delta u_t(1)& \Delta u_{t+1}(1)& \cdots & \Delta u_{t+M-1}(1)\end{array}\right]}^T$ and the nonlinear polynomial $f_1(\cdot )$;

(3)

Adding the sequence $\Delta \overline{U}(1)={\left[\begin{array}{cccc}\Delta {\overline{u}}_t(1)& \Delta {\overline{u}}_{t+1}(1)& \cdots & \Delta {\overline{u}}_{t+M-1}(1)\end{array}\right]}^T$ as feedforward information to Equation (15) gives an output prediction of chilled water temperature:

$$\begin{array}{cc}{\widehat{y}}_{t+j}(2)& =g_{j,0}^{\left(2\right)}\Delta u_{t+j-1}(2)+g_{j,1}^{\left(2\right)}\Delta u_{t+j-2}(2)+\cdots +g_{j,j-1}^{\left(2\right)}\Delta u_t(2)+l_j(t)\hfill \\ & +g_{j,0}^{\left(1\right)}\Delta u_{t+j-1}(1)+g_{j,1}^{\left(1\right)}\Delta u_{t+j-2}(1)+\cdots +g_{j,j-1}^{\left(1\right)}\Delta u_t(1)\hfill \\ & =g_{j,0}^{\left(2\right)}\Delta u_{t+j-1}(2)+g_{j,1}^{\left(2\right)}\Delta u_{t+j-2}(2)+\cdots +g_{j,j-1}^{\left(2\right)}\Delta u_t(2)+{\stackrel{⌢}{l}}_j(t)\hfill \end{array}$$

(22)

where ${\stackrel{⌢}{l}}_j(t)=g_{j,0}^{\left(1\right)}\Delta u_{t+j-1}(1)+g_{j,1}^{\left(1\right)}\Delta u_{t+j-2}(1)+\cdots +g_{j,j-1}^{\left(1\right)}\Delta u_t(1)+l_j(t)$. Note that ${\stackrel{⌢}{l}}_j(t)$ is known at time $t$;

(4)

Based on (22), using the linear GPC algorithm to calculate the intermediate control variable increment sequence $\Delta \overline{U}(2)={\left[\begin{array}{cccc}\Delta {\overline{u}}_t(2)& \Delta {\overline{u}}_{t+1}(2)& \cdots & \Delta {\overline{u}}_{t+M-1}(2)\end{array}\right]}^T$, and then obtain ${\overline{u}}_t(2)={\overline{u}}_{t-1}(2)+\Delta {\overline{u}}_t(2)$;

(5)

Using the intermediate control variable ${\overline{u}}_t(2)$ and the inverse of the nonlinear polynomial $f_2(\cdot )$ to obtain the refrigerant valve opening $u_t(2)$;

(6)

Next sampling time, return to step (1), repeat steps (1) to (5).

Remark 2: For the sake of simplicity, the proposed method uses the control model obtained by offline identification. If the actual MGT-CCP system has large time-varying characteristics, the performance of the proposed control algorithm will be attenuated. But this problem can be overcome by adding online identification of the control model.

Remark 3: This algorithm can only eliminate the influence of external disturbance of the step type. For time-varying disturbances, adding on-line identification of the disturbance model to the proposed algorithm can improve control performance [31].

5. Simulation Results

This section verifies the effectiveness of the multivariate Hammerstein model identification strategy and the coordinated predictive controller designed for the MGT-CCP system. First, the accuracy of the multivariate Hammerstein model was verified, and then the proposed coordinated predictive controller was tested and compared to the PI controller.

5.1. Verification of Identification Model

5.1.1. Power Generation System

The power generation system is affected only by the fuel, so it can be considered as a SISO system. According to the nonlinear analysis in Section 2, the nonlinearity of the power generation system is weak, so the conventional recursive least squares method is used to identify the model of the power generation system. The input sequence shown in Figure 6a is applied to the MGT-CCP system to obtain the output data shown in Figure 6b. The model of the power generation system was identified using the first 3000 s of data. The identified input and output models are as follows:

$$y_t-1.5691y_{t-1}+0.5895y_{t-2}=176940u_{t-1}-42160u_{t-2}$$

(23)

The accuracy of the identification model is verified using the last 3000 s of the data of Figure 6. The verification results are shown in Figure 7. It can be seen that the identified linear model has high accuracy, which also indicates that the nonlinearity of the power supply system is very weak, which is consistent with the nonlinear analysis results in the second section.

5.1.2. Cooling System

The cooling system is affected by both the fuel and the refrigerant valve opening, so the cooling system can be viewed as a two-input single-output system. According to the nonlinear analysis in Section 2, the nonlinearity of the cooling system is very strong, so the Hammerstein model identification method proposed in Section 3 is used to identify the model of the cooling system. The input sequence shown in Figure 8a is applied to the MGT-CCP system to obtain the output data shown in Figure 8b. The model of the power supply system was identified using the first 3000 s of data. The identified Hammerstein model is as follows:

$$\begin{array}{l}{y}_t-2.003y_{t-1}+0.9596y_{t-2}+0.1375y_{t-3}-0.09403y_{t-4}=[0.2851-1.846\times {10}^{-5}]{\overline{u}}_{t-1}\\ +[0.0185-3.288\times {10}^{-5}]{\overline{u}}_{t-2}+[-0.056765.084\times {10}^{-5}]{\overline{u}}_{t-3}+[-0.30847.234\times {10}^{-6}]{\overline{u}}_{t-4}.\end{array}$$

(24)

where $\overline{u}=\left[\begin{array}{c}{f}_1\left(u_1\right)\\ f_2\left(u_2\right)\end{array}\right]=\left[\begin{array}{c}{u}_1-0.0019u_1{}^2\\ u_2-0.0108u_2{}^2\end{array}\right]$.

The accuracy of the identification model is verified using the last 3000 s of the data of Figure 8. The verification results are shown in Figure 9. It can be seen that the identified Hammerstein model has high accuracy and can reflect the changing trend of the system.

5.2. Testing of the Proposed Control Strategy

This section tests the proposed GPC-HGPC coordinated controller for the MGT-CCP system. The parameters of the proposed controller are set as follows: for linear GPC controllers, prediction horizon $N_1=20$, control horizon $M_1=10$, error weight coefficient $q_1=1$, control weight coefficient $r_1=0.8$; for the Hammerstein-GPC controller, prediction horizon $N_2=600$, control horizon $M_2=20$, error weight coefficient $q_2=1$, control weight coefficient $r_2=0.1$.

The PID controller and MPC control are compared with the proposed controller. The parameters of the PID controller are as follows: for the power supply system, the proportional gain $k_{\mathrm{p}}=8\times {10}^{-8}$, the integral gain $k_{\mathrm{i}}=8.17\times {10}^{-9}$, the differential gain $k_{\mathrm{d}}=0$; for the cooling system, the proportional gain $k_{\mathrm{p}}=197.57$, the integral gain $k_{\mathrm{i}}=1.45$, the differential gain $k_{\mathrm{d}}=1183.39$. Note that the PID controller is optimally tuned using the Tuner function of the PID module provided by Simulink.

The parameters of the MPC controller are as follows: prediction horizon $N=39$, control horizon $M_1=5$, sampling time $t_s=20s$, control weight coefficient $Q_u=diag\left(0.1377,0.1377\right)$, control rate weight coefficient $Q_{\Delta u}=diag\left(0.0726,0.0726\right)$, output weight coefficient $Q_y=diag\left(0.0689,38.5596\right)$.

The control structure of the comparison methods is shown in Figure 10.

In order to quantitatively compare the load tracking performance of different controllers, the integrated absolute error (IAE) is used in the following simulations. IAE is given by:

$$\mathrm{IAE}={\displaystyle {\int }_0^{\infty }\left|ref\left(t\right)-y\left(t\right)\right|}dt$$

(25)

5.2.1. Case 1: Simultaneous Cooling and Power Load-Following Performance

Case 1 was designed to test the simultaneous cooling and power load-following performance of controllers. First, compare the proposed method with conventional MPC. Assume that the MGT-CCP system is operating at point #3. At 200 s, the MGT-CCP system changes from the operating point #3 ($73.5\mathrm{kW},8.5{}^{°}\mathrm{C}$) to the operating point #4 ($69.69\mathrm{kW},8.8{}^{°}\mathrm{C}$). The simulation results are shown in Figure 11 and

Table 2. Performance indices for Figure 11.

Method	$\mathbf{I}\mathbf{A}{\mathbf{E}}_{\mathbf{Power}}$	$\mathbf{IA}{\mathbf{E}}_{\mathbf{Temperature}}$
MPC	907	607
Proposed method	841	121

.

Both the proposed method and the conventional MPC can track the set points. But the proposed method is faster for load tracking. The power tracking performance under the IAE indicator increased by 7%, and the cold load tracking performance increased by 80%. The reason is that the proposed method uses coupling compensation and nonlinear compensation strategies.

Then we test the tracking performance of the controllers for a wide range of load changes. At 1500 s, the MGT-CCP system changes from the operating point #4 ($69.69\mathrm{kw},8.8{}^{°}\mathrm{C}$) to the operating point #5 ($65.18\mathrm{kw},9.1{}^{°}\mathrm{C}$), and then keeps running at the operating point #4, at 3000 s, the MGT-CCP system changes from operating point #5 to operating point #1 ($80\mathrm{kw},8{}^{°}\mathrm{C}$), and finally keeps running at operating point #1. The simulation results are shown in Figure 12 and

Table 3. Performance indices for Figure 12.

Method	${\mathbf{IAE}}_{\mathbf{Power}}$	${\mathbf{IAE}}_{\mathbf{Temperature}}$
PID	1636	364
Proposed method	846	292

. The simulation results show that the proposed method has smaller overshoot and faster tracking speed. The power tracking performance under the IAE indicator increased by 48%, and the cold load tracking performance increased by 20%.

5.2.2. Case 2: Performance of Unknown Input Disturbances Rejection

Case 2 is used to test the controller’s ability to suppress unknown input disturbances. Since the MGT-CCP system operates in a high temperature environment, the wear caused by the frequent operation of the valve will affect its stability. Therefore, the impact of input disturbances on the MGT-CCP system cannot be ignored.

Assume that the MGT-CCP system is running at working point #3. The input side of the system is subject to unknown input disturbances as shown in Figure 13, and the simulation results are shown in Figure 14. The simulation results show that the controllers can eliminate the influence of unknown input disturbances, but the proposed method suppresses the disturbance faster and the deviation from the operating point is smaller.

5.2.3. Case 3: Tracking Performance on Time-Varying Power Load Demands

Renewable energy sources such as wind power and photovoltaic power have time-varying characteristics. When they are connected to a micro-grid with a MGT-CCP system, it is necessary to adjust the power output of gas turbine to maintain the stability of the micro-grid. Therefore, Case 3 tests the controller’s ability to track time-varying power load demands. Figure 15 is a block diagram of the MGT-CCP system connected to renewable energy.

Assuming the MGT-CCP system is operating at operating point #3, the power required by the user is a fixed value and the power from the renewable sources is shown in Figure 16. In order to verify the effectiveness of the proposed method in depth, in Figure 16, the power changes rapidly in the first 3000 s, and the power changes slowly in the next 3000 s. The results are given in Figure 17 and

Table 4. Performance indices for Figure 17.

Method	${\mathbf{IAE}}_{\mathbf{Sum}\mathbf{of}\mathbf{power}\mathbf{output}}$	${\mathbf{IAE}}_{\mathbf{Temperature}}$
PID	4053	22
Proposed method	1525	20

. The simulation results show that the proposed method can track the electric load command better, the tracking speed is faster and the dynamic deviation is smaller. In addition, because of the implementation of coordination, the cold-water temperature deviates from the set value by a smaller amount. The power tracking performance under the IAE indicator increased by 62%, and the cold load tracking performance increased by 9%.

6. Conclusions

In order to achieve a sustainable future with more renewable energy integrated into distributed energy systems, MGT-CCP systems must frequently change loads to achieve grid balance. This requires more flexibility in the operation of the MGT-CCP. To provide guidance for controller design, this paper first studies the dynamic behavior and nonlinearity distribution of the MGT-CCP process. The results show that the MGT speed is affected only by the fuel, and the nonlinearity of the power generation subsystem is extremely weak. However, the cold-water temperature is affected by both the fuel and the refrigerant valve opening, and the cooling subsystem is extremely nonlinear. Therefore, based on these results, a coordinated predictive control strategy is proposed to simultaneously improve cooling-and-power load-following capability. The simulation results show that compared with PID and MPC, the proposed control method not only can improve cooling and power load-following capability at the same time, but also has the best control effect when accessing with renewable energy. In the future, we will strive to consider economics into the controller design of the MGT-CCP system, which is also an important issue for implementation of predictive controls.