Nonlinear PI controllers based on low-order empirical process models

4668

Ind. Eng. Chem. Res. 2003, 42, 4668-4677

Nonlinear PI Controllers Based on Low-Order Empirical Process Models Robert J. Brendel, Babatunde A. Ogunnaike*, and Prasad S. Dhurjati Department of Chemical Engineering, University of Delaware, Newark, Delaware 19716-3110

Although virtually all industrially relevant processes exhibit nonlinear behavior, a significant number are still operated using linear controllers. Such controllers, especially those of the classical single-loop PI(D) type, remain popular because they are relatively easy to designs they require very modest modeling efforts, and straightforward tuning rules aboundsand they are convenient to implement using standard hardware. In this paper, we present a technique for constructing nonlinear PI-type controllers using low-order empirical nonlinear models. These controllers, by design, retain the simplicity and structure of the familiar PI controllers but are more effective over wider operating regimes and are thus capable of better performance on nonlinear processes. As part of the technique, we present a procedure for process characterization from input/ouput data, as well as explicit extensions for implementation on multi-input multioutput processes. We demonstrate the controller design and implementation procedure on simulations of industrially relevant SISO and MIMO processes. 1. Introduction The linear PID controller has long been the mainstay of industrial control practice primarily because of its versatility, ease of application, and intuitive appeal. However, its performance deteriorates in cases where the process exhibits significant nonlinear characteristics. Numerous techniques have been proposed for modifying the algorithm to deal effectively with process nonlinearities while retaining its simplicity, intuitive appeal, and familiar structure. Some of these techniques involve making the controller parameters depend on such process variables as the set point,1 feedback error,2,3 input,4 output magnitude,5 or even time.6 Perhaps the most popular technique is gain scheduling, where either the process operating regime is decomposed into discrete regions and the appropriate controller gain for each region is chosen a priori,7 or an algorithm for scheduling the gain changes is derived from a process model.8,9 A relatively recent review of research results in gain scheduling can be found in the paper by Rugh and Sharma.10 The approach developed in this paper uses a specific empirical nonlinear process model form as the basis for a direct synthesis controller design, giving rise to a PItype controller with parameters that depend on the process output. (The approach can also be used in conjunction with any other standard linear controller design technique, and the method for doing so is also discussed.) The process model form arises naturally from a Taylor series approximation of the control-affine structure often encountered in first-principles models of chemical processes. As we will show, several characteristics of this model make it well-suited for use in developing industrially applicable nonlinear controllers. A similar approach for deriving nonlinear PI and PID controllers was recently presented,11 where the authors * To whom correspondence should be addressed. Tel.: (302)-831-4504. Fax: (302)-831-1048. E-mail: ogunnaike@ che.udel.edu.

showed how such controllers can be obtained from certain first- and second-order model equations. The process models required for this approach, however, must be obtained from first principles, and if the resulting model is more complicated than the desired first- or second-order form (as is quite likely for processes of industrial significance), the modeling effort must be augmented by an additional model reduction procedure. Thus, as discussed further in their paper,11 the technique requires (i) a first-principles process model and (ii) a procedure for reducing this model to the required form using first-principles knowledge. Given that broad application of nonlinear control in industry has been hindered in large part by (i) the lack of adequate first-principles knowledge for many complex processes and (ii) the significant effort required to build such models, even when such first-principles knowledge is available, one would expect that the “activation energy barrier” to routine application of this technique in industry will be quite high. The approach presented in this paper is predicated upon retaining the spirit of the “process reaction curve”s the response of the process to a single step change in the input variableswhich has been the vehicle of choice for process characterization in tuning industrial PID controllers. The methodology we discuss in this paper is therefore based on a sufficiently simple nonlinear model that can be identified from easily obtainable process data. The rest of the paper is organized as follows: section 2 introduces the empirical model form and its properties; the controller design procedure is discussed in section 3, important issues related to identification of the model from process data are then discussed in section 4, and a possible extension to multi-input multioutput (MIMO) systems is discussed in section 5. The performance of the proposed controller design methodology is illustrated for both SISO (single-input singleoutput) and MIMO systems using a polymerization reactor simulation.

10.1021/ie0300934 CCC: $25.00 © 2003 American Chemical Society Published on Web 08/01/2003

Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003 4669

2. Empirical Input/Output Model First-principles models for many chemical processes take the following control-affine form12

dX ) F(X) + G(X)U dt Y ) H(X)

dy ) (φ1y + φ2y2) + (γ0 + γ1y + γ2y2)u(t - R) (5) dt (1)

where X is the n-vector of states, U is the m-vector of inputs, and Y is the r-vector of outputs; F, G, and H are vector mappings of appropriate dimensions. (See section 6 for two example models of a reactor used for free-radical polymerization of methyl methacrylate; the set of eqs 37 for the special case of isothermal operation are particularly explicit.) In the special case of a single measured state and a single input (n ) m ) r ) 1), eq 1 becomes

dY ) f(Y) + g(Y)U dt

The model form in eq 4swhich contains no explicit time delayssis the basis of the work described below. Note that, for processes with time delays, this model is simply modified to

where R is the process time delay. 3. Nonlinear PI Controller Design 3.1. Direct Synthesis Approach. Observe that eq 4 can be rearranged to the form

τ(y)

τ(y) )

Y - Yss y) Y* U - Uss U*

(3)

where (Uss, Yss) represents a steady-state operating point and Y* and U* are scaling factors possibly equal to Yss and Uss, respectively, we can obtain a more useful approximation of eq 2 as

dy ) (φ1y + φ2y2) + (γ0 + γ1y + γ2y2)u dt

(6)

with

(2)

which is still too general for empirical modeling. By defining scaled deviation variables y and u as

u)

dy + y(t) ) K(y) u(t) dt

(4)

This choice of an approximation to eq 2 arises directly from Taylor series expansions of f(Y), and of the product g(Y)U, around (Uss, Yss), retaining only the y2u thirdorder term and truncating the remaining third- and higher-order terms.13 (Note that, when φ2, γ1, and γ2 are zero, eq 4 reverts to the familiar linear model as a special case.) The variable scaling implied in eq 4 is important for minimizing the ill-conditioning that might arise from the cross terms in the original variables, especially when the magnitudes of U and Y differ greatly.14 Thus, Y* and U* should be chosen such that y and u are similar in magnitude. Despite its simplicity, the model structure in eq 4 is capable of exhibiting several important nonlinear characteristics such as asymmetric responses to symmetric inputs, state- and input-dependent “gains” and “time constants”, and a wide variety of stability characteristics including finite escape times.13 However, for the purposes of this current paper, the most important characteristics of eq 4 are (i) that the model form is useful for designing nonlinear PI controllers because of its structural similarity to the familiar first-order ODE and (ii) that it can be identified from process data in a manner similar to that by which the linear first-order model is obtained from the classic process reaction curve. We now discuss, in sequence, how to design nonlinear PI controllers using eq 4 and then how to identify this model from process data.

-1 φ1 + φ2y

K(y) ) τ(y)(γ0 + γ1y + γ2y2)

(7)

where the structural similarity with the standard, linear first-order model is obvious. The implication is that, by approximating process behavior with eq 4, the “effective” time constant and gain are as given in eqs 7. This structural similarity with the linear form can then be exploited for controller design by analogy with established linear methods. One specific proceduresthe direct synthesis approachsis presented explicitly next, but other techniques that use linear first-order process characterizations for controller design can be applied as well. The direct synthesis controller15 required to cause the first-order system described by

τ

dy + y(t) ) Ku(t) dt

(8)

to follow a desired first-order closed-loop response trajectory with a time constant of τR is the PI controller

[

u(t) ) Kc (t) +

1 τI

∫0t(σ) dσ

]

(9)

with

Kc )

τ KτR

τI ) τ

(10)

When time delays are involved, the expression for Kc above is modified by introducing (τR + R) in place of τR. (See Table 15.5 of ref 15.) The corresponding direct synthesis controller for the nonlinear system in eq 4 is therefore the following PI controller

[

u(t) ) Kc(y) (t) + where

∫0t(σ) dσ

1 τI(y)

]

(11)

4670 Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003

Kc(y) )

1 (γ0 + γ1y + γ2y2)τR

τI )

The same procedure can be used to obtain equivalent Cohen-Coon or Ziegler-Nichols tuning parameters.

-1 φ1 + φ2y

(12)

with τR as the only tuning parameter. Again, by analogy, when significant time delays are involved, one simply replaces τR with (τR + R). Thus, as linear PI controller design requires that the process be characterized with the first-order-plus-deadtime model and then that these characteristic parameters be used to design the controller, the technique proposed here requires initial characterization of the process with the model in eq 4, followed by use of these model parameters in the nonlinear PI controller given in eqs 11 and 12. 3.2. Other Approaches. The direct synthesis approach is particularly effective when the process output is required to move from one set point to another. Such a technique will be appropriate, for example, for polymer reactors that must transition from one operating condition to another (see the illustrative example in section 6). On the other hand, many industrial processes operate at fixed set points for extended periods of time. Under these circumstances, the control objective is disturbance rejection rather than set-point tracking, and other controller design techniques might be more appropriate. The time-integral tuning rules (ISE, IAE, and ITAE) in particular are geared toward designing controllers for disturbance rejection. By appropriately substituting the process characterization in eqs 6 and 7 above for the respective linear counterparts, other linear controller design techniques could be similarly extended for nonlinear PI designs. For example, for minimum ITAE controller tuning, the recommendations are as follows (see ref 15, Table 15.4):

0.859 τ 0.977 K R

() τ R τ ) 0.674( τ )

Kc )

0.680

(13)

I

for disturbance rejection, where R is the associated process time delay. By substituting the expressions in eq 7 for τ and K and the appropriate process time delay (or the minimum allowable value of 0.1 for R/τ where there is no significant delay; cf p 537 of ref 15), one can then obtain the corresponding minimum ITAE nonlinear PI controller parameters as

Kc )

() ( )

1 0.859 (γ0 + γ1y + γ2y2) R τI )

0.977

1 τ(y)

τ(y)0.320R0.680 0.674

0.023

(14)

4. Parameter Identification A key factor in the enduring utility of the linear PID controller in industry is that its design can be based on the process reaction curvesthe process response to the easily implementable single step change in the input variable, whose shape is assumed to be well-approximated by a theoretical first-order-plus-dead-time response. The corresponding problem for the nonlinear PI controller design technique involves identifying eq 4 from appropriately generated process data. It is well-known that input sequence design is a critical aspect of nonlinear process identification, as the pertinent process nonlinearities must be excited and reflected in the process output data if they are to be represented adequately in the resulting model.16 Therefore, the input signals to be used for identifying of the model in eq 4 should ideally be of such magnitude and frequency as to generate process data that (a) cover the entire operating range of interest and (b) reflect the inherent process nonlinearities. Clearly, the single input step, often sufficient to identify the linear first-orderplus-dead-time-model, will no longer be sufficient to identify the model in eq 4, as a single step response is incapable of adequately reflecting many important nonlinear process characteristics. Nevertheless, as shown in ref 13, a simple three-step input sequence is sufficient for adequate process characterization using eq 4. This simple input sequence permits any inherent asymmetric nonlinear steady state and dynamic nonlinearites in the operating range of interest to be capyured, it is relatively easy to implement, and the estimation of model parameters from the acquired process data is reasonably straightforward. Strategies for identifying the model parameters from input/output process data were first presented in ref 13. The parameter estimation procedure is greatly simplified by the fact that the model is linear in the parameters; the required least-squares computation is therefore straightforward and can be carried out easily using standard software packages, even spreadsheets. 4.1. Basic Identification Procedure. The specific identification method we recommend is the vint method13 (for velocity integration) discussed in detail in the reference cited and summarized here. Consider a sequence of input and output data, represented by the vectors u(t) and y(t), respectively, collected from a process over a time interval (t0, tN). Next, we divide the total time interval into segments T1, T2, ..., TM, where each Ti refers to the interval tLi e t e tUi, of length δ(Ti) ) tUi - tLi. By integrating eq 4 between these boundaries, we obtain

˜ 0(Ti) + y(tUi) - y(tLi) ) φ1z˜ 1(Ti) + φ2z˜ 2(Ti) + γ0w γ1w ˜ 1(Ti) + γ2w ˜ 2(Ti) (16)

or, without delay where

8.146 Kc ) τ(y)(γ0 + γ1y + γ2y2) τI )

τ(y) 0.310

(15)

z˜ 1(Ti) )

∫tt

z˜ 2(Ti) )

∫tt

Ui

y(t) dt

(17)

y2(t) dt

(18)

Li Ui

Li

Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003 4671

∫tt

w ˜ 0(Ti) )

Ui

Li

u(t) dt

(19) (20)

w ˜ 1(Ti) )

∫tt

y(t) u(t) dt

w ˜ 2(Ti) )

∫t

y (t) u(t) dt

Ui

Li

(φ1 + γIu*)2 - 4γ0u*(φ2 + γ2u*) > 0

tUi 2

Li

˜ 0(Ti) + y(tUi) - y(tLi) ) φ1z˜ 1(Ti) + φ2z˜ 2(Ti) + γ0w γ1w ˜ 1(Ti) + γ2w ˜ 2(Ti) + (Ti) (22) where (Ti) is the model error. This is an equation of the form

y˜ ) X ˜Θ + 

[

z˜ 1(T1) z˜ (T ) X ˜ ) 1 2 l z˜ 1(TM)

z˜ 2(T1) z˜ 2(T2) l z˜ 2(TM)

w ˜ 0(T1) w ˜ 0(T2) l w ˜ 0(TM)

w ˜ 1(T1) w ˜ 1(T2) l w ˜ 1(TM)

(27)

(21)

The model identification problem can therefore be posed as follows: Given input/output data u(tk), y(tk), respectively, for k ) 1, 2, ..., i, ..., M, obtain the model parameters using the equation

where

When the input u is a constant u*, a similar analysis results in the stability criterion that the following discriminant must be positive

which reduces to the requirement that φ1 be real for u ) 0. The discriminant also appears in a calculation of the steady-state outputs y* predicted by the empirical model for a given constant input u*. Setting the derivative in eq 4 equal to zero results in two possible values for the steady-state output / ) y1,2

-(φ1 + γ1u*) ( x(φ1 + γ1u*)2 - 4γ0u*(φ2 + γ2u*) 2(φ2 + γ2u*)

(28)

(23)

w ˜ 2(T1) w ˜ 2(T2) l w ˜ 2(TM)

]

(24)

and Θ is the vector of desired model parameters. This is clearly a standard linear regression problem for which any number of standard software packages can be used. It should be noted that the cross terms in eq 4 require that careful consideration be given to variable scaling. If the numerical values of u and y differ by orders of magnitude, the parameter estimation problem shown above can become ill-conditioned. As mentioned earlier, this potential problem can be addressed by judicious choice of the scaling values Y* and U* in eq 3. The steady-state values Yss and Uss have been used with success in a number of process applications (see section 6 below for an example). For the remainder of this paper, the variables u and y are taken to be properly scaled. 4.2. Additional Model Parameter Considerations. The nature of the process model in eq 4 and the resulting controller parameter structures to which it gives rise, raise some additional subtle issues that we must now consider. The most critical of these, as we now show below, is that it is possible to identify an unstable process model inadvertently even though the physical process itself is stable. For an input u ) 0, the empirical model becomes

dy ) φ1y + φ2y2 dt

(25)

When integrated from the initial condition y(t)0) ) y0, the trajectory of the output y(t) is the solution to

y0 y eφ1t ) 0 φ1 + φ2y φ1 + φ2y0

A negative discriminant gives a nonzero imaginary component to the steady-state output, which results in an unstable empirical model. Ordinarily then, the potential exists that the model stability characteristics might not match those of the actual process if care is not taken in identifying model parameters from data. Thus, for processes known to be stable, it might be necessary to modify the parameter identification technique we have presented thus far if one wishes to guarantee that the identified model is stable over the range spanned by the open-loop identification data. As shown above, this involves no more than the imposition of eq 27 as a constraint on the standard least-squares estimation problem. Keep in mind, however, that, by imposing the nonlinear constraint required to guarantee model stability, one converts an otherwise straightforward procedure into one that requires nonlinear optimization. Thus, the standard procedure should be employed first and the identified model stability characteristics validated against that of the actual process; only if a discrepancy exists in this crucial aspect should the model identification be revisited as a constrained optimization problem. If one wants to avoid constrained optimization entirely, it is possible to employ an alternative, suboptimal procedure. This sequential procedure involves first estimating a time constant for each input step, from which one then determines the φi parameters. These parameters are then fixed, and the entire data set is used for identification of the remaning γi parameters via least squares. This procedure is illustrated in section 6 where it is applied to a multivariable problem.

(26)

Observe, therefore, that the only requirement for stability of the model when u ) 0 is that φ1 must be negative. This is equivalent to the stability requirement that τ must be positive in the first-order linear model in eq 8.

5. Extension to Multi-Input Multi-Output Systems The empirical model form in eq 4 can be extended to multivariable systems in several different ways. The most straightforward extension with the simplest structure for a 2 × 2 process is

dy1 ) (φ11y2 + φ12y22) + (γ110 + γ111y2 + γ112y22)u1 + dt (γ120 + γ121y2 + γ122y22)u2

4672 Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003

dy2 ) (φ21y2 + φ22y22) + (γ210 + γ211y2 + γ212y22)u1 + dt (γ220 + γ221y2 + γ222y22)u2 (29)

gI1 ) gI2 ) -

Γ12(y1) Γ11(y1) Γ21(y2)

(34)

where the two outputs are y1 and y2 and the two inputs are u1 and u2. It is quite straightforward to extend the least-squares technique presented for the scalar problem to this multivariable process model. In a vein similar to the preceding discussion, it is easy to show that the MIMO model is stable when the following discriminants are positive

as the expressions required of the functions gI1 and gI2 to acheive loop decoupling. In combination with eqs 31, these functions define a simplified decoupling controller based on the nonlinear empirical model form in eqs 29 (cf. eqs 22.20 and 22.21 of ref 15).

(φ11 + γ111u/1 + γ121u/2)2 - 4(γ110u/1 + γ120u/2)(φ12 +

6. Illustrative Example

γ112u/1

+

γ122u/2)

(φ21 + γ211u/1 + γ221u/2)2 - 4(γ210u/1 + γ220u/2)(φ22 + γ212u/1 + γ222u/2) (30) for constant u/1 and u/2. Distributed (single-loop) nonlinear PI controllers can now be constructed in precisely the same manner as described above for the SISO case. The controller for a y1-u1 pairing, for example, would use the parameters φ11, φ12, γ110, γ111, and γ112 to construct the gain and integral time values as functions of y1 precisely as was done for the SISO process, with each distributed controller assigned its own unique value for τR, the tuning parameter. When time delays are involved, u1 and u2 are replaced respectively by u1(t - R11) and u2(t - R12) in the first equation and by u1(t - R21) and u2(t - R22) in the second equation. If necessary, the empirical model in eq 29 can be used to design decouplers to compensate for loop interactions along the same lines as with linear systems, as we now illustrate. Let v1 and v2 be the outputs of a single-loop nonlinear PI controller while u1 and u2 are the actual process inputs (cf. Chapter 22 of ref 15). With interaction compensators gI1 and gI2, the control action is determined according to

Daoutidis and co-workers17 presented the following model for the free-radical polymerization, in a jacketed CSTR, of methyl methacrylate in toluene solvent with azo-bis-isobutyronitrile as the initiator

[

( )]

( )

-Efm -Ep dCm + Zfm exp CmP0(CI,T) + ) - Zp exp dt RT RT F(Cmin - Cm) V

( )

FICIin - FCI dCI -EI CI + ) -ZI exp dt RT V

[

( ) ( )]

-ETc dD0 + ZTd ) 0.5ZTc exp dt RT -ETd exp [P0(CI,T)]2 + Zfm RT -Efm FD0 exp CmP0(CI,T) RT V

[

( )

( )

-Ep dD1 + Zfm ) Mm Zp exp dt RT -Efm FD1 exp CmP0(CI,T) RT V

( )]

( )

-Ep (-∆Hp) dT ) Zp exp Cm P0(CI,T) dt RT Fcp

u1 ) v1 + gI1v2 u2 ) v2 + gI2v1

Γ22(y2)

F(Tin - T) UA (T - Tj) + FcpV V

(31)

If eqs 29 are rewritten in the more compact form

dTj Fcw UA ) (T - Tj) + (T - Tj) dt Vo w o FwcwVo

dy1 ) Φ1(y1) + Γ11(y1)u1 + Γ12(y1)u2 dt dy2 ) Φ2(y2) + Γ21(y2)u1 + Γ22(y2)u2 dt

where

(32) P0(CI,T) )

and combined with eqs 31, the result is

from which one obtains

( )

ZTd exp

dy1 ) Φ1 + (Γ11 + Γ12gI2)v1 + (Γ11gI1 + Γ12)v2 dt dy2 ) Φ2 + (Γ21 + Γ22gI2)v1 + (Γ21gI1 + Γ22)v2 dt

[

2f*CIZI exp

(33)

-ETd RT

( ) -EI RT

+ ZTc exp

(35)

)]

0.5

(

-ETc RT

(36)

Cm and CI are the monomer and initiator concentrations in the reactor, respectively; D0 and D1 are molar and mass concentrations respectively, of “dead” polymer chains; T and Tj are temperatures of the reaction mixture and jacket water, respectively; FI and Fcw are the flow rates of initiator and jacket cooling water,

Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003 4673 Table 1. Process Model Parameters parameter

value

units

ZTc ZTd ZI ZP Zfm R CMin CIin F U Tin Fw Vo

3.8223 × 3.1457 × 1011 3.7920 × 1018 1.7700 × 109 1.0067 × 1015 8.314 6.0 8.0 866 720 350 1,000 0.02

m-3

1010

h-1

kmol kmol m-3 h-1 h-1 kmol m-3 h-1 kmol m-3 h-1 kJ kmol-1 K-1 kmol m-3 kmol m-3 kg m-3 kJ h-1 K-1 m-2 K kg m-3 m3

parameter

value

units

ETc ETd EI EP Efm F V ∆H cP A Two cw

2.9442 × 2.9442 × 103 1.2550 × 105 1.8283 × 104 7.4478 × 104 1.00 0.1 -57 800 2.0 2.0 293.2 4.2 103

kJ kmol-1 kJ kmol-1 kJ kmol-1 kJ kmol-1 kJ kmol-1 m3 m3 kJ kmol-1 kJ kmol-1 K-1 m2 K kJ kmol-1 K-1

Table 2. Steady-State Operating Conditions and Scaling Factors parameter

value

units

parameter

value

units

Cm CI D0 D1 T Tj

5.5068 0.1329 0.001 975 2 49.382 335 297.2

kmol m-3 kmol m-3 kmol m-3 kg m-3 °C °C

FI Fcw Y /1 Y /2 U /1 U /2

0.016 79 3.2636 25 000 10 0.016 79 3.2636

m3 h-1 m3 h-1 kg kmol-1 °C m3 h-1 m3 h-1

respectively; V and Vo are the reactor and jacket volumes, respectively; U and A are the heat-transfer coefficient and area, respectively; and F and c are the density and heat capacity, respectively. The subscript in refers to inlet stream values. The heat-transfer fluid is water. Process parameter and steady-state values are listed in Tables 1 and 2. Additional details, including references and modeling assumptions, are available in ref 17. The state and input/output variables can be defined as follows: X1 ) Cm, X2 ) CI, X3 ) D0, X4 ) D1, X5 ) T, X6 ) Tj, Y1 ) X4/X3, Y2 ) X5, U1 ) FI, and U2 ) Fcw. 6.1. Single-Input Single-Output Case. For the special case of isothermal operation, the model reduces to the following form18

Figure 1. Open-loop step responses of the SISO polymerization reactor. Solid lines, (50% input step changes; dashed lines, (10% input step changes.

dX1 ) 10(6 - X1) - 2.4568X1xX2 dt dX2 ) 80U - 10.1022X2 dt dX3 ) 0.002 412 1X1xX2 + 0.112 191X2 - 10X3 dt dX4 ) 245.978X1xX2 - 10X4 dt Y)

X4 X3

U ) FI

(37)

Note that both eqs 35 and 37 are in the form of eq 1. The chosen base operating conditions are the steadystate conditions shown in Table 2. Scaled deviation variables are defined as

y)

Y - 25 000.5 25 000.5

u)

U - 0.016 79 0.016 79

(38)

Figure 2. Polymerization reactor model identification and validation. Solid line, first-principles kinetic model; dashed line, empirical model. Data from 0-3 h used for identification; data from 3-5 h used for validation.

The nonlinear character of this process is evident from Figure 1, which shows open-loop responses of the output (number-average molecular weight) to step changes in the input (initiator flow). Observe that the responses are nearly symmetric for small input steps but asymmetric for larger steps. Model Identification and Validation. Figure 2 shows the sequence of input steps used for empirical

4674 Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003

Figure 3. Identified model stability characteristics: (a) discriminant, (b) predicted number-average molecular weight loci.

Figure 4. Nonlinear PI controller parameter variations over operating region: (a) Kc(y)τR, (b) τI(y).

Table 3. Polymerization Reactor Empirical Model and Nonlinear PI Controller Parameters parameter

value

φ1 φ2 γ0 γ1 γ2 KcτR τI τR

-2.895 -2.355 -1.385 -2.943 -1.224 -3.4 × 10-7 0.3 h 0.3 h

model identification, along with the “true” process output and the empirical model prediction. Data from the first three of the five steps were used for model identification (using the technique outlined in section 4), whereas those from the last two steps were reserved for model validation. The operating range for this process is given as Y ) [15 000 35 000] or y ) [-0.4 0.4]; note that this range is spanned by the model identification experiments. The sample time was 30 s, giving 360 input/output data points available for model identification. The estimated empirical model parameter values are listed in Table 3. It is important now to confirm that the identified empirical model is stable over the input range of interest. Figure 3a shows that the discriminant from eq 27 is positive; also, Figure 3b shows the existence of two real, distinct steady-state solutions for the polymer number-average molecular weight for all initiator flow rates in the range of interest. Of course, one solution, even though real, is negative and therefore physically meaningless. This branch can thus be ignored, leaving the upper curve as the only relevant predicted steadystate number-average molecular weight. Controller Design and Closed-Loop Performance. Over the operating range of interest, the direct synthesis nonlinear PI controller constructed from the identified model exhibits the controller gain and integral time characteristics shown in Figure 4. Closed-loop simulation responses for set-point tracking are shown in Figure 5. Included for comparison are results for the same tests with a linear PI controller that was also tuned using the direct synthesis method. (The first-order linear process model parameters were determined graphically from the single step response process reaction curve.) The linear PI controller con-

Figure 5. Closed-loop polymerization reactor responses. Solid line, nonlinear PI; dashed line, linear PI; dotted line, set point.

stants Kc and τI are included in Table 3. The value of τR, the desired closed-loop response time constant, was 0.3 h for both linear and nonlinear controllers. The controller implementation was the digital “velocity” form to avoid reset windup. The simulation environment was Matlab v5.3 and 6.0. Observe that, in the operating range Y ) [25 000 35 000] (where the linear controller model was identified), the performances of the two controllers are virtually identical. Outside this region, however, the linear controller response is somewhat more sluggish, especially for low values of the process output. The nonlinear PI controller parameters, on the other hand, change with operating conditions, as shown in Figure 4, resulting in much faster closed-loop responses and, hence, better performance. The nonlinear PI controller also performed well for disturbance rejection, even though the direct synthesis technique is intended primarily for set-point tracking. At 15 h, the monomer concentration was decreased by one-third in a step, and at 20 h, it was restored to its original value. As seen in Figure 5, the nonlinear PI controller is able to return the output to the set point more quickly than the linear PI controller.

Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003 4675

Figure 6. Polymerization reactor robustness tests: -50% error in a critical process parameter. Solid lines, nonlinear PI; dashed lines, linear PI; dotted line, set point.

In principle, the question of a separate investigation of controller performance in the face of plant/model mismatch does not really arise here, even though this is a simulation study. This is because the controller that was implemented on the full first-principles process model was designed on the basis of a simplifed empirical model identified from data (albeit simulated). Nevertheless, we performed a more stringent test of the controller’s robustness by introducing a -50% error in ZI, the physical parameter corresponding to the preexponential factor of the rate constant associated with the initiation reaction. As is well-known, the dynamic behavior of polymer reactors of this type is heavily influenced by the initiator concentration, which, in turn, is very sensitive to the value of ZI. A -50% change in this parameter constitutes a severe parametric error. The closed-loop test was repeated, implementing the linear and nonlinear controllers as previously designed, but this time, on the “plant” with the altered value for ZI. The result is shown in Figure 6, where there is only a slight degradation in controller performance, with the nonlinear controller maintaining its performance edge over the linear controller. 6.2. Multi-Input Multi-Output Case. In designing effective control systems for multivariable processes, the engineer has several options. In the simplest cases, multiple single-loop controllers implemented with appropriate loop pairing sometimes provide adequate performance. The advantage of this choice is its simplicity. When interactions are significant, the performance degradation can be compensated through the introduction of decoupling. When time delays, constraints, and other complex dynamics are involved and the interactions are significant, model predictive control (MPC) has become the technique of choice;19,20 the well-documented performance improvement, however, comes at the cost on increased effort. The control engineer’s decision is therefore often one of trading simplicity for performance. In keeping with the spirit of this paper, therefore, we now illustrate how the techniques presented earlier can be employed either for designing multiple single-loop controllers or for augmenting them with decoupling. When performance is more important that simplicity, of course, MPC (linear or nonlinear) is recommended.

Figure 7. Open-loop identification data. Solid lines, y1 and u1; dashed lines, y2 and u2. Table 4. 2 × 2 System Empirical Model Parameters parameter

value

parameter

value

φ11 φ12 γ110 γ111 γ112 γ120 γ121 γ122

-2.541 -2.348 -3.084 -3.769 -0.7116 0.6499 1.410 -1.788

φ21 φ22 γ210 γ211 γ212 γ220 γ221 γ222

-2.951 0.3020 4.644 7.024 7.210 -1.299 -1.174 -0.4482

Simple MIMO Extensions. The model identification procedure of section 4 was extended for the multivariable case and modified somewhat to help ensure stability of the identified models. The input and output variables were defined as follows: Y1 ) X4/X3, Y2 ) X5, U1 ) FI, U2 ) Fcw, y1 ) (Y1 - 25 000.5)/25 000.5, y2 ) (Y2 - 335)/10, u1 ) (U1 - 0.016 79)/0.016 79, u2 ) (U2 - 3.2636)/3.2636. The scaling factors were chosen such that the scaled deviation variables were all of the same order of magnitude. In all cases except for y2, the choice was the respective steady-state value. The multivariable empirical model parameters were identified in two steps. First, the time constants were determined for each input step according to established techniques. Then, the φ parameters were identified by least squares, using the average output value for that step. For step i

-1 ) φ11 + φ12y1h ,i τ1,i -1 ) φ21 + φ22y2h ,i τ2,i

(39)

The γ parameters were then identified using leastsquares methods as discussed previously. Figure 7 shows the input and output trajectories used in identifying a model for the multivariable reactor. The identified model parameters are reported in Table 4. Figure 8 shows a comparison of the model predictions and the plant output, indicating a resonably good fit. The identified empirical model predicts a region of instability for output y2, as shown in Figure 9, where the square root of the discriminant has a nonzero imaginary part. Note that the predicted instability

4676 Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003

Figure 8. Model prediction vs actual outputs. Solid lines, kinetic model; dashed lines, empirical model.

Figure 10. SISO controller parameter variations over the operating range.

Figure 9. Imaginary part of y˘ 2 discriminant.

occurs at high initiator flows (u1) and low cooling water flows (u2), which corresponds to known unstable operating conditions for this process. Thus, the stability characteristics of the empirical model match those for the actual full-scale model. Closed-Loop Performance. The performance of distributed (single-loop) SISO controllers based on the identified empirical model was evaluated with the following configuration: The initiator flow (FI) was used to control product molecular weight (D1/D0), and the jacket cooling water flow (Fcw) was used to control the reactor temperature (T), for a u1-y1, u2-y2 pairing. How the nonlinear controller parameters change over the operating range is shown in Figure 10. Note that the changes in gain and effective time constant are significant only for loop 1, suggesting that closed-loop performance can likely be improved over linear PI controller performance only for this loop. Figure 11 shows a comparison of the closed-loop responses for distributed nonlinear and linear PI controllers with tuning parameters as listed in Table 5. The linear PI parameters are those suggested by Daoutidis and co-workers.17 At 0.5 h, the temperature set point was changed to 340 from 335 °C, and at 5.0 h, the number-average molecular weight set point was reduced to 20 000 from 25 000 kg/kmol. At 8.0 h, the inlet-stream monomer concentration was changed from 6.0 to 5.0

Figure 11. Closed-loop responses of 2 × 2 polymerization reactor. Solid line, distributed nonlinear PIs; dot-dashed lines, distributed linear PIs; dashed line, decoupled nonlinear PIs; dotted line, set points. Table 5. Distributed PI Controller Tuning Parameters parameter Kc1 τI1 Kc2 τI2 τR,1 τR,2

value 10-7

-1 × 0.075 -1 × 10-1 0.075 0.3 2.5

units kmol kg-1 h m3 h-1 K-1 h h h m3

h-1

kmol/m3, and at 10.0 h, the inlet temperature was decreased from 350 to 345 °C. The response of the nonlinear PI controllers is somewhat better than that of the linear PI controllers for both set-point tracking and disturbance rejection only for loop 1, as we had suspected from an inspection of Figure 10. To complete the simulation study, we applied simplified decoupling to the 2 × 2 polymerization reactor. We note first the following about the closed-loop system behavior: The number-average molecular weight (y1) response was significantly affected as the cooling water flow rate (u2) was being manipulated to control the temperature (y2). In contrast, the interactive effect of the initiator flow rate (u1) on temperature was not as

Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003 4677

significant. Thus, we considered implementing only partial decoupling, with the decoupling function gI1 implemented as in eq 34 while gI2 was set to zero. The resulting closed-loop responses are included in Figure 11, where the partial decoupling did not appear to offer any significant improvement in the overall performance. Decoupling, of course, is well-known to be sensitive to model errors; thus, as demonstrated by this example, closed-loop performance can be improved by decoupling only by upgrading model fidelity. If the distributed nonlinear controller performance is acceptable, decoupling might not be necessary. If, however, none of the simpler techniques illustrated here provides acceptable performance for the multivariable system and it becomes necessary to upgrade model fidelity, the control engineer might consider this to be sufficient justification for investing in the increased controller design effort required for implementing a full-scale nonlinear MPC. 7. Summary and Conclusion We have developed a nonlinear PI controller design technique based on a low-order, nonlinear empirical model. Parameter identification from a simple input sequence and appropriate nonlinear modification of linear controller designs was shown to be sufficient to provide improved closed-loop performance relative to the traditional linear PI controller. One process example was used to illustrate controller design and performance for both SISO and MIMO systems. The controller was demonstrated to be particularly effective in dealing with the important nonlinearity of asymmetric responses to symmetric inputs that is common in many chemical processes, especially for polymer reactors that are required to transition from one operating condition to another. The discussion of extensions to MIMO systems was restricted to the simplest of all possible multivariable model structures. Determining which multivariable model structure is most appropriate is still an open question. However, we believe that it is possible to obtain improved multivariable performance by basing controller design on a model structure that better fits the observed process behavior. Literature Cited (1) Rugh, W. J. Design of Nonlinear PID Controllers. AIChE J. 1987, 33, 1738-1742. (2) Cheung, T.-F.; Luyben, W. L. Nonlinear and Nonconventional Liquid Level Controllers. Ind. Eng. Chem. Fundam. 1980, 19, 93-98. (3) Jutan, A. A Nonlinear PI(D) Controller. Can. J. Chem. Eng. 1989, 67, 485-493.

(4) Wang, J. L.; Rugh, W. J. Parametrized Linear Systems and Linearization Families for Nonlinear Systems. IEEE Trans. Circuits Syst. 1987, 34, 650. (5) Marini, L.; Georgakis, C. Low-Density Polyethylene Vessel Reactors, Part II: A Novel Controller. AIChE J. 1984, 30, 409415. (6) Sung, S. W.; Lee, I.-B.; Lee, J. Modified ProportionalIntegral-Derivative (PID) Controller and a New Tuning Method for the PID Controller. Ind. Eng. Chem. Res. 1995, 34, 4127-4132. (7) Rugh, W. J. Analytical Framework for Gain Scheduling. In Proceedings of the American Control Conference; IEEE Press: Piscataway, NJ, 1990; pp 1688-1694. (8) Bequette, B. W. Nonlinear Control of Chemical Processes: A Review. Ind. Eng. Chem. Res. 1991, 30 (7), 1391-1413. (9) Bequette, B. W. Practical Approaches to Nonlinear Control: A Review of Process Applications. In Nonlinear Model Based Control; Berber, R., Kravaris, C., Eds.; NATO-ASI Series E.; Kluwer: Dordrecht, The Netherlands, 1998; Vol. 353. (10) Rugh, W. J.; Sharma, J. S. Research on gain scheduling. Automatica 2000, 36, 1401-1425. (11) Wright, R. A.; Kravaris, C.; Kazantzis, N. Model-Based Synthesis of Nonlinear PI and PID Controllers. AIChE J. 2001, 47 (8), 1805-1818. (12) Kantor, J. C. An Overview of Nonlinear Geometric Methods for Process Control. In Shell Process Control Workshop; Butterworths: Stoneham, MA, 1987; pp 225-250. (13) Ogunnaike, B. A.; Pearson, R. K.; Samardzija, N.; Bomberger, J. D. Low Order Empirical Modeling for Nonlinear Systems. In Advanced Control of Chemical Processes (ADCHEM ′94): IFAC Symposium; Pergamon Press: New York, 1994; pp 41-47. (14) Skogestad, S. Consistency of Steady-State Models Using Insight Extensive Variables. Ind. Eng. Chem. Res. 1991, 30, 654661. (15) Ogunnaike, B. A.; Ray, W. H. Process Dynamics, Modeling, and Control; Oxford University Press: New York, 1994. (16) Pearson, R. K. Nonlinear Input/Output Modeling. J. Process Control 1995, 5, 197-211. (17) Daoutidis, P.; Soroush, M.; Kravaris, C. Feedforward/ Feedback Control of Multivariable Nonlinear Processes. AIChE J. 1990, 36, 1471-1484. (18) Doyle, F. J.; Ogunnaike, B. A.; Pearson, R. K. Nonlinear Model Predictive Control Using Second-Order Volterra Models. Automatica 1995, 31, 697-714. (19) Qin, S. J.; Badgwell, T. A. An Overview of industrial model predictive control technology. In Chemical Process Control V; Kantor, J. C., Garcia, C. E., Carnahan, B., Eds.; CACHE: Austin, TX, 1997; pp 232-256. (20) Qin, S. J.; Badgwell, T. A. An Overview of nonlinear model predictive control applications. In International Symposium on Nonlinear Model Predictive Control: Assessment and Future Directions; Algo¨wer, F., Zheng, A., Eds.; Birkhauser Verlag: Basel, Switzerland, 1998; pp 128-145.

Received for review February 6, 2003 Revised manuscript received May 19, 2003 Accepted May 22, 2003 IE0300934