Abstract: The study of problems and inconveniences appearing in the Object-oriented (OO) implementation of continuous simulation systems based in Differential-algebraic Equations (DAE) is presented. It was found that the numerical treatment of the equations is crucial to find a design compatible with OO programming practices which lead naturally to implicit schemes. The family of Backward Differential Formulas (BDF) was found particularly appropriate to achieve high levels of software flexibility and reusability. A series of numerical studies were carried out comparing numerical performances with software quality metrics. It was found that BDF implementations improve substantially the software quality, although the computer costs also increase significatively which ultimately calls for pondering the importance of each software characteristic (i.e., modifiability-extensibility vs. calculation time). The equilibrium of this balance is determined by the size of the problem to solve. A utility function is proposed which can be used to determine the optimum choice.
INTRODUCTION
Object Oriented Programming (OOP) supported by proper framework architectures provide features that might significantly improve the design, implementation and maintenance of large codes (Meyer, 2000; Yang, 2008). This methodology has reached many branches of applied computation, such as data bases, software agents, expert systems, computer graphics, communications support, among others (Adetunde, 2009).
Every computer code has two fundamental elements: data and functions. Data are the carriers of the digital information by means of a proper codification, whereas functions (which actually also should be represented through a codification) are procedures responsible for data modifications. In the procedural programming methodology the programs are divided in modules defined by the functions (subroutines, procedures, functions, etc.). In turn, OOP changes radically the principle of program structure: the data become the fundamental elements, whereas functions (called methods) are defined and implemented according to the modularization naturally suggested by the data. In this way, each program module (called object) is responsible for their own data which can only be modified by their own methods and can only be “observed” by the other objects through communication protocols. In addition OOP introduces other modularity concepts (e.g., abstraction, inheritance and polymorphism) that help the programmer in structuring the systems ensuring the future extensibility and modifiability.
Most of the software tools currently available for scientific computation are still structured based in function-driven modularizations, requiring substantial additional work and often being practically impossible to maintain and extend the systems, even for minor changes. The main reason of this is that the scientific software is still implemented using procedural languages (e.g., in FORTRAN) or OOP languages but following essentially procedural designs. Typical symptoms of this kind of troubles is that often large scale numerical problems are performed employing several tools generating partial calculations which afterward are processed using other applications in order to obtain the final results. There are some interesting coding efforts of complex numerical methods applying OOP (Langtangen and Munthe, 2001; Machiels and Deville, 1997; Liu et al., 1996; Kees and Miller, 1999). Most of them conclude that the roll played by the mathematical approach to each particular problem is crucial to achieve good software designs.
Continuous simulation calls for interesting requirements from the programming language and the modeling environment (Aqel, 2006; Mohamed et al., 2008). The diversity of the simulated phenomena, often requiring interdisciplinary approaches, leads to increasingly complex designs covering progressively wider scopes.
The latter clashes with the demands of simplifications, modifications and multidimensional structure, not only from the point of view of the construction of models with structured methodologies (e.g., blocks) but also in the numerical solution of the underlying equations. Formal programming languages have the advantage of univocacy which is a problem of some block diagrams or textual descriptions (Hakman and Groth, 1999). For example, the mathematical expression z = y/3 + 2x, will always generate an executable code yielding the same results independently of the code generating tool. However, classical formal languages usually lead to structural, conceptual and topological information loss in the equation mass. In turn, block diagrams preserve the structural information but often darken the mathematical representation, hiding the local behavior and the conceptual information.
Nevertheless simulation software definitely goes toward OOP languages and Object Diagram Editors (Otter and Elmqvist, 1995), for they have demonstrated wider scopes and represents more naturally the physical reality. Therefore the application of OOP to continuous simulation problems by means of numerical solvers of differential equations combines the structural and mathematical information, without hiding one from the other.
In this article, a study of the problems arising in the application of OOP in continuous modeling and simulation by means of systems of combined Differential and Algebraic Equations (DAE) is presented. The purpose of the study is to assess an appropriate OO design for solving DAE initial-value problems in a generalized approach, in order to facilitate the task of the modeler that often involves a process of progressive model sophistication, modifying, extending and eliminating equations and system variables. It will be shown that in this type of problems involving a strong bias toward functionality relative to data structures it is important to start from mathematical approaches compatible with OO methodology. In the case of continuous simulation this leads naturally to implicit numerical schemes. Within this type of schemes, the family of Backward Differentiation Formulas (BDF) appears as a solid candidate compatible with software flexibility and modifiability requirements.
CONTINUOUS SIMULATION
Generally speaking, a continuous simulation is a numerical representation of temporal-evolving entities based in synchronized processes; therefore time is a single global variable serving as reference of all other variables. This concept which can look obvious to those familiarized to numerical simulation through differential equations is not the only possibility to simulate dynamical systems (Bhatti et al., 2006). There exists also the concept of local time, used in the so called discrete simulation, where each thread of the general process is associated to a different time line, all of which should be synchronized by means of certain defined criteria.
Physical systems involving time dependent variables are generally represented within continuous simulation frameworks by means of ordinary differential equations. Problems involving additionally the space as independent variable are more complex and should be represented with partial differential equations. The target of the present work is the first class of problems i.e., a single independent variable although some of the conclusions can support the treatment of partial differential equations.
Let us consider for example a simple system of two tanks with two pumps transferring water between them (Fig. 1).
The dynamics of this system can be described as a first approximation with two differential equations for the water volumes and two algebraic equations for the flow rates:
| (1) |
where a simple proportional control of the water levels is assumed, aiming to keep a reference volume, in both tanks.
The usual procedure to numerically solve such problem is to apply a classical integration scheme (i.e., Runge-Kutta), for the sake of which the set of equations is written as:
| (2) |
where, is the so called state vector. In this case:
| (3) |
| Fig. 1: | Tanks of regulation by pumps |
During regular modeling processes, the equations become more and more complicated as additional terms and variables are introduced in order to represent more details of the reality. The specialist task is to experiment with the successive sophistications following a series of hypotheses and comparisons against reality. The following are a list of typical modifications usually encountered during the development of dynamical models which complicate the modeler job:
Change a constant by a function of the state variables: During the development of differential dynamical models, it is usual to extend a system introducing dependencies of the control parameters on the state variables. This modification in turn can increase programming complications when the data structure is too rigid. For example, in the base case the model can be extended by postulating a volume dependence of the control parameters, pi (i.e., non-linear control). If change generality is intended, instead of changing the flow rate equations one should change the declaration of from constant to auxiliary variable and then include an additional line with the volume function. Moreover, depending on the way each data type (constant, auxiliary variable, state variable) is communicated, consistency should also be ensured during the calculations in the data transfer between program modules.
Change a constant by a state variable: This modification is similar to the precedent but here a differential equation is added to account for the variation of the parameter. The inconveniences are analogous than before.
Replace a differential equation by an algebraic equation: This change is a very common transformation when working with balance equations of the type:
| (4) |
where the increase and decrease rates of the variable, E φout, y φin, are functions of the state variables.
If the time constant of this equation is much smaller than the other system equations, a quasi-static approximation can be written by zeroing the temporal derivative, that is:
| (5) |
which is an algebraic equation. Modifications of this type change the order
of the system of equations, affecting among other things the time step control.
However, the main inconveniences is not this but the fact that the entire set
of equations should be rewritten all over again to produce the appropriate form
required by explicit schemes, that is
Replace an algebraic equation by a differential equation: This change is the transformation opposite to the precedent. In this case also inconveniences in the time step control and the reimplementation can appear.
Add a term depending on a temporal derivative: This change typically appears when a derivative control term is included. For example in the level control case:
| (6) |
This kind of modification is complicated for it requires rewriting the system
as
Generate a system of N instances: This type of change is usual in the modeling of processes with a large number of similar components. Typical examples are the chemical processes, piping transport and population models. Actually the spatial discretization of partial differential equations generally leads to sets of many similar ordinary differential equations coupled (cells). In the level control case, the extension to N instances is direct, yielding:
| (7) |
In such cases it is very useful for the modeler to make use of the inheritance property which enables encapsulate the common parts of the subsystems, saving implementation work. Clearly the inconveniences associated with all the preceding modifications increase by a factor N with this extension.
SYSTEMS OF DIFFERENTIAL AND ALGEBRAIC EQUATIONS (DAE)
The numerical solution of initial-valued problems modeled by DAE attracted the interest of the numerical community form the last 30 years (Brenan et al., 1996). Many engineering and scientific problems can be naturally modeled with DAE, yielding mostly to balance equations of the general form:
| (8) |
together with a set of auxiliary algebraic equations representing restrictions or constitutive relations:
| (9) |
In order to ensure that real initial-valued problems can be appropriately represented, the specific form of Eq. 8 and 9 should ensure that the temporal derivatives are uniquely determined for every valid set of values of and. In such cases one can always transform a DAE in a system of Ordinary Differential Equations (ODE) of the type:
| (10) |
which can be solved with some of the classical numerical solver schemes (Runge-Kutta, etc.).
The actual transformation of a DAE to the form given by Eq. 10 can be troublesome for various reasons. For instance, DAE often increase their stiffness when they are rewritten as ODE. However, even if this were not the case, the interest in keeping the DAE structure during the numerical calculation arises from the search of efficient software designs, particularly from the point of view of the modeler.
The general objective of the study can be therefore written as a totally implicit non-linear function of the form:
| (11) |
which are often called residual functions and are taken as a combination of
differential and algebraic equations. An Initial-valued Problem (IVP) consists
of determining the temporal functions
There are a variety of numerical methods for solving IVP of DAE which can be classified in three classes: single-step, multistep and extrapolation methods. The efficiency of each type of approximation depends on the specific problem. In what follows, the software-design aspects of a subclass of multistep method, BDF (Backward Difference Formulas) (Gear, 1971), having interesting flexibility properties, will be studied.
Generally speaking, BDF is based in transforming the DAE in a purely algebraic
system by means of numerical approximations. In order to do that let us define
a transformation replacing
| (12) |
This operation defines a set of estimated residuals:
| (13) |
Equation 13 are a set of algebraic equations whose unknown
variables are the state vector in the new time,
Within the described general procedure, different alternatives exist according to the specification of:
| • | The formulas |
| • | The control strategy of the time step |
| • | The numerical method for finding the roots of Eq. 13 |
| • | The data structure and classes to handle the different steps of the general algorithm |
The abstract BDF method presented previously yields to an advantageous calculation
scheme for the model developments, mainly because it solves naturally most of
the problems introduced by the typical modifications appearing during the modeling
of system dynamics. Actually, the basic functional form
ANALYSIS AND OO DESIGN OF DAE SOLVERS
From the point of view of programming, traditional subroutine libraries available to solve differential equations typically present two main weaknesses: rigid interfaces and complexity exponential with the variety of equations to solve. The first weakness affect the users for it is relatively difficult to elaborate models using these libraries, the modeler assuming a double roll: to represent its model mathematically and to adapt the resulting representation to the required interface. The second problem affects the library developers, since the number of programming and maintenance tasks explode.
The mentioned troubles are consequence of deficient abstractions in the design phases. A first problem is that in order to optimize the numerical performance the access to subroutine is often forced to rigid data structures which require the implementation of new subroutines when any of the data structures is changed during the modeling process. An addition inconvenience is that these subroutines often include input variables that describe the data structure.
Although better implementations of numerical methods can be produced using OO programming, one should remind that there always will be a penalty on the performance. Therefore, since numerical simulation often requires complex and highly efficient codes, the performance of high level programming languages will compete with the design improvements achieved from their application. In this section, an analysis of these issues and the results obtained during the design and implementation of a general DAE solver tool for continuous simulation are presented.
Analysis: Essentially a DAE is a set of scalar functions F = {F1 (A, B), F2 (A, B), ..., FN (A, B)}, where A is a set of N functions of the temporal coordinate (real and independent) and B are the corresponding temporal derivatives of A. The methodology BDF consists of transforming the set F in a system of N algebraic equations G {G1 ([X]), G2 ([X]), ..., GN ([X]) = 0, where [X] is a string of N elements representing the values of the functions A evaluated at current time. The transformation F → G is performed assigning to each A and B element an multistep estimator function f ([[X]]) , where [[X]] represents a list of strings [X] increasing its length as time advance (i.e., the elements of the list [X] are the values of A calculated at each time step).
Design: The purpose of the system to implement is the numerical solution of DAE sets, trying to maintain flexibility to support appropriately the development of mathematical models for continuous simulation. Two methods will be implemented and compared: explicit fourth-order Runge-Kutta and generalized BDF. Runge-Kutta methods and BDF method were implemented to solve the explicit and implicit schemes respectively. In what follows, the following design aspects will be omitted for the sake of clarity: constructors, destructors and methods of consultation of attributes. The syntax for attributes declaration, methods and pseudocodes is a relaxed syntax C++ but the schemes can be implemented in any OOP technology (e.g., JAVA, Small talk). If using another OOP technology to implement the proposed schemes, the comparative results will be similar.
Classes for the implementation of implicit schemes
DAE system: This class defines generically DAE sets of the type given
by Eq. 11. Dynamic matrices are used to represent the dependent
variables
| • | “Initial-guess()”: used to obtain the seed for the
solver of |
| • | “Posible-solution()”: stores the roots of |
| • | “Update-history()”: accept a possible solution that pass the consistency test |
| • | “Retrieve-history()”: retrieve the past values of the state vector |
Controller: This class represents the master control responsible for synchronizing the entire calculation process. Controller manages the current time and the time step. The method “Advance()” is used to advance the system a time step beginning from the present time. During the execution of “Advance()” the last calculated values of the state vector is stored in the history list and the current time step is calculated calling the method “Calculate-New-DT()” (Fig. 3).
| Fig. 2: | Temporary dynamic representation of the variables |
| Fig. 3: | Class diagram for implicit schemes |
| Fig. 4: | Class diagram for explicit scheme |
Estimator: This is an abstract class defining generically relations
of the variables of
Constructor implicit system: This is another subclass of “Estimator”
defining the method “G()” which represent the algebraic equations
associated to the residual functions
Root finder: This is an abstract class defining generically the numerical
methods that solve the roots of the algebraic equation
Implicit solver: This class constructs the sequence required to transform the user model equations to an implicit algebraic problem. The method “Calculate()” receive as a parameter an instance of the class “DAE System”.
Classes for the implementation of explicit schemes: Most of the classes designed for explicit schemes are similar to those corresponding to the implicit schemes (Fig. 4). The main differences are the following:
ODE system: This class specifies generically systems of ordinary differential
equations defined by Eq. 10. The functions of the derivatives,
Constructor explicit system: This is a subclass of “Estimator”
defining the method “H()” that represents the algebraic equations
associated to the functions
Explicit solver: This class builds the calculation sequence of the explicit scheme which is performed through the method “Calculate()” that receives as input an instance of the class “Explicit System”.
RESULTS
Numerical study of the model tanks: Consider a water level control of a series of tanks (Fig. 1). This case is appropriate to analyze the advantages of implementing BDF schemes in a modeling process, starting from a simple base case (proportional control) and progressively rising the model sophistication following a circular process of trial and error. Two schemes were implemented following the OO paradigm, BDF and explicit Runge-Kutta. Using both implementations, different modifications of the model were applied, assessing in each case the software flexibility.
The importance of each modification was characterized by the following parameters:
| • | Number of equations affected by a modification (r1) |
| • | Number of equations that should be added to the model after applying the required algebraic steps (r2) |
The software adaptability to each modification was characterized by the following indicators:
| • | Number of additional code lines required (r3) |
| • | Intermediate algebraic steps required to implement a modification to the model (r4) |
A series of two tanks with classical proportional level control was taken as the base case. The dynamics of this system can be described by two differential equations for each water volume and two algebraic equations for the flow rates:
| (14) |
where, vref is the reference volume.
The described base model was extended to series of 3 and 4 tanks and the following modifications were imposed to each extension:
| • | Include derivative control terms |
| • | Introduce a dependence of the control coefficients on the state variables |
| • | Transform the algebraic equations in differential equations |
The initial conditions and the values of the constant parameters in each particular
case are detailed in Table 1. The indicators ri
were calculated in each case, are detailed in Table 2. For
example, include one derivative control term
The general observation is that the BDF implementation reduces substantially the number of code lines and the number of intermediate algebraic steps. However, similarly to the previous cases, this advantage increase the computer costs which calls for pondering the importance of each software characteristic (i.e., modifiability-extensibility vs., calculation time). The equilibrium of this balance is determined by the size of the problem to solve. If the number of equations of the model is high modifiability and extensibility are the preponderant factors, whereas in models with few equations the gain in software implementation does not match the decrease in numerical performance.
To analyze the mentioned balance two metrics are defined that represent the competition between both elements of quality, software metric, S and numerical metric, C:
| S | = | Number of code lines plus number of algebraic steps (r3 + r4) |
| C | = | Real time/calculation time |
| Table 1: | Study of the model tanks - Initial conditions and values of parameters |
| Table 2: | Characterization parameters of the proposed modifications, software metric S and numerical metric C |
| 1BDF: Backsword diffential formulas, 2RK4: Runge-kutta 4 order under score | |
Table 2 shows the calculated values of S and C for each modification of the tank-series problem. Figure 5 shows the dependence of S and C with the number of equations which can be taken as the model size. It can be seen that C is lower in the BDF implementation, whereas the opposite occurs with S.
The conclusion is that the optimum scheme will ultimately depend on the relative importance that the modeler assigns to the metrics C and S. This type of decision is usually solved by means of a utility function representing the subjective usefulness perceived by the users for a given pair (C, S). In the present analysis, a utility function similar to those used in microeconomy is proposed, that is (Henderson and Quandt, 1980):
| (15) |
where, α represents the relative importance of the numerical performance respect to the software flexibility. The best alternative is determined by the scheme having the higher value of U.
Figure 6 shows the optimum decision map calculated for the modeling of the series of tanks, in the plane defined by the parameter α and the number of equations. It can be seen that for α values higher than 10 explicit Runge-Kutta implementations are more convenient for systems of few equations, whereas BDF schemes are better when the number of equations exceed a certain threshold. The critical number of equations separating both alternatives is lower for lower values of α. This means that as the importance software implementation quality increase, the minimum system size for convenience of implicit BDF implementations decrease. For α values lower than 10, BDF implementations is always recommended.
Similarly than in economy, the value of α in particular cases can be estimated by means of the relative impact of the factors C y S on the utility.
| Fig. 5: | Relations between software metric, S and numerical metric, C |
| Fig. 6: | Optimum decision map calculated for the modeling of the series of tanks |
In the software context, α can be interpreted as the number of S units producing the same impact than a C unit on the development work which depends on the developer skills and what type of work requires more efforts and resources.
Implementation of the core model of CAREM NPP: The CAREM Nuclear Power Plant (NPP) has an integrated reactor, the whole high-energy primary system, core, steam generators, primary coolant and steam dome is contained inside a single pressure vessel (Fig. 7). The flow rate in the reactor primary system is achieved by natural circulation. The driving forces obtained by the differences in the density along the circuit are balanced by friction and produce losses, resulting a flow rate in the core that allows for sufficient thermal margin to critical phenomena.
Mathematically, the core model is represented by differential algebraic equations, corresponding to the mass and energy. All equations are expressed with implicit DAE's.
| (16) |
| Fig. 7: | CAREM nuclear power plant-single pressure vessel |
In Fig. 8 it can see the classes implemented for this example. The core model has 6 differential equations and 18 algebraic equations, where several of them are nonlinear. The update process through the implicit scheme is very simple because the changes are local. On the other hand, these equations implemented with the explicit scheme can be difficult, especially in some equations that may have more derivative terms.
| Fig. 8: | Class diagram of core model |
A minimum change to a component linked to another may cause changes for all components of the reactor. Also, there are algebraic equations which the conversion to ODE’s can be causing a double error of approximation.
CONCLUSIONS
The application of OO designs to improve the development of continuous simulation systems was studied. The analysis has tackled the difficulties and complications that affect the work of modelers and library developers, either about the numerical performance or the software implementation.
A class structure providing rapid implementations of continuous dynamic simulation problems was presented, emphasizing the natural representation by means of DAE systems. This approach leads naturally to implicit numerical schemes among which BDF was found compatible with OOP. However, improvements in software flexibility and reusability introduce numerical penalties that should be taken into account in the general balance between flexibility and performance.
It was demonstrated than a metric established in some cases to determine the most convenient scheme, explicit or implicit, regarding the user requirements. The concepts of inherited and encapsulation play an important roll, producing programs much more modular than the procedural versions, thus promoting code reusability, extensibility and maintainability.
A recent research direction on modeling and simulation is to use Model Based Designs (MBD) (Bhatta and Goel, 1996) which helps to accomplish the type of design modifications mentioned with relatively little overhead compared to procedural methodologies. MBD is a mathematical and visual method of addressing the problems associated with designing complex control systems and is being used successfully in many motion control, industrial equipment, aerospace and automotive applications. However highly coupled non-linear systems (e.g., nuclear reactor neutronics, thermohydraulics and control dynamics) still can involve a strong overhead of design and testing, or otherwise increase the computational cost associated to the precision requirements of some engineering problems. Indeed, the use of DAE based models is complementary to MBD and in fact they can be implemented in MBD tools.