Multilinear Model Control of Hammerstein-like Systems Based on an

Mar 23, 2009 - angle dividing method is proposed to decompose the operating spaces and determine linear model banks for. SISO Hammerstein-like systems...
0 downloads 4 Views 272KB Size
Download PDF
3934

Ind. Eng. Chem. Res. 2009, 48, 3934–3943

Multilinear Model Control of Hammerstein-like Systems Based on an Included Angle Dividing Method and the MLD-MPC Strategy Jingjing Du, Chunyue Song,* and Ping Li State Key Laboratory of Industrial Control Technology, Institute of Industrial Process Control, Zhejiang UniVersity, Hangzhou 310027, China

This article introduces the multilinear model control approach to SISO Hammerstein-like systems. An included angle dividing method is proposed to decompose the operating spaces and determine linear model banks for SISO Hammerstein-like systems through measuring the nonlinearity of the static IO curves. On the basis of the linear model bank, the MLD-MPC strategy is employed for set-point tracking control. Two chemical processes are studied using the proposed included angle dividing method and MLD-MPC strategy. Closedloop simulations demonstrate that the proposed dividing method is useful and effective for the multilinear model approach and the MLD-MPC strategy is excellent for nonlinear systems. 1. Introduction Hammerstein models and Wiener models are popular nonlinear empirical modeling structures.1,2 As is depicted in Figure 1, a Hammerstein model is a block-oriented nonlinear model, which consists of the cascade structure of a static nonlinear function f( · ) followed by a linear dynamic block H(z),3 whereas a Wiener model contains the same elements in the reverse order.3,4 It has been shown that the Hammerstein and Wiener model structures can effectively represent and approximate many industrial processes, as they may account for nonlinear effects encountered in most chemical processes.5 For example, pH neutralization processes, distillation columns, heat exchangers, polymerization reactors, and dryer processes have been modeled by Hammerstein and Wiener models.6,7 When the Hammerstein model or the Wiener model is used for control purposes, a popular way is to remove the static nonlinearity from the original control problem via the inversion of the nonlinear element first and convert the nonlinear control problem into a new linear one, and then solve the linear control problem, and finally return the linear conversion to a nonlinear one by nonlinear inversion.4,8-10 However, to obtain a unique solution with the nonlinearity inversion method, the nonlinearity of the model needs to be bijective, which is generally not the case.9 And because nonlinearity is always incorporated in the original control problem, for it usually involves both the physical output and the physical input, separating the nonlinearity from the control problem therefore may lead to a decreased performance.10 The multilinear model control approach has attracted much attention and has been studied extensively in the past years11-18 for its potential ability to control chemical processes with highly nonlinear characteristics as a result of wide operating ranges and large set-point changes. The key concept is to represent the nonlinear system as a combination of linear systems to which classical controller design techniques can be easily applied.17 The multilinear model control approach seems attractive for nonlinear systems with a Hammerstein or Wiener model structure. However, the question of how many and which linear models are required to span the expected operating range of a nonlinear system remains largely unanswered,11,17,18 though it must be addressed in the first place when designing a multilinear * To whom correspondence should be addressed. Tel.: +86 571 87951442. E-mail: [email protected].

model controller.11 Gala´n et al.17 use gap metric to analyze the relationships among candidate models to get a reduced model bank that provides enough information for multilinear controller design. However, how to get those candidate models is not addressed explicitly. Tan et al.18 extend the method in Gala´n et al.17 to integrate the procedure of selecting operating points and local controller design. Whereas, their method needs extra compensators, and moreover, the combined global model might not cover the entire operating space, which gives rise to poor control performances and even instability when the multilinear model approach is used.11 In this article, an included angle dividing method is proposed, aiming to solve the problem of how many and which models are required to span an expected operating range primarily for SISO nonlinear systems with Hammerstein or Wiener model characteristics. For a given threshold value, the proposed included angle method divides the operating range and obtains a linear model bank without a prior set of linear models (which is needed in ref 17). This method is intuitive and easy to carry out, and, moreover, it works effectively. On the basis of the linear model bank, a global multilinear model-based controller is designed. However, general hard switching and soft switching multilinear model-based approaches may have difficulty in scheduling linear subsystems in coordination with one another and thus may incur oscillation when the system is transiting between different subsystems.19,20 Because any control method can be used within the multilinear model framework, as is pointed out in ref 11, in this article the MLD-MPC strategy is employed for the following reasons. First, multilinear models can be turned into an MLD form by HYSDEL directly and conveniently.21-23 Second, MLD-MPC can unite multiple linear models in a unified framework, optimize the control variable

Figure 1. a. Structure of a Hammerstein model. b. Structure of a Wiener model.

10.1021/ie8009395 CCC: $40.75  2009 American Chemical Society Published on Web 03/23/2009

Ind. Eng. Chem. Res., Vol. 48, No. 8, 2009 3935

systematically, avoid oscillation during switching, and augment the stability robustness of the whole system.19,20 Third, the toolbox MPT solves the MLD-MPC conveniently.24 The remaining part of the article is organized as follows: In section 2, Hammerstein-like systems are briefly discussed. Section 3 discusses the included angle method in detail. Section 4 is about the MLD-MPC strategy, and this strategy is applied to two examples later in section 5 where closed-loop simulation results are reported. In section 6, some conclusions and discussions are presented. 2. Hammerstein-like Systems Most chemical processes exhibit strong nonlinearity as a result of wide operating ranges and set-point changes. Generally speaking, these processes have both static nonlinearity and dynamic nonlinearity, whereas some may have primary static nonlinearity and some may have primary dynamical nonlinearity. For example, consider a system described by (dy)/(dt) + {y3 - u3}1/3 ) 0. It is obvious that this system is nonlinear as it is represented by a nonlinear differential equation. It can be further found that its static behavior is linear. So, the nonlinearity lies totally in its dynamical behavior. However, for a Hammerstein model described by eq 1 and a Wiener model eq 2, as is well-known that their nonlinear characteristics lie in the static gains not in the dynamics.25 That is to say, their dynamics remain unchanged, whereas their static gains vary with operating points. As is known to all, the static gain of a SISO system (no matter linear or nonlinear) at any operating point is the slope of the steady-state input-output curve.26 Then for a SISO Hammerstein model or a SISO Wiener model, the slope of its static IO curve varies with its operating point and reflects all of its nonlinearity. That is a pretty good feature.

{ {

x(k) ) f (u(k)) n

y(k) )

w(k) )

∑

m

ai y(k - i) +

i)1

(1)

i

i)1

n

∑

∑ b x(k - i) m

aiw(k - i) +

i)1

∑ b u(k - i) i

i)1

(2)

y(k) ) f (w(k))

where x and w are intermediate variables; and f ( · ) is analytic and not necessarily invertible in this article. In this article, a physical system that can be represented or approximated by a Hammerstein or Wiener model structure is defined as a Hammerstein-like system for notation simplicity. Namely, a Hammerstein-like system denotes a nonlinear process whose nonlinearity lies primarily in its static behavior, whereas dynamics are basically linear. Then most of the industrial processes where Wiener and Hammerstein model structures are employed, like pH neutralization processes, distillation columns, heat exchangers, polymerization reactors, and dryer processes,6,7 belong to Hammerstein-like systems. A Hammerstein model or a Wiener model denotes a system that can be modeled exactly (note: not approximately) by eq 1 or eq 2, which is an ideal mathematical model and its dynamic behavior is linear. Apparently, a Hammerstein-like system is not equivalent to a Hammerstein model or a Wiener model in this article. In brief, the Hammerstein model and the Wiener model are included in the Hammerstein-like system. From the previous analyses, we can get the conclusion that a Hammerstein-like nonlinear system exhibits primarily static

nonlinear behavior, that is its static gain varies as the operating point changes, whereas its dynamic parameters can be supposed to keep constant and can use the average values of the entire operating space26 or the approximate values instead. Then it is seemingly reasonable to measure the nonlinearity of a Hammerstein-like system through its static behavior. As has been pointed out in ref 13 that using physical information from the process static input-output map, the controller designer is able to select a piecewise approximation, which contains most of the information available, for the nonlinear system. Then if the multilinear model approach is applied to a Hammersteinlike system to design a controller, we can try to decompose the entire operating range according to the nonlinearity of its steadystate IO curve, and determine the number of submodels to approximate the nonlinear system. In this work, for a SISO Hammerstein-like system, an included angle method is proposed to decompose its operating range and then the multilinear model approach is employed to design a global controller. Thus, invertibility of the nonlinear element is not necessary and the design goal will not be influenced by the nonlinearity directly. Note that in this work only SISO systems are discussed and we focus on Hammerstein-like systems that contain Hammerstein models and Wiener models. 3. Included Angle Dividing Method in Operating Space Division For a Hammerstein-like system, the identification of the linear model and the static nonlinear function can be carried out independently.6,7,27-29 That is to say, the static and dynamic models of a Hammerstein-like system can be obtained separately. In particular, the steady-state process characterizations are often more readily available than dynamic ones, partly because steady-state models are simpler and easier to develop than detailed dynamic models.27,30 In this article, we always assume that the steady-state model of a nonlinear system is already aVailable (either from identification or from first principle modeling) and continuous analytic. Then for a nonlinear SISO Hammerstein-like system, we assess its nonlinearity through its static IO map and decompose its operating range into subregions via the following included angle method. Then for every subregion a local linear model can either be identified (for a black-box system and a gray-box system) or linearized around a proper operating point (for a white-box system). Finally, the local linear models are combined into a global model, that is a MLD model, based on which the MPC strategy is employed for set-point tracking control. 3.1. Included Angle Concept. As was pointed out previously, the static gain of a SISO nonlinear system equals the slope of its steady-state IO curve at any operating point and varies as its operating point changes. So, for a Hammersteinlike system with primarily static nonlinearity, the slope of its steady-state IO map can be used to measure its nonlinearity. However, the variation range of the slope is from negatively infinite to positively infinite, which is not convenient for computation or comparison. So, it is better to change the variable range into a finite one. As is known, the slope is a function of its slope angle, and the variation range of an angle is finite. So, the slope angle can be used instead of slope to describe the variation of a system’s static behavior. The slope angle at the ith operating point is computed via: θi ) arctan ki,i ) 1, 2,..., Na where ki is the slope at the ith operating point.

(3)

3936 Ind. Eng. Chem. Res., Vol. 48, No. 8, 2009

The difference between two slope angles at two operating points is defined as the absolute value of the difference: θij ) |θi - θj |,

i, j ) 1, 2,..., Na

(4)

We define θij as the included angle of the two operating points i, j. Its variation range is [0, π], which is finite. Here, we use the included angle of two operating points to measure the difference of the system’s behavior at the two operating points. When the included angle between two operating points is 0, the nonlinear system has similar static gains and similar system behaviors at the two points; when the included angle is bigger than π/2, the static gains of the nonlinear system are far apart and the system’s behaviors are rather different at the two points. 3.2. Procedure of the Included Angle Dividing Method. The procedure of the included angle dividing method consists of the following steps. Step 0. Prescribe a threshold value γa and the number of steady-state points Na according to the previous knowledge and performance requirement about the system. Step 1. Distribute the Na steady-state points along the steadystate IO curve evenly. Mark the ith steady-state point as Pi, i ) 1, 2,..., Na. Step 2. Compute the slope value ki at the ith steady-state point Pi and calculate its slope angle according to eq 3. Step 3. Then the included angles between any two of the Na steady-state points are computed via eq 4; and an Na × Na anglematrix [θij]Na × Na is computed. Step 4. Compare θij (i g j) with γa one by one in sequence from i ) j ) 1. If θij g γa, stop. Suppose nm (m means the times Step 4 has been executed) successive steady-state points are concerned during the comparison, and then these points are to be classified in one subregion and modeled by one submodel. The one in the middle of these nm points is chosen as the operating point, and, later on, a submodel can be identified or linearized with respect to this point. The range covered by these nm points is the submodel’s valid subregion. Step 5. Set i ) j ) i + 1, and repeat Step 4 until all the Na steady-state points are classified. Step 6. If some subregion, for example the mth subregion, acquired by the previous steps covers a rather small range with respect to the entire range, this subregion will be merged into its neighbor subregion(s) (m - 1th and/or m + 1th) because it is trivial to build a submodel for a rather small range. This subregion can be merged into its left (m - 1th) neighbor, right (m + 1th) neighbor, or both averagely according to its relationships to the neighbors in the included angle sense. Then the submodel of the neighbor is updated accordingly. Then the decomposition of the operating space is finished. Step 7. In each subregion, either identify the local linear model with respect to its operating point or linearize the nonlinear first-principle model around its operating point; and finally a linear model bank is set up. Remark 1. In this dividing procedure, the number of steadystate points Na and the threshold value γa are two designerassigned parameters, which should be chosen carefully according to a priori knowledge and performance requirements of the system. Generally speaking, to cover the system characteristic as much as possible, Na should be a large number, maybe 50, 100, and so forth, whereas a larger Na involves more computation complexity than a smaller one; and a smaller γa may result in more linear models, which may complicate the following controller design, whereas a bigger γa may produce less linear models that may not be enough to approximate the nonlinear

system. So, the choice of Na and γa is a tradeoff issue and should be addressed cautiously. Remark 2. The angle-matrix [θij]Na × Na is symmetric with respect to the main diagonal, so in Step 4 comparing the elements either below or above the main diagonal with γa is enough. Remark 3. A steady-state point denotes a point (x0, u0, y0) that satisfies the steady-state equation of the system. An operating point denotes the steady-state point around which a linear model is built (see Step 4). Remark 4. In this work, the proposed method can be applied to both black-box systems and white-box systems; and the controller design is based on a state-space model structure, so if IO models are identified at first, then they will be converted into state-space models immediately. Remark 5. Apparently, the dividing results have close relation to the boundary of the operating space, especially the steady-state point that starts the procedure. Different starting points usually have different dividing results. From the dividing procedure above, we can see that the efficiency of the proposed method depends on the calculation of the included angle heavily as the included angle is a function of the static gain. So, we would like to analyze the sensitivity of the proposed method to the static gain. The relation between the slope angle θ and the static gain k is described by eq 3. The derivative of θ with respect to k is (dθ)/(dk) ) (1)/(1 + k2). As we know 0 e |k| < ∞, so when |k| is small, θ is very sensitive to it, whereas when k f ∞, θ become insensitive. So, if the slope of the static IO curve is small, we should choose a larger Na to distribute more steady-state points in the operating space to capture the nonlinear behavior of the system; and if the IO curve slope is big, a small Na may be enough. In the following examples, the effectiveness of this method is illustrated by a SISO CSTR system and a SISO pH process. 3.3. Examples. In this part, two Hammerstein-like systems are presented to demonstrate the proposed included angle dividing method, a CSTR process, and a pH process. The former is a benchmark Hammerstein-type system and the latter is a typical Wiener-type system. For the CSTR system, we divide the operating space into subspaces and identify the linear local models for each subregion, whereas for the pH process we assume the first-principle model is given, so the linear models are obtained by linearizing the first-principle model around the chosen operating points. 3.3.1. CSTR Process. Consider a continuous stirred-tank reactor (CSTR) process, which has been modeled into a Hammerstein structure by many researchers.6,31,32 This CSTR system consists of an irreversible, exothermic reaction, A f B, in a constant volume reactor cooled by a singled coolant stream. It can be modeled by the following equations:

{

E q C˙A ) [CA0 - CA(t)] - k0CA(t)e- RT(t) V E q T˙(t) ) [T0 - T(t)] + k1CA(t)e- RT(t) + V

(5)

k3

k2qc(t)[1 - e- qc(t) ][Tc0 - T(t)] y(t) ) CA(t)

The output and input variables are CA and qc, respectively. The objective is to control the measured concentration of A, CA, by manipulating the coolant flow rate qc for set-point tracking control. The nominal parameter values appear in Table 1.

Ind. Eng. Chem. Res., Vol. 48, No. 8, 2009 3937 Table 1. Nominal CSTR Parameter Values measured product concentration coolant flow rate feed concentration inlet coolant temperature heat transfer term activation energy term liquid densities k1 ) -(∆Hk0)/(FCp)

CA qc CA0 TC0 hA E/R F,Fc

0.1 mol/L reactor temperature 103.41 Lmin-1 process flow rate 1mol/L feed temperature 350K CSTR volume 5 7 × 10 cal/min K reaction rate constant 1 × 104 K heat of reaction 1 × 103 g/L specific heats k2 )(FcCpc)/(FCpV)

T q T0 V k0 ∆H Cp,Cpc

438.51 K 100 L min-1 350 K 100 L 7.2 × 1010 min-1 -2 × 105 cal/mol 1 cal g-1 K-1 k3 ) (hA)/(FcCpc)

Table 2. Dividing Results of Case 1 for the CSTR Process subregion

1st

2nd

3rd

steady-state point-included operating point (CA, qc) subrange

1 f 128 65th (0.038122, 76) {y|0.02 e y e 0.06593}

129 f 176 153rd (0.082119, 98) {y|0.06593 < y e 0.10035}

177 f 200 189th (0.11481, 107) {y|0.10035 < y e 0.13}

Table 3. Dividing Results of Case 2 for the CSTR Process subregion

1st

2nd

3rd

steady-state point-included operating point (CA, qc) subrange

1 f 64 33rd (0.038122, 76) {y|0.02 e y e 0.06536}

65 f 88 77th (0.082119, 98) {y|0.06536 < y e 0.10035}

89 f 100 95th (0.11481, 107) {y|0.10035 < y e 0.13}

Table 4. Dividing Results of Case 3 for the CSTR Process subregion

1st

2nd

3rd

steady-state point-included operating point (CA, qc) subrange

1 f 32 17th (0.038122, 76) {y|0.02 e y e 0.06423}

33 f 44 39th (0.082119, 98) {y|0.06423 < y e 0.09849}

45 f 50 48th (0.11481, 107) {y|0.9849 < y e 0.13}

The steady-state input-output map of the CSTR process is shown in Figure 2, from which we can see the CSTR system exhibits strong static nonlinearity. A single linear controller is inadequate and a multilinear controller may be satisfactory. Apply the proposed included angle dividing method to the CSTR and decompose its operating space into subregions. We will give four dividing results with different designer-assigned parameters. Case 1. Na ) 200, γa ) 0.0015. The dividing results are summarized in Table 2. The entire operating space is decomposed into three subregions marked first, second, and third. The second column of Table 2 means first subregion includes first to 128th steady-state points, the operating point of the first subregion is the 65th steady-state point, and the subrange that first subregion covers is {y|0.02 e y e 0.06593}. The other columns can be similarly interpreted, and the following tables are interpreted like Table 2.

Figure 2. Steady-state input-output map of CSTR.

Table 5. Dividing Results of Case 4 for the CSTR Process subregion steady-state pointincluded operating point (CA, qc) subrange

1st

2nd

1 f 165

166 f 200

83rd (0.044632, 80.5)

183rd (0.10827, 105.5)

{y|0.02 e y e 0.0915}

{y|0.0915 < y e 0.13}

Case 2. Na ) 100, γa ) 0.0015. In this case, the operating space is also divided into threes subregions, and the results are listed in Table 3. Case 3. Na ) 50, γa ) 0.0015. Also, three subregions are produced. The results are summarized in Table 4. In the three cases above, Na is chosen to be 200, 100, and 50, respectively; and γa is chosen to be 0.0015. The whole operating range is decomposed into three subregions in all of these cases. Compare the three tables and it is clearly seen that the three cases above almost have the same decomposition: the same number of subregions, nearly the same subranges, and even the same operating points, respectively! Generally speaking, the bigger Na is, the more information the dividing procedure carries, and the more precise the dividing is, though a bigger Na also means a heavier computational burden. The unusual case of this CSTR process is supposed to be closely related to the monotonous and uniform rise of its steady-state curve. Although for this CSTR system the change of Na does not make big differences in the dividing results, it does make some differences. Look at the first subranges in the three cases: the bigger Na is, the wider the subrange is, because a bigger Na makes the dividing finer and more accurate. The other two subregions cannot compare in this way as they may not have the same initial point to start the dividing procedure. Case 4. Na ) 200, γa ) 0.0025. In this case, the operating space is decomposed into two subregions as is summarized in Table 5. Comparing Case 4 with Case 1, they have the same Na and different γa′s, whereas then they have totally different decompositions. So, it is important to choose a proper threshold value

3938 Ind. Eng. Chem. Res., Vol. 48, No. 8, 2009 Table 6. Linear Model Bank for Case 1 subregion

1st

2nd

3rd

transfer g1 ) 0.001 × g2 ) 0.001 × g3 ) 0.001 × function (0.13809z- 1 + (0.18165z- 1 + (0.18712z- 1 + 0.06849z- 2)/ 0.15272z- 2)/ 0.17708z- 2)/ (1 - 1.0429z- 1 + (1 - 1.5136z- 1 + (1 - 1.7857z- 1 + 0.19872z- 2) 0.62735z- 2) 0.86422z- 2) Table 7. Linear Model Bank for Case 4 subregion

1st

2nd

transfer function

g1 ) 0.001 × (0.15734z- 1 + 0.097429z- 2)/ (1 - 1.2384z- 1 + 0.3751z- 2)

g2 ) 0.001 × (0.18747z- 1 + 0.17369z- 2)/ (1 - 1.7755z- 1 + 0.8518z- 2)

according to the previous knowledge and control requirement to have a proper division. We identify the linear model banks in the second-order ARX form for Case 1 and Case 4, which are listed in Table 6 and 7, respectively. These models will be converted into state-space models, based on which MLD-MPC controllers will be designed for set-point tracking in section 5. Parameters of the state-space models (Ai, Bi, Ci, Di, i ) 1, 2, 3) are omitted for brevity. Simulation results demonstrate that Case 1 is a reasonable choice. 3.3.2. pH Process. The pH process is a typical Wiener type process where the material balance differential equations are almost linear and the equilibrium equation (titration curve) is a strong nonlinear static function.6,28 Consider the pH process in Figure 3, where neutralization reaction take place among the strong acid HNO3, the buffer agent NaHCO3, and the strong base NaOH (for details, refer to ref 6). The dynamic behavior of this reaction is described by the following differential equations.

{

1 h˙ ) (q1 + q2 + q3 - cvh0.5) A ˙ a ) 1 [(Wa - Wa )q1 + (Wa - Wa )q2 + (Wa - Wa )q3] W 4 1 4 2 4 3 4 Ah 1 ˙ Wb4 ) [(Wb1 - Wb4)q1 + (Wb2 - Wb4)q2 + (Wb3 - Wb4)q3] Ah 1 + 2 × 10pH4-pK2 Wa4 + 10pH4-14 + Wb4 × - 10-pH4 ) 0 1 + 10pK1-pH4 + 10pH4-pK2 (6)

where q1, Wa1, and Wb1 are the flow, the charge balance factor, and the material balance factor of the strong acid HNO3, respectively; and q2, Wa2, and Wb2 are those for the buffer agent NaHCO3, respectively; q3, Wa3, and Wb3 are those for the strong base NaOH, respectively; q4, Wa4 and, Wb4 are those for the effluent stream, respectively; and h is the level of the reactor; pH4 is the measured pH value of the effluent stream, also the output variable of the process. The control objective is to drive the system to different pH conditions (tracking control) by manipulating the input variable q3, the flow of NaOH. The nominal parameters are listed in Table 8. The steady state input-output map is displayed in Figure 4, from which we see that this pH process exhibits strong output nonlinearity, and five distinct operating subregions are visually

Figure 3. pH neutralization process.

discerned. We will apply the proposed included method to it for three different cases in the sequel. Case 1. Na ) 100, γa ) π/5. The entire operating space is divided into five subregions just as intuition. The decomposition is listed in Table 9. Case 2. Na ) 50, γa ) π/5. Four subregions are produced through the included angle dividing, which is summarized in Table 10. Compared to Case 1, Case 2 has the same γa, whereas a much smaller Na results in a completely different dividing result. That is because a smaller Na means less information about the system behavior, which makes the dividing rougher and poorer. Case 3. Na ) 100, γa ) π/4. As is listed in Table 11, the operating space is decomposed into four subregions, which does not agree with intuition. Compared with Case 1, Case 3 has the same Na and a different threshold value γa, which produce different numbers of subregions. An appropriate threshold value is indispensable for a reasonable division. We get the linear model banks for case 1 and case 3 by linearizing the nonlinear system around the operating points in the corresponding subregions. Parameters of the state-space models (Ai, Bi, Ci, Di, i ) 1, 2, 3, 4, 5) are omitted for brevity. MLD-MPC controllers will be designed based on the local linear models in Section 5, and closed-loop simulations approve of Case 1. 4. MLD-MPC for Nonlinear Systems 4.1. Mixed Logical Dynamical Model. The general form of a mixed logical dynamical (MLD) hybrid system is:33,34

{

x(k + 1) ) Ax(k) + B1u(k) + B2δ(k) + B3z(k) y(k) ) Cx(k) + D1u(k) + D2δ(k) + D3z(k) E2δ(k) + E3z(k) e E1u(k) + E4x(k) + E5

where, k ∈ Z, x ∈ Rnc × {0,1}nl denotes the states of the system, u ∈ Rmc × {0,1}ml denotes the inputs, and y ∈ Rpc × {0,1}pl denotes the outputs, with both real and binary components. Furthermore, δ ∈ {0,1}γl and z ∈ Rγc represent the binary and auxiliary continuous variables, respectively. The detail theory is referred to ref 33. 4.2. Convert a Hammerstein-like Nonlinear System into the MLD Formulation. Consider a continuous Hammerstein-like system described by a nonlinear state-space model in the following form:

Table 8. Nominal Parameter Values of the pH Process A ) 207 cm2 h ) 14.0 cm pK1 ) 6.35 pK2 ) 10.25

q1 ) 16.6 mL/s q2 ) 0.55 mL/s q3 ) 15.6 mL/s pH4 ) 7

(7)

Wa1 ) 3 × 10- 3 mol/L -2

Wa2 ) -3 × 10 mol/L Wa3 ) -3.06 × 10-3 mol/L cv ) 8.75 mL/(cm · s)

Wb1 ) 0 Wb2 ) 3 × 10-2 mol/L Wb3 ) 5 × 10-5 mol/L

Ind. Eng. Chem. Res., Vol. 48, No. 8, 2009 3939

{

x˙ ) f (x, u) y ) g(x, u)

(8)

where, x ∈ Rn is the state vector, u ∈ Rr is the control input vector, and y ∈ Rm is the output vector. f( · ) and g( · ) are nonlinear functions. Suppose the operating space of system (8) is decomposed into Nm subregions by the proposed included angle dividing method, and each subregion has an operating point in it. Then linearize eq 8 around the ith operating point (xoi, uoi, yoi) in the ith subregion, and we get the following linear state space equations.

{

δx˙ ) Aˆiδx + Bˆiδu ˆ iδu δy ) Cˆiδx + D

i ) 1,..., Nm

(9)

where δx ) x - xoi, δu ) u - uoi, δy ) y - yoi, Aˆi ) (∂f(x,u))/ (∂x)|(xoi,uoi), Bˆi ) (∂f(x,u))/(∂u)|(xoi,uoi), Cˆi ) (∂g(x,u))/ ˆ i ) (∂g(x,u))/(∂u)|(xoi,uoi). (∂x)|(xoi,uoi), and D After discretization with sampling period Ts, eq 9 is transformed as:

{

∆x(k + 1) ) Ai∆x(k) + Bi∆u(k) ∆y(k) ) Ci∆x(k) + Di∆u(k)

Figure 4. Titration curve (static input-output map) of the pH process.

i ) 1,..., Nm (10)

where, k ∈ Z, ∆x(k + 1) ) x(k + 1) - xoi, ∆x(k) ) x(k) - xoi, ∆y(k) ) y(k) - yoi, ∆u(k) ) u(k) - uoi, and Ai, Bi, Ci, and Di ˆ i. are the discretized matrices of Aˆi, Bˆi, Cˆi, and D Eq 10 can be rewritten as the following piecewise affine (PWA) form: x(k + 1) ) Aix(k) + Biu(k) + bi i ) 1,..., Nm (11) y(k) ) Cix(k) + Diu(k) + di

{

where bi ) xoi - Aixoi - Biuoi, and di ) yoi - Cixoi - Diuoi. Introduce logic variables δi(k), i ) 1,..., Nm, and define them as: δi(k) ) 1 T

{

x(k + 1) ) Aix(k) + Biu(k) + bi y(k) ) Cix(k) + Diu(k) + di

Figure 5. Closed-loop response and control input of CSTR using MLDMPC controller based on two linear models.

(12)

Then we have Nm

∑ δ (k) ) 1

(13)

i

i)1

{

The global model of the nonlinear system is combined as: Nm

x(k + 1) )

∑ [A x(k) + B u(k) + b ]δ (k) i

i

i

i

i)1

Nm

y(k) )

∑ [C x(k) + D u(k) + d ]δ (k) i

i

i

(14)

i

i)1

Introduce auxiliary vectors zxi(k), zyi(k), which have the same dimension of x(k). Define them as zxi(k) ) [Aix(k) + Biu(k) + bi]δi(k) zyi(k) ) [Cix(k) + Diu(k) + di]δi(k)

(15)

The following operating constraints may also exist: umin e u e umax, xmin e x e xmax

(16)

Transform eqs 12, 13, 15, and 16 into inequalities according to the rules in ref 33. Let δ(k) )[δ1′(k), · · · ,δ′Nm(k)]′, z(k) )

Figure 6. Closed-loop response and control input of CSTR using MLDMPC controller based on three linear models.

[z′x1(k), · · · , z′xNm(k),z′y1(k), · · · , z′yNm(k)]′. On the basis of these inequalities, eq 14 can be rewritten as:

{

x(k + 1) ) Amx(k) + Bm1u(k) + Bm2δ(k) + Bm3z(k) y(k) ) Cmx(k) + Dm1u(k) + Dm2δ(k) + Dm3z(k) Em2δ(k) + Em3z(k) e Em1u(k) + Em4x(k) + Em5 (17)

where Am, Bm1, Bm2, Bm3, Cm, Dm1, Dm2, Dm3, and Em1,..., Em5 can be obtained from the transformation procedure accordingly.

3940 Ind. Eng. Chem. Res., Vol. 48, No. 8, 2009 N-1

min J(ukN-1, x(k)) :)

{ukN-1}

∑ (|u(k + i) - u |

e Q1

+ |δ(i|k) - δe | Q2+

i)0

|z(i|k) - ze | Q3 + |x(i|k) - xe | Q4 + |y(i|k) - ye | Q3)

(18)

s. t.

{

Figure 7. Closed-loop response and control input of CSTR using MultiPID controller based on three linear models.

x(N|k) ) xe x(i + 1|k) ) Amx(i|k) + Bm1u(k + i) + Bm2δ(i|k) + Bm3z(i|k) y(i|k) ) Cmx(i|k) + Dm1u(k + i) + Dm2δ(i|k) + Dm3z(i|k) Em2δ(i|k) + Em3z(i|k) e Em1u(k + i) + Em4x(i|k) + Em5 (19)

where, k is the current time step, N is the prediction horizon, and x(k) is the current state. x(i|k) is the state predicted at time step k + i according to x(k) and input sequence uNk - 1: ){u(k), u(k + 1),..., u(k + N - 1)}, and δ(i|k), z(i|k), y(i|k) are similarly defined, i ) 1,..., N, Qj ) Q′j g 0, j ) 1,..., 5 are the weighted matrix, and xe, ue, ye are the values of a steady point that satisfies the nonlinear eq 8 and δe, ze are the corresponding auxiliary variables. Assume for the moment k that the optimal solution {u*k (i)}i ) 0,..., N-1 exists. According to the receding horizon philosophy, we get the control input, u(k) ) u*k (0)

(20)

At time step k + 1, the whole optimization procedure is repeated, and this online replanning provides the desired feedback control action on the premise of guaranteeing the stability of system (17). See the following lemma. Figure 8. Closed-loop response and control input of CSTR using nonlinear PID controller.

Figure 10. pH set-point tracking performance of MLD-MPC controller with five models. Figure 9. pH set-point tracking performance of MLD-MPC controller with four models.

Other physical constraints and logical rules can also be incorporated into the inequality of eq 17, and eq 17 is the MLD formulation of nonlinear system (8). The above procedure of obtaining an MLD form of eq 17 can be automatized by the HYSDEL software package effectively and conveniently.21-23 Remark 6. For systems where no first principle models is available, a local linear model can be identified in each subregion after operating space decomposition, and then turned into the PWA formulation (11) and finally converted into eq 17. 4.3. MPC Based on MLD Model. The MPC problem based on an MLD model is formulated as:33

Figure 11. pH set-point tracking performance of Multi-PID controller with five models.

Ind. Eng. Chem. Res., Vol. 48, No. 8, 2009 3941 Table 9. Dividing Results of Case 1 for the pH Process subregion

1st

2nd

3rd

4th

5th

steady-state point included operating point (q3, pH4) subrange

1 f 34 18th (5.4, 3.128) {y|2 e y e 4.10}

35 f 42 39th (11.7,5.60) {y|4.10 < y e 6.00}

43 f 54 49th (14.7,6.67) {y|6.00 < y e 7.52}

55 f 57 56th (16.8, 8.82) {y|7.52 < y e 9.13}

58 f 100 79th (23.7, 10.43) {y|9.13 < y e 11}

Table 10. Dividing Results of Case 2 for the pH Process subregion

1st

2nd

3rd

4th

steady-state point-included operating point (q3, pH4) subrange

1 f 17 9th (5.1, 3.10) {y|2 e y e 3.94}

18 f 27 23rd (13.5,6.29) {y|3.94 < y e 7.24}

28 f 29 29th (17.1,9.13) {y|7.24 < y e 9.13}

30 f 50 40th (23.7, 10.43) {y|9.13 < y e 11}

Lemma.33 Let (xe, ue) be an equilibrium pair of eq 8 and (δe, ze) definitely admissible. Assume that the initial state x(0) is such that a feasible solution of problem (18)-(19) exists at time t ) 0. Then ∀Q1 ) Q1′ > 0, Qi ) Q′i g 0 (i ) 2, 3, 4, 5) and the mixed integer programming control (MIPC) law (18)-(20) stabilizes the system (17) in that limk f∞x(k) ) xe, limk f∞u(k) ) ue, limk f∞|δ(k) - δe|Q2 ) 0, limk f∞|z(k) - ze|Q3 ) 0, limk f∞|y(k) - ye|Q5 ) 0, while fulfilling the dynamic/ relational constraints in eq 17. Detail about the stability proof is referred to ref 33. Note that the stability result is valid for nonlinear system (8) provided that the system is exactly described by a set of linear models eqs 9 or 10. However, in this MLD-MPC strategy, the linear models are only used to approximate the nonlinear system behavior, there always exist some modeling errors, and it is rather probable that the resulting response can be quite different from what has been predicted. Therefore, stability of the hybrid system (17) does not necessarily guarantee stability of the nonlinear system (8), a limitation common to many multilinear model control approaches applied to nonlinear systems.16 However, if a proper linear model bank is selected to approximate the nonlinear system, the MLD-MPC strategy provides an acceptable control for nonlinear systems while guaranteeing robust stability of the closed-loop system as is done in refs 11 and 16. Further study of robustness is beyond the scope of this work. 5. Simulations In this section, the MLD-MPC strategy presented in section 4 is applied to the CSTR and pH processes for set-point tracking control based on the linear model banks determined in section 3. For comparison, Multi-PID controllers35 based on the same linear model banks and nonlinear PID controllers based on nonlinearity inversion28 are designed. In this work, the MultiPID controllers are designed strictly according to the method in ref 35, where the Gaussian-weighted function is used to combine the local linear PID controllers into a global one; and the nonlinear PID technique is totally from ref 28, where the design procedure can be summarized as the following two steps: First, consider the system as a connection of a static nonlinear function and a linear dynamic element, and design a linear PID controller for the linear dynamic element. Second, convert the linear PID controller, via the inversion of the static nonlinear function, into a nonlinear PID controller, which is the final controller that acts on the controlled nonlinear system. 5.1. CSTR Process. In Figure 5, the MLD-MPC controller is designed based on the model bank from Table 7. It is seen that the tracking performance is poor, the input jumps frequently, and the output oscillates fiercely in most areas and the system cannot even keep stable. It can be seen from the static IO map in Figure 2 that this CSTR process is approaching unstable as the operating condition moves to a higher area, which leads to

a difficult control task. However, the main reason of poor performance is that two linear models are not able to approximate the CSTR process in the whole operating space and not able to provide enough information for the controller. Choosing γa as 0.0025 is not proper. Figure 6 shows the set-point tracking performance of the MLD-MPC controller based on the three linear models in Table 6. As can be seen from Figure 6, the input qc varies acceptably in its range [60, 115] according to changes of the set-point, and the output CA follows the reference signal swiftly, precisely, and smoothly. As the CSTR is not easy to control, to keep the output CA tracking the reference closely, every time the setpoint changes the control input qc chatters for the first few steps, which is acceptable. As a whole, the MLD-MPC controller based on three linear models performs an excellent tracking task, especially compared to the other two PID controllers. That is to say, three linear models can give sufficient information to the MLD-MPC controller for set-point tracking. Na ) 200 and γa ) 0.0025 is a nice choice. Figure 7 depicts the tracking performance of the Multi-PID controller based on the three linear models of Table 6. The output is oscillating and the tracking performance is rather poor. With the exception that this CSTR process is not easy to keep stable because of its approaching unstable as the operating conditions move to a higher area, the relatively poor coordination ability of the Multi-PID control also leads to the unsatisfactory closed-loop response. The tracking response of the nonlinear PID controller based on nonlinearity inversion is displayed in Figure 8. This tracking is still poorer than that of the Multi-PID controller. The nonlinear PID controller is designed by first getting a linear PID controller without considering the input nonlinearity, and then turning back the linear PID controller into a nonlinear one via the nonlinear inversion function. However, the physical input nonlinearity is integrated in the original control problem; handling the input nonlinearity separately therefore leads to a degraded performance. From the analyses above, it is seen that the CSTR system is not trivial to control, and the Multi-PID controller based on three linear models and the nonlinear PID controller based on nonlinearity inversion are not able to stabilize it. However, the MLD-MPC controller based on three linear models determined by the proposed included angle method excels the two PID controllers and performs satisfactorily. That is because, with a proper linear model bank, the MLD-MPC controller schedules the local linear models properly and optimizes the control variable systematically and synthetically in a unified MLD framework, thus it can both guarantee closed-loop stability and perform excellent set-point tracking control. 5.2. pH Process. Figure 9 shows the closed-loop response of the pH process, where the MLD-MPC controller is designed based on the four linear models determined by the proposed included angle dividing method with Na ) 100 and γa ) π/4.

3942 Ind. Eng. Chem. Res., Vol. 48, No. 8, 2009 Table 11. Dividing Results of Case 3 for the pH Process subregion

1st

2nd

3rd

4th

steady-state point-included operating point (q3, pH4) subrange

1 f 35 18th (5.4, 3.128) {y|2 e y e 4.10}

36 f 54 45th (13.5,6.29) {y|4.10 < y e 7.52}

55 f 58 57th (17.1,9.13) {y|7.52 < y e 9.33}

59 f 100 80th (24.0, 11.45) {y|9.33 < y e 11}

The tracking performance is poor, the input and output are both chattering frequently, and the system is even unstable. Four subregions with four local linear models cannot provide enough information for controller design. The choice of γa ) π/4 is not appropriate. The tracking performance of MLD-MPC based on five linear models for pH process is displayed in Figure 10. Input q3 varies properly in its range [0, 30] according to the variation of the reference signal, and output pH4 follows the reference signal fast, accurately, and smoothly. So, it is concluded that this pH process should be divided into five subsystems as visual intuition, and five linear models can provide sufficient information for controller design. The choice of Na ) 100, γa ) π/5 is appropriate. Figure 11 shows the performance of the Multi-PID controller. The response is not bad; the output tracks the reference quickly and smoothly. However, compared to Figure 10, the overshoot is a bit bigger and the response is a little sluggish due to relatively poor scheduling of subsystems. Figure 12 displays the tracking response of the nonlinear PID controller for the pH process. This nonlinear PID controller performs a little better than the Multi-PID controller, whereas it is eclipsed by the MLD-MPC controller because of the performance degradation as a result of nonlinearity inversion. From the previous analyses, it is concluded that the proposed included angle method is useful when dividing the operating space of a Hammerstein-like system. Choosing proper designerassigned parameters carefully, the proposed method is able to decompose the system appropriately into subsystems that provide sufficient information for multilinear model controller design. The MLD-MPC controller based on a proper linear model bank excels the Multi-PID controller based on the same model bank and the nonlinear PID controller based on nonlinearity inversion and can realize nice set-point tracking performance.

control engineers to understand; it is also easy to carry out; and moreover, it works effectively. It provides a new thought to solve the problem of linear model selection in multilinear model approach. On the basis of a proper linear model bank, the MLD-MPC strategy is employed for set-point tracking control. MLD-MPC can unite multiple linear models in a unified framework, optimize the control variable systematically, avoid oscillation during switching, and augment the stability robustness of the whole system. Closed-loop simulations demonstrate that the proposed dividing method is a useful and effective tool in decomposing the operating space for multilinear model approaches, and the MLD-MPC strategy is excellent for nonlinear systems. Although the included angle dividing method can only be used for SISO Hammerstein-like systems, and not all nonlinear systems, it is a useful method and provides a new idea for the multilinear model approach because a large number of chemical processes can be approximated by Hammerstein and Wiener models. In our further work, we may extend this method to MIMO systems to make this method more general.

6. Conclusions

(1) Pearson, R. K. Nonlinear empirical modeling techniques. Comput. Chem. Eng. 2006, 30, 1514. (2) Pearson, R. K. Selecting nonlinear model structures for computer control. J. Process Control 2001, 13, 1. (3) Pearson, R. K. Nonlinear input/output modeling. J. Process Control 1995, 5, 197. (4) Norquay, S. J.; Palazoglu, A.; Romagnoli, J. A. Model predictive control based on Wiener models. Chem. Eng. Sci. 1998, 53, 75. (5) Fruzzetti, K. P.; Palazoglu, A.; McDonald, K. A. Nonlinear model predictive control using Hammerstein models. J. Process Control 1997, 7, 31. (6) Yang, J. Study on nonlinear model predictive control algorithm based on combination model and its application, Ph.D. Thesis, Zhejiang University, China, 2007. (7) Jia, L.; Chiu, M. S.; Ge, S. S. A noniterative neruo-fuzzy based identification method for Hammerstein processes. J. Process Control 2005, 15, 749. (8) Norquay, S. J.; Palazoglu, A.; Romagnoli, J. A. Model predictive control based on Wiener models. Chem. Eng. Sci. 1998, 53, 75. (9) Harnischmacher, G.; Marquardt, W. Nonlinear model predictive control of multivariable processes using block-structured models. Contr. Eng. Pract. 2007, 15, 1238. (10) Bloemen, H. H. J.; van den Boom, T. J. J.; Verbruggen, H. B. Model-based predictiVe control for Hammerstein systems, Proc. of 39th IEEE Conference on Decision and Control, Sydney, Australia, December, 2000; pp 963-4966. (11) Murray-Smith, R.; Johansen, T. A. Multiple Model Approaches to Modelling and Control; Taylor & Francis: London, England, 1997.

An included angle dividing method is proposed to decompose the operating spaces and determine linear model banks for SISO Hammerstein-like systems through measuring the nonlinearity of the static IO curves. This method is intuitive and simple for

Figure 12. pH set-point tracking performance of nonlinear PID controller.

Acknowledgment This work has been supported by the National Natural Science Foundation of China (No. 60404018) and the State Key Development Program for Basic Research of China (No. 2009CB320603). Note Added after ASAP Publication: This paper was published ASAP on March 23, 2009. Equation 4, Tables 3 and 5, and the Acknowledgment section were each corrected. The revised paper was reposted on April 8, 2009. Literature Cited

Ind. Eng. Chem. Res., Vol. 48, No. 8, 2009 3943 (12) Chiu, M. S.; Cui, S. Internal model control design for transition control. AIChE J. 2000, 46, 309. (13) Gala´n, O.; Romagnoli, J. A.; Palazoglu, A. Robust H∞ control of nonlinear plants based on multi-linear models: an application to a benchscale pH neutralization reactor. Chem. Eng. Sci. 2000, 55, 4435. (14) Wang, F. Y.; Bahri, P.; Lee, P. L.; Cameron, I. T. A multiple model, state feedback strategy for robust control of non-linear processes. Comput. Chem. Eng. 2007, 31, 410. (15) Toscano, R. Robust synthesis of a PID controller by uncertain multimodel approach. Information Sciences 2007, 177, 1441. ¨ zkan, L.; Kothare, M. V. Stability analysis of a multi-model (16) O predictive control algorithm with application to control of chemical reactors. J. Process Control 2006, 16, 81. (17) Gala´n, O.; Romagnoli, J. A.; Palazoglu, A.; Arkun, Y. Gap metric concept and implications for multilinear model-based controller design. Ind. Eng. Chem. Res. 2003, 42, 2189. (18) Tan, W.; Marquez, H. J.; Chen, T. W.; Liu, J. Z. Multimodel analysis and controller design for nonlinear processes. Comput. Chem. Eng. 2004, 28, 2667. (19) Du, J.; Song, C.; Li, P. Modeling and control of a continuous stirred tank reactor based on a mixed logical dynamical model. Chinese J. Chem. Eng. 2007, 15, 533. (20) Li, X. Research on modeling and control of hybrid systems based on mixed logical dynamical, Ph.D. Thesis, Institute of Automation Chinese Academy of Sciences, China, 2003. (21) Heemels, W. P. M. H.; Schutter, B. D.; Bemporad, A. Equivalence of hybrid dynamical models. Automatica 2001, 37, 1085. (22) Torrisi, F. D.; Bemporad, A. et al. HYSDEL 2.0.5-User Manual; ETH, 2002. (23) Torrisi, F. D.; Bemporad, A. HYSDEL-A tool for generating computational hybrid models. IEEE Trans. Contr. Systems Technology 2004, 12, 235.

(24) Kvasnica, M.; Grieder, P.; Baotic, M.; Christophersen, F. J. MultiParametric Toolbox; ETH, December 11, 2006. (25) Hlaing, Y. M.; Chiu, M. S.; Lakshminarayanan, S. Modelling and control of multivariable processes using generalized Hammerstein model. Chem. Eng. Res. Des. 2007, 85, 445. (26) Eskinat, E.; Johnson, S. H. Use of Hammerstein models in identification of nonlinear systems. AIChE J. 1991, 37, 255. (27) Pearson, R. K.; Pottmann, M. Gray-box identification of blockoriented nonlinear models. J. Process Control 2000, 10, 301. (28) Sung, S.; Lee, J. Modeling and control of Wiener-type processes. Chem. Eng. Sci. 2004, 59, 1515. (29) Sung, S. System identification method for Hammerstein processes. Ind. Eng. Chem. Res. 2002, 41, 4295. (30) Pottmann, M.; Pearson, R. K. Block-oriented NARMAX models with output multiplicities. AIChE J. 1998, 44, 131. (31) Cervantes, A. H.; Agamennoni, O. E.; Figueroa, J. L. A nonlinear model predictive control system based on Wiener piecewise linear models. J. Process Control 2003, 13, 655. (32) Menold, P. H.; Allgower, F.; Pearson, R. K. Nonlinear structure identification of chemical processes. Comput. Chem. Eng. 1997, 21, 137. (33) Bemporad, A.; Morari, M. Control of systems integrating logic, dynamics, and constraints. Automatica 1999, 35, 407. (34) Morari, M.; Baric´, M. Recent developments in the control of constrained hybrid systems. Comput. Chem. Eng. 2006, 30, 1619. (35) Gala´n, O.; Romagnoli, J. A.; Palazoglu, A. Real-time implementation of multi-linear model-based control strategies- an application to a benchscale pH neutralization reactor. J. Process Control 2004, 14, 571.

ReceiVed for reView June 16, 2008 ReVised manuscript receiVed January 16, 2009 Accepted January 29, 2009 IE8009395