Daedalus Modeling Framework: Building First-Principle Dynamic Models

Article pubs.acs.org/IECR

Daedalus Modeling Framework: Building First-Principle Dynamic Models Joaõ R. Leal,*,† Andrey Romanenko,*,‡ and Lino O. Santos*,† †

CIEPQPF, Department of Chemical Engineering, Faculty of Sciences and Technology, University of Coimbra, 3030-790 Coimbra, Portugal ‡ Ciengis, SA, 3030-199 Coimbra, Portugal S Supporting Information *

ABSTRACT: First-principles-based dynamic models of industrial-scale chemical processes involve the definition of many equations and variables, which can demand a considerable effort to setup and maintain. Existing tools require the modeler to either write all of the equations or to select them from a repository. The Daedalus Modeling Framework is a novel software library that expedites the creation of nonlinear dynamic lumped parameter models by employing an automatic equation and property selection algorithm. It uses a phenomena-oriented modeling strategy to define a higher-level model description, which can generate different systems of equations, depending on the available model inputs and requested outputs. This paper describes the design decisions and implementation strategy for this framework, such as its basic building blocks, model structure selection algorithm, differentiation index reduction, algorithmic differentiation, and documentation generation. Two examples demonstrate the modeling capabilities of Daedalus and its integration with a nonlinear model predictive control framework.

1. INTRODUCTION AND MOTIVATION Dynamic modeling frameworks are important tools in aiding process engineers in translating their knowledge of physical, chemical, and biological principles or even empirical understanding of a process into a sometimes highly complex, mathematical model. In a significant number of areas in process system engineering, model development plays a key role. In the area of NMPC, for instance, the development and identification of a suitable model can be one of the most time-consuming tasks during commissioning, representing ∼80% of the total time and expense.1 This fact has also been pointed out as being one of the factors preventing the adoption of NMPC in many industrial processes.2 1.1. Equation-Oriented Modeling. In equation-oriented modeling, the mathematical representation of the model is usually decoupled from its usage. Moreover, in this approach, the equations of the system are solved simultaneously. This facilitates the use of the same model in several highly sophisticated algorithms for different purposes, such as numerical integration or optimization.3 Consequently, this approach is well-suited for dynamic simulation, optimization, and state estimation. This modeling strategy is also particularly useful when dealing with recycle streams, since the simultaneous resolution of the system of equations allows faster convergence rates.4 Equation-oriented modeling can be achieved using one or a combination of different tools, which can be classified as (1) general-purpose mathematical languages, (2) general dynamic modeling frameworks, and © XXXX American Chemical Society

(3) phenomenological dynamic modeling frameworks. Mathematical models can be developed with general-purpose mathematical languages where equations are manually derived by the modeler, who has full control over all the aspects of the model’s assumptions and implementation. However, for industrial-scale processes, where the process complexity is usually high, it may require typing and handling tens to hundreds of thousands of equations, each containing a few variables, resulting in an extremely large and sparse system of equations.5 In general-purpose mathematical languages, there are no equation-based engineering tools that aid in the development of the physical/chemical model, and the modeler must be aware of the computational impact of the implemented equation format. It may require high-level skills in modeling dynamic processes, mathematics, and computer science. Implementing large and complex dynamic models using this approach can be an extremely time-consuming and error-prone task.5−7 Some examples of widely used general-purpose mathematical languages for modeling are GNU Octave,8 Julia,9 Mathematica, MATLAB, Maple, and Scilab, to name but a few. Dynamic models can also be developed using general dynamic modeling frameworks, with their own simulation oriented language, which supports structured and declarative Received: Revised: Accepted: Published: A

August 14, 2016 February 17, 2017 March 2, 2017 March 2, 2017 DOI: 10.1021/acs.iecr.6b03110 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

1.2. Matching Algorithms. In some modeling systems, such as Simulink, the flow direction of information between model variables is established by the modeler through causal equations, i.e., expressions explicitly assigned to a variable. This is commonly called causal modeling.41 Alternatively, in an acausal modeling environment (also known as noncausal), such as Modelica and PSE gPROMS, there may be no assumptions on the direction of information flow. Equations in these environments can be reused for the evaluation of multiple properties in different contexts. Acausal models are commonly associated with very large and sparse differential DAE systems. In large causal models, it is possible to significantly reduce model size and sparsity resulting in an ordinary differential equation (ODE) system. Therefore, causal models are commonly easier to solve by standard numerical integration algorithms. This is why acausal models are typically transformed to causal models. Acausal equations are commonly symbolically manipulated and sorted using matching algorithms, in a process known as causalization, in order to generate the causal model. The task of matching algorithms is to associate each equation with a variable, commonly assuming that all equations are required in the final model, so that the system of equations is transformed into a block lower triangular form and, hence, solvable by a forward substitution process. Although it is not always possible to transform an acausal model into a fully causal one, the size of large acausal models can often be greatly reduced. See Frenkel et al.42 and references therein for more in-depth information on matching algorithms. 1.3. Automatic Index Reduction. The numerical and algebraic difficulty of solving DAEs can be measured by several concepts of indexes. The higher the index is, the harder it is to solve the system. The differentiation index is one of the most widely used index concepts. It is the minimum number of differentiations, with respect to time, that all parts of these DAE systems must undergo to reach an equivalent ODE system.43 A zero differentiation index system is an ODE system while an index−1 system requires one differentiation to produce an ODE system and contains a structurally nonsingular subset of algebraic equations. Common numerical integration methods for high-index systems may have convergence difficulties or may even provide wrong results.43 In order to overcome these problems, several methods have been developed by the scientific community to automatically reduce the index of high-index DAE systems.44−53 Further information on these index reduction algorithms can be found in Section S2 in the Supporting Information. 1.4. Algorithmic Differentiation. An important feature of any modeling framework is the ability to determine numerical values for the model and its derivative information. Most modeling environments allow the evaluation of differential information using automated strategies, such as finite-difference derivative approximation, symbolic differentiation, and algorithmic differentiation (AD). Daedalus, which is the modeling framework presented in this paper, uses AD, also known as automatic differentiation, which is a family of computational techniques that automatically performs differentiation of a computer call sequence. Additional information on AD can be found in Section S1 in the Supporting Information. 1.5. Scope. This paper presents a modeling library, implemented in C++, which was designed to streamline and accelerate the development process of dynamic lumped models. The Daedalus Modeling Framework builds on principles of earlier environments, such as the Foresee dynamic modeling

representation of process models. These languages are commonly designed to provide symbolic model manipulation, such as differentiation index reduction of differential algebraic equations (DAEs), elimination of variables/equations that can be made implicit, symbolic, or algorithmic differentiation, and an explicit support for dimensional analysis and units of measurement. Although these languages can be used for the dynamic modeling of chemical processes, they do not assist the modeler in the translation of phenomenological knowledge into a mathematical representation and, consequently, they still require the definition of all equations and variables. Nevertheless, implementation of the mathematical dynamic models can be easier than in general-purpose mathematical languages. These languages do not impose any specific ontology onto the modeler, and, therefore, they are typically not associated with any specific physical/chemical conceptualization. However, there can be special libraries or modules that provide useful modeling structures for some phenomena. Some examples of process modeling oriented frameworks are ASCEND,10,11 Aspen Custom Modeler, DAE Tools,12 Dymola,13 EMSO,14 ICAS-MoT,15,16 JModelica.org,17 Mosaic,18,19 Omola,20 PSE gPROMS, 21 the framework developed by Henao and Maravelias,22 and a framework using a Unit-Port-Conditioning Stream (UPCS) approach.23,24 Phenomenological dynamic modeling frameworks explicitly use first-principles to generate equation-oriented models that translate physical/chemical process knowledge into a mathematical representation in a semiautomated process. The focus shifts from manual definition of equations and variables to a theoretical description of the system by specifying relevant phenomena, such as mass and energy transport or chemical reactions. This allows models to be easily modified and reused by reorganization and concatenation of different phenomena. Furthermore, each phenomenological concept can significantly condense model information while also making it easily understood.25 These languages/environments are designed from the ground up for a particular application domain, such as process system engineering for chemical processes. They identify different concepts in their knowledge domain, such as chemical components and chemical reactions, and their taxonomic relationships.26 The modeler would then use those concepts to define the model by specifying information such as the system’s topological structure, underlying assumptions, and simplifications. Such frameworks can provide, for example, stoichiometric validation of chemical reactions and an automated system for the creation of mass and energy balance equations. Phenomenological modeling makes it easier to perform many of the subtasks associated with process modeling identified by Preisig.27 However, phenomenological modeling frameworks can have a limited range of applicability when compared with general-purpose mathematical languages and general modeling frameworks. Stephanopoulos et al.28 implemented the first process modeling language using these concepts in Model.la, which was later improved by Bieszczad.29 ModKit26,30 is another phenomena-oriented language which, just as Model.la, can generate PSE gPROMS source code. Fedorovas’ template-based modeling framework,31 ForeSee,7 Mobatec Modeller,32 ModDev,33 Modeller,34 OntoCAPE,35,36 Ontology Editor,37 OntoMODEL,38 TechTool,39 and the modeling framework of Preisig,40 are additional examples of phenomenological modeling environments. In order to simplify model development, many of these frameworks have specialized graphical user interfaces (GUIs). B

DOI: 10.1021/acs.iecr.6b03110 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Chart 1. Example of a C++ Code Excerpt for the Generation of an ODE System with the Daedalus Modeling Frameworka

a The data given in this chart depict a perfectly mixed tank with an inlet and an outlet stream containing a solution of a component A. Legend: m denotes the mass of the tank, cA is the concentration of A, inFlow is the inlet flow rate, cAin is the concentration of A in the inlet, and outFlow is the outlet flow rate. D__m and D__cA represent the time derivate of the mass in the tank and the time derivative of the concentration of A, respectively.

system7 and Modelica.13,17 It uses a new matching algorithm that allows the definition of equation pools for each model property. The pools represent alternative evaluation paths, from which the framework can automatically choose the equation to determine each variable that is required throughout the model. Each equation and property is a candidate for inclusion in the mathematical representation of the process. They will be present in the model if they are selected by the matching algorithm. This avoids manual selection of equations and variables by the modeler in order to create a well-defined system of equations. One should emphasize that this approach does not involve branches found in conditional modeling,54,55 where different equations are selected based on user-defined conditions at runtime. Equations and variables are selected during the model structure definition stage. The matching algorithm assigns costs to each mathematical operation and function in order to generate a model that is as simple as possible. For instance, when formulating the model to compute the mass fractions for several chemical components from the component masses, the last fraction calculation can be formulated as the difference between unity and the summation of the remaining fractions rather than through the division of masses. Since the generated system of equations is not overdefined, its solution does not require least-squares solvers. Furthermore, equations in Daedalus can also be classified as higher-order functions. Similar to model properties, they can be shared, overridden, and disabled.41 A hybrid AD library, CppADCodeGen,56 was developed to reduce derivative evaluation runtime for Daedalus models. This library also features three structural index reduction algorithmsthe Pantelides method,45 the Soares−Secchi method,53 and the Dummy Derivatives method47and it is capable of simple symbolic manipulations. Furthermore, it also aids in the creation of human-readable documents for the system of equations generated by Daedalus.

1.6. Outline. A generic overview of Daedalus, including its design concepts and implementation strategy, is presented in Section 2. The ontology defined in Daedalus for its phenomena-oriented modeling approach is presented in Section 3. Section 4 describes the graph-based model structure evaluation algorithm that aims to select the best equations to include in the model. Two implementation examples are discussed in Section 5. The paper ends with some conclusions and ideas for future work.

2. DAEDALUS MODELING FRAMEWORK The Daedalus Modeling Framework is implemented so that it would be relatively easy to use by a chemical engineer with some modeling and programming background, without the need of learning yet another domain-specific language. For a comparison between the use of domain-specific and general programming languages for process modeling, see the work of Nikolić.12 The framework aids chemical engineers to create their models in a C++ programming environment. See Chart 1 for a small example of a C++ model file that creates an ODE system that describes a perfectly mixed tank with an inlet and an outlet stream containing a solution of a component A. Daedalus has the following major objectives in mind: (1) to use a phenomena-oriented modeling strategy; (2) apply object-oriented modeling principles; (3) allow both causal and acausal equations; (4) allow the specification of multiple equation alternatives in a context-aware modeling environment; (5) generate distinct systems of equations for different input/ output variables using the same higher-level model description; (6) produce a highly efficient library for the model and its derivative evaluation; and (7) generate human-readable documents for model inspection. C

DOI: 10.1021/acs.iecr.6b03110 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

In addition to the documentation files used to display and navigate through the system of equations, Daedalus can also create a graph with the phenomenological components and their connections, which were used to define the model. When requested, Daedalus performs a structural index analysis and reduction of the generated DAEs model. It is noteworthy that the models generated by Daedalus have a fixed structure. Therefore, methods that determine near-index problems, 46 that would change the system along the integration, are not considered. Three algorithms using the structural information from the equations were implemented: the Pantelides method,45 the Soares−Secchi method,53 and the Dummy Derivatives method.47 The first two methods analyze the structural index of the DAE system and determine which equations must be differentiated and then added to the original system of equations in order to reduce the index, which results in an overdetermined system. The Dummy Derivatives method uses the Pantelides method to determine the index and equations to be differentiated and then generates a welldetermined system of equations. In Daedalus, the user can also select the Soares−Secchi method as a alternative to the Pantelides method used by the Dummy Derivatives method. Further reading on these index reduction algorithms and the algorithms implemented in modeling framework presented here can be found in Section S2. 2.1. Model Generation. The model generation process goes through several steps, which are depicted in Figure 1. First,

The first five objectives are meant to simplify the development of large and complex models. Like other frameworks using phenomena-oriented process modeling, Daedalus automatically generates basic equations, such as those resulting from mass/energy balances, reaction and transport rates, equilibrium equations, etc., which allows rapid prototyping of equation-based process models.7,26,28−30,32,39 The object-oriented programming paradigm can significantly help the modeling process by structuring large models with simpler modeling concepts. It also helps to improve the reuse and extensibility of existing source code. Acausal equations further simplify model development by allowing a single equation to be included in the equation pool of multiple properties through an automatic manipulation of expressions. This reduces the total number of causal equations that would have to be typed for each of those properties. Additionally, the flexibility to use the same phenomenological model definition to generate distinct systems of equations for different sets of output variables also reduces model maintenance. For instance, in the context of NMPC, the same model can be used to create a state-space DAE system and the observations model. Note that the Daedalus function is to generate source code of the models and that external/third party software is required in order to set up and run simulation, control or optimization tasks. For instance, it is currently used to generate models for the Plantegrity nonlinear model predictive control system but it does not include any NMPC ad hoc capabilities. Daedalus represents a system of equations through an acyclic operation graph that can be used for numerical evaluation purposes. Since Daedalus often produces models that contain a high number of temporary variables that are reused in several expressions, the AD presents itself as a well-suited approach to determine numerical derivative values. Differential information on Daedalus-generated models can be computed using one of the following AD packages: ADOL-C,57 CppAD,58 and CppADCodeGen.56 The first two packages use an Operator Overloading (OO) AD strategy, while the last package uses a hybrid AD strategy, which generates source code. Introductory information on AD and additional details on the CppADCodeGen package is given in Section S1 in the Supporting Information. All of these AD libraries use a special class for differentiable model variables, also known as active variables, so that it is possible to register all operations and intrinsic mathematical functions using the operator overloading and function overloading capabilities of C++. In order to use each package, this computer-aided modeling system retraces the operation graph using the required active variable types. In addition to generating C source code through CppADCodeGen for the numerical evaluation of models and their differential information, Daedalus also produces documentation-oriented sources. The package CppADCodeGen56 is used to create LATEX sources, which can be used to compile PDF files, and HTML documents with presentation MathML markup and Javascript. The idea of automatically generating LATEX documentation for a model is not new; it has been implemented, for instance, in BIMAP,59 while MathML has been used, for instance, by the MOSAIC system.18 Since the documentation expressions should leave no room for misinterpretation due to typesetting and variable naming, mathematical notation recommendations from Kuntsche et al.18 were closely taken into consideration. Section S4 in the Supporting Information presents an example of generated documentation.

Figure 1. Typical model generation procedure in the Daedalus Modeling Framework.

it involves the creation of a C++ file that defines the model. Chart 1 presents an example of a phenomenological model description of a perfectly mixed tank. A complete model definition for a CSTR model is presented in Listing S1 in the Supporting Information. The modeler starts by specifying the independent model variables that must be provided to the model in order to compute the desired dependent variables/ D

DOI: 10.1021/acs.iecr.6b03110 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research equation residuals. In Chart 1, the independent variables are m, cA, inFlow, cAin, and outFlow. The modeler then translates her/his phenomenological knowledge of the process into C++, using object-oriented modeling principles. No order is enforced for the definition of phenomena in the source code. These files are then used to request the set of dependent variables and/or equation residuals from the model that must be included in the final system of equations. In Chart 1, the dependent variables are D__m and D__cA. The next step in the model generation process is the compilation of the C++ files into a library, which is then loaded by the Daedalus Modeling Framework. This is followed by the creation of a set of interconnected nodes in the form of an equation/property digraph representing the internal model properties and equations, which are registered in repositories. These repositories are used by the framework to produce an operation graph that defines the mathematical structure of the system of equations used to create the model. The resulting model structure is also known as the flattened model.60 This flattening process is also used in general dynamic modeling frameworks such as JModelica.org. In Daedalus, the choice of the independent variables (exogenous model properties) and the model dependent variables has influence on the evaluation of the model structure. Note that, even for the same phenomenological model description, very distinct model structures may be created by Daedalus, based on the selection of dependent and independent variables made by the modeler. Daedalus only uses the equations and properties from the repository, which are strictly needed to evaluate the requested dependent variables and equation residuals without any intervention from the user. It is noteworthy that the framework generates a well-defined system of equations without the need for a least-squares method to solve the system. Optionally, the operation graph for the flat system of equations can still go through a structural index analysis and reduction in order to produce a new operation graph for a system of equations, which is easier to solve by conventional numerical integration and optimization algorithms. Of particular interest for NMPC formulations is the model computation time. Since the computational effort for the NLP minimization and the model integration can be dominated by the calculation of model derivatives, in Daedalus, a reduced index operation graph can be used with two existing OO AD libraries and the hybrid AD library CppADCodeGen. The OO AD libraries represent the model through a tape of operations while CppADCodeGen generates C source code, which is compiled at a later stage (see Section S1 in the Supporting Information). Daedalus generated models can also be used in MATLAB and GNU Octave via MEX-functions. Furthermore, operation graphs are also used to produce human readable documentation essential to the modeling process. 2.2. Object-Oriented Modeling. Daedalus’ phenomenological model description uses an object-oriented paradigm that allows an abstraction centered around data organization, access, and modification. The framework allows one to define process models via the specification of classes which are used to organize and to describe real-world concepts, such as a chemical substance, a reactor, or a heat capacity equation. Most general and phenomenological dynamic modeling languages have some features of the object-oriented programming, such as composition and inheritance.7,11,13,14,17,32,39

Two fundamental object-oriented programming concepts are encapsulation and polymorphism. As stated by Linninger and Krendl,39 however, it is difficult to take advantage of encapsulation features since most class attributes, such as model properties and equations, are commonly fully accessible. Nevertheless, encapsulation and polymorphism can improve the reusability of modeling components through the definition of interfaces. Just as in the GML language used by TechTool,39 Daedalus supports deep interfaces, where equations and properties defined internally by interface implementations are considered to be included in the full model. Other modeling components do not need to have access to internal implementation details beyond variables and equations defined in the interface. For example, a generic interface for the evaluation of activity coefficients can have multiple implementations. In the proposed modeling framework, a model class can be derived from multiple superclasses (multiple-inheritance), where each particular parent class focuses on a different modeling concept. For instance, the class “FlashStreamPhaseSplitElement”, which implements a zero-holdup flash splitter, is a subclass of “StreamPhaseSplitterElement” and “IntensivePhaseEquilibrium”. While the “StreamPhaseSplitterElement” introduces mass and energy conservation equations, the “IntensivePhaseEquilibrium” defines equilibria relations between chemical component compositions in the liquid and vapor stream phases. Class inheritance can also be used to change the original behavior of a superclass by overriding its methods. The class “FlashStreamPhaseSplitElement” overrides a method in order to indicate that the outlet stream phases are at the equilibrium conditions. 2.3. General Dynamic Modeling. In order to translate the phenomena-oriented model into a system of equations, concepts from general dynamic modeling languages are also included in Daedalus. In the following sections, an ontology used to define generic mathematical concepts is presented. 2.3.1. Properties. Internal model variables, such as temperatures and pressures, are denoted as properties. Each property instance has a unique identifying name and a dimension. Unique variable names can be generated for each property, according to the type of output document for which the model is generated. For instance, while it is important to have human readable names in LATEX and HTML documents, short names improve the compilation time of C files. In addition, each property also contains a pool of functions. Each of these functions constitutes an alternative evaluation path that can be used to determine the property. Functions are nothing more than a particular symbolic manipulation of a model equation and can be seen as directed edges from an equation to a property in the property/equation digraph. A node for the flat operation graph is also associated with a property once it is defined either by the modeler or by using one of the available functions. Properties can be shared across several model objects. For instance, a well-mixed element and its outlet streams can use a single temperature property, hence reducing the need to add additional trivial assignment equations. 2.3.2. Equations. Daedalus allows both causal and acausal equations. While acausal equations can be useful to reduce the total number of equations, the restriction on the flow of information in causal equations reduces the complexity of the property/equation digraph. An example of a definition of an acausal equation is given in Chart 2, and an example of the E

DOI: 10.1021/acs.iecr.6b03110 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research Chart 2. Acausal Equation Definition for Each Component Concentration in a Phase Class

Chart 3. Causal Equation Definition for Each Component Concentration in a Phase Class

dynamic frameworks. The closest modeling framework, which was also an inspiration to Daedalus, is the ForeSee dynamic modeling system.7 Therefore, this section will focus on differences but without neglecting a short description of the underlying concepts, such as phases and composition, streams, modeling components (namely Equilibrium), and tracking chemical components and phases. 3.1. Phases and Composition. One of the most commonly used classes in Daedalus represents a chemical or biological component. It can contain information such as name, molar mass, volumetric mass equations for each physical state, and diverse thermodynamic properties and equations. One or several components can constitute a substance in a given physical state, which is defined by the “IntensivePhase” class. In this class, component compositions are strictly defined by intensive variables (e.g., mass fraction, molar concentration). Moreover, several acausal equations can be used to convert between each of these properties. Other intensive properties can be defined, such as fugacity and activity coefficients for each component, mixture vaporization enthalpy, etc. As shown in Figure 2, the “IntensivePhase” class can be extended by classes introducing extensive properties, such as the “Phase” and “StreamPhase” classes. Note that, although related, extensive variables in these two classes have different dimensions.

definition of a causal equation is given in Chart 3. As previously stated, acausal equations are also used by other modeling environments, such as Modelica and PSE gPROMS. In addition, in this framework, equations can also be considered as higher-order functions,41 since they are objects that may be moved and shared among multiple instances of modeling classes. Moreover, equations can also be overridden and disabled. For instance, the default phase equilibria equations can be assigned different expressions in order to include Murphree efficiency coefficients for distillation column models. The same effect could be achieved by disabling the original equations and introducing new ones. Some equations are only considered during the matching algorithm if a given set of conditions are met. These conditions can include, for instance, the verification that all user-required information was provided. 2.3.3. Dimensional Analysis. Daedalus also performs a dimensional analysis through the “boost::units” library. The most widely used dimensions, units, and quantities are defined for the model variable types and constant scalar types, which allows the following to be written: Qa_pressure

p1 = 1.0*bar

Qa_pressure

p2 = 100000.0*pascal

Here, the variables p1 and p2 could be used interchangeably without any further model modification. This avoids the errorprone task of unit conversion and dimensional checking. Furthermore, note that, in many situations, dimensional checking is able to prevent physically or conceptually incorrect models. For instance, attempting to create a mass flow rate through the division of a quantity with units of mass by a quantity with units of time would be a valid expression whereas the same expression would produce a compilation error for a molar flow rate. Scaling of model variables and equations in order to improve the conditioning of matrices, such as the Jacobian and Hessian of the model, may still need to be performed. Plantegrity, for instance, scales Daedalus-generated models through userdefined normalization constants and/or by analyzing the Jacobian and Hessian matrices.

3. PHENOMENA-ORIENTED MODELING The ontology for the domain of chemical process modeling is implemented through object-oriented programming by the definition of classes. Many concepts implemented in Daedalus are common to multiple phenomenological modeling frameworks and even “sub-models” implemented using general

Figure 2. Phase classes in the Daedalus Modeling Framework. F

DOI: 10.1021/acs.iecr.6b03110 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

• constitutive equations describing the rates of chemical/ biological transformation or transport of mass and energy, and • algebraic equations, such as constraints associated with the behavior of geometric variables. Each model equation type can be represented by a distinct fundamental modeling component. These are process elements, connectors, coordinators, and equilibria. The ForeSee dynamic modeling system7 associates the first three modeling components with different equation types. For a more in-depth description of the first three modeling components, the reader is referred to the work of Tu and Rinard.7 The phenomena-oriented modeling concepts implemented in Daedalus are based on some of the concepts implemented in this modeling system. However, note that ForeSee only achieves the first two objectives laid out at the beginning of this section. Process elements in Daedalus are primarily associated with conservative equations and can be used to determine the dynamic behavior of a specific physical state. Figure 4 presents

The “Phase” class, representing a mixture of components with a single physical state, contains additional properties to specify a components’ composition: the amount (number of moles) and mass. The quantity of material of a specific component can be either specified directly by the modeler or evaluated using additional equations available in the “Phase” class. The mixture’s total amount of material may be defined using either the total mass, the total amount, or the total volume. The conversion between these different composition variables is a good example where the use of a pool of acausal equations can add significant flexibility to the modeling framework. In order to avoid defining the same information in multiple modeling objects, such as the same volumetric mass equation in multiple instances of the Phase class, each model implementation requires the creation of a single, globally accessible, object of the class “Model”. This allows for default values to be defined and accessible throughout the model. Furthermore, the Model object also holds references to all objects of many different classes, making it possible to perform global modifications and analysis. In Daedalus, several precoded process units, or part of them, may be available through the implementation of a “SubModel” class, just like in block-oriented modeling. These types of model blocks are also present in most modeling environments featuring aggregation. 3.2. Streams. A stream phase is considered to be a flow of material with a single physical state while a stream can contain one or more stream phases in a temperature and pressure equilibrium. In many parts of the library, however, both stream and stream phase objects can be easily handled by using the “BaseStream” as a uniform interface (see Figure 3). Similar to

Figure 4. Process elements in the Daedalus Modeling Framework.

some of the implemented process elements. The Well-Mixed Element (WME) is also a subclass of “Phase”, where all intensive properties are considered to be homogeneous throughout its volume. It can be connected to any number of inlet and outlet streams with the same physical state. A Stream Phase Splitter is a zero-holdup element that separates one stream into two or more single-phase streams. Notable subclasses are flash elements and splitters based on constant split ratios. In these Stream Phase Splitters, total flow rates can be determined using acausal equations for mass balances with the inlet and outlet streams. A similar element allowing multiphase outlets are Stream Splitters. A Stream Mixer is another zero-holdup element. It combines two or more multiphase streams into one multiphase stream. Mixing single-phase streams into another single-phase stream with the same state as the inlets is also possible through a Stream Phase Mixer. Connectors are used to describe constitutive equations with the exception of state equations. These include heat or mass transfer to/from phases and composition transformation in a single phase. Some examples include mass transport connectors, that provide mass-transfer fluxes between phases, chemical reaction connectors, that determine chemical/biological reaction rates within a phase, and heat transport connectors. Coordinators are used to define constraint equations that describe the dynamic behavior of geometric variables of a given phase, such as volume or surface area. Finally, Equilibrium components are used to account for extremely fast dynamic responses between two or more phases.

Figure 3. Relationship between Stream, PhaseStream, and Phase classes.

the “Phase” class, in addition to the intensive properties, component compositions can be defined as extensive partial flows (either mass or molar flow). The total flow of a stream phase can be specified in mass, amount (number of moles), or volume. Similar to the “Phase” class, many automatic conversions between different composition definitions are possible. 3.3. Modeling Components. Most model equations in chemical engineering can be divided into four categories. These include: • equations for conservation laws (based on mass, energy, and momentum), • equations for physical and chemical equilibrium among species and phases, G

DOI: 10.1021/acs.iecr.6b03110 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

product components to a well-mixed element when all reactants are present in the phase.

Modeling components can use many other auxiliary classes defined within this chemical engineering modeling library, such as the “Container” class. Containers describe the physical features of an equipment item, such as a distillation column tray or a reactor, and may contain one or more phases where dynamic process changes occur. Note that the current modeling framework implementation does not yet support the generation of models for distributed parameter systems. Many of the modeling component classes perform model definition validations. For instance, the class “DefaultChemicalReactionConnector” checks the chemical reaction stoichiometry and molar mass of the chemical components. The process element “WellMixedElement” only accepts inlet and outlet streams with the same physical state. 3.4. Equilibrium. Three different equilibria were considered in Daedalus: temperature, pressure, and phase equilibrium. A temperature equilibrium can be assumed when the heat transfer rate between phases is so high that any temperature change is considered to be instantaneous. When a temperature equilibrium is for two or more phases, a single energy balance can be performed globally for all phases, instead of each individual phase. Pressure equilibrium in a given volume with one or more phases can also be considered in the modeling formulation under the assumption that any pressure change is very fast, and the spatial differences due to gravity are negligible. These two types of equilibria typically do not require the addition of algebraic equations to the model. The phase equilibrium can be used when two or more phases transfer mass with each other at a rate fast enough to reach a chemical potential equilibrium almost instantaneously. Phase equilibria typically also implies temperature and pressure equilibria. The implementation considers that fugacity is the same in all phases, which can be determined by several methods (e.g., using equations of state to determine fugacity coefficients or using group contributions method for activity coefficients). While determining equilibrium relations using activity coefficients, it can optionally also determine whether the Poynting factor is relevant enough to be included in the expression. By specifying a phase equilibrium between a set of phases, only the chemical composition of a single phase should be defined since other compositions can be determined through the equilibrium. This type of equilibrium may require the addition of algebraic equations to the model. As an alternative, Newton methods can be employed to iteratively solve these equations within the model (e.g., solving bubble and dew points), without requiring the inclusion of algebraic equations in the model. Since conditional modeling is not yet supported by Daedalus, changing the number of phases while time progresses is not allowed. 3.5. Tracking Components. An important feature in any phenomenological modeling framework is the automatic propagation of chemical/biological components and physical states throughout the flowsheet. Automated definition of the chemical species topology has also been implemented, for instance, by Preisig and co-workers61−63 and in Modeller.34 Collections of distinct modeling objects such as these can be saved and tracked by storing references in Set objects. For instance, whenever a component is added to one of the inlet streams of a stream phase mixer, it is automatically added to the outlet stream. In stream phase splitters, however, the addition a component to an outlet will be dependent on the splitters configuration. A chemical reaction connector will only add

4. MODEL STRUCTURE EVALUATION As mentioned in Section 1.5, an acyclic operation graph, which constitutes the mathematical structure of the flat system of equations, is created using some of the equations and properties defined during the phenomenological model description. Although there could be many equations and properties defined within the framework, only the strictly required ones should be used to create the acyclic operation graph, based on the independent variables provided by the user and the requested model-dependent variables or equation residuals. For instance, there are several equations defined for the evaluation of a molar fraction, but only one of these equations needs to be used, if that molar fraction is required by the model, or none at all otherwise. The selection of an equation and property set and its evaluation order is performed through an analysis of a directed graph composed of nodes representing each equation and property (see Figure 5). While equations nodes (e) receive

Figure 5. Example of a property-equation digraph.

directed edges from every property present in the equation, properties (p) only receive directed edges from the equations for which they can be solved. In the following description, the latter edges will be denoted as functions. The set of functions connected to a property p constitutes a pool of alternative evaluation paths. Therefore, an attempt to solve acausal equations during an initial step is performed in order to create these edges. Because of their nature, causal equations require no symbolic manipulation. In Daedalus, no manual selection of equations and properties is required. Although the modeler has full control over which exogenous properties are model inputs and which properties and equation residuals must be model outputs, the framework can automatically select which additional properties must be evaluated in order to determine such outputs. Intermediate properties are determined only if they are considered as a dependency to a user-requested model output and whenever it is possible to evaluate them with the provided model inputs. Furthermore, considering that there can be several evaluation paths for the same model inputs and outputs, the final flat model structure should contain as few “costly” mathematical operations as possible. Therefore, the final flat model structure is selected such that the total operation cost is minimized for the evaluation of each property. To achieve this goal, individual costs are associated with each mathematical operation and function. For instance, while each addition operation has a cost H

DOI: 10.1021/acs.iecr.6b03110 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

the acyclic operation node is determined with that function and because there are no other unknown properties, it is now possible to determine the full evaluation cost for xa associated with this function. This value becomes the maximum allowed full cost, Cmax, for the evaluation of additional functions for xa. If, during the evaluation of these other alternatives, the cost rises above Cmax, that evaluation is immediately aborted. The matching procedure would then remove the last node from the stack and add the node associated with eq 1, since it has only one unknown property and the lowest shallow cost. In order to compute the full cost associated with using eq 1, property xb must also be evaluated. The property xb is added to the stack and Procedure S1 considers eqs 1, 3, and 6 for the evaluation of xb. Since eq 1 is already on the stack, it would be immediately ignored. The full cost for eq 6 would be computed first, because it does not have any unknown properties. Equation 3 is dependent on the property n which, because of the content of the stack, can only be computed through eq 4. Now it is possible to select the function with the lowest full cost for the evaluation of xb by comparing the full cost of eqs 6 and 3. The latter is considered to have the lower full cost (xb = nb/ (na + xb)) when compared to the former (xb = yb × p/pvap,b * ). While moving up in the stack back to the evaluation of xa, the acyclic operations resulting from using eq 1 is backed up. The full cost associated with eq 1 is also used to reduce the value of Cmax. The next alternative in the pool for xa, associated with eq 8, can now be computed. Equation 8 is dependent on the property n, which, because of the contents of the stack, can only be determined using eq 4 or eq 3. The former alternative is investigated first, leading to an evaluation cost for n related to a single addition. The latter equation requires property xb, which, because of the contents of the stack, must be determined using eq 6. This equation leads to a full cost that rises above Cmax, which also means that eq 3 is too expensive. In a situation where there would be additional unknown properties associated with eq 6 or eq 3, there would be no need to even attempt to compute them. The matching algorithm would backtrack earlier to stack node n. Therefore, the property n must be computed using a backed-up value computed with eq 4. Moving up in the stack to the property xa, it is now possible to select the alternative with the lowest full cost, which is the value computed using eq 8. During this evaluation, the property n was also associated with an acyclic operation node computed with eq 4. If the property xb is also requested as a dependent model variable, Procedure S1 orders eqs 1, 8, and 5 by their shallow costs. All properties required by these equations are defined. Equation 1 provides the evaluation alternative for xb with the lowest cost. The resulting system of causal equations contains only half of the initial set of equations:

of 2, each cosine function has a cost of 50. It is noteworthy that other, more-complex, user-defined strategies to determine the cost of each evaluation path can also be provided. These evaluation paths are explored in the matching algorithm presented in Procedure S1 (Section S3) in the Supporting Information. Procedure S1 is able to determine an acyclic operation node for a property node by assigning a cost to each mathematical operation and function, and searching every conceivable alternative in order to keep the cost as low as possible. This matching algorithm only associates equations to variables if the variable can be made implicit. This procedure does not yet identify algebraic loops. The algorithm saves the order of visited nodes in a stack in order to appropriately identify and avoid graph recursions, i.e., returning to a previously visited node. Note that a previously used equation for the evaluation of a property cannot be used to determine another property. A depth-first/breadth-first hybrid search algorithm was devised to select appropriate evaluation paths. A traditional depth-first search algorithm transverses the graph as far as possible from the model outputs toward model inputs without backtracking. In contrast, in a breadth-first search, each node is visited before going to a lower level. The hybrid algorithm, implemented in Procedure S1, uses a breadth-first search to preorder nodes to be visited in the depth-first search by determining shallow function evaluation costs (costs). A generic overview of Procedure S1 is presented in Section S3. To illustrate the application of Procedure S1, the following small system of acausal equations is considered:

1 = xa + xb na n n 0 = xb − b n 0 = xa −

(1) (2) (3)

0 = na + nb − n

(4)

⎛ p* ⎞ vap,a ⎟ 0 = ya − ⎜⎜ ⎟xa ⎝ p ⎠

(5)

⎛ p* ⎞ vap,b ⎟ 0 = yb − ⎜⎜ ⎟x b ⎝ p ⎠

(6)

where the variables xa, xb, n, na, nb, ya, yb, p*vap,a, p*vap,b, and p are model properties. Now assume that the properties na, nb, ya, yb, pvap,a * , pvap,b * , and p are defined as independent variables. If the property xa is requested as a dependent model variable, Procedure S1 would start by adding the node xa into the stack. It would then order the various alternatives that can be used to compute the property xa: the functions associated with eq 1 (xa = 1 − xb), eq 8 (xa = na/n), and eq 5 (xa = ya·p/pvap,a * ). Since the full cost of these functions, costf, is only known once all of their descendants have been visited, and hence, all mathematical operations are known, the shallow cost (costs) is determined without the evaluation of any additionally undefined property. While determining costs, undefined properties are assigned a high cost. Seeing that the function associated with eq 5 has no unknown variables, Procedure S1 would start by adding the node for that equation to the stack and then evaluating xa. Once

n = na + nb

(7)

na n

(8)

xb = 1 − xa

(9)

xa =

Because of the potentially high complexity that comes from transversing the same node multiple times in more-complex systems, each node keeps a list of previous evaluation failures due to recursion. This information can be used to improve performance by avoiding visitation to the same node under I

DOI: 10.1021/acs.iecr.6b03110 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 6. Schematic representation of the usage of Daedalus to create a model for a continuously stirred tank reactor (CSTR) heated by jacket. “TE” represents a temperature equilibrium and “WME” represents a Well Mixed Element.

and the temperature of the jacket fluid (Tj). The dynamic model has two input variables and 16 parameters. Average times for 100 runs taken to prepare the model for numerical evaluation using the full digraph search (greedy = false) with several AD packages are presented in Table 1. The

similar conditions, which would lead, once again, to an evaluation failure. Backtracking earlier due to evaluation costs being higher than Cmax has the drawback of not being possible to save evaluation failures due to recursion, since not all paths have been searched. It is possible to increase the speed of the matching algorithm by using a greedy selection strategy where any valid path is selected as soon as there are no undefined nodes in the path, thus ignoring selection based on function full costs. In the provided example, this would lead to property xa being evaluated with eq 5. The modeler can also manually select equations for the evaluation of any given property, by inspecting the available alternatives in the generated documentation.

Table 1. Preparation Time Comparison for the CSTR Modela Preparation Time (s)

5. EXAMPLES Two examples are presented to illustrate the application and capabilities of the modeling framework: a simple CSTR with an exothermal chemical reaction and a water−ethanol atmospheric distillation column. A comparison between Adol-C 2.6.0, CppAD 2016, and CppADCodegen 2.0 is also presented. The generated models are applied by simulation in the context of NMPC. The computational platform is an Intel Core i5− 4690K CPU (3.5 GHz) with 16 GiB (1600 MHz) of RAM and a solid-state drive running under GNU Linux (Kubuntu 15.10). The GNU gcc 5.2.1 compiler was used for both the modeling framework and the models. 5.1. Continuously Stirred Tank Reactor (CSTR). A model of a CSTR is used to demonstrate several phenomenological modeling components within the proposed framework. In this model, two aqueous liquid streams feed a reactor while only one stream exits. The flow of the outlet stream from the reactor (u2) is manipulated by a pump. An exothermic zeroorder reaction occurs inside the reactor where a chemical component A is transformed to B. The volumetric flow rate of a cooling fluid (u1) in a jacket is used to control the reactor’s temperature. The masses of the metallic reactor wall and the jacket wall are significant and, consequently, the walls are considered in the model dynamics. The reactor wall is in a temperature equilibrium with the contents of the reactor, while the jacket wall is also in a temperature equilibrium with the jacket fluid. The configuration of this system is adapted from Santos et al.64 Figure 6 shows how Daedalus can be used to model this system while the corresponding C++ source and generated equations can be found in Section S5 in the Supporting Information. The generated ODE system contains four differential state variables: the reactor liquid height (h), the concentration of A (cA), the temperature of the reactor (Tr),

a

Adol-C

CppAD

CppADCodeGen

model structure model conversion source-code generation compilation sparsity evaluation

1.50155 0.00008

1.50155 0.00100

1.50155 0.00110 0.00471 0.52450

0.00014

0.00001

total

1.50176

1.50255

2.03186

Average times obtained for 100 runs with a single thread.

model structure includes the property/equation repository generation, the model structure evaluation, and the DAE structural index reduction tasks presented in Figure 1. By using the greedy selection strategy, the time required to create the model structure can be reduced from 1.502 s to 0.170 s. Model conversion in Table 1 refers to the generation of the model representation required by each AD library from the operation graph created by Daedalus. The Adol-C and CppAD packages make a tape to represent the model, which is available until the end of each simulation run. At the model preparation stage, the Jacobian and Hessian sparsity patterns must be determined so that it is possible to compute numerical values for the sparse Jacobian of the DAE system and the sparse Hessian of the weighted sum of the several equations, henceforward simply referenced as Hessian. The Adol-C and CppAD packages determine these sparsities whenever the tape is created. Consequently, the sparsity patterns are computed once per simulation run for these packages. CppADCodeGen compiles the generated source into a library, which is reused, greatly reducing preparation times after the first simulation run. The source-code generation and compilation tasks are, consequently, performed only once for the model and not once per simulation run. The source-code generation requires the evaluation of sparsity patterns and, therefore, the time required to determine these patterns are included in the source-code generation task. The significant difference between the AD packages using the OO strategy and the hybrid AD is due to the compilation time. J

DOI: 10.1021/acs.iecr.6b03110 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

yi, ysp,i, ui, and uref,i represent predicted output variables, setpoints, predicted control variables, and control references at a time interval i, respectively. Qy,i and Qu,i are diagonal weight matrices. The output variables are the level and the temperature of the reactor. The manipulated variables are the outlet stream from the reactor and the coolant flow rate. Variable bounds, setpoints/references, and weights associated with each variable are presented in Table S1 (Section S6) in the Supporting Information. The length of the predictive (np) and control horizon (nm) is 60 and 40 time intervals, respectively. The time interval duration is 30 s. The resulting NLP is composed of 920 decision variables and 840 equality constraints. This NLP is solved using Ipopt, which is a primal-dual Interior Point solver,69 with the MUMPS linear solver. Each model variable and equation in the optimization problem is scaled according to the data in Table S2 in the Supporting Information. Model values and derivative information are determined with the CppADCodeGen library. Values for each time horizon i are determined in an independent thread with a total of three threads. Figure 7 shows the closed-loop response to a sequence of step changes in the setpoint of the reactor level. The controller takes an average computational time of 0.084 s at each time instant. This computational time also includes printing to screen and saving data in an SQL database. 5.2. Distillation Column. In Daedalus, an atmospheric water−ethanol distillation column with equilibrium stages can be easily modeled as an ODE or a DAE system. The distillation column contains 14 stages, a reboiler, and a condenser. The feed is added at the 11th stage from the top of the column with a liquid a mixture of 70% (mol/mol) water and 30% (mol/mol) ethanol at 82 °C (355.15 K). The column is considered to be a cylinder with a radius of 0.5 m. Both the condenser and the reboiler are cylindrical tanks with a radius of 0.5 m and a height of 1.0 m. Each column tray has a height of 0.3 m with a weir height of 0.2 m and a weir length of 0.6 m. The model for this column was implemented under the following assumptions: • liquid and vapor phases are uniformly mixed; • the vapor holdup is negligible; • the vapor phase is ideal; • there is no pressure drop; and • the stage efficiency is 100%. Figure 8 presents a schematic representation of a stage with uniformly mixed liquid and vapor represented through two well-mixed elements connected by a temperature, pressure, and phase equilibrium.

The times presented in Tables 1 and 2 should not be used to make direct performance comparisons between Adol-C and CppAD, since implementation details may have a relevant contribution to these results. Differences in sparsity evaluation, for instance, are partly related to the need to retrace the model in order to create a weighted sum of the several equations for the Hessians in the case of Adol-C, which is not required for CppAD, since it already contains appropriate drivers. In the context of an application of NMPC, this second-order differential information is used to evaluate the Hessian of the Lagrangian associated with the corresponding NLP. Results in Table 2 show that the runtime of CppADCodeGen is, at least, lower, by 1 order of magnitude. Table 2. Runtime Comparison for the CSTR Modela Runtime (μs) model sparse Jacobian sparse Hessian a

Adol-C

CppAD

CppADCodeGen

14.053 21.937 233.883

2.218 13.270 24.333

0.179 0.546 0.860

Average times obtained for 100 runs with a single thread.

Next, the use of the generated model by Daedalus is illustrated in the context of an NMPC application using the Plantegrity nonlinear model predictive control system, which is based on the NEWCON framework.65 In this work, the NMPC problem is solved using a simultaneous optimization approach. Output, control, and state variables are discretized/parametrized over several time intervals. Within each time interval, state variables are parametrized using third-order Gauss-Radau interpolation polynomials while the control variables are considered constant. Details on NMPC strategies and applications can be found in several review works and books.2,66−68 The objective function in the NMPC control problem of the CSTR is as follows: np

Ψ(u , y) =

nm − 1

∑ (Δysp,i )T Q y,iΔysp,i + ∑ i=1

(Δuref, i)T Q u, iΔ

i=0

uref, i

where Δysp, i = ysp, i − yi Δuref, i = uref, i − ui

Figure 7. Closed-loop response of the CSTR model under ideal NMPC control. K

DOI: 10.1021/acs.iecr.6b03110 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Table 3. Preparation Time Comparison for the Distillation Column Model with the Fast Energy Dynamics Assumptiona Preparation Time (s)

Figure 8. Schematic representation of the usage of Daedalus to create a model for a distillation column stage. “PhE”, “PE”, “TE”, “WME” represents a phase equilibrium, a pressure equilibrium, a temperature equilibrium, and a Well Mixed Element, respectively.

a

how +

CppAD

CppADCodeGen

100.83324 0.00247

100.83324 0.00282

100.83324 0.00651 1.06438 63.67807

0.10153

0.00062

total

100.93724

100.83668

165.58220

Average times obtained for 100 runs with a single thread.

Table 4. Runtime Comparison for the Distillation Column Model with the Fast Energy Dynamics Assumptiona

The distillation column is implemented as a submodel within the framework, being composed of additional submodel objects for each equilibrium stage. Liquid flow rates are determined through the Francis equation for segmental weirs. The liquid height over the weir (how) is approximated by the following smoothing function:70 2

Adol-C model structure model conversion source-code generation compilation sparsity evaluation

Runtime (μs) model sparse Jacobian sparse Hessian a

Adol-C

CppAD

CppADCodeGen

365.478 4063.273 82 434.633

59.181 2196.490 6234.845

20.332 106.553 125.929

Average times obtained for 100 runs with a single thread.

2

how + 0.01

information represents the most time-consuming model evaluations during the nonlinear optimization in the Plantegrity NMPC system. In the second version of the distillation column model, the fast energy dynamics assumption is not considered, leading to the formulation of index-2 DAE system. Stage temperatures become differential states while vapor flow rates and the condenser cooling heat are added as algebraic states. This leads to a high index DAE system with 48 differential states and 48 algebraic states. The Pantelides method differentiates the vapor−liquid equilibrium equations for both components and the vapor molar fraction sum equations, thus adding 48 new equations, and introduces new time derivatives for the vapor molar fractions. The Soares−Secchi method provides the same results as the Pantelides method. By applying the Dummy Derivatives method and symbolic manipulations, 48 equations are made implicit and the time derivatives for the vapor molar fractions and temperatures are transformed into dummy derivative variables. The resulting index-1 DAE system has 32 differential equations and 64 algebraic equations. The design of a model predictive control system requires both steady-state and dynamic model analysis. The former is presented in Section S7.4 in the Supporting Information. The dynamic model analysis of nonlinear models commonly requires the study of open- and closed-loop responses to multiple perturbations. The dynamic behavior of the second model version of the distillation column model is now illustrated in an NMPC application using the Plantegrity system. Just as in the CSTR example, the simultaneous optimization approach is used and state variables are parametrized using third-order Gauss-Radau polynomials. The same objective function is used. Table S3 (Section S8) in the Supporting Information presents parameters used to specify the NLP, such as variable bounds, setpoints, references, and weights in the objective function. The prediction and control horizons have 200 intervals of 10 min each. The resulting NLP has 59 200 decision variables and 58 400 equality constraints, and it is solved by the Ipopt optimizer69 with the MUMPS linear solver. Model information

2 Nonlinear expressions for heat capacities, volumetric masses, water vaporization enthalpy, and ethanol vapor pressure were taken from Perry et al.,71 while water vapor pressure and ethanol vaporization enthalpy expressions were retrieved from the work of Burgess72 and Acree et al.,73 respectively. Component activity coefficients are estimated using a group contribution method: the modified UNIFAC (Dortmund) method.74 Two versions of the column model are considered: a first version, assuming fast energy dynamics, which leads to an index-1 DAE system, and a second version without this assumption, which leads to an index-2 DAE system. A detailed list of the equations present in both versions are available in Section S7 in the Supporting Information. In the version with fast energy dynamics assumption, vapor flow rates are determined by neglecting molar enthalpy changes on each stage, with respect to time. The resulting index-1 DAE system contains 32 differential state variables (molar component holdup in each stage), 48 algebraic state variables (component vapor molar fractions and temperature in each stage), 4 controls (reboiler steam power, distillate molar flow rate, reflux ratio, and rectifier flow rate), and 4 additional process input variables (feed molar flow rate, water molar fraction in the feed, feed temperature, and column pressure). Table 3 presents times taken by Daedalus to prepare the model for numerical evaluation using the full digraph search (greedy = false). As in the CSTR example, the model structure evaluation consumes the largest amount of time when using OO AD packages, whereas compilation is the most expensive task when using the hybrid AD approach. The greedy selection strategy can reduce the model structure evaluation down from 100.833 s to 0.594 s, making OO AD approaches better suited for the initial model development tasks. Table 4 shows average runtimes for the distillation column model with the fast energy dynamics assumption for different AD packages. The hybrid AD leads to lower runtimes, especially for the evaluation of the Jacobian and the Hessian. The evaluation of this differential L

DOI: 10.1021/acs.iecr.6b03110 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 9. Closed-loop response of the binary distillation model under ideal NMPC control.

significantly lower runtimes, but at the cost of a longer preparation stage. As for future work, the performance of the digraph search algorithm can still be improved by caching successful property evaluations. The generation of source code for additional languages and interfaces, such as PSE gPROMS and Scilab, will also be pursued. The exploration of new index reduction algorithms, the exploitation of common equation patterns, support for conditional modeling, and the development of a graphical user interface (GUI) should also be considered for future versions of the modeling framework.

is determined with the CppADCodeGen library. Values for each time horizon are computed in an independent thread with a total of four threads of execution. Figure 9 shows the closed-loop response to a step change in the setpoint of ethanol distillate molar fraction from 0.80 to 0.85, bringing the system closer to the azeotrope. Since it is not possible to simultaneously have all output and control variables at the selected setpoints/references, the controller performs a compromise based on the objective function weights. In practice, setpoints and references should all be updated by an optimization layer using economic information; hence, not only would the molar fraction setpoint change but also all other outputs and controls. The ideal NMPC controller requires an average computational time of 19.162 s at each time instant.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.iecr.6b03110. Algorithmic differentiation; automatic index reduction; matching algorithm; documentation example; CSTR model; NMPC parameters for the CSTR example; distillation column models; NMPC parameters for the distillation example (PDF)

6. CONCLUSIONS AND FUTURE DIRECTIONS Daedalus is a new phenomenological modeling framework used to create lumped dynamic models for chemical engineering applications. It is written in C++ and takes advantage of the object-oriented programming paradigm. The ontology used by Daedalus attempts to define individual concepts for relatively simple and basic engineering laws. Each of these concepts is implemented in an individual class. Through different combinations of instances of these classes, it is possible to model several chemical engineering systems. Each modeling component provides a superset of equations and properties that are be used to generate the flat system of equations. The creation of a flat system is performed through a depth-first/breadth-first hybrid search algorithm of the equation−property digraph. This avoids manual selection of equations and implicit model variables in order to create a welldefined system of equations. The flat system allows one to generate source code such as C, for model execution, and HTML, to support the modeling process. Several algorithmic differentiation (AD) packages were compared. AD packages using the traditional operator-overloading approach exhibited smaller preparation times. The hybrid approach implemented in CppADCodegen achieves

■

AUTHOR INFORMATION

Corresponding Authors

*Tel.: +351 239 798 700. E-mail: [email protected] (J. R. Leal). *Tel.: +351 239 402 578. E-mail: andrey.romanenko@ciengis. com (A. Romanenko). *Tel.: +351 239 798 700. E-mail: [email protected] (L. O. Santos). ORCID

Joaõ R. Leal: 0000-0002-7441-122X Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors would like to thank the reviewers for their insightful comments. João Rui Leal would like to acknowledge M

DOI: 10.1021/acs.iecr.6b03110 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research the Ph.D. Grant SFRH/BDE/51185/2010 from Fundaçaõ para a Ciência e Tecnologia (FCT), Portugal, and the financial support from Ciengis.

Computer Aided Chemical Engineering, Vol. 16; Elsevier: Amsterdam, 2003; pp 209−249. (16) Heitzig, M.; Linninger, A. A.; Sin, G.; Gani, R. A computer-aided framework for development, identification and management of physiologically-based pharmacokinetic models. Comput. Chem. Eng. 2014, 71, 677−698. (17) Åkesson, J.; Årzén, K.-E.; Gäfvert, M.; Bergdahl, T.; Tummescheit, H. Modeling and Optimization with Optimica and JModelica.orgLanguages and Tools for Solving Large-Scale Dynamic Optimization Problems. Comput. Chem. Eng. 2010, 34, 1737−1749. (18) Kuntsche, S.; Barz, T.; Kraus, R.; Arellano-Garcia, H.; Wozny, G. MOSAIC: A web-based modeling environment for code generation. Comput. Chem. Eng. 2011, 35, 2257−2273. (19) Kraus, R.; Fillinger, S.; Tolksdorf, G.; Minh, D.; MerchanRestrepo, V.; Wozny, G. Improving Model and Data Integration Using MOSAIC as Central Data Management Platform. Chem. Ing. Tech. 2014, 86, 1130−1136. (20) Grundelius, M.; Mattsson, S.; Åström, E. Object-oriented components for simulation of adaptive controllers. In Proceedings of the 35th IEEE Conference on Decision and Control, Vols. 1−4; Institute of Electrical and Electronics Engineers (IEEE): New York, 1996; pp 2450−2455. (21) Barton, P. I.; Pantelides, C. C. Modeling of combined discrete/ continuous processes. AIChE J. 1994, 40, 966−979. (22) Henao, C. A.; Maravelias, C. T. Surrogate-based superstructure optimization framework. AIChE J. 2011, 57, 1216−1232. (23) Wu, W.; Henao, C. A.; Maravelias, C. T. A superstructure representation, generation, and modeling framework for chemical process synthesis. AIChE J. 2016, 62, 3199−3214. (24) Wu, W.; Yenkie, K.; Maravelias, C. T. A superstructure-based framework for bio-separation network synthesis. Comput. Chem. Eng. 2017, 96, 1−17. (25) Linninger, A. MetamodelingA novel approach for phenomena-oriented model generation. In Fifth International Conference on Foundations of Computer-Aided Process Design; American Institute of Chemical Engineers (AIChE): New York, 2000; pp 462−465. (26) Bogusch, R.; Lohmann, B.; Marquardt, W. Computer-aided process modeling with ModKit. Comput. Chem. Eng. 2001, 25, 963− 995. (27) Preisig, H. A. Constructing and maintaining proper process models. Comput. Chem. Eng. 2010, 34, 1543−1555. (Selected papers from the Seventh International Conference on the Foundations of Computer-Aided Process Design (FOCAPD), 2009, Breckenridge, CO, USA.) (28) Stephanopoulos, G.; Henning, G.; Leone, H. MODEL.LA. A modeling language for process engineeringI. The formal framework. Comput. Chem. Eng. 1990, 14, 813−846. (29) Bieszczad, J. A Framework for the Language and Logic of Computer-Aided Phenomena-Based Process Modeling. Ph.D. Thesis, Massachusetts Institute of Technology, Cambridge, MA, 2000. (30) Bogusch, R.; Marquardt, W. A formal representation of process model equations. Comput. Chem. Eng. 1997, 21, 1105−1115. (31) Fedorova, M.; Sin, G.; Gani, R. Computer-aided modelling template: Concept and application. Comput. Chem. Eng. 2015, 83, 232−247. (24th European Symposium on Computer Aided Process Engineering (ESCAPE).) (32) Westerweele, M. R.; Laurens, J. Mobatec ModellerA flexible and transparent tool for building dynamic process models. In 18th European Symposium on Computer Aided Process Engineering; Braunschweig, B., Joulia, X., Eds.; Elsevier: Amsterdam, 2008; pp 1045−1050. (33) Jensen, A. K.; Gani, R. A computer aided modeling system. Comput. Chem. Eng. 1999, 23, S673−S678. (34) Westerweele, M.; Preisig, H.; Weiss, M. Concept and design of Modeller, a computer-aided modelling tool. Comput. Chem. Eng. 1999, 23, S751−S754. (35) Morbach, J.; Yang, A.; Marquardt, W. OntoCAPEA largescale ontology for chemical process engineering. Engineering

■

NOMENCLATURE AD = algorithmic differentiation API = application programming interface CAS = computer algebra system DAE = differential algebraic equation MPC = model predictive control NMPC = nonlinear model predictive control NLP = nonlinear programming problem ODE = ordinary differential equation JIT = just in time OO = operator overloading SCT = source code transformation WME = well-mixed element CSTR = continuously stirred tank reactor GUI = graphical user interface

■

REFERENCES

(1) Agachi, P. S.; Nagy, Z.; Cristea, M.; Imre-Lucaci, Á . Model Based Control: Case Studies in Process Engineering; Wiley−VCH: Weinheim, Germany, 2007. (2) Qin, S. J.; Badgwell, T. A. An overview of nonlinear model predictive control applications. In Nonlinear Model Predictive Control; Allgower, F., Zheng, A., Eds.; Birkhauser Verlag: Basel, Switzerland, 2000; pp 369−392. (3) Pantelides, C.; Barton, P. European Symposium on Computer Aided Process Engineering−2 Equation-oriented dynamic simulation current status and future perspectives. Comput. Chem. Eng. 1993, 17, S263−S285. (4) Biegler, L. Simultaneous Modular Simulation and Optimization; Technical Report, Carnegie Mellon University, Pittsburgh, PA, 1983. (5) Tolsma, J. E.; Clabaugh, J. A.; Barton, P. I. Symbolic Incorporation of External Procedures into Process Modeling Environments. Ind. Eng. Chem. Res. 2002, 41, 3867−3876. (6) Linninger, A.; Chowdhry, S.; Bahl, V.; Krendl, H.; Pinger, H. A systems approach to mathematical modeling of industrial processes. Comput. Chem. Eng. 2000, 24, 591−598. (7) Tu, H.; Rinard, I. H. ForeSeeA hierarchical dynamic modeling and simulation system of complex processes. Comput. Chem. Eng. 2006, 30, 1324−1345. (8) Eaton, J. W.; Bateman, D.; Hauberg, S.; Wehbring, R. GNU Octave version 4.0.0 manual: A high-level interactive language for numerical computations; 2015. (9) Bezanson, J.; Edelman, A.; Karpinski, S.; Shah, V. B. Julia: A Fresh Approach to Numerical Computing. ArXiv e-prints 2014. (10) Locke, M. H.; Westerberg, A. W. The ascend-II systemA flowsheeting application of a successive quadratic programming methodology. Comput. Chem. Eng. 1983, 7, 615−630. (11) Piela, P.; Epperly, T.; Westerberg, K.; Westerberg, A. ASCEND: an object-oriented computer environment for modeling and analysis: The modeling language. Comput. Chem. Eng. 1991, 15, 53−72. (12) Nikolić, D. D. DAE Tools: equation-based object-oriented modelling, simulation and optimization software. PeerJ. Computer Science 2016, 2, e54. (13) Cellier, F. Continuous System Modeling; Springer: New York, 1991. (14) Soares, R. d. P.; Secchi, A. R. In European Symposium on Computer Aided Process Engineering−13: 36th European Symposium of the Working Party on Computer Aided Process Engineering (ESCAPE13); Kraslawski, A., Turunen, I., Eds.; Computer Aided Chemical Engineering, Vol. 14; Elsevier: Amsterdam, 2003; pp 947−952. (15) Sales-Cruz, M.; Gani, R. In Dynamic Model Development: Methods, Theory and Applications; Asprey, S., Macchietto, S., Eds.; N

DOI: 10.1021/acs.iecr.6b03110 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

(56) Leal, J. R. CppADCodeGen: Source Code Generation for Automatic Differentiation using Operator Overloading. http://github.com/joaoleal/ CppADCodeGen, 2016. (57) Walther, A.; Griewank, A. In Combinatorial Scientific Computing; Naumann, U., Schenk, O., Eds.; CRC Press: Boca Raton, FL, 2012; Chapter 7, pp 181−202. (58) Bell, B. M. CppAD: A package for C++ algorithmic differentiation. Computational Infrastructure for Operations Research (COIN-OR), 2012. (59) Dieterich, E.; Salden, A.; Schäfer, J.; Schmidt, J.; Eigenberger, G. Computer aided modeling with BIMAP. Comput. Chem. Eng. 1997, 21, 1191−1201. (60) Fritzson, P. Principles of Object-Oriented Modeling and Simulation with Modelica 2.1; Wiley−Interscience and IEEE Press: New York and Piscataway, NJ, 2003. (61) Preisig, H. A.; Lee, T. Y.; Makela, M.; Whittacker, A. D.; Little, F. On the Representation of Life-Support System Models. In 19th Intersociety Conference on Environmental Systems (ICES); SAE Technical Paper Series; Society of Automotive Engineers: Warrendale, PA, 1989; Paper No.89147910.4271/891479 (ISSN 0148-7191). (62) Preisig, H. A.; Lee, T. Y.; Little, F. A Prototype Computer-Aided Modelling Tool for Life-Support System Models. In 20th Intersociety Conference on Environmental Systems (ICES); SAE Technical Paper Series; Society of Automotive Engineers: Warrendale, PA, 1990; Paper No. 90126910.4271/901269 (ISSN 0148-7191). (63) Preisig, H. A. In Advanced Control of Chemical Processes 1994; BONVIN, D., Ed.; IFAC Postprint Vol.; Pergamon: Oxford, 1994; pp 131−136. (64) Santos, L. O.; Afonso, P. A.; Castro, J. A.; Oliveira, N. M.; Biegler, L. T. On-line implementation of nonlinear MPC: an experimental case study. Control Eng. Practice 2001, 9, 847−857. (Advanced Control of Chemical Processes.) (65) Romanenko, A.; Pedrosa, N.; Leal, J.; Santos, L. A Linux Based Nonlinear Model Predictive Control Framework. In II Seminario de Aplicaciones Industriales de Control Avanzado, 2007; p 229. (66) Camacho, E. F.; Bordons, C. Nonlinear model predictive control: An introductory review. In Assessment and Future Directions of Nonlinear Model Predictive Control, 2007; pp 1−16. (67) Diehl, M.; Ferreau, H. J.; Haverbeke, N. Efficient Numerical Methods for Nonlinear MPC and Moving Horizon Estimation. In Nonlinear Model Predictive Control: Towards New Challenging Applications; Magni, L., Raimondo, D. M., Allgöwer, F., Eds.; Springer: Berlin, Germany, 2009; pp 391−417. (68) Biegler, L. Nonlinear Programming; Society for Industrial and Applied Mathematics, 2010. (69) Wächter, A.; Biegler, T. L. On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming. Math. Progr. 2006, 106, 25−57. (70) Balakrishna, S.; Biegler, L. T. Targeting strategies for the synthesis and energy integration of nonisothermal reactor networks. Ind. Eng. Chem. Res. 1992, 31, 2152−2164. (71) Perry, R.; Green, D.; Maloney, J. Perry’s Chemical Engineers’ Handbook, 7th Edition; Chemical Engineering Series; McGraw−Hill: New York, 1997. (72) Burgess, D. In NIST Chemistry WebBook, NIST Standard Reference Database Number 69; Linstrom, P., Mallard, W., Eds.; National Institute of Standards and Technology (NIST): Gaithersburg, MD, 2016; retrieved 01/03/2016. (73) Acree, W., Jr.; Chickos, J. In NIST Chemistry WebBook, NIST Standard Reference Database Number 69; Linstrom, P., Mallard, W., Eds.; National Institute of Standards and Technology (NIST): Gaithersburg, MD, 2016; retrieved 01/03/2016. (74) Weidlich, U.; Gmehling, J. A modified UNIFAC model. 1. Prediction of VLE, hE, and γ∞. Ind. Eng. Chem. Res. 1987, 26, 1372− 1381.

Applications of Artificial Intelligence 2007, 20, 147−161. (Special Issue on Applications of Artificial Intelligence in Process Systems Engineering.) (36) Morbach, J.; Wiesner, A.; Marquardt, W. OntoCAPEA (re)usable ontology for computer-aided process engineering. Comput. Chem. Eng. 2009, 33, 1546−1556. (Selected Papers from the 18th European Symposium on Computer Aided Process Engineering (ESCAPE18).) (37) Elve, A. T. Ontology Design for Representation of mathematical models. M.Sc. Thesis, NTNU, Trondheim, Norway, 2015. (38) Suresh, P.; Hsu, S.-H.; Akkisetty, P.; Reklaitis, G. V.; Venkatasubramanian, V. OntoMODEL: Ontological Mathematical Modeling Knowledge Management in Pharmaceutical Product Development, 1: Conceptual Framework. Ind. Eng. Chem. Res. 2010, 49, 7758−7767. (39) Linninger, A.; Krendl, H. TechToolComputer-aided generation of process models (Part 1A generic mathematical language). Comput. Chem. Eng. 1999, 23, S703−S706. (40) Preisig, H. A. Constructing an ontology for physical-chemical processes. In 12th International Symposium on Process Systems Engineering (PSE) and 25th European Symposium on Computer Aided Process Engineering (ESCAPE), 2015; pp 001−1006. (41) Broman, D.; Fritzson, P. Higher-Order Acausal Models. In Proceedings of the 2nd International Workshop on Equation-Based ObjectOriented Languages and Tools (EOOLT), July 8, 2008, Paphos, Cyprus; Linköping Electronic Conference Proceedings, No. 29; Linköping University Electronic Press: Sweden, 2008; pp 59−69, ISBN 978-917519-823-1 (http://www.ep.liu.se/ecp/029/). (Also see Simulation News Europe 2009, Vol. 19, pp 5−16.) (42) Frenkel, J.; Kunze, G.; Fritzson, P. Survey of appropriate matching algorithms for large scale systems of differential algebraic equations. In Proceedings of the 9th International Modelica Conference, Münich, Germany, 2012. (43) Cash, J. R. Review Paper: Efficient numerical methods for the solution of stiff initial-value problems and differential algebraic equations. Proc. R. Soc. London, Ser. A 2003, 459, 797−815. (44) Gear, C. Differential-Algebraic Equation Index Transformations. SIAM J. Sci. Stat. Comput. 1988, 9, 39−47. (45) Pantelides, C.; Gritsis, D.; Morison, K.; Sargent, R. The mathematical modelling of transient systems using differentialalgebraic equations. Comput. Chem. Eng. 1988, 12, 449−454. (46) Chung, Y.; Westerberg, A. W. Solving stiff DAE systems as “near” index problems; Technical Report, Carnegie Mellon University, Pittsburgh, PA, 1991. (47) Mattsson, S.; Soderlind, G. Index Reduction in DifferentialAlgebraic Equations Using Dummy Derivatives. SIAM J. Sci. Comput. 1993, 14, 677−692. (48) Unger, J.; Kröner, A.; Marquardt, W. Structural analysis of differential-algebraic equation systemstheory and applications. Comput. Chem. Eng. 1995, 19, 867−882. (49) Reissig, G.; Martinson, W.; Barton, P. Differential-algebraic equations of index 1 may have an arbitrarily high structural index. SIAM J. Sci. Comput. 2000, 21, 1987−1990. (50) Pryce, J. A Simple Structural Analysis Method for DAEs. BIT Num. Math. 2001, 41, 364−394. (51) Chowdhry, S.; Krendl, H.; Linninger, A. A. Symbolic Numeric Index Analysis Algorithm for Differential Algebraic Equations. Ind. Eng. Chem. Res. 2004, 43, 3886−3894. (52) Kunkel, P.; Mehrmann, V. Index reduction for differentialalgebraic equations by minimal extension. Z. Angew. Math. Mech. 2004, 84, 579−597. (53) Soares, R. d. P.; Secchi, A. R. Structural analysis for static and dynamic models. Math. Comput. Modell. 2012, 55, 1051−1067. (54) Zaher, J. J. Conditional Modeling. Ph.D. Thesis, Carnegie Mellon University, Pittsburgh, PA, 1995. (55) Grossmann, I. E.; Türkay, M. Solution of algebraic systems of disjunctive equations. Comput. Chem. Eng. 1996, 20, S339−S344. O

DOI: 10.1021/acs.iecr.6b03110 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX