A novel robust MPC scheme established on LMI formulation for surge instability of uncertain compressor system with actuator constraint and piping acoustic

Abstract

This paper presents a new control method for surge instability of the compressor system. Due to the importance of the active control approach in improving the efficiency of the compressor system, a new scheme based on the model predictive control (MPC) technique has been considered to control surge which gives the benefits of control signal optimization by considering the constraints on states and input. Because of designing a novel MPC by using of linear matrix inequality scheme, the optimization problem is solved in less time and with less complexity. The proposed method is able to consider the limitation of the close coupled valve actuator and also covers the uncertainty on the model. In addition, the developed model describes the nonlinear behaviour of the compressor system and handles the effects of the pipe on surge instability. Using Lyapunov theory, the stability of the closed-loop system under the proposed robust active controller is shown. The simulation results show the high pressure and high efficiency operation of the system under study in different operating conditions.

KEYWORDS:

• Centrifugal compressor
• surge
• piping acoustic
• linear matrix inequality
• uncertainty
• model predictive control
• close coupled valve

Nomenclature

MPC	=	Model predictive control
NMPC	=	Nonlinear model predictive control
LMI MPC	=	Linear matrix inequality model predictive control
CCV	=	Close coupled valve
${\mathrm{\Phi }}_c$	=	Compressor mass flow
$\mathrm{\Psi }c$	=	Compressor pressure rise
${\mathrm{\Phi }}_p$	=	Plenum mass flow
${\mathrm{\Psi }}_p$	=	Plenum pressure rise
${\mathrm{\Phi }}_r$	=	Close coupled valve mass flow
${\mathrm{\Phi }}_{\mathrm{th}}$	=	Throttle valve mass flow
${\mathrm{\Psi }}_{\mathrm{th}}$	=	Throttle valve pressure rise
$U$	=	Impeller velocity
$\rho$	=	Air density (kg/m3)
$A_c$	=	The cross sectional area of the compressor duct
$P_{o1}$	=	The inlet absolute pressure
${\rho }_{o1}$	=	The density at the inlet
${\omega }_H$	=	Helmholtz frequency
$B$	=	The Greitzer stability parameter
$c_r$	=	Valve constant
$u_r$	=	The ccv opening percentage
$A_1$,$B_1$ and $D_1$	=	The characteristic curve correspond to the stable flow region of the compressor
$A_2$, $B_2$ and $D_2$	=	The characteristic curve correspond to t unstable flow region
$P_P$ and $Q_P$	=	Upstream (plenum) pressure and volumetric flow rate
$P_{\mathrm{th}}$ and $Q_{\mathrm{th}}$	=	Downstream (throttle duct) pressure and flow rate.
$u_{\mathrm{th}}$	=	The throttle valve opening percentage
$c_{\mathrm{th}}$	=	The valve constant

1. Introduction

Centrifugal compressors are used in turbochargers, petrochemicals, aerospace, power plants and turbo shaft gas turbine engines due to their high pressure ratio, extensive operational range and cost benefits. They are usually engine or turbine driven and their widespread use is mostly due to smooth operation, high tolerance of process fluctuations and higher reliability compared to other types of compressors. These devices are an integral part of the industry, and any malfunction in their operation can cause upstream and downstream equipment failure. Due to the special place of compressors in industrial energy consumption, many efforts are being made to ensure the stable and optimal performance of these devices simultaneously [1–3]. The performance of the compression system depends on the performance curve of compressor and its connected downstream and upstream equipment.

Various limitations affect the performance of the compressor system, the most important of which are: surge, stall and choke. Despite the negative effects of all these instabilities, surge, as the most destructive phenomenon, can cause severe mechanical damage to the compressor system and take it out of the circuit, and is therefore of particular importance. Compressor surge is a type of aerodynamic instability in axial or centrifugal compressors. The term refers to the intense airflow that oscillates in the axial direction of the compressor, indicating that the axial component of the fluid velocity changes intermittently and may even become negative. This is a fundamental dynamic instability that affects the all-inclusive compressor system and can cause high thermal and mechanical stresses in the system and cause severe damage to the compressor and its accessories. In addition, surge fluctuations can lead to unacceptable noise levels. Therefore, designing a suitable surge controller is essential for the compressor system.

Having a real model of the compressor system, considering the effect of different components on its dynamic behaviour, is the most important step in scheming a suitable controller. Over the years, several mathematical models have been proposed to describe the flow dynamics in compressor systems. Greitzer presented a dynamic model for the compressor system, compared its behaviour with real system, and evaluated the effect of different model parameters [4,5]. After Greitzer's worthy effort, the model and physics in [6] concluded that the geometric accuracy and characteristics of the compressor are the two main issues in model development for the main compressor system. To evaluate the influence of different parts on surge oscillation, the effect of pipe on the stability of the compressor system was investigated in [7,8]. Spark and Brown showed that the pipe system connected to the compressor has a very important role in the stability of the compressor system and can lead to strengthening and accelerating the occurrence of surge and, of course, limiting the stability range of the compressor. Accordingly, this paper presents a control-oriented model for the compressor system that, along with a description of the system state space, considers the effects of piping and the CCV actuator.

There are two main categories for compressor system surge controllers, which include passive and active control approaches [9,10]. Extensive efforts to increase efficiency, reliability and performance range have led to significant advances for controlling the surge instability of the compressor system. Passive control methods or surge avoidance methods, despite their high reliability and widespread use in existing compressor systems, are no longer the preferred method in designing new compressor systems, as they reduce system efficiency and sacrifice operating range. The active control method, despite reducing reliability, extends the compressor operating range to less mass flows. The use of this scheme has been highly regarded due to the importance of energy optimization and increasing efficiency, and various actuators have been used in this regard. Active surge control makes the compressor system much more flexible and allows the machine to perform more efficiently, usually near the surge line.

Among the various actuators used in the design of surge controllers, the close coupled valve is the most common actuator used in the compressor systems of the oil and gas industry due to its high speed and suitable performance [11,12]. Another aim of this study is to cover the physical limitations governing CCV actuator, while most active control techniques introduced to date [13,14] have not considered these important and fundamental restrictions. Given that there is a certain range for flow and pressure of different parts of the compressor system, considering the constraints in the states is another goal in the design of the controller. Also, the next goal is to consider the uncertainty that results from the uncertainty in the percentage of throttle valve opening and of course the uncertainty in modelling.

To the best of the authors’ knowledge to date, only a handful of surge controllers have been designed with pipe effects in mind. The controller designed in [15] is not able to capture the completely nonlinear dynamics of the compressor system and another disadvantage is the low operating range. Also, the control method designed in [16] despite the numerous advantages, suffers from a high computational volume and, of course, is not able to guarantee stability in the presence of uncertainty. In view of the above and considering the ability and application of predictive and robust methods in different applications [17–20], this paper presents a novel robust model predictive control approach by using of LMI to control the surge instability of the compressor system in the presence of uncertainty and includes the following innovations: (1) Considering the nonlinear model of compressor system. Due to the various components affecting the behaviour of the compressor, it is necessary to consider its nonlinear model; (2) Considering the effects of pipe on surge instability. The downstream pipe connected to the compressor has a significant effect on the operating range of the compressor and the displacement of the surge line; (3) Providing robust control scheme based on LMI. The LMI-based model predictive control technique, in addition to reducing computational volume, offers the ability to online executing simultaneously with optimization; (4) Ensuring stability in the presence of model and throttle valve uncertainties; The proposed technique is able to ensure stability and surge-free operation in the presence of uncertainty in the model and in addition has the ability to deal with changes in the throttle valve; (5) Considering the limitations of the CCV actuator and system states; An important advantage of the planned scheme is the ability to satisfy the constraints on the compressor system states and to fulfil the actuator limitations.

Therefore, this paper is organized as follows: The second part considers the dynamic model of a compressor system consisting of several compressors and considering their interactions with each other, offers the state space model that includes the effects of piping and nonlinear terms. In the third part, for the working class of the compressor system, a new MPC controller based on LMI is presented and then the stability of the closed loop system is guaranteed under the proposed method. In the fourth section, the simulation results are presented and the results are compared and evaluated. Finally, the fifth section includes conclusions and suggestions for future studies.

2. The model of compressor system

This research focuses on multi-compressor compression systems. To achieve the required pressure increase in a typical gas processing system in Iran, the fluid first enters the compression system duct. The compression system uses a bank of constant speed compressors to raise the fluid pressure to the required level and deliver it directly to the transmission areas. Figure 1 shows a schematic of a sample system with three compressors, including two parallel compressors and a series compressor with the same specifications. According to [11] in the compression system, the recycle valve is considered as the actuator for each compressor. The fluid is amplified in two parallel compressors and then sent to the third compressor. For more reliability, the Spill Back is envisaged as a safety valve for the system that causes the compressors to interact.

Figure 1. Schematic of a sample system with three compressors.

Figure 2 shows a block diagram of the various components of each compressor, while its corresponding dynamics are fully described in reference [15]. According to the characteristics of the Spill Back valve, the correlation between the flow of the Spill Back valve ${\mathrm{\Phi }}_{\mathrm{SB}}$ and the pressure of the third compressor throttle (${{\mathrm{\Psi }}_{\mathrm{th}}}_3$) in the compressor system of Figure 1 is as follows. (1) $\begin{array}{r}{\mathrm{\Phi }}_{\mathrm{SB}}={\gamma }_{\mathrm{SB}}\sqrt{{{\mathrm{\Psi }}_{\mathrm{th}}}_3}\end{array}$(1) Where, ${\gamma }_{\mathrm{SB}}$ indicates the valve’s gain.

Figure 2. Block diagram of each compressor with pipeline dynamics.

In the resulting nonlinear system, ${\mathrm{\Phi }}_c$ indicates the compressor mass flow,${\mathrm{\Psi }}_p$ specifies the plenum pressure rise, ${\mathrm{\Psi }}_{\mathrm{th}}$ states the throttle section pressure rise and ${\mathrm{\Phi }}_p$ shows the plenum mass flow rate. As stated by Figure 1, the state space equations for the compression system with two parallel compressors and one serial compressor are like this: (2) $\begin{array}{rl}{\left[\begin{array}{c}{\dot{\mathrm{\Phi }}}_c\\ {\dot{\mathrm{\Psi }}}_p\\ {\dot{\mathrm{\Psi }}}_{\mathrm{th}}\\ {\dot{\mathrm{\Phi }}}_p\end{array}\right]}_i& ={\left[\begin{array}{cccc}0& -B{\omega }_H& 0& 0\\ {\displaystyle \frac{{\omega }_H}{B}}& 0& 0& -{\displaystyle \frac{{\omega }_H}{B}}\\ 0& 0& 0& {\displaystyle \frac{2A_{12}{A}_c}{\mathrm{\rho U}}}\\ 0& {\displaystyle \frac{B_{21}\mathrm{\rho U}}{2A_c}}& {\displaystyle \frac{A_{21}\mathrm{\rho U}}{2A_c}}& A_{22}\end{array}\right]}_i\\ & \times {\left[\begin{array}{c}{\mathrm{\Phi }}_c\\ {\mathrm{\Psi }}_p\\ {\mathrm{\Psi }}_{\mathrm{th}}\\ {\mathrm{\Phi }}_p\end{array}\right]}_i+{\left[\begin{array}{c}B{\omega }_H(\mathrm{\Psi c},\mathrm{ss})\\ 0\\ 0\\ {\displaystyle \frac{\rho P_{o1}}{P_{o1}U{A}_c}}(A_{21}+B_{21})\end{array}\right]}_i\\ & +{\left[\begin{array}{cc}0& 0\\ -{\displaystyle \frac{{\omega }_H}{B}}& 0\\ 0& {\displaystyle \frac{2B_{12}{A}_c}{\mathrm{\rho U}}}\\ 0& B_{22}\end{array}\right]}_i{\left[\begin{array}{c}{\mathrm{\Phi }}_r\\ {\mathrm{\Phi }}_{\mathrm{th}}\end{array}\right]}_i+g_i(x)\end{array}$(2) Where $i=\mathrm{1.2.3}$ and $g_1=\left[\begin{array}{c}0\\ {\mathrm{\Phi }}_{\mathrm{SB}}\\ 0\\ 0\end{array}\right].g_2=\left[\begin{array}{c}0\\ {\mathrm{\Phi }}_{\mathrm{SB}}\\ 0\\ 0\end{array}\right].g_3=\left[\begin{array}{c}0\\ 0\\ 0\\ 0\end{array}\right]$. The corresponding matrix coefficients of the piping equation in (2) are found to be $A_{12}=3.7\times {10}^6$, $A_{21}=-1.92\times {10}^{-3}$, $A_{22}=-8$, $B_{12}=-3.7\times {10}^6$, $B_{21}=1.92\times {10}^{-3}$ and $B_{22}=7.98$. The parameters of the compressor system model used in this study are given in Table 1.

Table 1. Model parameters for the compression system.

Parameter	Symbol	Unit	Value
Comp. duct length	LC	m	1.86
Comp. duct cross area	AC	m^2	0.0082
Char. curve A1 coeff.	A1	–	−4.717
Char. curve A2 coeff.	A2	–	−172.6
Char. curve B1 coeff.	B1	–	−2.859
Char. curve B2 coeff.	B2	–	36.88
Char. curve D1 coeff.	D1	–	1.193
Char. curve D2 coeff.	D2	–	1.029
GreitzerStab. parameter	B	–	0.44
Helmholtz frequency	wH	Rad/s	80.1
Impeller Tip speed	U	m/s	213.24
Inlet pressure	P01	Pa	101,325
Inlet gas density	p01	Kg/m^3	1.165
Plenum volume	Vp	m^3	0.049
Pipeline length	L	m	6.5
Throttle constant	Cth	–	1.8569
CCV constant	Cr	–	0.5

An Euler approximation of system (2) with sampling time $\mathrm{Tc}$ equals to (3) $\begin{array}{rl}& {\left[\begin{array}{c}{\mathrm{\Phi }}_c\\ {\mathrm{\Psi }}_p\\ {\mathrm{\Psi }}_{\mathrm{th}}\\ {\mathrm{\Phi }}_p\end{array}\right]}_i(\text{k}+1)\\ & =\mathrm{Tc}.{\left[\begin{array}{ccccccccccccc}{\displaystyle \frac{1}{\mathrm{Tc}}}& -B{\omega }_H& 0& 0{\displaystyle \frac{{\omega }_H}{B}}& {\displaystyle \frac{1}{\mathrm{Tc}}}& 0& -{\displaystyle \frac{{\omega }_H}{B}}0& 0& {\displaystyle \frac{1}{\mathrm{Tc}}}& {\displaystyle \frac{2A_{12}{A}_c}{\mathrm{\rho U}}}0& {\displaystyle \frac{B_{21}\mathrm{\rho U}}{2A_c}}& {\displaystyle \frac{A_{21}\mathrm{\rho U}}{2A_c}}& {\displaystyle \frac{1}{\mathrm{Tc}}}+A_{22}\end{array}\right]}_i\\ & \times {\left[\begin{array}{c}{\mathrm{\Phi }}_c{\mathrm{\Psi }}_p{\mathrm{\Psi }}_{\mathrm{th}}{\mathrm{\Phi }}_p\end{array}\right]}_i(\text{k})+\mathrm{Tc}.{\left[\begin{array}{c}B{\omega }_H(\mathrm{\Psi c}.\mathrm{ss})00{\displaystyle \frac{\rho P_{o1}}{P_{o1}U{A}_c}}(A_{21}+B_{21})\end{array}\right]}_i& +\mathrm{Tc}.{\left[\begin{array}{ccccc}0& 0-{\displaystyle \frac{{\omega }_H}{B}}& 00& {\displaystyle \frac{2B_{12}{A}_c}{\mathrm{\rho U}}}0& B_{22}\end{array}\right]}_i{\left[\begin{array}{c}{\mathrm{\Phi }}_r{\mathrm{\Phi }}_{\mathrm{th}}\end{array}\right]}_i(\text{k})+\mathrm{Tc}.g_i(x(k))\end{array}$(3)

Which is a subset of discrete nonlinear systems as follows (4) $\begin{array}{r}{x}_i(k+1)=A_i{x}_i(k)+B_i{u}_i(k)+f_i(x_i(k).u_i(k)).\end{array}$(4)

The next section provides a robust LMI-based MPC for this compressor system. The parameters of the theoretical model can be found in .

3. The proposed robust model predictive control by using of linear matrix inequality

This section presents the robust model predictive control formulation for the work class introduced in the previous section. The basis for ensuring the stability of the proposed controller is based on the fact that first a proof of stability process is presented for one of the agents and the working class under study, and since this process can be generalized to all agents, the stability of the overall closed-loop system in the presence of all agents can be imagined despite all the interactions of agents on each other and other inconveniences.

Consider a discrete-time nonlinear system of one of the agents as (5) $\begin{array}{r}x(k+1)=\mathrm{Ax}(k)+\mathrm{Bu}(k)+f(x(k).u(k)).\end{array}$(5) Where $A$ and $B$ are coefficient matrices, $u(k)\in R^m$ indicates the control input conditional on $|u_r(k)|\le u_{r.\max }.r=\mathrm{1.2.}\dots .m.$ $x(k)\in R^n$ specifies the state of the system bound by $|x_r(k)|\le x_{r.\max }.r=\mathrm{1.2.}\dots .n$ and $f(\cdot .\cdot )$ satisfies a Lipschitz condition and $C^1$ function ($f(0.0)=0$ is the equilibrium of the system). The MPC scheme is used to find state feedback by minimizing an objective function to optimize closed-loop system performance over an optimal time horizon. Based on measurements obtained from the plant, the calculations are repeated at subsequent sampling times. To determine the control law $u(\cdot )$ of system (5), the objective function is defined as follows (6) $\begin{array}{rl}J(k)& =\sum _{i=0}^{\mathrm{\infty }}x(k+i|{k)}^T\mathrm{Qx}(k+i|k)\\ & +u(k+i|k{)}^T\mathrm{Ru}(k+i|k),\end{array}$(6) Where $Q>0.R>0.x(k+i|k)$ signifies state at time $k+i$ founded on the measurements at time $(k.x(k))$ and $u(k+i|k)$ indicates control signal at time $k+i$ calculated by minimizing $J(k)$ at time $k$.

To obtain optimization problem solution for the nonlinear system in (5) through linear matrix inequality, first, the upper limit of $J(k)$ is substituted in Equation (6) to obtain the desired inequality. By supposing a quadratic function $V(x)=x^T\mathrm{Px}$ with $P>0$ and $V(0)=0$, the following inequality holds at sampling time $k$ (7) $\begin{array}{rl}& V(x(k+i+1|k))-V(x(k+i|k))\\ & \le -x(k+i|k{)}^T\mathrm{Qx}(k+i|k)\\ & -u(k+i|k{)}^T\mathrm{Ru}(k+i|k).\end{array}$(7)

It is obtained by summing both sides of (7) from $i=0$ to $i\mathrm{\infty }$ (8) $\begin{array}{r}J(k)\le V(x(k|k))\le \gamma .\end{array}$(8) Where $\gamma$ indicates a positive scalar and it is an upper boundary for (6). Consequently, the robust model predictive control problem is defined in the following.

Theorem 3.1:

Suppose $x(k|k)$ be the measured state $x(k)$ at sample time $k$ in system (5) conditional on $|u_r(k+i|k)|\le u_{r.\max }.i\ge 0.r=\mathrm{1.2.}\dots .m$ and $|x_r(k+i|k)|\le x_{r.\max }.i\ge 0.r=\mathrm{1.2.}\dots .n$. Then, the state feedback $F(k)$ in the control law $u(k+i|k)=F(k)x(k+i|k)$ is gotten by $F=\mathrm{YG}-1$. Matrix $F(k)$ minimizes the upper boundary of $V(x(k|k))$ in the objective function $J(k)$ at time instant $k$, and $G$ and $Y$ are achieved from the solution of the subsequent optimization problem with variables $M.G.Y.X.W.\xi .\gamma .$ and $Z=[G;Y]$ and $\epsilon$ indicates a positive coefficient. (9) $\begin{array}{r}\underset{M.G.Y.X.W.\xi .\gamma }{min}\gamma \end{array}$(9) Conditional on (10) $\begin{array}{rl}& \left[\begin{array}{cc}I& x(k|k{)}^T\\ x(k|k)& M\end{array}\right]\ge 0,\end{array}$(10) (11) $\begin{array}{rl}& \left[\begin{array}{ccccc}G+G^T-M& \ast & \ast & \ast & \ast \\ \sqrt{(1+\epsilon )(\mathrm{AG}+\mathrm{BY})}& M& \ast & \ast & \ast \\ \sqrt{(1+{\epsilon }^{-1})\mathrm{LZ}}& 0& \mathrm{\xi I}& \ast & \ast Q^{1/2}G& 0& 0& \mathrm{\gamma I}& \ast R^{1/2}Y& 0& 0& 0& \mathrm{\gamma I}\end{array}\right]\ge 0\end{array}$(11) (12) $\begin{array}{rl}& M-\mathrm{\xi I}>0.\end{array}$(12) (13) $\begin{array}{rl}& X-M\ge 0\mathrm{with}{X}_{\mathrm{rr}}\le x_{r.\max }^2.r=\mathrm{1.2.}\dots .n.\end{array}$(13) (14) $\begin{array}{rl}& \left[\begin{array}{cc}\begin{array}{c}W\\ Y^T\end{array}& \begin{array}{c}Y\\ G+G^T-M\end{array}\end{array}\right]\ge 0\mathrm{with}\\ & W_{\mathrm{rr}}\le u_{r.\max }^2.r=\mathrm{1.2.}\dots .m.\end{array}$(14)

Proof:

the following lemma is used to establish the proof of Theorem 3.1.

Lemma 3.1:

Suppose $\overline{M}$ and $\overline{N}$ are real constant matrices and $P$ is a positive matrix of compatible dimensions. Then ${\overline{M}}^{T}P\overline{N}+{\overline{N}}^{T}P\overline{M}\le \epsilon {\overline{M}}^{T}P\overline{M}+{\epsilon }^{-1}{\overline{N}}^{T}P\overline{N}$ holds for any $\epsilon >0$.

Proof:

See . [21]

Applying Schur complements to $V(x(k|k))\le \gamma$ in (8) and replacing $P=\gamma M^{-1}$, relation (10) is derived. To obtain (11), relation (7) is rewritten as (15) $\begin{array}{rl}& x(k+i+1|k{)}^T\mathrm{Px}(k+i+1|k)\\ & -x(k+i|k{)}^T\mathrm{Px}(k+i|k)\\ & +x(k+i|k{)}^T\mathrm{Qx}(k+i|k)\\ & -u(k+i|k{)}^T\mathrm{Ru}(k+i|k)<0.\end{array}$(15)

Considering Lemma 3.1, the following inequality can be found for the first term in (15). (16) $\begin{array}{rl}& x(k+i+1|k{)}^T\mathrm{Px}(k+i+1|k)\\ & =(\mathrm{Ax}(k+i|k)+\mathrm{Bu}(k+i|k)\\ & +f(x(k+i|k+u(k+i)){)}^{T}P(\mathrm{Ax}(k+i|k)\\ & +\mathrm{Bu}(k+i|k)+f(x(k+i|k),u(k+i|k)))\\ & \le (1+\epsilon )(\mathrm{Ax}(k+i|k)+\mathrm{Bu}(k+i|k){)}^T\\ & \times P(\mathrm{Ax}(k+i|k)+\mathrm{Bu}(k+i|k))\\ & +(1+{\epsilon }^{-1})f(x(k+i|k),u(k+i|k){)}^T\\ & \times \mathrm{Pf}(x(k+i),u(k+i|k))\end{array}$(16)

According to $P<\mathrm{\mu I}$ ($\mu$ indicates the largest Eigen value of $P$) and by using of (16), relation (15) can be stated as. (17) $\begin{array}{rl}& (1+\epsilon )(\mathrm{Ax}(k+i|k)+\mathrm{Bu}(k+i|k){)}^T\\ & \times P(\mathrm{Ax}(k+i|k)+\mathrm{Bu}(k+i|k))\\ & +(1+{\epsilon }^{-1})\mu [x(k+i|k);u(k+i|k){]}^T{L}^T\\ & \times L[x(k+i|k);u(k+i|k)]\\ & -x(k+i|k{)}^T\mathrm{Px}(k+i|k)+x(k+i|k{)}^T\mathrm{Qx}(k+i|k)\\ & +u(k+i|k{)}^T\mathrm{Ru}(k+i|k)<0.\end{array}$(17)

Replacing $u(k+i|k)=\mathrm{Fx}(k+i|k)$, relation (17) holds for every $x(k+i|k)$ if and only if (18) $\begin{array}{rl}& (1+\epsilon )(A+\mathrm{BF}{)}^T+(A+\mathrm{BF})\\ & +(1+{\epsilon }^{-1})\mu [I;F{]}^T{L}^{T}L[I;F]\\ & -P+Q+F^T\mathrm{RF}<0\end{array}$(18) is fulfilled. Let $P=\gamma M^{-1}$ and $\mu =\gamma {\xi }^{-1}$ and by pre and post multiplying (18) by means of $M$ and then applying Schur complements, it can be obtained (19) $\begin{array}{r}\left[\begin{array}{ccccc}M& \ast & \ast & \ast & \ast \\ \sqrt{(1+\epsilon )(A+\mathrm{BF})M}& M& \ast & \ast & \ast \sqrt{(1+{\epsilon }^{-1})L[I;F]M}& 0& \mathrm{\xi I}& \ast & \ast Q^{1/2}M& 0& 0& \mathrm{\gamma I}& \ast R^{1/2}\mathrm{FM}& 0& 0& 0& \mathrm{\gamma I}\end{array}\right]\ge 0.\end{array}$(19)

According to , (19) is guaranteed every time (11) is fulfilled. Considering $F=Y{G}^{-1}$, relation (11) can be rewritten as follows. (20) $\begin{array}{r}\left[\begin{array}{ccccccccccccccccccccc}G+G^T-M& \ast & \ast & \ast & \ast \sqrt{(1+\epsilon )(A+\mathrm{BF})G}& M& \ast & \ast & \ast (1+{\epsilon }^{-1})L[I;F]G& 0& \mathrm{\xi I}& \ast & \ast Q^{1/2}G& 0& 0& \mathrm{\gamma I}& \ast R^{1/2}\mathrm{FG}& 0& 0& 0& \mathrm{\gamma I}\end{array}\right]\ge 0.\end{array}$(20)

The inequality denotes that component $G^T+G-M$ in (20) should be a positive semi-definite matrix. On the other hand, since $M>0$ then $(G^T-M)M^{-1}(G-M)\ge 0$ is true, equivalent to (21) $\begin{array}{r}0\le G+G^T-M\le G^T{M}^{-1}G.\end{array}$(21)

By considering (21), relation (20) can be written as. (22) $\begin{array}{r}\left[\begin{array}{ccccccccccccccccccccc}{G}^T{M}^{-1}G& \ast & \ast & \ast & \ast (1+\epsilon )(A+\mathrm{BF})G& M& \ast & \ast & \ast (1+{\epsilon }^{-1})L[I;F]G& 0& \mathrm{\xi I}& \ast & \ast Q^{1/2}G& 0& 0& \mathrm{\gamma I}& \ast R^{1/2}\mathrm{FG}& 0& 0& 0& \mathrm{\gamma I}\end{array}\right]\ge 0.\end{array}$(22)

The relation (19) is gotten from (22) by pre-multiplying it by diag $(M{G}^{-T}.I.I.I)$ and post-multiplying by diag $(G^{-1}M.I.I.I)$. The condition $P<\mathrm{\mu I}$ and replacing $P=\gamma M^{-1}$ and $\mu =\gamma {\xi }^{-1}$ results in (12).

In Theorem 3.1, relation (10) states that $x(k|k)M^{-1}x(k|k)\le 1$ and (11) indicates that $V(x)$ is a decreasing function or equivalently $0$ for $i\ge 0$. Thus, $x(k+i|k).i\ge 0$ fits in an invariant set $S=\{x\in R^n|x^T{M}^{-1}x\le 1\}$. If there exists a symmetric matrix $X$ in a way that $X_{\mathrm{rr}}\le x_{r.\max }^2.r=\mathrm{1.2..}...n$, then relation (13) assurances the state constraints fulfilment in Theorem 3.1.

For driving (14), consider (10) and the fact that $x(k+i).i\ge 0$ fit in an invariant set $S=\{x\in R^n|x^T{M}^{-1}x\le 1\}$. Then input limitation can be described as (23) $\begin{array}{rl}& \underset{i\ge 0}{max}{u}_r(k+i|k{)}^2=\underset{i\ge 0}{max}|{(\mathrm{Fx}(k+i|k))}_r{|}^2\\ & =\underset{i\ge 0}{max}|Y{G}^{-1}{M}^{1/2}(M^{1/2}x(k+i|k)){)}_r{|}^2\\ & \le {(Y{G}^{-1}{M}^{1/2}{)}_r}_2^2.\end{array}$(23)

It can be achieved by using of Cauchy–Schwarz inequality (24) $\begin{array}{r}{(Y{G}^{-1}{M}^{1/2}{)}_r}_2^2=(Y{G}^{-1}{M}^{-1}{Y}^T{)}_{\mathrm{rr}}.\end{array}$(24)

So if there exists a symmetric matrix $W$ such that $W_{\mathrm{rr}}\le u_{r.\max }^2.r=\mathrm{1.2..}...m$ and (25) $\begin{array}{r}\left[\begin{array}{cc}W& W\\ Y^T& G^T{M}^{-1}G\end{array}\right]\ge 0,\end{array}$(25)

Then the input limitations in Theorem 3.1 are guaranteed. Bearing in mind (19), (25) is assured if (14) holds.

In the same way, the proof process presented for one of the agents can be generalized to all agents and stability with the same design can be achieved for all agents.

Thus for the nonlinear system in (5), RMPC based on Theory 1 can asymptotically stabilize the closed-loop system with the estimated region $S=\{x\in R^n|x^T{R}^{-1}x\le 1\}.$ where $R$ is achieved from $\begin{array}{c}\text{max}\\ M.G.Y.Z.X.W.\xi .\gamma \end{array}$ $\mathrm{logdet}(M)$ conditional on (11)–(14).

4. Simulation results and discussion

In this section, simulation in MATLAB environment is used to show the ability of the proposed predictive approach to control surge instability, cover uncertainty effects as well as high efficiency performance. The scenario is as follows: the entire simulation occurs in one second. Up to $t=0.3s$, $u_{\mathrm{th}}$ is equal to $0.1$ and after this time $u_{\mathrm{th}}$ is considered equal to $0.05$. Also in this simulation up to $t=0.5s$, $c_{\mathrm{th}}$ is considered as a nominal value, while after this time a tolerance of $0.25$ is embedded on $c_{\mathrm{th}}$. Also, for better evaluation, the results obtained using the proposed robust scheme have been compared with the NMPC and Tube-MPC techniques [16,22]. These techniques are nominal model dependent and do not cover any structural uncertainties. Of course, the NMPC technique cannot guarantee stability in the event of uncertainty, but it has been used to better compare and demonstrate the controller capability.

The simulation results are described below and the comparative results are plotted. The diagrams of compressor mass flow, plenum pressure increment, throttle pressure increment and plenum mass flow are shown in Figures 3–6. Figure 7 also shows the control signal obtained from the two NMPC and proposed robust LMI MPC techniques.

Figure 3. Compressor Mass Flow.

Figure 4. Plenum pressure.

Figure 5. Throttle pressure.

Figure 6. Plenum mass flow.

Figure 7. Control Signal.

As shown in Figures 3–7, by changing the $u_{\mathrm{th}}$ at $t=0.3s$, the operating point of the compressor system starts to change and moves to the new situation, and at the same time, NMPC and LMI-MPC controllers apply control signals to the compressor system by performing optimization. This process leads to limited changes in the compressor system states under the NMPC technique and, of course, milder changes occur in the compression system states under the proposed LMI MPC scheme, while at this stage both controllers prevent surge instability. The next inconvenience considered in this simulation, i.e. 25% uncertainty on $c_{\mathrm{th}}$, leads to further change in the behaviour of the compressor system states. Examination of the simulation results in Figures 3–7 shows that from $t=0.3s$ the all compressor system states fluctuates sharply under NMPC approach, while under the proposed robust LMI MPC technique, the compressor system states reach a new equilibrium point with optimal operation without surge instability despite drastic changes, and from $t=0.3s$ onwards, no oscillations are observed on the system states. The changes in the control signal shown in Figure 7 indicate the possibility of control signal to ensure compressor stability and surge-free operation.

The results obtained by applying the proposed LMI MPC and Tube-MPC approaches are now compared and shown in Figures 8–12. Although the Tube-MPC compressor does not experience negative flow, but a sharp drop occurs and the compressor operating point rapidly tilts toward the surge area, as shown in Figure 8. However, the current under the proposed LMI MPC method is far from the surge area and in stable operation mode.

Figure 8. Compressor mass flow.

Figure 9. Plenum pressure.

Figure 10. Throttle pressure.

Figure 11. Plenum mass flow.

Figure 12. Control Signal.

Carefully in Figure 9, it can be seen that small amplitude fluctuations in plenum pressure occur using the LMI MPC technique, and of course, more compression and faster achievement of a stable state is achieved in the proposed scheme, while less compression is obtained in the Tube-MPC technique, the amplitude of fluctuations is high, and of course, the plenum pressure reaches its final value in a longer period of time than the proposed method. It is also clear from Figure 10 that the results obtained for plenum pressure are also valid for throttle pressure. More compression, less oscillation amplitude and less settling time are the advantages of the proposed LMI-MPC technique. The plenum flow shown in Figure 11 also shows that less fluctuations and less settling time are advantages of the proposed LMI MPC technique and although the plenum flow is negative for a moment, it recovers rapidly and reaches its stable working point. The results obtained in Figure 12 show that less fluctuations occur in the control signal obtained from the proposed LMI MPC technique and reach its final value faster. Also, the control signal is in the range between 0 and 1, i.e. the practical range of the CCV actuator is well observed and the obtained signal is executable and operational.

The time of performing control signal calculations in three methods: NMPC, Tube-MPC and LMI MPC is given in Table 2. The calculations are performed on a computer with Intel Core i5 processor, Intel UHD Graphics and 8 GB memory and it is obvious from Table 2 that the time required to perform calculations in proposed approach is less than the other two techniques, and of course, it is very important in the case of online implementation and control.

Table 2. Calculation Time for NMPC, Tube-MPC and LMI-MPC Approaches.

Surge Controller	Robust LMI MPC	Tube-MPC	Nonlinear MPC
Calculation Time	0.249	0.317	0.359

5. Conclusion

In this paper, an optimal robust scheme is presented for controlling the surge instability of a constant speed centrifugal compressor system. The effects of uncertainty, input signal limitations, and system states constraints considered in the nonlinear model of compressor system including pipe effects, and a novel robust model predictive control scheme presented by using of linear matrix inequality scheme. Due to the simplicity of the proposed scheme structure, the time required to perform control signal calculations and optimization is less than other techniques and can be implemented online. The simulation results also indicated the capability of the LMI MPC technique to control surge, fast flow stabilization with insignificant fluctuations and obviously extra compression.

Disclosure statement

No potential conflict of interest was reported by the author(s).