Model simplification and reduction for simulation and optimization of chemical processes

Computersthem. Engng,Vol. 17, No. 516, PP. 465483, Printedin Great Britain.All rightsreserved

MODEL FOR

0098-1354/93 $6.00 + 0.00 Copyright Q 1993 PergamonPressLtd

1993

SIMPLIFICATION AND REDUCTION SIMULATION AND OPTIMIZATION OF CHEMICAL PROCESSES J . PERREGAARD

Engineering Research Center, IVC-SEP, Institut for Kemiteknik, The Technical University of Denmark, DK-2800 Lyngby, Denmark (Received

1 September

1992; received for publication

30 September

1992)

Abstract-This paper presents a general model simplification procedure applicable for problems related to simulation (steady-state and dynamic), control, design/analysis and optimization of chemical processes. It also presents a partition technique which is applicable in equation-oriented simulations, and a model

reduction procedure for process optimization with modular based simulators. Several examples illustrating the applications of model simplification, partition techniques and model reduction are presented.

INTRODUCTION

Efficient solution of any simulation (steady-state and dynamic) problem and/or optimization problem, requires suitable numerical methods with special computational features. For example, for steady-state simulation and optimization using the modular approach an efficient coupling of the optimization algorithm with the simulator has to address questions related to the calculation order of the units, convergence of the flowsheet (full, partial or no convergence of the flowsheet at each optimization step), method of convergence, as well as the optimization strategy and location of this method as a unit in the flowsheet. For dynamic simulation using an equation-oriented simulator, issues related to the solution mode, sparsity or partitioning of the Jacobian matrix, consistent initialization and numerical stability have importance with respect to the development of robust and efficient methods of solution. By solution mode, it is meant whether the model equations represented by a set of differential-algebraic equations should be solved in the DAE mode or in the ODE plus procedures mode (solve the differential equations separately from the algebraic equations). Through model simplification, the complexity and nonlinearity of model equations are reduced. Model reduction on the other hand, reduces the number of equations that have to be solved. Model simplification and/or reduction therefore allow opportunities for reduction of the computational burden if and when their applications are feasible. The objective of this paper is to present a model simplification procedure that through simplification of the algebraic equations for phase equilibria calculations is capable of reducing the computing time to

solve the model equations without affecting the convergence characteristics of the numerical method. Also, this paper presents an algorithm for partitioned integration, which through using the inherited structure of the equations representing the chemical process, partitions the equations into subsets of equations that can be solved separately. For optimizusing the modular simulation ation problems approach an algorithm for determination of the unit calculation order that gives the minimum number of unit operation evaluations required in the solution of an optimization problem is presented. Also, for steady-state optimization, a model reduction procedure which reduces the number of algebraic equations based on the formulation of the optimization problem is presented. The application of the partitioned integration technique is illustrated through dynamic simulation of a process for esterification of acetic acid. The process consists of a mixer, a reactor converting ethanol and acetic acid to ethyl acetate and water, and a distillation tower for separation of the products from the unreacted feed components. This chemical process is also used fo illustration of the model reduction procedure presented in connection with optimization problems. The application of the model simplification procedure to problems related to optimization is illustrated for two chemical processes. The first process is a separation process where a feed of five normal alkanes is separated into four products using a distillation tower. The second proce&concerns the isomerization of n-butane. This process involves a reactor converting n-butane to isobutane and two distillation towers for separation of the products. The unreacted n -pentane is recycled.

J.

466

The application

of the model

simplification

PERREGAARD

pro-

cedure for dynamic simulation problems is also illustrated through simulation of two chemical processes. The first is a process where a four-component mixture of hydrocarbons is separated in a distillation tower and the second is a process involving two flash unit-operations with one recycle. DYNAMIC SIMULATION WITH STRUCTURED APPROACH

THE

Using the general framework of an equation oriented (EO) simulator, Gani et al. (1990) have demonstrated that increased flexibility and efficiency can be achieved by applying a structured approach. In the structured approach, the model equations are classified into three types, ordinary differential equations (ODES), explicit algebraic equations (EAEs) and implicit algebraic equations (IAEs). The different classes of equations are also grouped into those that are local to a single unit-operation only @ME, single model-equations) and those that are related to more than one unit-operation model (LE, linked-equations). Simulation strategies and therefore, methods of solution, can be derived for various simulation problems with different combinations of ODES, EAEs, IAEs and their grouping into SME and LE. The structured approach therefore allows a high degree of flexibility while maintaining the efficiency of any EO method of solution. Also, different combination and grouping of the model equations allow one to solve the equations in the DAE-mode (solving the differential-algebraic equations as one set), the ODE-mode or the AE-mode (solving only algebraic equations). In the structured approach, the generic form of the ODES, IAEs and EAEs representing a chemical process are given by equations (l-3) respectively:

dy

-

dt

= f(Y, p, t>,

(1)

a,

(2)

0 = g(y, x,

In the DAE-, type:

ODE-

P = h(x, Y, z).

(3)

or AE-mode,

an equation of the

D % = my*,

P, t)

is solved. In equation (4), the elements of y* determine the type of mode of solution. Due to the inherent modular structure of chemical process flowsheets, some of the ODES are more heavily coupled than others and the objective of a partitioned technique is to exploit this modular structure. In the

partitioned technique therefore, smaller ODES, IAEs and EAEs are solved.

subsets of

Partitioning In the predictor-corrector based integration methods, each step of integration employs a predictor-step, where an initial estimate for ~(t + h) is computed explicitly and a subsequent corrector-step where a set of equations of the form given in equation (5) is solved: o=~o~h*f(y)+D(!P

-y).

(5)

Equation (5) in this paper is called the generalized corrector equation (GCE) and the partitioned technique deals with how the set of equations represented by equation (5) are solved. If equation (5) is split into np (precedence-ordered) partitions, the set of GCEs for each partition takes the form: O = BO.h sfd(Yi, Pi, PY) + Di(Yi

- Yi),

(‘5)

where pi is a vector function of variables belonging to partitions with indices lower than i and pi”is a vector function of variables belonging to partitions with indices higher than i. The problem is then to determine a solution to equation (6) for all np partitions. For partitions belonging to a procedure order which does not contain recycle streams between the partitions, the vector p” will contain no elements and the partitions can be solved sequentially. However, if the recycle streams between the partitions are not “too sensitive” it is possible to solve the partitionequations sequentially by using an approxiation p,!‘, as shown in equation (7): O=~o~h~~(y,,pf,p,“‘)+Di(!P,-yi)=ri*

(7)

In the above equation p:’ = py(y’) where y” is the value of y after the predictor step. Sorensen (1991) has proposed a mathematical expression for a bound on the error committed by solving equation (7) instead of equation (6). According to SPrrensen (1991) this error can be made arbitrarily small by selecting a sufficiently small error bound for the integration. The step-by-step procedure for integration in the partitioned mode is discussed in detail by Sorensen (1991). In this paper an update version of the partition mode algorithm is given in Algorithm I. It should be noted that, since only the block diagonal elements (the derivatives of functions in partition i with respect to variables in partition i) of the total Jacobian matrix are required, it is possible to (artificially) decouple equations by keeping stream variables between partitions constant and increase the number of variables which can be perturbed simultaneously. When the number of variables perturbed simultaneously is increased a reduction of the

Model simplification and reduction computational

burden is obtained.

Since each par-

tition is solved separately, but simultaneously, the partition-mode also makes it feasible to apply different integration-methods for each partition. The “nonstiff partitions are to solved using a “nonstiff’ integration method and the “stiff’ partitions are solved using a “stiff’ integration method. Algorithm I. The partitioned-mode

algorithm:

1. Given u(r), t and the integration step size h. 2. Predict r(t + h), and set the partition index ip = 1. 3. Update Jacobian matrix if required. 4. Solve equation (7) for y,(t + h) and compute output stream variables from partition ip. 5. Set ip = ip + 1 and repeat from Step 4 unless ip > np. 6. In case of convergence failure, reduce h and restart from Step 2. Otherwise, continue. 7. Set t = t + h and update step size if required. 8. If t < tena repeat from Step 2, otherwise, stop. MODEL

SIMPLIFICATTON

Steady-state and dynamic simulation of chemical processes where thermodynamic models are employed to predict the physical properties, up to 30-80% (Grens, 1983) of the computational time is spent on these predictions. It is therefore obvious that a reduction in the thermodynamic model complexity may have a great impact on the overall computational burden. Indeed, earlier works by Johns and Vadhwans (1986), Macchietto et al. (1986) and Hillestad et al. (1990) have proposed methods for thermodynamic model simplification in process simulation. Issues related to the form of the simplifiedmodel (for example, what is the correlation capability of the simplified model?) and when and how to update the simplified model without causing any instability of the numerical method however, still need to be resolved. The technique developed in this work addresses the above questions and the important features of this technique is described below. Also, comparisons are made with model-simplification methods proposed earlier (Johns and Vadhwans, 1986; Macchietto er al., 1986; Hillestad et al., 1990). Simulation with combined models

simplified and rigorous

In the general framework of an equation-oriented simulator (steady-state or dynamic) a “Newton-like” method is usually employed to solve a set of noniinear equations. Each new estimate of the set of unknown variables are ‘obtained from the current

values

of

467

the variables

and the current

Jacobian

matrix and, function values of the set of nonlinear equations. In this paper, we propose to use the simplified model only to determine the Jacobian matrix but not for the evaluation of the function values. That is, for the function values, the rigorous model is used. If it can be guaranteed that the Jacobian matrix with the simplified model is correct, then the convergence characteristics of the method of solution (with the rigorous model) are not affected since the function values are not affected (by the simplified model). Therefore, it is possible to reduce the computational time without causing any numerical instability to the method of the solution. This can be seen from the fact that in any “Newton-like” method the updating of variables can in principle be written as: -F(Yk-‘);

n,-Yyk-‘].J(Yk-‘)=

&=g,

(8) .I

where Yk is the vector of the variables at iteration k, F is the vector of residuals and JI is the Jacobian matrix. Since the simplification procedure does not effect the evaluation of F but only .J, the norm of the updating of Y only becomes zero when the norm of F becomes zero. Since F is evaluated using the rigorous thermodynamic model the solution found will correspond exact@ to the solution obtained when the rigorous model is used in all calculations. It should be noted that this simplification procedure does not change the character of the corrector step of predictor-corrector integration methods such as BDF-methods (Gear, 1971). In the solution of the generalized corrector equation, usually, a Newtonlike method is employed. This results in the use of an approximated Jacobian matrix. Using the simplification procedure, a very “good” approximation to the Jacobian matrix is obtained. It should also be noted that the simplification procedure provides additional savings since the simplified thermodynamic model is only used when a full Jacobian evaluation is required. The same is true also for steady-state simulation with an equation-oriented approach and using a quasiNewton method. Form of the simplt$ed thermodynamic model The form of the simplified model is important because it is closely related to the correlation capabilities of the simplified model. Generally, the calculation of equilibrium ratios (Ki) in separation process models involves the solution of a set of nonlinear IAEs. This can be written as: Ki=f

(T, P,n).

(9)

J.

468

F'ERRFGAARD

In this equation, Ki is the equilibrium ratio for component i, T is the temperature, P is the pressure and n is a vector of mole numbers of each component (from which the mole fractions can be evaluated) in the unit-operation model (or part of a unit-operation model) under consideration. This general expression for the equilibrium ratio implies that the Ki is a function of temperature, pressure and composition. A simplified thermodynamic model is a simple analytical function, fitted to the desired property values predicted by the rigorous thermodynamic model within a limited range of temperature, pressure and composition. In the development of a simplified expression for the equilibrium ratio it is essential that the expression is explicit in order to ensure speed of calculation. It

The first system to be examined in this paper is the binary azeotropic system consisting of acetone and chloroform. Using the MHV2 model (Dahl et al., 1991) K-values were evaluated at a constant pressure of 1 atm over the entire liquid composition range for the binary system. This of course means that the temperatures are different at each liquid composition. The pressure term in the simplified model was excluded in all calculations using the simplified model (Di= 0). The effect of the composition-dependent term can be evaluated through calculations excluding (Cj = 0) and including this term. Therefore the generated K-values were correlated using a two- (A and B) and a three-parameter (A, B and C) version of the simplified model. The results are illustrated in Fig. 1. The difference in the ability of the two simplified

is also very important that the parameters to be fitted to the rigorous model are linear in the expression, in order to ensure a fast, robust and consistent evaluation of the parameters. Most important, however, is the physical meaning of the terms in the expression, since this determines the predictive capabilities of the simplified model. Earlier workers have considered the logarithmic of the K-value as an explicit expression of temperature, pressure, and composition as the simplified thermodynamic-models [see, for example, Grens (1983) Johns and Vadhwans (1986), Macchietto et al. (1986), Hillestad et al. (1990)] Simplified models using reference compositions and reductions of the number of components, through use of pseudo-components (Chimowitz et al.,1983) have also been suggested. The form of the simplified model used in this work is given in equation (10) and is based on the Wilson expression (Wilson, 1968) for K-values. However, in addition to the Wilson expression a composition dependent term is added:

models to represent the K-values is dramatic. The model using only temperature dependence shows significant errors as well as the wrong curvature for acetone. While the model including both temperature and composition dependence represents the K-values quite well. The fact that the two-parameter model for a binary azeotropic mixture cannot give the correct curvature of the K-values for both components is to be expected, since minimum (or maximum) in K-values at the azeotropic temperature is present and the simplified expression using only the temperature dependence is monotonic in the temperature. Therefore, the nonideal composition term is required. The azeotropic system described above has also been studied by Chimowitz et al.(1983) in the context of application of simplified thermodynamic models.

Ln(Ki)=A,+$+Ci[l

-zT+Q.Ln[$].

(10)

In equation (10) A,, B,, C, and D, are constants fitted to a rigorous model and superscript 0 indicates the reference pressure. Since the Wilson expression is known to provide a initial good guess for rigorous equilibrium calculations, it therefore also constitutes a good basis for developing simplified models. Since in many cases it is necessary to take the cross-component interaction into account, an additional term for incorporating this cross-component interaction is introduced. This is the term with the parameter C in equation (10). The necessity of using a compositiondependent term along with the general correlation capability of the simplified thermodynamic model is explained below.

0.00

.

uq”i~wlno~:w~ac~:n or0 @Let::

Fig. 1. Comparison of simplified thermodynamic models for representing vapor-liquid equilibria in the binary azeotropic system consisting of acetone and chloroform at 1 atm.

Model simplification and reduction

469

His conclusions are in agreement with the ones made here and it should be mentioned that the composition-dependent simplified model used by Chimowitz er al. (1983) gives an even better representation of the et al. system. However, the model of Chimowitz (1983) uses a reference composition and through this introduces three additional parameters (50% more parameters). Also, the use of reference compositions for fitting parameters does not make it feasible to employ the simplified model in regions other than the one where the parameters are obtained. This is illustrated in the next example. This quaternary system ethanol( 1), 1-propanol(2), cyclohexane(3) and toluene(4) was also studied by Chimowitz et al. (1983). Here the same calculations as those performed by Chimowitz et al. (1983) is performed but using the simplified thermodynamicmodel represented by equation (10). At a fixed total composition (given in Table 1 as mixture I) and a constant pressure of 10atm the TXY-diagram is determined using a rigorous thermodynamic model. These TXY-data points are then used to fit the parameters Ai, Bi and Ci in equation (10) for each component and the fitted simplified model is then used to predict the K-values for mixtures II and III (also given in Table 1). Two different rigorous models, the MHVZ-model and the original UNIFAC-model (Fredenslund et al., 1977) were used in the generation of the data-points. The difference between the two models at this pressure (10 atm) is mainly the fact that the MHV2 model uses three temperature-dependent parameters whereas the original UNIFAC uses only one temperature-independent parameter. The results are illustrated in Figs 2 and 3 for the MHVZmodel and in Figs 4 and 5 for the original UNIFAC-model. It can be seen from these figures that the values and the dependency of the K-values for all components are different for the two rigorous models. An analysis of the applicability of the models (rigorous models) for the evaluation of equilibrium for this system therefore requires experimental data. An analysis of the capability of a simplified model can, however, still be performed using this system. The Figs 4 and 5 correspond nearly exactly to the results of Chimowitz et al. (1983) suggesting that the

Table 1. Composition

of mixtures 1, 11 and 111

I

II 111

1 0.50

1

0.45 0.45 0.50 0.53 0.55 0.55 0.50 0.43 0.55 0.55 Mole fraction of tolueae In the liquid phua

Fig. 2. Vapor-liquid equilibrium ratios for mixture II in the quaternary system consisting of ethanol, l-propanol, cyclohexane and toluene at 10atm. The full lines are predicted values using the MHV2 model and the dotted lines are predicted using the simplified model with parameters evaluated based on mixture I. rigorous model used by Chimowitz et al. (1983) has no temperature-dependent parameters. The simplified-model with parameters obtained based on mixture I is for both the rigorous-model, MHV2 and the rigorous model, original UNIFAC, capable of prediction of the K-values for mixture II with a very 3.0

_ --

Rigorous

2.5

4 i

2.0

%

r;

81.5 B 8 1.0

l&propanol

--. ___--4 C+hexane

1

Toluene ____----

0.25

Z(2)

213)

214)

TD

TB

0.25 0.18 0.10

0.15 0.12 0.30

0.20 0.22 0.32

0.40 0.48 0.28

460.8 467.3 458.3

443.6 450.2 448.2

__.

,,,,,,,,,,,,,,,,,,,,,,,,,,I,,(, 0.50

Yale fraction

Z(l)

value0

Simplified prediottons

z

0.5

Total composition of the mixtures for the system ethanol(l), I-propanol(2). cyclohexam(3) and toluene(4) TD is the dew point temperature and Ts is the bubble ooint temDerature Mixture

1 .QO

ol

025 tolueae

I ,,,,,, in

0.40 the Uqtid

0.45 pham

Fig. 3. Vapor-liquid equilibrium ratios for mixture III in the quaternary system consisting of ethanol, l-propanol, cyclohexane and toluene at 10 atm. The full lines are predicted values using the MHV2 model and the dotted lines are predicted using the simplifkd model with parameters evaluated based on mixture I.

470

J. 12

J

Rigorous

-

-0

PERREGAARD

Slmplifiad

veluce

predictions

Toluene ,““,““,““,““,““,““,“‘~

0.45 Mole

0x4

0.56

0.60

0.68

0.70

0.76

fraction of toluenc in the liquid

a

case

Fig. 4. Vapor-liquid equilibrium ratios for mixture II in the quatemary system consisting of ethanol, 1-propanol, cyclohexane and toluene at 10 atm. The full lines are predicted values using the original UNIFAC model and the dotted lines are predicted using the simplified model with parameters evaluated based on mixture I.

high accuracy. The compositions of mixture II differ from mixture I by lO--20% but the dominating component (toluene) is the same. In contrast, it can be seen from Figs 3 and 5 that the errors with the simplified-model are greater compared to the

--

.

.

rigorous-models for mixture III. This mixture is very different from mixture I and the dominant components are no longer toluene. From this example it can be seen that even with composition-dependent simplified models evaluating the parameters in a simplified model based on data in one region and applying the simplified model in another region can result in very significant errors. In the work of Chimowitz et al. (1983) the simplified model incorporates a reference composition using “extra” parameters fitted at this reference composition. This, however, did not result in an improvement in the predictions for these systems. From our experience with many mixtures that we have examined, the general conclusions are that for mixtures that are only moderately nonideal (hydrocarbon mixtures, mixtures with only one polar component, etc.) the simplified model is capable of correlating the equilibrium ratios to a very high accuracy (less than 1% deviation) over a large range of compositions, whereas the correlation of highly nonideal mixtures cannot be made with the same accuracy. For small composition ranges however, highly nonideal systems can also be correlated to a high degree of accuracy by incorporating a composition term. It should be noted that any procedure for applying simplified-models in the context of accurate process simulation problems must ensure that the region in which the simplified model is adapted is also the region where it is applied.

Updating Rlgomu6 v6luoa Simplttted predictions

.

.

.

.

.

0.m 0.m 0.00 0.a 0.46 a.60 0.00 Eclotmoticnoltolu.aehthollguldph6oe

. ada

Fig. 5. Vapor-liquid equilibrium ratios for mixture III in the quaternary system consisting of ethanol, I-propanol, cyclohexane and toluene at 10 atm. The full lines are predicted values of the original UNIFAC model and the dotted lines are predicted using the simplified model with parameters evaluated based on mixture I.

of the parameters

in the simpltjied

model

Application of any simplified model requires a procedure for updating of the parameters in the simplified model. Also, it is necessary to determine when to update the parameters and what criteria to apply for updating the parameters. If all parameters are linear in the simplified model, it is possible to update the parameters using the concept of least squares approximation (Hansen, 1985). Recently Macchietto et al. (1986) and Hillestad et al. (1990) have suggested the use of recursive least-squares procedures. The idea behind a recursive least-squares method is that the information of past points can be preserved in the covariance matrix for the parameters. When a new point is known, the covariance matrix can be updated and new values of the parameters can be obtained. One of the advantages of the recursive least-squares methods is that points from previous calculations using the rigorous model can be included in the parameters of the simplified model without storing these points but only storing the covariance matrix. However, this also results in one of the difficulties in applying a recursive leastsquares method. The problem is that the region in

471

Model simplification and reduction which the simulation problem is being solved may change during the iteration/time of solution of the

Criteria for up&ting

simulation problem. Thus, the simplified model parameters obtained at one point in the solution of the simulation problem (one region) may have no relevance at another point. Therefore the recursive least squares methods need to incorporate some way of “forgetting” past information. If, however, no new information is introduced due to little change in the can lead to “blow-up situregion this “forgetting” ations” where the parameters in the simplified model become extremely sensitive to even small deviations between the rigorous and simplified models. The procedure applied in this work for updating parameters in the simplified model involves equations (lo), (11), (12) and (13). It is possible to evaluate all parameters in the simplified model based only on one evaluation of the rigorous model which also supplies the derivatives with respect to temperature, pressure and comparisons. Since there are four equations and four unknown variables (parameters), the simplified model parameters can easily be determined. If, however, more than one rigorous model evaluation is performed, a least-squares approximation is applied. This is also the case if all the derivative information is unavailable and more than one rigorous model evaluation is required:

is important to know when to update the simplified model. Macchietto et al. (1986) and Hillestad et al. (1990) have proposed updating policies based on the estimation of error between the simplified- and the rigorous-models. According to these policies, the simplified model parameters are updated when the errror estimates, determined through a truncated Taylor series expansion around the last updated point, is greater than a user-specified tolerance (usually set to 1%). In steady-state and dynamic simulations, much of the computational burden is associated with the generation of the Jacobian-matrices related to the solution of the set of model equations. If the simplified model is used only for the generation of the Jacobian-matrices (as in this work), then the simplified model parameters must be updated every time a new Jacobian-matrix is needed by the numerical method. Since the proposed simplified model, in addition to matching the reference conditions also matches the derivatives with respect to temperature, pressure and composition, the Jacobian-matrix calculated with the simplified model can be expected to have correct values. Whether the Jacobian matrix is determined (with the simplified model) numerically or analytically, the computational time is greatly reduced (compared to use of the rigorous model). It should be noted that in this work the simplified model is being used only for the generation of Jacobian-matrices.

Zn,

a Ln(K) ’

8%

2

=-2ci

1-z

[

aLn(K,) a* a Ln(&)

aP

Bi =-?’ D,

=P’

1 ,

(11)

(12)

(13)

From the rigorous model [equation (911 the left-hand side of the equations (1 1), (12), (13) and (10) are known and used to evaluate parameters C,, Bi, Di and Ai, respectively. The evaluation of parameters in the simplified model by means of derivative information from the rigorous model ensure that the simplified model at least locally represents the derivatives of the rigorous-model correctly. This is important since for example the Newton step in a Newton method where the Jacobian-matrix calculated using the simplified model must result in the same directions as when the rigorous model is used. The implementation of the procedure for updating parameters also allows the use of one set of parameters for each component in a full or part of a unit-operation model. For example one set of parameters may be used for an entire distillation column.

When applying

the simplified model parameters

model simplification

techniques it

OPTIMIZATION USING THE FRAMEWORK SEQUENTIAL MODULAR METHOD

OF

THE

Advances in numerical optimization techniques in the last decade have made it possible to consider larger and more complex problems. One of the most widely used optimization methods is the successive quadratic programming method which is capable of incorporating constraints in the solution of the optimization problem [see for example Biegler and Cuthrell(l985)]. In order to solve optimization problems involving complex chemical processes it is necessary to couple the optimization method with simulation tools. In this work the modular-based simulator [SEPSIM, Andersen et al. (1991)] has been coupled with the successive quadratic programming method of Bielger and Cuthrell (1985). The coupling has been reported earlier by Perregaard et al. (1991) who proposed two different strategies for solution of the optimization problems (black-box and simultaneous). In the black-box strategy the simulation problem is solved in an inner-loop and the

J. PJZRREGAARD

472

optimization

problem

in an outer-loop,

while in the

simultaneous strategy (infeasible path) the convergence of tear-streams from the recycle loops are included as equality constraints in the optimization problem. For an ethcient solution of the optimization problem, the minimum number of unit-operation evaluations should not be related to flowsheet “tearing” with respect to process simulation. A rearrangement of the calculation order as well as identification of the unit operations required for the optimization loop is necessary. In this work an algorithm for determination of a calculation order that requires a minimum number of unit evaluations has been developed. The algorithm (see Algorithm II) also determines whether to use the black-box or simultaneous strategies. The criterion for using the black-box approach is that no partition can be found having both design variables and tear-stream(s). In this case no optimization variables directly affect the partitions with tear-streams and the initial guesses (taken from the previous sohition of the simulation problem) are very close to the true solution. Hence, only one or two iterations in the convergence of the tear-stream(s) is required making the black-box more efficient than the simultaneous strategy. The application of Algorithm II is illustrated with the help of two flowsheets shown in Figs 6 and 7. In the first example, a flowsheet with four units and one recycle loop is considered and the flowsheet can be partitioned into three partitions (see Fig. 6). In the second example (see Fig. 7), a flowsheet consisting of six units with three recycle loops is considered and the flowsheet can be partitioned into three partitions (see Fig. 7). For different definition of objective functions and design variables, unit calculation orders for a single iteration are given in Table 2, where also the unit caIculation orders based only on flowsheet information are given. From Table 2, it can be seen that the developed algorithm for determination of the minimum number of units required in the solution of the optimization problem is capable of reducing the number of units in many cases. For example, in the definition of the optimization problem represented by

P2

Pl u2

Eirl

3

u3

P3 88

3

/ \ Fig. 6. Flowsheet Example 1, P means partition, LJunit and S stream.

Pl

1

Fig. 7. Flowsheet Example 2, P means partition, U unit and S stream.

example 1.4 in Table 2, the proposed algorithm for each iteration in the optimization loop reduces the required unit-operation evaluations by 75%. It should also be mentioned that the algorithm is also capable of handling situations with more than one tear-stream per partition. For example, if stream S7 in the flowsheet given in Fig. 6 was recycled back to unit number 1 (Ul) then the total flowsheet would be treated as one partition but with two recycle streams in the same partition. Note that the optimization algorithm must be added to the flowsheet as a unit. Algorithm II. Determination of the sequence for optimization problems:

calculation

Determine the number of partitions and the calculation sequences of each partition along with the number of tear-streams in each partition. For each partition make a list of units where design variables are present and a list of the units which have output streams or equipment parameters needed for the evaluation of the objective function or the constraints included in the definition of the optimization problem. If the number of tear-streams or the number of partitions with both tear-streams and design variables are zero, use the black-box approach and place the optimization unit at the end of the flowsheet. Solve the resulting simulation problem. In the solution of the optimization problem, exclude partitions with no design variables if the partition does not follow a partition with design variables. Also exclude partitions which do not have design variables and are not needed for the evaluation of the objective function or the constraints The check can be based on the two lists generated in Step 2. Add to the list of partitions required in the optimization loop, all partitions following a partition with design variables

Model simplification aad reduction

473

Table 2. Unit calculation sequence Unit and streams used in the objective function of constraints Example 1.1 1.2 1.3 1.4 2.1

2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9

CvalUatiOU

Ul,S4, s2 u3,ss Ul, S6 U3, U4, S3, S4, S6 Ul, u2, u3, u4, Sl, s2, s3, S4, s5, S6, s7, SlO SlO SlO SIO Ul, U4, SlO, S6 Ul, U4, SIO, 56 SIO Ul, s3, s7 Ul, U3, U6,S9

Units where design variables

Optimization strategy BB

=-e-t

u2 u3 Ul.U2

ss X

X X X

z

X

Ul

X

z U6 Ul, u3 Ul, U6

:: X X X

::

::

unit calclllatioo order for the optimization using Algorithm II

Conventional unit calculation order for the simultaneous strategy

Ul,UZ u3 Ul, U2, U3, U4 u4 US, U6

Ul. u2 Ul. u2, u3 Ul, u2, u3. u4 Ul. u2, u3. u4 Ul-U6

ULU6 U3, U4, US, U6 US, U6 US, U6 UlLU6 U&U6 u3, u4 US. U6

Ul-U6 Ul-U6 Ul-U6 Ul-U6 W-U6 W-U6 Ul-U6 W-U6

Exampleswith

the first index equal 1 corresponds to tlowsheet given in Fig. 6. Examples with the first index equal 2 corresponds to Bowsheet 2 given in Fig. 7. U = unit; S = stream; X = use this strategy; UlLU4 = all units between Ul and U4 (Ul, U2, U3 and U4): BB = black-box strategy; SS = simultaneous strategy.

which has linked streams to partitions already in this list. For each iteration in the optimization loop, start the solution of the simulation prob-

variable has been changed. If more than one design variable has been changed start with the first unit in the flowsheet where a design

lem with evaluation of the unit where the last design variable has been changed. If more than one design variable has been changed, start the solution of the simulation problem with the first unit in the flowsheet where a design variable has been changed. Stop when the Kuhn-Tucker condition is satisfied. 4. If the number of partitions with both tearstreams and design variables are greater than zero, use the simultaneous strategy and place the optimization unit directly after the first tom stream in the first partition including design variables. Solve the resulting simulation problem. In the solution of the optimization problem, exclude partitions with no design variables if the partition does not follow a partition with design variables. Also exclude partitions which do not have design variables and are not needed for the evaluation of the objective function or the constraints. The check can be based on the two lists generated in Step 2. Add to the list of partitions required in the optimization loop all partitions following a partition with design variables which has linked streams to partitions already in the list of partition required in the optimization loop. Expand the partition with the optimization unit to include all following partitions not excluded above. In the optimization problem, start with the evaluation of the unit where the last design

variable has been changed. Make one evaluation of all units in the expanded optimization loop. Stop when the Kuhn-Tucker condition is satisfied. Model reduction in connection with steady-state timization

op-

The algorithm for determination of a calculation sequence of optimization problem can be thought of as a model reduction scheme, since by reducing the number of unit evaluations in the optimization problem, it also reduces the number of equations. A more traditional scheme for model reduction applicable for the simultaneous optimization strategy has also been developed_ Consider a flowsheet consisting of two partitions (pl and p2) where all design variables belong to the first partition but the evaluation of the objective function involves unit-operation and/or stream variables belonging to both partitions (for example, Problem 1.3 in Table 2). In this case the evaluation of the objective function requires an evaluation of all units and streams in the two partitions. It is, however, possible to apply model reduction for the second partition. The objective function can be calculated in two parts (Fl and F2). One based on the first partition pl, and one based on the second partition p2. Since no design variables are present in the second partition model reduction for this partition is feasible. Algorithm III, given below, transforms the rigorous evaluation of partition p2 into a simple evaluation based only on the stream variables entering p2.

J. PERRBOAARD

474 Algorithm III. Model optimization:

reduction in connection

with

a rigorous model If II@ -SW, > Cl Perform evaluation of partition p2, evaluate F2 and set S”=S. If e1 > Il(S - W’>II, > Q Perform a rigorous model evaluation of partition ~2, evaluate F2. Set F2rrf = F2 and set S& = S. For i = 1,2..nv. Perturb stream variable No, i and perform a rigorous model evaluation of p2 and evaluate F2. Calculate

AFp2.i GZ - FEZ x = and reset S.I-- SF’ I.

I

SF’ - si

If c2 > Il(S - S=p 11 Z Apply the reduced model evaluation of F according to:

In the above algorithm, S is a vector of stream variables entering ~2, c1 and + are predefined tolerances indicating when the reduced model has to be updated and applied. F2 is the part of the objective function due to partition ~2. The algorithm can be applied for all calculations in connection with the solution of an optimization problem of the nature stated above. The first step in the algorithm corresponds to a rigorous evaluation of the second partition and is used far from the optimal solution. The second step in the algorithm is the evaluation of the derivatives of F2 with respect to the stream variables entering the second partition. The aim is to use this step only once. In the third step it can be seen that the equation representing partition p2 is reduced to a simple summation. This step is used when the stream variables entering this partition are not changing too dramatically.

SimplijScation design

in

connection

with

optimization

/

Solving optimization problems using a successive quadratic programming method coupled with a modular-based steady-state simulator (in this work represented by SEPSIM) also makes it feasible to apply the proposed thermodynamic model simplification procedure. The solution of an optimization problem using a successive quadratic programming method consisting of three parts: 1. Gradient evaluation of objective function and constraints 2. Finding a search direction through solution of a quadratic programming problem. 3. Line search along the search direction.

The optimal solution is considered found when the Kuhn-Tucker criterion is satisfied. This criterion is based on zero gradient. In the solution of the quadratic programmi ng problem no flowsheet (and no part of the flowsheet) evaluation is required. Therefore no thermodynamic property evaluations are required_ But in both the gradient evaluation and the line search, flowsheet evaluation and therefore thermodynamic property evaluation are involved. Simplification of thermodynamic model equations is therefore feasible when the gradients are evaluated or a line search is performed. Since, a simplified thermodynamic model is only valid in a limited region, applying a simplified thermodynamic model in a line search where the region of operation may change drastically is dangerous. However, in the evaluation of the gradient where the region of operation is not changing (small perturbation around the last obtained point is performed when numerical derivatives are used) it is feasible and logical to apply simplified models. The proposed procedure for model simplification related to optimization is therefore applied in this work, only when a gradient evaluation is performed. Since the criteria for terminating the optimization calculations is based on gradients, it is not adequate to let the termination of the calculations be based on the simplified model, as it may result in a slightly different gradient vector. The differences in the gradient vector between the rigorous and the simplified model is only important close to the solution. Since the difference between the rigorous gradients and the simplified gradients evaluation is small this difference only becomes important when the norm of the gradients is small. This is exactly the case close to the optimal solution where it is required that the norm of the gradient is lower than an error tolerance. The behavior of the gradients in the 2-D space can be graphically illustrated as shown in Fig. 8. At the starting point and the first gradient evaluations the simplified and rigorous evaluation of the gradient vector results in practice in the same

Start

point

f Gradient Full

and

1 almpllfled

Fig. 8. The gradient in the solution of an optimization

problem.

Model simplification and reduction direction whereas close to the solution the small difference becomes important, since the norm of the vector is small. To overcome the problem of differences in the gradient vector the procedure for application of simplified model first determines the optimal solution using the simplified model in gradient evaluations, then an evaluation of the gradients using the rigorous model are made and the iteration for solution to the optimization problem is continued only using the rigorous model. The step-by-step procedure for applying the simplified thermodynamic model in connection with optimization problems using a successive quadratic programming method is outlined in Algorithm IV. The numerical experience using this procedure has shown that no more than two rigorous gradient evaluations are required. Algorithm IV. Algorithm for applying the simplified thermodynamic model in connection with optimization: (I) (II)

(III)

(Iv)

0 WI)

Set flag I indicating the use of simplified model to one. If I is equal to one update the A,, Bi, Ci and Dj parameters from the last rigorous model evaluation and evaluate the gradients applying equations (10). Otherwise, use equation (9) for evaluation of the gradient. If the Kuhn-Tucker criterion is satisfied and I is equal to one, set I to zero and reevaluate the gradient. If the Kuhn-Tucker criteria is satified then the optimum is found and the calculations stopped. Solve the quadratic programming problem to obtain a search direction. Perform a line search using the rigorousmodel [use equation (9)] and go to II for a new evaluation of the gradient. RESULTS

AND

DISCUSSION

475

of the bottom product from the fh-st column is recycled. A second distillation column (used only for steady-state optimization) consisting of 10 trays, separates ethyl acetate from water. Details of this process is given in Table 3 and the process flowsheet is shown in Fig. 9. Steady-state optimization (case study la). All results for steady-state optimization of the esterification process has been obtained using the modular steady-state simulator SEPSIM (Andersen, 1986). The reactor section has been simulated as a multiphase equilibrium reactor and for the distillation columns, constant equimolar overflow was assumed. An optimization of the full flowsheet partitioned into two subsections was performed. The pressure and temperature in the reactor were selected as design variables. The objective function was in this case defined as: F”@ = -(flow of ethyl acetate in top product DIST2-flow of acetic acid in bottom product of DISTI)

of

Therefore, variables from the second partition (the second column) enter into the evaluation of the objective function while the tear-stream variables belong only to the first partition. The calculation order based only on flowsheet information would be to solve the first partition until convergence and then to solve the second partition. If simultaneous updating of tear-streams and design variables are to be performed (Perregaard et al., 1991) this procedure is inadequate, since the evaluation of the objective function requires the values of streams from the second partition (which are not known when the first partition is solved). Therefore, using Algorithm II, a new calculation order is determined where one evaluation of the model equations for the second partition is combined with the solution of the first partition and the optimization problem. The definition of this optimization problem makes it feasible to apply the model-reduction algorithm (Algorithm III). When applying Algorithm III and using the values of the tolerances (defined in the

Case study 1. Esterz#cation of acetic acid An equimolar mixture of acetic acid and ethanol is being fed to a mixing tank in which recycled unreacted acetic acid and water also enter. The mixture is passed to a reactor and reacted at 300 K according the reaction scheme: C, H, OH + CH, COOH

--+ H, COOC,

The reactor products are distillation column, where taken out at the top of the and acetic acid is obtained

to

H, + H, 0.

separated in a IO tray acid-free ethyl acetate is column and mainly water as bottom product. 80%

Table 3. Design of the cstcrification

pr-ss

Feed stream Flow 100 kmol h-’

Temperature 297.0 K

Pressure l.OatHl

Composition (0.5,0.5,0.0,0.0)

Parameters for the tank Mixer: Vol = 18.4 m’, area = 6.62 m2, R, = 204 kmol h-’ II-’ Parameters for the CSTR Vol = 16.45 m3, area = 5.92 m2, R, = 204 kmol h-’ II-’ ks-8x 10Srn’kmol-‘h-‘, E,-59.45MJkmol-’ Design of the distillation tower I Tray configuration: IO sieve trays Weir length: 2.65111, Weir height: 0.19m Plate spacing: 0.16 m. Active area: 9.15 m’ Vapor boilup: 296.0 krnol h-‘, Distillate rate: 64.0 kmol h-’

476

J. P~REGAARD

Mixer

Reac

Fig. 9. Flowsheet of the esterification of acetic acid.

algorithm) of 6, = lo-* and e2 = low3 in the solution of the optimization problem a reduction of the computational time of 38% was achieved. In general it cannot be guaranteed that the use of a reduced-model finds the same optimal solution as a rigorous model, however for this example the deviation between the rigorous model solution and the reduced model solution in the outgoing streams of the process were all lower than 0.1%. Dynamic simulation (case study lb). In this case, the reactor was modeled as a CSTR. The pre-exponential factor, k, was determined by forcing the value of k, to match the conversion factor obtained from steady-state simulations of the process (see “steadystate optimization” section). The activation energy of the reaction EA is taken from Venkataraman et al. (1990). A linear pressure profile with atmospheric pressure at the top tray and 1.1 atm at the bottom tray was employed in the calculations. The second distillation tower was not considered (see Fig. 9). The strategy for the simulations is as follows: first the steady-state for the process at a reactor temperature of 300 K was determined. At the steady-state, the reactor temperature was suddenly increased to 310 K and the transient responses of the process were computed. Even for this simple chemical process, a fairly large set of equations (approx. 1000 equations) have to be solved simultaneously to simulate the process. Thus, for a purely equation-oriented approach quite a large Jacobian-matrix needs to be handled. In the case of the structured approach, however, only 49 ODES (5 ODES for the mixer unit, 4 ODES for the CSTR unit and 40 ODES for the distillation tower) have to be solved simultaneously in the ODE-mode and no iterations on stream values are required. In the DAE-mode the number of equations solved simultaneously is 101 (49 ODES and 52 IAEs). From flowsheet analysis the structure of the Jacobianpattern shown in Fig. 10 was obtained by employing the simple but efficient method of Sorensen (1991). Note that each of the white blocks in Fig. 10 contains only zero elements and that the Jacobian pattern does not take into account the equations describing each

wm

Oist l Fl

m

)

Fig. 10. Jacobian pattern for the esterification process.

module. This is the reason why the fraction of nonzero elements in the Jacobian-pattern is as high as 84% and therefore, the sparse mode cannot be expected to provide much saving compared to the dense mode. In Table 4 the efficiency of the partitioned technique outlined in Algorithm I is compared with the dense mode. In Table 4 dense corresponds to an evaluation of all elements in the Jacobian matrix, Pl corresponds to applying Algorithm I and P2 corresponds to applying Algorithm I with the extra feature that different integration methods are applied for each partition. For all simulation modes, no error failures and/or convergence failures were observed. The computed transient responses are identical within 0.1% for the different integration methods. Note that the step-size and necessity of Jacobian evaluations are determined on the basis of the computed data. In terms of efficiency, the two most important numbers NUE (the number of times a unit-operation model is evaluated) and CPU (the CPU-time in seconds needed for the simulation). These numbers favor the partitioned-mode. In fact, the reduction in the number of evaluations of unitoperation models (22%) shows the potential of this technique. The reduction in the number of evaluations of unit operation models is, however, mainly obtained in the first partition which involves the mixer unit and since the computational burden for the mixer unit is much smaller than that for the distillation unit, the reduction in terms of computational time cannot be expected to be too great for this example. However, as the size of the equation set increases, the efficiency of the partitioned technique should increase. Also for simulation of larger processes with similar dynamic characteristics, the partitioned-mode should be even more efficient. It can therefore be concluded that simulation in the Table 4. Integrationstatistics Method Dense

Pl P2

CPU

Mixer

NUE CSTR

DIST

Total

821 809 806

1287 485 491

1287 1285 1251

1287 1285 1251

3861 3055 2993

Model simplification and reduction partitioned-mode

a means for reduction

can provide

of the computational

time even for flowsheets

Four

with

recycle streams. Also, for the flowsheet of Fig. 9, the feasibility of solving different partitions using different integration methods is illustrated. Method P2 in Table 4 indicates the use of two different integration methods. The first partition (mixer) is integrated using a “nonstiff’ method and a “stiff’ method is used for the other partitions. As it can be seen from Table 4, this results in a higher efficiency.

1. 2. 3. 4.

477

optimization

(design)

in

steady-state

optimization

In Case studies 2 and 3 results for application of the model simplification procedure in connection with steady-state optimization problems are presented. Case study 2. Separation of hydrocarbons In a distillation tower a mixture of n-pentane, n-hexane, n-heptane, n-octane and n-nonane is separated into four products. The unit-operation model is a constant-molar overflow distillation model. Two different rigorous thermodynamic-models were used. In one case the SRK equation (Soave, 1972) of state was used and in another the MHV2 model was used. The details of the process is given in Fig. 11. The objective of the optimization is to obtain the best possible separation and the objective function (Fob) to be minimized is defined as: F”b’ = c x [&* + (G,

x s, + &-*

x s,

+ J&Q) x S.i x Xcg x

&I,

(14)

X, is the mole fraction of component i and S, is the total flow of stream j in kmol h-l - C is the cost scaling factor (C = -0.002204585).

Beside the normal upper and lower bounds on the optimization variables, a further inequality constraint is introduced: s3

+

s4

25

s2 s3

?isztEF

Sl

198’6Nl4 1414.ka

T~~~~K

:.

s3

(

()g5.

(15)

.

has been applied using both the SRK- and MHVZmodels as the reference rigorous model. The form of the simplified thermodynamic-model is given in equation (10) but since the pressure is constant in the column, parameter D is set equal to zero for all components. Also parameter C is set equal to zero for all components since the mixture is close to ideal. Only one set of parameters (A and B) per component was used for the entire column. The efficiency of a model simplification procedure can be measured in many ways. It is common to use the reduction in the number of rigorous thermodynamic evaluations as a measure of the efficiency. Also the ratio between the computational time for a rigorous and a simplified thermodynamic evaluation is often employed as a measure of efficiency. However, since the use of model simplification procedures also require administration and regression of the simplifled model parameters a more realistic measure of the efficiency is the reduction in overall computational time. Another important aspect is the measure of the error introduced as a result of the

Fed-S1

zi E

+

The initial values of the optimization variables together with the optimal values using both the SRKand the MHVZmodel are given in Table 5. The model simplification procedure given in Algorithm IV

. P-l.

are chosen:

Sidestream 4 to total feed ratio (S,/S, ). Sidestream 3 to total feed ratio (Z&/S,). Overhead stream 2 to total feed ratio (S,/S, ). Reflux ratio.

s,

Model-simpl$?cation problems

variables

. 54

Fig. 11. Flowsheet for the separation of hydrocarbons.

J. PERReGAARD

478

Table 5. Results of the optimiition

of the separation of hydrocarbons

optimal values SRK

SRK & SIM IO 0.237 0.140 0.315

RCAUX ws, .%I-? W%

QP

Rigorous Simp. Time (6) Objective function

22,8Z 0 50

Start

MHV

MHV & SIM 10 0.237 0.140 0.314

12 6400 7075 34

14 23.650 0 98

12 6650 6900 51

- 8.78

All 07.2 0.2 0.3 -8.27

SRK = The SRK equation of state is used in calculations involving phase equilibria. MHV = The MHV2 equation of state is used in calculations involving phase equilibria. SIM = The simplification procedure given in Algorithm IV is used. The simplified equilibrium expression is given in equation (i0) with the C and D are set equal to zero for all components. Rigorous = Number of calls to the rigorous thermodynamic model. Simp. = Number of calls to the simplified modal. Time = Time for the solution of the optimization problem in seconds. The calculations were performed on an IBM 486 machine (PSZ, Model 70). QP = Number of quadratic subproblems.

model simplification procedure. The model simplification procedure developed in this work does not introduce any additional error and the optimal solution obtained both with and without the use of the simplification procedure is identical. The overall saving in computational time when the simplification procedure is applied is 32% for the case where the SRK equation of state is used and 48% for the case where the MHV2 equation of state is used. It should be noted that only two gradient evaluations using the rigorous model were needed. This corresponds to the minimum possible, since once the optimal solution using the simplified model is obtained, the gradients are reevaluated using the rigoroous model, a line search is performed and a second evaluation of the gradients are made in order to evaluate the Kuhn-Tucker criteria. Case study 3. Isomerization of n-butane The second process in connection with optimization where the model simplification procedure is

I--

applied concerning the isomerization of n-butane. This process was simulated by Andersen (1985). A feed of n-butane polluted with some n pentane is fed to a reactor converting n-butane to isobutane. Andersen (1985) assumed that n-pentane did not participate in the reaction and the same assumption is used in this work. The isobutane product is recovered in a two-distillation column system. In the 6rst column (20 ideal stages) n-pentane is removed and in the second column (35 ideal stages) isobutane and n-butane are separated. The bottom products from both columns contain some unreacted n-butane and are therefore recycled. Five percent of the total recycle stream is purged to prevent build-up of n-pentane and n-butane. A 46% conversion of nbutane is assumed in the reactor. The process flowsheet is shown in Fig. 12. An optimization problem was formulated and the objective was to maximize the amount of isobutane in the product stream. Also a purity constraint requiring that the amount of n-butane in the product stream should be lower than

Reactor

Purge

Fig. 12. Flowsheet of the process for isomerizatlon of n-butane.

Model simplification and reduction 1% was imposed. Two optimization variables were selected. These were the distillate-to-feed ratio in the two distillation columns. The optimization problem was solved using the SRK equation of state for the calculation of equilibrium constants. A second run of the optimizaton problem including the model simplification procedure of Algorithm IV was performed and the results are given in Table 6. Two sets of A and B parameters for each column (that is four sets in total) were required in order to describe the equilibrium constants predicted by the SRK equation of state. One set for each column described the part of the column above the feed tray and one set for each column described the part below the feed tray. From the results given in Table 6 it can be seen that a reduction in the number of rigorous thermodynamic evaluations of 45% was obtained using the model simplification procedure. The reduction in the overall computing time was 20%. Still the results using the simplified procedure match the rigorous solution perfectly. Model simplzjication examples in dynamic simulation In Case studies 4 and 5, results for application of simplified thermodynamic models in dynamic simulation problems are presented. Case study 4. Distillation tower example In this example a 4-component mixture of hydrocarbons (n-butane, n -pentane, n -hexane and n-heptane) is separated in a distillation tower. Details of the process and the flowsheet are given in Fig. 13. The SRK equation of state is used for the prediction of the thermodynamic properties. All simulations were Table

6.

Results

for

the

optimization

of

the

isomer&&on

of

n-butane Optimal SRK

Rigorous Simp. Time (s) Objective function

start

SRK

Dist 1: D/F Dist 2: D/F

QP

values UCSIM

0.68 0.50 6 36,210 14,900 153

5

66,010 19: -90.37

All 0.60 0.41 -88.40

D/F = Distillate-to-feed ratio. SRK = the SRK equation of state is used in calculations involving phase equilibria. SIM = the simplification procedure given in Algorithm IV is used. The simplified expression is given in equation (IO) with the C and D set equal to zero for all components. Two sets of A and B parameters per distillation unit were used. Rigorous = number of calls to the rigorous thermodynamic model. Simp. = number of calls to the simplified thermodynamic model. Time = time for the solution of the optimization problem in seconds. The calculations were performed on an IBM 486 machine (PSZ, Model 70). QP = Number of quadratic subproblems.

479

P-latm 15

8

RdUXnb8-15OKllW

P-1.14cmn

m

Fig.

13.

Qb-6020Mjlhr

Flowsheet and specifications for the distillation tower example.

performed using the mixed-mode simulator DYNSIM (Sorensen, 1991). The simulation strategy is as follows: from the default initial conditions generated by DYNSIM dynamic simulation is performed until steady-state is obtained. At the steady-state, a step change of 4% in the reboiler heat duty is made and the corresponding transient response is simulated. The model simplification procedure was applied and the pressure and composition dependence for Ki was neglected [that is, assuming all Cj and Di in equations (9-12) are zero]. Only one set of A and B parameters per component was needed for the entire column. The computational times for the different simulation modes are given in Table 7 and it can be seen that a significant reduction in the required computational time can be obtained by using the proposed model simplification procedure. It is important to note that the proposed model simplification procedure can be applied together with quasi-Newton updating of the Jacobian and with the use of sparse methods. Hence, the model simplification procedure provides an additional reduction in computation time. In Table 7 it can be seen that applying the model simplification procedure together with a sparse method gives a further reduction in the required total computational time of 15%. It should also be noted that these reductions are compared with what is normally considered a very fast rigorous model (SRK equation of state). Also it is seen from Table 7 that for the denseand sparse-modes, the bubble point calculations take up 32 and 41% of the total simulation time, respectively. Therefore it is of course not possible to save more than these percentages when model simplification is employed. The savings obtained using the model simplification procedure with respect to the

480

“possible’*

J. PERREOAAIW savings are 59 and 40% for the dense- and

sparse-modes, respectively. Another interesting aspect to be noted from Table 7 is that the size of the saving in the computational time depends on the phase of the simulation problem, that is, the time from initialization of the problem. After a step-change or start of a new simulation an evaluation of the full Jacobian matrix is required and since the model simplification procedure saves computational time on the evaluation of the Jacobian matrix, a greater reduction in the computational time is obtained. However, due to the overhead of reinitialization of the integration method (step-size, order selection, etc.) which is not affected by the model simplification procedure, the percentage saving is smaller. The process has also been simulated using the simultaneous solution of algebraic and differential equations (DAE mode). It was found that the computational time for the DAE-mode was three times higher than for the ODE-mode (dense). This is in agreement with what has been reported earlier by Perregaard et al. (1992) where it was also shown that the efficiency of the DAE-mode is problem specific. Case study 5. Two-flash example This example was first reported by Macchietto et al. (1986). The example is a dynamic simulation of two flash units with a recycle. The mixture in the process consists of acetone, acetonitrile and water and the MHV2 model is used for phase equilibrium calculations in this work. The process is shown in Fig. 14 and the equipment parameters and feed specifications are given in Table 8. The difference between the process simulated in this work and the one by Macchietto et al. (1986) is first of all that

Macchietto er al. (1986) neglected the vapor holdup and therefore only bubble-point calculations were required and, a volume constraint and the energy balance were avoided. This assumption may be valid as a first approximation but in this work a more rigorous unit operation model including both dynamic energy balance equations and vapor holdup is used. Also Macchietto et al. (1986) calculated the K-values each time a Jacobian matrix was calculated and kept the K-values constant for the rest of the integration step. This assumption is not made here and explains why more thermodynamic property evaluations are required in this work. The dynamic response was simulated using the simulator DYNSIM. A mixture of acetone, acetonitrile and water was fed to the first flash unit, with overall composition z = (0.3,0.3,0.4). The initial holdups (moles) in the two flash units were set to (0.3,0.3,0.4). Dynamic simulation was performed for 80 h when the steady-state was reached. At this time the feed concentration was switched to a new compoand the transient resition Z* = (0.90,0.05,0.05) sponse was obtained through dynamic simulation for an additional 100 h. In order to compare with the results of Macchietto et al. (1986) the dynamic behavior of the energy balance was neglected at the time for the change in feed composition. Since, the response in the composition is one of interest, and the calculation of the enthalpy is a function of composition, temperature and pressure, either the temperature or pressure must be fixed when the energy balance dynamics are neglected. It was chosen to fix the temperature at a steady-state value using the first feed composition. The liquid phase composition in flash 1 is shown as a function of time in Fig. 15 and for flash 2 in Fig. 16. 5

F I-

4

A S H

2 I:

Fig. 14. The two-flash process.

Model

simplffrcationand reduction

481

Table 7. Simulation time for different mode of simulation for the distillation tower example (Cast study 4). Method of solution Dense Jacobian no simplification

Dense Jacobian + simplification

Sparse Jacobian no simplification

Sparse Jacobian +simpliiication

Time interval 1 (O-20 x IO-’ h) Function evaluation Jacobian evaluation Bubble point Total

9.4 9.1 3.1 13.4

6.3 6.4 1.1 8.5

Time interval 2 [2.0 x IO-‘-

1.

6G#

0.3002 h (the steady-state is reached)] 22.8 13.7 11.9 26.4

27.1 20.0 8.4 33.3

33.4 27.1 15.2 39.2

Function evaluation Jacobian evaluation Bubble point Total

3.4 2.1

4.1 3.1 1.6 7.3

17.6

8.1 2:::

Time interval 3 [0.3002-0.3382 h (at 0.3002 the step change is introduced)] Function evaluation Jacobian evaluation Bubble point Total

34.1 30.6 13.6 39.5

Function evaluation Jacobian evaluation Bubble point Total

30.1 30.4 11.3 34.6

Function evaluation Jacobian evaluation Bubble point Total

25.1 26.2 9.6 29.1

Function evaluation Jacobian evaluation Bubble point Total

132 123 53 156

26.9 22.2 6.5 32.0

17.7 11.0 9.5 24-4

13.9 7.1 6.0 21.8

9.8 8.4 4.8 12.0

7.5 5.4 2.1 9.8

9.6 8.2 4.5 11.8

7.6 5.4 2.3 9.7

Time interval 4 (0.3382-I .3382 h) 22.6 22.2 4.0 27.1 Time interval 5 (1.3382-2.3382 h) 18.8 19.1 2.4 22.9 Total simulation 102 90 22 124

64 45 32 78

50 28 tz

Time for the solution of the simulation problem &en in Case study 4 in seconds. The calculations were performed on an IBM 486 machine (PSZ. Model 701. Function evaluation = time used in the &al&on of thk right-hand sides. Jacobian evaluation = time used in the evaluation of the fuI1 Jacobian matrix. Bubble point = time used in the evaluation of bubble-points. Total = time used in the total solution of the simulation problem. Dense = a dense Jacobian matrix is us&. Sparse = a sparse Jacobian matrix is used.

Two runs were performed using the rigorous thermodynamic model in all calculations and one where the model simplification procedure was used. Since, the temperature is not changing, the B parameter was set to zero in the simplified model. One set of A, C and D parameters for each component in each flash unit was used. The transient responses in the two simulation runs were, as expected, identical. The use of the model simplification procedure resulted in a re-

duction in the overall computation time of 48%. Even though the model simplification procedure ensures that no loss of accuracy and is employed only in the evaluation of the Jacobian-matrix, this reduction in computing time is higher than the best obtained by Macchietto et al. (31%). This is due to the fact that a more rigorous unit operation model for the flash drums which required more thermodynamic property evaluations is used in this work.

Table 8. Equipment of the two-flash example Feed stream Flow

Tcmpcrature.

1 kmol h-’

297 K

Pressure. latm

Composition 0.30,0.30,0.40 (time = O-80) 0.90,0.05,0.05 (time = 80-200)

Parameters for the Bash I Vol = 0.75 m3, Rv = 0.5 kmol h-’ atm-’

Area = 0.3 m2, Q=lOOMJh-’

R, = I kmol h-’ m-‘,

Parameters for flesh 2 Vol = 0.50 ms, Rv==0.3kmolh-‘atm-‘.

Area = 0.3 m2, Q = -0.1 MJ h-’

R,=

1 kmolh-‘m-’

J.

482 1.00

PJBREGMRD

1.00

1 E

.z 0.60 & 2

ifi 0.60

Acetone _. ___. __ Acetonitril ---_ water

_______ 0.40

z

in flash in flash in flash

1 1 1

~---.___________

-

0.60 0)

5

Acetonitril Water

____

E 0.40

in flash in flash in flash

2 2 2

s

.s

3

%20 3‘

0.20

.--___________.

‘_________

0.00 0

Fig. 15. Dynamic concentration profile for flash unit 1.

CONCLUSION

An algorithm for integration with “partitioning” of the model equations has been presented. The algorithm allows the use of different integration techniques for different partitions and the efficiency and flexibility of the method is illustrated through the simulation of the esterification of acetic acid process. An algorithm for model reduction in optimization problems solved with a modular steady-state simulator has also been developed, implemented in the simulator and tested. The algorithm can be applied for optimization of flowsheets consisting of two or more partitions and is based on the observation that one or more partition does not include optimization/ design-variables. For optimization problems of this nature, it is demonstrated that a reduction in the computational time can be achieved without loss of accuracy. Also in connection with steady-state optimization an algorithm for detection of the minimum number of unit operation models required in the solution of an optimization problem is presented. The algorithm leads to a substantial reduction in the number of unit operation evaluations for the solution of an optimization problem. A general model simplification procedure has been developed. The method simplifies the IAEs of the system of equations and is applicable for problems in simulation (steady-state and dynamic), control problems, design/analysis and, for steady-state optimization. Test examples have indicated the reductions in overall computational times of the order of 2040% compared to the reference rigorous model at no loss of accuracy and with no side-effects on the convergence characteristics of the numerical method.7

zs

80

76

,a

LII

I_

170

me

Time in hr Fig. 16. Dynamic concentration profile for

flash unit 2.

Present work is applying the algorithms developed and presented in this paper on more complex chemical processes and future work will involve the development of model-reduction techniques based on the ideas of Hendirks and Van Bergen (1992).

NOMENCLATURE

A = Parameter in the simplified model B = Parameter in the simplified model C = Parameter in the simplified model

c = cost coefficient

D = Parameter in the simplified model D = Diagonal matrix containing ones for ODEs and zero for IAEs E,, = Activation energy (kJ mol - ‘) f = Vector of differential residuals F = Vector of residuals g = Vector of (implicit) algebraic residuals h = Step-size h = Vector of explicit algebraic functions J = Jacobian matrix k, = he-exponential factor for the Arrhenius expression (mole rnm3)’-” h-’ (n is the order of the reaction) K = Equilibrium ratio n = Mole number (mol) p = Vector of explicit algebraic variables P = Pressure (atm) Q = Heat supplied (MJ h-t) r = Vector of residuals in the partition integration method R = Valve resistant (kmol h-’ m-r) S = Vector of stream variables t = Time (h) T = Temperature (K) Vol = Volume (m3) x = Vector of implicit algebraic variables X = Mole fraction y = Vector of state (differential) variables Y = Vector of variables z = Vector of design variables Z = Mole fraction Greek

tNote that the computational efficiency of prediction of physical property is high and in most cases it is only between 30 and 50% of the total computation time.

PO= ;z;ent

for

the predictor+orrector

integration

P = Past information for the integration method

Model simplification and reduction

i= i= L = V =

Partition number Component number Liquid Vapor

B = Bubble-point D = Dew-point k = Iteration number k 1 = Lower u = Upper l = Approximated 0 = Reference value obj = Objective function ref = Reference value

REFERF.NCES

An&men P. M., Steady-state simulation of chemical processes. Ph.D. Thesis, Institut for Kemitelcnik, The Technical University of Denmark (1985). Andersen P. M.. F. Genovese and J. Perregaard, Manualfor the Steady-state SimJator SEPSIM. Institut for Kemiteknik, The Technical University of Denmark (1991). Biegler L. and J. E. Cuthrell, Improved infeasible path optimixation for sequential modular simulators--II. The optimization algorithm. Computers them. Engng 9, 257-267 (1985). Chimowitz E. H., T. F. Anderson, S. Macchietto and L. F. Stutxman, Local models for representing phase equilibria in multicomponent, nonideal vapor-liquid and liquidliquid systems. 1. Thermodynamic approximation functions. Iud. Emme Chem. Process Des. Dev. 22. 217-225 (1983). - Dahl S.. Aa. Fredenslund and P. Rasmussen. The MHV2 model-a UNIFAC based equation of s&e model for prediction of gas solubility and vapor-liquid equilibria at low and high pressures. Ind. Engng Chem. Res. 30, 1936-1945 (1991). and P. Rasmussen, J. Gmehling Fredenslund Aa., Vapor-Liquid E?pdIibria Using VNIFAC. Elsevier, New York (1977). Gani R., J. Perregaard and H. Johansen, Simulation strategies for design and analysis of complex chemical processes. Trans. IChemE 68, (Part A), 407-417 (1990).

483

Gear C. W., Numerical Initial Value Problems in Ordinary D@rential Equations. Prentice-Hall, Englewood Cliffs, NJ (1971). Gums E. A., Efficient use of thermodynamic-models in process calculations. Proc. 2M Int. Con3 on Foundations of Computer-Aided Process Design, pp. 249-302 (1983). Hansen P. C., Hsefte 53-Kurvetilpasning, Numerisk Institut, The Technical University of Denmark (1985). Hendirks E. M. and R. D. Van Bergen, Application of a deduction method to phase equilibria calculations. Sixth Int. Conf. Fluid Properties & Phase Equilibria for Chemical Process Design, Cortina d’Ampexxo, Italy (1992). Hillestad M., C. Sarlie, T. F. Anderson, I. Olsen and T. Hertzberg, On estimating the error of local thermodynamic-models. A general approach. Model. Ident. Control 11, 7348 (1990). Johns W. R.. and V. Vadhwans, Convergence studies in dual-level flowsheeting. Chem. Enann Res. Des. 64, -332-334 (1986). Macchietto S., E. H. Chimowitx, T. F. Anderson and L. F. Stutxman, Local models for representing phase eqiulibria in multicomponent nonideal vapor-liquid and liquid-liquid systems 3. Parameter estimation and update. Znd. Engng Chem. Process Des. Dev. 25, 674-682 (1986). Perregaard J., F. Genovese and R. Gani, Efficient integration of simulation and optimization for study of complex chemical processes. Presented at Z&d European Symp. on Use of Computers in Chemical Engineering, Barcelona, Spain (1991). Perreaaard J.. B. S. Pedersen and R. Gani, Steady-state and dynamic simulation of complex chemical _processes. Trans. IChemE 70, 99-109 (1992). Soave G., Equilibrium constants from a modified Redlich-Kwong equation of state. Chem. Engng Sci. 27, 1197-1203 (1972). Sorensen E. L., Dynamic simulation of chemical processes using accurate thermodynamic-modelling. Ph.D. Thesis, Institut for Kemiteknik. The Technical Universitv of Denmark (1991). Venkataraman S., W. C. Chan and 3. F. Boston, Reactive distillation using ASPEN PLUS. Chem. Engng Prog. 45-54 (1990). Wilson G., A modified Redlich-Kwong equation of state. Application to general physical data calculations. Paper 15c, Presented at the American Institute of Chemical Engineers Nafional Meeting, Cleveland, OH (1968).