Dynamic Modeling and Nonlinear Predictive ... - ACS Publications

May 16, 2011 - PFC optimization based new PID design for fluidized catalytic cracking unit ... State space model predictive fault-tolerant control for...
1 downloads 0 Views 2MB Size
ARTICLE pubs.acs.org/IECR

Dynamic Modeling and Nonlinear Predictive Control Based on Partitioned Model and Nonlinear Optimization Ridong Zhang,*,† Anke Xue,† and Shuqing Wang‡ † ‡

Information and Control Institute, Hangzhou Dianzi University, Hangzhou 310018, People's Republic of China National Key Laboratory of Industrial Control Technology, Institute of Cyber-Systems & Control, Zhejiang University, Hangzhou 310027, People's Republic of China ABSTRACT: The paper presents a combination modeling procedure and the implementation of a nonlinear predictive control scheme for the optimization of industrial chemical processes. The model structure is first based on a simple step response method. This provides a way to use prior knowledge about the dynamics, which has a general validity, while additional information about the process behavior is derived from measured plant data-model error. This data error driven model framework is applicable for a wide range of chemical operating units under a certain control policy. The same idea is also used to solve the online optimization problem in the predictive controller. The efficiency and effectiveness of the modeling training algorithm and the nonlinear predictive control approach are demonstrated through a coke furnace case study. A good model fitting for the nonlinear plant is obtained by using the new method. A comparison with traditional approaches shows that the new algorithm can considerably reduce modeling error and improve control accuracy.

1. INTRODUCTION Dynamic modeling and control for chemical process optimization requires that models describe the essential dynamic characteristics of the process and control strategies are model relevant and easy to implement. For complex nonlinear chemical processes, to derive a simple and effective model may be difficult.13 A simple data-driven approach may be limited in an industrial environment where the process is subject to noise, while the development of a first principles model may be time-consuming and expensive.48 These drawbacks all pose the difficulty for modeling and nonlinear predictive control. Obviously, a fully mechanistic model can help one to have a better insight into the nonlinear process mechanism. However, developing rigorous models is a difficult and time-consuming task. Moreover, when the final goal is to simulate, monitor, or control the system, such a detailed nonlinear understanding of the system is not suitable for an easy controller design. However, a simple linear control technique often fails to achieve a satisfactory system response if a process exhibits strong nonlinearities. This contradiction leads to an inevitable nonlinear system identification and control.9,10 In general, two main sources of information are available during model building: first principles and process data.1113 And a lot of work has been done. First principles models are said to grasp the dynamics of the process well but are less effective for further model predictive control (MPC) controller designs.14 What’s more, they are generally not readily available in practice due to a lack of detailed physicochemical knowledge of the process. Some examples are in refs 48. For the MPC controllers design, the direct method of developing models based on inputoutput data is recommended.15,16 Some examples are: Greco et al.17 proposed one of the earliest research works reported in this area. In their paper, they proposed an identification method of using separate models for each prediction horizon r 2011 American Chemical Society

for adaptive control to deal with the convergence properties and modeling errors. Rossiter and Kouvaritakis18 proposed another identification algorithm for MPC control by minimizing multisteps ahead of prediction errors with separate models for different prediction horizons. They also proposed a least-squares algorithm for the autoregressive with exogenous (ARX)-type models. One of the disadvantages of these methods is that a large number of model parameters need to be identified. Another is that noise may have a negative impact on the model when the prediction is large. Another category of inputoutput data based MPC-relevant identification can be seen in refs 19 and 20, where they proposed the methods of data prefiltering. Experimental results showed that by prefiltering the inputs and the outputs through a noise model, the minimization algorithm based on multistep prediction error can be reduced to that of an one-step prediction error method. The disadvantage is, as they said, the accuracy of the prediction depends both on the deterministic and the noise mode of the estimated process model. However, in many industrial applications, neither of the above two sources is the sufficient information to build a simple model. This shows that all the available information has to be combined. A good way to accomplish this is by designing hybrid models, in which first principles are combined with data-driven approaches. Or combining both black-box and phenomenological ideas to build a model. These insights are often important to understand the key dynamics in a given process. Some research results have been done in this area.2128 However, the underlying physics is often too complex to derive a suitable phenomenological model. What’s Received: November 1, 2010 Accepted: May 15, 2011 Revised: May 2, 2011 Published: May 16, 2011 8110

dx.doi.org/10.1021/ie102211x | Ind. Eng. Chem. Res. 2011, 50, 8110–8121

Industrial & Engineering Chemistry Research more, the successive control design based on the model is also a difficult task.2934 Some examples of hybrid or partitioned models based MPC can also be seen in representative research papers, such as refs 3540, in which a partitioned model is first derived, and then corresponding MPC controllers are designed. These design methods are reviewed as follows: Doyle et al.35 selected a Volterra model-based predictive control. The controller design is through the local expansion, such as Carleman linearization. This retains the merit of conventional linear controller design and extends it to nonlinear processes. However, it may be erroneous when a global nonlinear behavior is required.36 Maner et al.37 proposed another nonlinear MPC based on the second-order Volterra model, and the control law requires a nonlinear optimization. The main deficiency of Volterra models is that time histories need to be incorporated in the models which sacrifice the computational requirements. Besides, the models may be grossly overparameterized, and the reliability may be questionable.28 Harris and Palazoglu38 used a functional expansion model to design an internal model controller (IMC). This model is limited to fading memory systems, and the controller only gives satisfactory results for a small operation range. Maksumov et al.39 proposed in their paper a global model of neural network-based IMC design, the controller design is through a local linear ARX model. Kalmukale et al.36 and Cheng and Chiu40 proposed the IMC and proportional integral derivative (PID) controllers design based on just-in-time learning (JITL) technique, respectively. These strategies select a linear model plus a nonlinear part to design the control law. The nonlinear part is identified through the JITL technique. The JITL technique selects a fixed first- or second-order model. In this article, a partitioned model with sep-response identification plus a nonlinear optimization and a nonlinear predictive control method are presented. The identification procedure consists of two steps: The first is a step response method to get the general information of the process, and the second step is based on the error between the plant data and the derived step response model. This second error-driven step can be implemented through universal approximation properties of powerful tools such as neural networks, support vector machines, etc.41,42 Unlike a traditional data-driven method, it is to ensure a better accuracy of the identification of nonlinear processes. A further nonlinear predictive control is proposed for the above partitioned model. This predictive control selects the same errordriven strategy and is implemented in a simple linear way. Compared with other nonlinear predictive methods, this allows the easy optimization of nonlinear predictive control based on linear theory for practice.4346 The paper is organized as follows: Section 2 presents the identification approach. Section 3 deals with the nonlinear predictive control method. Section 4 details a case study on the plant data of an industrial coke furnace. Section 5 shows the convergence of the proposed nonlinear predictive control strategy. Conclusion is in Section 6.

2. THE MODELING PROCEDURE The modeling framework to be introduced is called combination of step response and nonlinear optimization. This section will demonstrate the general ideas. The case study in Section4 will discuss this idea in detail.

ARTICLE

2.1. Process Model. The process is composed of a stepresponse linear part and a nonlinear optimization part:

ð1Þ

ym ðkÞ ¼ yL ðkÞ þ yNL ðkÞ T

yL ðkÞ ¼ Φ X,

Φ ¼ ½a1 , a2 , ..., an , b0 , b1 , ..., bm1 

T

X ¼ ½yðk  1Þ , ..., yðk  nÞ, uðt  d  1Þ , ..., uðt  d  mÞT

ð2Þ

yNL ðkÞ ¼ f ðX, EÞ

where X is the inputoutput data sample set, Φh is the identified parameters, yL(k) is the local linear model, yNL(k) is the nonlinear part derived through optimizing E, which is the deviation between yL(k) and the real plant output data y(k), f is a nonlinear function, k is the sampling time instant, and n, m, and d þ 1 are the output, input order, and dead time, respectively. 2.2. Identification Algorithm. This section discusses the twostep modeling of the nonlinear process. 2.2.1. Identification of the Linear Part. The linear model is achieved through the “step-response” method.47 The accurate step-response identification requires that the magnitude of the input signal be reasonably chosen and that several repeated tests be done under the same conditions. After obtaining enough data, certain data-processing methods are chosen to get the model of the process. In this case, several groups of step tests were performed. Wang47 shows the simple method to get the linear first-order, second-order, or higher models which are generally thought as the behavior of industrial processes. First-order plus dead time model : GðsÞ ¼

Keτs Ts þ 1

Let y(t) denote the real-time process output, y(¥) denote the steady value of the process output, and U0 the amplitude of input step signal. The output is first treated as: y*(t) = y(t)/y(¥), then the model gain is calculated as K = y(¥)/U0. Typically for a first-order plus dead time model, y*(t) represents the following form: ( 0 t t1 > τ, and especially when choosing y*(t1) = 0.39, y*(t2) = 0.63, T and τ are calculated as follows: T τ

¼ 2ðt2  t1 Þ ¼ 2t1  t2

Second-order plus dead time model : GðsÞ ¼

ð4Þ Keτs ðT1 s þ 1ÞðT2 s þ 1Þ

For the linear part of the process that is represented as a second-order model, the model gain can also be calculated through K = y(¥)/U0, while τ is directly obtained through the no-response part of the process’s step response. Then remove the no-response part and treat the step response of the process into y*(t) = y(t)/y(¥). Thus the following transfer function will be used to obtain the parameters T1 and T2: GðsÞ ¼ 8111

1 , ðT1 s þ 1ÞðT2 s þ 1Þ

T1 g T2

ð5Þ

dx.doi.org/10.1021/ie102211x |Ind. Eng. Chem. Res. 2011, 50, 8110–8121

Industrial & Engineering Chemistry Research

ARTICLE

Table 1. Relations of n and t1/t2 for Higher Order Processes 1/(Ts þ 1)n

The corresponding response of the above model is 1  yðtÞ ¼

T1 T2 et=T1  et=T2 T1  T2 T1  T2

ð6Þ

Select two points [t1, y*(t1)] and [t2, y*(t2)], and especially when choosing y*(t1) = 0.4, y*(t2) = 0.8, T1 and T2 are calculated through the following equations: 8 T1 T2 > > et1 =T1  et1 =T2 ¼ 0:6 < T1  T2 T1  T2 ð7Þ T1 T2 > > et2 =T1  et2 =T2 ¼ 0:2 : T1  T2 T1  T2 with

8 1 > > ðt1 þ t2 Þ < T1 þ T2  2:16 T1 T2 t1 > > : ðT þ T Þ2  ð1:74t  0:55Þ 2 1 2

ð8Þ

Higher order plus dead time model with model order n : Keτs GðsÞ ¼ ðTs þ 1Þn Wang47 states that when t1/t2 > 0.46, it is recommended to choose a higher-order model. The corresponding parameters t1 and t2 are first decided for y*(t1) = 0.4 and y*(t2) = 0.8, respectively, which is the same as the second-order model derivation case. Then n and T are calculated by the following equation (see Table 1): nT 

t1 þ t2 2:16

ð9Þ

Discussions: • The above procedure describes the linear model derivation. However, how to appropriately select a linear model structure (or order) is a major problem because one- or multistep ahead predictions for the given process will partly be based on it. A typical criterion is provided in ref 45, where t1/t2 is selected as: t1/t2 > 0.32 for first-order plus dead time model, 0.32 < t1/t2 < 0.46 for second-order plus dead time model, and t1/t2 > 0.46 for higher-order plus dead time model. • Thus the typical linear model derivation can be summarized as follows: (1) The model gain K and the dead time τ are decided through K = y(¥)/U0 and the step response of the process. (2) Two time points t1 and t2 corresponding to y*(t1) = 0.4 and y*(t2) = 0.8 are selected. (3) Through the criterion in Table 1, t1/t2 will be calculated to determine the linear model order n. (4) Parameters of the linear part will be calculated through the algorithms in Section 2.2.1. • When the sampling time Ts is selected, the corresponding discrete form of the above continuous models is derived as the linear form yL(k) in eq 1. 2.2.2. Identification of the Nonlinear Part. Define the output error as: e(k) = y(k)  yL(k), the model error is further identified

n

1

2

3

4

5

6

7

8

t1/t2

0.32

0.46

0.53

0.58

0.62

0.665

0.67

0.685

by minimizing the following optimization problem: min E ¼

N

½eðkÞ2 ∑ k¼1

ð10Þ

where N is chosen to ensure the step response of the process reaches its steady state. The optimization problem eq 10 leads to the nonlinear modeling part yNL(k). Denote if NL(X) = yNL(k), then model eq 1 can be rewritten as ym ðkÞ ¼ a1 yðk  1Þ þ 3 3 3 þ an yðk  nÞ þ b0 uðk  d  1Þ þ 3 3 3 þ bm1 uðk  d  mÞ þ NLðXÞ ð11Þ The nonlinear model eq 11 will be used to design an overall nonlinear predictive control strategy. Remarks: • Under the assumption that the real process consists of a linear step response part and a nonlinear optimization part, the quality of the multistep prediction for MPC will depend on both the linear and the nonlinear parts. However, it is not very easy to know the structure of the real process and to obtain an optimal linear structure. The best we can do is to search for a model that is to some degree “optimal” enough for process prediction and controller design. • When a process exhibits nonlinear behavior, the ability of a linear model to grasp its dynamic is limited. What we can do is try to find a relatively suitable linear model. However, for the whole nonlinear identification strategy, linear model mismatch is not very strict since a nonlinear optimization will compensate for its modeling errors. Especially when the final goal is for MPC controller design, the feedback correction in MPC will also help overcome the model mismatch.

3. PREDICTIVE CONTROLLER DESIGN 3.1. Output Prediction. Add the back shift operator Δ = 1  z1

to eq 11, where Δ is an operator to obtain an incremental value of a variable. For example, Δy(k) = (1  z1)y(k) = y(k)  y(k  1), then it can be rewritten as ym ðkÞ ¼ A1 yðk  1Þ þ 3 3 3 þ Anþ1 yðk  n  1Þ þ B1, 0 Δuðk  d  1Þ þ 3 3 3 þ B1, m1 Δuðk  d  mÞ þ ΔNLðXÞ

ð12Þ

where A1 = 1 þ a1, Ai = ai  ai1 (i = 2, 3, ..., n), Anþ1 = an, B1,i = bi(i = 0, 1,..., m  1). The predictive controller will be designed through eq 12. The optimal prediction is composed of three parts: the first part Ypast is determined by past inputs and outputs; the second part GU is determined by present and future inputs; the third part is the prediction error feedback correction, consisting of the nonlinear error E1 and the error E2 caused by external disturbances. Here the elements of E2 are the same. They are the difference between 8112

dx.doi.org/10.1021/ie102211x |Ind. Eng. Chem. Res. 2011, 50, 8110–8121

Industrial & Engineering Chemistry Research

ARTICLE

the real process output and the model output. Note that the prediction error feedback correction here is separated into two parts: E1 and E2. Generally, in the MPC strategy, the error feedback correction consists of only the error between the real process and the model outputs. And this is enough as the correction for model mismatch and impact of unmeasured disturbance.48 However, for the sake of designing the proposed controller through linear theory, the linear part of the nonlinear model will be used for the controller design. Thus a model error will be resulted. This is why another error E1 is introduced into the prediction, and the controller serves a two-fold function of calculating both the correct E1 for the nonlinear process and the optimal control action. In the meantime, the predictive horizon is chosen as p and define Xi ¼ ½yðk þ i  1Þ , ..., yðk þ i  nÞ, uðk þ i  d  1Þ , ..., uðk þ i  d  mÞT ði ¼ 1, 2 , ..., pÞ Then the prediction based on eq 12 is ^ ¼ Y past þ GU þ E1 þ E2 Y

ð13Þ

where ^ = (^y(k þ d þ 1/k), ^y(k þ d þ 2/k), ..., ^y(k þ d þ p/k))T Y Ypast = (ypast(k þ d þ 1), ypast(k þ d þ 2), ..., ypast(k þ d þ p))T U = (Δu(k), Δu(k þ 1), ..., Δu(k þ p  1))T E1 = (ΔNL(X1), ΔNL(X2), ..., ΔNL(Xp))T E2 = (y(k)  ym(k), y(k)  ym(k), ..., y(k)  ym(k))T 2 3 B1, 0 6 7 6 B2, 0 B1, 0 7 0 7 G ¼ 6 6 ::: 7 ::: 4 5 Bp, 0 Bp1, 0 ::: B1, 0 Ypast is the free output of the process which can be calculated by its model eq 12, the elements of G are calculated as follows: B1, 0 ¼ b0

Bi, 0 ¼ bi1 þ

i1

∑ Aj Bij, 0 ,

j¼1

i ¼ 2 , ..., p

E01 ¼ 0 E0 ¼ ðGT G þ λ2~IÞ1 GT ðY r  Y past  E2  E01 Þ ^kÞ Ek1 þ 1 ¼ Ek1 þ δðkÞðY km  Y Uk þ 1 ¼ ðGT G þ λ2~IÞ1 GT ðY r  Y past  E2  Ek1 þ 1 Þ ð19Þ where E01 is the initial value of E1,U0 is the initial value of U, superscript k stands for the value of k step. Ykm is derived by ^ k is derived by substituting Uk into substitutingUk into eq 12, Y eq 13, δ(k) is the convergent factor, it is a positive scalar decreasing sequence, δ(k) ∈ (0,1), limkf¥ δ(k) = 0. When δ(k) is properly chosen, the above recurrence formula is convergent, here δ(k) is chosen as δ(k) = 1/k to ensure the convergence of ^ = Ym, the control law Uk is the optimal the control law. When Y kþ1 one, actually, if ||E1  Ek1|| < a desired tolerance, Uk can be thought as the optimal control law.

4. CONVERGENCE ANALYSIS This section analyzes the performance of the proposed predictive control strategy. It will show that the proposed control strategy is an optimal one. From eq 19, it is seen that when k is becoming larger and larger, δ(k) is becoming smaller and smaller, and the changes of E1k and Uk are becoming smaller and smaller. When Ek1 f E1, where E1 is the true value of Ek1, ^ k f Ykm is derived, thus the optimal control law is derived. So the Y next step is to analyze the convergence of Ek1, if Ek1 f E1, then the control law is convergent. Lemma.49 Consider the following recurrence algorithm: RðkÞ ¼ ΦðkÞ ¼

ð15Þ

Y r ¼ ðyr ðk þ d þ 1Þ, yr ðk þ d þ 2Þ , ..., yr ðk þ d þ pÞÞT ð16Þ

where R(k) is the variable to be estimated, define Ds = {R|all the eigenvalues of A(R) are within the unit circle}, DR is a connected open subset of Ds, and in DR, the functions of eq 20 satisfy the following regular conditions C1∼C5: • C1: Q(k,R,Φ) is Lipschitz continuous of R and Φ near (R, Φ), where R ∈ DR, Φh is arbitrary, and Q(k,R,Φ) is also a continuous differentiable function of R and Φ. • C2: For R ∈ DR, functions A(R) and B(R) are Lipschitz continuous of R. • C3: {e(k)}is the independent random vector sequence, which makes limxf¥ E{Q(k,R,Φ(k,R))} f f* (R); "R ∈ DR, where Φ(k,R) is defined as follows: n Φðk, R Þ ¼ AðRÞΦðk  1, RÞ þ BðRÞeðkÞΦð0, RÞ ¼ 0

where λ2 is the weighting part, from (∂J)/(∂U) = 0 derives U ¼ ðGT G þ λ2~IÞ1 GT ðY r  Y past  E1  E2 Þ

Rðk  1Þ þ δðkÞQ ðk, Rðk  1Þ, ΦðkÞÞ AðRðk  1ÞÞΦðk  1Þ þ BðRðk  1ÞÞeðkÞ ð20Þ

where μ is the smoothing factor and ys is the set point. The vector forms of the reference trajectory and the cost function are set as

^ ÞT ðY r  Y ^ Þ þ λ2 UT Ug J ¼ minfðY r  Y

The above cannot be realized since E1 is not known, so a revised method is added to get the optimal control law:

ð14Þ

The reference trajectory is set as yr ðk þ dÞ ¼ yðkÞ yr ðk þ d þ iÞ ¼ μi yðkÞ þ ð1  μi Þys i ¼ 1, 2 , ..., p

where ~I is matrix of appropriate dimension. Denote qT as the first row of(GTG þ λ2~I)1GT, then the optimized control signal is ð18Þ uðkÞ ¼ uðk  1Þ þ qT ðY r  Y past  E1  E2 Þ

• C4: ∑k=1¥δ(k) = ¥, ∑k=1¥[δ(k)]m < ¥ (m is a positive number, m > 1)

ð17Þ 8113

dx.doi.org/10.1021/ie102211x |Ind. Eng. Chem. Res. 2011, 50, 8110–8121

Industrial & Engineering Chemistry Research

ARTICLE

Figure 1. Process flow of coking furnace (101/3).

• C5: δ(k)is a positive scalar decreasing sequence, with   1 1  lim δðkÞ ¼ 0, lim sup < V ðRD Þ e 0, "RD ∈ DA dτ ð22Þ d > : V ðRD Þ ¼ 0, "RD ∈ Dc , Dc ∈ DA dτ w:p:1

designed by eq 19, and the controller parameters satisfy the following conditions: δ(k) is a positive scalar decreasing sequence δ(k) ∈ (0,1), limkf¥ δ(k) = 0, here it is selected as δ(k) = 1/k, and then the predictive control law eq 19 is convergent. Proof. From eqs 19 and 13:

Then, when k f ¥,RðkÞ fk f¥ Dc . In the meantime, if R* is the overall asymptotic stable equilibrium point, then when k f w:p:1  ¥, R(k) f R* with 1 probability, that is, RðkÞ fk f¥ R Theorem. For nonlinear process eq 1, if its model is treated into the form of eq 12, then the predictive control law is

ð23Þ where, D and F are constant matrixes at each recurrence step: D ¼ δðkÞ½Y km  Y P  GUk  E2  F ¼ diagð1  δðkÞ, 1  δðkÞ , ..., 1  δðkÞÞ Define DS = {F| all the eigenvalues of F are within the unit circle}. It is seen that λk ðFÞ ¼ 1  δðkÞ < 1

k ¼ 1 , ..., p

ð24Þ

which meansDS is the whole plane, that is,DR = DS = Rn. From eq 23 and δ(k) = 1/k, all conditions C1∼C5 can be satisfied. 8114

dx.doi.org/10.1021/ie102211x |Ind. Eng. Chem. Res. 2011, 50, 8110–8121

Industrial & Engineering Chemistry Research

ARTICLE

Thus the adjoint differential equation is dE1D ðτÞ ¼ dτ ¼

f / ðE1D ðτÞÞ lim EfY km  Y past  GUk  E1D ðτÞ  E2 g

k f¥

ð25Þ While Ykm = Ypast þ GUk þ E1 þ E2, then 

f ðE1D ðτÞ ¼ E1  E1D ðτÞ

ð26Þ

Define E1D* as the equilibrium point of differencing equation eq 26, then 

E1D ¼ E1

ð27Þ

Figure 2. Real-time response of the process.

That is, the equilibrium point is at the true value. The Lyapunov function of eq 26 is 1 V ðE1D Þ ¼ ðE1  E1D Þ T ðE1  E1D Þ > 0 2

ð28Þ

then d V ðE1D Þ ¼ dτ ¼

d ðE1  E1D Þ dτ ðE1  E1D ÞT ðE1  E1D Þ e 0 " E1D

ðE1  E1D ÞT

ð29Þ with d V ðE1D Þ ¼ 0, dτ



E1D ¼ E1D ¼ E1

ð30Þ

which means that there exists an invariant set DC = {E1}for the differential equation, and the attraction domain is the whole plane. Thus from the Lemma, E1 f E1. Remarks. • The selection of δ(k) = 1/k ensures the insight into the proposed controller. It can be seen that as k increases, δ(k)decreases. And at each calculation step in a control cycle, a value of E1k is first calculated, and then it is for the successive calculation of control action Uk. For the ideal case that k f þ¥, then δ(k) = 1/k f 0, leading to E1kþ1 = E1k; this means that E1k converges to a fixed value, then the control action Uk is also a fixed one. • The conventional GPC is the special case of the proposed GPC when we select only one step, i.e., just let E10 = 0 and calculate the control action U0 = (GTG þ λ2~I)1 GT (Yr  Ypast  E2  E01) for the process without the following iteration part in eq 19.

5. CASE STUDY This section describes a case study of the above identification and nonlinear predictive control on the real-time data of a coke furnace. Generally speaking, the most common way is to design a linear controller and ignore the nonlinearities. However, there are reasons that nonlinear controllers should be used: (1) when the process to be controlled exhibits a strong degree of nonlinearity and suffers frequent disturbances; and (2) servo problems where operating points change frequently and span a large

Figure 3. Combination modeling results.

range of nonlinear dynamics. The process of industrial coke furnace is a nonlinear process due to the above two main reasons. 5.1. Coke Furnace Process. See Figure 1, it mainly consists of a fractionating tower (T102), a process furnace (101/3), and two coke towers (T101/5,6). The main job of the equipment is to coke residues oil. The process flow is as follows: A flow of residual oil is separated into two branches (FRC8103, FRC8104) and flows into the convection room of the process furnace (101/3). After being heated to about 330 C, the two branches join together and flow out of the radiation room of the furnace and go to the fractionating tower (T102) for heat exchange with oil gas from the top of the coke towers (T101/5,6). After heat exchanging, the heavy part of both residual oil and the oil gas join together (now called circulating oil), the circulating oil is divided into two branches (FRC8107, FRC8108) by pumps (102/1,2,3) and are sent back to enter the radiation room of the furnace (101/3) to be heated to about 495 C. Then the two branches join together and go to the coke towers (T101/5,6) to remove coke. This process is called coking of residues. In this study, one of the control objectives is to keep the output temperature of circulating oil of each branch as close as possible to its desired value, here it is 495 C. The manipulated variables are two fuel flows that are controlled by two valves, respectively. The reasons that generate the nonlinearities of this process are: the flames are not steady, which results in the model uncertainty and the measurement error of the controlled variables. The system is affected by load changes, the amount of input fuel, the 8115

dx.doi.org/10.1021/ie102211x |Ind. Eng. Chem. Res. 2011, 50, 8110–8121

Industrial & Engineering Chemistry Research

ARTICLE

Figure 4. (A) Output signal of proposed under no model mismatch. (B) Input signal of proposed under no model mismatch.

Figure 5. (A) Output signal of GPC under no model mismatch. (B) Input signal of GPC under no model mismatch.

outer environment, and the valves might saturate or stick; these changes bring nonlinearities in the operation. 5.2. Modeling Based on Partitioned Strategy. This section discusses the combination modeling procedure of one of the output temperature loops (as an example, TRC8103 ). For a relatively accurate step response, the magnitude of input signal should be reasonably chosen, and several repeated tests should be done under the same condition. In this case, several groups of step responses are selected that include forward and inverse step tests (a typical group including the set point change cycle is, for example: 495 f 497 f 495 f 493 f 495 C). However, to illustrate the modeling and controller design, just one of them is given for an example. The real-time process data when the set point is changed from 493 to 495 C can be seen in Figure 2. The temperature loop is first viewed as a rough first-order model. Generally, this is acceptable for an unknown industrial process, and also it will be easy to design the successive controller. However, if the first-order model cannot get a satisfactory accuracy, then a second-order model will be considered. Thus the linear part based on real-time data is as described in Section2 and eqs 3 and 4. When the sampling time Ts is selected, here it is chosen 5 s, the corresponding linear part is

E = ∑k=1N[e(k)]2, where N = 450, which derives the following: " #1 " # " # 0 β 0 IT ð32Þ ¼ I Ω þ γ1 I E W

yL ðkÞ ¼ a1 yðk  1Þ þ b1 uðk  50Þ

where W = [w1, w2, ..., wN]T, I = [1, 1, ..., 1]T, E = [e(1), e(2), ..., e(N)]T, Ω is a square matrix, and Ωij = j(Xi)Tj(Xj) = K(Xi,Xj), K( 3 , 3 ) is the kernel function. In this case, a radial basis function (RBF) function is adopted, K (X,Xi) = exp{||X  Xi||22/σ2}, here the parameters of LSSVM are selected as σ2 = 10 and γ = 10, and I is the unit matrix of appropriate dimension. Then yNL(k) is derived as follows: 2    N   yNL ðkÞ ¼ wi expf  X  X i  =σ2 g þ β ð33Þ   i¼1



2

where {Xi,y(i)} i=1 is the inputoutput data sample set, X is the input vector, Xi is the input vector of sampling time i, and wi, σ2, and β are the parameters of LSSVM. Thus the combination model of the process is N

ym ðkÞ ¼ a1 yðk  1Þ þ b1 uðk  50Þ 2    N   þ wi expf  X  Xi  =σ 2 g þ β   i¼1



ð31Þ

with a1 = 0.9908 and b1 = 0.0092. A further error identification between the real-time data and the above linear model e(k) = y(k)  yL(k) is through a least-squares support vector machine (LSSVM) by minimizing

ð34Þ

2

Figure 3 shows the modeling result using the above combination method. It can be seen that the modeling result is satisfactory. Note that from Figure 2, one may think that the process exhibits kind of reverse response. And both the linear model eq 31 and the nonlinear model eq 34 do not include it in them for 8116

dx.doi.org/10.1021/ie102211x |Ind. Eng. Chem. Res. 2011, 50, 8110–8121

Industrial & Engineering Chemistry Research

ARTICLE

Figure 6. (A) Output signal of proposed under model mismatch. (B) Input signal of proposed under model mismatch.

Figure 7. (A) Output signal of GPC under model mismatch. (B) Input signal of GPC under model mismatch.

their controllers’ designs, respectively. This is true because the reverse behavior also contains the process noises. After filtering the real-time response, the reverse response is not too dominant, and thus it is ignored in both the linear and the nonlinear models. 5.3. Control Strategy Design. The nonlinear predictive controller design for the derived model eq 34 is as described in Section 3. Define: NL(X) = yNL(k), then model eq 34 is rewritten as

associated with the set point. A value of μ that is larger than 0.9 but smaller than 1 is often preferred because the set point will be a very smooth one. Here it is chosen as 0.97. Experimental results are conducted under two cases: there is no model mismatch and there exists a model mismatch. The process initially operates at 491 C. From time t = 300 min, successive step changes in the set point are conducted, followed by a step disturbance of 2 t/h introduced to the process output at t = 500 min. Figures 47 show the comparison results of the two controllers. It is seen that the proposed nonlinear predictive controller gives a relatively good performance when there is no model mismatch. It manages to bring the output to the desired range as close as possible and to reject unmeasured disturbances as it can. Conventional GPC does not yield an ideal result because of the nonlinearities of the process. This is because the proposed GPC adopts an iterative strategy that manages to bring the error between the linear and nonlinear model as small as possible. And step by step, an optimal control action will be found. Conventional GPC only uses the linear model for the prediction of the process’ future behavior and calculation of control signal, thus a relatively large error is introduced to the optimization of control law which results in an unsatisfactory performance. Figures 6 and 7 show the results when there exists a model mismatch. At time t = 300 min, the model parameters are changed to a1 þ 0.008 and b1  0.008. This will cause the relatively big changes of the process gain and time constant. The control parameters for the two controllers remain the same:p = 56, λ2 = 1, and μ = 0.97. The control results in Figures 6 and 7 show that the proposed nonlinear predictive controller still yields satisfactory results. However, results of conventional GPC are not very satisfactory.

ym ðkÞ ¼ a1 yðk  1Þ þ b1 uðk  50Þ þ NLðXÞ

ð35Þ

This is the first-order case of eq 11, the nonlinear predictive controller design based on it is through eqs 1219. 5.4. Experimental Results. This section gives the experimental results of the proposed identification and nonlinear predictive control strategy compared with conventional predictive control method (GPC).48 The two controllers are both designed based on the real-time data model eq 35, they are compared under the same condition. The control parameters for the two controllers are the same: p = 56, λ2 = 1, and μ = 0.97. Though there has not been any quantitative criterion of control parameters for predictive control, some qualitative rules are already known for researchers. The basis for the choice of the above control parameters is based on these rules. That is: p is associated with the stability and closed-loop response speed of the control system. It should be chosen to ensure a stable and a fast response of closed-loop system. To guarantee this, it must incorporate the dead time of the process and not be very large. For the studied process, since the dead time is 50, p is chosen as 56. λ2 is associated with the control action. The control action is expected not to violate too much. In fact, a small value of λ2 can meet this requirement. Here λ2 is chosen as 1, and μ is actually

8117

dx.doi.org/10.1021/ie102211x |Ind. Eng. Chem. Res. 2011, 50, 8110–8121

Industrial & Engineering Chemistry Research

ARTICLE

6. CONCLUSION In this paper, a combination of a step response and nonlinear optimization based identification method and a nonlinear predictive control strategy is proposed. This combination model combines general prior knowledge and nonlinear optimization in the form of a dynamic model structure with specific information about the behavior of the process under study, which is derived from experimental process data. This reduces the modeling error. The performances of the proposed nonlinear predictive controller are verified as the optimal one and tested by real-time experiments. The proposed method provides a better response when compared with conventional GPC controller. ’ APPENDIX A: DESIGN WITH CONSTRAINTS Constraint dealing is another important issue for MPC. This section will introduce a simplified numerical linear programming solution50,51 to the constrained predictive control problem. Typically, the constraints are imposed on all future sampling instants. For unconstrained predictive control eq 16, the optimal control is U ¼ ðGT G þ λ2~IÞ1 GT ðY r  Y past  E1  E2 Þ

Figure A1. Output signal without constraints.

ðA1Þ

When the constraints are imposed on the inputs and outputs as follows: ymin e yðkÞ e ymax ,

umin e uðkÞ e umax

ðA2Þ

where ymin, ymax, umin, and umax are the constraints. Then the above eq A1 will not be the optimal one. Introduce a vector Z = [Z1, Z2, ..., Zp]T, where Z ¼ ðGT G þ λ2~IÞU  GT ðY r  Y past  E1  E2 Þ

ðA3Þ

Define a new cost function as J½U ¼ min

p

∑ jZi j i¼1

ðA4Þ

where the definition of eq A4 means that the goal is to find a control law that drives the output of the process to the set-point, while at the same time insures the steady control input.50 Define:  Zi ¼ Zþ i  Zi

ðA5Þ

 where Zþ i > 0, Zi > 0. Then eq A4 is expressed as follows:

J½U ¼ min

p

∑ jZþi  Zi j i¼1

ðA6Þ

Then, the constrained predictive control is converted to the following linear programming problem: p

min

 jZþ ∑ i  Zi j i¼1

Figure A2. Output signal with constraints.

The constraints are all expressed in terms of the manipulated variablesU. For the manipulated variables, 2 3 1 6 7 6 7 6 7 7 6 617 7 6 uðk þ 1Þ 6 7 7 6 7 ¼ 6 7uðk  1Þ 6 6 7 7 6 617 7 6l 6 7 7 6 4 5 5 4 1 uðk þ p  1Þ 2

e yM ði ¼ 1, 2 , ..., pÞ e uM ði ¼ 1, 2 , ..., pÞ

3

2

ðA7Þ

1 6 61 þ6 6l 4 1

subject to ym e ^yðk þ d þ i=kÞ um e uðk þ i  1Þ

uðkÞ

ðA8Þ

0 1

333 333

1

1

32 3 0 ΔuðkÞ 76 7 6 7 07 7 6 Δuðk þ 1Þ 7 7 6 7 05 4l 5 1 Δuðk þ p  1Þ ðA9Þ

8118

dx.doi.org/10.1021/ie102211x |Ind. Eng. Chem. Res. 2011, 50, 8110–8121

Industrial & Engineering Chemistry Research

ARTICLE

Chart 1.

Then combine eqs A2 and A9, the constraints of manipulated variables can be expressed as C1 ðumin  uðk  1ÞÞ e CU e C2 ðumax  uðk  1ÞÞ ðA10Þ where C, C1, and C2 are the corresponding matrices. For the output constraints, since E1 in eq 13 is not ^ andU is done by linearizing known, the linear relation of Y E1 through the Taylor expansion at Uk. Then the following

is derived _ ^ ¼ Y þ GU Y

ðA11Þ

where Y is the part that does not include any relation toU, and _ G is the corresponding matrix. Then combine eqs A2 and A11, the constraints of output variables can be expressed as _ C3 ymin  Y e G e C4 ymax  Y ðA12Þ 8119

dx.doi.org/10.1021/ie102211x |Ind. Eng. Chem. Res. 2011, 50, 8110–8121

Industrial & Engineering Chemistry Research where C3 and C4 are the corresponding matrices. Combine eqs 24, A3, and A12, the constraints on Z can be calculated. Then the linear programming problem eq A7 can be solved with the simplex method. Figures A1 and A2 show the constrained cases, where in Figure A1, no constraints are imposed on the process. It is seen that the output has kind of an overshoot. To test the constrained predictive control, the constraints are imposed as 490 e y(k) e 494.5 and 497 e u(k) e 500, that is, no overshoot is allowed. In Figure A2, it can be seen that the constraints are well met. And this shows that even though Talor expansion is not an accurate one, the constraints can be well treated when the linear expansion errors are small (Chart 1).

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected]. Telephone: þ86-0571-86878566.

’ ACKNOWLEDGMENT Part of this project was supported by National Natural Science Foundation of China (Grant no. 60804010 and 60974138), National Basic Research Program of China (National 973 Program: Grant no. 2009CB320602) and Key Innovation Team Program of Zhejiang Province of China (Grant no. 2009R50019). ’ REFERENCES (1) Pearson, R. K. Nonlinear empirical modeling techniques. Comput. Chem. Eng. 2006, 30, 1514. (2) Van Lith, P. F.; Betlem, B. H. L.; Roffel., B. Combining prior knowledge with data driven modeling of a batch distillation column including start-up. Comput. Chem. Eng. 2003, 27, 1021. (3) Maciej, L. Modelling and nonlinear predictive control of a yeast fermentation biochemical reactor using neural networks. Chem. Eng. J. 2008, 45, 290. (4) Dalaouti, N.; Seferlis, P. A unified modeling framework for the optimal design and dynamic simulation of staged reactive separation processes. Comput. Chem. Eng. 2006, 30, 1264. (5) Flaus, J. M.; Thevenon, L. Data flow modeling for batch and hybrid processes. ISA Trans. 2003, 42, 361. (6) Kishalay, M.; Mahesh, G. Modeling of an industrial wet grinding operation using data-driven techniques. Comput. Chem. Eng. 2006, 30, 508. (7) Oliveira, R. Combining first principles modelling and artificial neural networks: a general framework. Comput.-Aided Chem. Eng. 2004, 28, 755. (8) Van Lith, P. F.; Betlem, B. H. L.; Roffel, B. A structured modeling approach for dynamic hybrid fuzzy-first principles models. J. Process Control 2002, 12, 605. (9) Al Seyab, R. K.; Cao, Y. Differential recurrent neural network based predictive control. Comput. Chem. Eng. 2008, 32, 1533. (10) Peng, H.; Wu, J.; Inoussa, G.; Deng, Q. L.; Nakano, K. Nonlinear system modeling and predictive control using the RBF nets-based quasi-linear ARX model. Contr. Eng. Pract. 2009, 17, 59. (11) Ahmed, S.; Huang, B.; Shah, S. L. Identification from step responses with transient initial conditions. J. Process Control 2008, 18, 121. (12) Liu, M.; Wang, Q. G.; Huang, B.; Hang, C. C. Improved identification of continuous-time delay processes from piecewise step tests. J. Process Control 2007, 17, 51. (13) Zorzetto, L. F. M.; Maciel Filho, R.; Wolf-Maciel, M. R. Process modelling development through artificial neural networks and hybrid models. Comput. Chem. Eng. 2000, 24, 1355.

ARTICLE

(14) Ljung, L. System identification-theory for the user, 2nd ed; Prentice-Hall: Upper Saddle River, NJ, 1999. (15) Rivera, D. E.; Flores, M. E. Plant-friendly closed-loop identification using model predictive control, paper 227d. In Proceedings of 1999 AIChE Annual Meeting, Dallas, TX, October 31November 5, 199, AIChE: New York, 1999. (16) Zhu, Y.; Butoyi, F. Case studies on closed-loop identification for MPC. Contr. Eng. Pract. 2002, 10, 403. (17) Greco, C.; Menga, G.; Mosca, E.; Zappa, G. The performance improvements of self tuning controllers by multi-step horizons: The MUSMAR approach. Automatica 1987, 20, 681. (18) Rossiter, J. A.; Kouvaritakis, B. Modeling and implicit modeling for predictive control. Int. J. Control 2001, 74, 1085. (19) Shook, D. S.; Mohtadi, S.; Shah, S. L. A controlrelevant identification strategy for GPC. IEEE Trans. Autom. Control 1992, 37, 975. (20) Huang, B.; Wang, Z. The role of data prefiltering for integrated identification and model predictive control. In Proceedings of IFAC World Congress, July 59, 1999. (21) Agarwal, A.; Grossmann, I. E. Linear coupled component automata for MILP modeling of hybrid systems. Comput. Chem. Eng. 2009, 33, 162. (22) Georgieva, P.; Meirelesb, M.; Feyo de Azevedo, J. S. Knowledgebased hybrid modelling of a batch crystallisation when accounting for nucleation, growth and agglomeration phenomena. Chem. Eng. Sci. 2003, 58, 3699. (23) Potocnik, B.; Bemporad, A.; Torrisi, F. D.; Music, G.; Zupancic., B. Hybrid modelling and optimal control of a Multiproduct Batch Plant. Contr. Eng. Pract. 2004, 12, 1127. (24) Romijn, R.; Ozkan, Leyla.; Weiland, S.; Ludlage, J.; Marquardt, W. A grey-box modeling approach for the reduction of nonlinear systems. J. Process Control 2008, 18, 906. (25) Safavi, A. A.; Nooraii, A.; Romagnoli, J. A. A hybrid model formulation for a distillation column and the on-line optimization study. J. Process Control 1999, 9, 125. (26) Sanchez, A.; Parra, L. F. Hybrid modeling and dynamic simulation of automated batch plants. ISA Trans. 2003, 42, 401. (27) Sanchez, A.; Parra, L. F.; Baird, R.; Macchietto, S. Hybrid modeling and dynamic simulation of automated batch plants. ISA Trans. 2003, 42, 401. (28) Shaw, A. M.; Doyle, F. J.; Schwaber, J. S. A dynamic neural network approach to nonlinear process modeling. Comput. Chem. Eng. 1997, 21, 371. (29) Grosman, B.; Lewin, D. R. Automated nonlinear model predictive control using genetic programming. Comput. Chem. Eng. 2002, 26, 631. (30) Huang, B.; Malhotra, A.; Tamayo, E. C. Model predictive control relevant identification and validation. Chem. Eng. Sci. 2003, 58, 2389. (31) Kadali, R.; Huang, B.; Rossiter, A. A data driven subspace approach to predictive controller design. Contr. Eng. Pract. 2003, 11, 261. (32) Shafiee, G.; Arefi, M. M.; Jahed-Motlagh, M. R.; Jalali, A. A. Nonlinear predictive control of a polymerization reactor based on piecewise linear Wiener model. Chem. Eng. J. 2008, 143, 282. (33) Wang, X. R.; Huang, B.; Chen, T. W. Data-driven predictive control for solid oxide fuel cells. J. Process Control 2007, 17, 103. (34) Yuzgec, U.; Becerikli, Y.; Turker, M. Nonlinear predictive control of a drying process using genetic algorithms. ISA Trans. 2006, 45, 589. (35) Doyle, F. J.; Ogunnaike, B. A.; Pearson, R. K. Nonlinear modelbased control using second-order Volterra models. Automatica 1995, 31, 697. (36) Kalmukale, A. G.; Chiu, M. S.; Wang, Q. G. Partitioned modelbased IMC design using JITL modeling technique. J. Process Control 2007, 17, 757. (37) Maner, B. R.; Doyle, F. J.; Ogunnaike, B. A.; Pearson, R. K. Nonlinear model predictive control of a simulated multivariable polymerization reactor using second-order Volterra models. Automatica 1996, 32, 1285. 8120

dx.doi.org/10.1021/ie102211x |Ind. Eng. Chem. Res. 2011, 50, 8110–8121

Industrial & Engineering Chemistry Research

ARTICLE

(38) Harris, K. R.; Palazoglu, A. Studies on the analysis of nonlinear processes via functional expansions  III: controller design. Chem. Eng. Sci. 1998, 53, 4005. (39) Maksumov, A.; Mulder, D. J.; Harris, K. R.; Palazoglu, A. Experimental application of partitioned model-based control to pH neutralization. Ind. Eng. Chem. Res. 2002, 41, 744. (40) Cheng, C.; Chiu, M. S. Robust PID controller design for nonlinear processes using JITL technique. Chem. Eng. Sci. 2008, 63, 5141. (41) Al Seyab, R. K.; Cao, Y. Nonlinear system identification for predictive control using continuous time recurrent neural networks and automatic differentiation. J. Process Control 2008, 18, 568. (42) Zhang, R. D.; Wang, S. Q. Support vector machine based predictive functional control design for output temperature of coking furnace. J. Process Control 2008, 18, 439. (43) Balaji, S.; Fuxman, A.; Lakshminarayanan, S.; Forbes, J. F.; Hayes, R. E. Repetitive model predictive control of a reverse flowreactor. Chem. Eng. Sci. 2007, 62, 2154. (44) Causa, J.; Karer, G.; Nunez, A.; Saez, D.; Skrjanc, I.; Zupancic, B. Hybrid fuzzy predictive control based on genetic algorithms for the temperature control of a batch reactor. Comput. Chem. Eng. 2008, 32, 3254. (45) DeHaan, D.; Guay, M. A new real-time perspective on nonlinear model predictive control. J. Process Control 2006, 16, 615. (46) Karer, G.; Music, G.; Skrjanc, I.; Zupancic, B. Hybrid fuzzy model-based predictive control of temperature in a batch reactor. Comput. Chem. Eng. 2007, 31, 1552. (47) Wang, S. Q. Control engineering in industrial processes; Chemical Engineering Press: Beijing, China, 2003. (48) Clarke, D. W.; Mohtadi, C.; Tuffs, P. S. Generalized predictive control—part I: the basic algorithm. Automatica 1987, 23, 137. (49) Fang, C. Z.; Xiao, D. Y. System identification; Tsinghua University Press: Beijing, China, 1998. (50) Zou, Z. Y.; Liu, J. Y.; Yu, D. H.; Liu, X. H.; Zhao, D. D.; Wang, Z. Z.; Guo, N. A novel predictive control algorithm of linearly constrained system. J. Chem. Ind. Eng. (China) 2010, 61, 413. (51) Luenberger, D. G.; Ye, Y. Y. Linear and Nonlinear Programming; Springer: New York, 2008.

8121

dx.doi.org/10.1021/ie102211x |Ind. Eng. Chem. Res. 2011, 50, 8110–8121