versión impresa ISSN 0327-0793
Lat. Am. appl. res. vol.39 no.3 Bahía Blanca jul. 2009
Linearly implicit discrete event methods for stiff ODE's.
Abstract This paper introduces two new numerical methods for integration of stiff ordinary differential equations. Following the idea of quantization based integration, i.e., replacing the time discretization by state quantization, the new methods perform first and second order backward approximations allowing to simulate stiff systems. It is shown that the new algorithms satisfy the same theoretical properties of previous quantization-based integration methods. The translation of the new algorithms into a discrete event (DEVS) specification and its implementation in a DEVS simulation tool is discussed. The efficience of the methods is illustrated comparing the simulation of two examples with the classic methods implemented by Matlab/Simulink.
Keywords Stiff System Simulation. Quantization Based Integration. DEVS
The use of traditional methods (Hairer et al., 1993; Hairer and Wanner, 1991; Cellier and Kofman, 2006) based on time discretization to integrate stiff systems require the use of implicit algorithms since the required step size used by explicit methods is limited by the stability region and the resulting step size becomes inadmissibly small (Cellier and Kofman, 2006).
In fact, numerical integration methods that include in their numerically stable region the entire left half (λ·h) plane (or at least a large portion of it) are necessary for stiff systems integration (Cellier and Kofman, 2006). Only some implicit methods have this type of stability region. Explicit algorithms showing that feature do not exist.
The problem with implicit methods is that they are computationally expensive because in each step they need to use iterative algorithms to determine the next value (usually with the Newton Iteration). The problem becomes critical in applications related to real time simulation, where in many cases performing iterations becomes unacceptable.
An alternative approach to classic time slicing started to develop since the end of the 90's, where time discretization is replaced by state variables quantizacrete time but discrete event systems. The origin of this idea can be found in the definition of Quantized Systems (Zeigler et al., 2000).
This idea was then reformulated with the addition of hysteresis -to avoid the appearance of infinitely fast oscillations- and formalized as the Quantized State Systems (QSS) method for ODE integration in (Kofman and Junco, 2001). This was followed by the definition of the second order QSS2 method (Kofman, 2002), the third order QSS3 method (Kofman, 2006).
The QSS methods showed some important advantages with respect to classic discrete time methods in the integration of discontinuous ODEs (Kofman, 2004), sparsity exploitation (Kofman, 2002), the property of absolute stability, and the existence of a global error bound (Cellier and Kofman, 2006).
In spite of these properties, QSS, QSS2 and QSS3 fail when applied to stiff systems due to the appearance of fast oscillations. They remain numerically stable at the expense of utilizing excruciatingly small step sizes, as any explicit algorithm is expected to require in this situation.
To solve this problem, a first order backward QSS method (called BQSS, after Backward QSS) was proposedinMigoni et al. (2007). The method was able to efficiently integrate many stiff systems. The basic idea of the BQSS method is to use a future value of the states to obtain the quantized value. Yet, whereas the algorithm is implicit, it still can be implemented without Newton iterations. The reason is that each state has only two possible future values. If a state variable currently assumes a value of qj , the next value of that state variable must be either qj + Δqj, or qj − Δqj. Hence all possible next values can be enumerated without an open-ended search. In other words BQSS was the first explicit method for stiff ODEs.
The main drawback of BQSS is that it performs only a first order approximation and accurate results cannot be obtained. Another problem is that BQSS introduced an extra perturbation term that increases the error bound and, in some nonlinear systems, might provoke the appearance of spurious equilibrium points.
This paper presents first a new method that combines the idea of BQSS and linearly implicit integration.
The first order accurate linearly implicit QSS (LIQSS) method follows the idea of BQSS, but avoids the presence of the mentioned perturbation term and the appearance of spurious equilibrium point using a linearly implicit idea to find the state values where some derivatives cross by zero. Then, a second order accurate LIQSS method (called LIQSS2) is defined combining the principles of LIQSS and the second order accurate QSS method (QSS2). This new method solves, up to certain limits, the problem of accuracy.
The work is organized in the following way: Section II.recalls the principles of Quantization-Based Integration and introduces the problems of QSS methods to deal with stiff systems. Then, Section III.introduces the new methods, namely, LIQSS and LIQSS2 and discuss the main implementation issues of them. Section IV.studies the properties of the new methods (legitimacy and error bound). Finally, Section V. introduces some simulation examples and Section VI. presents conclusions and future work.
II. QUANTIZATION-BASED INTEGRATION
This section introduces the principles of QSS methods, shows the problems they have in the simulation of stiff systems, and presents a first approximation to solve these problems (the BQSS method). Finally, a brief introduction to the DEVS formalism is provided and the DEVS implementation of QSS methods is discussed.
A. QSS Method
Consider a time invariant ODE in its State Equation System (SES) representation:
where x(t) ∈ is the state vector and u(t) ∈ is an input vector, which is a known piecewise constant function.
The QSS method simulates an approximate system, which is called Quantized State System:
where q(t) is a vector of quantized variables which are quantized versions of the state variables x(t). Each component of q(t) follows a piecewise constant trajectory, related with the corresponding component of x(t) by a hysteric quantization function so that:
and qj(t0) = xj(t0). Thus, qj(t) only changes when it differs from xj(t) in ±ΔQj. The magnitude ΔQj is called quantum.
B. QSS and Stiff Systems
has eigenvalues λ1 ≈ −0.01 and λ2 ≈ −99.99 which means that the system is stiff.
The QSS method approximates this system as
Considering initial conditions x1(0) = 0, x2(0) = 20, and quanta ΔQ1 = ΔQ2 = 1, the QSS integration performs the following steps:
In t = 0 we set q1(0) = 0 and q2(0) = 20. Then, = 0.2 and = 20. This situation remains until |qi − xi| = ΔQi = 1.
The next change in q1 is then scheduled at t = 1/0.2 = 5 while the change in q2 is scheduled at t = 1/20 = 0.05.
Thus, a new step is performed in t = 0.05. After this step it results q1(0.05) = 0, q2(0.05) = 21, x1(0.05) = 0.01, x2(0.05) = 21. The derivatives are (0.05) = 0.21 and (0.05) = −80.
The next change in q1 is rescheduled at 0.05 + (1 − 0.01)/0.21 = 4.764 while the next change in q2 is scheduled at 0.05 + 1/80 = 0.0625. Hence, the next step is performed at t = 0.0625.
In t = 0.0625 it results q1(0.0625) = 0, q2(0.0625) = x2(0.0625) = 20, x1(0.0625) ≈ 0.0126 and the derivatives coincide with those of t = 0.
This is cyclically repeated until a change in q1 occurs. That change occurs at t ≈ 4.95, after 158 changes in q2 (which oscillates between 20 and 21).
The simulation continues in the same way. Fig.1 shows the evolution of q1(t) and q2(t) through 500 units of simulated time.
Figure 1: QSS Simulation
The fast oscillations of q2 provoke a total of 15995 transitions in that variable, while q1 only changes 21 times. Consequently, the total number of steps to complete the simulation is greater than 16000 (this number is of the order of the 25000 steps needed by Forward Euler method to obtain a stable result).
Evidently, the QSS method is unable to efficiently integrate System (4).
C. QSS2 Method
The second order QBI method uses first order quantization. As it is shown in Figure 2, a first order quantizer produces a piecewise linear output trajectory. Each section of that trajectory starts with the value and slope of the input and finishes when it differs from the input in ΔQi. A formal definition of a first order quantization function can be found in (Kofman, 2002).
First Order Quantizer
Figure 2: Trajectories of a first-order quantizer.
The QSS2 method then approximates a system like (1) by (2) but now, the quantized variables qi(t) follow piecewise linear trajectories and the state variables xi(t) are piecewise parabolic functions of the time.
Like the first order QSS, QSS2 can be represented by a DEVS model. The advantage of QSS2 is that it permits using a small quantum -i.e., setting a small error tolerance- without increasing considerably the number of calculations. In QSS, the number of steps is inversely proportional with the quantum, while in QSS2 it is only inversely proportional with its square root.
QSS2 also exhibits the problem of fast oscillations in the simulation of stiff systems. For instance, if we use the QSS2 method to simulate System (4), it performs 65467 steps (19 changes in q1 and 65448 changes in q2).
D. BQSS Method
The BQSS method is similar to QSS, but qi is always chosen so that xi(t) goes to qi(t).
Basically, given a state variable xj(t), BQSS uses two hysteretic quantization functions: one from below (q(t) ≤ xj(t)) and the other from above ( ≥ xj(t)). Both quantization functions are defined so that they never differ from xj in more than 2ΔQj.
The quantized variable qj is chosen equal to either qj or , according to the direction of (t). When > 0 we use qj = , and viceversa.
When = fj(q, u) depends on qj, it could happen that the sign of the derivative changes when we evaluate fj using each possibility, i.e., > 0 and < 0. Thus, we cannot find a correct value for qj.
However, in that situation, continuity in fj ensures that a value exists, with qj < < , such that = 0.
Thus, the BQSS method sets either qj = or qj = qj and enforces the derivative = 0 adding an extra perturbation term Δfj = .
Then, given the system Eq.(1), instead of simulating a system like Eq.(2), BQSS simulates a system of the form:
where Δf(t) is normally zero, except when the situation described above is found (i.e., when we cannot find a correct value for qj).
BQSS works fine with most stiff systems. Anyway, the presence of the perturbation term Δf(t) increases the error and can cause the appearence of spurious equilibrium points in some nonlinear systems.
Another limitation of BQSS is that it is only first order accurate. We could not find, based on that idea, a second order accurate method.
E. DEVS Formalism
A DEVS model (Zeigler et al., 2000) processes an input event trajectory, and, based on that trajectory and the initial state, produces an output event trajectory.
The behavior of an atomic DEVS model is formally defined by the structure:
- X is input event set, i.e., the set of all possible input event values.
- Y is the output event set.
- S is the state value set.
- δint, δext, λ and ta are the functions that define the model dynamics.
Each possible state s (s ∈ S) has an associated time advance given by the time advance function (ta(s) → ). ta(s) is a non-negative number indicating how long the system remains in a given state in absence of input events.
Then, if at time t1 the state takes the value s1, after ta(s1) time units (i.e., at time t1 + ta(s1)) the system performs an internal transition going to a new state s2 = δint(s1). Function δint (δint : S → S) is called internal transition function.
When the state goes from s1 to s2 an output event is produced, with value y1 = λ(s1). Function λ (λ : S → Y) is called output function. Functions ta, δint and λ defined the autonomous behavior of a DEVS model.
When an input event arrives, the state changes instantaneously. The new state depends not only on the input event value but also on the previous state and the elapsed time since the last transition. If the model arrives to state s3 at t3 and an input event arrives at time t3 + e with a value x1, the new state is calculated as s4 = δext(s3,e,x1) (notice that e < ta(s3)). In this case, we say that the system performs an external transition. Function δext (δext : S × × X → S) is called external transition function. During an external transition no output event is produced.
DEVS models can be coupled and the result of the coupling defines an equivalent atomic DEVS model.
F. DEVS and QSS
Each component of Eq.(2) can be thought as the coupling of two elementary subsystems. A static one,
and a dynamical one
where Qj is the hysteretic quantization function (it is not a function of the instantaneous value xj(t), but a functional of the trajectory xj(·)).
Since the components uj(t) and qj(t) are piecewise constant, the output of Subsystem (8), i.e., dj(t), will be piecewise constant. In this way, both subsystems have input and output piecewise constant trajectories that can be represented by sequences of events.
Then, Subsystems (8) and (9) define a relation between their input and output sequences of events. Consequently, equivalent DEVS models can be found for these systems, called static functions and quantized integrators respectively (Kofman and Junco, 2001).
The second order accurate QSS2 method can be implemented in the same way of QSS. However, the trajectories are now piecewise linear instead of piecewise constant. Thus, the events carry two numbers that indicate the initial value and the slope of each segment. Also, the static functions and quantized integrators are modified with respect to those of QSS so they can take into account the slopes.
III. LINEARLY IMPLICIT QSS METHODS
In this section, we introduce the new Linearly Implicit QSS (LIQSS) methods of first and second order and we discusses the implementation of them as DEVS models.
A. First Order LIQSS Method
The idea of LIQSS is very similar to BQSS. However, when a value qj so that xj goes to it cannot be found, LIQSS tries to find the value for which = 0 instead of adding the perturbation term Δfj to enforce that situation.
In order to illustrate this idea, we shall simulate System (4) from the same initial conditions and quantum than before.
In t = 0, we can choose q2 = 19 or q2 = 21 according to the sign of (t). In both cases, (0) > 0 so the quantized future value of x1 will be q1(0) = 1. On the other hand, if we choose q2(0) = 21 then (0) = −180 < 0 and if we choose q2(0) = 19, it results (0) = 20 > 0 so there exists a point 19 < (0) < 21 in which (0) = 0. The value for (0) can be calculated (exploiting the linear dependence of with q2) as
Then, the state derivatives result: (0) = 0.192, (0) = 0.
The next change in q1 is scheduled for t = 1/0.192 ≈ 5.2083 while the corresponding in q2 is scheduled for t = ∞
Then, the next step takes place in t = 1/0.192 ≈ 5.2083. At this time, x1 = 1 and x2 = 0. After that, q1(5.2083) = 2 (because (5.2083) > 0). If we reevaluate for q2 = 19 and q2 = 21 it results lower than zero in both cases, so the correct value is q2(5.2083) = 19 because in this way x2 goes to q2.
With these values of q1 and q2 the following state derivatives are obtained: (5.2083) = 0.19 and (5.2083) = −80. The next change in q1 is scheduled to t = 1/0.192 + 1/0.19 ≈ 10.47149, while the one in q2 is scheduled to t = 1/0.192 + 1/80 ≈ 5.22083. So, the next step is given in t = 5.22083, when x2 reaches q2.
Calculations follow in this way. Figure 3 show the evolution of q1(t) and q2(t) through 500 units of simulated time. As it can be seen, in this method the fast oscillations of q2 are not present. In this way, q1 changes 21 times and q2 changes 25 times, which totalizes 46 steps (this constitutes a rather decent result for a stiff system).
Figure 3: LIQSS Simulation
B. LIQSS Definition
Given System (1), the LIQSS method approximates it by Eq.(2), where each qj is defined by the following function
where is equal to q(t−) except for the j-th component, where it is equal to and Aj,j is the j,j component of the Jacobian matrix evaluated in , i.e.,
As we shall see now, when Aj,j 0, setting qi = provokes (in the linear case) the situation = 0.
C. Calculation of
We take qj(t) equal to , except that the j-th component is qj. Notice that we use when fj changes the sign between qj and . Thus, an intermediate point exists where fj = 0.
We take (t) equal to , except that the j—th component is .
and the residual
calling the point where fj = 0, we can write
Taking into account that (because and only differ in the j-th component) and considering that fj(,u) = 0, we can solve the previous equations for , obtaining
In order to estimate Aj,j we can use the expression
If fj depends linearly on qj, Eq.(16) gives the exact value of the Jacobian. Moreover, in that case the least term in Eq.(15) results zero and = , i.e., the value of provokes = 0.
In a nonlinear case, we will have ≈ , and ≈ 0. Although the oscillations will not disappear in this case, they will have a low frequency.
This is analogous to what linearly implicit discrete time methods do by solving the implicit equation only for the linear part of the problem, and this is the reason for calling LIQSS to this method.
D. DEVS Implementation of LIQSS
The main difference between LIQSS and QSS is the way in which q is obtained from x. Eq.(10) says that qj not only depends on xj but also on q.
In QSS, the changes in qj are produced when xj differs from ir in ΔQi. Similarly, in LIQSS qi changes when xj reaches qj. However, changes in qi can trigger changes in other quantized variables because of Eq.(10). In the same way, changes in some component of ui can also change some quantized variable.
Following the idea of the DEVS implementation of QSS, i.e., by the coupling of quantized integrators and static functions, but taking into account the mentioned differences, we can obtain a DEVS model for the LIQSS algorithm.
The structure of this implementation is shown in Fig. 4.
Figure 4: Block Diagram of LIQSS
Since the trajectories of uj(t) and qj(t) are piecewise contant, the static funcions are the same than those of QSS.
Then, we only need to define the LIQSS integrators so that they calculate qj according to the definition of LIQSS.
In order to build the DEVS model of the LIQSS quantized integrator, we shall first analyze its behavior.
Let us suppose that in time t the state xj reaches the value of qj with a positive slope ((t−) > 0). Then, the upper and lower quantization functions and qj must be updated (increasing them by ΔQj). In this situation we first try with an output value qj(t) = = qj(t−) + ΔQj.
If, due to the feedback, we receive an input event with dj = (t) < 0, then we are in the situation where we need to calculate the value that provokes = 0, i.e., we estimate using Eq.(13), that in this case takes the form
and we set qj(t) = (supposing that Aj,j = 0, otherwise we just use qj = ).
The value of Aj,j can be easily estimated using the values of (t−), (t), qj(t−) and (t).
It could happen that the integrator then receives (also at time t) another event with a nonzero value (because of the error in the calculation of ). In this case we shall not calculate any further value for qj at time t.
In the opposite case (when xj reaches qj from above) we proceed in an analogous way.
The other case in which qj changes and the quantized integrator must provoke output events is when a change in the sign of the derivative is received. Suppose that xj was going up towards qj = and at time t (due to the change in some quantized variable or input), an event with dj = (t) < 0 is received.
Then, the integrator must send the new output value qj(t) = qj so that xj goes towards qj. In this case, it can also happen that, due to the feedback, a new event with < 0 is received at time t and we are again in the situation where we need to calculate . In this case, we proceed exactly as before, calculating and ignoring any further change of sign at time t.
As before, the case where xj is initially going down to qj is completely analogous to the one described above.
In any situation, after calculating qj, it results easy to schedule the next output event time:
The behavior described for the quantized integrator can be easily translated into a DEVS model.
E. Second Order LIQSS: Basic Idea
The second order linearly implicit method, called LIQSS2, was developed combining the ideas of QSS2 and LIQSS.
The quantized variables of this new method are piecewise linear instead of being piecewise constant and they are chosen in order to verify , i.e., so that xj goes towards qj.
As an example, in Fig. 5 a general trajectory is shown following this idea.
Figure 5: LIQSS2 Trajectories
Analogously to the case of LIQSS and BQSS, it can happen that the sign of changes when we start a new segment of qj(t). Then, an intermediate slope mj exists that makes = 0. In this case, we can also select the initial value qj of the new segment so that = mj, i.e., we can make the state trajectory to run parallel to the quantized trajectory so that no events are generated. Both values, qj and mj, can be easily obtained when depends linearly on qj.
If we simulate System (4) from the same initial conditions but with quanta Δq1 = 0.1 and Δq2 = 0.1 (10 times smaller than before), the LIQSS2 method only performs 59 steps (20 changes in q1 and 39 in q2). The simulation results can be seen in Fig. 6
Figure 6: LIQSS2 Simulation
F. LIQSS2 Definition
Given the system (1), the LIQSS2 method approximates it by (2), where each component qj is defined as:
Note that in Eq.(21), the condition (t) · (t−) < 0 means that an intermediate value mj exists so that = 0. In a linear system this value can be calculated analogously to LIQSS with the expression mj(t) = mj(t−) − (t−)/Aj,j. In a nonlinear case, we shall also obtain an approximate value.
G. DEVS Implementation of LIQSS2
The simulation scheme for LIQSS2 is the same than before (Fig. 4), but now the trajectories are piecewise linear an parabolic.
Since the quantized variable trajectories of LIQSS2 are, as in QSS2, piecewise linear, the static functions are the same of QSS2.
The main difference between QSS2 and LIQSS2 is the way in which the quantized state variable trajectories are calculated in the integrator.
Each segment of the quantized variable trajectory can be characterized by an initial point qj and a slope mqj. In the same way, the state derivative can be characterized by the pair (dj,mdj). Thus, each input and output event of the quantized integrator will carry two numbers.
Let us then analyze the behavior of the resulting DEVS model. Suppose that at time t the state xj reaches either or qj with (t−) > 0. Using the fact that xj is piecewise parabolic, we know both (t−) and (t−). Thus, we set the new segments of and qj with slope mq = (t−) and initial values xj + ΔQj and xj − ΔQj respectively. Then, as (t−) > 0 we select qj(t) = (t).
It could happen that, due to the feedback, at time t we then receive an event with slope mdj(t) = (t) < 0. Thus, an intermediate output slope between the old and the new one exists that provokes the situation mdj = = 0. If Aj,j = 0, this slope can be estimated as
We also calculate the value (t) that makes mdj = (t) = (t):
Thus, we update the slopes of qj and to the value and we output an event with the pair (,.
Like the case of LIQSS, if we receive another event at time t we do not produce any further output event.
Another situation in which qj changes is when an event with a change in the sign of the second derivative is received. Suppose that xj was moving towards qj with > 0 and at time t an event is received with mdj = < 0 (due to the change is some other quantized or input variable).
Thus, we first try with qj = qj(t) and we update the slope mqj = (t−). If, due to feedback, we receive another event at time t with mdj > 0, we must calculate the value that makes mdj = 0, and we repeat what we did from Eq.(22).
Once qj is calculated, we must schedule the next event time. The time to the next event is given by the first crossing of xj with either or qj.This can be calculated as the minimum positive solution σj of the following equations
Similarly to the case of LIQSS, Aj,j can be estimated as:
All these ideas can be easily translated into the DEVS model of the LIQSS2 quantized integrator.
H. PowerDEVS Implementation
PowerDEVS (Pagliero et al., 2003) is a software tool for DEVS simulation. It has a graphical editor that permits building block diagrams of DEVS models. PowerDEVS libraries contain all the blocks needed to implement the QSS methods, including quantized integrators, static functions, source terms and blocks for discontinuity handling.
The DEVS models corresponding to both LIQSS methods were added to the generic quantized integrator that implementes the QSS methods. Now, this block permits selecting among the following methods: QSS, QSS2, QSS3, BQSS, CQSS, LIQSS and LIQSS2.
This block also permits selecting the quantum and the initial state value.
Section V.illustrates the use of these new methods in PowerDEVS.
IV. THEORETICAL PROPERTIES
We shall treat here the most important properties of the LIQSS methods. We shall show first that the methods perform a finite number of steps in a finite interval of time (this guarantees that the simulation time will always advance). Then we shall analyze the stability and accuracy properties.
A. Trajectories and Legitimacy of LIQSS
A crucial requirement of QSS methods is the legitimacy condition, which ensures that a finite number of events occurs in any finite interval of time. The following theorem proves this property for the first order LIQSS method.
Theorem 1. Suppose that function f in (1) is bounded in a domain D × Du, where D ⊂ , Du ⊂ and assume that the trajectory u(t) ∈ Du is piecewise constant. Then,
- Any solution x(t) of (2) is continuous while q(t) remains in D.
- The trajectory q(t) is piecewise constant while it remains in D.
Proof. The proof of (1) is straightforward since, according to (2), the derivative of x is bounded.
For the item (2), in order to prove that q is piecewise constant it is necessary to ensure that it only experiences a finite number of changes in any finite interval of time.
Let (t1,t2) be an arbitrary interval of time in which q(t) remains in D. We shall prove that, within this interval, q(t) has a finite number of changes.
The assumptions of the theorem ensure that f(q,u) is bounded and, taking into account the relation between xj and qj, positive constants exist so that, for t ∈ (t1,t2)
Let tc ∈ (t1,t2) and suppose that . According to (12), this situation cannot be repeated until |xj(t) − xj(tc)| ≥ ΔQj. Thus, the minimum time interval between two discontinuities in (t) is
Then, calling the number of changes of (t) in the interval (t1,t2), it results that
Since u(t) is piecewise constant, it will perform a finite number of changes nu in the interval (t1,t2).
The definition of qj ensures that it can only change when (t) changes or when there is a change in some other quantized or input variable (qi(t) or ui(t)) that inverts the sign of .
In conclusion, changes in qj(t) are linked to changes in some (t) or ui(t). Thus, the total number of changes will be equal or less than the sum of all the changes in those variables, i.e.,
which is a finite number.
Although this theorem is only valid for the first order LIQSS method, an analog result for LIQSS2 can be obtained combining this proof with that of QSS2 (Kofman, 2002).
B. Perturbed representation
The theorical properties of the QSS methods are based in a perturbed representation of the original system (1) that is equivalent to the approximation Eq.(2). Defining Δx(t) = q(t) − x(t) each row of System (2) can be rewritten as:
From (10), (11), (12) and (13) in the definition, it can be ensured that each component Δxi(t) of Δx(t) is bounded by1
where ΔQi is the quantization adopted for xj(t). Thus, the LIQSS methods simulate an approximate system which only differ from the original SES(1) due to the presence of the bounded state perturbation.
C. Global Error Bound and Stability
Given the LTI system
were A is a Hurwitz matrix with Jordan canonical form Λ = V−1AV, the LIQSS(2) approximation simulates the system
Defining the error as e(t) = x(t) − xa(t), and following the procedure of Kofman (2002) and Cellier and Kofman (2006), it results that
where ΔQ is the vector of quantum adopted.
The error bound is twice the error bound of QSS, QSS2 and QSS3.
D. Equilibrium Points
One of the drawbacks of BQSS was the appearence of spurious equilibrium points. Due to the term Δf, Eq.(6) admits equilibrium points also when f(q,u) 0.
However, the only possibility in which LIQSS or LIQSS2 arrive to an equilibrium point is when f(q,u) = 0, i.e., when the quantized variables reach an equilibrium point.
This section is aimed to introduce some examples which show the advantages of the LIQSS methods.
A. Second Order Linear System
In Section II the following system was presented:
with initial conditions x1(0) = 0 and x2(0) = 20
The system was first simulated using LIQSS and LIQSS2 methods with quantum ΔQi = 1. Then, the simulations were repeated decreasing the quantum 10, 100 and 1000 times. The following table shows the number of steps performed for each method using the mentioned quantization:
This table shows that the number of steps performed by LIQSS linearly varies with the quantization, while the number of steps in LIQSS2 grows approximately with the square root of the quantum reduction.
Figure 7 shows the simulation results using Simulink ODE15s with tolerance 10−3 and PowerDEVS with ΔQi = 0.001 (the difference between both methods cannot be appreciated with the naked eye).
Figure 7: Linear Stiff System Simulator
The simulation time could not be evaluated under PowerDEVS using ΔQi = 0.001 (it was too smal in order to be accurately measure). Thus, we compared the simulation times using quantum ΔQi = 0.0001. The simulation in Simulink takes 0.078 seconds (ODE15s, in accelerated mode, with tolerance 10−3), while in PowerDEVS it takes 0.015 seconds.
The global error bound -Eq.(27)- ensures that the error in the LIQSS2 simulation is always less than 2 · 10−4 in x1 and 6 · 10−4 in x2.
Figure (8) shows the simulation error in PowerDEVS with quanta ΔQi = 0.0001. Comparing it with the theoretical bound, we see that the last one is a bit conservative.
Figure 8: Error
B. Van der Pol Oscillator
The problem consists of a second order differential equation proposed by B. van der Pol in the 1920"s, that describes the behavior of nonlinear vacuum tube circuits. It has two periodic solutions, the constant solution, z(t) = 0, which is unstable, and a nontrivial periodic solution that corresponds to an attractive limit cycle. The equation depends on a parameter that weights the importance of the nonlinear part of the equation. The corresponding state equations are:
We fixed the parameter μ = 1000, what gives rise to a stiff problem that is often used as a test problem for stiff ODEs solvers (Enright and Pryce, 1987).
The model was then built in PowerDEVS (Fig.9)
Figure 9: PowerDEVS model of the Van der Pol oscillator
For the simulation, we used initial conditions x1(0) = 2, x2(0) = 0 and quantization ΔQ1 = 0.001, ΔQ2 = 1. We simulated the system with LIQSS2 until tf = 4000. The results are shown in Figure 10. The total number of steps was 2159 (838 in x1 and 1321 in x2). The simulation takes 0.031 seconds.
Figure 10: Van der Pol oscillator simulation
The same system was simulated using different Matlab/Simulink methods. The best results were obtained with ODE15s. In order to obtain similar results to the ones given by PowerDEVS the tolerance error must be set to 10−14 . Using larger tolerances the results were qualitatively similar but with a significant phase error. In this case, the number of steps was 2697 and the simulation takes 0.056 seconds.
It must be taken into account as we compare the results, that the order of LIQSS2 is smaller that the one of ODE15s.
Finally, the system was simulated again using LIQSS2 but with quanta ΔQ1 = 0.0001 and ΔQ2 = 0.1 (ten times smaller). The number of steps performed result 4148(1975 in x1 and 2173 in x2) and the simulation takes 0.047 seconds. Comparing this result with the previous one, it can be seen that the number of steps was increased just in a factor of 2 while the quantized levels were reduced in a factor of ten. This put in evidence the second order nature of the LIQS2 method
We presented two novel QBI methods that are able to integrate stiff systems based on linearly implicit principles. The fact that the methods do not call for iterations permits achieving an important reduction of computational costs when compared with traditional implicit discrete time methods.
We showed that these methods, satisfy the global error bound property of QSS method.
We showed that these methods, satisfy the global error bound property of QSS method.
The methods improve the performance of BQSS by solving the problem of the spurious equilibrium points and by increasing to 2 the order of accuracy.
Future work must be done in order to develop higher order methods (following the idea of QSS3 for instance). It is also important to establish which kind of stiff system can be efficiently simulated with LIQSS methods.
1 The symbols |·| and "≤" denote the componentwise module and inequallity, respectively.
1. Cellier, F. and E. Kofman, Continuous System Simulation, Springer, New York (2006). [ Links ]
2. Enright, W. H. and J. D. Pryce, "Two (fortran) packages for assessing initial value methods," (ACM) Transactions on Mathematical Software, 13, 1&-27 (1987). [ Links ]
3. Hairer, E., S. Norsett, and G. Wanner, Solving Ordinary Differential Equations I. Nonstiff Problems, Springer, 2nd edition (1993). [ Links ]
4. Hairer, E. and G. Wanner, Solving Ordinary Differential Equations II. Stiff and Differential-Algebraic Problems, Springer, 1st edition (1991). [ Links ]
5. Kofman, E., "A Second Order Approximation for DEVS Simulation of Continuous Systems," Simulation, 78, 76-89 (2002). [ Links ]
6. Kofman, E., "Discrete Event Simulation of Hybrid Systems," SIAM Journal on Scientific Computing, 25, 1771-1797 (2004). [ Links ]
7. Kofman, E., "A Third Order Discrete Event Simulation Method for Continuous System Simulation," Latin American Applied Research, 36, 101-108 (2006). [ Links ]
8. Kofman, E. and S. Junco, "Quantized State Systems. A DEVS Approach for Continuous System Simulation," Transactions of SCS, 18, 123-132 (2001). [ Links ]
9. Migoni, G., E. Kofman, and F. Cellier, "Integración por Cuantificación de Sistemas Stiff," Revista Iberoam. de Autom. e Inf. Industrial, 4, 97-106 (2007). [ Links ]
10. Pagliero, E., M. Lapadula, and E. Kofman, "PowerDEVS. Una Herramienta Integrada de Simulación por Eventos Discretos," Proceedings of RPIC'03,, San Nicolas, Argentina 1, 316-321 (2003). [ Links ]
11. Zeigler, B., T. Kim, and H. Praehofer, Theory of Modeling and Simulation, Academic Press, New York, 2nd edition (2000). [ Links ]
Received: October 18, 2007.
Accepted: March 1, 2009.
Recommended by Guest Editors D. Alonso, J. Figueroa, E. Paolini and J. Solsona.