Equation-oriented dynamic simulation current status and future perspectives

EQUATION-ORIENTED DYNAMIC SIMULATION CURRENT STATUS AND FUTURE PERSPECfIVES C.c. Pantelides and P.I. Bartont Centre for Process Systems Engineering, Imperial College of Science, Technology and Medicine, London SW7 2BY, United Kingdom

ABSTRACf Some of the key issues associated with the equation-oriented approach to process modelling and simulation are examined. It is argued that most processes and operations of industrial interest have combined discrete and continuous characteristics. Furthermore, effective tools for describing both the physical behaviour of process systems, and the external control actions and disturbances imposed on them are required. In both areas, major concerns include managing the complexity of real processes, and supporting the development of general components that can be re-used within a number of different applications. Recent work in this direction is reviewed. A general mathematical statement for the combined discrete and continuous dynamic simulation problem is presented, and the use of symbolic, structural and numerical solution techniques is discussed. Methods for the diagnosis of problems that are either badly-posed or difficult to solve with currently available methods are also considered.

The areas of modelling distributed parameter systems, and supporting the formal description and solution of general optimisation problems are identified as being of strategic importance for the future evolution of general-purpose process modelling environments.

KEYWORDS Dynamic simulation; combined simulation; process modelling; numerical methods; differentialalgebraic equations.

INTRODUCTION Dynamic simulation is well recognised as a valuable tool at all stages of plant design and operation (see, e.g. Perkins and Barton, 1987, and Naess et al., 1992). Traditional applications in plant operations include the study of the response of continuous plants to disturbances causing deviations from steady-state, and the tuning of controllers. More recently, increasing safety and environmental concerns and stricter regulation have provided strong incentives for the study of plant operation during abnormal conditions (often involving significant excursions from the nominal steady-state), for performing hazard and operability studies, and for validating automatic emergency control procedures. The use of dynamic simulation for training plant operators to deal with such relatively rare situations, as well as to carry out their normal duties, is also well documented (Mani et al., 1990, Macchietto and Kassianides, 1991). In addition to serving as an aid to plant operation, dynamic simulation is increasingly gaining a role as a design tool. In order to deal effectively with the growing complexity of process plants, many applications of dynamic simulation such as the analysis of intrinsic controllability of the process, the t Current address: Department of Chemical Engineering, Massachusetts Institute of Technology. Cambridge. MA 02139, U.S.A.

S263

S264

European Symposium on Computer Aided Process Enginc:c:ring-2

design of the control system, and the study of appropriate start-up and shut-down procedures now form part of the basic design process. Another factor leading towards the increased use of dynamic simulation is the re-emergence of batch plants as a convenient vehicle for flexible manufacturing in the chemical industry, in particular for the production of diverse small-volume but high-value added chemicals and foodstuffs. Such processes always operate in the transient regime, and dynamic simulation is required to study their behaviour in detail. The evolution cf dynamic simulation tools mirrored to a large extent the earlier one of their steadystate counterparts, starting with the construction of programs tailored to simulate specific operations, followed first by the development of generic simulators for individual unit operations, and then by the emergence of general-purpose dynamic simulation packages for entire continuous processes. Commercial examples of the latter type of package, such as GEPURS (Shinohara, 1987) and SpeedUp (Aspen Technology, 1991), are now available and are increasingly being used in industrial practice. The development of general-purpose dynamic simulation packages raises a number of issues. One of them is the need for extensive support for model construction activities. In contrast to general-purpose steady-state packages for which it has been possible to build comprehensive libraries containing a relatively small number of models, experience with dynamic simulation so far indicates that different applications may require significantly different models of ostensibly similar unit operations. This may be because, in fact, the unit operations behave differently in the dynamic domain despite their identical steady-state behaviour. For instance, a mixing tank with holdup requires a different model to a pipe junction, although they have the same steady-state functionality. In other cases, the difference in the models may be dictated by the different nature of the application, e.g, much more rigorous and detailed models are often required for safety studies involving significant depanures from normal operation than for studying the regulatory control of a process about its steady-state. For all these reasons, a "comprehensive" library of dynamic unit operations models is probably an unattainable goal. We therefore believe that the provision of powerful model development facilities is an essential requirement for dynamic simulation packages. In attempting to satisfy the above requirement, equation-oriented modelling packages separate the description of the physical system from the mathematical solution method. The former is then provided by the model developer, while the system takes responsibility for the latter. It is now well recognised that a further advantage of the separation of model definition from model solution is the scope that it allows for the development of multi-purpose models. Effectively, through the use of appropriate solution methods, the same model can be used for a number of different tasks, such as steady-state simulation and design, dynamic simulation, steady-state and dynamic optimisation, data reconciliation etc. Overall, this feature has resulted in a shift away from purely dynamic simulation packages towards general-purpose modelling environments, in which a common process model is used in a number of different modes. The benefits in tenus of avoiding duplication of model developmcnt effon and guaranteeing the consistency of the results obtained are quite evident. The equation-oriented approach is not new: it was clearly described, for instance, in early discussions of computer-aided design by Sargent (1967), its application to steady-state flowsheeting was reviewed by Shacham et at. (1982) a decade ago, and a detailed analysis of its merits was presented by Perkins (1983). Yet, it is only in the last few years that general-purpose process modelling tools based on this approach have reached the stage of being used widely in industrial practice. In retrospect, much of this delay can be attributed to the lack of, on one hand, sufficient computing power, and, on the other, reliable algorithms and software for the solution of the underlying mathematical problems. Most current process modelling packages (equation-oriented or otherwise) are mainly concerned with the description of predominantly continuous processes (see, e.g., the recent review by Marquardt, 1991). However, a large proportion of processes of industrial interest also involve complex discrete manipulations. Examples include batch operations, periodic separation processes, the start-up and shut-down of (predominantly) continuous processes, and the application of digital control strategies. Specialised packages for modelling and simulating some of the above types of processes (e.g. batch operations) have been described in the literature (Joglekar and Reklaitis, 1984, Czulek 1988). While these are certainly useful within their specific domain of interest, their flexibility in terms of supporting general model development is rather limited, especially when compared with the capabilities of state-of-the-an modelling packages for continuous processes. Overall, it is becoming increasingly apparent that a systematic and general framework for describing processes with combined discrete and continuous characteristics is required. In the particular context of dynamic simulation, it must be recognised that the description of the complex sets of external

European Symposium on Computer Aided Process Engineering-2

S265

actions and disturbances imposed on a plant is no less important than the modelling of the physical system characteristics. The next two sections of this paper are concerned with modelling, first of the physical behaviour of process systems, and then of the external actions that may be imposed on them. In both cases, particular emphasis is placed on effective mechanisms for handling modelling complexity, and facilitating the development of re-usable components. This is followed by a brief consideration of some of the issues associated with building complete dynamic simulation descriptions from the combination of models of physical systems and external actions. The following section is concerned with the definition and solution of the mathematical problems posed by combined discrete and continuous dynamic simulation. The paper then proceeds to discuss techniques for problem diagnosis, an important issue on which relatively little specific information has been published to date. We conclude by identifying some areas that we consider to be of considerable significance for the future evolution of equation-oriented modelling packages.

MODELLING PHYSICAL BEHAVIOUR Much attention has been focused in recent years on tools for the formal definition of models describing the physical behaviour of process systems. In reviewing the current state of the art, this section concentrates on three aspects of such model descriptions. First, formalisms for the construction of models for simple unit operations are considered. This is followed by a discussion of approaches for the description of discontinuities in such models. Finally, mechanisms for managing model complexity are presented. At the lowest level, process models are often described by sets of variables and the equations that relate them. The latter typically represent conservation laws, equilibrium conditions, thermodynamic relations etc. This is the approach adopted by SpeedUp (perkins and Sargent, 1982, Pantelides, 1988a), ASCEND (Piela, 1989, Piela et al., 1991), and gPROMS (Barton and Pantelides, 1991). All these packages provide high level declarative languages for describing models in a fairly natural fashion. Alternatively, model equations may be coded in a procedural computer programming language such as FORTRAN. An example of this approach is provided by the DIVA package (Holl et al., 1988, Kroner et al., 1990). The correct, complete and non-redundant definition of models as sets of equations is not always easy. This realisation has recently led several researchers to investigate the definition of models in terms of high-level descriptions of the physicochemical mechanisms that take place in each unit operation (see Henning et al., 1990, Vazquez-Roman, 1992). This approach potentially provides a much more convenient and powerful means for model building as the mathematical model construction is now completely automated. Although in their simplest form models for unit operations are described in terms of continuous equations, many such models also involve one or more discontinuities. These typically arise from thermodynamic (e.g. phase) or flow (e.g. from laminar to turbulent regime) transitions, or from irregularities in the geometry of process vessels (e.g. non-uniform cross-sections. overflow pipes or weirs). Effective representation mechanisms for these model discontinuities are clearly very important. One such mechanism is the IF-THEN-ELSE construct of SpeedUp described by Pamelides (l988a):

IF Logical Condition THEN Equation-I ELSE

Equation-Z ENDIF where both the logical condition and the equations in each clause are expressed in terms of the model variables and/or their time derivatives. Such a construct defines two system states corresponding to the IF and the ELSE clauses respectively. Discontinuities involving multiple states may be described by nesting severaIIF-THEN-ELSE statements within each other, or using equivalent constructs, such as the CASE statement of ASCEND (Piela, 1989).

5266

European Symposium on Computer Aided Process Engineering-2

It is worth emphasising that, although the above constructs are syntactically similar to the corresponding statements of procedural programming languages (e.g. FORTRAN), they represent discontinuous equations that are purely declarative in nature - just like all other equations in the model. However, only a subset of the discontinuities that typically occur in process systems can be represented in this fashion. In particular, these are discontinuities in which transitions between the two system states are, first, allowed in both directions, and, second, triggerred by the same logical condition switching from TRUE to FALSE and vice versa. We call these reversible and symmetric discontinuities. An example is the discontinuity describing the flow over an overflow weir in a process vessel: transitions between the two system states, "Flow" and "No Flow", can occur in both directions, as determined by whether the single logical condition

Liquid Level < Weir Height is true or false. On the other hand, many process discontinuities are neither symmetric nor reversible. An example of an irreversible discontinuity is that of a vessel fitted with a bursting disc. The transition between the two disc states Intact and Burst is triggered by the condition

Pressure

~

Bursting Pressure

and is irreversible: once burst, the disc will not be restored to its former state irrespective of the value of the vessel pressure. Finally, some discontinuities are reversible but asymmetric. These typically occur in systems with hysteresis, such as safety valves that open and close at different pressures. I. Tank with a weir

Level > We ir_Heighl

NOT Level > Weir_Heighl 2 Vessel filled with a bursting disc

3. Vessel filled with a safety relief valve Press > SCI_Press

Fig. 1. General Representation of Discontinuities in Process Systems All of the above types of discontinuities can be represented as finite automata, as shown in Fig. I. An appropriate language construct for defining such discontinuous equations in their most general form is included in the model description of the gPROMS package (Barton and Pantelides, 1991, Barton, 1992). An example of the model of a safety valve expressed in this language is shown in Fig. 2. It should be noted that the system states are enumerated explicitly by the special SELECTOR variable, Valve_State. A CASE construct, similar to that used by Cellier and Bongulielmi (1980) in their COSY language, is used to define both the appropriate modelling equations in each state, and the logical conditions for transitions between states. Clearly the constructs used by earlier packages for the description of reversible and symmetric discontinuities are only special cases of this CASE construct. However, they may be retained as convenient short-hand for use when appropriate. The declarative model language enhanced by powerful mechanisms for discontinuity representation provides an adequate tool for the description of models of simple unit operations. Models for more complex systems are also characterised by sets of variables and equations, and could, in principle, be

European Symposium on Computer Aided Process Engineering-2

5267

MODEL Vessel_With_Safety_Valve PARAMETER Vessel Volume, R SecPress, ReseatPress Valve_Const

AS REAL AS REAL AS REAL

VARIABLE Flow In, Flow_Out, RelieT Flow Holdup Temp Press, Press_In

AS Molar Flowrate AS MolesAS Temperature AS Pressure

STREAM Inlet Outlet Relief

AS MainStream AS MainStream AS MainStream

: Flow_In, Press_In : Flow Out, Press : RelieT_Flow, Press

SELECTOR Valve_State

AS (Closed, Open)

EQUATION $Holdup :::; Flow_In - Flow_Out - ReliefFlow ; Press"Vessel_Volume e R*Holdup*Temp ; CASE Valve State OF WHEN Closed : WHEN Open

END

Relief Flow =0 ; SWITCH TO Open IF Press >= Set Press ; RelieCFlow = Valve_Const*Press/SQRT(Temp) ; SWITCHTO Closed IF Press <= Reseat Press ;

END

Fig. 2.

Example of Representing General Model Discontinuities in gPROMS

described in a similar way. However, the efficient and error-free construction of such a description would be an overwhelming task in most cases, so appropriate decomposition mechanisms are necessary for handling this additional complexity. One such mechanism is provided by the traditional process flowsheet represent ation as a network of interconnected unit operations . Most equation-oriented process simulation packages support such a facility through "stream" connections used to build higher level models from lower level ones. This can be done by recognising a fixed maximum number of modelling levels, such as the MODEL/MACRO/FWWSHEET hierarchy of SpeedUp. Alternatively, the concept of hierarchical submodel decomposition (Elmquist, 1978) can be introduced to allow the formal definition of models in a recursive manner, thus permitting an effectively unlimited number of hierarchical levels. The latter approach has been used in OMOLA (Nilsson, 1989, Andersson, 1990), ASCEND. and gPROMS . Furthermore, general operators for declaring that two subsets of the attributes of two separate models are, in fact, identical can be used as alternatives to stream connections (cf, the ARE_THE_SAME operator of ASCEND). An example of a hierarchical definition of a plant model is shown schematically in Fig.3, and a declaration of one level of this hierarchy is illustrated in the gPROMS language in Fig. 4. The implementation of hierarchies of unlimited depth is facilitated by effectively making no syntactical distinction between models at different hierarchical levels. In particular, any model may contain sub-models declared as instances of other models defined earlier (UN/D. The stream connections in any model (STREAM) may either be identified with existing stream connections of its sub-models (as in Fig. 4), or be declared as subsets of the variables in the model and/or its sub-models (as in Fig. 2). Stream connections between sub-models are treated as equations within the model (EQUATION) .

S268

European Symposium on Computer Aided Process Engineering-2

nAY

I

i

I

i

~t '" ~ J., l Fig. 3.

Hierarchical Submodel Decomposition

MODEL Distillation_Column PARAMElER No_Trays, Feed_Position

AS INlEGER

UNIT Condenser Top_Section, Bottom_Section Feed Reboiler

AS AS AS AS

Total Condenser Linked_Trays Feed3ray Partial_Reboiler

VARIABLE Net_Energy-Requirement STREAM Feed Stream Top_Product Bottom Product SET Top_Section.NoTrays Bottom_Section.NoTrays

IS Feed.Feed Stream IS Condenser~Liquid Product IS Reboiler.LiquidProduct := NoTrays - Feed Position ; := Feed_Position - 1 ;

EQUATION # Define the net energy requirement for the column Net_Energy-Requirement = Reboiler.Heat Load + Condenser.Hcat_Load ;

# Continuous connections between sub-models Condenser.Reflux IS Top_Section.Liquid_In Condenser.Vapour_Feed Top_Section.Vapour_In Top_Section.Liquid_Out

Feed.VapourIn

Feed.Liquid_Out

Bouom Section.VapcurIn BottomSection.LiquidOut

IS Top_Section.Vapour_Out IS Feed.Vapour_Out IS Feed.Liquid_In IS BOllom_Section.Vapour_Out IS BOllom_Section.Liquid_In IS Reboiler.VapourOut IS Reboiler.Liquid_Feed

END

Fig. 4.

Example of Higher-level Model Definition in gPROMS

A different mechanism for building model hierarchies is through the use of inheritance, aconcept popularised by the object-oriented programming languages (Meyer, 1988). Although much has been

European Symposium on Computer Aided Process Engineering-2

S269

published in this area in the literature, the main idea of relevance to the present paper is that a detailed model can be constructed from a simpler one by "inheriting" and "refining" it (i.e. adding more detail to it). By repeating this type of operation, a hierarchy of models organised in a tree structure can be built, as shown schematically and syntactically in Fig. 5 and Fig. 6 respectively'. The model for a second-order CSTR in Fig. 6 inherits that of a generic CSTR which it then refines by adding an equation defining the reaction rate. All identifiers not declared explicitly in this model are inherited from its parent.

TAo"'"

Fig. 5.

Model Inheritance

MODEL SecondOrderCstr INHERITS Cstr PARAMElER Key_Component AS INlEGER VARIABLE Rate_Const EQUATION Reaction_Rate = Rate_Const*Concentration(Key_Component)"2 ; END Fig. 6.

Example of Model Inheritance in gPROMS

Whilst the attractions of model hierarchies built through the combination of models in flowsheets are quite obvious, those of constructing model hierarchies through inheritance are rather less so. Inheritance certainly obviates the need for specifying the same information in a number of similar models, thus reducing the possibility of error. Also changes in a lower level model automatically propagate throughout all its descendants in the hierarchy. However, beyond these practically useful but rather mundane functions, inheritance serves as a mechanism for grouping diverse models into sets, each member of which is guaranteed to have a minimum set of attributes (variables, equations, streams etc.), namely those of the root model in the hierarchy. This establishes a concept of model type (or "class") - an important contribution given the complexity of model entities. One application of such model types will be demonstrated in the next section.

t In addition to the single inheritance illustrated in Fig. 5, one can envisage multiple inheritance hierarchies, in which each model is allowed to inherit attributes from more than one parent model, and also selective inheritance, in which not all the attributes of the parent models are inherited. Both may pose problems with respect to ensuring and maintaining model consistency (Meyer, 1988).

S270

European Symposium on Computer Aided Process Engineering-2

MODELLING EXTERNAL ACfIONS The accurate analysis of the effects of external control actions and disturbances imposed on the physical system is one of the prime motivations for developing dynamic models and using dynamic simulation. However, in contrast to the attention lavished on modelling the physical behaviour of process systems, much less emphasis has been devoted to providing effective representations for such actions and disturbances, in particular when these are of a discrete nature. For instance, the SpeedUp language (Pantelides, 1988a) only allows the declaration of discrete changes to the system inputs. This modest capability is further restricted by the requirement that the time of the occurrence of these changes must be specified explicitly. The facilities provided by most other general-purpose process simulation packages (at least as documented in the literature) are no more powerful. Yet even relatively simple simulations may involve much more sophisticated discrete manipulations than those described above. For instance. although many discrete control actions do indeed correspond to changes in the values of the input variables, they are often triggered by logical conditions expressed in terms of the values of the system unknowns. Consequently, the precise timing of such events, or even whether or not they ever occur during a given simulation, cannot be predicted a

priori. Moreover, other types of discrete change to the system model are also common. For instance. switching a controller from manual to automatic mode can be viewed as a change in which an input variable specification (that fixing the value of the controller output) is replaced by a general equation (that representing the control law). Further complications are introduced by the modelling of impulsive discrete actions. such as the " instantaneous" introduction of a quantity of a certain reactant in a reactor. or resetting the integral error of a proportional/integral controller to zero. Although in reality some of these events occur over a finite length of time, the latter is much shorter than the rest of the simulation. and the overall action is best approximated as an impulse. Clearly. a formal representation of external discrete actions in dynamic process simulation is required. This was one of the main motivations for the development of the gPROMS package (Barton and Pantelides, 1991. Barton, 1992), and in particular. the introduction of the concept of a TASK to describe a set of discrete actions that may operate on a model describing a physical system. At the lowest level. elementary tasks are required to operate directly on the basic characteristics of the model used for the simulation. An analysis of the underlying mathematical problem has revealed that, in fact. only two such fundamental tasks of a discrete nature are strictly necessary. These characterise respectively changes in the describing equations (the REPLACE task), and the effects of impulsive inputs on the system (the REINITIAL task). The equations in a particular simulation comprise symbolic relations between variables and stream connections declared in models. as well as input variable specifications. Therefore the REPLACE task allows the declaration of a very wide range of external actions, such as. for instance, the dynamic alteration of the plant connectivity. or the replacement of one input specification (e.g. an upstream pressure) by another (e.g. a downstream pressure) during the course of the simulation. However, a frequently occurring case is that of the specification of an input variable as a given function of time (or constant value) being replaced by a specification of the same variable as a different function of time (or constant value). The RESET task, a special case of REPLACE, is introduced as a convenient short-hand for defining such simple changes. The REINITIAL task allows the effects of impulses to be described as general sets of equations that relate the condition of the system immediately following the impulse to that before it. For instance. in the case of the instantaneous addition of a quantity of a component to a reactor, such an equation is given by the instantaneous mass and energy balances

Mass-before-addition + Amount-added = Mass-after-addition Internal-energy-before-addition + Internal-energy-added = Internal-energy-after-addition In most simulations. discrete changes of the type introduced above are separated by periods of continuous2 plant operation, the duration of which may be specified explicitly (as a given length of time), or implicitly (as a logical condition expressed in terms of the values of the system variables). Both of these cases are accommodated through the CONTINUE elementary task. which is of the form 2 Further discrete changes may occur during these periods as a result of discontinuities within the physical models itself, as discussed in the previous section.

European Symposium on Computer Aided Process Engineering-2

S271

CONTINUE FOR time-duration or CONTINUE UNTIL logical-condition. More complex tasks can be defined as sets of several elementary tasks. One important issue here is the specification of the relative ordering of the latter. Although tasks can be thought of as procedural entities similar in many respects to statements in conventional programming languages. it must be remembered that. in contrast to the strictly sequential nature of the latter, in the real world multiple actions can occur in parallel. In gPROMS, two distinct constructs, SEQUENCE and PARALLEL, are provided to express the fact that a set of tasks occur in sequence or in parallel respectively (if the SEQ and PAR constructs of the OCCAM parallel programming language (INMOS, 1984)). Each such construct is itself considered to be a task, and therefore they can be nested within each other to any level. TASK Perform_Reaction PARAMETER Reactor SR StartTem perature Termination_Condition

AS MODEL Stirred Reactor AS REAL AS REAL AS LOGICAL_EXPRESSION

SCHEDULE SEQUENCE RESET Reactor.Stearn_Rate := SR ; END; CONTINUE UNTIL Reactor.Temperature > Start Temperature RESET Reactor.Stearn_Rate:= 0.0; END; CONTINUE UNTIL Terminationcondition ; END END Fig. 7.

Example of gPROMS TASK

A simple example of a gPROMS TASK for carrying out a reaction operation in a reactor vessel is shown in Fig. 7. In this case, the steam rate in the heating jacket is maintained at a non-zero value until the temperature reaches the level necessary for starting up the reaction. The steam supply is then stopped, and the reaction is allowed to proceed until a certain termination condition is satisfied. It is enlightening at this point to draw some analogies between the mechanisms used for modelling the physical behaviour of a processing system, and those employed for modelling the external actions imposed on it. One very important issue in modelling is that of model re-usability. As far as models of physical systems are concerned, it has been well understood for several decades that a general model, once developed and verified, can be re-used within either the same or different simulations with minimal additional effort. It was indeed recognition of this benefit that provided the impetus for the development of the first general-purpose sequential modular packages for steady-state simulation.

Similarly. it is equally advantageous to ensure that models of external actions are as general as possible so as to be applicable to multiple situations without modification. One way of achieving this reusability is through extensive task parameterisation, as illustrated in the example of Fig. 7. Once again, these parameters may be thought of as analogous to the arguments of SUBROUTINEs or PROCEDUREs in conventional programming languages. However. much more generality is required in the case of tasks. For example. since most tasks are designed to operate on process units of a particular type, one or more of their parameters may be instances of entire models. Also, it may be preferable to have the termination condition of a task passed to it as a symbolic logical expression parameter, rather than defining a whole range of overall similar tasks that diffcr only with respect to this condition.

S272

European Symposium on Computer Aided Process Engineering-2

TASK StareUp_2_Reactors

SCHEDULE PARALLEL Perform Reaction

(Reactor SR StartTemperature

Tennination_Condition Perform Reaction (Reactor SR StartTemperature

END

Tennination_Condition

IS PIant.Rl.

IS 35.3. IS 70.0. IS PlantR1.Conversion > 0.95) IS Plant.R2. IS 10.0, IS 40.0, IS Plant.R2.Temperature > 90.0)

END Fig. 8.

Example of Higher-Level gPROMS TASK

Describing very complex tasks. such as the start-up or shut-down of an entire plant. entirely in terms of elementary ones is an overwhelming and error-prone undertaking. It is in managing task complexity that a second. particularly useful, analogy can be drawn between modelling physical systems on one hand. and external actions on the other. In the previous section, it was seen that complex plant models can be constructed by the hierarchical combination of lower level models in ftowsheets. Similarly. describing a complex task can be done by assembling lower-level (but not necessarily elementary) ones in sequence or in parallel. For instance. the stan-up of the plant shown in Fig. 3 can be described simply in terms of the start-up in parallel of the reaction and separation tasks. Each one of the latter can then be described in terms of lower level tasks, such as starting individual reactors or columns. At the lowest level, tasks such as opening or closing individual valves can be described in terms of elementary operations such as RESETting the values of input variables representing the stem positions to one or zero respectively. Of course. task parameterisation is crucial in rendering such a hierarchical approach to task description both feasible and useful. An illustration of the ideas of task parameterisation and hierarchical decomposition. and the close relation of the latter to the corresponding concept for models of physical systems. is given in Fig. 8. The task here operates on a plant model comprising two reactor sub-models, RI and R2. It uses the task of Fig. 7 to perform the reaction in both Rl and R2 in parallel. It can be seen that the parameterisation of the lower level task significantly expands the scope of using it in different modes. More complex illustrations of these concepts are provided by Barton et al. (1991). A further enhancement to task re-usability is provided by the use of inheritance hierarchies for building physical system models . As was argued at the end of the previous section, probably the key contribution of inheritance is in grouping a set of diverse models into a single type with a guaranteed minimal set of attributes. This implies that a set of external actions that is entirely defined in terms of these attributes can be applied to all models in this set. For instance. the Perform.Reaction task of Fig. 7 can be applied not only to reactor vessels described by the StirredReactor model, but also to vessels described by any model inherited (directly or indirectly) from StirredReactor. Thus. in the terminology of the object-oriented programming languages. model parameters in tasks (such as the Reactor parameter in the task of Fig. 7) arc polymorphic.

BUILDING COMPLETE SIMULATION DESCRIPTIONS In the previous two sections. powerful mechanisms for the modelling of physical systems and the external actions imposed on them were discussed. The two must be brought together in order to build descriptions of entire simulations. However. some further information is usually required before such descriptions are complete. First. most system models involve more variables than equations . An appropriate number of such variables must therefore be specified (either to fixed numerical values. or to explicit expressions of time) so as to leave a well-posed square system. These then take the role of system inputs. Of course. the subset of variables that can be designated as inputs is not normally unique. On the other hand, not every subset of the appropriate cardinality leads to a problem, the solution of which can be obtained given the current state of numerical solution methods. or even one that possesses a solution! We return to these issues in later sections of this paper.

European Symposium on Computer Aided Process Engineering-2

S273

Another type of information that must be specified to complete the simulation description is that of the initial condition of the system. A fairly general mechanism for specifying these conditions is implemented in the SpeedUp package, as described by Pantelides (1988a). This involves the specification of initial values for a subset of the system differential and/or algebraic variables, and/or of the time derivatives of the former. More general initial condition specifications expressed in terms of equations relating the system variables at the initial time are allowed by gPROMS (Barton, 1992). The well-posedness, consistency and completeness of initial conditions is an issue closely related to the properties of the underlying mathematical system of equations. As such, its treatment is also postponed for later sections of the paper. It should be noted that for system models involving irreversible, or reversible but asymmetric discontinuities, the specification of the initial condition of the system must also include the initial state of each discontinuous equation, as this cannot always be deduced from the values of the system variables. Consider, for instance, the model of the safety valve shown in Fig. 2. If the initial vessel pressure is between the opening and reseat values, then it is impossible to deduce automatically whether the valve is initially open or closed. An explicit specification must therefore be provided. On the other hand, the initial state thus specified may be unstable. For instance, if the state of a bursting disc (illustrated in Fig. 1) is initially set to "intact" but the initial pressure in the vessel exceeds the disc bursting value, then the disc will switch instantaneously to the "burst" state at the very beginning of the simulation. Finally, most attention to date has been focused on the simulation of deterministic systems. In fact, an interesting application of dynamic simulation is in establishing the variability in the behaviour and performance of complex plants subject to stochastic variability in time-independent process parameters and/or initial conditions. Plant behaviour under noisy, time-varying inputs is also of interest, especially for evaluating control system performance. However, despite the obvious attraction of stochastic simulation, very few general-purpose process simulation packages provide facilities for it. A recent example is provided by the work of Diwekar and Rubin (1991) on the public-domain version of the ASPEN steady-state simulator. In the dynamic process simulation field, gPROMS allows the specification of values for both input variables and initial conditions in terms of samples from a number of different probability density functions. Some applications making use of this facility were recently reported by Naf (1992).

SOLUTION METHODS The transient behaviour of most lumped parameter systems can be described mathematically by mixed sets of differential and algebraic equations (DAEs). A consequence of the presence of model discontinuities and/or the imposition of external discrete actions is the partitioning of the time domain of interest [t(O),t(f)] into NCD continuous sub-domains [t(k-l),t(k»), k=1..NCD . Normally the initial time teO) is assumed to be given, while the final time t(f)=.t(NCD) may either be given or be specified implicitly through one or more termination logical conditions. The sub-domain boundaries t(k), k=1..NCD -1 may also be fixed or determined through logical conditions. The dynamic simulation problem is defined as the solution of a sequence of DAE systems of the form j(k) (X(k), i(k), y(k), U(k), t) = 0 t e [t(k-l), t(k») k=1..NCD (1) Here X(k) and y(k) are respectively the system differential and algebraic variables over the kth subdomain, and i(k) the time derivatives of X(k>Ct), all of which must be determined by the solution of the equations. On the other hand, u (k) are the input variables, the time variation of which over the k th sub-domain is known. In general, the transition from the kth sub-domain (k
S274

European Symposium on Computer Aided Process Engineering-2

Before we proceed to consider each of these issues in tum, it is worth noting that one of the advantages of the high-level model representation afforded by many packages is the fact that much more than numerical information is available to the solution process. For instance, the sparsity pattern of the underlying mathematical system can be generated automatically, while symbolic differentiation can be used for generating analytical partial derivative expressions for use by numerical codes (see, e.g., Pantelides, 1988a). The combined application of numerical, structural and symbolic techniques to problem solution is, in our opinion, an important ingredient for the success of any equation-oriented simulation system. Consistent Initialisation The initial condition of the DAE system at the start of the first time sub-domain can be expressed as a set of general nonlinear relationships of the form C(l) (x(l)(t(O», i(1)(t(O\ y(l){t(O», u(l)(t(O~, teo»~ = 0 (3) The issue of the consistency and completeness of a set of initial conditions is particularly complex, and is closely related to the categorisation of DAE system& according to their index (see, e.g., Brenan et al., 1989). A set of consistent initial conditions (x(l)(t~ \ y(l)(t(°l), i(l)(t(O») for any DAE system must certainly satisfy equations (1). However, for DAE systems of index greater than one, and even for some index-l systems, consistency demands that the initial conditions satisfy not only the original equations (I), but also the first, and possibly higher-order, time differentials of some of these equations. The existence of these "hidden" constraints on the initial conditions imply that, unlike the ODE case, the number of additional relations (3) that need to be specified may be smaller than the number of differential variables in the system. Unfortunately as noted by Pantelides et al. (1988), high index problems occur relatively frequently in process engineering applications. Although sometimes they can be avoided through careful modelling (see, e.g., Ponton and Gawthrop, 1991), in other cases they need to be detected and the full set of consistency relations for them must be derived. Pantelides (l988b) presented a graph theoretical algorithm that performs this task by analysing the occurrence structure of an arbitrary DAE system. Modifications to this algorithm, and their implementation have recently been reported by Unger and Marquardt (1991). In theory, because of its structural nature, the algorithm only establishes a set of necessary consistency conditions. However, in practice, experience with many process modelling applications indicates that it usually detects all necessary and sufficient conditions. The problem of initialising the DAE system (1) at the start of sub-domains k=2..NCD is slightly different to that for the first sub-domain. Most discontinuities do not alter the number or the identity of the differential variables in the system''. The usual practice is then to assume that these variables are continuous across the discontinuity, i.e., that x(k)(t(k-I)) = x(k-l)(t(k-I)_)

(4)

the physical reasoning being that the differential variables either represent conserved quantities (e.g. mass, energy), or are directly related to them. However, the above continuity assumption may be problematic in at least two cases. First, if the DAE system /(k)(.) involves "hidden" consistency relations of the type mentioned earlier, some of the continuity relations (4) may be redundant or inconsistent. This is the case when the index of the DAE system increases at the discontinuity. A simple illustration is provided by a vessel that is partially full of an incompressible liquid mixture in the (k-l)th sub-domain becoming completely full at the start of the kth sub-domain. The index of the DAE system increases from one to two at this point. as the previously independent component molar hold-ups arc henceforth constrained by the fixed volume of the vessel and the incompressibility of the liquid. In this case, one of the component hold-up continuity relations is redundant. It is perhaps worth mentioning that we have not, so far, encountered any practical situation in which one or more of the continuity relations becomes inconsistent in this manner. A second situation in which the use of the continuity conditions (4) is not justified is that of the discontinuity at t(k-I) being caused by an impulsive input to the system. For instance, if a certain amount of one component is added instantaneously to a reactor vessel, then neither the hold-up of that component, nor the internal energy of the reactor contents are continuous across this discontinuity. The condition of the reactor following the discontinuity is determined by the instantaneous 3 A discussion of the handling of cases in which this is not true is provided by Barton (1992).

European Symposium on Computer Aided Process Engineering-2

S275

material balance on the component being added, an instantaneous energy balance, and the continuity of the hold-ups of all other components in the system.

In summary, the consistent condition of the system immediately after the discontinuity at time t(k) is determined by the simultaneous solution of (a) the system equations j(k)=:J.J, (b) the hidden consistency relations (if any) implied by these equations, and (c) additional conditions of the general form: C(k)

~(k)(t(k-l»), i(k)(t(k-l»), y
i(k-I)(t(k-lL),

y(k-I)(t(k-lL), u(k)(t(k-I»), t(k-I)]

=0

(5)

In most, but not all, cases, equations (5) takes the simpler form (4), while for the first sub-domain (k=l) they are replaced by (3).

In general, all the above equations form a set of nonlinear algebraic relations that have to be solved numerically. This is, in fact, a far from trivial undertaking, especially for the very first initialisation (k=l), for which good estimates for the values of the variables are not always available. In our experience, the very simple modified Newton algorithms used by most integration algorithms for converging the corrector iterations (see below) are rather inadequate for this purpose. More reliable behaviour has been obtained by using a quasi-Newton type algorithm that has been specially modified for dealing with local Jacobian singularities, and for enhancing convergence from very poor initial estimates (pantelides, 1988c). A priori re-arrangement of the system of nonlinear equations into block triangular form, followed by sequential solution of these blocks is also preferable to simultaneous solution of the entire system, As observed by Kroner et at. (1992), even more powerful algorithms, based on globally convergent (e.g. continuation) methods may be needed for the initialisation of some problems. Further utilisation

of symbolic manipulation may also be of benefit. For instance, symbolic equation transformations could be used for maximising system linearity, and reliable scaling factors for the equations could be derived from an analysis of the relative magnitudes of the various terms occurring in them. The initialisation of problems involving reversible and symmetric discontinuities poses particular problems, as in general, the branch of the discontinuity that is active at the start of the simulation is not known a priori. Some work on the solution of nonlinear algebraic systems that contain such conditional equations has recently been reported by Zaher and Westerberg (1991). DAE System Integration Once the condition of the system at the start of the kth sub-domain has been computed, the next task is the integration of the system of equations from t(k-I) to t(k)4, We consider this problem below. Elimination of Algebraic Equations vs. Simultaneous Approach. Early methods for the solution of index-l DAE systems in the common semi-explicit form gl (x, x, y, t) = 0

g2 (x, y,

t) = 0

(6a)

(6b)

relied on solving the algebraic equations (6b) numerically to obtain the algebraic variables y in terms of the differential ones x S. The algebraic variables can then be substituted into the differential equations (6a), thus effectively converting them to a set of implicit ODEs in x. The latter is then solved using standard ODE methods. However, most sophisticated ODE solution algorithms (especially implicit ones) require multiple evaluations of the residuals of the differential equations, with each such evaluation necessitating the iterative solution of the algebraic equations (6b). Thus, the elimination approach can often be quite inefficient. To remedy this problem, Gear (1971) proposed that the differential and algebraic equations be solved simultaneously at each step. A backward difference formula (BDF) is used to approximate the time 4 Of course, the end of the sub-domain t(k) may not be known a priori; its determination is considered later in this section. 5 Note that for index-I semi-explicit systems, the matrices dg IldX and dg21dY are nonsingular everywhere on the solution trajectory. CACE 17 Suppl-S

S276

European Symposium on Computer Aided Process Engineering-2

derivatives x in terms of the values of x at the current and past steps. This method has the advantage of requiring the solution of only one set of nonlinear algebraic equations per step. A modified Newton iteration is normally used for this purpose (the so called "corrector iterations"). The efficiency of the iteration is further enhanced by starting it from good initial approximations for both differential and algebraic variables obtained via polynomial extrapolation of past solution values (the "predictor step''). Several codes implementing this method are available, the most widely used being the DASSL code developed by Petzold (1983). Recently Perregaard et at. (1992) used several process engineering problems to compare the performance of the elimination and the simultaneous approaches, which they call the "ODE mode" and "DAE mode" respectively. As expected, the simultaneous approach was found to be more efficient when the numerical elimination requires many iterations, or when the evaluation of the residuals of the algebraic equations is very expensive (e.g. due to an internal iterative calculation). Interestingly. the ODE mode appeared to be significantly more efficient for some problems, notably those involving distillation column calculations. However, it should be noted that the problems examined in this study were rather small (not exceeding 135 equations in total), and the numbers of differential and algebraic equations in them were of comparable magnitude. In our experience, most realistic models tend to be much larger than these (by one or two orders of magnitude), and the numberof algebraic equations in them by far exceeds that of the differential equations. For such systems, and provided that currently available efficient sparse linear algebra codes are used, the cost of the simultaneous solution of the two types of equations is usually not much higher than that of a single solution of the algebraic equations. Ease of initialisation is sometimes cited as another advantage of the ODE mode, as the algebraic equations (6b) can be solved independently to yield values of the algebraic variables y that are consistent with given initial values of the differential variables x. However, as discussed earlier in this paper, much more general specifications of initial conditions often need to be accommodated. In any case. the application of block triangularisation to the consistent initialisation equations automatically decomposes the problem into a sequence of the smallest possible sub-problems. Consequently, we agree with Marquardt's (1991) observation that the simultaneous approach is preferable overall. Solution of High Index Problems. As already mentioned, the consistent initialisation of DAE systems of index exceeding unity is particularly difficult. However. even if these problems are overcome, the integration of such systems using standard algorithms of the type described above is problematic. Index 2 systems can still be solved using BDF-type methods starting from consistent initial conditions provided that some attention is paid to the estimation and control of the error of integration. Gupta et al. (1985) showed that, for semi-explicit index-2 systems, the local error in the algebraic variables is of order one less than that in the differential variables. It can therefore be argued that, for the efficient integration of such systems, the algebraic variables should be omitted from the error control mechanisms of the integrator. However, this may be undesirable if one indeed wishes to control the error in some or all of these variables. Gritsis et al. (1988) considered a sub-class of these systems in which the high index is effectively caused by a small subset of the algebraic variables occurring in differential equations only. In this case, it can be shown that it is only these variables that suffer from error order reduction, and therefore need to be excluded from error control. In the case of systems of index exceeding two, the predominant solution approach is one of reducing them to index one or zero by repeated differentiations of the system equations with respect to time. In principle, one could differentiate all equations repeatedly until the system is reduced to a set of ODEs. Whilst this approach is closely related to the definition of index itself, it is both impractical and unnecessary for large systems, in which the high index problem may be caused by a very small subset of the equations. Gear and Petzold (1984) proposed an index reduction algorithm for linear time-varying coefficient systems that only differentiates purely algebraic equations. In general, a set of numerical row manipulations is required to identify the latter. The process is repeated until the system is reduced to a set of ODEs. This technique suffers from the disadvantage that, on a finite precision computer, rounding error accumulating during the linear algebra operations may obscure the distinction between differential and algebraic equations, a problem that tends to be especially severe for large systems. Furthermore, complete reduction to ODEs is unnecessary as most index one systems can also be solved with effectively the same degree of difficulty. The process of index reduction through equation differentiation is also closely related to the task of deriving a complete set of consistency relations for the initial conditions of a DAE system. In particular, it can be shown that the index of the system formed by the highest order differentials of each equation in these relations (or the original equation, if this has not been differentiated) does not

European Symposium on Computer Aided Process Engineering-2

S277

exceed unity. Therefore algorithms for consistent initialisation can also be employed for index reduction. One advantage of such an approach is that these algorithms determine the minimal number of consistency relations, and therefore are selective as to which equations they differentiate and how many differentiations they apply (Pantelides, 1988b, Unger and Marquardt, 1991). Mathematically, it is easy to show that a solution that (a) is consistent with respect to all the equations and their time differentials at the initial time, and (b) satisfies the highest-order time differential of each equation over the time domain of interest, also satisfies the original high index system equations over the same time domain. Therefore, it would appear that all one has to do to solve any high index system is simply to derive the consistency relations, solve them to determine a set of consistent initial conditions, and then integrate the index-l (or 0) system formed from the highest order time differentials of the equations, discarding all others. Unfortunately, because of the effects of integration and rounding errors, this may not be true for numerical implementations of such an algorithm on a finite-precision computer. The problem is that after a number of steps, the solution obtained does not necessarily satisfy the original equations within the required accuracy , despite satisfying their highest-order differentials. The above difficulty can be overcome by retaining both the original and the differentiated forms of each equation in the system. This, however, leads to a redundant set involving more equations than variables. In an attempt to remedy this problem, Bachmann et al. (1990) use numerical elimination to produce algebraic equations that only involve the differential variables x. If this is not possible, the index of the system is at most 1. Otherwise these equations are differentiated with respect to time, and the resulting differentials are converted into algebraic relations by using the original differential equations to eliminate all occurrences of the time derivatives i from them. The redundancy created by the introduction of these new relations is removed by dropping an equal number of the differential equations used in the elimination. The procedure is repeated until the index problem is removed. Whilst the Bachmann et al approach results in a system that has the desired numerical characteristics, the use of elimination at two of its steps limits its applicability, especially for large nonlinear systems. Even for linear systems , rounding errors may also create problems (see comments on the Gear and Petzold (1984) method above). A simpler approach has very recently been proposed by Mattsson and Soderlind (1992). Their "dummy derivative" method effectively removes the redundancy created by the introduction of the differentiated equations by treating some of the time derivatives Xi as new, independent algebraic variables ignoring their relationship to the corresponding differential variables Xi ' The method does not require any symbolic or numerical eliminations, relying entirely on the equation differentiations determ ined by the application of Pantelid es's (1988b) algorithm . Other approaches to the solution of high-index problems have been proposed by Gear (1988) and Chung and Westerberg (1990). Beyond Index Problems. Notwithstanding the importance of problems posed by high-index DAE systems, it should be recognised that these are not the only major outstanding ones. In fact, the increasing use of simulation for the study of process engineering systems under conditions that are far from steady-state has revealed new challenges for numerical integration codes. Jarvis and Pantclides (1991) identified two such classes of problems. The first involves so-called partially determinable systems. the equations of which do not determine unique values for all variables over certain parts of the time domain. Although, in many cases , this indeterminacy can be traced back to the lack of modelling detail, it can be argued that integration of such systems should still be possible provided that unique values are obtainable for all variables of primary interest. However, most existing DAE solvers simply fail the moment the system becomes partially indeterminate . A second class of difficult problems are those involving fast transients. These occur very frequently in practice (see, e.g., Heyen et al., 1992), creating problems for existing solvers which often tend to take unnecessarily small steps in the region of the transients. A possible remedy. based on a novel measure of the integration error has been proposed by Jarvis and Pantelides (1991). This. together with algorithms for the automatic detection and integration of partially determinable systems, has been implemented in the DASOLV code (Jarvis and Pantelides, 1992). Finally, despite their overall efficiency, sophisticated codes based on BDF approximations may behave quite inefficiently in the presence of frequent discontinuities such as those arising from the interaction of dynamic simulations with real-time control systems (e.g. for operator training applications). The presence of discontinuities prevents the codes from using previous step information to achieve high integration orders; in the worst case, the entire integration is effectively carried out with an implicit Euler method and a very small step size.

5278

European Symposium on Computer Aided Process Engineering-2

A possible remedy to the above problem may be the use of high-order single-step (e.g. Runge-Kutta) methods, as these would allow a relatively large step to be taken immediately following the discontinuity. Furthermore, although the stiffness of most process engineering problems usually dictates the use of implicit integration methods, explicit methods could be advantageous for problems in which the time interval between discontinuities does not exceed the stability limit Adaptive Runge-Kutta methods that automatically switch from implicit to explicit form and vice versa (Cameron and Gani, 1988), may also be particularly suitable in these situations. A different approach, suggested by Gear (1980), is to use the Runge-Kutta method only for the initial steps following the discontinuities before switching to the usual BDF approach. Discontinuity Location As described in the introduction to this section, the end of each interval of continuous operation is determined by the logical condition (2) becoming true. In the general case, such logical conditions may involve the values of the unknown variables themselves, and therefore define implicit discontinuities, the time of occurrence of which cannot be determined a priori. A significant body of literature exists on the problem of detecting and locating implicit discontinuities. One class of approaches is based on formulating an appropriate "switching function", which changes sign across the discontinuity. A sign change detected during the integration then triggers a search for the precise location of the discontinuity (Smith and Morton, 1988, Preston and Berzins, 1991). Implementations of this technique differ as to whether checks for sign changes are carried out at any stage during the integration (e.g. after the predictor step, or during the corrector iterations), or only at points determined as lying on the solution trajectory. Both suffer from certain disadvantages. A sign change detected at a point not on the solution trajectory may be entirely spurious. On the other hand, postponing the search for the discontinuity location until a "true" solution point past it has been obtained implies incurring all the cost associated with integrating across discontinuities (e.g. failed corrector iterations, and/or error tests). In either case , the task of actually locating the discontinuity following its detection can itself be a difficult, and often costly, undertaking. Fortunately, an alternative approach that seeks to address the above problems, and one that is particularly suited to the mathematical problem described by equations (1) and (2) does exist. It is simply implemented by checking the status of the logical conditions (2) only at the end of each successful integration step. If one or more of these are found to be true, then interpolation is used over the last step to determine the earliest time at which this switching occurred . The main advantage of this approach is that every integration step is taken on a continuous DAE system, and this normally minimises the cost involved in terms of failed steps. Furthermore. because of the continuity of the step across the discontinuity, efficient interpolation procedures can be used to locate the discontinuity within the required accuracy at minimal cost. Finally, no ad hoc switching functions are required - all necessary logical checks are carried out using the logical conditions naturally associated with the equation definition. Joglegar and Reklaitis (1984) employed this technique in their batch process simulator, with each DAE system (1) corresponding to a different batch task, the end of which is determined by logical conditions (2). Pantel ides (1988a) used it for handling conditional equations representing reversible and symmetric discontinuities within the SpeedUp package . In this case, the techn ique can be interpreted as determining the branch of the equation that is active at the stan of each integration step. and then "locking" the equation in this branch for the durat ion of the step - hence the term discon-

tinuity locking. One potential disadvantage of the above approach is that it assumes that a solution point can still be found past the discontinuity at t(l:) even though the wrong set of describing equations (namely J (t)=:{] instead of J (I:+1 )={)) is used for the last step. As this is not generally true, the method may fail to establish a point past the discontinuity, which inevitably leads to repeated corrector failures and severe reductions in step-length. However, it may be worth mentioning that, despite the very extensive use of the technique in process dynamic simulation over the past five years, to our knowledge there has been little evidence of serious problems, provided some modest care is exercised in coding conditional equations.

European Symposium on Computer Aided Process Engineering-2

S279

PROBLEM DIAGNOSIS It is important for any equation-oriented package to be able to detect and diagnose problems that either are inherently badly-posed, or cannot readily be solved using currently available techniques. Problems in the first category include DAE systems that either do not possess any solution, or do not possess a solution satisfying the given initial conditions. On the other hand, most high index problems fall in the second category. In this section. three potential sources of problems are examined, namely lack of DAE solvability, high index, and inconsistent initial condition specification. In all of these, the difficulty of diagnosis is compounded by the fact that problems may be caused not only by the general form of the modelling equations, but also by the particular selection of input variables. i.e., the partitioning of the variables in (1) into known inputs u"(lc), and unknowns x(kl and y(k l. Therefore. the task of aiding the user in arriving at a valid selection of input variables is closely related to the problem diagnosis. Solvability of DAE Systems In contrast to the well-known classical results on the solvability of systems of ODEs. very little theory exists on the existence and uniqueness of solutions of general nonlinear DAE systems. For instance, the recent monograph by Brenan et at. (1989) limits itself to defining solvability for nonlinear. and even linear, time-varying coefficient. DAE systems . Theorems actually establishing conditions under which solutions exist are given only for the linear. constant coefficient system A

x + B x = v (t )

(7)

In particular. it can be shown that the above system is solvable if and only if the matrix pencil (M + B) is regular. i.e. the function

(A.) ;; det(M + B) is not identically zero for all values of the scalar A.. D

(8)

In the absence of theoretical results for nonlinear systems, and by analogy to (8). we relate the existence of a solution to the DAE system (1) to the nonsingularity of the following matrix on the solution trajectory'':

[

A.~ + af af] ax

ax ay

(9)

A further practical justification for this assumption is provided by the use of the above matrix in numerical integration algorithms. As already described , the application of implicit integration algorithms to (1) requires the solution of a set of nonlinear equations at every step . It can easily be shown that the Jacobian of these equations is of the form (9), with A. being inversely proportional to the step-length of integration. Thus. if (9) were identic ally singular. numerical methods would probably fail , and (1) would be categorised as "unsolvable" for all practical purposes. Accepting (9) as a solvability criterion, we are still faced with the problem of checking for singul arity of a potentially large, sparse matrix . and in particular one. the elements of which depend on the precise point at which the various partial derivatives in (9) are evaluated. Although this is a difficult undertaking in general, often the singularity of (9) is independent of the symboli c form or the numerical values of its elements. but only depends on its structure. The latter is easily obtained from the sparsity pattern of the Jacobian of ( I) [f x f:i f y] by merging the sparsity patterns of f i. and f x ' i .e. replacing each occurrence of the time derivative Xi by one of the corresponding differential var iable Xj, and then removing any duplicate elements that may arise. Structural singularity in large matrices can be detected very efficiently using graph -theoretical algorithms for output set assignment. in particular one proposed by Duff (1980). If matrix (9) is indeed structurally singular. Duff's algorithm detects at least one equation subset. the cardinality of which exceeds by exactly one the number of distinct variables x and/or y occurring in it One of the equations in this subset must therefore be either inconsistent or redundant. If the equations actually involve one or more input variables u , it may be possible to attribute the problem to the 6 In the interests of simplicity, we drop the superscript (k) from the notation in this section.

S280

European Symposium on Computer Aided Process Engineering-2

corresponding input variable specification(s). On the other hand. if this is not the case, then an error in the construction of the physical model is detected. In either case. the user may be advised accordingly. Structural analysis may also offer valuable help in the converse problem of guiding the user in selecting input variable specifications so as to arrive at a well-posed DAE system. In this case, the matrix (9) is initially rectangular, with the number of columns exceeding that of the rows by the number of outstanding input specifications. Sargent (1978) showed how such rectangular matrices can easily be rearranged to a block triangular form, such that (a) all but the last block are square and structurally nonsingular, and (b) the last block is rectangular and structurally nonsingular. It follows that further input variable specifications can only be selected from this last block. Of course. the latter may decompose further once an additional specification is made. Thus the block triangularisation procedure must be applied after each specification until all blocks are square. Fortunately. extremely efficient algorithms are already available for all necessary structural manipulations (Duff, 1980, Tarjan. 1972). Of course, it should be borne in mind that structural nonsingularity of a matrix does not guarantee that the matrix is not identically or locally singular. Consequently, the structural techniques outlined above and in the rest of this section are useful in eliminating some, but not necessarily all, badlyposed problems. Avoiding High Index DAE Problems Despite the recent advances in understanding high index DAE systems, robust, reliable and efficient codes for the solution of such systems are still largely unavailable. Furthermore, in our experience, many - although by no means all - high index problems are caused by assumptions that may not be essential, or even simply by modelling errors. It may therefore be important for the modeller to be provided with some understanding regarding the causes of high index, even in those cases that it is possible to solve the associated mathematical problem. It can be shown that all systems of the general form (1), the index of which exceeds unity, are such that the matrix

[1f ~ ]

(10)

is singular. The above matrix is also singular for some index-I systems, but even these pose special difficulties with respect to the specification of consistent initial conditions. Hence, the singularity of the above matrix can be used as an overall indication of index-related difficulties. Despite the difference in the forms of (9) and (10), the techniques that can be employed for their analysis are quite similar. In particular, structural singularity can be used to identify input specifications and/or equation subsets that may be responsible for the high index problem. Conversely, the user may be guided towards selecting a set of input specifications that lead to a DAE system that is not only solvable in principle, but of index-I as well. If a complete characterisation of the index of a given DAE system (rather than simply an indication of whether it is of index 1) is required, then one may use Pantelides' (1988b) structural algorithm, the first step of which is indeed the analysis of the structural singularity of (10). Interestingly, it has been shown that the algorithm terminates if and only if (9) is structurally nonsingular. Thus, the structural analysis for system solvability must always be carried out before attempting to determine the system index. Well-posed Initial Conditions In general. the specification by the user of a consistent set of initial condition specifications of type (3) may be a non-trivial undertaking, even for index-l systems for which matrix ~1O) is nonsingular. For this class of problems, a set of consistent values for x(t(O), x(t(O) and yet ») is simply determined by the solution of conditions (3) (provided by the user) together with the system equations (1) at the initial point, i.e., there are no "hidden" initial conditions of the type mentioned earlier in this paper. Although. strictly speaking. nonsingularity of the Jacobian of a set of nonlinear algebraic relations is not a necessary condition for the existence of real solutions, it can be considered to be so for most practical problems. We therefore require that the following matrix be nonsingular:

European Symposium on Computer Aided Process Engineering-2 fJj(l}

fJj(l)

fJf(l)

ai(l)

fJX(I}

fJy(l)

fJC(l}

fJC(l}

fJC(l)

Cli(l}

fJx(l)

fJy(l)

S281

(1 I)

Once again, structural singularity concepts may be used for detecting inconsistencies in the user specified relations (3). Furthermore, when (3) take the special form of simply specifying initial values for a subset of the variables, some guidance can be provided to the user as to the correct selection of these variables.

FUTURE PERSPECfIVES Despite the recent progress in equation-oriented process modelling and simulation, there are still many outstanding issues to be resolved. Some of these have already been pointed out earlier in this paper, their resolution, especially where solution techniques are concerned, will greatly increase the reliability of the available tools. Other areas, such as the development of better user interfaces (e.g. Stephanopoulos et ai., 1987, Bar and Zeitz, 1989, Westerberg at al., 1991), and improved data organisation and representation (e.g. Marquardt, 1992) are also particularly important in ensuring the most widespread use of this technology. Beyond the kind of developments mentioned above, however, there are those that are of a more strategic nature, in the sense that fulfilling them will allow the definition and solution of entirely new classes of problems within the unified framework of general-purpose process modelling environments. The rest of this section is devoted to two such potential developments. Declarative Modelling and Simulation of Distributed Parameter Systems A substantial proportion of all process unit (such as tubular reactors, packed separation columns and pipelines) are distributed in nature, i.e. the conditions within them vary with respect to both time and position determined in terms of one or more spatial coordinates. The conditions in other unit operations (e.g. crystallisers, polymerisation reactors and driers for particulate solids) are described by distributions of properties (such as crystal size, product molecular weight, or particle moisture content). From a mathematical viewpoint, the detailed modelling of such systems generally leads to mixed systems of partial and ordinary differential and algebraic equations (PDAEs), or even integral PDAEs (IPDAEs). As noted by Marquardt (1991), no current general-purpose process modelling environment supports either the direct modelling or the numerical treatment of such systems. In the area of modelling, one solution is to devise model description languages that allow the user to express relationships between variables directly in symbolic IPDAE form. An analysis of the issues involved, and a design for such a language as an extension of gPROMS are the objects of on-going research at Imperial College (Oh, 1992). On the numerical front, the well-known method of lines that discrctises all independent variables except time, thus reducing the IPDAE system to one of DAEs, seems to be the most appropriate as it

naturally allows the simulation of systems involving both distributed and lumped unit operations using well-established time integration techniques (see, e.g. Heydweiller et al., 1977). Of course, this kind of transformation could be done automatically from a high-level representation of the modelling equations, perhaps with some guidance from the user as to the most suitable type of discretisauon (e.g. finite differences, finite clements, orthogonal collocation on finite clements etc.) for each spatial dimension. The main issue, however, is that of controlling the spatial discretisation error to the required accuracy while minimising the size of the resulting DAE system. The use of non-uniform discretisation grids is particularly important for the accurate solution of problems involving steep spatial variations in one or more dimensions. Furthermore, for problems with moving fronts, it is often essential for the position of the grid to change with time, either continuously (as in the moving finite element method first proposed by Miller and Miller, 1981, and its derivatives), or at discrete time instances (as in the class of adaptive grid methods, Flaherty et al., 1989).

S282

European Symposium on Computer Aided Process Engineering-2

Of course, much work has already been done in related areas (such as computational fluid dynamics).

in analysing individual classes of problems of PDE systems and developing the most efficient methods for each of them (see. for instance, Fletcher, 1988). However, one could argue that the key requirement in the context of general-purpose process modelling environments, is primarily one for methods that are robust and reliable for the solution of general PDAE systems subject to general boundary conditions - even if these attributes are achieved, to a certain extent, at the expense of solution efficiency. A class of general higher-order moving finite element methods fulfilling some of these requirements were proposed by Pipilis (1990), but much .remains to be done to achieve uniformly acceptable solution efficiencies. Provision of General Optimisation Capabilities The possibility of performing both simulation and optimisation calculations within a unified environment, potentially utilising a common system model, is one of the main attractions of equation-oriented modelling packages. Thus. several of the packages mentioned in this paper (such as SpeedUp and ASCEND) already offer facilities for steady-state optimisation alongside their simulation facilities. However, much less attention has been devoted to dynamic optimisation problems, despite their importance in diverse practical applications, such as the optimal control of batch distillation and reaction operations, the determination of optimal feedstock changeover policies etc. We believe that this is a particularly important area for the evolution of equation-oriented packages. As in the case of distributed parameter systems, two issues have to be addressed. First, a natural, high-level representation of general optimisation problems in this class must be devised . Secondly. general numerical algorithms and codes for the solution of large-scale problems are needed. The latter problem has received much attention in recent years with the development of algorithms based on both the control vector parameterisation approach (Sargent and Sullivan, 1977. Morison, 1984. Vassiliadis, 1992), and the complete parameterisation approach (Renfro et al.• 1987. Logsdon and Biegler, 1989, 1991). Despite their still rather ad hoc treatment of inequality path constraints, these algorithms have already \cd to the development of codes that can solve many optimisation problems of practical significance. On the contrary, little work has been done on devising languages for the formal high-level description of these problems, comparable in power and generality to those now available for defining simulation problems. For applications that seek to derive an optimal set of external control actions to be imposed on a system. a generalisation of the concept of a TASK, as used in gPROMS, might provide an appropriate mechanism for attaining this aim.

CONCLUDING REMARKS As hopefully this paper indicates, substantial progress has already been achieved in all aspects of equation-oriented dynamic simulation. Effective tools are now available for modelling both complex physical systems and the complex actions imposed on them. thus covering a wide spectrum of processes ranging from purely continuous to batch operations . Whilst many applications can adequately be served by relatively simple models, we believe that the main thrust in the future will be in the development of detailed, high-fidelity models. These will both benefit from. and support, the on-going progress in other areas such as chemical reaction engineering and separation processes. Clearly, the introduction of facilities for the direct modelling of distributed unit operations would be a significant step in this direction. Given the difficulty of building correct mathematical descriptions of very complex models, the current research into tools for model construction from purely physical descriptions is of great importance. One particularly significant aspect of the construction of high-fidelity models is taking account of the detailed geometry of plant equipment items and their interconnections. The natural and effective representation of the latter perhaps deserves to receive more attention to complement that already directed towards the representation of basic physicochemical mechanisms and their interactions in unit operations. In the area of solution methods, notable achievements have been recorded over the past decade. Through the combined application of symbolic. structural and numerical techniques, and the continuing advances in computer hardware, it is now possible to solve practical problems of significant size

European Symposium on Computer Aided Process Engineering-2

S283

and complexity. Of course, with the ever widening range of applications of dynamic simulation, and the increasing flexibility in problem specification allowed by the available software, new classes of mathematical and numerical problems will continue to be identified.

REFERENCES Andersson, M. (1990). Omola - An Object-Oriented Language for Model Representation. PhD Thesis, Lund Instritute of Technology. Aspen Technology (1991), SpeedUp User Manual. Cambridge, United Kingdom. Bachmann, R, L., Brilll, Th. Mrziglod, and U. Pallaske (1990). On Methods for Reducing the Index of Differential Algebraic Equations. Comput. chem. Engng., 14, 1271-1273. Bar, M, and M. Zeitz (1989). A Knowledge-Based Flowsheet-Oriented User Interface for a Dynamic Process Simulator. Comput. chem. Engng., 14 , 1275-1280. Barton, P.I. (1992). The Modelling and Simulation of Combined Discrete/Continuous Processes. PhD Thesis, University of London. Barton, P.I. and C.C. Pantelides (1991). The Modelling and Simulation of Combined Discrete/Continuous Processes. Proc. PSE'91 Conf', Montebello. Canada. Barton, P.I., E.M.B. Smith. and C.c. Pantelides (1991). Combined Discrete/Continuous Process Modelling Using gPROMS. Paper 152c, AIChE Annual Meeting, Los Angeles. Brenan, K.E., S.L. Campbell, and L.R. Petzold (1989). Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations. North-Holland. Cameron, LT., and R. Gani (1988). Adaptive Runge-Kutta Algorithms for Dynamic Simulation. Comput. chem. Engng., 12 , 705-717. Cellier, F.E., and A.P. BonguIielmi (1980). The COSY Simulation Language. In: Simulation of Systems'79 (L. Dekker, G. Savastano, and G.c. Vansteenkiste, eds.), pp. 271-281. North-Holland. Czulek, A.J. (1988). An Experimental Simulator for Batch Chemical Processes. Comput. chern. Engng., 12 , 253-259. Chung Y., and A.W. Westerberg (1990). A Proposed Numerical Algorithm for Solving Nonlinear Index Problems. Ind. Eng. Chem. Res., 29 , 1234-1239. Diwekar, U.M., and E.S. Rubin (1991). Stochastic Modeling of Chemical Processes. Comput. chem. Engng., 15 , 105-114. Duff, LS. (1980). On Algorithms for Obtaining a Maximum Transversal. ACM Trans. Math. Softw., L. 315-330. Elmquist, H. (1978). A Structured Model Language for Large Continuous Systems. PhD Thesis, Lund Institute of Technology. Flaherty, J.E., PJ. Paslow, M.S. Shephard, and J.D. Vasilakis (eds.) (1989). Adaptive Methods for Partial Differential Equations, SIAM Publications. Fletcher, CAJ. (1988). Computational Techniques for Fluid Dynamics - Fundamental and General Techniques. Springer-Verlag. Gear, C.W. (1971). Simultaneous Numerical Solution of Differential-Algebraic Equations. iEEE Trans. Circ. Theory, CT-18 , 89-95. Gear, C.W. (1980). Runge-Kutta Starters for Multistep Methods. ACM TOMS, !5. , 263-279. Gear, C.W. (1988). Differential-Algebraic Equation Index Transformation. SIAM J. Sci. Stat. Comput., 2 , 39-47. Gear, C.W., and L.R. Petzold (1984). ODE Methods for the Solution of Differential/Algebraic Systems. SIAM J. Numer. Ana!., 2.l , 716-728. Gritsis, D., C.C. Pantelides, and R.W.H. Sargent (1988). The Dynamic Simulation of Transient Systems Described by Index Two Differential-Algebraic Equations. Proc. Third lnt. Symp. Proc, Syst. Engng., Sydney, Australia. Gupta, G.K., C.W. Gear, and B. Leimkuhler (1985). Implementing Linear Multistep Formulas for Solving DAEs. Technical Report UIUCDCS-R-85-1205, Dept. of Computer Science, Univ. of Illinois at Urbana-Champaign. Henning, G., G. Stephanopoulos, and H. Leone (1990). MODEL.LA. A Modeling Language for Process Engineering. Part 1: The Formal Framework. Comput. chern. Engng., 14 , 813-869. Heydweiller, J.C, R.F. Sincovec, and L.T. Fan (1977). Dynamic Simulation of Processes Described by Distributed and Lumped Parameter Models. Comput. chern. Engng., 1 , 125-131. Heyen, G, B. Kalitventzeff. P. Hutchinson, H. Yeung, and M. Gill (1992). Simulation of Fast Transients in Fluid Transport Equipments and Utility Networks. In: European Symposium on Computer Aided Process Engineering -I, (R. Gani, ed.), , pp. Sl09-Sl16. Pergamon Press. Holl, P., W. Marquardt, and E.D. Gilles (1988). DIVA - a Powerful Tool for Dynamic Process Simulation. Comput. chem. Engng., 12,421-426. INMOS (1984). OCCAM Programming Manual. Prentice Hall. Jarvis, R.B., and c.c. Pantelides (1991). Numerical Methods for Differential-Algebraic Equations New Problems and New Solutions. Paper 158c, AIChE Annual Meeting, Los Angeles.

S284

European Symposium on Computer Aided Process Engineering-2

Jarvis, R.B., and C.C. Pantelides (1992). DASOLV: A Differential-Algebraic Equation Solver. Technical Report, Centre for Process Systems Engineering, Imperial College, London. Joglegar G.S, and G.V. Reklaitis (1984). A Simulator for Batch and Semi-continuous Processes. Compo chem. Engng., 315-327. Kroner, A., P. Holl, W. Marquardt, and E.D. Gilles (1990). DIVA - an Open Architecture for Dynamic Process Simulation. Comput. chem. Engng., 14, 1289-1295. Kroner. A., W. Marquardt, and E.D. Gilles (1992). Computing Consistent Initial Conditions for Differential-Algebraic Equations. In: European Symposium on Computer Aided Process Engineering -I, (R. Gani, ed.), pp. S131-S138. Pergamon Press. Logsdon, J.S., and L.T. Biegler (1989). On the Accurate Solution of Differential-Algebraic Optimization Problems. Ind. Eng. Chem. Res., 28 , 1628-1639. Logsdon, 1.S., and L.T. Biegler (1991). Decomposition Strategies for Large-Scale DifferentialAlgebraic Optimization Problems. Chem. Eng. Sci., 47 , 851-864. Macchietto, S., and S. Kassianides (1991). Use of Dynamic Simulation in Operator Training. Paper presented at European Simulation Multicon/erence'91, Copenhagen, Denmark. Mani, S., S.K. Shoor, and H.S. Pedersen (1990). Experience with a Simulator for Training Ammonia Plant Operators. Plant/Operations Progress, 2 , 6-11. Marquardt, W. (1991). Dynamic Process Simulation - Recent Progress and Future Challenges. In: Chemical Process Control (Y. Arkun and W.H. Ray, eds.), pp. 131-180. CACHE-AIChE Publications. Marquardt, W. (1992). An Object-Oriented Representation of Structured Process Models. In: European Symposium on Computer Aided Process Engineering - I, (R. Gani, ed.), pp. S329-S336. Pergamon Press. . Mattsson, S.E., and G. SOderlind (1992). Index Reduction in Differential-Algebraic Equations Using Dummy Derivatives. SIAM 1. Sci. Stat. Comput. (to appear). Meyer. B. (1988). Object-Oriented Software Construction. Prentice Hall. Miller, K., and R.N. Miller (1981). Moving Finite Elements. SIAM J. Numer. Anal.• li. 1019-1032. Morison, K.R. (1984). Optimal Control of Processes Described by Systems of Differential-Algebraic Systems. PhD Thesis, University of London. Naess, L., A. Mjaavatten, and J-O. Li (1992). Using Dynamic Process Simulation from Conception to Normal Operation of Process Plants. In: European Symposium on Computer Aided Process Engineering -I. (R. Gani, ed.), pp. S119-S130. Pergamon Press. Niif, U. (1992). Stochastic Process Simulation using gPROMS. Technical Report. Centre for Process Systems Engineering, Imperial College. London. Nilsson, B. (1989). Structured Modeling of Chemical Processes with Control Systems. Paper 27g. AIChE Annual Meeting. San Francisco. Oh, M. (1992). A General Process Modelling Environmentfor the Dynamic Simulation of Distributed Parameter Systems. Technical Report. Centre for Process Systems Engineering. Imperial College. London. Pantelides, C.C. (1988a). SpeedUp - Recent Advances. in Process Simulation. Comput. chem. Engng., 12 , 745-755. Pantelides, c.c. (l988b). The Consistent Initialization of Differential-Algebraic Systems. SIAM J. Sci. Stat. Comput., 2 . 213-231. Pantelides, C.C. (1988c). Symbolic and Numerical Techniques for the Solution of Large Systems of Nonlinear Algebraic Equations. PhD Thesis. University of London. Pantelides, c.c., D. Gritsis. K.R Morison, and R.W.H. Sargent (1988). The Mathematical Modeling of Transient Systems Using Differential-Algebraic Equations. Comput. chem. Engng.• .ll ,449454. Perkins, J.D. (1983). Equation-Oriented Flowsheeting. In: Foundations of Computer Aided Chemical Process Design (A.W. Westerberg and H.H. Chien. eds.), pp. 309-367. CACHE Publications. Perkins. J.D .• and G.W. Barton (I987). Modelling and Simulation in Plant Operation. In: Computer Aided Process Operations (G. V. Reklaitis and H.D. Spriggs. eds.), pp. 287-316. CACHEElsevier. Perkins, J.D., and RW.H. Sargent (1982). SPEEDUP - a Computer Program for Steady-State and Dynamic Simulation and Design of Chemical Processes. AlChE Symp. Series. 214 • I-II. Perregaard, J., B.S. Pedersen. and R. Gani (1992). Steady-state and Dynamic Simulation of Complex Chemical Processes. Chem. Engng. Res. Des.• 70 • 99-109. Petzold, L.R. (1983). A Description of DASSL: a Differential/Algebraic System Solver. In: Scientific Computing (R.S. Steplernan et al., cds.), pp. 65-68. North-Holland. Piela, P. (1989). ASCEND: an Object-Oriented Computer Environment for Modeling and Analysis. PhD Thesis. Carnegie Mellon University. Piela, P., T.G. Epperly, K.M. Westerberg. and A.W. Westerberg (1991). ASCEND: an ObjectOriented Computer Environment for Modeling and Analysis: The Modeling Language. Comput.

.a,

chem. Engng., U. 53-72.

Pipilis, K.G. (1990). Higher-Order Moving Finite Element Methods for Systems Described by Partial Differential-Algebraic Equations. PhD Thesis. University of London.

European Symposium on Computer Aided Process Engineering-2

S285

Ponton, J.W., and P.J. Gawthrop (1991). Systematic Construction of Dynamic Models for Phase Equilibrium. Comput. chem. Engng., li , 803-808. Preston, AJ.• and M. Berzins (1991). Algorithms for the Location of Discontinuities in Dynamic Simulation Problems. Comput. chem. Engng.• li .701-713. Renfro. J.G., AM. Morshedi, and O.A. Absjomsen (1987). Simultaneous Optimization and Solution of Systems Described by Differential/Algebraic Equations. Comput. chem. Engng.. 11 • 503· 517. Sargent, RW.H. (1967). Integrated Design and Optimization of Processes. Chern. Eng. Progr.• Q.3. (9), 71-78. Sargent, RW.H. (1978). The Decomposition of Systems of Procedures and Algebraic Equations. In: Lecture Notes in Mathematics, Q3.Q • pp. 158-178. Springer-Verlag. Sargent. RW.H., and G.R Sullivan (1977). The Development of an Efficient Optimal Control Package. Proc. 8th IFfP ConJ. Optim. Tech. Pt. 2. Shacham, M., S. Macchietto, L.F. Stutzman. and P. Babcock (1982). Equation Oriented Approach to Process Flowsheeting. Comput. chem. Engng., Q• 79-95. Shinohara, T. (1987). A Block Structured Approach to Process Simulation. In: Computer Aided Process Operations (G.V. Reklaitis and H.D. Spriggs. eds.), pp. 317-331. CACHE-Elsevier. Smith, G.l, and W. Morton (1988). Dynamic Simulation Using an Equation-Oriented Flowsheeting Package. Comput. chem, Engng.• 12 .469-473. Stephanopoulos, G.. J. Johnston, T. Kriticos, R. Lakshmanan, M. Mavrovouniotis, and C. Siletti (1987). DESIGN-KIT - an Object-Oriented Environment for Process Engineering. Comput. chem. Engng., 11 , 655-674. Tarjan, RE. (1972). Depth-First Search and Linear Graph Algorithms. SIAM 1. Comput., 1 . 146160. Unger, J.• and W. Marquardt (1991). Structural Analysis of Differential-Algebraic Equation Systems. In: Computer-Oriented Process Engineering (L. Puigjaner and A Espuna, eds.), pp. 241-246. Elsevier. Vassiliadis, V.S. (1992). DAEOPT- A Differential-Algebraic Optimisation Solver. Technical Report. Centre for Process Systems Engineering, Imperial College, London. Vazquez-Roman, R. (1992). Computer Aids for Process Model Building. PhD Thesis. University of London. Westerberg, A.W., P.C. Piela, RD. McKelvey, T.G. Epperly (1991). The ASCEND Modeling Environment and its Implications. In: Proc. PSE'91 Conf.; Montebello, Canada. Zaher, ].J., and AW. Westerberg (1991). Conditional Modeling in an Equation-Based Modeling Environment. Paper l38d, AIChE Annual Meeting, Los Angeles.