Stochastic Timed Discrete-Event Systems: Modular Modeling and Performance Evaluation Through Markovian Jumps

We are interested in a scalable, flexible, and modular methodology, for modeling and performance analysis of stochastic discrete-event systems (SDES). In this sense, we propose a modular approach for timing non-markovian SDES expressed as a parallel composition of modules that interacts with each other through events. We show how general distribution for event lifetimes can be implemented systematically by coupling timing modules to the system model. As a result, this coupling mechanism preserves modularity, leading to a compact markovian model expressed in terms of flexible modules. Therefore the methodology allows us to write the whole SDES model as a composition of the system model and the timing one, giving flexibility and scalability in modeling design, as we can modify the modules individually according to the designer’s interests. In addition, from the whole markovian SDES model, we show how to perform the model analysis through the analytic approach, as well as through Monte Carlo computer simulation. As an application, we present a numerical example of computing the abandonment rate for a service network with general service time employing both analytical and computer-simulation models.


I. INTRODUCTION
Deterministic Timed Discrete-Event Systems (DTDES) applies well in situations where the uncertainties associated with the description of the variables involved are negligible, in other words, from a probabilistic point of view, the variance is very small compared to the expected value. Deterministic models can be used, for example, in the modeling and control of automated manufacturing systems [10], [14], [26], [27], transportation networks [17], diagnosing [5], learning [2], [9] and real-time scheduling of tasks on processors [6], [7].
In the present paper, we are concerned with Stochastic Discrete-Event Systems [4], which are used for modeling, controlling, diagnosing, and analyzing the performance of uncertain systems whose dynamics are driven by the occurrence of random events [1], [16], [24]. Applications include The associate editor coordinating the review of this manuscript and approving it for publication was Tao Wang . industry and intelligent systems, networked systems, service networks, material handling systems [15], [22], [23] and maintenance [21]. In particular, we are interested in taking into account explicitly time in the model construction, leading to what we denote as a Stochastic Timed Discrete-Event Systems (STDES).
An important point to be highlighted is that the models in STDES are in, general, complex due to the discrete nature of the entities, with the number of states that grow in a combinatorial way as more information about the studied process is added to the models. In this sense, the construction of models in a modular way simplifies and systematizes the task of modeling complex systems since the whole model is obtained by composing simpler sub-models.
The utilization of Modular Models is a kind of ''divideand-conquer'' strategy used to solve problems in Discrete-Event Systems in different ways. In this context, we can cite, for instance: developments in analysis of untimed models [11], [18], [28], [12], and [13]; approaches for reducing the computational complexity of synthesizing controllers in the supervisory control theory for timed models [25] and for simulation purposes of complex systems [8].
Unlike previous papers found so far in the literature, we are interested in the study of a modular approach for timing non-markovian STDES expressed as a parallel composition of modules that interacts with each other through events. We show how the general distribution for event lifetimes can be implemented by coupling timing modules to the system model. As a result, this coupling leads to a compact markovian system expressed in terms of flexible modules. This modularity is interesting for modeling purposes since it divides the whole model into system sub-models and the timing sub-models for events that have a non-markovian lifetime, giving flexibility and scalability in modeling design, since the modules can be modified individually according to the designer's interests. In the paper, we show how to use the resulting modular markovian STDES methodology to perform the model analysis through analytic and Monte Carlo computer simulation approaches.
Paper Contributions: The model construction for a Stochastic Timed Discrete Event -System (STDES) is done through the parallel composition of modules based on Stochastic Timed Automaton (STA). STA is a very compact structure with few entities allowing solid mathematical developments. As a result, the approach simplifies and reduces the errors in the modeling task, since the whole model is constructed through the composition of simpler sub-models. Moreover, it allows us to incorporate timing mechanism as modules described in terms of STA's, as well. Therefore, we have a complete modular description for the STDES in terms of STA's. Besides the flexibility in changing the modules, or parameters, according to the designer's interests, the approach gives us flexibility in developing equivalent analytical or Monte Carlo computer simulation models, as we show in the paper.
The sequence of this paper is organized as follows: in Section II, we define the basic elements of STDES that we are interested in, as well as the operation of parallel composition of sub-models. We also discuss how to analyze the dynamics of a STDES using analytic and Monte Carlo computer simulation approaches. In section III, we show how to construct and couple markovian timing modules to markovianize a STDES model, in particular, we develop timing modules for hypoexponential and hyperexponential distributions. In Section IV we discuss the space complexity for the model representation, as well as time complexity for Analytic and Computer Simulation approaches. In Section V, we show a numerical-application example for a service network with customers' abandonment. The conclusion and perspectives for this work are presented in Section VI.

II. STOCHASTIC TIMED AUTOMATA MODELS
It is known that the stochastic automaton is a very general structure to represent STDES, and it can even be used for simulation via discrete computing. In this paper, we are interested in STDES that can be modeled through the stochastic timed automaton defined as follows: is a stochastic set of clock structure for generating lifetime of events, being V i a nonnegative random variable. In this work, the transition function is not total, that is, quite often some events are not allowed, or not feasible in a given state. This fact leads to the Definition 2.
Definition 2: For a STA G = (X , E, f , x 0 , V), we define by the function of feasible events: As a result, we can observe that : X → 2 E , being 2 E the power set, that is, the set of all sub-set of E.
To computer simulate the behavior of a STA, we need a numerical representation for the transition function f . To this end, the transition functions f of an automaton is represented by matrix according to Definition 3, in a general way, assuming that E is a subset of larger event sets denoted by E tot . Definition 3 (Matrix Representation): Given a STA G = (X , E, f , x 0 , V), an event set expressed as a sequence E tot = {1, 2, . . . , L} such that E ⊆ E tot and the state set given by X = {1, 2, . . . , N }, the transition function f is presented by the matrix A, N × L such that: In addition, we need to define the input-event set for a state, given a precedent state. This is achieved with the input function as defined in Definition 4.
Definition 4: Consider a STA G = (X , E, f , x 0 , V), the input-event set of a state x ∈ X , given a state y ∈ X , is defined as: To deal with huge models, we use a ''divide and conquer strategy'' that consists of a decomposition of the whole model into sub-models that interact with each other using common events. This can be done systematically by using the parallelcomposition operation as presented in Definition 5.
Definition 5 (Paralel Composition of STA): Parallel composition is a binary operation between two STA's G 1 = X 1 , E 1 , f 1 , x 1,0 , V 1 and G 2 = X 2 , E 2 , f 2 , x 2,0 , V 2 defined by: if e ∈ 2 (x 2 ) This operation is associative, that is: Example 1: We consider a simple manufacturing cell, with a waiting space F with capacity for two pieces and a machine S with capacity for processing one piece at a time. We consider the following events a, which represents the arrival of pieces, c, representing piece admission by the machine, and d representing that a piece was processed and has left the system. This system is then divided into two subsystems representing by the modules F and S, whose automaton model can be obtained easily focusing only on such subsystems, as show in Figure 2 As a result, the complete model is represented, in a compact way, by the parallel composition F||S. It's important to remark that all states of the system are ''embedded'' in the automaton resulting from F||S and we do not need to explicit the states. However, if desired, we can explicitly show all the states using the Definition 5. For the present example F||S results in the single automaton shown in Figure 3.

A. ANALYTIC APPROACH
Since exponential distribution is a practical way of generating non-negative random variables, as event lifetime, we can π 0 is the vector probability 2 of initial state x 0 ∈ X G . As a result, the dynamical evolution of a MCA if expressed by the following autonomous system: for which π (t) = [π 1 (t) . . . π n (t)] is the state probability vector whose π 1 (t) = P[X (t) = i], being X (t) a random variable that indicates the state of the system at time t. Quite often the system converge to a steady-state behavior, whose probabilities are constant, leading to: In addition, we need to add the probability constraint Since the rows of Q have linear dependence, we can eliminate, for instance, the first row, resulting in a system of the form:Q beingQ identical to Q, except for the first row, whose entries are 1, and b the null vector except for b(1) = 1.
There are several methods to solve Equation 7, among them the simplest one is based on inverse ofQ.
Example 2: Returning to Example 1, we consider the single automaton depicted in Figure 3. Given the temporal evolution of the system, the state space on the date t is given by: If the event lifetimes are random variables, the state of the system at date t is uncertain. Thus, we define the following probabilities: We suppose that the lifetimes are generated to respect exponential distributions as follows: where V e is the random variables that define the lifetimes for the events e. So in the present example e ∈ {a, c, d}.
In this case, we obtain an exact representation of the stochastic automaton as a Markov chain, as shown in Figure 4.

B. MONTE CARLO COMPUTER SIMULATION APPROACH
First of all, we recall that in this case, in general, it is assumed that the system enters into a steady state, more specifically, it is assumed that the associated stochastic process is stationary on average and ergodic. This being the case, the expected value can be estimated by the equation: where X i is the i-th sample of a random variable X . However, in view of the limitation of discrete computing resources, in terms of time and memory, a finite size sample must be considered. In this sense, in general, an estimatorˆ is built, that is, a finite approximation for Equation 10 is given by the the following estimator:ˆ Note thatˆ is a random variable that depends on the X i samples. For this estimator to be an acceptable approximation for θ = E[X ], we need an evaluation of the variance ofˆ .
Performance evaluation of a STDES through Monte Carlo Computer Simulation demand a lot of human and computational resources, since the models are, in general, quite complex. However the implementation of a computer simulation of a Markov STA is direct due to the memoryless property of the exponential distribution. Indeed, it can be implemented for a given STA G through a simple Monte Carlo simulation based on generation of uniform random numbers U ∈ [0 1] as presented in the Algorithm 1.

III. MARKOVIANIZATION WITH MARKOV JUMP MODULES
In Section II, we defined the basic elements of STDES that we are interested in, as well as the operation of parallel composition of sub-models. We also show how to analyze the dynamics of STDES using analytic or Monte Carlo computer simulation approaches. Based on those concepts, we show In this section how to construct Markov jump modules in order to transform a non-markovian STDES into a markovian one.
First, in order to have a general automaton representation for the timing mechanism, we construct an automaton T e as depicted in Figure 5. For this automaton all event lifetime are exponentially generated as: • the jump rate off state X into I being λ; • the jump rate from state I to (k, 1) being p k λ; • from state (k, j) given by λ jk ; • and from state E given by λ. So if we have an event that causes a ''non-markovian jump'', let's say from state ''X'' to ''Y'', the operation of parallel composition insert automatically several intermediary states whose transitions operate with markovian jumps in such a way that the '' big jump'' between ''X'' and ''Y'' respects the desired distribution. This operation must be done for each event e that does not respect an exponential distribution. In practical situations, minimal topology can be derived by matching the mean and variance of the collected data as we show in the following. VOLUME 10, 2022 FIGURE 5. Automaton T e that implements general distribution as a timing mechanism, in which states X and E represent respectively the origin and the destination of an event. As a result, if we denote the corresponding timing STA for each event e i by T i , the complete timing mechanism is therefore given by the parallel composition: In addition, if the starting of the jumping process is indicated with an event, let's say s e , we must include this event in the alphabet of the system by making it feasible whenever event e is. This can be done by inserting an appropriate trigger mechanism as shown in Figure 6.
Therefore given an STA, denoted as G, with general lifetime distribution for some events, lets say e ∈ E t = {1, . . . p}, we create and associate starting events s e ∈ E s for those events, resulting the set of pairs: Then we construct a modified STA by adding self-loops of these starting event e s in all states x of G for which event ''e'' is feasible, i.e. e ∈ (x). This procedure of insertion self-loops is illustrated in in Figure 6 and results in a new automaton, which we denote as T (G, P). This transformed STA is presented formally in Definition 7 Definition 7: Given G = (X , E, f , x 0 , V) and the set of pairs P, we define the Triggered STA as: for which: x if e = s e and e ∈ (x).
• V t = V ∪ V s ; being V s the lifetime structure for the starting events and (x) the set of feasible events for a given state x of the automaton G.
As a result the the complete STA, with markovian jumps, denoted as G Markov , is given by the parallel composition of triggered STA and the corresponding STA for time mechanism as: It is important to remark that the corresponding Markov chain is systematically obtained directly from G Markov by replacing events with their rates. An important advantage of the compact representation of the Markov Chain as given by the Equation 15 is the reduced space complexity, since the operation of parallel composition allows us to represent the whole system through modules. Another interesting feature of this modularity is the flexibility to change the distributions by replacing modules.

A. TIMING MODULES FOR HYPOEXPONENTIAL AND HYPEREXPONENTIAL DISTRIBUTIONS
In a practical situation, hypoexponential and hyperexponential are quite a general way of representing unimodal distribution for time random variables [3]. Important specifications to fit those distributions are the well-known mean (µ) and variance (σ 2 ) of the input data. In this sense, we will discuss in sequence how can we deal with a practical situation with data described as a hypoexponential distribution, that is σ < µ or a hyperexponential one, that is σ > µ. In particular, we will see how to obtain minimal structures ensuring mean and variance matching for those distributions FIGURE 7. Sequence of markov jumps for hypoexponential with N steps. We assume that λ i the jump rate from step i to step i + 1.

1) CONVENIENT MINIMAL STRUCTURES FOR MEAN AND VARIANCE MATCHING
The simplified diagram for the states of an Hypoexponential distribution is depicted as shown in Figure 7.
Therefore, our matching problem is achieved by solving the following pair of equation for a minimum N : whose m i = 1/λ i , being λ i the jump rate from state i to state i + 1. First, we must observe that triangular and internal product inequalities (Cauchy-Schwartz Inequality) lead to: As a result, provided that a solutions exists, they are all such that: So the number of states N must satisfy: Since N is an integer, its smallest possible value is 3 N = µ 2 σ 2 . Keeping in mind those observations, we present a solution for the system of equations 16. First we observe that if µ 2 σ 2 = 1, the solution is trivial, with only one state, that is N = 1. So let us concentrate our attention in situations whose µ 2 σ 2 ≥ 2.

Proposition 1 (Minimal Hypoexponential Distribution):
If µ 2 σ 2 ≥ 2, a solution for the system of equations 16 is given by: and being 3 x stands for the ceil of x, i.e. smallest integer greater than or equal to x. . STA that implements a hypoexponential distribution as an event timing mechanism for σ < µ, whose µ 2 σ 2 ≤ N < µ 2 σ 2 + 1.

FIGURE 9.
Sequence of markov jumps possibilities for hyper-exponential distributions. We denote by p the probability of initially routing to state ''1'', and by λ i the jump rate to get out from state i .
Proof: By denoting m j = y (1 ≤ j ≤ N −1) and x N = z, The system 16 is written as: being N = µ 2 σ 2 . As a result, we can check that: are indeed solutions for the pair of Equations 22. As a result, we obtain the minimum STA as shown in Figure 8.
To couple the time module into the system for a given event e, we do the following: • Create a Triggered STA according to Definition 7 using the procedure illustrated previously in Figure 6; • The time scheme is presented by the STA T e depicted in Figure 8, with jump rates given off state j ∈ {1, . . . , N } by λ j = 1/m j , being m j given by Equations 20 and 21. So far, we have established that if σ 2 ≤ µ 2 , the minimum number of states is given by N = µ 2 σ 2 , being the means between states provided by Proposition 1. It remains to solve the cases for which σ is strictly greater than µ, that is σ > µ. We solve this problem by considering two states and two routing probabilities, as depicted in Figure 9. This configuration leads us to an Hyperexponential distribution.
Denoting m i = 1 λ i , we can write that: VOLUME 10, 2022  FIGURE 10. Visualization of the structure that implements a hyper-exponential distribution. Considering, without loss of generality that m 1 ≥ m 2 , we solve the system of equation, whose solutions are given by: and C v the coefficient of variation, i.e C v = σ µ . So the simplest topology to ensure the desired results is given by choosing and m 2 = 0. This topology, which has only one state, is shown in Figure 10. Remark 1: An interpretation of the resulting distribution for this topology, in terms of service time, is the following: with probability p, some clients are served with exponentially distributed service time with mean m 1 , while others are served with negligible service mean with probability (1 − p). As a result, we obtain the corresponding minimal automaton implementation for this distribution as depicted respectively in Figure 11.
In order to couple the time mechanism into the system, given event e, we follow a similar procedure as we did for hypoexponential distributions: • Create a Triggered STA according to Definition 7 using the procedure illustrated previously in Figure 6; • The time scheme is presented by the STA T e depicted in Figure 8, with jump rates given off state j ∈ {1, . . . , N } by λ j = 1/m j , being m j given by Equations 20 and 21.

IV. COMPLEXITY ANALYSIS OF THE MODULAR MARKOV REPRESENTATION
First, we analyze the memory requirements complexity of the modular model representation, then we analyze the time complexity for performance evaluation using analytic approach as well as computer simulation one.

A. MODEL REPRESENTATION
The complexity Analysis is performed for a model given by a parallel composition of sub-models expressed by Equation 15. Explicitly: being G s = G 1 . . . G L the modular STA model for the system, P the set of pairs for the non-markovian events, as indicated in Equation 13, and T m = T 1 . . . T p the modular STA representation for the timing mechanism.
If we denote |H| the number of states of an arbitrary STA H, the space complexity for the implicit model expressed by the parallel composition is linear in terms of the number of the sub-models: On the other hand, for the explicit (monolithic) model, for which we explicitly represent all states, the space complexity is exponential: Therefore, comparing Equations 28 and 29, we can see that parallel representation saves much more memory than the explicit one.
Regarding time complexity for performance evaluation, in the following we compare the worst-case computational complexity for analytic as well as Monte Carlo computer simulation approaches.

1) ANALYTIC APPROACH
In this case the method is base on matrix inversion. So the worst-case complexity is given by: However, in practice, the matrices are very sparse for large instance problems, since the number of events is quite small in comparison with the number of states. Therefore the computational complexity can be much more smaller depending on the numerical method.
Numerical and Symbolic solution methods are presented in [3]. Numerical Methods for solving large sparse linear equation systems can be found in [19]. Symbolic representations and analysis of large system can be found in [20].

2) COMPUTER SIMULATION APPROACH
This method requires intensive computer iteration as we can observe in Algorithm 1. The number of iterations depends FIGURE 12. Service network with customers abandonment: Event ''a'' represents that a customer arrives to the queue F; ''c'' customer admission by the server S; ''d'' departure of customer from the server after receiving the service. The customer abandonment when the queue F has length i is represented by event ''q i ''.
on the desired precision. For instance, if we are interested in estimating a given expected value, we observe in Equation 10 and 11 that the large the number of iterations the better is the approximation. In practice, given a desired confidence interval, the worst-case complexity is written as: being N ev , N sp ×N rn respectively the number of events, computer iteration and the number of runs. We can observe that the complexity depends on sample number N sp , as well as on the number of runs N rn . These parameters depend on a previously established confidence interval. The sample number is in general very large (thousands and even millions of samples) and the number of runs is usually small (some units). Therefore, the approach requires intensive use of computer time, but, for large systems, it is simpler to implement than analytic approach.

V. NUMERICAL RESULTS: SERVICE NETWORK WITH CUSTOMERS ABANDONMENT
We evaluate the methodology with a study of a Service Network with Customers Abandonment whose service time is non-markovian but specified in terms of mean and variance. For the representation of this system, we consider two main subsystems: a queue, which is a space where the customers wait for the service, and a server that can process one customer at a time. In this system, customers conditions are sensitive to queue waiting time, in the sense that the longer the queue length the higher the abandonment rate. A block diagram of this system with the event descriptions is depicted in Figure 12.
The correspondent modular automaton model for the system is shown in Figure 13.
To numerically evaluate the performance of this system in terms of the overall abandonment rate of parts, we consider that all event lifetimes follow an exponential distribution, except for the service time, which is specified in terms of a hypoexponential distribution. We consider that the abandonment rates that are proportional to the queue length.
In the present example, we consider that the maximum queue length is M = 10 and the mean lifetime for  The corresponding STA's that implement S s = T (G, P) and T d are depicted in Figure 14 The complete STA model is the result of the parallel composition of the modules F, S d , and Td presented in the Figure 15.
In the sequence, we are interested in the performance evaluation of the system. We show how this can be achieved through analytic approach, as well, as through

3) ANALYTIC APPROACH
First, from Figure15, we check that the resulting STA model has 11 × 2 × 4 = 88 states, leading to an equivalent Markov Chain with the same number of states.
In order to evaluate the abandonment rate for the system, we need to compute the throughput rate, that is, the rate of the number of customers that receive a service from the server per time unit. Since we can say that a customer is served whenever an event ''d'' happens, we compute the number of customers served by the number of times the event ''d'' occurs. In this sense, we observe that the event of departure, ''d'' is feasible whenever the server is on state ''2'' and the Markov module is on state ''4'' no matter the queue state. As result, the throughput rate Tr is computed using the conditional probability formula: So the overall abandon rate, Tq is computed by the difference between arrival and throughput rates: whose λ a = 1 E[V a ] is the arrival rate. As a consequence, we obtain an abandonment percentage rate 100 × Tq λ a = 33, 50%. The computation time was 0.4516 s in a DELL Computer Intel I7, 3.10Ghz, 8 GB RAM, Windows 8 Pro with 64 bits.
For the sake of comparison, in the sequence, we do the same analysis through Computer Simulation Approach, using ScicosLab scicoslab.org, a variant of Scilab, to develop the programs.
Using the Algorithm1, with the same STA used previously in analytic approach, as depicted in Figure15, we do the computer simulation.
We performed the simulation until the conclusion of the service for 10000 pieces. To eliminate the transient behavior, the throughput was computed only after the 400th piece departure. We ran the simulation 10 times, to compute a confidence interval (IC) of 95% using t-Student distribution.
As a result, we obtain the IC for the abandonment percentage rate as [33.37% 33.60%], being the computation time 104.49 s in the same DELL Computer Intel I7, 3.10Ghz, 8 GB RAM, Windows 8.1 Pro with 64 bits.
Remark 2: Comparing the results obtained through the two approaches, that is, the abandonment percentage rate of 33, 50% for analytic and an IC of [33.37% 33.60%] for the computer simulation, we observe that they are quite coherent. However, in terms of computation time, we obtained 0.4516 s for analytic approach versus 104.49 s for computer simulation one, that is computer simulation is more than 200 times slower.

A. DISCUSSION
The model used in the numerical example had a moderate complexity, but we think that it was enough to show the potential of the methodology since we could observe the flexibility in the model construction by changing the whole model by simply replacing modules. For instance, changing the queue capacity, or timing mechanism, is a matter of replacing the respective module for the system. Moreover, despite its apparent simplicity, the total number of states for the example is 88, that is, it is difficult to visualize (or to even draw) the whole model resulting from the compositions of the modules. On the other hand, modular models for the queue, server and timing mechanism had respectively only 11, 2 and 4 states. We remark that the analytic approach gives us real values instead of interval estimates. For systems, with a small number of modules, it is faster, since the number of samples required in the computer simulation approach to generate a small confidence interval is usually very large (thousands or even millions of samples). However, computer simulation algorithms are simpler to implement.

VI. CONCLUSION
In this paper, we presented a complete methodology to analyze the behavior of stochastic discrete-event systems represented using STA, being the whole system expressed as a parallel composition of sub-systems or modules. We have shown how to incorporate timing modules to describe general non-exponential distributions for event lifetimes. In particular, we have shown how to design timing modules for hypoexponential and hyperexponential distributions. Therefore, we obtained a complete description of the system dynamics by the parallel composition of individual modules that interacts with each other through events, being dynamics operating through markovian jumps among states. With the aid of a numerical example we show how the task of obtaining the whole model of a system could be drastically simplified through the utilization of modular models for the subsystems, as well as for the timing mechanism. We can observe the versatility of the representation as we could analyze the dynamic behavior of the system by analytic or by Monte Carlo computer simulation methods. Finally, we remark that other application perspectives of the methodology, in terms of performance evaluation of STDES, include computing probabilities of accessing states, blocking, failures, or executing a given sequence of events, which is useful in the analysis and decision of system issues in areas such as control, diagnosing, maintenance and security. Another interesting works to be exploited in future works are an automated tool to convert from STDES to STA or formal system verification.