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 functiondriven 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 initialvalue 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
temporalevolving 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:
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., RungeKutta), for the sake of which the set of equations is written as:
where, is the so called state vector. In this case:

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, p_{i} (i.e., nonlinear 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:
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 quasistatic approximation can be written by zeroing the temporal derivative, that is:
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 .
In effect, after Eq. 4 is eliminated the variable E is no
longer a state variable, although it will continue to be called to calculate
,
an therefore it should be calculated from the auxiliary algebraic equations
including the new one. These annoyances imply reprogramming of the model, with
the additional charge of associated implementation errors.
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:
This kind of modification is complicated for it requires rewriting the system
as
which imply as in the case 3 the reprogramming of the model.
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:
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 initialvalued 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:
together with a set of auxiliary algebraic equations representing restrictions or constitutive relations:
In order to ensure that real initialvalued 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:
which can be solved with some of the classical numerical solver schemes (RungeKutta, 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 nonlinear function of the form:
which are often called residual functions and are taken as a combination of
differential and algebraic equations. An Initialvalued Problem (IVP) consists
of determining the temporal functions
that satisfy Eq. 11 given the values
in t = 0.
There are a variety of numerical methods for solving IVP of DAE which can be
classified in three classes: singlestep, multistep and extrapolation methods.
The efficiency of each type of approximation depends on the specific problem.
In what follows, the softwaredesign 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
and
in the residual Eq. 11 by functions
and ,
representing multistep estimators of the state vector and its derivatives (by
convention ).
For example, a semiimplicit 2step scheme would be:
This operation defines a set of estimated residuals:
Equation 13 are a set of algebraic equations whose unknown
variables are the state vector in the new time,
which is the calculated knowing its past history.
Within the described general procedure, different alternatives exist according to the specification of:
• 
The formulas
and

• 
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
is invariant under any of the mentioned model changes mentioned in the previous
section. From this perspective, one is tempted to say that implicit schemes
are a “natural mathematics” for OO, the “natural” objects
being the functions .
The implementation of the numerical solution can be easily encapsulated, leading
the modeler to “see” a solver “black box” in which he or
she introduces, adds, changes, eliminates, functions ,
without worrying about rewriting the equations or data communication problems.
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 = {F_{1} (A, B), F_{2} (A, B), ..., F_{N} (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 {G_{1} ([X]), G_{2} ([X]), ..., G_{N} ([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 fourthorder RungeKutta and generalized BDF. RungeKutta 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
and
, the columns being associated to the system variables and the rows representing
the discretized time coordinate (Fig. 2). The functions
are declared on the method “CalculatF()” which receives through parameter
the estimation of and
,
returning the corresponding values of .
Since a single problem can imply several systems of equations linked through
the variables, an association relation 0..* was included between objects of
the same class. Moreover, the following methods were defined to complete the
functionality of this class:
• 
“Initialguess()”: used to obtain the seed for the
solver of
roots,

• 
“Posiblesolution()”: stores the roots of ,
taken as a possible solution of the current state vector which should be
tested for consistency before being accepted 
• 
“Updatehistory()”: accept a possible solution that pass the
consistency test 
• 
“Retrievehistory()”: 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 “CalculateNewDT()” (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
with the state vector and the current unknown state. The abstraction guarantees
the possibility of user implementation as well as the use of libraries of estimators.
The current time step is included as class property. The subclasses “Estimator
DXDT” and “Estimator X” define the relation of the temporal derivative
and the state vector with the past values of the state vector and the corresponding
unknown current value.
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 .
All class that inherits from “Root Finder” should implement the method
“FindRoots()” having the algebraic function “G()” and the
guess values
as input parameters and returning the root vector .
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,
and the auxiliary algebraic equations are declared in the method “Calculate_DXDT()”
which receives as input the values of
and delivers the values of .
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 RungeKutta. 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:
where, v_{ref} 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 r_{i}
were calculated in each case, are detailed in Table 2. For
example, include one derivative control term
(case a in Table 1) for tank 1 in the explicit scheme, the
number of equations affected by a modification is r1 = 2 (
and q_{1} = p_{1} (v_{ref}v_{1}), the number
of equations that should be added to the model after applying the required algebraic
steps is r2 = 0 (replace the original equations), the number of additional code
lines required is r3 = 2 (declaration and initialization of δ_{1})
and the intermediate algebraic steps required to implement a modification to
the model r4 = 3 (steps to obtain the expression of ).
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., modifiabilityextensibility 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 

^{1}BDF: Backsword diffential formulas, ^{2}RK4:
Rungekutta 4 order under score 
Table 2 shows the calculated values of S and C for each modification of the tankseries 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):
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 RungeKutta 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 highenergy 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.

Fig. 7: 
CAREM nuclear power plantsingle 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 nonlinear 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.