Symbolic Incorporation of External Procedures into Process Modeling

Despite the widespread availability of sophisticated, user-friendly process modeling environments, ... proach, the user constructs the process flowshe...
0 downloads 0 Views 183KB Size
Download PDF
Ind. Eng. Chem. Res. 2002, 41, 3867-3876

3867

Symbolic Incorporation of External Procedures into Process Modeling Environments John E. Tolsma, Jerry A. Clabaugh, and Paul I. Barton* Department of Chemical Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139

Despite the widespread availability of sophisticated, user-friendly process modeling environments, the use of external procedures within these software packages will be necessary for quite some time. This paper illustrates the importance of properly handling these external procedures and describes a novel automated, symbolic approach for incorporating them correctly into an equationoriented modeling environment. 1. Introduction The value of process modeling within the chemical and biological processing industries today is indisputable. Its use can be found in an increasing number and variety of applications. For example, modeling is used in activities such as experimental design, process feasibility studies, process synthesis6 and design,21 debottlenecking and optimization,3,18 operator training,13,15 and real-time optimization. The information gained from proper modeling includes improved understanding of the process; safer, cleaner, and more profitable designs and operations; and reduced time to market for new products. For these reasons alone, process modeling, simulation, and optimization play a crucial role within the process industries. Two main approaches exist for process modeling: the sequential (and simultaneous) modular approach and the equation-oriented approach. In the modular approach, the user constructs the process flowsheet model, typically with a user-friendly graphical interface, by connecting blocks corresponding to members of a library of unit operation model subroutines. The simulator then analyzes the flowsheet structure and solves the problem with an appropriate algorithm. Substantial research has been performed over the past several decades in flowsheet analysis and solution strategies. The unit operation model libraries typically contain collections of subroutines implemented in lower-level programming languages such as C or Fortran. These subroutines take as arguments a relatively small number of parameters (e.g., flow rate, composition, temperature, and pressure of input streams; unit operation parameters; etc.) and return the variables characterizing a relatively small number of output streams. The implementation of these subroutines typically includes sophisticated tailored solution strategies to ensure convergence to correct results. Information moves through the flowsheet model in a manner similar to the flow of material through the actual process, from the output of one unit operation to the input of the next, and so on. This unidirectional flow of information from inputs to outputs makes the modular simulator ideally suited for the steady-state solution of flowsheets with relatively few recycle streams. Flow* To whom correspondence should be addressed. Tel.: 617253-6526. Fax: 617-258-5042. E-mail: [email protected].

sheets containing many recycle streams and/or design constraints (e.g., specifications on outputs) are often much more difficult to converge. In contrast to the modular approach, modern equation-oriented modeling environments attempt to solve the entire system of equations involved in the process flowsheet model simultaneously. The overall system model is typically very large (tens to hundreds of thousands of equations and variables) but extremely sparse. For example, each equation typically involves only five to ten variables. Having simultaneous access to the entire system of equations allows for the application of very sophisticated large-scale numerical integration and optimization codes to the process flowsheet model. Thus, the equation-oriented approach is more suited for dynamic calculations, such as dynamic simulation, parametric sensitivity analysis, and dynamic optimization, and for steady-state optimization. When the entire system of equations is in this open form, that is, all equations and variables are visible to the numerical algorithms, it is not typically possible to combine subsets of the model equations with specifically tailored solution algorithms as is the case in the modular approach. Consequently, equation-oriented process modeling environments tend to be far less robust for steadystate simulations than modular environments. Many equation-oriented process modeling environments provide high-level declarative input languages that allow users to develop their own models in a very natural and intuitive way, much as the modeler would with a pencil and paper using standard mathematical notation. This model prototyping capability is extremely important from the standpoint of dynamic calculations. To capture properly the dynamic behavior of a process, very detailed models must be constructed that include descriptions of vessel geometry and internal arrangement, possible flow reversals and transitions, phase changes, etc. Because it is obviously not possible to provide a library of unit operation models with sufficient detail for all processes, the availability of a flexible input language is crucial for dynamic calculations. Furthermore, most input languages support the ability to describe complex hybrid phenomena such as discrete control actions, safety interlock systems, and process disturbances. This capability enables activities such as the development of startup and shutdown procedures,19

10.1021/ie0107946 CCC: $22.00 © 2002 American Chemical Society Published on Web 05/07/2002

3868

Ind. Eng. Chem. Res., Vol. 41, No. 16, 2002

disturbance studies, operator training, and development of batch process recipes. The ability to prototype rapidly detailed models is important within many emerging fields where new and innovative processes require custom models. Some of these fields are biotechnological processes, biomedical devices, specialty polymer and other advanced materials design, microprocess systems, and modeling of biological systems (e.g., cell signaling). The model development involved in these emerging areas must be performed by experts in the field, who are often not experts in process modeling and the subsequent calculations. Consequently, the availability of sophisticated, user-friendly model prototyping tools offers several benefits such as rapid model development, model analysis tools for debugging, automated application of advanced numerical algorithms, sensitivity analysis, and visualization of results. 2. Motivation Although the input languages of modern equationoriented modeling environments enable quite general models to be constructed, limitations exist. For example, the input languages of these modeling environments are ideal for describing large-scale systems of differentialalgebraic equations (DAEs) or partial differentialalgebraic equations (PDAEs) where every equation and variable should be made accessible symbolically to the numerical engine. As alluded to above, this eliminates the possibility of harnessing sophisticated tailored solution algorithms for specific portions of the model that are difficult to converge, are expensive to compute, or exhibit multiple solutions. For example, a subset of the model equations might compute vapor-phase partial molar volumes from a cubic equation of state. A specific numerical algorithm should be applied for this subtask to ensure that the calculation always converges to the correct vapor-phase root. Also, a user might wish to apply an inside-out flash calculation to compute VLE robustly. For these reasons and others, most modern equation-oriented input languages provide the user with the ability to incorporate external procedures coded in programming languages such as C or Fortran into the overall process model. There are a number of reasons why a user would want to employ external procedures within a process model in addition to being able to incorporate custom numerical algorithms. One is to make use of third-party libraries such as physical property or chemical kinetics packages, which are usually available as subroutine or procedure libraries. Another is the fact that many organizations have vast amounts of legacy code (typically written in Fortran) that contain proprietary or classified information. Requiring the user to rewrite the equations encapsulated in this existing code in the form of the input language is tedious, error-prone, and often not possible, because of the input language limitations mentioned above. Furthermore, many of these existing codes have been well-tested and validated and are trusted. Consequently, they should be used “as is” whenever possible. When model equations are written in the input language of the modeling environment, the simulator is able to analyze these equations and construct analytical derivatives, sparsity patterns, and essentially any other symbolic information that can be exploited during

the subsequent numerical calculation. However, with current technology, external procedures are evaluated as “black boxes” within the simulator. That is, given a set of independent variables, the external procedure is called by the simulator merely to calculate values for the corresponding dependent variables. The implication of this is that, in general, necessary partial derivatives are approximated using finite differences, sparsity internal to the procedure is not exploited, and any discontinuities in the equations are not handled explicitly. Using finite differences to compute derivatives can be costly, and the values can be inaccurate. These effects might not be detrimental during dynamic simulation; however, if a parametric sensitivity analysis is performed, the use of inaccurate derivative information can result in the computation of highly inaccurate sensitivities. These problems associated with treating external procedures as black boxes are illustrated in section 5. This paper describes a novel approach in which automated code-analysis and code-transformation techniques are used to incorporate external procedures automatically and properly into an equation-oriented modeling environment. In this approach, the external source code is analyzed automatically, and new code is generated that is compiled and dynamically linked into the process simulator, providing all of the necessary symbolic information that would otherwise be neglected. In fact, all of the symbolic information that is available to the numerical algorithms for the subset of equations written in the input language of the simulator also becomes readily available for the equations corresponding to external procedures, provided source code is present. This enables the modeler to use external procedures with the confidence that subsequent calculations will be performed efficiently, robustly, and correctly. Furthermore, these techniques, which will be described below, enable the user to apply a much broader class of numerical algorithms to models embedding external codes. For example, the automated codeanalysis and code-generation techniques can be used to generate code for evaluating convex relaxations of nonlinear functions7 and interval extensions of the process model,16 allowing activities such as robust location of all roots of nonlinear systems of equations, global optimization, and nonconvex MINLP to be performed. Another advantage of the correct incorporation of external procedures is in applications where speed is crucial (e.g., online applications). Most modern process modeling environments employ an interpretive architecture for evaluating model equations and derivatives. That is, the model equations and some symbolic form of the partial derivatives are held in computer memory as data structures that are “interpreted” to provide values. It is well-known that interpreted evaluation can be as much as an order of magnitude slower than compiled evaluation. Although this does not typically imply that the overall calculation is an order of magnitude slower, the speedup can be significant if optimized, compiled external code is used, particularly in applications such as parametric sensitivity analysis, dynamic optimization, and stochastic optimization. The ideas developed in this paper have been implemented in two software packages, ABACUSS II and DAEPACK. The next two sections elaborate on the features of these software packages related to the symbolic incorporation of external procedures. They are

Ind. Eng. Chem. Res., Vol. 41, No. 16, 2002 3869

followed by example problems illustrating the importance of the proper incorporation of external procedures into an equation-oriented modeling environment. 3. Equation-Oriented Modeling Environment The equation-oriented modeling environment employed to demonstrate the ideas presented in this paper is ABACUSS II.27 Similarly to other equation-oriented modeling environments, such as gPROMS,2 Aspen Custom Modeler, DIVA,10 and ASCEND,20 ABACUSS II provides an intuitive, high-level declarative input language with which the user can describe the process model (see syntax manual in ref 27 for details). ABACUSS II also allows the user to formulate some (or all) of the process model using external subroutines or procedures. The algorithmic (or automatic) differentiation (AD) literature (e.g., ref 9) assumes that the code to be differentiated operates in the following manner. A collection of independent variable values are taken as arguments of the code, and the code evaluates the values of a collection of dependent variables as a function of these independent variables. This very general model for the operation of a code is also adopted for external procedures interfaced to ABACUSS II. If we denote the set of independent variables by x ∈ Rnx and the set of dependent variables by y ∈ Rny, a reference to an external procedure within an ABACUSS II model amounts to inserting the following ny equations into the overall process model

y ) f(x)

(1)

and the last of which is a dependent variable. Because any argument can be an open-dimension array, the actual number of equations implied by the external procedure is not fixed at this point. Once the EXTERNAL block has been introduced, it is possible to employ it to define equations within ABACUSS II MODEL blocks. For example

MODEL Example1 VARIABLE U,V,W AS NOTYPE EQUATION MULT(U,V,W) ; END introduces the equation

w ) uv

(2)

into the process model. The example above is equivalent to the model

MODEL Example1b VARIABLE U,V,W AS NOTYPE EQUATION W)U*V ; END

where f : Rnx f Rny is the function evaluated by the external code. Note that the number of equations implied by the external procedure is determined by the number of dependent variables. We will use the following simple Fortran subroutine that just multiplies two numbers and returns the result to illustrate the utility of this formulation.

where the equation is written in the ABACUSS II input language. Of course, in MODEL Example1, the equation residual is evaluated by calling compiled source code, whereas in MODEL Example1b, the equation residual is evaluated from the data structure holding the equation in memory. The ABACUSS II input

SUBROUTINE MULT(A,B,C)

VARIABLE

DOUBLE PRECISION A,B,C C)A*B

MODEL Example2 U,V AS NOTYPE EQUATION MULT(U∧2,V,0) ;

RETURN END

END

In this example, A and B are the independent variables, and C is the dependent variable. It is first necessary for the user to communicate this information in the ABACUSS II input language, which is done with the following code segment.

introduces the equation

EXTERNAL MULT(INDEPENDENT, INDEPENDENT, DEPENDENT) ; END This segment of code declares that there is an external procedure identified by MULT that has three arguments, the first two of which are independent variables

u2v ) 0

(3)

into the process model. This constraint could be written in the ABACUSS II input language as

U∧2 * V ) 0; It should be noted that any of the arguments in the reference to the external procedure can be expressed as an arbitrary function of the MODEL variables, including constants, as illustrated in Example2 above. Often, an external procedure will correspond to the residual evaluator for a system of equations, the

3870

Ind. Eng. Chem. Res., Vol. 41, No. 16, 2002

Figure 1. Layered architecture of ABACUSS II.

residual being zero at a root of the equations. For example, the subroutine

SUBROUTINE RES(NZ,T,Z,ZPRIME,DELTA) INTEGER NZ DOUBLE PRECISION T,Z(NZ),ZPRIME(NZ), DELTA(NZ) DELTA(1) ) ... etc. RETURN END would be introduced by the EXTERNAL block

EXTERNAL RES(INTEGER, INDEPENDENT, INDEPENDENT, INDEPENDENT, DEPENDENT) ; END and could be used in the MODEL block

MODEL Example3 PARAMETER NX AS INTEGER VARIABLE X AS ARRAY(NX) OF NOTYPE EQUATION RES(NX,TIME,X,$X,0(1:NX)) ; END to introduce the equations

f(t,x,x˘ ) ) 0

(4)

into the overall process model where f : R × Rnx × Rnx f Rnx. Note that the number of equations introduced is inferred from the number of dependent variables (the final argument of this external procedure). In this case, the dependent variables are the zero vector, so it is necessary to associate a dimensionality with this constant to infer the number of equations. This is denoted by the slice syntax 0(1:NX). The external subroutine referenced in ABACUSS II can be arbitrarily complex and can call any number of additional subroutines and/or functions. All that is important is that, given values for the independent variables, the code returns the corresponding values for the dependent variables. The interface to external procedures used in ABACUSS II models is completely general, provided that all arguments are native types (e.g., integers, reals, and characters) and that the external procedure can be implemented in any compiled programming language (e.g., C/C++, Fortran, and Pascal). The implementation inside the external code can be completely general. The key novelty in this paper is that, if source code is available for this subroutine, ABACUSS II will use DAEPACK to construct automatically all of the additional symbolic information it needs to perform the subsequent numerical calculation. DAEPACK is described in the following section. The remainder of this section outlines the architecture of the ABACUSS II software. The architecture of ABACUSS II was designed to enable high levels of flexibility in several areas, including how the software is interfaced, how the software is deployed and in what environment, how the process model is described, and what numerical algorithms are applied to the process model. Figure 1 contains a schematic of the three-layer architecture of ABACUSS II. The top layer is the interface level, where the user has access to the functionality of ABACUSS II to perform the desired analyses and calculations. The ABACUSS II distribution provides both a commandline interface (CLI) and a graphical user interface (GUI). In addition, ABACUSS II exports a set of documented interfaces that can be accessed in any number of third-party applications, including Microsoft Excel and

Ind. Eng. Chem. Res., Vol. 41, No. 16, 2002 3871

Matlab. The exported interfaces also allow the user to construct readily custom GUIs using Microsoft Visual Basic or Java, for example. Finally, the full functionality of ABACUSS II is available within programs written by the user in programming languages such as C, C++, or Fortran. The middle layer of ABACUSS II is the input translator and calculation executive. This layer implements the interfaces exported to the top layer. This level includes features such as input file translation, model symbolic analysis, and the calculation executive. This level also exports a number of visualization tools for examining information such as model structure and numerical results. This layer is available in several formats, which allows ABACUSS II to be deployed in a variety of environments. The software is available as a shared library (.so in UNIX/Linux and .dll in Windows) and as a C++ object for using the software locally. The software can also be executed remotely as a CORBA or DCOM component or using XML-RPC, an XML-based remote procedure call specification.28 The bottom layer provides numerical algorithms and the ability to incorporate external procedures symbolically. The external procedures incorporated can be either portions of an overall model or user-supplied numerical components. The ability to incorporate external procedures properly, the focus of this paper, is provided by the software package DAEPACK, described below. 4. Symbolic Incorporation of External Procedures The symbolic techniques used to incorporate external procedures into an equation-oriented modeling environment have been implemented in DAEPACK, a generalpurpose software library for numerical calculations.24,25 DAEPACK is divided into two main sublibraries, one containing a collection of numerical components for performing calculations such as root finding, numerical integration and parametric sensitivity analysis and the other containing a set of symbolic components that automatically construct all of the symbolic information required by the numerical components. DAEPACK provides both the standard numerical functionality of ABACUSS II and the ability to incorporate external codes into an ABACUSS II model properly. Currently, the symbolic components of DAEPACK support only Fortran source code; however, the ideas presented here can be readily extended to any procedural programming language. The symbolic components of DAEPACK have emerged from ideas originally developed in the algorithmic, or automatic, differentiation (AD) community. AD is a technique for computing exact derivative values (accurate to within roundoff error) for functions implemented in the form of codes written in some programming language. Derivative values are obtained by simply decomposing the computer code into sequences of elementary operations for which partial derivatives are known symbolically and applying the chain rule. It is important to note that AD does not furnish symbolic expressions for the derivatives; rather, it furnishes values for the partial derivatives of the dependent variables with respect to any desired values for the independent variables. How the chain rule is applied gives rise to several variants of AD, each of which is appropriate under different circumstances. A full description of AD is beyond the scope of this paper, but

excellent descriptions of AD can be found in refs 8, 9, 11, and 12, and a description of chemical engineering applications can be found in ref 22. What is important to emphasize, however, is that AD can be quite efficient, general, and automated. By automated, we mean that the user simply supplies the source code of the function evaluator and minimal additional information, such as the identities of the independent and dependent variables, and the derivative code is automatically constructed without user intervention. However, it should be noted that significant improvements can sometimes be realized by tailored applications of AD techniques. With the appropriate variant of AD, the cost of evaluating gradients and general vector-Jacobian and Jacobian-vector products is only a small multiple of the cost of evaluating the underlying function code, independent of the number of variables or functions involved. If the full Jacobian is desired, then sparsity can be exploited to reduce the cost of Jacobian evaluation in a number of ways.1,4,23,24 AD is also quite general. The underlying functions need not simply be implemented as a sequence of assignments but can be arbitrarily complex code including common blocks, loops, IF statements, and complex hierarchies of subroutines and functions. Moreover, AD is more properly termed algorithmic differentiation (see ref 9) because the code can contain complex iterative solution algorithms (e.g., computing partial molar volumes from an equation of state using an iterative algorithm). Two main approaches for generating derivative code with AD are available. The first employs the operator overloading capabilities of many modern programming languages (e.g., C++, Ada, and Fortran 90) and relies on the compiler to generate in the object code the additional instructions necessary for derivative evaluation. Other AD tools, including DAEPACK, operate as source-to-source translators where the additional instructions are added directly to the transformed source code, which can be compiled and linked into applications to provide derivative values. The underlying ideas of AD are relatively easy to understand once one realizes that any function that can be implemented on a computer can be viewed as a composite function in terms of several elementary operations, including binary operators such as + and * and intrinsic functions such as sin and exp. Derivatives of these elementary operations are usually trivial and known symbolically. An AD tool will simply identify the underlying elementary operations of a function evaluation and insert instructions into the new code corresponding to the derivatives of the elementary operations (or small composites of elementary operations). Derivatives of the original function are obtained by simply propagating the numerical values of these elementary derivatives through the code using the chain rule. The symbolic components of DAEPACK extend the sourceto-source transformation ideas of AD described above to generate from a function evaluator code a wide variety of additional codes required in the application of state-of-the-art numerical algorithms to the model. DAEPACK provides a component that constructs derivative code for evaluating gradients of scalar functions and vector-Jacobian and Jacobian-vector products of vector functions. In addition, if the model is sparse, DAEPACK can be used to generate code exploiting this sparsity to compute sparse Jacobian matrices with surprising efficiency (see section 5). DAEPACK can also

3872

Ind. Eng. Chem. Res., Vol. 41, No. 16, 2002

Figure 2. Automatic generation of code with DAEPACK.

be used to generate code that determines the sparsity pattern of the code for use by sparse linear solvers, in block triangularization, and in structural diagnosis. Similar to the DAEPACK AD component, which augments the original code with instructions for computing derivative values, the sparsity pattern component will augment the original code with instructions that simply keep track of the relationships between the user specified independent and dependent variables. If the original code contains nonsmooth intrinsic functions or IF statements, then DAEPACK can generate new code that allows these discontinuities to be handled properly during numerical integration17 and parametric sensitivity analysis.26 This component works as follows. For a given set of input variable values, the path through the conditional statements and nonsmooth functions defines a particular discrete mode, that is, a particular set of statements that are executed to compute the output variables from the input variables. DAEPACK augments the original code with instructions that record the current discrete mode in a bit vector. A single bit in the bit vector will indicate whether a particular IF statement evaluates to true or false. On subsequent calls to the locked model evaluator, the bit vector is used to define the conditional branching, rather than the actual values of the logical expressions in the conditional statements evaluated at the current values of the independent variables. Additional information, including values for the discontinuity functions, whose zero crossings are necessary conditions for discrete events, are also returned so that they can be utilized by an appropriate state event detection algorithm. DAEPACK can also be used to generate new code that evaluates the interval extension of the original code and convex relaxations of nonlinear functions in the original code. This enables the user to apply algorithms such as interval Newton/generalized bisection16 and global optimization and nonconvex MINLP.7 Figure 2 contains a schematic showing the steps performed by DAEPACK to transform the user’s original code into a set of new codes providing a wide variety of information. To use DAEPACK, the user must simply

provide the collection of Fortran source codes defining the external model and a DAEPACK specification file. The specification file contains information such as what code is to be generated, which are the independent and dependent variables, etc. The first task performed by DAEPACK in constructing new code is translation. At this step, the original source code is converted into intermediate data structures that can be analyzed and manipulated in the subsequent steps. Given information about the dependent and independent variables, a dependency analysis is performed to determine which sections of the overall code are active, that is, compute intermediate quantities that depend on independent variables and influence dependent variables. The next step is code transformation, where the manipulations required for performing the desired computation are identified. Finally, code generation is performed, during which the augmented model is written to a source file. The intermediate data structures mentioned above have been designed such that they can represent a wide variety of different procedural programming languages. Thus, to extend DAEPACK to support other programming languages, a translator must be written that maps the new target source onto the intermediate data structures, and a new code generator must be written that writes the transformed code in the desired language. ABACUSS II is able to generate automatically the specification file from the ABACUSS II input block containing the declaration of the external code. In Figure 2, six new Fortran codes are generated from the original code, providing all of the information for performing hybrid parametric sensitivity analysis.26 Figure 3 shows how DAEPACK is used in conjunction with ABACUSS II to incorporate external codes into an overall process model properly. In this example, the ABACUSS II input file contains a reference to an external subroutine EXT, for which the user has source code. Given the location of this source code and a description of the independent and dependent variables (from the declaration of EXT in the ABACUSS II input file), ABACUSS II will call DAEPACK to generate all

Ind. Eng. Chem. Res., Vol. 41, No. 16, 2002 3873

Figure 3. Symbolic incorporation of external procedures with DAEPACK.

of the code it needs to perform the subsequent calculation robustly, efficiently, and correctly. The code generated by DAEPACK will be compiled and placed into a shared library that can be dynamically loaded by ABACUSS II prior to executing the calculation. To summarize this section, all of the additional symbolic information required by ABACUSS II for proper numerical calculation (e.g., sparsity pattern, analytical derivative values, discontinuity information, etc.) can be obtained for external code with DAEPACK (provided the source code is available and written in Fortran). Thus, all equations, regardless of whether they are written in the ABACUSS II input language or are available as external code, are treated in a consistent and correct manner within ABACUSS II. 5. Examples This section contains several example problems that illustrate the importance of our novel methods for proper incorporation of external procedures into an equation-oriented modeling environment. 5.1. Hidden Discontinuities. The first example consists of three (very) small DAEs containing discontinuities. These DAE systems were coded into a Fortran subroutine and incorporated into ABACUSS II in two ways: (1) as a black box and (2) with additional symbolic information provided by DAEPACK. The first example, case A, is

x˘ ) 1 y)

{

(5)

xk1

if x e 0

xk2

otherwise

(see top diagram in Figure 4). Two numerical integrations were performed from t ) 0 to t ) 1.5 with integration tolerances (relative and absolute) set at

Table 1. Integration Statistics for Three Simple Discontinuity Examples case A

case B

case C

black black black box DAEPACK box DAEPACK box DAEPACK stepsa RESb JACc ETFd CTFe

118 215 81 19 0

110 157 69 2 0

103 177 131 33 0

79 89 74 0 0

F A I L S

71 87 63 0 0

a Steps ) number of integration steps performed. b RES ) number of model residual evaluations required. c JAC ) number of Jacobian evaluations required. d ETF ) number of integration error test failures. e CTF ) number of integration corrector test failures.

10-10, k1 ) 2, k2 ) 5, and x(0) ) -1. The integration statistics are shown in the first two columns of Table 1. As might be expected by the differentiability at the event, there is no significant difference between the two cases, although explicit handling of the event does reduce the number of error test failures. The second example, case B, is

x˘ ) 1 y)

{

(6)

xk1

if x e 0

Rx

otherwise

(see center diagram in Figure 4). Two numerical integrations were performed from t ) 0 to t ) 1.5 with integration tolerances (relative and absolute) set at 10-10, k1 ) 2, R ) 1000, and x(0) ) -1. The integration statistics for this example are shown in the second two columns of Table 1. In this example, there is a significant difference between the case where the external procedure is incorporated as a black box and when it is not. The black-box mode required nearly twice as many residual evaluations, Jacobian evaluations, and LU factorizations. If these equations were part of a much

3874

Ind. Eng. Chem. Res., Vol. 41, No. 16, 2002 Table 2. Timing Information for Constant-Pressure Batch Reactor Examplea timing (s) black box

DAEPACK

454.8

20.6

overall simulation time a

Calculations performed on a 1.4-GHz PC with 512 MB of RAM.

quently, they should be used with minimal modification whenever possible. The first example is a model of a perfectly stirred batch reactor operating adiabatically and at constant pressure. The model consists of the following Nc + 2 equations

F

dyi ) Wiwi dt FCp

dT

i ) 1, ..., Nc Nc

)

dt

∑ Wkhkwk

larger overall model, then this additional work could be quite significant. The third example, case C, is

x˘ ) 1 y)

{

(7)

xk1

if x e 0

R + xk2

otherwise

(see bottom diagram in Figure 4). Two numerical integrations were performed from t ) 0 to t ) 1.5 with integration tolerances (relative and absolute) set at 10-10, k1 ) 2, k2 ) 2, R ) 1, and x(0) ) -1. In this example, the numerical integration fails at the event when the external procedure is treated as a black box. If these equations were part of a much larger overall model, they would probably also cause a numerical integration failure, and it would be very difficult to isolate the cause of the failure. Although these three simple DAEs are small, they illustrate the importance of properly handling external procedures containing discontinuities. 5.2. Correct Use of Chemical Kinetics Libraries. The remaining two examples illustrate the importance of properly handling external procedures for physical property and chemical kinetics calculations. These computations are typically very complex, costly, and difficult to code correctly and efficiently. In addition, many existing high-quality codes are available that have been extensively validated and are trusted. Conse-

(9)

k)1

F ) F(T,y)

Figure 4. Simple discontinuities: example A (top), example B (middle), example C (bottom).

(8)

(10)

where F is the mass density, yi is the mass fraction of component i, Nc is the number of chemical species, Wi is the molecular weight of species i, T is the temperature, Cp is the constant-pressure heat capacity, wi is the molar production rate of species i per unit volume, and hi is the specific enthalpy of species i. The molar production rates, heat capacity, mass density, molecular weights, and enthalpies were computed with external Fortran subroutines from the CHEMKIN-II library.14 The chemical mechanism for the reaction of oxygen, nitrogen, and n-heptane used in this example involves 544 chemical species (i.e., Nc ) 544) and 2446 reactions. This mechanism was obtained from Curran et al.5 Two simulations were performed with this model, one in which the external procedures were treated as black boxes and the other in which they were incorporated properly using DAEPACK. The numerical integration was performed for a simulation time of 5.0 s on a 1.4GHz PC with 512 MB of RAM. The initial conditions were yO2 ) 0.0252, yN2 ) 0.9734, yC7 ) 0.0014, yO ) yH ) 10-16, and T ) 800 K. Table 2 contains the timing information for this example. Clearly, significant improvements can be realized by exploiting the additional symbolic information available when the external procedures are incorporated properly. This benefit is examined in more detail in the following example. The next example is a reacting flow simulation. A gaseous mixture of oxygen, nitrogen, and n-heptane is injected into a tubular reactor. In this model, the reactor is assumed to be isothermal and isobaric, and the gas is assumed to be ideal. Also, it is assumed that there are only variations in time and the axial direction and that diffusion is negligible. These assumptions result in the following system of equations

∂xi ∂(xivz) RT + w )0 ∂t ∂z P i ∂vz ∂z

-

RT P

i ) 1, ..., Nc - 1

(11)

Nc

∑ wk ) 0

(12)

k)1

Nc

∑ xk ) 1

k)1

(13)

Ind. Eng. Chem. Res., Vol. 41, No. 16, 2002 3875 Table 3. Timing Information for Unsteady PFR Examplea timing (s) overall simulation time single Jacobian evaluation time single LU factorization time single LU backsubstitution time a

Figure 5. Sparsity pattern of the unsteady PFR model.

where xi(t,z) is the mole fraction of species i, Nc is the number of species, t is the time, z is the axial coordinate, R is the gas constant, T is the temperature, P is the pressure, vz(t,z) is the gas velocity in the z direction (all other velocity components are assumed to be negligible), and wi(T,P,x(t,z)) is the molar production rate of species i per unit volume. The PDAE above was discretized using upwind finite differences and coded into an ABACUSS II input file. As in the previous example, the Nc , were computed with molar production rates, {wk}k)1 external Fortran subroutines from the CHEMKIN-II library, and the same n-heptane mechanism was used. Using 10 grid points in the discretization, the overall model consisted of 10 890 variables (Nc mole fractions, Nc molar production rates, and one velocity on each of the 10 grid points) and equations. Figure 5 contains a diagram of the exact sparsity pattern for the unsteady PFR model. The ordering of variables in this sparsity pattern is the molar production rates on grid point 1, followed by the molar production rates on grid point 2, and so on. The molar production rates are followed by the mole fractions on grid point 1, followed by the mole fractions on grid point 2, and so on. The velocities on each grid point are the last 10 columns of the sparsity pattern. The order of the equations (rows of the sparsity pattern) are the calls to the external Fortran code to compute the molar production rates on each grid point, followed by the discretization of the species balance equations, followed by the discretization of the velocity relationships, followed last by the summation of mole fraction constraints on each grid point. The sparsity pattern for the portions of the overall model corresponding to external procedures was obtained with DAEPACK, and these blocks can be seen in the upper right corner of the sparsity pattern in Figure 5. Note that each of these subblocks of the Jacobian matrix is 544 rows by 544 columns; however, each one contains only 12 518 nonzero entries. Although not clearly evident because of the limited resolution of the figure, these blocks are actually approximately 96% sparse. When these external procedures are treated as black boxes, these 10 blocks of the overall Jacobian matrix would contain 295 936 entries each. Two numerical integrations were performed with this model using ABACUSS II, one in which the external

black box

DAEPACK

8252 160 2.0 0.12

280 3.3 0.1 0.01

Calculations performed on a 1.4-GHz PC with 512 MB of RAM.

procedure was treated as a black box and another in which symbolic information was obtained with DAEPACK. Both simulations were performed at an absolute pressure of 12.5 atm and a temperature of 900 K. The reactor length was 5 cm. The initial condition was the tubular reactor filled with pure nitrogen. The boundary condition was a fixed composition and velocity at the inlet of the reactor (mole fractions of O2, N2, and n-heptane equal to 0.0252, 0.9734, and 0.0014, respectively, with trace quantities of oxygen and hydrogen free radicals at mole fractions of 10-16 each). The inlet velocity was fixed at 1 cm/s. The numerical integration was performed for a simulated time of 0.1 s on a 1.4GHz PC with 512 MB of RAM. Table 3 contains the timing information for this example. This example clearly illustrates the benefit of properly incorporating external procedures during numerical calculations. Using this additional symbolic information, obtained automatically by ABACUSS II with DAEPACK, the simulation time was reduced from approximately 2.3 h to 4.7 min, a 30-fold speed improvement. The performance improvement can essentially be attributed to two factors. First, by efficient accumulation of (structurally) nonzero derivative values in the sparse blocks of the Jacobian matrix corresponding to the external code, the overall Jacobian matrix evaluation time was reduced from 160 to 3.3 s. In both scenarios, all derivative values other than those associated with the external code were computed in the same manner. Second, by exploiting sparsity in the linear solver, the time for a single LU factorization was reduced from 2 to 0.1 s, and the time for a single backsubstitution was reduced from 0.12 to 0.01 s. Again note that, in the black-box example, it is only the subblocks associated with the external code that are dense and sparsity is exploited for all other portions of the Jacobian matrix in both scenarios. This example highlights an interesting observation regarding simulations involving complex physical properties or kinetic mechanism calculations. In many numerical integration calculations, it is the cost of the LU factorization that tends to dominate the cost of the overall calculation (which is why most modern numerical integration codes attempt to reduce the number of LU factorizations performed). However, if the residual and Jacobian evaluations are costly, as is the case in the complex chemistry example above, the Jacobian evaluation time can significantly exceed the cost of the linear algebra, even in quite large-scale models. Thus, the proper incorporation of external procedures can substantially reduce the cost of the calculation even if matrix-free linear solvers are applied. Also of significance is the amount of memory saved by exploiting sparsity. For example, when the external procedures were treated as black boxes, each dense block of the overall Jacobian matrix corresponding to these equations contained 295 936 entries. Because the matrix is stored in sparse triplet form (i.e., a double-precision array containing the values of the Jacobian matrix and

3876

Ind. Eng. Chem. Res., Vol. 41, No. 16, 2002

two integer arrays containing the row and column indices), the external blocks alone require approximately 47 MB to store the derivative information. This is compared to only 2 MB of storage required for the external procedure blocks when sparsity is exploited. 6. Conclusion This paper describes how source-to-source codetransformation techniques can be used to incorporate external code into an equation-oriented process modeling environment properly. This enables the user to write complex models described partly in the input language of the process simulator and partly with new or legacy external codes. By properly handling these external procedures, the user can be confident that the subsequent calculation will be performed robustly, efficiently, and correctly. The ideas described in this paper have been implemented with the equation-oriented process simulator ABACUSS II and the numerical and symbolic software library DAEPACK. DAEPACK currently works with Fortran but can be readily extended to other procedural programming languages. Comparable capabilities can be achieved with object-oriented languages such as C++ using operator overloading features. Although the focus of this paper is on incorporating external procedures into equation-oriented modeling environments, the techniques described are quite useful for incorporating external procedures into modular simulators for steady-state simulation and optimization. In particular, the ability to generate fast and accurate analytical derivative values can often substantially improve the performance of steady-state calculations, in particular, steady-state optimization. Acknowledgment The authors acknowledge support from the EPA Center for Airborne Organics at MIT and Mitsubishi Chemical Corporation. Literature Cited (1) Averick, B. M.; More´, J. J.; Bischof, C. H.; Carle, A.; Griewank, A. Computing large sparse Jacobian matrices using automatic differentiation. SIAM J. Sci. Stat. Comput. 1994, 15, 285-294. (2) Barton, P. I. The modeling and simulation of combined discrete/continuous processes. PhD Dissertation, University of London, London, U.K., May 1992. (3) Biegler, L. T.; Hughes, R. R. Process optimization: A comparative case study. Comput. Chem. Eng. 1983, 7, 645-661. (4) Bischof, C. H.; Khademi, P.; Bouaricha, A.; Carle, A. Efficient computation of gradients and Jacobians by transparent exploitation of sparsity in automatic differentiation. Optim. Methods Software 1996, 7, 1-39. (5) Curran, H. J.; Gaffuri, P.; Pitz, W. J.; Westbrook, C. K. A comprehensive modeling study of n-heptane oxidation. Combust. Flame 1988, 114, 149-177. (6) Douglas, J. M. Conceptual Design of Chemical Processes; McGraw-Hill: New York, 1988. (7) Gatzke, E. P.; Tolsma, J. E.; Barton, P. I. Construction of convex function relaxations using automated code generation techniques. Optim. Eng. 2002, in press.

(8) Griewank, A. On automatic differentiation. In Mathematical Programming: Recent Developments and Applications; Iri, M., Tanabe, K., Eds.; Kluwer Academic Publishers: Dordrecht, The Netherlands, 1989; pp 83-108. (9) Griewank, A. Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation; SIAM: Philadelphia, PA, 2000. (10) Holl, P. W.; Marquardt, W.; Gilles, E. D. DIVAsA powerful tool for dynamic process simulation. Comput. Chem. Eng. 1988, 12, 421. (11) Iri, M.; Tsuchiya, T.; Hoshi, M. Automatic computation of partial derivatives and rounding error estimates with applications to large-scale systems of nonlinear equations. J. Comput. Appl. Math. 1988, 24, 365-392; Original Japanese version appeared in J. Inf. Process. 1985, 26, 1411-1420. (12) Iri, M. History of automatic differentiation and rounding estimation. In Automatic Differentiation of Algorithms: Theory, Implementation, and Application; Griewank, A., Corliss, G. F., Eds.; SIAM: Philadelphia, PA, 1991; pp 1-16. (13) Kassianides, S. C. An integrated system for computer based training of process operators. PhD Dissertation, University of London, London, U.K., 1991. (14) Kee, R. J.; Rupley, F. M.; Miller, J. A. CHEMKIN-II: A FORTRAN Chemical Kinetics Package for the Analysis of GasPhase Chemical Kinetics; Technical Report SAND89-8009; Sandia National Laboratory: Albuquerque, NM, 1980. (15) Mani, S.; Shoor, S. K.; Pederson, H. S. Experience with simulator training for ammonia plant operators. Plant/Oper. Prog. 1990, 10, 6-10. (16) Moore, R. E. Methods and Applications of Interval Analysis; SIAM: Philadelphia, PA, 1979. (17) Park, T.; Barton, P. I. State event location in differential algebraic models. ACM Trans. Modell. Comput. Simul. 1996, 6, 137-165. (18) Parker, A. L.; Hughes, R. R. Approximation programming of chemical processes. 1: Optimization of FLOWTRAN models. Comput. Chem. Eng. 1981, 5, 123-133. (19) Perris, F. A. The growing importance of dynamic simulation for process engineers; In Dynamic Simulation in the Process Industries; UMIST: Manchester, U.K., 1990. (20) Piela, P. C. ASCEND: An object-oriented computer environment for modeling and analysis. PhD Dissertation, Carnegie Mellon University, Pittsburgh, PA, 1989. (21) Seader, J. D. Computer Modeling of Chemical Processes; AIChE Symposium Series 81; American Institute of Chemical Engineers: New York, 1985. (22) Tolsma, J. E.; Barton, P. I. On computational differentiation. Comput. Chem. Eng. 1998, 22, 475-490. (23) Tolsma, J. E.; Barton, P. I. Efficient calculation of sparse Jacobians. SIAM J. Sci. Comput. 1999, 20, 2282-2296. (24) Tolsma, J. E.; Barton, P. I. DAEPACK: A Combined Symbolic and Numeric Library for General Numerical Calculations; Massachusetts Institute of Technology: Cambridge, MA, 2000; http://yoric.mit.edu/daepack/daepack.html. (25) Tolsma, J. E.; Barton, P. I. DAEPACK: An open modeling environment for legacy models. Ind. Eng. Chem. Res. 2000, 39, 1826-1839. (26) Tolsma, J. E.; Barton, P. I. Hidden discontinuities and parametric sensitivity calculations. SIAM J. Sci. Comput. 2002, 23 (6), 1862-1875. (27) Tolsma, J. E.; Clabaugh, J.; Barton, P. I. ABACUSS II: Advanced Modeling Environment and Embedded Process Simulator; Massachusetts Institute of Technology: Cambridge, MA, 2000; http://yoric.mit.edu/abacuss2/abacuss2.html. (28) XML-RPC Web Page; UserLand Software, Inc.: Burlingame, CA, 2001; http://www.xmlrpc.com.

Received for review September 25, 2001 Revised manuscript received March 7, 2002 Accepted March 7, 2002 IE0107946