Latin American applied research - A third order discrete event method for continuous system simulation

## Articulo

• Similares en SciELO

## 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. Kofman1

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 = {Q0, Q1, ..., Qr} be a set of real numbers where Qk-1 < Qk with 1 ≤ kr. Let Ω be the set of piecewise continuous real valued trajectories and let xi ∈ Ω be a continuous trajectory. Let b : Ω → Ω be a mapping and let qi = b(xi) where the trajectory qi satisfies:

 (3)

and

Then, the map b is a hysteretic quantization function.

The discrete values Qk are called quantization levels and the distance Qk+1 - Qk 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 et al., 2001).

Kofman and Junco (2001) proved that the quantized variable trajectories qi(t) and the state derivatives i(t) are piecewise constant and the state variables xi(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 fi, 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 parabolic1.

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-1q (5)

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

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 xi(t) and qi(t) are related by a second order quantization function if qi(t0) = xi(t0) and

 (6)

with tttj+1, and the sequence t0, ..., tj, ... defined so that tj+1 is the minimum t > tj where

(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 tj 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 Dn.

Proof. Let qi(t) be an arbitrary component of q(t). Taking into account (6) and the fact that tj+1 > tj, ensuring that qi(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 Ut 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 xi ensures that

|xi(tj) - xi(tj+1)| < Fi(tj+1 - tj)

and then, using (12) and (13) we obtain

Δqi < 2Fi(tj+1 - tj) + Pi(tj+1 - tj)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 xi(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 qi(t) from di(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, dk), 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(t0) = x(t0) (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(tj), and . The time succession tj can be calculated by

tj+1 = min[t/|q(t-) - x(t)| = Δqi] (14)

with t > tj. Notice that the calculation of tj+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 dk), and , we can calculate the trajectory q(t), expressed as another sequence of events with values q(tj), 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 fi(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 dk), 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) = fi(v1(t), ..., vN(t)) = fi(v(t))

where

Using Taylor's formulae, we have

and then

Discarding the higher order terms (starting from t3), we have an expression for dk), and .

Notice that we need to evaluate and at vk). When we have the expression of fi 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 ui(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

|qi(t) - xi(t)| ≤ Δqi (15)

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 tj+1 - tj in (14). Immediately after a step, 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 · pd in x and 0 in q. Provided that d(t) does not change between tj and tj+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 × 106 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 ia(t) and ω(t) are the armature current and the angular speed of the motor.

The inputs Ua(t) and τ(t) are the armature voltage and load torque. The parameters La, Ra, km, J and bm are the armature inductance and resistance, the motor constant, the inertia, and the friction coefficient respectively.

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: La = 0.003, Ra = 0.05, km = 6.783, Jm = 15, and bm = 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 t = 3. The triangular waveform was set with a frequency of 1000Hz 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 ia 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 ia 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 ia(t). The 27572 steps of the integrator of ia(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 sw 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 1m by 1m 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, vx(0) = 0.5 and vy(0) = 0 and the first 10 seconds of the system evolution were simulated.

The QSS3 method performed 464, 346, 10 and 6 steps in the integrators which calculate y, vy, x and vx 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.

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).