## Servicios Personalizados

## Articulo

## Indicadores

- Citado por SciELO

## Links relacionados

- Citado por Google
- Similares en SciELO
- Similares en Google

## Bookmark

## Latin American applied research

*versión impresa* ISSN 0327-0793

### Lat. Am. appl. res. v.36 n.2 Bahía Blanca abr./jun. 2006

**A third order discrete event method for continuous system simulation**

**E. Kofman ^{1}**

^{1} *Laboratorio de Sistemas Dinámicos. FCEIA - UNR - CONICET. Riobamba 245 bis - (2000) Rosario - Argentina. kofman@fceia.unr.edu.ar*

*Abstract* — This paper introduces a new numerical method for integration of ordinary differential equations. Following the idea of quantization based integration, i.e., replacing the time discretization by state quantization, this new method performs a third order approximation allowing to achieve better accuracy than their first and second order predecessors. It is shown that the new algorithm satisfies the same theoretical properties of the previous methods and also shares their main advantages in the integration of discontinuous systems.

*Keywords* — Hybrid Systems. ODE Integration. Discrete Event Systems.

**I. INTRODUCTION**

Numerical integration of ordinary differential equations (ODEs) is a topic which has advanced significantly with the appearance of modern computers. Based on classic methods like Euler, Runge-Kutta, Adams, etc., several variable-step and implicit ODE solver methods were developed (Hairer *et al.*, 1993; Hairer and Wanner, 1991). Simultaneously, several software simulation packages have been developed implementing these algorithms. Among them, Matlab/Simulink (Shampine and Reichelt, 1997) is probably the most popular and one of the most efficient.

In spite of the several differences between the mentioned ODE solver algorithms, all of them share a property: they are based on time discretization, giving a solution obtained from a difference equation system (i.e., a discrete-time model).

A completely different approach for ODE numerical integration started to develop since the end of the 90's, in which time discretization is replaced by state variables quantization. As a result, the simulation models are not discrete time but discrete event.

The origin of this idea can be found in the definition of Quantized Systems (Zeigler and Lee, 1998). Quantized Systems were reformulated with the addition of hysteresis -to avoid the appearance of infinitely fast oscillations- and formalized as a numerical algorithm for ODE's by Kofman and Junco (2001), where the Quantized State Systems (QSS) and the QSS method were defined.

The following step was the definition of the method of second order quantized state systems (QSS2) (Kofman, 2002), and then both methods were extended to the simulation of differential algebraic equations (DAEs) (Kofman, 2003) and discontinuous systems (Kofman, 2004).

The discrete event nature of these methods make them particularly efficient in the last case, and a considerable reduction of computational costs with respect to the most sophisticated discrete time methods can be observed.

Despite their simplicity, the QSS and QSS2 methods satisfy some strong stability, convergence and error bound properties, and they intrinsically exploit spar-sity in a very efficient fashion.

This paper continues the previous works by formulating the method of third order quantized state systems (QSS3) which permits improving the accuracy of QSS and QSS2 conserving their main theoretical and practical advantages. An additional advantage of QSS3 is that the choice of the quantum becomes less critical than in the lower order methods since it can be adopted in a conservative fashion without affecting considerably the number of calculations.

After a brief introduction recalling the principles of quantization based integration, the definition of the QSS3 method will be introduced. Then, we shall prove that it is legitimate, i.e., that it cannot produce a Zeno-like behavior, and we shall deduce the input-output relationships of the basic components of QSS3 (quantized integrators and static functions). Then, after a brief discussion about the theoretical properties of QSS3, two relatively complex simulation examples will be introduced.

**II. QUANTIZATION BASED INTEGRATION**

**A. QSS-Method**

Consider a time invariant ODE in its State Equation System (SES) representation:

(*t*) = *f*(*x*(*t*), *u*(*t*)) (1)

where *x*(*t*) ∈ ^{n} is the state vector and *u*(*t*) ∈ ^{m} is an input vector, which is a known piecewise constant function.

The QSS-method (Kofman and Junco, 2001) simulates an approximate system, which is called Quantized State System:

(*t*) = *f*(*q*(*t*), *u*(*t*)) (2)

where *q*(*t*) is a vector of quantized variables which are quantized versions of the state variables *x*(*t*). Each component of *q*(*t*) is related with the corresponding component of *x*(*t*) by a hysteretic quantization function, which is defined as follows:

**Definition 1.** *Let Q* = {*Q*_{0}, *Q*_{1}, ..., *Q _{r}*}

*be a set of real numbers where Q*

_{k}_{-1}<

*Q*1 ≤

_{k}with*k*≤

*r*.

*Let*Ω

*be the set of piecewise continuous real valued trajectories and let x*∈ Ω

_{i}*be a continuous trajectory. Let b*: Ω → Ω

*be a mapping and let q*(

_{i}= b*x*)

_{i}*where the trajectory qi satisfies*:

(3) |

*and*

*Then, the map b is a* hysteretic quantization function.

The discrete values *Q _{k}* are called quantization levels and the distance

*Q*

_{k}_{+1}-

*Q*is defined as the quantum, which is usually constant. The width of the hysteresis window is ε and the best choice is taking it equal to the quantum (Kofman

_{k}*et al.*, 2001).

Kofman and Junco (2001) proved that the quantized variable trajectories *q _{i}*(

*t*) and the state derivatives

_{i}(

*t*) are piecewise constant and the state variables

*x*(

_{i}*t*) are piecewise linear. As a consequence, those trajectories can be represented by sequences of events and then the QSS can be simulated by a DEVS model. DEVS (Zeigler

*et al.*, 2000) is a formalism whose acronym stands for Discrete Event System specification.

The mapping of a QSS like (2) into a DEVS model can be done in several ways and one of the easiest is based on coupling principles. A generic QSS can be represented by the block diagram of Fig. 1. That block diagram is composed by static functions *f _{i}*, integrators and quantizers.

Each pair formed by an integrator and a quantizer is called *quantized integrator* and it is equivalent to a simple DEVS model. Similarly, the static functions have DEVS equivalents and consequently, the entire block diagram has an equivalent coupled DEVS which represents it (Kofman and Junco, 2001). These DEVS models can be found in the cited reference.

Figure 1. Block Diagram Representation of a QSS

Some simulation programs -PowerDEVS (Pagliero and Lapadula, 2002) for instance- have libraries with DEVS blocks representing quantized integrators and static functions. Thus, the implementation of the QSS-method consists in building the block diagram in the same way that it could be done in Simulink.

**B. QSS2-Method**

QSS only performs a first order approximation. Due to accuracy reasons, a second order method was proposed which also shares the main properties and advantages of QSS (Kofman, 2002).

The basic idea of the new method, (called QSS2) is the use of first-order quantization functions instead of the quantization function given by (3). Then, the simulation model can be still represented by (2) but now *q*(*t*) and *x*(*t*) have a different relationship. This new system is called Second Order Quantized State System or QSS2 for short.

A first-order quantization function can be seen as a function which gives a piecewise linear output trajectory, whose value and slope change when the difference between this output and the input becomes bigger than certain threshold (Fig. 2)

Figure 2. I/O trajectories in a *First Order* quantizer

In that way, the quantized variable trajectories are piecewise linear and the state trajectories are piecewise parabolic^{1}.

As before, the system can be divided into quantized integrators and static functions like in Fig. 1. However, the DEVS models of the QSS2 quantized integrators are different due to the new behavior of the quantizers. Similarly, the DEVS models of the QSS2 static functions are also different since they should take into account the slopes of the piecewise linear trajectories.

The formal definition of first order quantization functions and the DEVS models associated to the QSS2 integrators and static functions were also developed by the author (Kofman, 2002).

Like QSS, PowerDEVS also implements the QSS2-method. Thus, a continuous or hybrid system described by a block diagram can be simulated with QSS2 by selecting the corresponding method in the integrators.

**C. Properties of QSS and QSS2**

There are properties -proven by Kofman and Juco (2001) and Kofman (2002)- relating the solutions of Systems (1) and (2). These properties not only show theoretical features but also allow deriving rules for the choice of the quantization according to the desired accuracy.

The mentioned properties are stability, convergence and error bound and the corresponding proofs were built based on perturbation studies. In fact, defining Δ*x*(*t*) = *q*(*t*) - *x*(*t*), System (2) can be rewritten as

(*t*) = *f*(*x*(*t*) + Δ*x*(*t*), *u*(*t*)) (4)

From the definition of the hysteretic and the first order quantization functions, it can be ensured that each component of Δ*x* is bounded by the corresponding quantum adopted. Thus, the QSS and QSS2 methods simulate an approximate system which only differs from the original SES (1) due to the presence of the bounded state perturbation Δ*x*(*t*).

The Convergence Property ensures that an arbitrarily small error can be achieved by using a sufficiently small quantization. A sufficient condition which guarantees this property is that the function *f* is locally Lipschitz.

The Stability Property relates the quantum adopted with the final error. An algorithm can be derived from the proof of this property which allows the choice of the quantum to be used in the different state variables.

Finally, the Global Error Bound is probably the most important property of quantization based methods. Consider a LTI system (*t*) = *Ax*(*t*) + *Bu*(*t*) where *A* is Hurwitz and diagonalizable. Let us call (*t*) and (*t*) to the exact solution and the approximate obtained with the QSS or QSS2 method respectively, with (0) = (0). Then, using Theorem 3 of (Kofman, 2005b), it can be seen that the error in the QSS or QSS2 simulation is bounded by

|(*t*) - (*t*)| ≤ |*V*||e(Λ)^{-1}Λ||*V*^{-1}|Δ*q* (5)

where Λ and *V* are the matrices of eigenvalues and eigenvectors of *A* (Λ is diagonal), that is, *V*^{-1}*AV* = Λ and Δ*q* is the vector of quantum adopted at each component^{2}.

Inequality (5) holds for all *t*, for any input trajectory and for any initial condition.

**III. QSS3 METHOD**

In order to obtain a third order approximation, we need to consider not only the first but also the second derivative of the system trajectories. Thus, we can redefine the first order quantizer shown in Fig. 2 so that the output is piecewise parabolic.

Then, given a system of state equations like (1), the QSS3 method will approximate it by (2) where x and q are related component-wise by *second order quantization functions*.

**A. Second order quantization**

Formally, we say that the trajectories *x _{i}*(

*t*) and

*q*(

_{i}*t*) are related by a second order quantization function if

*q*(

_{i}*t*

_{0}) =

*x*(

_{i}*t*

_{0}) and

(6) |

with *t* ≤ *t* ≤ *t _{j}*

_{+1}, and the sequence

*t*

_{0}, ...,

*t*, ... defined so that

_{j}*t*

_{j}_{+1}is the minimum

*t*>

*t*where

_{j}(7)

with and defined as

(8)

(9)

Figure 3 illustrates the behaviour of a second order quantizer.

Figure 3. Input and output of a second order quantizer

**B. Trajectories in QSS3**

The basis of QSS and QSS2 are the trajectory forms that allow the discrete event representation. A crucial requirement is the legitimacy condition which requires that only a finite number of events occurs in any finite interval of time.

Although the definition of a second order quantization function (6) suggests that the components of *q*(*t*) are piecewise parabolic, it should be proven that the sequence *t _{j}* does not have infinite components in a finite interval of time. The following theorem gives sufficient conditions to this property.

**Theorem 1.** *Consider system* (2). *Assume that the input u*(*t*) *is bounded and left differentiable and the function f is continuously differentiable. Then, the components of q*(*t*) *are piecewise parabolic while it remains inside any bounded set D* ⊂ ^{n}.

*Proof.* Let *q _{i}*(

*t*) be an arbitrary component of

*q*(

*t*). Taking into account (6) and the fact that

*t*

_{j}_{+1}>

*t*, ensuring that

_{j}*q*(

_{i}*t*) is piecewise parabolic is equivalent to prove that .

We shall start assuming the opposite, i.e.,

(10) |

which implies that . Then, given > 0, a natural *N* exists so that

Continuously differentiability of *f* and left differentiability of *u*(*t*) with the conditions *q*(*t*) ∈ *D* and ||*u*(*t*)|| < *U* imply that positive constants and *U _{t}* exist so that

Then, from (8) it results that

(12)

and, from (9), we have

Using (11) it results that, when *j* > *N*, it is true that

Choosing we get

Thus, when *j* > *N*, the succession can only grow if . Then,

(13)

where

From (7) it follows that

Continuity of *x _{i}* ensures that

|*x _{i}*(

*t*) -

_{j}*x*(

_{i}*t*

_{j}_{+1})| <

*F*(

_{i}*t*

_{j}_{+1}-

*t*)

_{j}and then, using (12) and (13) we obtain

Δ*q _{i}* < 2

*F*(

_{i}*t*

_{j}_{+1}-

*t*) +

_{j}*P*(

_{i}*t*

_{j}_{+1}-

*t*)

_{j}^{2}

and

which contradicts (10) and completes the proof.

**Corollary 1.***When f is (piecewise) linear w.r.t. x and u, and u*(*t*) *is piecewise parabolic, the state derivatives _{i}*(

*t*)

*are piecewise parabolic and the state variables x*(

_{i}*t*)

*follow piecewise cubic trajectories.*

In the nonlinear cases, we cannot say anything about the form of the state trajectories. However, we can approximate function *f* by a piecewise linear one. Similarly, when *u*(*t*) is not piecewise parabolic, we can use a piecewise parabolic approximation. Thus, we shall consider that the state trajectories and their derivatives are piecewise cubic and parabolic respectively.

Thus, as we did in QSS and QSS2, we can divide the system (2) into quantized integrators and static functions as shown in Fig. 1 so that each subsystem has piecewise parabolic input and output trajectories.

**C. Third order quantized integrator**

Quantized integrators in QSS3 will have piecewise parabolic input and output trajectories. They calculate *q _{i}*(

*t*) from

*d*(

_{i}*t*)

_{i}(

*t*). We shall deduce here the relationship between these trajectories in order to build an equivalent DEVS model. In order to simplify the notation, we shall eliminate the sub-index

*i*.

Let *d*(*t*) be a known piecewise parabolic trajectory

where τ_{k} ≤ t < τ_{k+1}. Notice that *d*(*t*) is defined by the sequences τ_{k}, *d*(τ_{k}), and .

Let *x*(*t*) be its integral, then

and let *q*(*t*) be a second order quantized version of *x*(*t*). Taking into account that it is piecewise parabolic, *q*(*t*) can be written as

with *q*(*t*_{0}) = *x*(*t*_{0}) (the initial condition) and = = 0. After that, whenever *q*(*t*) differs from *x*(*t*) in Δ*q*, *q*(*t*) starts a new segment. In those instants the quantized variable can be computed as

and their first and second derivatives are

This allows to calculate the output sequences *q*(*t _{j}*), and . The time succession

*t*can be calculated by

_{j} *t _{j}*

_{+1}= min[

*t*/|

*q*(

*t*

^{-}) -

*x*(

*t*)| = Δ

*q*] (14)

_{i} with *t* > *t _{j}*. Notice that the calculation of

*t*

_{j}_{+1}requires solving a cubic equation to determine when the difference between

*q*(

*t*) and

*x*(

*t*) becomes equal to the quantum.

Then, given a piecewise parabolic trajectory *d*(*t*) expressed by a sequence of events carrying the values *d*(τ_{k}), and , we can calculate the trajectory *q*(*t*), expressed as another sequence of events with values *q*(*t _{j}*), and .

This behavior can be easily represented by an atomic DEVS model that we shall call *third order quantized integrator* (Kofman, 2005a). The model was also implemented in PowerDEVS (the integrator block permits choosing the method, i.e., QSS, QSS2 or QSS3).

**D. Static functions**

The right hand side of (2) is composed by *n* functions *f _{i}*(

*q*,

*u*). Since

*q*and

*u*follow piecewise parabolic trajectories, in the linear case, the output

_{i}(

*t*) will be also piecewise linear.

In order to obtain a DEVS model relating the input and output trajectories of such a function we shall consider a linear function

where

so that

This expression defines *d*(*t*), and the sequences *d*(τ_{k}), and which allow expressing the trajectory as a sequence of events.

A DEVS model of this behavior is almost straightforward (Kofman, 2005a). It was also implemented in PowerDEVS in a block called *WSum* (which stands for weighted sum), that is common to the three methods.

The nonlinear case is a bit more complicated. Here, *d*(*t*) is no longer piecewise parabolic. However, we can approximate it by a piecewise parabolic trajectory, discarding the higher order terms.

We have

*d*(*t*) = *f _{i}*(

*v*1(

*t*), ...,

*v*(

_{N}*t*)) =

*f*(

_{i}*v*(

*t*))

where

Using Taylor's formulae, we have

and then

Discarding the higher order terms (starting from *t*^{3}), we have an expression for *d*(τ_{k}), and .

Notice that we need to evaluate and at *v*(τ_{k}). When we have the expression of *f _{i}* in closed form, it can be easily done. In more general cases, we might need to evaluate it numerically.

Some particular cases of nonlinear static functions are already programmed in PowerDEVS for the QSS3 method. For instance, there are blocks that calculate the product, trigonometric and power functions.

**E. Input signals**

The incorporation of input signals in QSS3 does not differ from QSS and QSS2. The only difference is that now we are allowed to consider piecewise parabolic approximations, which can reduce the number sections (and the number of events) with respect to the piece-wise constant and linear signals of QSS and QSS2.

Thus, the corresponding *signal sources* of QSS3 will be just DEVS generators which provoke events with the successive values of *u _{i}*(

*t*) (and their first and second derivatives).

A wide variety of input sources for QSS3 were included in PowerDEVS, and some of them are used in the examples of Section IV.

**F. Discontinuity handling**

Discontinuities can be managed in a similar way to QSS2. There, the discontinuity conditions were predicted by looking at the piecewise parabolic evolution

of the state *x*(*t*) (Kofman, 2004). Here, taking into account that a smaller quantization will be used, it is convenient to observe directly the quantized variable *q*(*t*) in order to avoid solving a cubic equation.

Some examples including blocks with discontinuity detection capabilities will be discussed in Section IV.

**G. Theoretical properties of QSS3**

From the definition of the second order quantizer -see (7)- it results that

|*q _{i}*(

*t*) -

*x*(

_{i}*t*)| ≤ Δ

*q*(15)

_{i}and then, all the properties mentioned in Sec.II.C are also satisfied by QSS3.

As in the lower order methods, the strongest property is that QSS3 has a calculateable global error bound in the simulation of LTI systems. Although this bound is the same for the three methods, in QSS3 we can use smaller quanta.

The time between successive steps on each quantized integrator is determined by the distance *t _{j}*

_{+1}-

*t*in (14). Immediately after a step,

_{j}*x*(

*t*) and

*q*(

*t*) have the same value and the same first and second derivatives. The only difference is in the third derivative, which is 2 ·

*p*in

_{d}*x*and 0 in

*q*. Provided that

*d*(

*t*) does not change between

*t*and

_{j}*t*

_{j}_{+1}, the time in which the difference between

*x*and

*q*becomes equal to Δ

*q*can be calculated as

Thus, we can conclude that the step size in a single integrator is proportional to the cubic root of the quantum. In QSS the step size was proportional to the quantum and in QSS2 it was proportional to the square root.

For example, if we want to increase the accuracy by a factor of 1 × 10^{6} we will have to reduce the quantum by the same factor (due to the global error bound property). In QSS it is equivalent to multiply by 1000000 the number of calculations. In QSS2 we would have to multiply it by 1000 and in QSS3 only by 100.

This fact not only permits using smaller quantization achieving better accuracy, but also makes less critical the election of the quantum since it can be chosen conservatively small without affecting considerably the computational costs.

**IV. EXAMPLES**

**A. DC motor with PWM control**

We consider in this example a DC motor with constant field, described by the equations

where *i _{a}*(

*t*) and ω(

*t*) are the armature current and the angular speed of the motor.

The inputs *U _{a}*(

*t*) and τ(

*t*) are the armature voltage and load torque. The parameters

*L*,

_{a}*R*,

_{a}*k*,

_{m}*J*and

*b*are the armature inductance and resistance, the motor constant, the inertia, and the friction coefficient respectively.

_{m} A typical strategy to control the motor speed is called pulse width modulation. The motor speed is compared with a desired reference and then the armature voltage switches from a positive value (+*V*) to a negative value (-*V*) so that the duration of the resulting pulses is proportional to the error.

A way of achieving this is comparing the error with a triangular waveform (carrier wave) applying +*V* or -*V* according to the sign of the difference.

Using the parameters corresponding to a real DC motor: *L _{a}* = 0.003,

*R*= 0.05,

_{a}*k*= 6.783,

_{m}*J*= 15, and

_{m}*b*= 0.005 we simulated the response of the system to a speed reference which goes from 0 to 60 with a rising time of 2 seconds. The DC motor is initially unloaded (τ(0) = 0) and a step of 2500 is applied in

_{m}*t*= 3. The triangular waveform was set with a frequency of 1000

*Hz*and an amplitude of 1.1

The PowerDEVS model is shown in Fig. 4. There, the blocks *QSS3 Integrator and QSS3 Linear* correspond to the third order quantized integrator and the linear static function described in the previous section. The *Triangular* and *Step* blocks are simple DEVS generators and the *Saturation* block bounds the output between 1 and -1 (so that the error does not becomes greater than the triangular wave, limiting in that way the maximum duty cycle).

The *SwitchTraj* block compares the error and the triangular wave, provoking events with values *V* or -*V* when they become equal to each other. This block predicts the trajectory crossings using the fact that they are piece wise parabolic. Thus, the discontinuities are exactly detected and handled.

Figure 4. PowerDEVS model of a PWM control.

For the simulation, a quantum of 0.001 was used in both state variables in order to appreciate the oscillations of the speed and current.

The QSS3 method completed the simulation of the first 5 seconds of the evolution after 7029 and 27572 steps in the integrators which calculate ω and *i _{a}* respectively. Additionally, the block which produces the triangular waveform provoked 10000 events (i.e. 2 events by period during 5 seconds of simulation) and consequently, the switch produced 10000 commutations. In that way, there were a total of 54600 events. 34600 events corresponded to the continuous part and 20000 to the discrete part.

The simulation in Power DEVS took 1.37 seconds on 450MHz PC under Windows 98. The results are shown in Figures 5-7.

Figure 5. Motor speed.

Figure 6. Motor speed (transient after the torque step).

The experiment was repeated using the QSS2 method with the same quanta. Now, the number of steps in the integrators was 17644 and 124021 in ω and *i _{a}* respectively. The addition of the 10000 events in the discrete subsystem gives a total of 161665 events. Power DEVS needed 3.18 seconds to simulate the system on the same computer than before.

We tried to simulate the same system with Simulink. The best results were obtained with ode23s. In order to get a qualitatively good result we needed to set

Figure 7. Armature current (final oscillations).

the relative tolerance to 2 · 10^{-10} and the absolute tolerance to 10^{-7}. The number of steps was 186593 and the simulation took 21.75 seconds.

The number of steps performed by QSS3 was less than the third part of the steps of ode23s. Besides this, each step in QSS3 only involves a few calculations. The 10000 events given by the triangular wave generator are only seen by the switch model which decides when to apply +*V* or -*V*.

Similarly, the 10000 events given by the switch are only seen by the integrator that computes *i _{a}*(

*t*). The 27572 steps of the integrator of

*i*(

_{a}*t*) are only seen by itself and the other integrator. The only expensive events are the 7029 given by the integrator of ω(

*t*) which are seen by both integrators and the discrete subsystem.

This clearly explains the fact that the simulation with Power DEVS is much faster -about 15 times- than the simulation with Matlab.

**B. A ball bouncing downstairs**

Consider a ball moving in two dimensions (*x* and *y*) bouncing downstairs. It will be assumed that the ball has a model in the air -with the presence of friction-and a different model in the floor (spring-damper).

A possible model is given by the set of equations

where *s _{w}* is equal to 1 in the floor and 0 in the air. Function int(

*h*+ 1 -

*x*) gives the height of the floor at a given position (h is the height of the first step and steps of 1

*m*by 1

*m*are considered).

The commutations are produced when the equation *y* = int(*h* + 1 - *x*) is verified.

The system was simulated in PowerDEVS with the QSS3 method. The block diagram was built in a similar way to the previous example.

In this case, a quantum of 0.00001 was chosen for the vertical position and 0.001 in the other variables. The initial conditions where *x*(0) = 0.575, *y*(0) = 10.5, *v _{x}*(0) = 0.5 and

*v*(0) = 0 and the first 10 seconds of the system evolution were simulated.

_{y} The QSS3 method performed 464, 346, 10 and 6 steps in the integrators which calculate *y*, *v _{y}*,

*x*and

*v*respectively. Thus, there were a total of 826 steps. PowerDEVS took about 0.033 seconds to complete the simulation on a 450 MHz PC running under Windows 98. The results are shown in Fig. 8.

_{x}

Figure 8. Ball bouncing downstairs (*y* vs. *x*).

The simulation with QSS2 and the same parameters takes a total of 9346 steps which are performed in 0.11 seconds on the same computer.

The simulation with Matlab variable step methods requires using a relative and absolute tolerance under 7 × 10^{-7}. Otherwise, the methods skip events. The best results are obtained with the ode23s method which performs 2960 steps and takes 0.27 seconds (also on the same computer). Again, a noticeable reduction of simulation time is observed.

**V. CONCLUSIONS**

We introduced a new numerical method called QSS3 which performs a third order discrete event approximation of ordinary differential equations.

We proved that QSS3 cannot produce an infinite sequence of events in a finite interval of time and hence the simulation with QSS3 will never be stuck.

The input-output relationship of quantized integrators and static functions was deduced so that DEVS models of these elementary blocks can be easily built.

Then, we saw that QSS3 has the same theoretical and practical properties than its lower order predecessors (QSS and QSS2). However, we showed that the reduction of the global error bound provokes a smaller increment of the computational costs.

The examples analyzed showed that the new method can achieve a noticeable reduction of the computational costs with respect to sophisticated discrete time methods in the simulation of discontinuous systems.

Future work should generalize the usage of the method to general nonlinear systems, considering also the DAE case.

^{1} In nonlinear systems this is only approximated.

^{2} Symbol | · | denotes the component-wise module of a complex matrix or vector and symbol "≤" in (5) also denotes a component-wise inequality.

**REFERENCES**

1. Hairer, E., S. Norsett, and G. Wanner, *Solving Ordinary Differential Equations I. Nonstiff Problems*, Springer, 2nd edition (1993). [ Links ]

2. Hairer, E. and G. Wanner, *Solving Ordinary Differential Equations II. Stiff and Differential-Algebraic Problems*, Springer, 1st edition (1991). [ Links ]

3. Kofman, E., "A Second Order Approximation for DEVS Simulation of Continuous Systems", *Simulation*, **78**, 76-89 (2002). [ Links ]

4. Kofman, E., "Quantization Based Simulation of Differential Algebraic Equation Systems", *Simulation*, **79**, 363-376 (2003). [ Links ]

5. Kofman, E., "Discrete Event Simulation of Hybrid Systems", *SIAM Journal on Scientific Computing*, **25**, 1771-1797 (2004). [ Links ]

6. Kofman, E., "A Third Order Discrete Event Simulation Method for Continuous System Simulation. Part II: Applications", *Proceedings of RPIC'05* (2005a). [ Links ]

7. Kofman, E., "Non conservative ultimate bound estimation in lti perturbed systems", *Automatica*, **41**, 1835-1838 (2005b). [ 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. Kofman, E., J. Lee, and B. Zeigler, "DEVS Representation of Differential Equation Systems. Review of Recent Advances", *Proceedings of ESS'01*, 591-595 (2001). [ Links ]

10. Pagliero, E. and M. Lapadula, "Herramienta Integrada de Modelado y Simulación de Sistemas de Eventos Discretos", Diploma Work. FCEIA, UNR, Argentina (2002). [ Links ]

11. Shampine, L. and M. Reichelt, "The MATLAB ODE Suite", *SIAM Journal on Scientific Computing*, **18**, 1-22 (1997). [ Links ]

12. Zeigler, B., T. Kim, and H. Praehofer, *Theory of Modeling and Simulation. Second edition*, Academic Press, New York (2000). [ Links ]

13. Zeigler, B. and J. Lee, "Theory of quantized systems: formal basis for DEVS/HLA distributed simulation environment", *SPIE Proceedings*, 49-58 (1998). [ Links ]

**Received: September 21, 2005. Accepted for publication: February 10, 2006. Recommended by Guest Editors C. De Angelo, J. Figueroa, G. García and J. Solsona. **