Robust adaptive control of surge instability in constant speed centrifugal compressors using tube-MPC

Abstract

Surge is one of the two dynamic instabilities that occur in centrifugal compressors and can pose severe safety issues to the compression system. This paper tries to investigate the active surge control methods in centrifugal compressors. To this end, Moore–Greitzer Compressor Model is used, which considers a close couple valve as the independent variable. Most of the existing active controllers require information about system characteristics and upper bound of disturbances. In designing a controller, it is assumed that some bounded external disturbances are applied to the compressor, but its upper bound is unknown. As such, the throttle valve and the compressor’s characteristic are uncertain. In the proposed new active controller, a fuzzy system is used to approximate the compressor characteristic. Also, the adaptive control is used to estimate the openings of the throttle valve and the parameters of fuzzy system. Then, a robust controller is employed to find the appropriate gain in Tube-MPC to avoid the system disturbances. Finally, by simulating the compressor’s behavior after applying the designed active controller, one can prove its efficiency and robustness to uncertainty in model and disturbance rejection.

Keywords:

• compressor
• surge
• active control
• close couple valve
• robust control
• adaptive control
• fuzzy system
• tube-MPC

Public Interest Statement

In this paper, an active controller is designed using Tube-MPC, a model predictive controller that can control system even when the model and the process are not matched. As such, in designing the controller it is assumed that the compressor characteristic is not known. Then, it is approximated by a fuzzy system. Next, the requirement of information regarding throttle valve is removed. This is why an adaptive controller has been utilized. Finally, a robust controller is used to find the ancillary controller gain in Tube-MPC. This leads to a controller robust to disturbances.

1. Introduction

Compressors play a key role in petroleum and energy industry as the heart of gas transformation systems, injection systems for petroleum, and gas fields. They can help to increase the efficiency of pickup, liquidized natural gas plants, and converting gas to liquid. The main limiting factor of a centrifugal gas compressor is aerodynamic instability that appears in the form of surge or rotational stall. However, avoiding surge or rotational stall can restrict the optimal performance and consequently reduces the plant yield. The surge control system in gas compressors is a basic and important issue in dynamic stability and determines the operating area of the machine.

There are two control methods to avoid surge: passive control and active control. The passive surge control is a general surge control method (Willems, 1999). This strategy, conducted by defining a surge control line that is chosen at a distance from the surge line, tries to keep the compressor’s equilibrium point to the right side of the surge line by recycling or blowing off the gas compressor. As a result of lag in control valve response times, the surge control line cannot be close to the surge line and it is usually selected at a distance of 10–20%. The weakness of this method is in proper selection of the distance of the working point from the surge line. If the working point is far away from the surge line, the efficiency is reduced and the compressor cannot produce a high pressure.

Active surge control is the opposite point of the above-mentioned strategy, since it extends the surge activity area. This method was introduced for the first time by Epstein, Williams, and Greitzer (1989) where the compressor was allowed to work in the instability area using some control feedback. In this method, when the working point becomes closer to the surge line, the compressor efficiency increases and the compressor can produce a higher pressure.

Many research papers have focused on active surge control, some of which are mentioned in the following. First, different operators have been used for active surge control, including piston damping mechanism, close-couple valve (CCV), throttle control valve (TCV) (Wang, Zhao, & Luo, 2006), hydraulic (Arnulfi, Giannattasio, Micheli, & Pinamonti, 2001) and motor rotation control (Bøhagen & Gravdahl, 2008).

Willems, Heemels, de Jager, and Stoorvogel (2002) presented a new controller, positive feedback, using pole placement method to actively control the surge. This controller has been conducted on an experimental sample and its efficiency has been proved, but there are four problems when using positive feedback control: complex open-loop poles, the domain of attraction of a stabilized equilibrium point, the system’s robustness, the disturbance and noise rejection properties (Zillinger Molenaar, 2007). Zillinger Molenaar (2007) improved the positive-feedback controller (Willems et al., 2002) by combining it with an MPC in order to solve the problems (2), (3) and (4). The obtained results from this combined controller were better than the original positive feedback controller. However, problem (1) is not solved yet, and does not consider all the nonlinear system dynamics because the system model is linear. Since surge is in fact a limit cycle, it is necessary to investigate the nonlinear model of the system. Employing nonlinear controllers is a better approach in active surge control.

According to the bifurcation theory, a classic nonlinear dynamics (a simple active controller) can be used for surge control (Abed, Houpt, & Hosny, 1993; Wang, Li, & Zhang, 2007). In (Cortinovis, Pareschi, Mercangoez, & Besselmann, 2012), MPC (an active controller) was used for surge control by controlling motor rotation and the opening of the control valve. This can stabilize the compressor working point in every point in the performance curve of the compressor. However, it was assumed in these references that the system model is known. Shehata, Abdullah, and Areed (2009) investigated the application of active surge control in constant speed centrifugal compressors according to Moore–Greitzer (MG) model. The control operators in this method were CCV and TCV. The applied control methods in this work were fuzzy, anti-windup PID and Gain scheduling, sliding mode and backstepping; however, the effect of disturbances and uncertainly were not considered in this study. MPC approaches for nonlinear centrifugal compression systems, controlled by torque and recycle valve actuation, are discussed in Torrisi et al. (in press). These systems are usually used in critical processes, where the compressed fluid can impose an irreversible damage to the compressor. Bentaleb, Cacitti, De Franciscis, and Garulli (2015) present a new multivariable predictive control model for centrifugal compressors in gas transportation plants. The results of the numerical simulations indicate that, under different types of variations in the upstream and downstream conditions of the pipeline, the proposed scheme can meet the requirements of discharge pressure regulation as well as surge prevention. Cortinovis, Ferreau, Lewandowski, and Mercangöz (2015) executes an experimental setup to test the compressor control. In addition to that, the proposed model, the parameter identification, and the model validation are also presented in Cortinovis et al. (2012). This model is based on a linear MPC controller, which regulates the recycle valve opening and the driver torque. To avoid disturbances, the torque action is increased in the first transient phase of the disturbance, which helps the anti-surge controller.

Most of the applied methods in active control require a complete knowledge about the compressor model. The uncertain dynamics and disturbances may prevent the active control from presenting a good performance. This is why the topic of active control, along with an uncertainty existing in the system, is very controversial (Chetate, Zamoum, Fegriche, & Boumdin, 2013; Javadi Moghaddam, Farahani, & Amanifard, 2011; Nayfeh & Abed, 2001; Sheng, Huang, Zhang, & Huang, 2014a). In (Sheng, Huang, Zhang, & Huang, 2014b) an adaptive robust controller, using backstepping method, has been presented for surge control, in which the compressor characteristic is identified by a fuzzy system. This controller is robust to unknown, bounded disturbances. However, the weakness of the controller (Sheng et al., 2014b) is its negative control signal, which is not acceptable for the CCV operator output.

As it was mentioned above, the most significant weakness of active controllers is their requirement of the system knowledge. In industry, there is a need to have an active controller that can avoid surge despite the existence of uncertainty in system model and compressor characteristic. In this paper, an active controller is designed using Tube-MPC, a model predictive controller that can control system even when the model and the process are not matched. As such, in designing the controller it is assumed that the compressor characteristic is not known. Then, it is approximated by a fuzzy system. Next, the requirement of information regarding throttle valve is removed. This is why an adaptive controller has been utilized. Finally, a robust controller is used to find the ancillary controller gain in Tube-MPC. This leads to a controller robust to disturbances.

The rest of this paper is organized as the following:

Section 2 presents the new robust adaptive Tube-MPC. Section 3 designs a robust adaptive fuzzy Tube-MPC controller to stabilize the compressor system. Section 4, conducts a simulation to prove the efficiency of the controller. Finally, Section 5 concludes the paper.

2. Constant speed centrifugal compressor model

Given the instable nature of the surge phenomenon in compressor and in order to design a suitable controller, a mathematical model is required to describe the surge and dynamics of the compressor flow. Greitzer in 1976 presented the first nonlinear model for axial compressors (Greitzer, 1976a, 1976b). However, Hansen et al. in 1981 (Hansen, Jo̸rgensen, & Larsen, 1981) showed that Greitzer’s model can be applied to centrifugal compressors as well. Greitzer’s model shows surge cycle completely, but does not incorporate rotational stall that leads to pressure drop. In 1986, Moore and Greitzer (Moore & Greitzer, 1986) proposed a new model (MG Model) that was able to describe surge and rotational stall simultaneously. Moore and Greitzer model equations for the centrifugal compressor are as the followings$$\frac{\text{d}\mathit{\psi }}{\text{d}\mathit{\xi }}=\frac{1}{4B^2{l}_{\text{c}}}\left(\mathit{\varphi }-{\mathit{\varphi }}_{\text{T}}\left(\mathit{\psi }\right)\right)$$

$$\frac{\text{d}\mathit{\varphi }}{\text{d}\mathit{\xi }}=\frac{1}{l_{\text{c}}}\left({\mathit{\psi }}_{\text{c}}\left(\mathit{\varphi }\right)-\mathit{\psi }-\frac{3H}{4}\left(\frac{\mathit{\varphi }}{W}-1\right)J\right)$$

(1) $$\frac{\text{d}J}{\text{d}\mathit{\xi }}=J\left(1-{\left(\frac{\mathit{\varphi }}{W}-1\right)}^2-\frac{J}{4}\right)\mathit{\delta }$$(1)

where $\mathit{\psi }$ is the increase coefficient of compressor pressure, $\mathit{\varphi }$ is the coefficient of compressor mass flow, J denotes the square of rotational stall range, and ξ is the non-dimensionalization time. Also, ${\mathit{\varphi }}_{\text{T}}\left(\mathit{\psi }\right)$ is the characteristic of throttle valve and ${\mathit{\psi }}_{\text{c}}\left(\mathit{\phi }\right)$ is the characteristic of the compressor. B is the Greitzer’s parameter and l_c shows the length of canals (ducts) and compressor and δ is a constant coefficient that is obtained from the following relation(2) $$B=\frac{U}{2a_{\text{s}}}\sqrt{\frac{V_{\text{p}}}{A_{\text{c}}{l}_{\text{c}}}}$$(2)

where U is the constant of compressor tangential speed, a_s is the sonic speed, V_p is the plenum volume, and A_c is the cross-section of the compressor.

A compressor with constant speed, flow rate m, pressure difference ΔP, and real time t is normalized according to the following formula(3) $$\mathit{\varphi }=\frac{m}{\mathit{\rho }{A}_{\text{c}}U}\mathit{\psi }=\frac{\mathrm{\Delta }P}{\mathit{\rho }{U}^2}\mathit{\xi }=\frac{U}{R}t$$(3)

This normalization or non-dimensionalization converts a class of curves for the compressor map into a single characteristic in every speed. This characteristic is a non-linear relationship between the compressor pressure $\mathit{\psi }$ and compressor flow $\mathit{\varphi }$. Various statements have been used for this characteristic, but Moore and Greitzer’s (1986) cube characteristic in 1986 is one the most well-known, which is defined as(4) $${\mathit{\psi }}_{\text{c}}\left(\mathit{\varphi }\right)={\mathit{\psi }}_{\text{c0}}+H\left(1+\frac{3}{2}\left(\frac{\mathit{\varphi }}{W}-1\right)-\frac{1}{2}{\left(\frac{\mathit{\varphi }}{W}-1\right)}^3\right)$$(4)

where ${\mathit{\psi }}_{\text{c0}}$ is the value of characteristic curve in zero dB, H is half of the height of the characteristic curve, and W is the half of the width of the characteristic curve (Sheng et al., 2014b).

The equation for throttle valve characteristic is also derived from (Gravdahl & Egeland, 1999) and is as follows(5) $${\mathit{\varphi }}_{\text{T}}\left(\mathit{\psi }\right)={\mathit{\gamma }}_{\text{T}}\sqrt{\mathit{\psi }}$$(5)

where γ_T is also the valve’s yield.

In (1), J = 0 is considered to obtain the surge equations. So, the equations are rewritten as$$\dot{\mathit{\psi }}=\frac{1}{4B^2{l}_{\text{c}}}\left(\mathit{\varphi }-{\mathit{\varphi }}_{\text{T}}\left(\mathit{\psi }\right)\right)$$

(6) $$\dot{\mathit{\varphi }}=\frac{1}{l_{\text{c}}}\left({\mathit{\psi }}_{\text{c}}\left(\mathit{\varphi }\right)-\mathit{\psi }\right)$$(6)

Like other physical systems, this system also has some disturbance inherent in it. This disturbance leads to instability and causes the active controller function to work improperly. The disturbances of flow and pressure are applied to the compressor equations as follows$$\dot{\mathit{\psi }}=\frac{1}{4B^2{l}_{\text{c}}}\left(\mathit{\varphi }-{\mathit{\varphi }}_{\text{T}}\left(\mathit{\psi }\right)+{\text{d}}_{\mathit{\varphi }}\left(\mathit{\xi }\right)\right)$$

(7) $$\dot{\mathit{\varphi }}=\frac{1}{l_{\text{c}}}\left({\mathit{\psi }}_{\text{c}}\left(\mathit{\varphi }\right)-\mathit{\psi }+{\text{d}}_{\mathit{\psi }}\left(\mathit{\xi }\right)\right)$$(7)

Figure 1 shows the diagram of compression system with CCV.

Figure 1. The compressor system with CCV.

Source: Sheng et al. (2014b).

Figure 2. Control system block diagram.

The system model equations, considering a CCV, are$$\dot{\mathit{\psi }}=\frac{1}{4B^2{l}_{\text{c}}}\left(\mathit{\varphi }-{\mathit{\varphi }}_{\text{T}}\left(\mathit{\psi }\right)+{\text{d}}_{\mathit{\varphi }}\left(\mathit{\xi }\right)\right)$$

(8) $$\dot{\mathit{\varphi }}=\frac{1}{l_{\text{c}}}\left({\mathit{\psi }}_{\text{c}}\left(\mathit{\varphi }\right)-\mathit{\psi }-{\mathit{\psi }}_{\text{V}}\left(\mathit{\varphi }\right)+{\text{d}}_{\mathit{\psi }}\left(\mathit{\xi }\right)\right)$$(8)

Considering ${\mathit{\psi }}_{\text{V}}\left(\mathit{\varphi }\right)$ as the input for system control and $x_1=\mathit{\psi },x_2=\mathit{\varphi }$, the equations of compressor state space are$${\dot{x}}_1=\frac{1}{4B^2{l}_{\text{c}}}\left(x_2-{\mathit{\varphi }}_{\text{T}}\left(x_1\right)+{\text{d}}_{\mathit{\varphi }}\left(\mathit{\xi }\right)\right)$$

(9) $${\dot{x}}_2=\frac{1}{l_{\text{c}}}\left({\mathit{\psi }}_{\text{c}}\left(x_2\right)-x_1-u+{\text{d}}_{\mathit{\psi }}\left(\mathit{\xi }\right)\right)$$(9)

Compressor parameters are listed in Table 1.

Table 1. Description of symbols and parameters

Symbols	Parameters
$\mathit{\varphi }$	Mass flow coefficient
$\mathit{\psi }$	Pressure rise coefficient
$l_{\text{I}}$	Inlet duct length
$l_{\text{c}}$	Dimensionless compressor length
B	Greitzer’s parameter
as	Speed of sound
lE	Exit duct length
Ac	Compressor duct area
U	Mean rotor speed
R	Rotor mean radius
VP	Plenum volume
ξ	Dimensionless time
t	Time
u	Valve opening
$L_{\text{c}}$	Compressor and ducts length
a	Reciprocal time lag of compressor passage
${\mathit{\psi }}_{\text{C0}}$	Compressor characteristic shut-off value
W	Compressor characteristic semi width
H	Compressor characteristic semi height
${\mathit{\varphi }}_{\text{T}}$	Throttle valve
${\mathit{\psi }}_{\text{V}}$	Close-couple valve
$d_{\mathit{\varphi }}$	Disturbance
$d_{\mathit{\psi }}$	Disturbance

Source: Shehata et al. (2009).

3. Robust adaptive fuzzy model predictive control on surge

3.1. Description of fuzzy system

A Fuzzy system is defined as a system that presents an outline showing the input vector to output vector scheme: x → y, where $x=\left[x_1,\dots ,x_n\right]\in X_1\times \dots \times X_n\subseteq R^n,y\in R$ given the fact that fuzzy system has an integrated fuzzification, Gauss membership function, product inference, and principal mean defuzzification, the ith law of fuzzy logic system looks like this:

Rule $i:\text{if}{x}_i\text{is}{F}_{i1},x_n\text{is}{F}_{in}\text{then}y=w_i$, where $i=1,\dots ,m$, m shows how many fuzzy logic rules are there, w_i displays the ith fuzzy law, $F_{ij}\left(j=1,\dots ,n\right)$ shows a fuzzy collection in the world of discourse X_i, which employs Gauss function as membership function(10) $$\mathit{\mu }{F}_{ij}\left(x_j\right)=e^{-{\left(\frac{x_j-a_{ij}}{b_{ij}}\right)}^2}$$(10)

a_ij and b_ij are design components.

To achieve the fuzzy system, the fuzzy laws should be combined one by one, i.e.(11) $$y\left(x\right)=\frac{\sum _{i=1}^m{w}_i\left(\prod _{j=1}^n\mathit{\mu }{F}_{ij}\left(x_j\right)\right)}{\sum _{i=1}^m\left(\prod _{j=1}^n\mathit{\mu }{F}_{ij}\left(x_j\right)\right)}$$(11) (12) $$W^{\text{T}}=\left[w_1,\dots ,w_m\right]$$(12) (13) $$P\left(x\right)={\left[p_1\left(x\right),p_2\left(x\right),\dots ,p_m\left(x\right)\right]}^{\text{T}}$$(13) (14) $$p_i\left(x\right)=\frac{\prod _{j=1}^n\mathit{\mu }{F}_{ij}\left(x_j\right)}{\sum _{i=1}^m\left(\prod _{j=1}^n\mathit{\mu }{F}_{ij}\left(x_j\right)\right)}$$(14)

Then, a constant value should be specified for the membership function, (i.e. a_ij and b_ij constant), and the fuzzy law w_i must be described as a variable component. Then we have(15) $$y\left(x\right)=W^{\text{T}}P\left(x\right)$$(15)

In which, $P\left(x\right)$ is fuzzy basis function vector, W component vector.

There exists a fuzzy system $y^{\ast }\left(x\right)$ for each real continuous function $y\left(x\right)$ in the set X ⊆ Rⁿ and also for each real number ɛ > 0, that meets sup $\left|y^{\ast }\left(x\right)-y\left(x\right)\right|<\mathit{\epsilon }$; so, a fuzzy system can be used to estimate a continuous function $y\left(x\right)$: $f\left(x\right)=W^{\ast \text{T}}P\left(x\right)+\mathrm{\Delta }f\left(x\right)$(16) $$f\left(x\right)=W^{\ast \text{T}}P\left(x\right)+\mathrm{\Delta }f\left(x\right)$$(16)

where $\mathrm{\Delta }f\left(x\right)$ satisfies $\mathrm{\Delta }f\left(x\right)<\mathit{\epsilon }$ (Wang, 1994).

3.2. Controller design

In designing a surge controller for the compressor system (9), it is assumed that the throttle valve opening, as well as the compressor characteristic, are not known. So, the fuzzy system is used to approximate it. To do this, a 9-membership is used, with the following equations(17) $$p_i\left(x\right)=e^{-{\left(x-0.1i\right)}^2},i=1,2,\dots ,9$$(17)

Having used (17), the system characteristic curve is obtained as(18) $${\mathit{\psi }}_{\text{c}}\left(x_2\right)=W^{\text{*T}}P\left(x_2\right)+\mathrm{\Delta }{\mathit{\psi }}_{\text{c}}\left(x_2\right)$$(18)

where W^*T is the weight of fuzzy function members. Also, $\mathrm{\Delta }{\mathit{\psi }}_{\text{c}}\left(x_2\right)$ is the approximate system error. Then the compressor equations are obtained$${\dot{x}}_1=\frac{1}{4B^2{l}_{\text{c}}}\left(x_2-{\mathit{\gamma }}_T\sqrt{x_1}+{\text{d}}_{\mathit{\varphi }}\left(\mathit{\xi }\right)\right)$$

(19) $${\dot{x}}_2=\frac{1}{l_{\text{c}}}\left(-x_1-u+W^{\text{*T}}P\left(x_2\right)+\mathrm{\Delta }{\mathit{\psi }}_{\text{c}}\left(x_2\right)+{\text{d}}_{\mathit{\psi }}\left(\mathit{\xi }\right)\right)$$(19)

Next, the compressor equations are rewritten as follows(20) $$\dot{x}\left(\mathit{\xi }\right)=Ax\left(\mathit{\xi }\right)+Bu\left(\mathit{\xi }\right)+f\left(x\right)\mathit{\theta }\left(\mathit{\xi }\right)+\text{d}\left(\mathit{\xi }\right)$$(20)

where $x\left(\mathit{\xi }\right)=\left[\begin{array}{c}{x}_1\hfill \\ x_2\hfill \end{array}\right]$ and(21) $$\begin{array}{cc}A=\left[\begin{array}{cc}0\hfill & \frac{1}{4B^2{l}_{\text{c}}}\hfill \\ -\frac{1}{l_{\text{c}}}\hfill & 0\hfill \end{array}\right],\hfill & B=\left[\begin{array}{c}0\hfill \\ -\frac{1}{l_{\text{c}}}\hfill \end{array}\right]\hfill \end{array}$$(21) (22) $$\begin{array}{cc}f\left(x\right)=\left[\begin{array}{cc}-\frac{1}{4B^2{l}_{\text{c}}}\sqrt{x_1}\hfill & 0\hfill \\ 0\hfill & \frac{1}{l_{\text{c}}}P\left(x_2\right)\hfill \end{array}\right],\hfill & \mathit{\theta }\left(\mathit{\xi }\right)=\left[\begin{array}{c}{\mathit{\gamma }}_{\text{T}}\hfill \\ W^{\ast }\hfill \end{array}\right]\hfill \end{array}$$(22) (23) $$d\left(\mathit{\xi }\right)=\left[\begin{array}{c}\frac{1}{4B^2{l}_{\text{c}}}{d}_{\mathit{\varphi }}\left(\mathit{\xi }\right)\\ \frac{1}{l_{\text{c}}}\left(\mathrm{\Delta }{\mathit{\psi }}_{\text{c}}\left(x_2\right)+{\text{d}}_{\mathit{\psi }}\left(\mathit{\xi }\right)\right)\end{array}\right]$$(23)

Now, the nominal model of the compressor system is given as(24) $$\dot{\overline{x}}\left(\mathit{\xi }\right)=A\overline{x}\left(\mathit{\xi }\right)+B\overline{u}\left(\mathit{\xi }\right)+f\left(\overline{x}\right)\widehat{\mathit{\theta }}$$(24)

where $\overline{x}\left(\mathit{\xi }\right)\in R^{n_x}$ shows the nominal model states, $\overline{u}(\mathit{\xi })\in R^{n_u}$ is the control input of nominal model and $\widehat{\mathit{\theta }}(\mathit{\xi })\in R^{n_{\mathit{\theta }}}$ is the estimator of θ uncertainty in the compressor system.

The cost function is defined as(25) $$J\left(\overline{x},\overline{u}\right)={\int }_0^{T_{\text{p}}}\left(\overline{x}{\left(\mathit{\tau }\right)}^{\text{T}}Q\overline{x}\left(\mathit{\tau }\right)+\overline{u}{\left(\mathit{\tau }\right)}^{\text{T}}R\overline{u}\left(\mathit{\tau }\right)\right)\text{d}\mathit{\tau }$$(25)

where T_p is the prediction horizon. In the objective function (25), the values of weight matrixes Q and R are selected as(26) $$Q=\left[\begin{array}{cc}100\hfill & 0\hfill \\ 0\hfill & 100\hfill \end{array}\right],R=1$$(26)

Next, the existing constraints in the compressor system should be incorporated in the optimization problem. The first existing constraint is on the control signal. Since the control signal has a CCV output, so we have(27) $$\overline{u}\left(\mathit{\xi }\right)>0$$(27)

The next constraint and limitation is that the flow has some maximum and minimum values. This constraint should also be considered.(28) $$-{\mathit{\varphi }}_{\text{m}}\le \mathit{\varphi }\left(\mathit{\xi }\right)\le {\mathit{\varphi }}_{\text{Choke}}$$(28)

The following problem is solved to find $\overline{u}(\mathit{\xi })$(29) $$\begin{array}{c}\text{minimize}J\left(\overline{x}\left(\mathit{\xi }\right),\overline{u}\left(\mathit{\xi }\right)\right)\hfill \\ \overline{u}\left(\mathit{\xi }\right)\hfill \\ s.t\hfill \\ \dot{\overline{x}}\left(\mathit{\xi }\right)=A\overline{x}\left(\mathit{\xi }\right)+B\overline{u}\left(\mathit{\xi }\right)+f\left(\overline{x}\right)\widehat{\mathit{\theta }}\hfill \\ \overline{u}\left(\mathit{\xi }\right)>0\hfill \\ -{\mathit{\varphi }}_{\text{m}}\le \mathit{\varphi }\left(\mathit{\xi }\right)\le {\mathit{\varphi }}_{\text{Choke}}\hfill \end{array}$$(29)

Consequently, the overall applied control input for the actual system (19) is(30) $$u\left(\mathit{\xi }\right)=\overline{u}\left(\mathit{\xi }\right)+Ke\left(\mathit{\xi }\right)$$(30)

where(31) $$e\left(\mathit{\xi }\right)=x\left(\mathit{\xi }\right)-\overline{x}\left(\mathit{\xi }\right)$$(31)

According to (20), (24) and (30), the dynamic error equations are given as(32) $$\dot{e}\left(\mathit{\xi }\right)=\left(A+BK\right)e\left(\mathit{\xi }\right)+f\left(x\left(\mathit{\xi }\right)\right)\mathit{\theta }-f\left(\overline{x}\left(\mathit{\xi }\right)\widehat{\mathit{\theta }}\right)+d\left(\mathit{\xi }\right)$$(32)

By defining $\tilde{\mathit{\theta }}=\mathit{\theta }-\widehat{\mathit{\theta }}$, we have(33) $$\dot{e}\left(\mathit{\xi }\right)=\left(A+BK\right)e\left(\mathit{\xi }\right)+B_{\text{w}}w\left(\mathit{\xi }\right)$$(33)

where(34) $$g\left(\mathit{\xi }\right)=\left(f\left(x\left(\mathit{\xi }\right)\right)-f\left(\overline{x}\left(\mathit{\xi }\right)\right)\right)\widehat{\mathit{\theta }}+\text{d}\left(\mathit{\xi }\right),w\left(\mathit{\xi }\right)=\left[\begin{array}{c}\tilde{\mathit{\theta }}\hfill \\ g\left(\mathit{\xi }\right)\hfill \end{array}\right],B_{\text{w}}=\left[f\left(x\left(\mathit{\xi }\right)\right)I_{n_x\times n_x}\right]$$(34)

Lemma 1

(Yu, Bohm, Chen, & Allgower, 2010): Let $S:R^{n_x}\to \left[0,\infty \right)$ be a continuously differentiable function and ${\mathit{\alpha }}_1\left(\Vert z\Vert \right)<S\left(z\right)<{\mathit{\alpha }}_2\left(\Vert z\Vert \right)$, where α₁, α₂ are class k_∞ functions. Suppose $u:R\to R^{n_u}$ is chosen, and there exists $\mathit{\lambda }>0$ and μ > 0 such that(35) $$\dot{S}\left(z\right)+\mathit{\lambda }S\left(z\right)-\mathit{\mu }{w}^{\text{T}}\left(t\right)w\left(t\right)\le 0$$(35)

With $z\in X,d\in D$. Then, the system trajectory starting from $z\left(t_0\right)\in \mathrm{\Omega }\subseteq X$, will remain in the set $\mathrm{\Omega }$, where(36) $$\mathrm{\Omega }=\left\{z\in R^{n_x}|S\left(z\right)\le \frac{\mathit{\mu }{w}_{max}^2}{\mathit{\lambda }}\right\}$$(36)

where w_max is the upper bound of $w\left(t\right)$.

Lemma 2

(Poursafar, Taghirad, & Haeri, 2010). Let M, N be real constant matrices and P be a positive matrix of compatible dimensions. Then(37) $$M^{T}PN+N^{T}PM\le \mathit{\epsilon }{M}^{T}PM+{\mathit{\epsilon }}^{-1}{N}^{T}PN$$(37)

Holds for any ɛ > 0.

Lemma 3

Suppose that there exist a positive definite matrix $X\in R^{n_x\times n_x}$, a non-square matrix $Y\in R^{n_u\times n_x}$, and scalars $\mathit{\alpha }>0,\mathit{\beta }>0,\mathit{\lambda }>0,\mathit{\eta }>0,\mathit{\epsilon }=\frac{\mathit{\alpha }\mathit{\lambda }}{\mathit{\eta }\mathit{\beta }}$ and $\mathit{\mu }=\frac{\mathit{\lambda }}{\mathit{\eta }\left(1+\mathit{\beta }\right)}$ such that(38) $${\left(AX+BY\right)}^T+AX+BY+\left(\mathit{\alpha }+\mathit{\lambda }\right)X\le 0$$(38) (39) $$X\le \frac{1}{\mathit{\epsilon }}I$$(39)

Then, the set $\mathrm{\Omega }=\left\{e\in R^{n_x}|S\left(e,\widehat{\mathit{\theta }}\right)\le \frac{\mathit{\mu }{w}_{max}^2}{\mathit{\lambda }}\right\}$ is a robust invariant set for the error system (33), where(40) $$u\left(\mathit{\xi }\right)=\overline{u}\left(\mathit{\xi }\right)+Ke\left(\mathit{\xi }\right)$$(40) (41) $$S\left(e\left(\mathit{\xi }\right),\tilde{\mathit{\theta }}\right)=e^T\left(\mathit{\xi }\right)Pe\left(\mathit{\xi }\right)+\frac{1}{\mathit{\eta }}{\tilde{\mathit{\theta }}}^T\tilde{\mathit{\theta }}$$(41)

And(42) $$P=X^{-1},K=Y{X}^{-1}$$(42)

Proof

According to Lemma 1, for system (33) we have(43) $$\dot{S}\left(e\left(\mathit{\xi }\right),\tilde{\mathit{\theta }}\right)+\mathit{\lambda }S\left(e\left(\mathit{\xi }\right),\tilde{\mathit{\theta }}\right)-\mathit{\mu }{w}^T\left(\mathit{\xi }\right)w\left(\mathit{\xi }\right)\le 0$$(43)

Then, according to (41) we have(44) $${\dot{e}}^T\left(\mathit{\xi }\right)Pe\left(\mathit{\xi }\right)+e^T\left(\mathit{\xi }\right)P\dot{e}\left(\mathit{\xi }\right)-\frac{1}{\mathit{\eta }}\left({\dot{\widehat{\mathit{\theta }}}}^T\tilde{\mathit{\theta }}+{\dot{\widehat{\mathit{\theta }}}}^T\right)+\mathit{\lambda }\left(e^T\left(\mathit{\xi }\right)Pe\left(\mathit{\xi }\right)+\frac{1}{\mathit{\eta }}{\tilde{\mathit{\theta }}}^T\tilde{\mathit{\theta }}\right)-\mathit{\mu }{w}^T\left(\mathit{\xi }\right)w\left(\mathit{\xi }\right)\le 0$$(44) (45) $$e^T\left(\mathit{\xi }\right)\left({\left(A+BK\right)}^{T}P+P\left(A+BK\right)+\mathit{\lambda }P\right)e\left(\mathit{\xi }\right)+\frac{\mathit{\lambda }}{\mathit{\eta }}{\tilde{\mathit{\theta }}}^T\tilde{\mathit{\theta }}+g^T\left(\mathit{\xi }\right)Pe\left(\mathit{\xi }\right)+e^T\left(\mathit{\xi }\right)Pg\left(\mathit{\xi }\right)-\mathit{\mu }\left({\tilde{\mathit{\theta }}}^T\tilde{\mathit{\theta }}+{\tilde{\mathit{\theta }}}^T\left(\mathit{\xi }\right)g\left(\mathit{\xi }\right)+g^T\left(\mathit{\xi }\right)\tilde{\mathit{\theta }}\left(\mathit{\xi }\right)+g^T\left(\mathit{\xi }\right)g\left(\mathit{\xi }\right)\right)+{\tilde{\mathit{\theta }}}^T{f}^T\left(x\left(\mathit{\xi }\right)\right)Pe\left(\mathit{\xi }\right)+e^T\left(t\right)Pf\left(x\left(\mathit{\xi }\right)\right)\tilde{\mathit{\theta }}-\frac{1}{\mathit{\eta }}\left({\dot{\widehat{\mathit{\theta }}}}^T\tilde{\mathit{\theta }}+{\tilde{\mathit{\theta }}}^T\dot{\widehat{\mathit{\theta }}}\right)\le 0$$(45)

By choosing(46) $$\dot{\widehat{\mathit{\theta }}}=\mathit{\eta }{f}^T\left(x\left(\mathit{\xi }\right)\right)Pe\left(\mathit{\xi }\right)$$(46)

We have(47) $$g^T\left(\mathit{\xi }\right)Pe\left(\mathit{\xi }\right)+e^T\left(\mathit{\xi }\right)Pg\left(\mathit{\xi }\right)\le \mathit{\alpha }{e}^T\left(\mathit{\xi }\right)Pe\left(\mathit{\xi }\right)+{\mathit{\alpha }}^{-1}{g}^T\left(\mathit{\xi }\right)Pg\left(\mathit{\xi }\right){\tilde{\mathit{\theta }}}^T\left(\mathit{\xi }\right)g\left(\mathit{\xi }\right)+g^T\left(\mathit{\xi }\right)\tilde{\mathit{\theta }}\left(\mathit{\xi }\right)\le \mathit{\beta }{\tilde{\mathit{\theta }}}^T\left(\mathit{\xi }\right)\tilde{\mathit{\theta }}\left(\mathit{\xi }\right)+{\mathit{\beta }}^{-1}{g}^T\left(\mathit{\xi }\right)g\left(\mathit{\xi }\right)$$(47)

According to Lemma 2, we have(48) $$g^T\left(\mathit{\xi }\right)Pe\left(\mathit{\xi }\right)+e^T\left(\mathit{\xi }\right)Pg\left(\mathit{\xi }\right)\le \mathit{\alpha }{e}^T\left(\mathit{\xi }\right)Pe\left(\mathit{\xi }\right)+{\mathit{\alpha }}^{-1}{g}^T\left(\mathit{\xi }\right)Pg\left(\mathit{\xi }\right){\tilde{\mathit{\theta }}}^T\left(\mathit{\xi }\right)g\left(\mathit{\xi }\right)+g^T\left(\mathit{\xi }\right)\tilde{\mathit{\theta }}\left(\mathit{\xi }\right)\le \mathit{\beta }{\tilde{\mathit{\theta }}}^T\left(\mathit{\xi }\right)\tilde{\mathit{\theta }}\left(\mathit{\xi }\right)+{\mathit{\beta }}^{-1}{g}^T\left(\mathit{\xi }\right)g\left(\mathit{\xi }\right)$$(48)

By substituting (48) in (47), it is obtained that(49) $$e^T\left(\mathit{\xi }\right)\left({\left(A+BK\right)}^{T}P+P\left(A+BK\right)+\left(\mathit{\alpha }+\mathit{\lambda }\right)P\right)e\left(\mathit{\xi }\right)+\left(\frac{\mathit{\lambda }}{\mathit{\eta }}-\mathit{\mu }\left(1+\mathit{\beta }\right)\right){\tilde{\mathit{\theta }}}^T\tilde{\mathit{\theta }}+{\mathit{\alpha }}^{-1}{g}^T\left(\mathit{\xi }\right)Pg\left(\mathit{\xi }\right)-\mathit{\mu }\left(1+{\mathit{\beta }}^{-1}\right)g^T\left(\mathit{\xi }\right)g\left(\mathit{\xi }\right)\le 0$$(49)

Consider(50) $$P\le {\mathit{\lambda }}_{\text{max}}I\le \mathit{\epsilon }I$$(50)

where ${\mathit{\lambda }}_{max}$ is the maximum eigenvalue of P and ɛI is the corresponding upper bound (Poursafar et al., 2010), then(51) $$e^T\left(\mathit{\xi }\right)\left({\left(A+BK\right)}^{T}P+P\left(A+BK\right)+\left(\mathit{\alpha }+\mathit{\lambda }\right)P\right)e\left(\mathit{\xi }\right)+\left(\frac{\mathit{\lambda }}{\mathit{\eta }}-\mathit{\mu }\left(1+\mathit{\beta }\right)\right){\tilde{\mathit{\theta }}}^T\tilde{\mathit{\theta }}+\left({\mathit{\alpha }}^{-1}\mathit{\epsilon }-\mathit{\mu }\left(1+{\mathit{\beta }}^{-1}\right)\right)g^T\left(\mathit{\xi }\right)g\left(\mathit{\xi }\right)\le 0$$(51)

By choosing(52) $$\begin{array}{c}\mathit{\mu }=\frac{\mathit{\lambda }}{\mathit{\eta }\left(1+\mathit{\beta }\right)}\hfill \\ \mathit{\epsilon }=\frac{\mathit{\alpha }\mathit{\lambda }}{\mathit{\eta }\mathit{\beta }}\hfill \end{array}$$(52)

Equation (51) is reduced to(53) $$e^T\left(\mathit{\xi }\right)\left({\left(A+BK\right)}^{T}P+P\left(A+BK\right)+\left(\mathit{\alpha }+\mathit{\lambda }\right)P\right)e\left(\mathit{\xi }\right)\le 0$$(53)

By multiplying (38) from left and right by $diag\left\{P\right\}$ and substituting P = X⁻¹ and K = YX⁻¹ we have(54) $$\left({\left(A+BK\right)}^{T}P+P\left(A+BK\right)+\left(\mathit{\alpha }+\mathit{\lambda }\right)P\right)\le 0$$(54)

By multiplying inequity (54) from left by $e^T\left(t\right)$ and from right by $e\left(t\right)$, (51) is obtained. Also, Equation (39) is obtained from (50).

According to Lemma 1, there is a set $\mathrm{\Omega }=\left\{e\in R^{n_x}|S\left(e,\tilde{\mathit{\theta }}\right)\le \frac{\mathit{\mu }{w}_{max}^2}{\mathit{\lambda }}\right\}$ so that it is a robust invariant set for system (33).

By choosing(55) $$\mathit{\alpha }={10}^{-3},\mathit{\beta }={10}^{-7},\mathit{\lambda }={10}^{-3},\mathit{\eta }={10}^3,\mathit{\mu }={10}^{-6},\mathit{\epsilon }={10}^{-2}$$(55)

And solving LMI for (38) and (39) relations, P and K are obtained as(56) $$K=\left[\begin{array}{cc}-0.34\hfill & 1.7\hfill \end{array}\right]$$(56) (57) $$P=\left[\begin{array}{cc}0.58\hfill & 0.18\hfill \\ 0.18\hfill & 0.56\hfill \end{array}\right]\times {10}^{-3}$$(57)

Also, the following controlling algorithm is employed to stabilize the system (Yu, Maier, Chen, & Allgöwer, 2013).

Algorithm

Step 0. At time t₀, set $\overline{x}(t_0)$ in which $x\left(t_0\right)$ is the current state.

Step 1. At time t_k and current state $\left(\overline{x}(t_k),x(t_k)\right)$, solve problem 1 to obtain the nominal control action $\overline{u}(t_k)$ and the actual control action $u(t_k)=\overline{u}(t_k)+Ke(t_k)$.

Step 2. Apply the control $u\left(t_k\right)$ to the system (20), during the sampling interval [t_k, t_k+1], where t_k+1 = t_k + δ.

Step 3. Measure the state x(t_k+1) at the next time instant t_k+1 of the system (20) and compute the successor state $\overline{x}(t_{k+1})$ of the nominal system (24) under the nominal control $\overline{u}(t_k)$.

Step 4. Set $(\overline{x}(t_k),x(t_k))=(\overline{x}(t_{k+1}),x(t_{k+1})),t_k=t_{k+1}$, and go to step 1.

Block diagram of the control system is shown in the following Figure 2.

Finally, suppose that K and Ω are given, and problem (29) is feasible at time ξ₀. Then,

(1)	Problem (29) is feasible for all $\mathit{\xi }>{\mathit{\xi }}_0$,
(2)	According to Algorithm, the trajectory of the system (9) under MPC control law is asymptotically ultimately bounded,
(3)	The closed-loop system is input-to-state stable.

Proof

(Yu et al., 2013), (1) the online optimization is completely independent of the disturbances, because at the next time instant, the only required elements for Problem (29) solution are the “measured” state of the nominal system and the nominal system dynamics. Therefore, if one solution for Problem (29) can be applied at the initial time instant, then one can guarantee the recursive feasibility of the system (Chen & Allgöwer, 1998).

(2) and (3) Since the nominal system is asymptotically stable (Chen & Allgöwer, 1998), there exists a class κ function $\mathrm{\Gamma }(\Vert \overline{x}\Vert ,\mathit{\xi })$ (Khalil, 2002) such that(58) $$\Vert \overline{x}\left(\mathit{\xi }\right)\Vert \le \mathrm{\Gamma }\left(\Vert \overline{x}\left({\mathit{\xi }}_0\right)\Vert ,\mathit{\xi }\right)\forall \mathit{\xi }\ge {\mathit{\xi }}_0$$(58)

Because $S(e,\tilde{\mathit{\theta }})\le \frac{\mathit{\mu }{w}_{max}^2}{\mathit{\lambda }}$ for all ξ ≥ ξ₀ and $e\left(\mathit{\xi }\right)\in \mathrm{\Omega }$, there exists a class function such that(59) $$\Vert e\left(\mathit{\xi }\right)\Vert \le \mathit{\gamma }\left(\begin{array}{c}\underset{{\mathit{\xi }}_0\le \mathit{\tau }\le \mathit{\xi }}{sup}\hfill \end{array}\Vert w(\mathit{\tau })\Vert \right)\forall \mathit{\xi }\ge {\mathit{\xi }}_0$$(59)

Since $\overline{x}(\mathit{\xi })=\overline{x}+e(\mathit{\xi })$ and $\overline{x}({\mathit{\xi }}_0)=x({\mathit{\xi }}_0)$,(60) $$\Vert x\left(\mathit{\xi }\right)\Vert \le \mathrm{\Gamma }(\Vert x\left({\mathit{\xi }}_0\right)\Vert ,\mathit{\xi })+\mathit{\gamma }\left(\begin{array}{c}\underset{{\mathit{\xi }}_0\le \mathit{\tau }\le \mathit{\xi }}{sup}\hfill \end{array}\Vert w(\mathit{\tau })\Vert \right)\forall \mathit{\xi }\ge {\mathit{\xi }}_0$$(60)

So, according to the Algorithm, the system (9) solution, exposed to the MPC control, is asymptotically ultimately bounded and the closed-loop system is input-to-state stable (Khalil, 2002).

4. Simulation

In this section a simulation is conducted to prove the robustness and effectiveness of the presented controller. To do this, the presented controller is compared to an active controller (Sheng et al., 2014b) in order to confirm its efficiency. The simulation scenario: before t = 150 s, the throttle valve output is 0.7. From t = 150 s to t = 300 s the throttle valve output experiences a 20% reduction and reaches to 0.5, which leads to surge. Finally, after t = 300 s, the throttle valve output increases to 0.65. Equations (62–64) show the system disturbances in t = 500 s. Values of the compressor parameters used in this simulation are according to Greitzer (1976a).(61) $$B=1.8,l_{\text{c}}=3,H=0.18,W=0.25,{\mathit{\psi }}_{\text{c0}}=0.3$$(61)

The initial points of process were $\left(x_1\left(0\right),x_2\left(0\right)\right)=\left(0.6,0.6\right)$ in the left side of the surge line and the initial points of model were $\left(x_1\left(0\right),x_2\left(0\right)\right)=\left(1,1\right)$. The optimization problem described by problem (29) is solved in discrete time with a sample time of δ = 0.1 time units and prediction horizon of T_p = 0.3 time units.

In the first mode, it is assumed that some stable disturbances are applied to the system.(62) $$\begin{array}{c}{\text{d}}_{\mathit{\varphi }}\left(\mathit{\xi }\right)=0.5\mathit{\varphi }\hfill \\ {\text{d}}_{\mathit{\psi }}\left(\mathit{\xi }\right)=0.5\mathit{\psi }\hfill \end{array}$$(62)

The system states after applying the active surge controller for compressor flow and pressure are shown in Figures 3 and 4, respectively. It can be seen that the presented controller, in comparison to the controller in Sheng et al. (2014b), can easily stabilize the system.

Figure 3. Compressor pressure.

Figure 4. Compressor flow.

Figure 5 presents the control signal. Since CCV is considered as the controller operator, its output have upper and lower bounds, while the controller (Sheng et al., 2014b) fluctuates continuously during this time. On the other hand, the controller presented in this paper has a smooth control signal.

Figure 5. Control signal.

Figure 6 shows the trajectory of the compressor in a performance curve. It shows how the controller prevents the compressor from entering the surge area.

Figure 6. Compression of the system trajectories.

In the second mode, it is assumed that some transitory disturbances are applied to the system. These disturbances are modeled as(63) $$\begin{array}{c}{\text{d}}_{\mathit{\varphi }}\left(\mathit{\xi }\right)=0.15e^{-0.015\mathit{\xi }}cos\left(0.2\mathit{\xi }\right)\hfill \\ {\text{d}}_{\mathit{\psi }}\left(\mathit{\xi }\right)=0.1e^{-0.005\mathit{\xi }}sin\left(0.3\mathit{\xi }\right)\hfill \end{array}$$(63)

Figures 7 and 8 present compressor flow and pressure, respectively. It is obvious that the presented controller has a better performance in surge avoidance than the controller in (Sheng et al., 2014b).

Figure 7. Compressor pressure.

Figure 8. Compressor flow.

The control signal for compressor surge avoidance is also presented in Figure 9, where the control signal for Sheng et al. (2014b) controller is still a fluctuation value.

Figure 9. Control signal.

The trajectory of the compressor is also shown in Figure 10 and it conveys the ability of the controller for disturbance avoidance.

Figure 10. Compression of the system trajectories.

In the third mode, it is assumed that some stable disturbances are applied to the system. These disturbances are modeled as follows(64) $$\begin{array}{c}{\text{d}}_{\mathit{\varphi }}\left(\mathit{\xi }\right)=0.02sin\left(0.1\mathit{\xi }\right)+0.02cos\left(0.4\mathit{\xi }\right)\hfill \\ {\text{d}}_{\mathit{\psi }}\left(\mathit{\xi }\right)=0.02sin\left(0.1\mathit{\xi }\right)+0.02cos\left(0.4\mathit{\xi }\right)\hfill \end{array}$$(64)

As it can be seen in Figures 11 and 12, the presented controller can stabilize the compressor flow better than the controller in Sheng et al. (2014b).

Figure 11. Compressor pressure.

Figure 12. Compressor flow.

The CCV output is also presented in Figure 13. It can be observed that this output is fluctuation in some time points for (Sheng et al., 2014b) controller.

Figure 13. Control signal.

The performance curve for the compressor is also shown in Figure 14. It can be seen that, despite the existence of disturbance, the controller has been able to avoid surge.

Figure 14. Compression of the system trajectories.

It should be noted that in the controller presented in this paper, like the controller in Sheng et al. (2014b), the compressor characteristic, as well as the disturbance bound are not known. However, the presented controller, unlike the controller in Sheng et al. (2014b), does not have access to throttle valve values. Nevertheless, the presented controller provides a better performance in surge avoidance and compared to the controller presented in Sheng et al. (2014b), it is more robust to uncertainties. Controller Sheng et al. (2014b) does not need Greitzer parameter in the designing phase, but the controller presented in this paper requires this parameter, which leads the controller to lose its robustness against B parameter.

5. Conclusion

Active surge control is one of the most important ways to increase the operating area of a compressor, but most of the active control methods need to recognize the compressor characteristic. In this research a new active surge controller was designed through a robust adaptive fuzzy method with Tube-MPC. Without requiring any information about throttle valve opening and compressor characteristic, this new surge controller avoided surge in an effective way. As such, this controller was robust against flow and pressure bounded disturbances, without considering their upper bounds. The results obtained from simulation show the efficiency and robustness of this controller.

Additional information

Funding

Funding. The authors received no direct funding for this research.

Notes on contributors

Masoud Taleb Ziabari

Masoud Taleb Ziabari is a PhD student in Control and System Engineering in Science and Research branch, Islamic Azad University, Tehran, Iran. His research interests include Nonlinear Control, MPC and Robust Control.

Mohammad Reza Jahed-Motlagh

Mohammad Reza Jahed-Motlagh is a professor at the Department of Electrical Engineering, Iran University of Science and Technology, Tehran, Iran. His research interests include Control Engineering, Complex Systems, Nonlinear System and Artificial Intelligence.

Karim Salahshoor

Karim Salahshoor is a professor at the Department of Automation and Instrumentation, Petroleum University of Technology (PUT), Tehran, Iran. His research interests include Intelligent Control Applications and Advanced Process Control and Intelligent Monitoring, Fault Detection and Diagnosis.

Amin Ramezani

Amin Ramazani is an assistant professor at the Department of Electrical Engineering, Tarbiat Modares University Tehran, Iran. His major research interests are Process Control and Automation and Instrumentation Analysis and Design.

Ali Moarefianpur

Ali Moarefianpur is an assistant professor of Electrical Engineering Department, Science and Research branch, Islamic Azad University, Tehran, Iran. His research interests are Convex Optimization, Industrial Control and Artificial Intelligence.