4506
Ind. Eng. Chem. Res. 2003, 42, 4506-4511
A Steady-State Process Resiliency Index for Nonlinear Processes: 1. Analysis Boris M. Solovyev and Daniel R. Lewin* PSE Research Group, Wolfson Department of Chemical Engineering, Technion-Israel Institute of Technology (IIT), Haifa 32000, Israel
Preliminary process resiliency analysis is a necessary part of any modern process design methodology. This paper extends the disturbance cost (DC) resiliency index, previously developed for linear plants, to nonlinear systems. As the nonlinear index can handle such nonlinear features as input and output multiplicities, it can accommodate the analysis of the potential of a partial control strategy to improve process resiliency. The index reliably identifies potential resiliency problems associated with process design. This paper deals with the DC analysis of processes that can be described by closed-form models and demonstrates its application to the design of a continuous stirred tank reactor (CSTR). 1. Introduction Traditionally, continuous chemical process design has been carried out at steady state for a given operating range, driven by the desire to minimize costs, largely ignoring controllability and resiliency. However, design on the basis of steady-state economics alone is risky because the resulting plants are often difficult to control, leading to off-specification product, excessive use of fuel, and associated profitability losses. Controllability can be defined as the ease with which a continuous plant can be held at a specific steady state. Resiliency measures the degree to which a processing system can meet its design objectives despite external disturbances and uncertainties in its design parameters. Because adequate process resiliency is a necessary part of optimal process design (Luyben et al.,1 Seider et al.2), it is desirable to consider process resiliency assessment when determining the process structure and establishing the operating range. A number of attempts to quantify process resiliency have been reported previously, most of which are based on linear process models (e.g., Skogestad and Morari,3 Stanley et al.,4 Wolff,5 Lewin,6 Cao et al.,7 and Solovyev and Lewin8). However, the actual process might exhibit highly nonlinear behavior (Solovyev and Lewin8), making the linear analysis less effective. Several approaches intended to overcome the limitations of linear analysis have been reported recently (e.g., Seferlis and Grievink9 and Subramanian and Georgakis10). This paper proposes an alternative approach to perform resiliency evaluation by using the “disturbance cost” (DC) philosophy (Lewin,6 Solovyev and Lewin8). This approach defines the resiliency measure in terms of the control effort required to perfectly reject the worst-case disturbance, normalized with respect to the available range of control action afforded by the designed process. Consequently, a DC value of unity indicates control system saturation, and thus, a process * To whom correspondence should be addressed. Tel.: +972-4-8292006. Fax: +972-4-8295672. E-mail: dlewin@ tx.technion.ac.il. URL: http://tx.technion.ac.il/∼dlewin/ pse.htm.
design that is characterized by DC values in excess of unity should be modified or discarded in favor of a more readily controlled alternative. Here, an extension of the DC methodology to nonlinear cases is proposed and demonstrated in the resiliency analysis of a continuous stirred tank reactor (CSTR) based on a first-principles model (Subramanian and Georgakis10). Transient nonlinear simulations verify the obtained results. 2. Disturbance Cost Definition The disturbance cost (DC) is defined as the control effort required to reject the worst disturbance. Mathematically, it can be defined as DC ) max min max {max[R ) (ui - uLi )/∆umax , 1 - R]} i d∈D
u∈U i)1,...,nu
subject to
F(x,u,d,v) ) 0 y ) G(x,u,d,v) y∈Y
(1)
In eq 1, D, U, and Y are the spaces of allowable disturbances/uncertainties, inputs, and outputs, respectively; y is a vector of process outputs; x is the state vector; u is a vector of nu manipulated (input) variables; d is a vector of disturbances; v is a vector of design variables; F is a vector function describing the process; and G is a vector function giving the outputs as a function of the states, manipulated variables, and represents the maximum allowed disturbances. ∆umax i change in input i, superscript L indicates the lower limit, and R is a temporary variable. The symbol max indicates the maximum value of a function, and max and min stand for the maximization and minimization procedures, respectively. The expression max[R ) (ui , 1 - R] in eq 1 represents the normalized uLi )/∆umax i control effort. When the input set that returns the process to the design operating region Y is not unique, it must be minimized over the U domain. It is important to note that multiple input sets that adequately return the process to the desired operating region can exist as a result of process nonlinearity (input multiplicity),
10.1021/ie020926a CCC: $25.00 © 2003 American Chemical Society Published on Web 05/06/2003
Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003 4507
unequal numbers of inputs and outputs (MPC), or partial control requirements. In this context, partial control ideas are of great practical use in keeping the process in the desired operating region rather than at the desired (nominal) set point, especially when there are more process outputs than inputs (Arbel et al.11). The maximization of the control effort on the D domain provides a condition for the worst-case disturbance. A consequence of the normalization in eq 1 is that a DC value of unity indicates control system saturation. Thus, a process design characterized by DC values in excess of unity should be modified or replaced by a more resilient alternative, noting that values of DC below unity imply designs that are easier to control. Note that mathematical programming could solve the inner minimization problem. The infinity norm of an array of linear expressions produces a concave shell whose edges are defined by arrays members. The solution lies on the intersection of the shells ribs with the space defined by the active constraints. In contrast, the outer maximization problem is less convenient to solve by programming, as it requires substituting the original min-max problem by a dual min or max problem over the united U ∪ D domain. A more attractive alternative is to discretize the D domain and then to solve the inner minimization problem for discrete values of disturbances. This is the method implemented in this study, which provides a measure of the process nonlinearity and monotonicity with respect to the disturbances. Using this approach, it is possible to distinguish most critical regions of the D domain and to reduce significantly the dimensionality of the problem. The aforementioned DC definition can be reduced to that used for the linear case (Solovyev and Lewin8), assuming that the following three conditions hold: (1) the processes are linear and square (i.e., number of inputs equal to the number of outputs), (2) the nominal values of all disturbances are located in the middle of their ranges (i.e., symmetry), and (3) the process is maintained at the nominal operating point (set-point control). Conditions 1 and 3 imply that the input set is unique; thus, optimization on the U domain in eq 1 is replaced by the direct computation of the set of inputs
y˜ ) Pu˜ + Pdd ˜ ) 0 w u˜ ) -P-1Pdd ˜
d∈D
{[|(P P ) d|/∆u P ) d|/∆u , ..., |(P
|(P
-1
d 2
-1
d 1
max 2
max 1 ,
-1
|
Pd)nud /∆unmax u
|
(|
|
| )
-1
|
..., P Pd
|
|
nu∆d
max
) max P-1Pd ∆dmax/∆umax
(|
/∆unmax u
)
(4)
The evaluation of the DC relies, first, on the availability of a process model and, second, on the ability to solve this model. Considering nondifferentiability and the possibly nonconvex structure of the goal function, the evaluation of DC might be a complex computational task. However, the recent growth of computational power on one hand and the progress in algorithmic development on the other allow us to attack this complex problem. As an illustration, the proposed resiliency measure is demonstrated on a first-principles model of a CSTR. It is noted that modeling and optimization are carried out using Mathematica (Wolfram12). The solution of algebraic equations is performed using the modified Newton method, available in Mathematica, while optimization is carried out using the multistart hill-climbing method, available in the Global Optimization package written for Mathematica by Loehle Enterprise (Loehle13). 3. Resiliency Analysis of a CSTR (Subramanian and Georgakis10) A three-state CSTR model (perfectly mixed exothermic CSTR with the single irreversible reaction AfB) is described by the following material and energy balances
d(VCA) ) Q0CA0 - QCA - Vk0e-Ea/RT dt dV ) Q0 - Q dt (5) d(VT) ∆H UA ) Q0T0 - QT Vk e-Ea/RT (T - TC) dt FCp 0 FCp dTC QC UA ) (T - TC) + (T - TC) dt VC C0 VCFCCpC
(2)
Here, x˜ indicates a deviation variable, i.e., x˜ ) x - xnom, and P and Pd are process transfer matrices relating process outputs to process inputs and disturbances, respectively. In this case, the computation of the DC is reduced to following maximization problem
DC ) max max
DC ) -1 max max P-1Pd 1∆dmax/∆umax /∆umax 1 , P Pd 2∆d 2 ,
]}
In the above equations, A is the heat-transfer area, CA is the concentration at the reactor outlet, Cp is the heat capacity, Ea is the activation energy, ∆H is the heat of reaction, k0 is the preexponential factor, Q is the volumetric flow rate, R is the gas constant, T is the reactor temperature, TC is the cooling-jacket temperature, U is the heat-transfer coefficient, V is the volume, and F is the density. Note that the indices 0 and C represent feed and coolant properties, respectively. The following dimensionless variables are defined
(3)
The indices 1, 2, ..., nu symbolize appropriate rows of the matrix (P-1Pd), and |x| denotes the absolute value of x. The assumption of symmetry leads to the simplification of eq 3 to the form defined by Solovyev and Lewin8
x1 )
CA CA,ref
x10 )
CA0 CA,ref
x2 )
T - Tref Tref
x20 )
T0 - Tref Tref
4508 Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003
TC - Tref x3 ) Tref
TC0 - Tref x30 ) Tref
Ea γ) RTTref
k0Vrefe-γ φ) Qref
β)
∆HCA,ref TrefF Cp
Vref δ1 ) VC
Table 1. CSTR Model Parameters (from Subramanian and Georgakis10)
β δ2 γ φ δS δB
V UA δS + δB ) Vref FCpQref FCp δ2 ) FCCpC
u2 )
V Vref
q)
Q Qref
q0 )
Q0 Qref
u1 )
qC )
QC,ref Qref
τ)
QC QC,ref
tQref Vref
Substituting the dimensionless variables into eqs 5 and expressing the model in the steady state gives
q0(x10 - x1) - u2φ x1f(x2) ) 0 q0(x20 - x2) - (x2 - x3)(δB + u2δS) + u2x1βφf(x2) ) 0 δ1[qCu1(x30 - x3) + (x2 - x3)(δB + u2δS)δ2] ) 0
(6)
where f(x2) ) eγx2/(x2+1). The process outputs are the conversion, y1 ) 1 - x1, and the reactor temperature, y2 ) x2; the disturbances are the feed composition, d1 ) x10, and the feed temperature, d2 ) x20; and the manipulated variables are the coolant flow rate, u1, and the reactor hold-up, u2. Following Subramanian and Georgakis,10 two distinct parameter sets were investigated: set 1, studied by Luyben,14 which features high conversion, and set 2, studied by Ray,15 which features intermediate conversion and potential steady-state output multiplicity. The model parameters relating to both sets are summarized in Table 1. Set 1 (Luyben13). The same seven designs investigated by Subramanian and Georgakis10 are treated here, featuring an increase in the operating temperature, while the conversion is maintained constant at 95%. Five alternative design approaches are considered. In the first case analyzed, it is assumed that the holdup and feed composition are constant (d1 ) 1, u2 ) 1), with a disturbance range of -0.15 e d2 e -0.083 and an allowable input range of 0.3u1nom e u1 e 1.05u1nom. The remaining cases investigate the impact of a larger disturbance range, a modified range of manipulated variables, and alternative specifications, as summarized in Table 2. In this case, the evaluation of inputs is trivial, given that, with the definition of y2, the system is linear with respect to u1 and u2 and has a unique solution. The results of the nonlinear DC analysis are reported in Table 3 and are in good agreement with those of Subramanian and Georgakis.10 As the reactor volume is decreased, the results show the expected decreasing resiliency, as indicated by increasing values of the DC, as pointed out by Luyben.14 Case 2 indicates that a feed composition disturbance in the range of 0.8 < d1 < 1.2 cannot be rejected if the
Luyben14
Ray15
1.33 0.60 25.18 19.26 91.42 11.43
0.35 1.00 20.00 0.11 0.44 0.06
q0 x10 x20 x30 qC Tref
Luyben14
Ray15
1.00 1.00 -0.12 -0.12 6.56 599.67 °R
1.00 1.00 0.00 -0.05 1.00 -
reactor level is not permitted to exceed its nominal value, which is in good agreement with previous results (Subramanian and Georgakis10). In contrast, the resiliency can be completely recovered by allowing a 25% additional hold-up (case 3). However, if the allowed range of u1 is reduced (case 4), the overall reactor volume needs to be increased significantly to provide adequate disturbance rejection. In this case, for the same composition disturbance as in cases 2 and 3, u2nom needs to be greater than 6. Finally, the nonlinear DC analysis indicates the expected increased resiliency afforded by adopting a partial control strategy (case 5), rather than attempting to control the reactor precisely at its set points; significantly smaller reactor hold-ups can be accommodated by adopting a partial control strategy. Set 2 (Ray,15 Russo and Bequette16). Following Subramanian and Georgakis,10 four different designs to maintain conversion at a constant value of 80% were investigated, again with increasingly smaller reactor hold-ups, leading to increasing values in the nominal operating temperature. Two distinct cases are examined, with the following ranges in variables: For case 1, 0.85 e d1 e 1.15, -0.083 e d2 e 0.083, 0.3 e u1/u1nom e 4, 0.25 e u2/u2nom e 1.3, y1 ) 0.8, and y2 ) y2nom. For case 2, 0.797 e y1 e 0.803, y2 ) y2nom ( 0.003, with the remaining variables in the same ranges as in case 1. It is noted that case 1 implies a set-point control strategy, whereas case 2 is an implementation of partial control, where outputs in the steady state are required to be maintained within (3% of the set-point values. The analysis results are summarized in Table 4, where the last entry is that associated with the high-temperature steady state. As in set 1, the tendency for larger reactors to have improved resiliency is preserved. For low temperatures and large reactor hold-ups, there is significantly improved resiliency when moving from a set-point to a partial control strategy. However, not even the partial control strategy can make the system work for excessively small hold-ups. It is noted that the proposed resiliency index can handle both output and input multiplicities. For example, design D of set 2 features output multiplicity for the nominal inputs (Russo and Bequette16). The partial control assumption can be interpreted as a problem featuring input multiplicity because the outputs are not fixed and there is an infinite number of inputs sets reliably keeping the process inside the desired bound. 4. Verification of Results by Simulation To verify the reliability of the predictions obtained by the DC analysis, a number of closed-loop simulations were carried out, using designs A and E in cases 4 and 5 of set 1. Recall that the DC predicts that the set-point
Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003 4509
Figure 1. Nonlinear dynamic simulation: set 1, design A-4 (set-point control).
Figure 2. Nonlinear dynamic simulation: set 1, design A-5 (partial control). Table 2. Disturbance and Manipulated Variable Ranges and Process Requirements for Set 1 case
1
remark
base case
d1 d2 u1/u1nom u2/u2nom y1 y2
1 -0.15 to -0.083 0.3-1.05 1 0.95 y2nom
2
3
4
5
feed composition disturbance 0.8-1.2 as in case 1 0.3-3 0.25-1 as in case 1 as in case 1
as in case 2 but with increased u2 range as in case 2 as in case 1 as in case 2 0.25-1.25 as in case 1 as in case 1
as in case 3 but with decreased u1 range as in case 2 as in case 1 0.3-1.25 as in case 3 as in case 1 as in case 1
as in case 4 but specifying partial control as in case 2 as in case 1 as in case 4 as in case 3 0.947-0.953 y2nom ( 0.003
design A-4 is close to its saturation point (DC ) 0.98), whereas E-4 is infeasible in that it cannot reject the disturbance (DC > 1). In contrast, from the DC analysis, one expects that both of the corresponding partial control designs, namely, A-5 and E-5, will successfully
reject the designed disturbance (both cases have DC < 1). To check these predictions, linear decentralized control was designed for all four cases, with pairings determined using the relative gain array (RGA) approach (Bristol17), with inversed-based controllers using
4510 Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003
Figure 3. Nonlinear dynamic simulation: set 1, design E-4, (set-point control).
Figure 4. Nonlinear dynamic simulation: set 1, design E-5 (partial control). Table 3. Optimization Results for Analysis of Set 1
Table 4. Optimization Results for Analysis of Set 2
DC design u1nom u2nom y1nom A B C D E F G
y2nom
DC
case 1 case 2 case 3 case 4 case 5
5.95 16.21 0.95 -0.100 0.615 3.18 9.75 0.95 -0.083 0.62 2.17 5.97 0.95 -0.067 0.63 1.32 2.35 0.95 -0.033 0.66 1.00 1.00 0.95 0.000 0.71 0.77 0.44 0.95 0.033 0.77 0.65 0.20 0.95 0.067 0.85
∞ ∞ ∞ ∞ ∞ ∞ ∞
0.84 0.84 0.84 0.84 0.84 0.84 0.84
0.98 0.99 1.00 1.04 1.09 1.17 1.29
0.32 0.40 0.74 0.88 0.97 1.07 1.19
linearized process models. Note that the RGA approach indicates that off-diagonal pairing is required, with the conversion (y1) controlled by the reactor hold-up (u2),
design
u1nom
u2nom
y1nom
y2nom
case 1
case 2
A B C D
2.30 1.64 0.80 0.25
6.48 3.52 1.52 0.80
0.8 0.8 0.8 0.8
0.094 0.132 0.189 0.236
0.80 0.99 2.55 39.1
0.37 0.51 2.21 25.57
and the temperature (y2) by the coolant flow rate (u1). The obtained transfer functions for the controllers were interpreted as a set of ordinary differential equations (ODEs) that have been added to the process dynamic model, which also accounts for possible saturation of the inputs. The obtained set of nonlinear ODEs was inte-
Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003 4511
grated numerically using the NDSolve function in Mathematica. The disturbances were simulated as step functions, with the step values obtained directly from the DC solution (the worst-case disturbance). The partial control was implemented by adaptation of the reference signals, again, obtained directly from the DC solution. The obtained results, presented in Figures 1-4, match the DC predictions precisely; note the failure of design E-4 and the successful disturbance rejection in all other cases. 5. Conclusions This paper defines a nonlinear resiliency measure, which is an extension of the previously presented disturbance cost concept. The nonlinear DC evaluation is demonstrated in a resiliency analysis of a CSTR based on a first-principles model. It is shown that the index highlights potential resiliency problems associated with design to be identified reliably. The approach also accommodates the analysis of the potential of a partial control strategy to improve process resiliency. The current results assume that the closed-form solution of the process model is readily available, which is appropriate for the simple CSTR model analyzed here, but is generally not the case. An extension of the approach to the case of more practical problems is research in progress and will be presented as continuation to this study. Acknowledgment The support of the Fund for the Promotion of Research at Technion is acknowledged. Nomenclature CA ) concentration of component A Cp ) heat capacity d ) vector of disturbances D ) disturbance domain DC ) disturbance cost F ) state equation vector G ) output equation vector ∆H ) heat of reaction k0 ) reaction rate constant nu ) number of inputs P ) matrix of transfer functions q ) dimensionless flow rate Q ) volumetric flow rate R ) gas constant RGA ) relative gain array s ) laplace variable t ) time T ) temperature u ) manipulated variable vector U ) domain of inputs UA ) specific heat-transfer coefficient V ) volume x ) state variable vector xmax ) maximum value of x
xnom ) nominal value of x x0 ) reactor feed attribute of x xC ) cooling jacket attribute of x Greek Symbols R ) temporary variable β ) heat generation/yield ratio γ ) dimensionless activation energy δ ) dimensionless volume φ ) material consumption/yield ratio F ) density τ ) dimensionless time
Literature Cited (1) Luyben, W. L.; Tyreus, B. D.; Luyben, M. L. Plantwide Process Control; McGraw-Hill: New York, 1999. (2) Seider, W. D.; Seader, J. D.; Lewin, D. R. Product and Process Design Principles: Synthesis Analysis and Evaluation; John Wiley and Sons: New York, 2003. (3) Skogestad, S.; Morari, M. Effect of Disturbance Direction on Closed-Loop Performance. Ind. Eng. Chem. Res. 1987, 27, 1848. (4) Stanley, G.; Marino-Galarraga, M.; McAvoy, T. J. Shortcut Operability Analysis. 1. The Relative Disturbance Gain. Ind. Eng. Chem. Process Des. Dev. 1985, 24, 1181. (5) Wolff, E. A. Studies on Control of Integrated Plants. Ph.D. Thesis, University of Trondheim-NTH, Trondheim, Norway, 1994. (6) Lewin, D. R. A Simple Tool for Disturbance Resiliency Diagnosis and Feedforward Control Design. Comput. Chem. Eng. 1996, 20, 13. (7) Cao, Y.; Rossiter, D.; Owens, D. Input Selection for Disturbance Rejection under Manipulated Variable Constraints. Comput. Chem. Eng. 1997, 21, S403. (8) Solovyev, B. M.; Lewin, D. R. Controllability and Resiliency Analysis for Homogeneous Azeotropic Distillation Columns. In Proceedings of ADCHEM 2000; Elsevier Science: New York, 2000. (9) Seferlis, P.; Grievink, J. Plant Design Based on Economic and Static Controllability Criteria. In Proceedings of the 5th International Conference on Foundations of Computer-Aided Process Design; American Institute of Chemical Engineers: New York, 1999; p 346. (10) Subramanian, S.; Georgakis, C. Steady-state Operability Characteristics of Reactors. Comput. Chem. Eng. 2000, 24, 1563. (11) Arbel, A.; Rinard, I. H.; Shinnar, R. Dynamics and Control of Fluidized Catalytic Crackers. 4. The Impact of Design on Partial Control. Ind. Eng. Chem. Res. 1997, 36, 747. (12) Wolfram, S. Mathematica 4: User Manuals; Wolfram Research, Inc.: Champaign, IL, 1999. (13) Loehle, C. Global Optimization V2.1s: User Manuals; Loehle Enterprises: Naperville, IL, 1998. (14) Luyben, W. L. Tradeoffs Between Design and Control in Chemical Reactor Systems. J. Process Control 1993, 3, 17. (15) Ray, W. H. New Approaches to the Dynamics of Nonlinear Systems with Implications for Process Control System Design. In Proceedings of Chemical Process Control 2; United Engineering Trustees: New York, 1982; p 245. (16) Russo, L. P.; Bequette, B. W. Impact of Process Design on the Multiplicity Behavior of a Jacketed Exothermic CSTR. AIChE J. 1995, 41, 135. (17) Bristol, E. H. On a New Measure of Interactions for Multivariable Process Control. IEEE Trans. Autom. Control 1966, AC -11, 133.
Received for review November 18, 2002 Revised manuscript received March 3, 2003 Accepted March 6, 2003 IE020926A
