Sensitivity analysis for chemical process optimization

Computers chem. Engng Vol. 20, No. 10, pp. 1177-1200, 1996 Copyright© 1996 ElsevierScienceLtd Printed in Great Britain. All fights reserved 0098-1354(95)00000-Q 0098-1354/96 $15.00+ 0.00

Pergamon

SENSITIVITY ANALYSIS F O R C H E M I C A L PROCESS OPTIMIZATION P. SEFERLIS and A. N. HRYMAK~ Department of Chemical Engineering, McMaster University, Hamilton, Ontario L8S 4L7, Canada (Received 19 September 1994; final revision received 1 June 1995)

Abstraet--A sensitivity analysis methodology based on continuation techniques is applied to chemical engineering processes. The method is based on the trace of the optimal solution path of a process model under multiple model parameter variations of finite magnitude. The first-order stationarity conditions form a system of parameterized nonlinear equations. The solution of the equation set for different values of the parameters is obtained using continuation methods (PITCON). The behavior of the optimal solution path is studied around singular points. Simple tests are used to detect singularities due to violation of the strict complementarity, linear independence constraint qualification or second-order optimality conditions. The example cases involve single and multiple unit models with equation sets of moderate to large size. The method can handle active constraint set changes and can determine the range of model parameter variations in a given direction for which the model is feasible. It can be used to determine the most sensitive unit in a plant for a given disturbance and calculate multiple solutions to the stationarity conditions for the same parameter values. Copyright © 1996 Published by Elsevier Science Ltd

INTRODUCTION The optimal solution obtained in model based optimization is subject to uncertainty associated with the estimates of the model parameters. The parameter uncertainty may be due to the random error associated with the plant measurements and the model mismatch. In addition, the effects of measured or unmeasured disturbances may cause the optimal solution to deviate from the nominal operating conditions in which the parameter estimates are valid and therefore increase the modeling error, Sensitivity analysis aims to estimate the change in the optimal solution given an estimated error (but not a statistical distribution for the error) in the model parameters. A very sensitive solution, where large changes in the optimal solution result from small changes in the parameter estimates, will not

for use in real-time optimization (Forbes and Marlin, 1994). Gazi et al. (1994) studied the structural (nonparametric) sensitivity analysis of systems described by ordinary differential equations (ODE) based on qualitative reasoning techniques and applied to the robustness assessment of controllers. Sensitivity analysis methods can be separated into two main categories: (a) local sensitivity methods that calculate the local gradients of the solution with respect to infinitesimal parameter variations; and (b) global sensitivity methods that determine the behavior of the solution under simultaneous parameter perturbations of arbitrary magnitude. Numerous methods have been developed for the study of the effects of model parameter variations to the solution trajectory of a system of differential equations (Tilden et al., 1981; Caracotsios and

generally be desirable. This implies that in a realtime optimization framework the set-points to the controllers calculated by the optimizer may have a large range of possible values even for small model

Stewart, 1985; Leis and Kramer, 1985). The sensitivity of the solution of a reaction-diffusion system described by partial differential equations has been examined by Koda et al. (1979a) and Koda and Seinfeld (1982).

mismatch (Koninckx, 1988). Model parameters that have a relatively large influence on the optimal solution would require a more accurate estimate, Sensitivity information has been utilized in order to develop process models and designs that are robust to model error (Uber and Brill, 1990) and to determine a criterion for the adequacy of a process model

The local sensitivity information is either obtained by the solution of the original differential equations coupled with the sensitivity equations (Tilden et al., 1981) or by the Green's function theorem (Hwang et al., 1978; Wacholder and Dayan, 1984). Global sensitivity analysis for systems of differential equations can be performed by direct search pattern

t To whom all correspondence should be addressed,

methods that require the solution of the original

1177

1178

P. SEFERLISand A. N. HRYMAK

system for different parameter values. Another method that requires less sampling points than the direct search method is the Fourier amplitude sensitivity test (FAST). The mathematical basis of the method was defined by Cukier et al. (1978) and has been implemented in the analysis of chemical kinetic mechanisms by Koda et al. (1979b). The FAST method transforms the multi-dimensional parameter space into a single scalar space. The sensitivity measure is the "partial variance" of the state vailable due to the variation of a given parameter and reflects the relative contribution of the given parameter to the total variance of the variable due to variations in all parameters. Stochastic sensitivity analysis methods attempt to evaluate the variance of the solution of the system given the probability distributions of the model parameters (Atherton et al., 1975). The main theoretical results for the stability and the local sensitivity of nonlinear programming problems were introduced by Fiacco (1983) with numerical procedures for the calculation of the sensitivity of the optimal solution that utilize the information available from different optimization algorithms, Shapiro (1985) studied the second-order sensitivity in nonlinear programming and investigated the asymptotic behavior of the optimal solution given the asymptotic probability distribution of the parameters. Ganesh and Biegler (1987) and Wolbert et al. (1993) proposed efficient computational approaches by working in the reduced space of the independent variables. The method provides local sensitivity information and the variations in the set of the active constraints require the solution of a quadratic program for changes in each parameter of interest, The sensitivity analysis of chemical process flowsheets was studied by Volin and Ostrovskiy (1981a,b), where they presented an adjoint method for a sequential modular process simulation environmerit. Leis et al. (1987) utilized the structure of the system and the sparsity pattern of the Jacobian matrix in order to improve the efficiency of the sensitivity calculations, Nonlinear parametric programming (NLPP) was used for the global sensitivity analysis of the firstorder optimality conditions with respect to changes in the model parameters. Jongen and Weber (1990) reviewed the main theoretical results of the parametric analysis of nonlinear programs (NLP). The parameterized, nonlinear first-order stationarity conditions of the NLP can be solved by continuation methods (Rheinboldt, 1986), which trace humerically the solution of a set of nonlinear equations parameterized by some artificial variable. Seider et

al. (1991) reviewed procedures based on continua-

tion methods that are used for the solution of NLP problems. A number of applications of the homotopy-continuation methods are available including: process design (Wayburn and Seader, 1987), heterogeneous azeotropic distillation (Kovach and Seider, 1987), continuous reactors (Seader et al., 1990), countercurrent separation processes (Salgovic et al., 1981) and multi-objective optimization (Rakowska et al., 1991). The main objective of this work is to develop a systematic methodology for sensitivity analysis in on-line process optimization. The main interest lies in the investigation of the changes in the optimal variable values and the optimal objective function value under multiple simultaneous parameter perturbations over an expected range of uncertainty. This will provide useful information about the relationship between the optimal solution and combinations of parameter perturbations. Nonlinear effects may appear in the optimal solution under multiple parameter variations that cannot be predicted by a local sensitivity analysis. Sensitivity analysis will expose the importance of the model parameters and combination of model parameters to the optimal solution which can be utilized in the parameter estimation stage by updating the estimates of those parameters whose impact on the optimal solution is the most significant. Furthermore, sensitivity analysis information can be used to determine the largest magnitude for single or multiple parameter variations for which the process becomes infeasible or loses optimality. In process optimization, with many process units, sensitivity analysis can be utilized to investigate the interaction between the different units under the influence of disturbances or parameter uncertainty. The solution of the set of stationarity conditions with the addition of the strict complementarity condition for different values of parameters will result in a sequence of Karush-Kuhn-Tucker (KKT) points (Guddat etal., 1990). Violations of either the strict complementarity condition or the feasibility requirement during the calculation of the optimal solution path will indicate an active set modification when inequality constrants or variable bounds are present. The method can determine the range of parameter perturbations in a given direction for which the active set remains unchanged. Active set changes will provide useful information about the importance of process constrants, such as product specifications and capacity constraints, in the presence of model uncertainty. Singularities in the set of equations due to violation of either the linear independence constraint qualification or the second-

Sensitivity analysis for chemical process optimization order optimality conditions may lead to multiple optima, transitions from local minima to local maxima or saddle points, or boundary points of the feasible region of the problem. Tests to detect different types of singularities during the course of the numerical solution are proposed and the structure of the solution path is examined around such points. The proposed methodology is applied to chemical processes that involve single and multiple units. Nonlinear sensitivity analysis of an optimal solution is performed for the modified Williams-Otto plant, a single multi-component distillation column and a train of distillation units, NONLINEAR PARAMETRIC PROGRAMMING

Background Nonlinear parametric programming tracks the behavior of the parameterized solution of the corresponding nonlinear program. The s~lution of the NLP is expressed in terms of the first order optimality conditions for a Karush-Kuhn-Tucker point. Kojima and Hirabayashi (1984) developed the theoretical background for the study of the continuous deformation of the optimal solution due to perturbations of one free parameter. The behavior of the optimal solution path around singular points was studied by Jongen et al. (1986). Jongen and coworkers classified the singular points into five types and investigated the characteristics of the stationary points near the singular points. Tiahrt and Poore (1990) analyzed the behavior of the optimal path near a singular point using bifurcation theory. Guddat et al. (1990) and Lundberg and Poore (1993) proposed efficient path-following algorithms that handle singularities in the solution path of the NLPP. Hirabayashi et al. (1993) studied the structure of the KKT solution set with two free parameters. The standard form of a parametric NLP is as follows: Minf(x, t)

(P1)

x

s.t. hi(x, e) = 0

ie I

I = {1 . . . . . m}

gj(x, t) ~<0 j ~ J J = {1 . . . . . 1} where x is the decision variable vector of dimension n and ~ is the parameter vector of dimension p. The Lagrangian of (P1)is given as: L(x,~,A,i.t)=f(x,e)+~lihi(x,~.)+E~igi(x,e.). ~E~

jEJ (1)

1179

A point $ that is a local minimizer for (P1) satisfies the first-order KKT stationary conditions: VxL(.f,,~,g,[)=Vxf(.~,E)+~,~.iVxhi(~,[ ) ~ + ~/~jV~gi($, t) = 0 J~J0 ~igi(~, t) = 0 j ~ J

(2a) (2b)

/zi>0 jeJo hi(2, ~) = 0 i ~ I

(2c) (2d)

gi($, t) ~ 0 j e J.

(2e)

J0 is the index set of the active inequalities at point .~ defined as: Jo = {j~J/gj(~, ~)=0}. In addition to the KKT conditions, the gradients of the active constraints at point ~ are assumed to form a set of linearly independent vectors which is an assumption known as the linear independence constraint qualification (LICQ) (Fletcher, 1987) which guarantees the uniqueness of the Lagrange multipliers at point 2. The stationarity conditions (2a) and (2b), combined with the feasibility relation (2d), form a set of parameterized nonlinear equations. An active set index strategy that eliminates inactive inequality constraints is used to reduce the size of the system. Since the Lagrange multipliers that correspond to the inactive inequalities are equal to zero, then the equation set takes the following form: -VrL(x, 2,#, ~)l

F(x, 2, g, ~) =

hi(x, t)

~ = 0 i ¢ 1, j e J0. (3)

gj(x, ~) The number of equations is equal to (n + m + r), where r<~l is the number of the active inequality constraints, including (n + m + r + p ) unknowns. Thus the system has p degrees of freedom, equal to the dimension of the parameter space. The set of the KKT points for (P1) for different values of the parameter vector is obtained by solving (3) with the addition of the strict complementarity condition (2c), that requires the Lagrange multipliers which correspond to the active inequality constraints to be positive. A point z0 = (x0, i0, ~0), that is a solution of (3) for some parameter values t0, is called a generalized critical point. When (3), LICQ and (2c) are satisfled, then the point is called a critical and KKT point. A point z which satisfies (3) and has a nonsingular Jacobian of F with respect to z (VzF) is called a regular point of F. If VzF is singular, then the point is called a singular point of F.

1180

P. SEFERLISand A. N. HRYMAK

The implicit function theorem (Fiacco, 1983) states that the solution of (3), (z0, to), can be parameterized by means of the parameter e in a neighborhood of (z0, ~o) provided that the Jacobian of system (3) with respect to z is nonsingular at the point (z0, ¢-0). More specifically, there exists a once continuously differentiable function ~0(e): such that F(9(e), e ) = 0 for all e near e0 and 9(Co)= z0. The solution set of system (3) is a p-dimensional manifold. The entire analysis of the behavior of the optimal solution path as the parameter vector t varies is based upon the following theorem (McCormick, 1976; Tiahrt and Poore, 1990). Let (z0,eo) be a solution of (3) and assume that f, h, g are twice continuously differentiable in a neighborhood of (x0, z0). Then a necessary and suffident condition that VzFis nonsingular is that each of the following three conditions hold: --Cl.

Strict complementarity condition (SC); /~j > 0 for j E J0. - - C2. {Vxhi(x,e), i ~ I O V~gj(x, e), j ~ J0} is a set of linearly independent vectors; (LICQ). - - C3. The reduced Hessian of the Lagrangian ZrV2xLZ is nonsingular at (z0, t0), where Z is a matrix whose columns form a basis for the null space of the gradients of the active constraints.

A direct extension of the previous statement is that the characteristic features of a stationary point may change only when a singular point is crossed in the optimal path trajectory (Guddat et al., 1990). Kojima and Hirabayashi (1984) and Guddat et al. (1990) introduced indices in order to describe fully the characteristics of a stationary point (local minimum, local maximum, saddle point). They used the linear index (linear co-index) as the number of the negative (positive) Lagrange multipliers gj that correspond to the active inequality constraints (jeJo) and the quadratic index (quadratic co-index)which denotes the number of the negative (positive) eigenvalues of the reduced Hessian of the Lagrangian. For instance, a local minimum has the linear and quadratic indices equal to zero (all gi are positive and the reduced Hessian is positive definite) and a local maximum has the linear and quadratic coindices equal to zero. The local characteristics of the optimal solution can be fully determined by these four indices and they are invariant under parameter variations until a singular point is reached. The type of singularity in the optimal solution path, which depends upon the violated conditions C1-C3, defines the way that the indices will be altered (Jongen et al., 1986; Poore and Tiahrt, 1987). A

detailed analysis of the solution curve around singular points is presented in the following sections. Multiple parameter perturbations

The parameter vector e may contain a large number of physical parameters of the process model. The senstivity of the optimal solution to multiple parameter perturbations may differ from the sum of the effects of individual parameter variations (Hearne, 1985). Most numerical continuation algorithms can handle problems with only one degree of freedom; i.e. one independently varying parameter. However, the equation system can be modified to accommodate the requirement for multiple parameter changes. A new scalar parameter (a) is introduced which may vary between zero and any defined target value. Additional equations are introduced that express the functional relation between the changes of the physical model parameters and the variation of the artificial free parameter a. Then the value of each model parameter at a certain point in the optimal solution curve satisfies the following relation: Ae(e, a)i = (e.ew- t,ef)i Oi(an~w)= 0 (e,~t)i

i = 1..... p

(4) where ~ t and r~,w are the reference (initial) and the current parameter values vectors, respectively. The first term corresponds to the relative change of the ith parameter at the current point. The relative rate of change is chosen for scaling purposes. The vector 0 represents a direction of perturbation in the parameter space which is chosen arbitrarily, by utilizing engineering knowledge about the most frequently occurring variations, or by using the correlation between the parameters. The product Oiot at the final point, where a is at its target value, reflects the final relative change of the ith parameter from the initial point. The error associated with the estimate of every parameter may be incorporated in the assignment of the value for 0. Therefore, larger moves are allowed for parameters that have larger error in their estimates and more conservative moves are allowed for parameters that have an accurate estimate, anew is the relative perturbation size at the current point in the optimal solution curve. Equation (4) implies that the parameter values vary linearly between the initial and final value. However, relation (4) may take the form of any nonlinear function that describes the change in the parameters (e.g. catalyst deactivation curve with time, kinetic parameter changes with temperature) as follows:

Sensitivity analysis for chemical process optimization (end,,-~!~ Ae(t, a ) = \ ~ /i-d(a)

(5)

where d(a) is a nonlinear function that relates the perturbation size with the relative change of the ith parameter. An example of nonlinear behavior for two model parameters is shown in Fig. 1.

Continuation methods for the path following of the optimal solution The solution of the parameterized equation set (3), with the addition of the relations that govern the variation of the physical parameters of the process model, (5), is obtained by continuation methods,

1181

Rheinboldt (1986) described a continuation method to determine the solution of nonlinear equations. Allgower and Georg (1990) reviewed a large variety of algorithms for continuation methods which can be classified into two types: predictor-corrector methods and piecewise linear methods that approximate the solution curve by a sequence of linear segments. A predictor-corrector continuation method is used in the present work as implemented in PITCON (Rheinboldt and Burkardt, 1983a,b). PITCON allows one parameter to be varied independently. The equation system that is to be solved and allows multiple model parameters to vary as a

slope=edO~

o . . 7 . . . . / D,,

0

1~2

(s-s~),/~,~.~

(a)

C~arget

I ~--0 m

(b)

Fig. 1. Relation between the relative changes of parameters for multiple model parameter variations (a) linear relation [equation (4); (b) nonlinear relation [equation (5)].

1182

P. SEV'ERLISand A. N. HRYMAK

function form:

of the

independent parameter a has the

In this section, the discussion involves the detection and study of the behavior of the optimal solution path around singular points. Without any loss of generality, the analysis is performed for system F, equation (3), and it is assumed that e is onedimensional. The extension of the results to system FMp is straightforward.

-VrL(x, ,~, iz, ~)

hi(x, 0

FMp(X,Adz, e, a)=

=0

gj(x, ~) Ae(e, a) i~I,j~Jo.

(6)

The predictor stage evaluates the tangent vector to the optimal solution path and then a step is taken along the tangent direction to calculate the predictor point. PITCON automatically adjusts the step length of the predictor step. A corrector stage based on Newton's method calculates the new continuation point of the optimal solution path. Since system (6) is under-determined, PITCON augments the equation set. This additional equation forces the variable that has the largest entry in magnitude in the tangent vector, calculated in the predictor step, to be fixed at its predictor point value during the Newton's iterations in the corrector step. The solution of the linear systems is performed by block elimination using LINPACK package (Dongarra et al., 1979). The system of equations that appears during the course of the continuation method includes the Jacobian of FMp, (9), which has the following form: IV20HL V r H VwF~=

V

V2xtL ' VrH I

0 0

0 i 0 e~f

(7)

where VxH =

SINGULARITIES IN THE OPTIMAL SOLUTION PATH

vxhiEl ] LVxgJ~J0J .

I is the (p x p ) identity matrix and er~f is a (p x 1) column vector with entries the reference values of the parameters. The matrix in (7) contains the full Hessian of the Lagrangian function which is evaluated numerically; a procedure that requires great computational effort. In the present work, a modifled Newton's procedure is followed that does not update the Jacobian matrix at every correction step iteration. In large-scale problems that are mildly nonlinear, the Jacobian may be updated at the continuation steps that the corrector stage fails to converge. Then the corrector stage is retried using the updated Jacobian. A smaller predictor step is tried if the corrector stage continues to fail with an updated Jacobian.

Violation of the strict complementarity (SC) condition In this case, the SC condition is violated at point (z0, ~0) causing VzFto become singular. It is assumed that conditions (C2) and (C3) hold. There are two distinct cases for violation of the SC condition. The first case involves the situation where one multiplier which corresponds to an active inequality crosses zero (SC violation). The second case occurs when one of the inactive inequalities becomes infeasible (feasibility loss). Jongen et al. (1986), Guddat et al. (1990) and Tiahrt and Poore (1990) studied extensively the behavior of the characteristic indices around a singular point that SC is violated. The check of the satisfaction of the strict complementarity and the feasibility conditions is carried out at the end of every continuation step before advancing to the calculation of the new point. In the case of SC violation, where a Lagrange multiplier that corresponds to an active inequality constraint becomes zero, thecorrespondingconstraintisremovedfrom system (6) and the path following continues on the modified set. In the case of feasibility violation, where an inactive inequality becomes infeasible for some parameter values in the optimal solution path, the violated inequality is added to the equation system and a corresponding Lagrange multiplier is introduced as an additional variable. Path following is then performed on the altered equation system. As mentioned before, the check for the violation of the SC and the feasibility condition is performed after the evaluation of the continuation point. In the case of SC loss as soon as a Lagrange multiplier of an active inequality becomes negative, the continuation direction is reversed temporarily to calculate the point where the respective Lagrange multiplier is equal to zero. Then, the corresponding constraint is removed from the equation set and the continuation algorithm proceeds on the modified set. This procedure is followed in order to calculate accurately the point in the optimal solution path at which the SC condition is violated. In the feasibility loss case, when an infeasibility is detected, the corresponding inequality constrant is added to the equation set with its associated Lagrange multiplier and

Sensitivity analysis for chemical process optimization the optimality conditions (2a) are properly modifled. The continuation algorithm performs a correctot step to adjust the point so that the Lagrange multiplier of the added constraint is equal to zero. Then, the continuation procedure begins for the altered set of equations, A simple example taken from Fiacco (1983) that illustrates the behavior of the optimal solution path around a singular point of this kind is the following: Min F(xl, X2) : (XI -- [)2 + (X 2 j¢ 1)2

identify the set of constraints that first violate any of the two conditions for the exact same parameter values. Then the corresponding constraints are added or removed from the equation set accordingly and the continuation procedure is applied in the new system. The violation of the SC condition may lead to bifurcations in the optimal solution path, when there are different possible active set selections that result in branches consisting of local minimizers for

(P2)

problem (P1).

Linear independence loss

s.t. gl(xl, x2) = xl -x2~<0

gE(X1,X 2 ) =

--X 1 --X 2~

Let (z0, eo) be a solution of F = 0 and the rank of at that point is equal to (m + r-l) while around this point the rank of V,F is equal to n + m + r, which implies that V,F does not belong to the range

0.

VxH

The objective function is the minimization of the distance from the feasible region to point (e, - 1). For e = - 2 . 0 (point A in Fig. 2), constraint g2 is

space of VzF. It can then be shown that the optimal solution path of (P1) exhibits a quadratic turning

active. At point B where e = - 1.0, constraint gl becomes active and the active set is properly modi-

point (Jongen et al., 1986; Keller, 1987). The Lagrange multipliers, that correspond to the

fled. Continuation in the dotted curve causes gl to be infeasible. In the segment B-C, both inequality constraints are active at the solution with both x~ and x2 equal to zero. At point C (~= 1.0) the Lagrange multiplier that corresponds to g2 becomes

constraints that cause the rank deficiency, tend to infinity when the LI condition is violated. As a consequence, the optimal solution path would be impossible to calculate using the described continuation method near the singular point. Guddat et al. (1990) suggested that a constraint should be

equal to zero, which results in a SC violation case. The constraint is removed from the active set and the optimal solution path is traced until point D (e= 2.0).

imposed on the magnitude of the multipliers in order to approach the singularity. The stationarity conditions are thus modified by introducing a

The proposed formulation can handle simultaneously multiple changes in the active set index,

multiplier for the gradent of the objective function.

meaning that more than o n e constraint may enter o r leave the basis. In a situation where multiple violadons of the SC and feasibility conditions occur, the last continuation step is repeated using a smaller step length in the predictor stage. The objective is to 0.6

1183

Vx L(x, 2, IZ, ~) vVxf (X, ~) + E 2iVxhi(x, E) ,~t =

+ E/*iVxgi(x' ~) =0. J, J0

i

0.5 A

D

0.4

gl

a ~

0.3 0.2 x2

0.1

/

o -oA -0.2 -2

gl,g2active

B'.,

C

'.,gl infeasible ~ -1.5

' -1

~ -0.5

~ 0

' 0.5

' 1

1.5

parameter, s Fig. 2. Optimal solution path for problem (P2).

2

(8)

1184

P. SEFERLISand A. N. HRYMAK 3

i

i

2

i

i

B( "",,

1

"'"",,

'""--,. max 0

"'"-

X

""'",.

-1

~

-2

C

~

-3

i

i

i

I

0

0.5

1

1.5

2

parameter, Fig. 3. Optimal solution path for problem (P3).

At KKT points, v is assumed to be equal to one. The imposed constraint has the following form: I~2"I'~T~'~-~.gT/u=fl 2. (9) The term f12 is a constant scalar that is equal to the sum of squares of Lagrange multipliers including v2 at the point where v is to be introduced, The system, for which the solution set is to be calculated, thus becomes:

straints. An example of LI loss is given by considering the following simple minimization problem (Guddat et

F u ( X , r , A , / z , 8) = -1

I VrL(x' 1/, 2, /t, 8) / hi(x, 8) I, g j ( x , 8)

=0

iel, jeJo.

(10)

/

Lv2 +2r2 +l~rla-321 The solution of system (10) is equivalent to the solution set of (3) with the Lagrange multipliers

givenby]q=(;ti/v)ieI, and#j=(l~j/v)jeJo(Guddat et al., 1990). As v approaches zero, then /[ and /i---* oo. If v = 0, then matrix VxH is singular (the Lagrangian becomes

~z

will the eigenvalues of the reduced Hessian (Tiahrt and Poore, 1990). A direct consequence will be that a path consisting of local minimizers will become a path consisting of local maximizers when a singular point due to LI loss is crossed. The branch of the local maximizers can be traced numerically using continuation but in the opposite direction. In practical applications LI loss may occur when the varying parameters appear in the set of nonlinear con-

~iVxhi + ~ #jVxgj = 0 J~o

with nonzero multipliers). If V~His nonsingular then v=/=0 (because if v = 0, then V~H would be singular), This simple test is used to detect the LI violation in the optimal solution path. v is normalized by dividing the stat;.onarity conditions by the numerical value of v, when it crosses zero and takes negative values. The variables x will remain unchanged, but all the Lagrange multipliers will reverse sign and so

al., 1990): Min f(x) = (x - 3) 2 s.t.

h(x) =

1 / 3 x 3 - 1 / 2 x 2 - 2 x -t- 5 -

(P3) 3.58

8 varies from one to zero. For z = 1 the minimum is at x = 3 (Fig. 3). The Lagrange multiplier for h approaches infinity at x = 2 and 8 = 0.48, since the gradient of the constraint becomes zero at that point. This is a case that the LI condition does not hold, the optimal solution path exhibits a quadratic turning point and the local minima branch turns to a local maxima branch. At x = - 1, the gradient of the constraint becomes zero again and there is another turning point in the solution path with the characteristics of the optimal solution changing to those of a local minimum. Figure 3 indicates the existence of multiple minima for the same parameter value which are calculated by the continuation method.

Specialcase of LI loss A special case of LI loss will be investigated in this section. In many practical applications of NLPP,

Sensitivity analysis for chemical process optimization

1185

where inequality constraints are involved, it is possible that the total number of active constraints (equality plus binding inequality constraints) exceeds the number of process variables. Such a situation may occur when the solution of a fully determined system hits another process constraint due to variation in the model parameters. More specifically, the case where m + • = n + 1 will be examined based on the work by Jongen et al. (1986) and Guddat et al. (1990). The analysis requires that the number of active constraints is n + 1 with r~>2 and V,H~R(V~H) (where R denotes the range space of V~H) which implies that {V~x.,)h~, i ~ I, V~x.,)gi,J • J0} is a set of linearly independent vectors (V(x.~)H has rank equal to n + 1). Since LI does not hold, then the set of

these branches will necessarily be comprised of local minimizers. Jongen etal. (1986)have proved that if MFCQ hold at point (z0, to), then there is a solution branch consisting of local minimizers for t > to. In contrast, if MFCQ does not hold at (z0, to), then this point is a boundary point for the KKT set. This implies that there are no local minimizers for e > to, thus defining the boundary of the feasible region for problem (P1). The boundary point will determine the range of parameter changes that the system reaches the limits of the feasibility region. The Lagrange multipliers of the active inequalities will tend to infinity and the numerical procedure described in the previous section will be followed.

Lagrange multipliers does not have a unique set of values, A weaker constraint qualification than the LICQ was introduced by Mangasarian and Fromowitz (1967) denoted as Mangasarian-Fromowitz constraint qualification (MFCQ) which requires that:

In this section it is assumed that SC and LI conditions hold and exactly one eigenvalue of the reduced Hessian of the Lagrangian vanishes at point (z0, to). The optimal solution path exhibits a quadratic turning point at (z0, to), with one eigenvalue of the reduced Hessian changing sign (Jongen et al.,

(i) V.h~ i • ! is a set of linearly independent vectors; and (ii) =1~•~" such that Vxhi~=0 i • l , Vxgj~<0

1986). Therefore, a solution branch that consists of a sequence of local minimizers will turn to a branch of saddle points. This is equivalent to loss of the second-order optimality condition. The active set remains unchanged and so does the linear index when the optimal solution path crosses such a point. As an example consider the following onedimensional problem:

j eJo. The Lagrange multipliers are bounded and unique up to a common multiple when the MFCQ is satisfled (Kojima and Hirabayashi, 1984). The system of equations (3) is feasible only at one point (z0, to), when n + 1 active constraints are present. In the neighborhood of the point (z0, to), the stationarity conditions are satisfied if one of the active inequalities is removed. The set of generalized critical points will consist of the union of solution branches that correspond to the different problems in which one active inequality has been removed from the basis (Fig. 4). However, not all of

Singularity of the •educed Hessian

Minf(x) = 1/5x 5 - 1/4~r*- 1 / 3 t ( x - 2 )

3. (P4)

Figure 5 shows the optimal solution path as the parameter ~ varies fron one to zero. At ~ = 0, the second derivative of the objective becomes zero and the second-order optimality condition is violated. The optimal solution has a turning point and the second derivative changes sign (negative) which results in a branch consisting of local maxima.

MFCQ holds

MFCQ violation

X

X

g

E

Fig. 4. Typical behavior of the optimal solution path around a singular point due to special case of LI violation (sofid curves: local minimizers; dashed curves: generalized critical points).

1186

P. SE~RLIS and A. N. HRVMAK 1.5

,

,

~--A

0.5'

0 x

B~

-0.5

-1

-1.5

max

................................ C

-2 -0.2

i 0.2

0

t 0.4

i 0.6

, 0.8

i 1

1.2

parameter, s Fig. 5. Optimal solution path for problem (P4).

The presence of this kind of singular point can be detected by inspection of the sign of the determinant of the reduced Hessian (det{ZTV2LZ}). Matrix Z has size n x ( n - m - r ) and columns that form a basis for the null space of VxH. Z can be easily calculated by performing a QR decomposition of VxH. However, for a large-scale problem, a QR decomposition is computationally intensive. An alternative way to evaluate Z would be to find a matrix that is orthogonal to VxH. If

r I ] Z = [ - (V~pH) -~ (V~pH)J '

the requirement that VxHZ=0 is satisfied for the partition of the variables into independent (Xi,d)and dependent (Xdep). The reduced Hessian is then calculated numerically by pertubing simultaneously all the variables along the directions defined by the columns of Z (Wolbert et al., 1993).

SP product Distillation

SA

SB

SE TUXW

__

[

[

Heat exchanger

SR

Reactor

[

/

FUXW

SX

Decanter

[

Ss

I FURW TURW

SG SD

SL, ST

Recycle stream Reactions a+b~c

b+c--~p+e p+c~g

Fig. 6. Williams-Otto plant flowslieet.

Sensitivity analysis for chemical process optimization

1187

Table 1. Model parameters for the Williams-Otto plant Reference value

Model parameter Reactor cooling jacket heat transfer coefficient,kJ/(s K) Heat exchanger heat transfer coefficient, kJ/(s K) Arrhenius gain in reaction 1, (s. wt. fract.) -1 Arrhenius gain in reaction 2, (s. wt. fract.) -1 Arrhenius gain in reaction 3, (s. wt. fract.) ~ Arrhenius exponent in reaction 1, K Arrhenius exponent in reaction 2, K Arrhenius exponent in reaction 3, K Heat of reaction 1, kJ/kg Heat of reaction 2, kJ/kg Heat of reaction 3, kJ/kg Reactor's cooling water inlet temperature, °C Heat exchanger's cooling water inlet temperature, °C Flow rate of stream SA (reactant a), kg/s Temperature of stream SA, °C Temperature of stream SB, °C Temperature of stream SL, °C Mass holdup of the CSTR reactor, kg

VxL(x + ~/e~Z)- VxL(x)

2.6 × 10-5 4.457 x 10-4 2.7778 4.1667 5.5556 6,666.7 8,333.3 11,111.1 0.523 0.209 0.599 15.7 15.7 1.827 21.3 21.3 38.0 2,104.7

0 in first 0 in second eigenvector eigenvector directon direction -0.1885 - 0.0674 0.0015 0.0004 -0.0003 -0.0066 -0.0005 0.0004 0.2868 0.3431 0.0592 0.0441 0.5807 0.2144 0.0938 0.1880 0.5719 0.0003

- 0.2584 0.0108 0.0027 -0.0002 0.0003 -0.0130 0.0013 - 0.0016 0.3869 0.4603 0.0742 0.0635 -0.0916 0.2045 0.1348 0.2524 0.6555 0.0006

0 in scenario2 direction -----0.5773 -0.5773 0.5773 -----------

(11)

multiple parameter changes in the optimal solution are investigated. A three-reaction model is used in

where 7 is the perturbation size and ei is the ith vector of the canonical basis. This requires only n - m - r evaluations of the constraint gradients which allows significant savings in the computation of the reduced Hessian. In cases where the n u m b e r of active constraints is equal to the number of variables (no degrees of freedom), the reduced Hessian cannot be defined since the null space of

the reactor with 'p' being the product and 'g' the undesired by-product. Pure components 'a' and 'b' (streams S A and SB) and a recycle stream (SL) are fed into the reactor, which is assumed to be perfectly mixed. All three reactions are exothermic with Arrhenius temperature-dependent kinetic parameters. There is a heat removal system from the reactor, and the reactor's effluent is further cooled at a target temperature in the heat exchanger. A

(er,ZrVEL)-

9,

V~H has dimension equal to zero. P I T C O N has the ability to trace solution curves that exhibit turning points (limit points) in some of the variables. A limit point for a variable occurs when the tangent element that correspond to this variable changes sign. A simple inspection of the tangent elements in the predictor stage will detect such a situation. The eigenvalues of the reduced Hessian do not change sign unless a singular point is encountered in the optimal solution path. H e n c e the

logarithmic m e a n temperature difference is used for the heat exchanger's energy balance. There is no constraint imposed on the temperature of the heat exchanger's outlet stream. The decanter removes all the amount of c o m p o n e n t 'g' from the system. The distillation column column is modeled as a simple split with the overhead stream consisted of pure product 'p'. The objective function to be maximized is: nstr

reduced Hessian will be examined for positive deftniteness only when the optimal solution crosses singular points such as limit points in the free parameter a, or points where either the SC or LI con-

%return = ~ (pricei FLOWi) i=t -A(FUXW+

ditions are violated. PARAMETRICSENSITIVITYANALYSISFOR CHEMICAL PROCESSES

Williams-Otto plant example (scenario 1) The optimal solution path of the Williams-Otto plant is calculated for model parameter variations along specified directions. The plant shown in Fig. 6 consists of a continuous stirred tank reactor, a heat exchanger, a decanter and a distillation column (Williams and Otto, 1960). The combined effects of

FURW)

- B ( T U R W - Tt~r)2.

(12)

denotes the flowrate of the ith stream, F U X W and F U R W are the heat exchanger and the reactor cooling water flowrates, T U R W is the outlet ternperature of the reactor's cooling water flowrate, Tt~ is a target temperature for the cooling water, pricei is the price per kg of each stream (negative for feed streams, positive for product streams, zero for intermediate streams), A is the cost factor for the cooling water requirements and B is the cost factor for the temperature difference in the reactor's coolant. Flow i

1188

P. SEFERLIS and A. N. HRYMAK

The values of the decision variables and the model parameters are scaled so that their values lie between zero and one. The optimal solution for the reference parameter values are given in Table 1 has three degrees of freedom. A local sensitivity analysis has been performed for the 18 parameters. The resulting sensitivity matrix X (with elements the

partial derivatives of the optimal variable values with respect to each parameter), reveals that the most sensitive variable is the heat exchanger cooling water flowrate (FUXW). The heat exchanger cooling water inlet temperature and the target temperature of the outlet process stream have the largest effect on FUXW. Another group of variables that

% return 15

i

lO(~ ~

i

i

J

X

-

X 5

X

0 o

X

0

X -s

0 X

-10

0

-20

I 0.02

] 0.04

X

I 0.06

I 0.08

0.1

perturbation size, ct

(a)

Heat exchanger cooling water flow rate kg/s 3o

O

O

'

'

X 2s

X

O O O O

20

X

X

X 0

0 X ~° ( x

0

X

O

O

is

X

O

X

O

0

'

0

0

~

I 0.02

I 0,04

I 0.06

perturbation size, a

I 0.08

0.1

(b)

Fig. 7(a-b). Optimal solution path for (a) objective function value; (b) heat exchanger cooling water flowrate.

Sensitivity analysis for chemical process optimization

1189

Reactor temperature deg C 80

7s _

i

i

i

i

OX~×

70

X ×

0 65

0

× 0

×

6O

0

×

0

50

I

I

0.02

I

0.04

I

0.06

0.08

0.1

perturbation size, a

(c)

Flow rate of reactant 'b', kg/s 3.5

3

×

2.5

× X

0 0

X 2

0 X 1.5

0

I

0,5

X

0

I 0.02

I 0.04

X

I 0.06

perturbation size, a

I 0,08

0,1

(d)

Fig. 7(c-d). (c) Reactor's temperature; and (d) flowrate of reactant 'b' (~: perturbation along first eigenvector direction; x : perturbation along the second eigenvector direction).

have the second largest, in absolute values, entries in the sensitivity matrix are the process stream flowrates (streams SL, SR, SX, SE, SS and ST). The definition of perturbation in the parameter space that causes the largest variability in the optimal variable value space is determined by the eigenvector direction that corresponds to the largest CA~ ZO:lO-a

eigenvalue of the matrix xTx (Hearne, 1985). The eigenvectors and eigenvalues are calculated using a singular value decomposition of matrix XTX. The relative importance of the model parameters to the optimal solution may be assessed by the contrihution of each parameter to the eigenvector direction. The eigenvector that corresponds to the second

1190

P. SEFERLIS and A. N. HRYMAK

Temperature, deg C 80

i

J

i

I

,0

mperature 60

A,t~ 50

..................................... ./:"

40

30

20 -0.25

/ /

A

:

Cooling water t e r n p ~ . ~ . s . . . . ~

~

I

I

I

i

-0.2

-0.t5

-0.1

-0.05

0

perturbation size, a Pig. 8. Optimal solution path for changes in the activation energies of the reactions (solid curves: local minimizers; dotted curves: saddle points),

largest (in magnitude) eigenvalue of the sensitivity matrix defines the second most important direction of variability. The optimal solution path is traced for changes along the first and second eigenvector directions which will cause the largest changes in the optimal variable values. It should be noted that the eigenvectors are calculated using local sensitivity information and thus represent the major directions of variation near the reference (base) optimal point, The eigenvectors are normalized and their entries that define the values of 0 for each parameter are shown in Table 1. There are 150 unknowns (67 process variables, 64 Lagrange multipliers, 18 model parameters and the perturbation size a) and 149 equations (67 gradients of the Lagrangian function, 64 equality constraints and 18 relations for the parameter moves). The optimal solution paths for the profit function and some of the most important variables are shown in Fig. 7a--d. The cooling water flow rate rapidly increases as the perturbation size increased until it

reaches a maximum value at a=0.03 and 0.054 for the first and second eigenvector direction respectively; then decreases due to quenching of the reaction. The sensitivity analysis using continuation reveals that a 1.7% increase in the heat exchanger's cooling water inlet temperature and in the target temperature for stream SL when combined with smaller variations in the values of the remaining parameters may cause the optimal value for FUXW to triple. The profit decreases for perturbations in both eigenvector directions, but in the first case the slope of the solution path is steeper. The marginal parameter values for which the plant remains profitable can be determined from Fig. 7a. This may be more beneficial in a case where the perturbed paremeters are cost terms or market values in the objective function. The optimal solution path obtained by PITCON is compared with the optimal solution calculated by MINOS 5.3 (Murtagh and Saunders, 1975) for the parameter values at every continuation point. Both

Table 2. Variable bounds for the DIB Upper bound Bottoms product flowrate, B (Mmol/d) Overhead product flowrate, D (Mmol/d) Reflux ratio, RF Reboiler's heat duty, QB (MJ/d) Condenser's heat duty, QD (MJ/d) iC4 molar fraction in bottoms product, xB(iC4) nC4 molar fraction in overhead product, XDnC4)

3.0 5.2 12.0 800.0 -650.0 0.040 0.035

Lower bound 1.5 3.0 8.0 700.0 -800.0 0.0 0.0

Sensitivity analysis for chemical process optimization

1191

MINOS and the continuation method give identical solutions for same parameter values. Each case in MINOS uses the previous basis and variable values as a starting point. The main advantage of the

arbitrarily and this may result in failure of the optimizer to converge.

continuation method over the case study approach for sensitivity calculations is that it systematizes the procedure by adaptively controlling the step size, so that convergence is facilitated. Furthermore, the continuation approach identifies the parameter values at the points where the active set varies or a singularity occurs. In the case study approach to sensitivity analysis, the sampling points are chosen

The optimal solution path of the WiUiams-Otto plant is calculated for changes in the activation energy of the three reactions. It is assumed that there is a negative bias for the activation energy of the first and third reaction and a positive bias for the activation energy of the second reaction. The values of the 0 coefficient, shown in Table 1, are equal for all three activation energies which implies that the

Profit $/d 5800

~

f

~ -r"

5600

~ ~3,,~

5400

~',~

',

'

I

~1 ~"

i

I

~,

J

,

,

I

I

r~

I

I ~

i

5000

~~

~

~

, I

~ r ~"

I

I""

I

I

,

t

~

,,

I

"~

~

~,~

I

~ ~r

"~

I

"~

~ ~<

.~

"--.

--.

-3.5

~ nC4 feed flow rate Mmol/d

13

1.5

.......~....

..........

i

....

i......

~

....... F

10 9

r i

i

~

8

......... ~"..

,.k.

.......

i

.........

i ..... i........ .............. '~....... i

~..

',

................ i ......... i ...... L

..i ......... c~

..----I

iCA feed flow rate Mrnol/d (a)

"........... -

. . . . . .

..... ."i ..... 11

3.4 3.2

_~

i 12

I ........ i a ....... L

'"-...

...::;.<::

....

1

i ........... !

........ ',

....... J : ....... '-.. i ........

..... ~:'

[

t

~

~"~

Reflux ratio

,

~

'~

- ~

I

"~ ~-

~t--..

l

~

5200

Williams-Ottoplant example (scenario2)

i

"" ...... i i ........ ', ".,.., "' i ......... i ...... !

'i

'::>~.:

...........

i

3 ~

4.2 •

""-,

.:'-~:.""

3.4 1.5

nCA feed flow rate Mmol/d

"-

4

-'-

3.2

iC4 feed flow rate Mmol/d

(b)

Fig. 9(a-b). Optimal solution path for changes in the feed composition of the D I B (a) objective

function value; (b) reflux ratio•

1192

P. SF_J~RLIS and A. N. HRYMAK

R e b o i l e r h e a t duty MJ/d

620

JfT

800

_-~,

1-~I

:

~ ~'i~

,,o

~-

, ~-~--

I+t

I

~

i

~--

,"

,it

'I

: ~--.. : 'I "-:I ~

J

~ i

"-...

~'. I ~....

/_.+ ~

760

"-.. "~m...

-2

I

'~- ,~.

'

I

~.~

I

=

~-...

I

I~

t

"1 "~.

'

"

740 I

~-~.~

"" -..

~

~

~"-~

I

~"-~.~ ~

n C A f e e d flow rate Mmol/d

1.5

n C 4 m o l e fraction in o v e r h e a d

" " - ..

3.2

~"1

~

"~

4.2

""

iC4 f e e d flow rate M m o l / d

- . _ ~ m .I .

(c)

~'~.

0.036 ~ u

I ~

I

I

0.034

.~

.__t

i

0.032

.

_.~'~

~

0.03

~

0.028

0.026

~

i

,

I

r~

~

'

I

~--.

~ "

2

Fig. 9 ( c - d ) .

~

i

i

i

~ -J~,,~

~

~~

"-

"~.

~

~<~

"-.

I "-.

~

"4...

~i'~

~,'~---. "~,~

~'-.

-~

~ ~-.

"~

"-

"~

3.2

I

,

'

'

~..

~~

~'--

"-~"~ - 3.6

1.5

I

"-~

~

~'--

"xl

I

I

~ v

'

"-

`+

n~

-~,~

~-~m,~

n C 4 feed flow rate Mmol/d

~.l-~

~,~ ~

l

,

4.2

3.8

i C 4 feed flow rate M m o l / d

(c) R e b o i l e r ' s h e a t d u t y ; a n d ( d ) n C 4 m o l a r f r a c t i o n in o v e r h e a d local m i n i m i z e r s , d a s h e d curves: local m a x i m i z e r s ) .

(d)

p r o d u c t (solid curves:

Table 3. Optimal solution path for DIB iC4 feedflow rate (Mmol/d)

nC4 feedflow rate (Mmol/d)

Continuation points

A B C D E F G

Min Min Min ]din Min

4.09594 3.6687 3.4816 3.3133 3.3133 3.2015 3.2015

1.81285 2.2401 2.4272 2.5955 2.5955 2.7073 2.7073

21 24 45 4 74 109

H I J K

Max Max "Max Max

3.2015 3.3351 3.3351 4.1088

2.7073 2.5737 2.57373 1.8000

183 9 3 31

Singular point type Starting point Upper bound violation for XD (nC4) = 0.035 Upper bound violation for Qe = 800 MJ/d Upper bound violation for reflux RF= 12.0 8C violation for QB Id Upper bound violation for xa(nC,) = 0.040 Linear independence loss Transition min ---, max SC violation for xD(nC4) Upper bound violation for Qs=800, MJ/d SC violation for RF Target point reached

Active set upper level {xD(nC4)}

{xv(nC4), QB} {xD[nC4), Qe, RF} {xn(nC4), RF} {xD(nC4), RF, xe(iC4)}

{xv(nC4), RF, xB(iC4)}

{xD(nC4), RF, xe(iC4)} {RF, xa(iC4)} {RF, xe(iC4), Qe} {ze(iC~), Qa}

Sensitivity analysis for chemical process optimization

1193

C3product l

C3-C4Splitter

I

im

I

~

Debutanizer

I

Feed =

i-C4 product m. Deisobutanizer

L~

l

Cs, heavier

/

= n-C4 product

Fig. 10. Schematic of the distillation column train.

absolute relative change for every parameter is the same at every point in the optimal solution curve, The optimal solution curve for the reactor and cooling water in the heat exchanger temperatures are shown in Fig. 8. As the kinetic parameters change, the reactor's optimal temperature decreases in order to favor the desired second reaction that produces the valuable component 'p'. For the perturbation size of -0.232 (13.4% relative change in the kinetic parameters) the temperature difference at the two sides of the countercurrent heat exchanger becomes equal (40°C), which causes the logarithmic mean temperature difference model for the heat exchanger to fail. In order to avoid such

numerical instabilities, an approximation is used for the logarithmic mean temperature as suggested by Prasad (1988). The new modeling alternative allows the continuation method to trace the curve around the turning, point at a = - 0 . 2 3 2 (AE1 =5773, AE2 = 9450, AE3 = 9622 K). A calculation of the (3 x 3) reduced Hessian reveals that one of the eigenvalues vanishes at this point. This results in a transition from a branch of local maximizers to a branch of saddle points that satisfy the parameterized equation system (6) (see dashed curves in Fig. 8). The direction of the predictor step in the continuation procedure is reversed and the solution curve is tracked until the temperature of the reactor

Table 4. Variable bounds for the column train Upper bound Debutanizer (DB) Bottoms product flowrate (Mmol/d) Reflux ratio Reflux to feed ratio Reboiler's heat duty (Mild) Condenser's heat duty (MJ/d) iC4 molar fraction in bottoms product nC4 molar fraction in bottoms product

Cr-C, Splitter (SP) Overhead product flowrate (Mmol/d) Reflux ratio Rebofler's heat duty (Mild) Condenser's heat duty (Mild) iC4 molar fraction in overhead product C3 molar fraction in bottoms product

Deisobutanizer (DIB) Bottoms product flowrate (Mmol/d) Reflux ratio Reboiler's heat duty (Mild)

Condenser's heat duty (M_J/d) nC4 molar fraction in overhead product iC4 molar fraction in bottoms product iCs + nCs + C6 molar fraction in bottoms product

Lower bound

1.50 1.050 1.20 90.0 - 130.0 0.011 0.022

0.90 0.90 0.80 74.0 - 170.0 0.0 0.0

2.00 8.0 260.0 - 180.0 0.015 0.015

1.00 11.0 210.0 -220.0 0.0 0.0

2.30 10.0 840.0 - 630.0 0.027 0.032 0.170

1.70 8.0 680.0 - 790.0 0.0 0.0 0.0

1194

P. SEFERLIS and A. N. HRYMAK Table 5a. Optimal solution path for the column train

iC4 feedflow nC4 feedflow DOF rate (Mmol/d) rate (Mmol/d) A B C D E F G H I J K L M N

2 2 1 0 1 2 1 0 1 0 1 0 -0

3.9750 4.0634 4.0713 4.0715 4.1284 4.1744 4.1984 4.1984 4.2234 4.2234 4.2987 4.2987 4.2987 4.2987

1.7559 1.6675 1.6594 1.6594 1.6025 1.5565 1.5325 1.5325 1.5075 1.4322 1.4322 1.4322 1.4322 1.4322

Singular point type

Active set

Starting point {1, 2, 3, 4} Lower bound violation for bottoms product flowrate in DIB (1.7 mMol/d) {1, 2, 3, 4, 5} Upper bound violation for purity level in DIB bottoms product (xB= 0.032) {1, 2, 3, 4, 5, 6} SC violation for bottoms product flowrate upper bound in DB {2, 3, 4, 5, 6} SC violation for reflux ratio upper bound in DB {3, 4, 5, 6} Lower bound violation for reboiler's heat duty in DB (74.0 MJ/d) {3, 4, 5, 6, 7l Lower bound violation for reflux to feed ratio in DB (0.80) {3, 4, 5, 6, 7, 8} SC violation for reboiler's heat duty lower bound in DB {3, 4, 5, 6, 8} Upper bound violation for reboiler's heat duty in DIB (840.0 MJ/d) {3, 4, 5, 6, 8, 9} SC violation for reflux ratio upper bound in DIB {3, 5, 6, 8, 9} Upper bound violation for pentane and heavier in DIB bottoms product (0.17) {3, 5, 6, 8, 9, 10} Upper bound violation for iC4molar fraction in SP overhead product (0.015) {3, 5, 6, 8, 9, 10, 11} Specialcase of LI loss SC violation for reflux to feed ratio lower gound in DB{3, 5, 6, 9, 10, 11} Upper bound violation for reboiler's heat duty in DB (90.0 mild) {3, 5, 6, 9, 10, 11, 12}

Table 5b. Constraint numbers for Table 5a Constraint No.

Description

Numerical value

1 2 3 4 5 6 7 8 9 10 11 12

Upper bound on bottoms product flowrate in DB Upper bound on reflux ratio in DB Upper bound on reboiler's heat duty in SP Upper bound on reflux ratio in DIB Lower bound on bottoms product flowrate in DIB Upper bound on iC4 molar fraction in DIB's bottom product Lower bound on reboiler's heat duty in DB Lower bound on reflux to feed ratio in DB Upper bound on reboiler's heat duty in DIB Upper bound on molar fraction of pentanes and heavier in DIB's bottoms product Upper bound on iC~molar fraction in SP's overhead product Upper bound on reboiler's heat duty in DB

1.50 Mmol/d 1.05 260.0 Mj/d 10.0 1.70 Mmoi/d 0.032 74.0 Mild 0.80 840.0 Mild 0.17 0.015 90.0 Mild

a n d t h e h e a t e x c h a n g e r cooling w a t e r b e c o m e a l m o s t identical which c o r r e s p o n d to a n infeasible situation. T h e sensitivity analysis reveals some significant aspects of t h e solution curve of t h e stationarity conditions. T h e m a g n i t u d e of v a r i a t i o n of the activation energies in t h e specified direction are d e t e r m i n e d for w h i c h the system lost its optimality due to t h e violation of t h e s e c o n d - o r d e r optimality condition. Multiple solutions of the stationarity con-

p r o p e r t i e s of t h e c o m p o n e n t s a n d the p h a s e equilibria d a t a are c o m p u t e d b y regressed e q u a t i o n s (Bailey et al., 1993). A m o d e l b a s e d o n t h e H u g h m a r k correlations (Lockett, 1988) is used to estim a t e the M u r p h r e e tray efficiencies. T h e size of t h e m o d e l is 374 equality c o n s t r a i n t s with 372 process variables. T h e r e f o r e , t h e p a r a m e t e r i z e d set of optimality conditions includes 746 equations, T h e o b j e c t i v e function t h a t is m i n i m i z e d is the total cost for t h e s e p a r a t i o n , which includes the

ditions for t h e same set of p a r a m e t e r values were f o u n d , b u t n o t all of t h e m c o r r e s p o n d to local optima,

utility costs a n d t h e differential p r o d u c t cost for the iC4 a n d nC4 lost in the b o t t o m s a n d o v e r h e a d product streams, respectively.

Deisobutanizer (DIB) example

COStDIB = 2.4 X QB + 0.4 x ( - Qo) + 15,000

This e x a m p l e involves the study of t h e o p t i m a l solution p a t h u n d e r feed c o m p o s i t i o n v a r i a t i o n s for a single m u l t i - c o m p o n e n t distillation column. T h e efficiency o f t h e m e t h o d o l o g y to h a n d l e active set c h a n g e s in t h e o p t i m a l solution p a t h of t h e c o l u m n is e x a m i n e d . A c t i v e set variations alter drastically the b e h a v i o r o f t h e process variables to p a r a m e t e r changes. T h e r a n g e of feed c o m p o s i t i o n v a r i a t i o n for which n o feasible o p t i m a l solution exist is determined, A tray-by-tray m o d e l b a s e d o n t h e M E S H e q u a tions ( H o l l a n d , 1981) is used to m o d e l the fourc o m p o n e n t D I B column. T h e t h e r m o d y n a m i c

x xB(iC4) x B + 15,000 x xt)(nC4) x D. (21) T h e base o p t i m a l solution o b t a i n e d for iC4 = 4.09594 M m o l / d a n d nC4 = 1.81285 M m o l / d with M I N O S 5.3 has two degrees of f r e e d o m a n d n o n e of the variables are at any of t h e i r u p p e r or lower b o u n d s . T h e o p t i m a l solution p a t h is traced for a d i s t u r b a n c e in t h e feed composition. T h e feed flowrate of iC4 (light key) decreases linearly f r o m its r e f e r e n c e value to a final value. Simultaneously, t h e feed flowrate for nC4 is c h a n g i n g in such a way t h a t t h e total feed to t h e c o l u m n is k e p t constant. T h e

Sensitivity analysis for chemical process optimization upper and lower bounds imposed to some of the variable values are shown in Table 2. Figures 9a--d show the optimal solution path for the objective function value, the reflux ratio, the reboiler heat duty and the iC4 molar fraction in bottoms product and Table 3 provides information about the type of singularity encountered in the optimal solution path. Starting from the base case (point A) with the feed composition varying as described before, the reflux ratio and the reboiler's heat duty are increasing in order to anticipate for the

increased losses of nC4in the overhead product that are penalized in the objective function. At point B of the optimal solution path (relative change in feed flowrates, iC~ = - 10.4% and nC4)= +23.6%), the upper purity level of the overhead product (3.5% mole in nC4) is reached and one degree of freedom is removed from the optimal solution. Beyond point B, the reflux ratio and the reboiler heat duty are increasing at a higher rate, because of the lost degree of freedom, in order to anticipate for the continuing change in the feed composition. At point

...................!............

Profit $/d 5

",...

x 10

2.35

1195

.....i.......

..... i................

i

C ~

......... i . . . . .

" ....

~:

2.25

2.2

........iiiii>.,
nC4 feed flow rate Mmol/d

1.4

3.9

DB bottoms product flow rate Mmol/d

1.5

1.25 1.8

................. ~

"..........................

"...........

i

iC4 feed flow rate Mmol/d

(a)

................i ..............

.......c;o ..................................

.......................

".......................................... ..................................... ::"=:::::............ ................................ "",

1.6 nC,4 feed flow rate Mmol/d

".............

.......... 1.4

" ...........

"......................

3.9

i

"......... i 4.2

4.3

"" iC4 feed flow rate Mmol/d

(b)

Fig. 11(a-b). Optimal solution path for changes in the feed composition in the column train (a) profit; (b) DB bottoms product flowrate.

1196

P. SEV'ERLXSand A. N. HRYMAK ............ i"'-.,

OB reflux to feed flow rate ratio .........

o9

.-'- "

. ............

......

........................... •.......................... E?

0.86

-..

~

............i...........

0.88

'"

i ................ ................ ~ ............. "...........

...... elt~ ....................

......................i

"i........... ...........

O,o, o, 0.78

.................. ~:";...... ..,

"........ ".,,.

.,-

...................

nC4 feed flow rate Mmol/d

"................ Qi-1.............. '-.,.. ~".-

............-"',, ,.

"..................

1.4

'...............

3.9

4.2

4.3

iG4 feed flow rate Mmol/d (c)

.......... •"?",.....,. ............. "i ........

DIB reflux ratio

10.02 10 9.95

"'"-.,.,.

. . . . . . . . . . . . .

...........~.....

i .....

........

i .......... ....

.......

"....

i........ii.. i ....... i..

9.9 9.85 9.8 1.8 4.3

nC4 feed flow rate Mmol/d Fig. l l ( c - d ) .

(c)

DB

1.4

3.9

iC4 feed flow rate Mmol/d (d)

reflux to f e e d f l o w r a t e ratio; (d) D I B reflux ratio.

C (iC4= - 15.0% and nC4= +33.9%), the reboiler reaches its highest heat duty level (800 MJ/d)resulting in a system with no degrees of freedom. The optimal solution at this point is simply the solution of the fully determined set of the nonlinear equatiffs for different values of the feed composition, As the path following procedure progresses, the increasing reflux ratio violates its maximum allowable level of 12.0 at point D in the optimal solution path (iC4= - 19.1% and n C 4 = +43.2%). The value of the reflux ratio is fixed at its upper level, which results in an overconstralned problem with three

fixed variables. Therefore, the system can be satisfled only for a unique feed composition. This is a situation where the path encounters a point at which a special case of LI loss occurs. The Lagrange multipliers associated with the active inequalities are changing dramatically for infinitesimal perturbations in the feed stream. According to the analysis described previously, it is suggested that the stationarity condition system is augmented as in (10), so that the multipliers are bounded and possible LI violation is detected. At point E, which almost coincides with point D, the Lagrange multiplier associated with the

Sensitivity analysis for chemical process optimization

1197

iC4 mole fraction in SP overhead product .......................~,............. .................i

Ii

.........

..,. . . . . . . . . . . . . . . . . •............ ....................... ii" --..... "

001,

............ i

..............i...........

i ...................i............

q~mm

b~

"-......

.............i

.............L

0,0135

0.0125 0.012 0.0115 1.8

............................

D

... 4.3 4.1

nCA feed flow rate Mmol/d

1.4

3.9

iC4 feed flow rate Mmol/d

(e)

Fig. ll(e) iC4 Molar fraction in SP overhead product. active upper bound on the reboiler's heat duty crosses zero, violating the SC condition, which implies that MFCQ is satisfied at this point. The constraint that forces the reboiler's heat duty to be

meaning or any practical application, but the path is tracked to show the optimal solution type change.

at its upper bound is removed from the basis thus leaving two fixed variables in the optimal solution. At point F (iC4 = - 21.8% and nC4 = + 49.3%)

This example involves the simultaneous optimiza-

the upper bound for the purity level in the bottoms product (4.0% mole of iC4)is reached. When the corresponding constraint, that fixes the bottoms purity level at its upper bound, is added to the basis all the Lagrange multipliers of the binding inequalities increase rapidly for infinitesimal changes in the feed composition. The augmented system (10) is used again and allows the continuation method to calculate the solution very close to point G where v = 0 and the LI condition is violated. This point also defines the borders of the feasible region for the column for the given parameter variation. No optimal solution can be obtaned for feed composition beyond this point. At point G, the direction that the

Column train example

tion of three multi-component stagewise distillation columns in series (Fig. 10). The behavior of the optimal solution path is studied for disturbances in the feed composition. It is of great interest to investigate the interaction between the different units under the influence of disturbances. The feed to the column train is a 12-component mixture of hydrocarbons. The debutanizer (DB) separates pentane and heavier components from the mixture which are used as fuel gas. The Cr-C4 splitter (SP) separates propane and lighter components, with the overhead product priced as pure propane. The overhead and bottoms products from the deisobutanizer (DIB) are priced as pure iC4 and nC4 respectively. The economic data are taken from Bailey et al. (1993). The columns are modeled using tray-by-tray models and the thermodynamic proper-

continuation method traverses the optimal solution path has to be reversed. The path following continues for negative values of v in order to avoid division by zero and then switches back to system (3) (point G). The transition from system (10) to (3) is

ties are calculated using regressed equations. The system originally has six degrees of freedom and the setpoints to the multivariable controller are:

made by dividing the optimality conditions by the current value of ~,, thus normalizing v and switching back to the KKT conditions. The optimal solution path consists of local maximizers since all Lagrange multipliers and all eigenvalues of the reduced Hessan changed sign (v is negative). Apparently, the local maximizers do not have any significant

(a) the reflux to feed ratio in the DB; (b) the nC4 molar fraction in the bottoms product in the DB; (e) the iC4 molar fraction in the overhead product in the SP; (d) the C3 molar fraction in the bottoms product in the SP;

1198

P. SEt~RLISand A. N. HRYMAK

(e) the nC4 molar fraction in the overhead product in the DIB; (f) the iC4 molar fraction in the bottoms product in the DIB. In this case, the /C4 molar flowrate in the feed increases from a base case value of 3.975 Mmol/d to

CONCLUSIONS In this work, the application of a parametric

sensitivity analysis based on continuation methods in chemical processes is studied. The procedure is a final target value of 4.3 Mmol/d linearly, while the applied to single and multiple unit processes under nC decreases at the same rate from 1.rrjV2 multiple model parameter variations of finite magni'1"/2 1/2 af72, to 1.4322 Mmol/d, so that the total feed flowrate in the column train is kept constant. The upper and lower bounds for the variables are shown in Table 4. The optimal solution for the initial values of the feed composition (iC4=3.975, nC4= 1.7559 Mmol/d), which is obtained using MINOS 5.3, has two degrees of freedom with four variables found at their bounds. The total number of equality constraints in the model are 855 and the equation set (3) has 1715 equations, The light and heavy keys for the DIB are iC4 and nC4, respectively. As a result, the main effects of the variation in the feed composition are initially noticed in the DIB where the bottoms flowrate reaches its lower bound and the bottoms purity level its upper bound for very small changes in the feed composition (0.2% increase in iC4). The system at this point has no degrees of freedom with the DIB overconstrained (three active variable bounds). In order to compensate for the continuing changes in the feed the active bound on the upper level of the bottoms flow rate in the DB becomes inactive. The profit starts to decrease, since the DB bottoms product has the highest molar value. As a result, heavier components (mainly pentanes) are sent into the SP and DIB so that the level of the iC4 in the DIB bottoms could satisfy the product specification, The DIB column is sensitive to this particular kind of feed composition disturbance and forces the degrees of freedom in the other columns to adjust so that the entire system remains feasible, At point K, that correspond to a relative change for iC4 equal to +8.1% and for nC4 equal to 18.4%, the upper bound on the molar fraction of pentanes and higher components in the DIB bottoms product are reached. The system could not absorb any further changes in the feed stream cornposition and there is a sequence of active set changes for infinitesimal changes in the feed. Table 5 shows the details of the type of singniarities encountered in the optimal solution path for the column train. Figures 11a-e show the optimal solution path for the profit, the bottoms product flowrate and the reflux to feed flowrate ratio in the DB, the reflux ratio in the DIB and the iC4 molar fraction in the overhead product of the SP, resvectivelv,

tude. The method is based on the solution of the parameterized set of optimality conditions for different values of the parameter set. The addition of the strict complementarity condition results in a sequence of KKT points. Multiple parameter perturbations along specified directions are possible by proper addition of equations that express the relation between the moves of the model parameters and an independently varying artificial parameter. These equations may reflect any frequently occurring disturbance or correlation between the parameters. The solution of the equation system is obtained using a predictor-corrector continuation algorithm (PITCON). The procedure handles active set changes by checking the feasibility relations and the SC condition at every continuation point in the solution curve. The equation set of the stationarity conditions is properly modified when the linear independence constraint qualification is violated so that the solution path is traced close to the singular point. Satisfaction of the second-order optimality conditions is checked at a point where there is evidence that the reduced Hessian of the Lagrangian is singular. Such a situation occurs when there is a limit point for the independent artificial parameter in the optimal solution curve or a singularity due to violation of either the SC and LI conditions. The evaluation of the second-order information is based on the calculation of a matrix whose columns form a basis for the null space of the active constraint gradients. A partition of the decision variables into independent and dependent variables is required for this purpose. The method allows the study of the optimal solution behavior for large-scale process systems. The parameter value ranges for which a different basis is required can be calculated. The feasibility region of a chemical process system is determined under parameter uncertainty or disturbances.

Acknowledgements--The authors are grateful for financial supportfor the research from an Imperial Off University Research Grant and the McMaster Advanced Control Consortium (MACC). NOMENCLATURE Nonlinear parameter variation function f= Obieetive f i m e t i n n

d =

Sensitivity analysis for chemical process optimization F = Parameterized set of stationarity conditions Fu = Parameterized set of stationarity conditions around a point that LI does not hold F ~ = Parameterized set of stationarity conditions with multiple parameter variations g = Inequality constraints h = Equality constraints I= Equality constraints index set J=Inequality constraints index set 1 = Number of inequality constraints in (P1) L = Lagrangian function of NLP m = Number of equality constraints in (P1) n = Number of decision variables in (P1) Q~ = Reboiler's heat duty (MJ/d) QD = Condenser's heat duty (MJ/d) r = Number of active inequality constraints RF = Reflux ratio x = Decision variables xB = Component molar fraction in bottoms product XD= Component molar fraction in overhead product z = Solution point of system F = 0 Z = Reduced gradient

Greek symbols a = Artificial parameter for multiple parameter changes fl =Constant ~,= Perturbation size for the numerical evaluation of the reduced Hessian 6 = Predictor stage step size e=Parametervector 0=Parameter perturbation vector 2 = Lagrange multipliers for equality constraints = Lagrange multipliers for inequality constraints v = Multiplier for objective gradient

Superscripts k = Continuation point index

Subscripts dep = Dependent decision variables i = Equality constraint index ind=Independent decision variables j = Inequality constraint index

Operators f = scalar x = n-dimensional vector V~f= (1 x n) row vector with elements Of/Oxi, i = 1, n .... llEI~llENCES Allgower E. L. and K. Georg, Numerical Continuation Methods. Springer-Verlag, New York (1990). Atherton R, W., R. B. Schainker and E. R. Ducot, On the statistical sensitivity analysis of models for chemical kinetics. A. 1. Ch.E. J 21,441-448 (1975). Bailey J. K., A. N. Hrymak, S. S. Treiber and R . B . Hawkins, Nonlinear optimization of a hydrocracker fractionation plant. Computers chem. Engng 17, 123-138 (1993). Caracotsios M. and W. E. Stewart, Sensitivity analysis of initial value problems with mixed ODEs and algebraic equations. Computers chem. Engng 9, 359-365 (1985). Cukier R. I., H. B. Levine and K. E. Sehuler, Nonlinear sensitivity analysis of multiparameter model systems. J. Comp. Phys. 26, 1-42 (1978). Dongarra J. J., J. R. Bunch, C. B. Moler and G . W . Stewart, LINPACK User's Guide. Society for Industrial and Applied Mathematics, Philadelphia (1979). Fiacco A. V., Introduction to Sensitivity and Stability Analysis in Nonlinear Programming. Academic Press, New York (1983).

1199

Fletcher R., Practical Methods of Optimization, 2nd edn. Wiley, New York (1987). Forbes J.. F. and T. E. Marlin, Model accuracy for economic optimizing controllers: the bias update case. Ind. Engng Chem. 33, 1919-1929 (1994). Ganesh N. and L. T. Biegler, A reduced Hessian strategy for sensitivity analysis of optimal flowsheets. A. I. Ch. E. J 33, 282-296 (1987). Gazi E., W. D. Seider and L. H. Ungar, Controller verification using qualitative reasoning. Proc. 2nd IFAC

Workshop on Comp. Soft. Struc. Integ, AI/KBS Sys. in Proc. Control, Lund, Sweden (1994). Guddat J., F. G. Vasqucz and H. T. Jongen, Parametric Optimization: Singularities, Pathfollowing and Jumps. Wiley, New York (1990). Hearne J. W,, Sensitivity analysis of parameter combinations. Appl. Math. Modeling 9, 106-109 (1985). Hirabayashi R., M. Shida and S. Shindoh, Manifold structure of the Karush-Kuhn-Tucker stationary solution set with two parameters. SIAM J. Opt. 3, 564-581 (1993). Holland C. D., Fundamentals of Multicomponent Distillation. McGraw-Hill, New York (1981). Hwang J. T., E. P. Dougherty, S. Rabitz and H. Rabitz, The Green's function method of sensitivity analysis in chemical kinetics. J. Chem. Phys. 69, 5180--5190 (1978). Jongen H. Th., P. Jonker and F. Twilt, Critical sets in parametric optimization. Math. Program. 34, 333-353 (1986). Jongen H. Th. and G.-W. Weber, On parametric nonlinear programming. Annl. Op. Res. 27,253-284 (1990). Keller H. B., Lectures on Numerical Methods in Bifurcation Problems. Springer-Verlag, New York (1987). Koda M., A. H. Dogru and J. H. Seinfeld, Sensitivity analysis of partial differential equations with application to reaction and diffusion processes. J. Comput. Phys. 30, 259-282 (1979a). Koda M., G, J. McRae and J. H. Seinfeid, Automatic sensitivity analysis of kinetic mechanisms. Int. J. Chem. Kinetics 11,427 444 (1979b). Koda M. and J. H. Seinfeld, Sensitivity analysis of distributed parameter systems. IEEE Trans. Automatic Control AC-27(4), 951-955 (1982). Kojima M. and R. Hirabayashi, Continuous deformation of nonlinear programs. Math. Program. Study 21,150198 (1984). Koninckx J., On-line optimization of chemical plants using steady-state models. Ph.D. Thesis, University of Maryland (1988). Kovach III, J. W. and W. D. Seider, Heterogeneous azeotropic distillation: homotopy-continuation methods. Computers chem. Engng II, 593-605 (1987). Leis J. R., S. A. Gailagher and M. A. Kramer, Parametric sensitivity analysis of complex process flowsheets using sequential modular simulation. Computers chem. Engng 11,409-421 (1987). Leis J. R. and M. A. Kramer, Sensitivity analysis of systems of differential and algebraic equations. Computers chem. Engng 9, 93-96 (1985). Lockett M. J., Distillation Tray Fundamentals. Cambridge University Press, U,K. (1988). Lundberg B. N. and A. B. Poore, Numerical continuation and singularity detection methods for parametric nonlinear programming. SIAM J. Opt. 3, 134-154 (1993). Mangasarian O. L. and S. Fromovitz, The Fritz-John necessary optimality conditions in the presence of equaiity and inequality constraints. J. Math. Anal. Appl. 17, 37-47 (1967). McCormick G. P., Optimality criteria in nonlinear programming. In Nonlinear Programming SIAM-AMS Proc. Vol. IX (1976). Murtagh B. A. and M. A. Saunders, MINOS 5.3 User's Guide. Technical Report SOL 83-20R (1975).

1200

P. SEVEmZSand A. N. HRYMAK

Poore A. B. and C. A. Tiahrt, Bifurcation problems in nonlinear parametric programming. Math. Program. 39, 189-205 (1987). Prasad R., An improved approximation for the log-mean temperature difference. Chem. Engng 95(14), 110 (1988). Rakowska J., R. T. Haftka and L. T. Watson, Tracing the efficient curve for multi-objective control-structure optimization. Comput. Sys. Engng 2, 461-471 (1991). Rheinboldt W. C., Numerical Analysis of Parametrized Nonlinear Equations. Wiley, New York (1986). Rheinboldt W. C. and J. V. Burkardt, A locally parameterized continuation process. ACM Tram. Math. Software 9, 215-235 (1983a). Rheinboldt W. C. and J. V. Burkardt, Algorithm 596. A program for a locally parameterized continuation process. ACM Trans. Math. Software 9, 236-241 (1983b). Salgovic A., V. Hlavacek and J. Ilavsky, Global simulation of countercurrent separation processes via oneparameter imbedding techniques. Chem. Engng Sci. 36, 1599-1604 (1981). Seader J. D., M. Kuno, W.-J. Lin, S. A. Johnson, K. Unsworth and J. W. Wiskin, Mapped continuation methods for computing all solutions to general systems of nonlinear equations. Computers chem. Engng 14, 7185 (1990). Seider W. D., D. D., Brengel and S. Widagdo, Nonlinear analysis in process design. A. I. Ch.E. J 37, 1-38 (1991). Shapiro A., Second order sensitivity analysis and asymptotic theory of parametrized nonlinear programs. Math.

Program.33, 280-299 (1985). Tiahrt C. A. and A. B. Poore, A bifurcation analysis of the nonlinear parametric programming problem. Math. Program.47, 117-141 (1990). Tilden J. W., V. Costanza, G. J. McRae and J. H. Seinfeid, Sensitivity analysis of chemically reacting systems. In K. H. Ebert, P. Deuflhand and W. Jager (eds), Modeling of Chemical Reaction Systems. Springer-Verlag, Berlin (1981). Uber J. G. and E. D. Brill Jr, Design optimization with sensitivity constraints. Eng Opt. 16, 15--28 (1990). Volin Y. M. and G. M. Ostrovskiy, Sensitivity calculation methods for complex chemical systems - - 1. Algorithmization. Computers chem. Engng 5, 21-30 (1981a). Volin Y. M. and G. M. Ostrovskiy, Sensitivity calculation methods for complex chemical systems - - 2. Comparison. Computerschem. Engng 5, 31-40 (1981b). Wacholder E. and J. Dayan, Appfication of the adjoint sensitivity method to the analysis of a supersonic ejector. Trans. ASME J. Fluids Engng 106, 425-429 (1984). Wayburn T. L. and J. D. Seader, Homotopy continuation methods for computer-aided process design. Computers chem. Engng 11, 7-25 (1987). Williams T. J. and R. E. Otto, A generalized chemical processing model for the investigation of computer control. AIEE Trans. 78, 458 (1960). Wolbert D., X. Joulia, B., Koehret and L. T. Biegler, Optimal flowsheet sensitivity in a sensitivity oriented environment. Computers Chem. Engng lg, 1083-1095 (1993).