Improved predictive functional control for ethylene cracking furnace

Abstract

The objective of this paper is to show the design and application of pass temperature balance control system using an improved predictive functional control method in eight 800 tone/year USC ethylene cracking furnaces. The advanced pass temperature balance controller is developed using the proposed method and implemented in proprietary APC-ISYS software, which is connected to Yokogawa distributed control system via an OPC server. The advantage of it lies in the fact that the dynamics of pass temperature with nonlinearity and time delay are described by Takagi–Sugeno model and transformed into time-varying extended state space model, and thus, the proposed controller can regulate pass temperature based on the extended state space formulation. In addition, the control law with a linear iterative form, easily applied to industrial process, is derived. The robust analysis for the set point, input disturbance and output disturbance to the output verifies the ability of tracking and disturbance rejection of the proposed method. Application results from an industrial furnace are shown to be markedly better in terms of lower variability in the outlet temperature of both the passes compared to the current proportional–integral–derivative control scheme.

Introduction

Ethylene cracking furnace is a vital equipment to petro-chemical industry in which a lot of important chemical products, such as propylene, ethylene, butadiene, benzene, dimethylbenzene, and so on, can be obtained by the high-temperature cracking reaction. ¹ Because the pass temperature of the radiation room in the cracking furnace has a direct impact on the cracking effect and the ethylene yield, it must be precisely controlled in order to ensure the stable operation of subsequent equipments for high quality and yield of the expected products.

In general, the traditional pass temperature system of ethylene cracking furnace based on proportional–integral–derivative (PID) ² control strategy is mainly to adjust the pass feed-in flow so that the pass temperature is as identical as possible. Because the total feed-in flow is required to remain constant in industrial production, there is strong dynamic coupling between the pass feed-in flow and pass temperature that leads to the poor effect on pass temperature. Wang and Zheng ³ proposed a difference control technique (DCT), that is, difference method,^4,5 applied in an actual heating furnace. The main idea of it is that the difference of pass temperature between two streams is adjusted to be zero using designed controllers. Based on the above control idea of two-passes balance, more passes can be controlled by reiterative use of the DCT method. However, it seems too boring to control such system at lots of iterative times due to many passes. To solve the existing problem of DCT method, Li et al., ⁶ Wang and Zheng,^7,8 Wang et al. ⁹ and Qian et al. ¹⁰ proposed a differences control technique (DsCT), that is, average method. Its main idea is to balance the pass temperature by the average value of all pass temperature, without the need for the DCT’s reiterative application frequently. But this strategy using average method would lead the controller to be large delay and insensitive. In addition, Wang et al.^11,12 also proposed a switching difference control strategy (SDCT) for parallel streams temperature. Its main ideal is to control two worst passes for all passes which can avoid the feed-in flow valves quite frequently being regulated. Luan and colleagues^13–15 proposed a distributed deviation-based uniform temperature control method on feed heater. This scheme requires the usage of the temperature information of adjacent passes, which is more effective for such cases with too many passes. But this method needs many computations of matrices.

As we know, the conventional PID control is used in the above researches. However, the pass temperature is complex with a number of problems, which is difficult to be solved by conventional PID control, such as time delay, multi-input and multi-output behavior and strong nonlinear performance. Therefore, some advanced control methods^16–20 are proposed to overcome these problems. Zhang et al. ¹⁶ designed a steady state energy pass balance control strategy to maintain the pass outlet temperature to be same for multi-pass fired heater. Rivas-Perez et al. ¹⁷ used a modified Smith predictor method to balance the pass outlet temperature and applied it in a crude oil preheating furnace. It solved the problem of the real dominant time delay and disturbance rejection on temperature uniformity of the furnace. Li et al.^18,19 proposed a pole placement control method combined with a humanoid intelligent multi-modality to balance the pass temperature of an ethylene cracking furnace. In addition, some other advanced control methods have also been used to control the pass temperature uniformity, such as adaptive and predictive method,^1,20–23 coordinate control strategy of the pass temperature and pass liquid level, ²⁴ multivariate and fuzzy control method^25–27 and so on. Ojasvi and Kaistha ²⁸ proposed pass balancing control scheme applying DMC as well as simple PID control. In addition, there are some other advanced control methods such as sliding mode control, ²⁹ adaptive control^30,31 and robust model predictive control.^32,33

Different from previous study, in order to balance pass temperature of the radiation room in the cracking furnace, an improved predictive functional control (IPFC) approach is proposed. The advantage of the pass temperature balance control method is that the dynamics of pass temperature with strong nonlinear time delay, multi-input and multi-output behavior and external disturbance is simultaneously solved. The pass temperature with multivariable nonlinear character is described by Takagi–Sugeno (T-S) model. Moreover, the linear model is acquired by weighting a set of T-S sub-model and then transformed into the extended time-varying state space model that can compensate steady state error where the state variables are constructed by the real-life output variables and current and past incremental input and output variables. Finally, the control law with a simple linear iterative form is derived. The robust analysis for the set point, input disturbance and output disturbance to the output verifies that the proposed approach has ability of disturbance rejection and tracking.

Using the IPFC approach, the advanced pass temperature balance controller is developed and implemented in proprietary APC-ISYS software, which is connected to Yokogawa DCS system via an OPC server. Application results of the pass temperature balance in the ethylene cracking furnace using proposed method compared with recent literature and PID control are presented. The fluctuation of the pass temperature has been decreased from $\pm 5°\mathrm{C}$ to $\pm 1°\mathrm{C}$ . Therefore, the developed pass temperature balance controller can reach the improved control performance.

The rest of this paper is arranged as follows: IPFC algorithm and its corresponding robust analysis are detailed in “Improved predictive functional control” section. “Implementation case” section presents the implementation case. Conclusions are presented in “Conclusion” section.

IPFC

Process description

In general, an industrial process can be described by the following nonlinear equation with p inputs, $\mathit{u}\in U\subset R^p$ , and q outputs, $\mathit{y}\in Y\subset R^q$

$$y_h(k+1)=R_h(\mathit{\nu }(\mathit{k}),\mathit{u}\left(\mathit{k}\right))$$

(1)

where h = 1, 2, …, q and $\mathit{\nu }(\mathit{k})\in {\mathit{R}}^{\mathit{q}}$ consists of current and past outputs as well as past inputs as follows

$$\mathit{u}\left(k\right)={[u_1\left(k\right),u_2\left(k\right),···,u_p\left(k\right)]}^{\mathrm{T}}$$

(2a)

$$\mathit{\nu }(\mathit{k})={[{\mathit{y}}_1\left(\mathit{k}\right),\dots ,{\mathit{y}}_{\mathit{q}}\left(\mathit{k}\right),{\mathit{u}}_1(\mathit{k}-1),\dots ,{\mathit{u}}_{\mathit{p}}(\mathit{k}-1)]}^{\mathrm{T}}$$

(2b)

where ${\mathit{y}}_h(k)=[y_h(k),···,y_h(k-n_{y,h}+1),y_h(k-n_{y,h})]$ and ${\mathit{u}}_l(k-1)=[u_l(k-1),···,u_l(k-n_{u,l}+1),u_l(k-n_{u,l})]$ , $l=1,\dots ,p$ , where n_{u, l} and n_{y, h} refer to the order of the lth input and hth output, separately.

Here, the multivariable nonlinear process described by Equation (1) can be approximated by many of T-S sub-models. The fuzzy rule R_hj is expressed as follows.

R_hj: If ${\nu }_1(k)$ is $Z_{\mathrm{hj},1}$ , ${\nu }_2(k)$ is $Z_{\mathrm{hj},2}$ ,…, ${\nu }_m(k)$ is $Z_{\mathrm{hj},m}$ , and $u_1(k)$ is $Z_{\mathrm{hj},m+1}$ , $u_2(k)$ is $Z_{\mathrm{hj},m+2}$ ,…, $u_p(k)$ is $Z_{\mathrm{hj},m+p}$ , then

$$y_{\mathrm{hj}}(k+1)={\mathit{\mu }}_{\mathit{hj}}\mathit{\nu }(\mathit{k})+{\mathit{\zeta }}_{\mathit{hj}}\mathit{u}\left(\mathit{k}\right)+{\mathit{\theta }}_{\mathit{hj}},\mathit{j}=1,2,···,{\mathit{R}}_{\mathit{h}}$$

(3)

where Z_hj refer to the antecedent fuzzy set of the jth rule and the hth output, ${\mathit{\mu }}_{\mathit{hj}}$ and ${\mathit{\zeta }}_{\mathit{hj}}$ stand for the vectors of consequent parameters, ${\theta }_{\mathrm{hj}}$ is a corrected item R_h, and $m=\sum _{h=1}^q(n_{y,h}+1)+\sum _{l=1}^p{n}_{u,l}$ stands for the number of fuzzy rules of the hth output and antecedent variables, separately. Equation (1) can be weighted by the following linear combination of a group of sub-model (3).

$$y_h(k+1)=\sum _{j=1}^{R_h}{\omega }_{\mathrm{hj}}(\mathit{\nu },\mathit{u})({\mathit{\mu }}_{\mathit{hj}}\mathit{\nu }(\mathit{k})+{\mathit{\zeta }}_{\mathit{hj}}\mathit{u}\left(\mathit{k}\right)+{\mathit{\theta }}_{\mathit{hj}})$$

(4a)

$${\omega }_{\mathrm{hj}}(\mathit{\nu },\mathit{u})=\frac{{\beta }_{\mathrm{hj}}(\mathit{\nu },\mathit{u})}{\sum _{l=1}^{R_h}{\beta }_{\mathrm{hj}}(\mathit{\nu },\mathit{u})}$$

(4b)

$$\sum _{j=1}^{R_h}{\omega }_{\mathrm{hj}}(\mathit{\nu },\mathit{u})=1,0\le {\omega }_{\mathrm{hj}}(\mathit{\nu },\mathit{u})<1$$

(4c)

where ${\beta }_{\mathrm{hl}}(\mathit{\nu },\mathit{u})=\underset{\mathit{\kappa }=1}{\overset{\mathit{m}}{\Pi }}{\mathit{\lambda }}_{{\mathit{Z}}_{\mathit{hl},\mathit{\kappa }}}({\mathit{\zeta }}_{\mathit{\kappa }})\underset{\mathit{\iota }=1}{\overset{\mathit{p}}{\Pi }}{\mathit{\lambda }}_{{\mathit{Z}}_{\mathit{hj},\mathit{\iota }}}({\mathit{u}}_{\mathit{\iota }})$ is the Gaussian membership function of the antecedent variables and $\Pi$ stands for the fuzzy operator.

Online identification of consequent parameters

To eliminate the effect of unknown disturbances and model mismatch, the consequent parameters of T-S models are identified online using the following recursive weighting least square approach.^34–36 Equation (4) is expressed as follows

$$y_h\left(k\right)=\mathit{\Phi }{(\mathit{k}-1)}^{\mathrm{T}}\mathit{\theta }$$

(5)

where

$$\mathit{\Phi }(\mathit{k}-1)={[{\mathit{\omega }}_{\mathit{h}1}\mathit{\phi }(\mathit{k}-1),···,{\mathit{\omega }}_{\mathit{h}{\mathit{R}}_{\mathit{h}}}\mathit{\phi }(\mathit{k}-1)]}^{\mathrm{T}}$$

$$\mathit{\theta }={[{\mathit{\theta }}_1,{\mathit{\theta }}_2,···,{\mathit{\theta }}_{{\mathit{R}}_{\mathit{h}}}]}^{\mathrm{T}}$$

$$\begin{array}{c}\mathit{\phi }(\mathit{k}-1)=\hfill \\ \left[\begin{array}{c}{\mathit{y}}_1(\mathit{k}-1),\dots ,{\mathit{y}}_1(\mathit{k}-{\mathit{n}}_{\mathit{y},1}-1),\dots ,{\mathit{y}}_{\mathit{q}}(\mathit{k}-{\mathit{n}}_{\mathit{y},\mathit{q}}-1),\hfill \\ \multicolumn{1}{c}{{\mathit{u}}_1(\mathit{k}-2),\dots ,{\mathit{u}}_1(\mathit{k}-{\mathit{n}}_{\mathit{u},1}-1),{\mathit{u}}_2(\mathit{k}-2),\dots ,}\\ \multicolumn{1}{c}{{\mathit{u}}_{\mathit{p}}(\mathit{k}-{\mathit{n}}_{\mathit{u},\mathit{p}}-1),{\mathit{u}}_1(\mathit{k}-1),\dots ,{\mathit{u}}_{\mathit{p}}(\mathit{k}-1),1]}\end{array}\right]\hfill \end{array}$$

$${\mathit{\theta }}_{\mathit{j}}={[{\mu }_{\mathrm{hj},1},\dots ,{\mu }_{\mathrm{hj},n_{y,1}},{\mu }_{\mathrm{hj},n_{y,1}+1},\dots ,{\mu }_{\mathrm{hj},\sum _{h=1}^q(n_{y,h}+1)},\dots ,{\mu }_{\mathrm{hj},m},{\zeta }_{\mathrm{hj},1},\dots ,{\zeta }_{\mathrm{hj},p},{\theta }_{\mathrm{hj}}]}^{\mathrm{T}}$$

where $\mathit{\Phi }(\mathit{k}-1)$ is the regression vector and $\mathit{\theta }$ is referred to the vector of consequent parameters that can be identified as follows

$${\mathit{\theta }}_{\mathit{j}}\left(k\right)={\mathit{\theta }}_{\mathit{j}}(k-1)+\frac{W(k-1){\mathit{P}}_j(k-1)\mathit{\Phi }(\mathit{k}-1)}{1+W(k-1)\mathit{\Phi }{(\mathit{k}-1)}^{\mathrm{T}}\mathit{P}(\mathit{k}-1)\mathit{\Phi }(\mathit{k}-1)}\sigma \left(k\right)$$

(6a)

$${\mathit{P}}_j\left(k\right)={\mathit{P}}_j(k-1)-\frac{W(k-1){\mathit{P}}_j(k-1)\mathit{\Phi }(\mathit{k}-1)\mathit{\Phi }{(\mathit{k}-1)}^{\mathrm{T}}\mathit{P}(\mathit{k}-1)}{1+W(k-1)\mathit{\Phi }{(\mathit{k}-1)}^{\mathrm{T}}\mathit{P}(\mathit{k}-1)\mathit{\Phi }(\mathit{k}-1)}$$

(6b)

where the initial values of it is set as $\mathit{\theta }(0{)}^{\mathrm{T}}=\mathit{\sigma }\mathit{I}$ and $\mathit{P}(0)={\gamma }^2\mathit{I}$ , where $\sigma$ is quite small positive real number, γ is the vector of quite a big positive real number, and I is an identity matrix. ${\mathit{P}}_j(k-1)$ is a portion of covariance matrix $\mathit{P}(k-1)=[{\mathit{P}}_1(k-1),···,{\mathit{P}}_{R_h}(k-1){]}^{\mathrm{T}}$ , $\epsilon (k)=[y_h(k)-\mathit{\Phi }(\mathit{k}-1{)}^{\mathrm{T}}\mathit{\theta }(\mathit{k}-1)]$ is the residual error, and $W(k-1)$ is a non-negative weighted coefficient.

Multivariable predictive functional control

Equation (4a) is rewritten as follows

$$y_h(k+1)={\hat{\mathit{\mu }}}_h\mathit{\nu }(\mathit{k})+{\hat{\mathit{\zeta }}}_{\mathit{h}}\mathit{u}\left(\mathit{k}\right)+{\hat{\mathit{\theta }}}_{\mathit{h}}$$

(7)

where

$${\hat{\mathit{\mu }}}_h\left(k\right)=\sum _{j=1}^{R_h}{\omega }_{\mathrm{hj}}(\mathit{\nu },\mathit{u}){\mathit{\mu }}_{\mathit{hj}}$$

(8a)

$${\hat{\mathit{\eta }}}_h\left(k\right)=\sum _{j=1}^{R_h}{\omega }_{\mathrm{hj}}(\mathit{\nu },\mathit{u}){\mathit{\zeta }}_{\mathit{hj}}$$

(8b)

$${\hat{\theta }}_h\left(k\right)=\sum _{j=1}^{R_h}{\omega }_{\mathrm{hj}}(\mathit{\nu },\mathit{u}){\theta }_{\mathrm{hj}}$$

(8c)

Taking a difference operation $\Delta =1-z^{-1}$ on both sides of Equation (7), it can obtain

$$\Delta y_h(k+1)={\hat{\mathit{\mu }}}_h\Delta \mathit{\nu }(\mathit{k})+{\hat{\mathit{\zeta }}}_{\mathit{h}}\mathit{\Delta }\mathit{u}\left(\mathit{k}\right)+\mathit{\Delta }{\hat{\mathit{\theta }}}_{\mathit{h}}$$

(9)

To eliminate steady state error, the incremental state space vector $\Delta \hat{\mathit{X}}(k)$ is given as

$$\begin{array}{c}\Delta \hat{\mathit{X}}\left(k\right)=\hfill \\ {[\Delta y_1\left(k\right),···,\Delta y_1(k-n_{y,1}),···,\Delta y_q(k-n_{y,q}),\Delta u_1(k-1),···,\Delta u_1(k-n_{u,1}),···,\Delta u_m(k-n_{u,m})]}^{\mathrm{T}}\hfill \end{array}$$

Then, the fuzzy models can be transformed into the linear time-varying augmented state space model as follows

$$\begin{array}{c}\Delta \hat{\mathit{X}}(k+1)=\hat{\mathit{A}}\left(k\right)\Delta \hat{\mathit{X}}\left(k\right)+\hat{\mathit{B}}\left(k\right)\Delta \mathit{u}\left(k\right)+\Delta \mathit{\Theta }\left(\mathit{k}\right)\hfill \\ \mathit{\Delta }\mathit{y}\left(\mathit{k}\right)=\hat{\mathit{C}}\left(\mathit{k}\right)\mathit{\Delta }\hat{\mathit{X}}\left(\mathit{k}\right)\hfill \end{array}$$

(10)

where

$$\hat{\mathit{A}}={\left[\begin{array}{cccccccc}{\hat{\mu }}_{1,1}\hfill & {\hat{\mu }}_{1,2}\hfill & \cdots \hfill & \cdots \hfill & \cdots \hfill & \cdots \hfill & {\hat{\mu }}_{1,m-1}\hfill & {\hat{\mu }}_{1,m}\hfill \\ \multicolumn{1}{c}{1}& 0\hfill & \cdots \hfill & \cdots \hfill & \cdots \hfill & \cdots \hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & \cdots \hfill & \cdots \hfill & \cdots \hfill & \cdots \hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \cdots \hfill & \cdots \hfill & \cdots \hfill & \cdots \hfill & ⋮\hfill & ⋮\hfill \\ \multicolumn{1}{c}{{\hat{\mu }}_{q,1}}& {\hat{\mu }}_{q,2}\hfill & \cdots \hfill & \cdots \hfill & \cdots \hfill & \cdots \hfill & {\hat{\mu }}_{q,m-1}\hfill & {\hat{\mu }}_{q,m}\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & \cdots \hfill & 1\hfill & 0\hfill & \cdots \hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & \cdots \hfill & 0\hfill & 1\hfill & \cdots \hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \cdots \hfill & \cdots \hfill & \cdots \hfill & \cdots \hfill & ⋮\hfill & ⋮\hfill \end{array}\right]}_{m\times m}$$

$$\hat{\mathit{B}}={\left[\begin{array}{ccccc}{\hat{\zeta }}_{1,1}\hfill & {\hat{\zeta }}_{1,2}\hfill & \cdots \hfill & \cdots \hfill & {\hat{\zeta }}_{1,p}\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & \cdots \hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \hfill & \hfill & ⋮\hfill \\ \multicolumn{1}{c}{{\hat{\zeta }}_{q,1}}& {\hat{\zeta }}_{q,2}\hfill & \cdots \hfill & \cdots \hfill & {\hat{\zeta }}_{q,p}\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & \cdots \hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \hfill & \hfill & ⋮\hfill \\ \multicolumn{1}{c}{1}& 0\hfill & \cdots \hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & \cdots \hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \hfill & \hfill & ⋮\hfill \end{array}\right]}_{m\times p}$$

$$\hat{\mathit{C}}={\left[\begin{array}{ccccccccc}1\hfill & 0\hfill & ···\hfill & ···\hfill & ···\hfill & ···\hfill & ···\hfill & ···\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & ···\hfill & 1\hfill & 0\hfill & ···\hfill & ···\hfill & ···\hfill & 0\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & ···\hfill & \hfill & \hfill & \hfill & \hfill & \hfill & ⋮\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & ···\hfill & ···\hfill & ···\hfill & 1\hfill & 0\hfill & ···\hfill & 0\hfill \end{array}\right]}_{q\times m}$$

$\mathit{\Theta }={{\left[\begin{array}{ccccccccc}{\hat{\mathit{\theta }}}_1\hfill & 0\hfill & ···\hfill & ···\hfill & {\hat{\mathit{\theta }}}_{\mathit{q}}\hfill & 0\hfill & ···\hfill & ···\hfill & 0\hfill \end{array}\right]}^{\mathrm{T}}}_{\mathit{m}\times 1}$

The output variable $\mathit{y}(k)$ is extended into the state variable $\Delta \overline{\mathit{X}}(k)$ as follows

$$\hat{\mathit{X}}\left(k\right)={[\Delta \hat{\mathit{X}}\left(k\right),\mathit{y}\left(k\right)]}^{\mathrm{T}}$$

where

$$\mathit{y}\left(k\right)={\left[\begin{array}{ccc}{y}_1\left(k\right)\hfill & ···\hfill & y_q\left(k\right)\hfill \end{array}\right]}^{\mathrm{T}}$$

Note that

$$\begin{array}{c}\mathit{y}(k+1)-\mathit{y}\left(k\right)=\hat{\mathit{C}}\left(k\right)\hat{\mathit{A}}\left(k\right)\Delta \hat{\mathit{X}}\left(k\right)\hfill \\ +\hat{\mathit{C}}\left(k\right)\hat{\mathit{B}}\left(k\right)\Delta \mathit{u}\left(k\right)+\hat{\mathit{C}}\left(k\right)\Delta \mathit{\theta }\left(\mathit{k}\right)\hfill \end{array}$$

(11)

Putting together Equations (10) and (11) leads to the extended state space model as follows

$$\begin{array}{c}\left[\begin{array}{c}\Delta \hat{\mathit{X}}(k+1)\hfill \\ \multicolumn{1}{c}{\mathit{y}(k+1)}\end{array}\right]=\left[\begin{array}{cc}\hat{\mathit{A}}\left(k\right)\hfill & {\mathit{0}}_{m\times q}\hfill \\ \multicolumn{1}{c}{\hat{\mathit{C}}\left(k\right)\hat{\mathit{A}}\left(k\right)}& {\mathit{I}}_q\hfill \end{array}\right]\left[\begin{array}{c}\Delta \hat{\mathit{X}}\left(k\right)\hfill \\ \multicolumn{1}{c}{\mathit{y}\left(k\right)}\end{array}\right]\hfill \\ +\left[\begin{array}{c}\hat{\mathit{B}}\left(k\right)\hfill \\ \multicolumn{1}{c}{\hat{\mathit{C}}\left(k\right)\hat{\mathit{B}}\left(k\right)}\end{array}\right]\Delta \mathit{u}\left(k\right)+\left[\begin{array}{c}\Delta \mathit{\Theta }\left(\mathit{k}\right)\hfill \\ \multicolumn{1}{c}{\hat{\mathit{C}}\left(k\right)\Delta \mathit{\Theta }\left(\mathit{k}\right)}\end{array}\right]\hfill \\ \mathit{y}\left(k\right)=[{\mathit{0}}_{q\times m},{\mathit{I}}_q]\left[\begin{array}{c}\Delta \hat{\mathit{X}}\left(k\right)\hfill \\ \multicolumn{1}{c}{\mathit{y}\left(k\right)}\end{array}\right]\hfill \end{array}$$

(12)

Moreover, Equation (12) is transformed as

$$\begin{array}{c}\hat{\mathit{X}}(k+1)={\hat{\mathit{A}}}_1\left(k\right)\hat{\mathit{X}}\left(k\right)+{\hat{\mathit{B}}}_1\left(k\right)\Delta \mathit{u}\left(k\right)+{\hat{\mathit{\Theta }}}_1\left(k\right)\hfill \\ \mathit{y}\left(k\right)={\hat{\mathit{C}}}_1\left(k\right)\hat{\mathit{X}}\left(k\right)\hfill \end{array}$$

(13)

where ${\hat{\mathit{A}}}_1(k)=\left[\begin{array}{cc}\hat{\mathit{A}}(k)\hfill & {\mathit{0}}_{m\times q}\hfill \\ \multicolumn{1}{c}{\hat{\mathit{C}}(k)\hat{\mathit{A}}(k)}& {\mathit{I}}_q\hfill \end{array}\right]$ , ${\hat{\mathit{B}}}_1(k)=\left[\begin{array}{c}\hat{\mathit{B}}(k)\hfill \\ \multicolumn{1}{c}{\hat{\mathit{C}}(k)\hat{\mathit{B}}(k)}\end{array}\right]$ , ${\hat{\mathit{\Theta }}}_1(k)=\left[\begin{array}{c}\Delta \mathit{\Theta }(\mathit{k})\hfill \\ \multicolumn{1}{c}{\hat{\mathit{C}}(k)\Delta \mathit{\Theta }(\mathit{k})}\end{array}\right]$ , and ${\hat{\mathit{C}}}_1(k)=[{\mathit{0}}_{q\times m},{\mathit{I}}_q]$ .

The extended state space model (13) used the predictive model to predict the output of process. $\hat{\mathit{X}}(k+\overline{h}|k)$ is the $\overline{h}-\mathrm{step}$ predictive state output. A multiple-step linearization strategy is adapted in the whole predictive horizon. Based on Equation (13), the $\overline{h}-\mathrm{step}$ ahead output predictive vector is presented as follows

$$\begin{array}{c}\mathit{y}(k+\overline{h}|k)={\hat{\mathit{C}}}_1{{\hat{\mathit{A}}}_1}^{\overline{h}}\hat{\mathit{X}}\left(k\right)+{\hat{\mathit{C}}}_1{{\hat{\mathit{A}}}_1}^{\overline{h}-1}{\hat{\mathit{B}}}_1\Delta \mathit{u}\left(k\right)+···\hfill \\ {\hat{\mathit{C}}}_1{\hat{\mathit{A}}}_1{\hat{\mathit{B}}}_1\Delta \mathit{u}(k+\overline{h}-2)+{\hat{\mathit{C}}}_1{\hat{\mathit{B}}}_1\Delta \mathit{u}(k+\overline{h}-1)\hfill \\ +{\hat{\mathit{C}}}_1{{\hat{\mathit{A}}}_1}^{\overline{h}-1}{\hat{\mathit{\Theta }}}_1+···+{\hat{\mathit{C}}}_1{\hat{\mathit{A}}}_1{\hat{\mathit{\Theta }}}_1+{\hat{\mathit{C}}}_1{\hat{\mathit{A}}}_1{\hat{\mathit{\Theta }}}_1\hfill \end{array}$$

(14)

Using the idea of PFC, the control law can be acquired by linearly weighing a number of basic functions. Here, we adopt the step basic function. The control law is

$$\mathit{u}(k+\overline{h})=\mu =\mathit{u}\left(k\right)$$

(15)

Equations (14) and (15) can be transformed into

$$\begin{array}{c}\multicolumn{1}{c}{\mathit{y}(k+\overline{h}|k)={\hat{\mathit{C}}}_1{{\hat{\mathit{A}}}_1}^{\overline{h}}\hat{\mathit{X}}\left(k\right)+{\hat{\mathit{C}}}_1({{\hat{\mathit{A}}}_1}^{\overline{h}-1}+···+{\hat{\mathit{A}}}_1+\mathit{I})}\\ \times ({\hat{\mathit{B}}}_1\Delta \mathit{u}\left(k\right)+{\hat{\mathit{\Theta }}}_1)\hfill \\ ={\hat{\mathit{C}}}_1{{\hat{\mathit{A}}}_1}^{\overline{h}}\hat{\mathit{X}}\left(k\right)+{\hat{\mathit{C}}}_1({{\hat{\mathit{A}}}_1}^{\overline{h}}-\mathit{I}){({\hat{\mathit{A}}}_1-\mathit{I})}^{-1}\hfill \\ \times ({\hat{\mathit{B}}}_1\Delta \mathit{u}\left(k\right)+{\hat{\mathit{\Theta }}}_1)\hfill \end{array}$$

(16)

To handle the problem of imprecision of model and some unknown disturbances, the predictive error is usually employed to modify the vector of predictive output as follows

$$\mathit{e}\left(k\right)=\overline{\mathit{y}}\left(k\right)-\mathit{y}\left(k\right)$$

(17)

where $\overline{\mathit{y}}(k)$ and $\mathit{y}(k)$ are the actual and predictive output vector of model, separately.

Then, the vector of modified predictive output at time $k+\overline{h}$ is as follows

$${\mathit{y}}_{\mathrm{c}}(k+\overline{h}|k)=\mathit{y}(k+\overline{h}|k)+\mathit{e}\left(k\right)$$

(18)

The reference trajectory, which is a predefined path we hope to track the expected set-value, is defined as

$${\mathit{y}}_r(k+\overline{h})=\mathit{c}\left(k\right)-{\mathit{\alpha }}^{\overline{h}}(\mathit{c}\left(k\right)-\overline{\mathit{y}}\left(k\right))$$

(19)

where ${\mathit{y}}_r(k+\overline{h})=\left[\begin{array}{c}{y}_{r,1}(k+\overline{h})\hfill \\ \multicolumn{1}{c}{y_{r,2}(k+\overline{h})}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{y_{r,q}(k+\overline{h})}\end{array}\right]$ , $\mathit{c}(k)=\left[\begin{array}{c}{c}_1(k)\hfill \\ \multicolumn{1}{c}{c_2(k)}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{c_q(k)}\end{array}\right]$ , $\overline{\mathit{y}}\left(k\right)={\left[\begin{array}{cccc}{\overline{y}}_1\left(k\right)\hfill & {\overline{y}}_2\left(k\right)\hfill & ···\hfill & {\overline{y}}_q\left(k\right)\hfill \end{array}\right]}^{\mathrm{T}}$ , and ${\mathit{\alpha }}^{\overline{h}}=\mathrm{diag}[{{\alpha }_1}^{\overline{h}}{{\alpha }_2}^{\overline{h}}···{{\alpha }_q}^{\overline{h}}]$ , where c_h(k) is the expected set-value, ${\overline{y}}_h(k)$ is the $\overline{h}\mathrm{th}$ actual process output at time instant k, ${\alpha }_i=e^{-T_{\mathrm{s}}/T_{r,\overline{h}}}$ , T_s is the sampling time, and $T_{r,\overline{h}}$ is the time constant of the reference trajectory.

Generally speaking, the following performance index is used to obtain the predictive control law of system

$$\begin{array}{c}\multicolumn{1}{c}{J_p=\sum _{\overline{h}=H_1}^{H_2}{{\Vert {\mathit{y}}_r(k+\overline{h})-{\mathit{y}}_c(k+\overline{h}|k)\Vert }_{\mathit{Q}}}^2}\\ +\sum _{\overline{h}=1}^{H_u}{{\Vert \Delta \mathit{u}(k+\overline{h}-1)\Vert }_{\mathit{R}}}^2\hfill \end{array}$$

(20)

where Q  > 0, R  > 0, H₁ and H₂ are the lower and upper limit of prediction horizon separately, and H_u is control horizon.

Based on the principle of the single-value optimization, the performance (Equation (20)) is simplified

$$J_p=\Vert {\mathit{y}}_r(k+H)-{\mathit{y}}_c(k+H|k){\Vert }_{\mathit{Q}}^2+\Vert \Delta \mathit{u}\left(k\right){\Vert }_{\mathit{R}}^2$$

(21)

The incremental control input at current time is derived by substituting Equations (17)–(19) into Equation (21) as follows

$$\begin{array}{c}\multicolumn{1}{c}{\Delta \mathit{u}\left(k\right)={(\mathit{R}+{{\hat{\mathit{B}}}_1}^{\mathrm{T}}{\mathit{F}}_2^{\mathrm{T}}\mathit{Q}{\mathit{F}}_2{\hat{\mathit{B}}}_1)}^{-1}{{\hat{\mathit{B}}}_1}^{\mathrm{T}}{\mathit{F}}_2^{\mathrm{T}}\mathit{Q}·}\\ [(\mathit{I}-{\mathit{\alpha }}^H)(\overline{\mathit{y}}\left(k\right)+\mathit{c}\left(k\right))+\mathit{y}\left(k\right)-{\mathit{F}}_1\hat{\mathit{X}}\left(k\right)-{\mathit{F}}_2{\hat{\mathit{\Theta }}}_1]\hfill \end{array}$$

(22)

where ${\mathit{F}}_1={\hat{\mathit{C}}}_1{{\hat{\mathit{A}}}_1}^H$ and ${\mathit{F}}_2={\hat{\mathit{C}}}_1({{\hat{\mathit{A}}}_1}^H-\mathit{I})({\hat{\mathit{A}}}_1-\mathit{I}{)}^{-1}$ .

The optimal control law at current time is

$$\mathit{u}\left(k\right)=\Delta \mathit{u}\left(k\right)+\mathit{u}(k-1)$$

(23)

where $\mathit{u}(k)=[u_1(k),···,u_p(k){]}^{\mathrm{T}}$ .

The detailed algorithm flow of IPFC using T-S model is as follows:

1. Set the initial parameters: the definite matrix $\mathit{P}(0)$ , the initial consequent parameters $\mathit{\theta }(0)$ of T-S model, the initial control input $\mathit{u}(0)$ , the parameters $T_s$ , $T_r$ , and $H$ of IPFC;

2. Read the process output $\overline{\mathit{y}}(k)$ and the desired output $\mathit{c}(k)$ ;

3. According to Equation (6), estimate the consequence parameters of T-S model: ${\mu }_{\mathrm{hj},1},···,{\mu }_{\mathrm{hj},m},{\zeta }_{\mathrm{hj},1},···,{\zeta }_{\mathrm{hj},p},{\theta }_{\mathrm{hj}}$ , $j=1,···,R_l$ , $h=1,···,q$ .

4. According to Equation (8), calculate ${\hat{\mu }}_{h,1},···,{\hat{\mu }}_{h,m}$ , ${\hat{\zeta }}_{h,1},···,{\hat{\zeta }}_{h,p}$ , and ${\hat{\theta }}_1,···,{\hat{\theta }}_q$ , then calculate the matrix $\hat{\mathit{A}}$ , $\hat{\mathit{B}}$ , $\hat{\mathit{C}}$ , and $\hat{\mathit{\Theta }}$ ;

5. If $k\le 10$ , then loop to step 6, or loop to step 7;

6. Set $\mathit{u}(k)=\mathit{u}(0)$ then loop 8;

7. Calculate the optimal control input $\mathit{u}(k)$ using Equation (23), loop to step 8;

8. Calculate ${\mathit{y}}_c(k+1)$ and implement $\mathit{u}(k)$ to the process;

9. Let $k=k+1$ , return to step 2.

Robust analysis

Setting the noise ${\hat{\mathit{\Theta }}}_1=0$ , Equation (22) is rewritten as

$$\Delta \mathit{u}\left(k\right)={(\mathit{R}+{{\hat{\mathit{B}}}_1}^{\mathrm{T}}{\mathit{F}}_2^{\mathrm{T}}\mathit{Q}{\mathit{F}}_2{\hat{\mathit{B}}}_1)}^{-1}{{\hat{\mathit{B}}}_1}^{\mathrm{T}}{\mathit{F}}_2^{\mathrm{T}}\mathit{Q}·[{\mathit{y}}_r-\mathit{e}-{\mathit{F}}_1\hat{\mathit{X}}\left(k\right)]$$

(24)

The control law is

$$\Delta \mathit{u}\left(k\right)=K_r{\mathit{Y}}_r-K_x\hat{\mathit{X}}\left(k\right)-K_e\mathit{e}$$

(25)

where

$$K_r={(\mathit{R}+{{\hat{\mathit{B}}}_1}^{\mathrm{T}}{\mathit{F}}_2^{\mathrm{T}}\mathit{Q}{\mathit{F}}_2{\hat{\mathit{B}}}_1)}^{-1}{{\hat{\mathit{B}}}_1}^{\mathrm{T}}{\mathit{F}}_2^{\mathrm{T}}\mathit{Q}$$

$$K_x={(\mathit{R}+{{\hat{\mathit{B}}}_1}^{\mathrm{T}}{\mathit{F}}_2^{\mathrm{T}}\mathit{Q}{\mathit{F}}_2{\hat{\mathit{B}}}_1)}^{-1}{{\hat{\mathit{B}}}_1}^{\mathrm{T}}{\mathit{F}}_2^{\mathrm{T}}\mathit{Q}{\mathit{F}}_1$$

$$K_e={(\mathit{R}+{{\hat{\mathit{B}}}_1}^{\mathrm{T}}{\mathit{F}}_2^{\mathrm{T}}\mathit{Q}{\mathit{F}}_2{\hat{\mathit{B}}}_1)}^{-1}{{\hat{\mathit{B}}}_1}^{\mathrm{T}}{\mathit{F}}_2^{\mathrm{T}}\mathit{Q}$$

For notational simplicity, the robustness of the system is verified for single-input and single-output (SISO) case. The process is described by Equation (1) and then is transformed into Equation (7) based on T-S models. The model (Equation (7)) can be expressed as the transfer function interpretation. The transfer function of the actual process is $B_n(z)/A_d(z)$ and the process model is ${\hat{B}}_n(z)/{\hat{A}}_d(z)$ . The model order is n and the extended state variable $\Delta \hat{X}(k)$ is 2n. The control gain vector of state feedback is

$$\begin{array}{c}{K}_x=[k_1,k_2,···,k_n,k_{n+1},k_{n+2},···k_{2n}]\hfill \\ K_r=[k_{s1},k_{s2},···,k_{\mathrm{sH}}]\hfill \\ K_e=[k_{s1},k_{s2},···,k_{\mathrm{sH}}]\hfill \end{array}$$

(26)

Define the two corresponding polynomial functions

$$\begin{array}{c}N\left(z\right)=k_1+k_2{z}^{-1}+···+k_n{z}^{-(n-1)}\hfill \\ M\left(z\right)=1+k_{n+1}{z}^{-1}+···+k_{2n-1}{z}^{-(n-1)}\hfill \end{array}$$

(27)

Then from Equation (25), the corresponding polynomial form in Figure 1 is viewed as follows

$$\begin{array}{c}(1-z^{-1})M\left(z\right)u\left(z\right)=W_s{y}_r\left(z\right)-(1-z^{-1})N\left(z\right)Y\left(z\right)\hfill \\ -k_{2n}\hat{y}\left(z\right)-W_{s}e\left(z\right)\hfill \end{array}$$

(28)

Figure 1.

The block diagram of the polynomial form.

The block diagram with respect to Equation (28) is depicted in Figure 1, in which $y_r(z)$ is the set-value, $D_i(z)$ is the input disturbance, and $D_o(z)$ is the output disturbance. The tracking of the set point, the output disturbance, and the input disturbance to the output is discussed to verify the robustness. First, the transfer function is presented from the set point to the output

$$T\left(z\right)=\frac{\overline{y}\left(z\right)}{y_r\left(z\right)}=\frac{W_s{\hat{B}}_n\left(z\right)A_d\left(z\right)}{\Gamma \left(z\right)}$$

(29)

where

$$\begin{array}{c}\multicolumn{1}{c}{\Gamma \left(z\right)=(1-z^{-1})[M\left(z\right)A_d\left(z\right){\hat{A}}_d\left(z\right)+{\hat{A}}_d\left(z\right)B_n\mathrm{N}\left(z\right)]}\\ +{\hat{A}}_d\left(z\right)B_n{k}_{2n}-{\hat{A}}_d\left(z\right)B_n{W}_s+W_s{\hat{B}}_n\left(z\right)A_d\left(z\right)\hfill \end{array}$$

From Equation (24) together with Equation (12), we can easily acquire that $W_s=k_{2n}$ . Therefore, the result is

$$\underset{z\to 1}{lim}\mathrm{T}\left(z\right)=1$$

(30)

From Equation (30), the system can track a constant set value without steady error.

Second, the transfer function is presented from the output disturbance to the output

$$S_{\mathrm{output}}\left(z\right)=\frac{\overline{y}\left(z\right)}{D_o\left(z\right)}=\frac{(1-z^{-1})M\left(z\right)A_d\left(z\right){\hat{A}}_d\left(z\right)}{\Gamma \left(z\right)}$$

(31)

It is

$$\underset{z\to 1}{lim}{S}_{\mathrm{output}}\left(z\right)=0$$

(32)

Thus, the constant output disturbance can be rejected.

Third, the transfer function from the input disturbance to the output is

$$S_{\mathrm{input}}\left(z\right)=\frac{\overline{y}\left(z\right)}{D_i\left(z\right)}=\frac{(1-z^{-1}){\hat{B}}_n\left(z\right)A_d\left(z\right)M\left(z\right)}{\Gamma \left(z\right)}$$

(33)

It can be obtained as follows

$$\underset{z\to 1}{lim}{S}_{\mathrm{input}}=0$$

(34)

Therefore, the constant input disturbance can be rejected.

Through the above roust analysis, we can conclude that the proposed approach has a better capability of tracking and disturbance rejection.

Implementation case

Process description

The process of the USC ethylene cracking furnace is presented in Figure 2. There are eight 800 tone/year USC ethylene cracking furnace (F1110-F1180) in practice. The structure of radiation room for USC cracking furnace is designed as two chambers labeled A and B, respectively. TI112018 and TI112019 are referred to temperature measurements of chambers A and B. Each chamber consists of two groups of pipes. Each pipe is composed of 44 furnace tubes. For chamber A of the cracking furnace (F1120), the feedstock (NAP) is divided into two passes (FIC112102 and FIC112302) and sent into the convection room of F1120 to be preheated. The first preheated temperature (HC PREHEAT1) is about 183°C and the second (HC PREHEAT2) about 396°C. The preheated dilution steam (DS) is injected into each pass and then the mixed substance is sent to radiation section of F1120 where a number of important chemical products, such as propylene, ethylene, butadiene, benzene, dimethylbenzene, and so on, can be obtained by the high-temperature cracking reaction. For the pass flow system, DS and pass flow are mixed using a double closed-loop ratio controller where controlled variable (CV) is the set-point of feed flow. They adopt PID control method. Likewise, there is the same process flow for other two branches (FIC112202 and FIC112402) of chamber B. The final objective is to adjust the pass feed flow such that the pass temperature (TIC112126A and TIC112326A in the south, and TIC112226A and TIC112426 in the north) is as identical as possible.

Figure 2.

Process of the USC ethylene cracking furnace.

Pass temperature balance control

To reduce the influence of the ethylene cracking furnace to the consequence units, the total feed flow has to remain unchanged. Also, there are two groups of tubular pipes (two passes), and the pass temperature is balanced by the change of the pass feed flow. The pass temperature balance control between the two passes becomes a multi-input multi-output control problem. Figure 3 presents a method where two outlet temperatures serve as the CVs, and the set-values of two pass flows are considered as the manipulate variables (MVs). $\mathrm{FI}{\mathrm{C}}_i,i=1,2$ are two feed flow controllers (PI control) of two inlet flows, and $\mathrm{T}{\mathrm{I}}_i,i=1,2$ are two thermocouples. The expression of the pass temperature balance is presented as follows

$$\underset{\Delta F_i}{minJ}=\sum _{i=1}^2|T_i-T|$$

(35)

$$\sum _{i=1}^2\Delta F_i=0$$

(36)

$$F_{i,\min }\le F_i\le F_{i,\max }$$

(37)

$$\Delta F_{i,\min }\le \Delta F_i\le \Delta F_{i,\max }$$

(38)

where T_i (i = 1,2) is the pass temperature; T is referred to the average temperature; and the constraints of the two pass flows, that is, the lower and upper limits of the adjustments ( $\Delta F_{i,\min }$ , $\Delta F_{i,\max }$ ), are applied in the pass temperature balance control. Also the fluctuation of the total feed flow should remain constant. The adjustments of two pass flows ( $\Delta F_i,i=1,2$ ) are the outputs of the pass temperature balance controller (IPFC).

Figure 3.

Pass temperature balance control.

Control strategy

The pass temperature of the ethylene cracking furnace is a complex process dynamic with multi-input and multi-output, nonlinearity and large delay as well as many of unknown disturbances. The traditional PID controllers have a good control performance, which can rapidly restrain the unknown disturbances for process, but PID controller has been proven incapable and insufficient because of the aforementioned complex performance. Therefore, a transparent architecture, that is, IPFC-PID, is used to achieve the temperature balance control as shown in Figure 4. The advanced control level is the IPFC controller which guarantees the capability of tracking and robustness of process, which the conventional control level in the traditional PID controller that rapidly eliminates the unknown disturbances. The advanced strategy implements better performance than the traditional PID control.

Figure 4.

Principle diagram of IPFC-PID transparent control.

The control block diagram is shown in Figure 4, in which $y_{r1}\mathrm{and}{y}_{r2}$ are the outputs of the reference trajectory and ${\overline{y}}_1\mathrm{and}{\overline{y}}_2$ are actual outputs for pass temperature, respectively. The actual process outputs are measured by the field instruments. These output signals are sent to the IPFC controller through transmitters. Then IPFC gives the set-value of feed flow PID ( $\mathrm{PID}1\mathrm{and}\mathrm{PID}2$ ). The outputs of IPFC are computed by Equation (23) and PID controller adjusts the control parameters to track the expected value. The output of PID control is sent to each actuator (Actuator1 and Actuator2) to control the controlled process of each loop ( $\mathrm{Process}1\mathrm{and}\mathrm{Process}2$ ). The IPFC controller proposed in this paper adopts a generalized process, including PID controller, actuators and processes.

Hardware system architecture of advanced control system

The control system for the ethylene cracking furnace in the PetroChina Sichuan Petrochemical Company is CENTUM VP distributed control system (DCS) made by Yokogawa Company. ¹ A control loop is constructed by field control station (FCS). The basic level uses DCS to achieve traditional PID control because of its capability of unknown disturbance rejection. The hardware architecture of system is presented in Figure 5. The upper computer adopts proprietary APC-ISYS software, while the lower computer is CENTUM VP DCS system. The ethylene cracking furnace is connected with the DCS system through a cable. The communication between the upper computer and the lower computer is achieved by OPC server. The proposed IPFC method is written in VB-Script of the upper computer. Here, the PID control is used in conventional control mode. In other words, the advanced control loop is open. The advanced controller is joined when the switching logic between the upper computer and the lower computer is triggered.

Figure 5.

Hardware system architecture.

Industrial implementation results

The pass temperature balance control system using the proposed control algorithm IPFC has been in operation since 2014 in eight 800 k tone/year ethylene cracking furnaces in PetroChina Sichuan Petrochemical Company.

To reach a good control performance, T-S fuzzy model of pass temperature is structured through Equation (3) using the actual process data. G-K clustering approach ³⁷ is used to acquire the premise fuzzy sets and its corresponding parameters. The consequence parameters are online identified by Equation (6). The identified T-S model of pass temperature is

R ₁: If y₁(k) is Z₁₁, y₂(k) is Z₁₂, u₁(k) is Z₁₃, and u₂(k) is Z₁₄, then

$$y_1(k+1)=a_{11}^1{y}_1\left(k\right)+a_{12}^1{y}_2\left(k\right)+b_{11}^1{u}_1\left(k\right)+b_{12}^1{u}_2\left(k\right)$$

$$y_2(k+1)=a_{21}^1{y}_1\left(k\right)+a_{22}^1{y}_2\left(k\right)+b_{21}^1{u}_1\left(k\right)+b_{22}^1{u}_2\left(k\right)$$

R ₂: If y₁(k) is Z₂₁, y₂(k) is Z₂₂, u₁(k) is Z₂₃, and u₂(k) is Z₂₄, then

$$y_1(k+1)=a_{11}^2{y}_1\left(k\right)+a_{12}^2{y}_2\left(k\right)+b_{11}^2{u}_1\left(k\right)+b_{12}^2{u}_2\left(k\right)$$

$$y_2(k+1)=a_{21}^2{y}_1\left(k\right)+a_{22}^2{y}_2\left(k\right)+b_{21}^2{u}_1\left(k\right)+b_{22}^2{u}_2\left(k\right)$$

R ₃: If y₁(k) is Z₃₁, y₂(k) is Z₃₂, u₁(k) is Z₃₃, u₂(k) is Z₃₄, then

$$y_1(k+1)=a_{11}^3{y}_1\left(k\right)+a_{12}^3{y}_2\left(k\right)+b_{11}^3{u}_1\left(k\right)+b_{12}^3{u}_2\left(k\right)$$

$$y_2(k+1)=a_{21}^3{y}_1\left(k\right)+a_{22}^3{y}_2\left(k\right)+b_{21}^3{u}_1\left(k\right)+b_{22}^3{u}_2\left(k\right)$$

(39)

where the Gaussian membership functions of the fuzzy set $Z_{1j}$ , $Z_{2j}$ , and $Z_{3j}$ , $j=1,2,3$ , for the pass temperature is shown in Figure 6.

Figure 6.

Degree of membership of T-S model for the pass temperature: (a) pass temperature for the first branch, (b) pass temperature for the second branch, (c) feed flow for the first branch and (d) feed flow for the second branch.

The fuzzy model (Equation (39)) is employed in the advanced pass temperature controller (IPFC). The parameters of it are $H=15$ , $T_r=0.5\mathrm{s}$ , $T_s=1\mathrm{s}$ , $\mathit{Q}=\mathrm{diag}(0.9,0.9,0.9,0.9,1,1)$ , and $\mathit{R}=\mathrm{diag}(1,1)$ that are acquired by repeated offline test.

The parameters of the corresponding Su and Li ³⁸ are $H=15$ , $T_r=0.5\mathrm{s}$ , $T_s=1\mathrm{s}$ , $\mathit{Q}=\mathrm{diag}(0.9,0.9,0.9,0.9)$ , and $\mathit{R}=\mathrm{diag}(1,1)$ .

Before the advanced controller is operated, traditional PID controller of each loop had to be adjusted separately. The corresponding parameters are acquired by Cohen-Coon tuning approach. ³⁹ Based on these tuning parameters, the engineer adjusts the parameters again by his abundant operational experience. The parameters of feed flow controllers of two passes are $K_{p1}=0.2,K_{i1}=0.013,K_{d1}=0;K_{p2}=0.2,K_{i2}=0.013,K_{d2}=0$ .

The results of the pass temperature before and after the introduction of the proposed IPFC control are shown in Figures 7–14. The control cycle is set to 15 s. The method has been applied in eight 800 k tone/year USC ethylene cracking furnace of PetroChina Sichuan Petrochemical Company. The results of the pass temperature (TIC112126A and TIC112326A in the south) in chamber A are presented in Figures 7, 9, 11 and 13. From these figures, a better control performance can be easily seen with the IPFC compared to the PID, Su and Li ³⁸ and standard PFC. The fluctuation of the pass temperature has been decreased from $\pm 5°\mathrm{C}$ to $\pm 1°\mathrm{C}$ . Hence, the IPFC method surpasses the PID, Su and Li ³⁸ and standard PFC. The corresponding changes of set point of pass feed flow (FIC112102 and FIC112302) are shown in Figures 8, 10, 12 and 14. The similar results on pass temperature (TIC112226A and TIC112426 in the north) for chamber B are obtained and omitted here for the sake of simplicity. Due to the decrease of the fluctuation of the pass temperature, the coking rate in the furnace will be decreased, leading to the interval cycle of decoking prolonged. Therefore, the goals of extending life of the furnace tube have been achieved.

Figure 7.

Pass temperature for chamber A under PID control.

Figure 8.

Set point of pass flow for chamber A under PID control.

Figure 9.

Pass temperature for chamber A under Su and Li. ³⁸

Figure 10.

Set point of pass flow for chamber A under Su and Li. ³⁸

Figure 11.

Pass temperature for chamber A under standard PFC control.

Figure 12.

Set point of pass flow for chamber A under standard PFC control.

Figure 13.

Pass temperature for chamber A under IPFC control.

Figure 14.

Set point of pass flow for chamber A under IPFC control.

Tables 1 and 2 show the results of pass temperature (TIC112126A and TIC112326A in the south) for chamber A, in which the standard deviations of them decrease by 42.73% and 41.79% compared to Su and Li, ³⁸ by 80.45% and 79.25% compared to traditional PID control and by 32.56% and 40.09% compared to standard PFC control. Therefore, compared to Su and Li, ³⁸ PID controller and standard PFC, the developed pass temperature controller can achieve the improved control performance with regard to the set point tracking and disturbance rejection.

Table 1.

Comparison results of pass temperature (TIC112126A).

Method	Average	Max	Min	Extreme deviation	Standard deviation
IPFC	848.0047	849.8350	846.2189	3.6161	0.6656
Wu et al. 32	848.0617	851.7710	842.5104	9.2605	1.1623
Standard PFC	848.0208	852.0275	844.5004	7.5272	0.9869
PID	849.9919	860.2997	843.3713	16.9284	3.4043

Table 2.

Comparison results of pass temperature (TIC112326A).

Method	Average	Max	Min	Extreme deviation	Standard deviation
IPFC	847.9827	850.0922	845.7072	4.3850	0.7455
Wu et al. 32	848.3626	851.8709	842.7717	9.0992	1.2807
Standard PFC	848.0097	852.7457	844.1626	8.5831	1.2444
PID	849.5174	860.1621	843.0208	17.1413	3.5934

Conclusion

The design and application of pass temperature balance control using proposed IPFC algorithm on the ethylene cracking furnace are presented in this paper. The pass temperature system is structured by T-S fuzzy models and then transformed into the extended time-varying state space model. The linear control law is computed by the predictive functional control method. The pass temperature balance controller is implemented in proprietary APC-ISYS software, which is connected to Yokogawa DCS system via an OPC server. The industrial application shows the developed pass temperature balance control system has a better control performance.

Because the uncertainties and unknown disturbance problems are inevitable in many industrial processes, it may lead to the instability of system and poor control performance. How to solve such problems is a great challenge for the control of pass temperature balance. Recent works have demonstrated that the disturbance observer (DOB) is good at addressing the above problems leading to the motivation of integrating to achieve a new model predictive control (MPC) scheme for pass temperature balance, which will be some interesting issues as our future works.