Efficient Handling of Phase-type Distributions in Generalized Stochastic Petri Nets

We study the introduction of transitions with Phase-type distribution firing time in (bounded) generalized stochastic Petri nets. Such transitions produce large increases of both space and time complexity for the computation of the steady state probabilities of the underlying Markov chain. We propose a new approach to limit this phenomenon while keeping full stochastic semantics of previous works. The method is based on a structural decomposition of the net. We establish conditions under which this decomposition leads to a tensor expression of the generator of the chain. The tensor expression is used to solve the chain with an iterative method.


Introduction
Generalized Stochastic Petri Nets (GSPNs) integrate two kinds of transitions, namely immediate and exponential. Although this allows to model a large range of systems, the availability of more general distributions of ring time would be very useful for many real life phenomena. Two general directions have been taken in this area: the rst one 5,9] introduces general distributions of ring time with speci c conditions, which allows one to extract so called subordinated Markov chains and compute steady state and transient probabilities of the states. The second approach 16,12,7] introduces Phase-type (PH) distribution ring time without further restrictions, and compute the resulting embedded Continous Time Markov Chain (CTMC). This needs to precisely de ne the stochastic semantics of these new transitions and to cope with the "bi-dimensional" complexity induced by the number of states generated: exponential with respect to the size of the net, exponential under speci c hypotheses with respect to the number of stages of the Coxian distributions.
In the framework of this last approach, previous works have taken into account the stochastic semantics problem either at the net level 7], or at the underlying CTMC level directly 12]. Unfortunately, the complexity problem has not been tackled.
Our approach deals not only with the stochastic semantics allowing the most general situations but also with the complexity of the Markov chain. We exploit the concurrency structure The outline of the paper is as follows: in Sect. 1 we investigate the problems arising with the introduction of PH and Coxian distribution ring time of transitions (PH and Coxian transitions) in GSPNs. We then formally de ne the syntax and semantics of bounded GSPNs with such transitions in Sect. 2. In Sect. 3 we explain how to decompose the net to get a tensor expression of the generator of its underlying CTMC. We study the impact of our method on the computational complexity in Sect. 4 and we introduce immediate transitions in Sect. 5 ending with concluding remarks. 1 Dealing with PH transitions in stochastic Petri nets 1.1

PH and Coxian distributions
Coxian distributions 11] generalize Erlang and hyper-exponential distributions which are serial and parallel network of exponential servers. Their model (Fig. 1) is a serial network of s exponential servers (or stages) with service rate i . At the end of service i, clients may leave the global service with probability i or enter next stage (if i < s) with probability i = 1? i . It is also possible to leave immediately the service with probability 0 . PH distributions 17] are another family of distributions which may be used to approximate many distributions. A PH distribution is the time until absorption of a CTMC with s + 1 states and only one absorbing state. It may be also seen as an exponential servers network (Fig. 2), the initial distribution of the chain being (c O;j ) 1 j s and the absorbing state being the state 0. The Laplace transform of a PH distribution is a rational function and it is known that any probability distribution with rational Laplace transform may be (non uniquely) interpreted as a PH distribution, so that many real life phenomena may be modelled with such distributions and their introduction in GSPNs enlarges the usability of this model. In the rest of this paper, we shall only deal with Coxian distributions which are a special class of PH distributions, but the proposed method can be extended straightforward to the general case.

Stochastic semantics of SPNs
The stochastic semantics of SPNs, also called the execution policy 1], rst de nes which is the transition to be red in a given marking: in this paper we follow the race policy 2, 1]: the transition to re among the enabled transitions is the one with the smallest remaining ring time. The two other main components of the execution policy are the memory policy and

Memory policy
As Coxian distributions are no more memoryless, it is necessary to specify how the remaining ring time of a non selected enabled transition is computed. There are three main possibilities.
Resampling: elapsed time is lost and a new sampling of ring delay will be done the next time the transition will be enabled. Resampling allows to model con icting activities one of which fully disables the others like failures.
Enabling Memory: elapsed time is kept if the transition stays enabled. Time-out transition for example may have this memory policy: as long as it is enabled, its timer is decreased. If watch-dog conditions are no more true, the time-out is stopped.
Age Memory: elapsed time is kept whatever the evolution of the system, that is to say: the next time the transition will be enabled, the remaining ring time will be the residual time from the last disabling and not a new sampling. Time-sharing service may be modelled by age memory policies. From a practical point of view, resampling policy has only meaning for con icting transitions 2], and in this case, it may be modelled with enabling policy and immediate transitions, so that we shall not deal further with resampling in the paper.

Service semantics
The service semantics de nes what happens about the service, when the Enabling Degree (ED) of the transition is greater than 1, i.e., when several clients may be served. In standard GSPNs, the ring rate of the transition is used to translate the service policy. If we denote by (t) the rate of t and (t; m) its ring rate in the marking m, there are three usual situations: Multiple (K) Servers: at most K clients may be served simultaneously. The ring rate is now (t; m) = min(K; ED(t; m)): (t) (the Single Server policy may be seen as the special case K = 1 of the Multiple Servers policy).
In nite Servers: all clients able to be served are served. We have here (t; m) = ED(t; m): (t).
Note the very important fact that, for Coxian transitions with service policy not Single Server, we shall have to retain in which stages are the clients being served and to know what stage ends its service to get the ring rate of the transition in a given marking.

Interruption/resume policy
If several clients are served in the stages of a Coxian transition, we have to de ne how to select the clients to be interrupted if necessary. Two choices, at least, correspond to standard solutions: the most recently entered client (Last engaged policy, denoted by F for its connection with the FIFO policy) or the client nearest from the rst stage (Least engaged policy, denoted by L) is interrupted. Furthermore, if the memory policy is Age Memory, the resume policy which tells which interrupted client is resumed has to be de ned. It is clear that it must be chosen accordingly with the interrupt policy (e.g., resume the oldest client or the most engaged one) so that we name interrupt policy both interruption and resume (if needed) strategies. We summarize the di erent execution policies in table 1: each line gives the three parameters needed (sp(t); mp(t) and ip(t)) for every Coxian transition. 1.2.4 Complexity of the stochastic model The size of the Markovian model of GSPNs is roughly the same as the one of the reachability graph of the net, that is to say exponential with respect to cardinalities of the sets of places and transitions. As soon as we introduce Coxian transitions with the di erent semantics above, the size of the new CTMC will now be the product of the previous model size by the cardinalities of each of the state space sizes of the Coxian transitions. With Multiple or In nite Servers policies (and far more with Age Memory policy and FIFO strategy) this increases hugely the size of the model which leads to nd e cient method to solve it.

Previous works
Since the introduction of stochastic Petri nets, several works have tried to overcome the exponential restriction allowing more general distributions. Two main lines were proposed. 1. General distributions are allowed 9]. If two transitions with general distribution are enabled in a marking and if the ring of one of them does not reset the service time of the other, (like with Age Memory policy or Deterministic ring time for instance) the stochastic process is no more Markovian. Under speci c hypotheses, it is however possible to de ne regeneration points of the process which is then semi-regenerative 10].
The steady state resolution is based on the extraction of a so called embedded Markov Chain and the study of the process between two regeneration points 5]. 2. Restrictions about new distributions are such that the underlying stochastic process remains a CTMC. In the second case, which is the framework of the present paper, two methods have been studied.
The rst one alters the initial GSPN with Coxian transitions 16,7]: each Coxian transition is replaced in the net with a sub-net which translates its stochastic semantics. In 7] the authors propose replacement sub-nets of such transitions with a xed number of places and transitions whatever the number of stages of the Coxian distribution in contrast with 16]: for instance, a Coxian transition with Age Memory and Single Server policies is expanded with a 8 places, 7 transitions sub-net. This expansion may be automated and integrated in a tool. Note however that only the Single Server policy is studied and that moreover, each Coxian transition must have only one input and one output place.
The second method expandes the embedded Markov chain 12] (the author deals with PH distributions but the principle is the same for Coxian distributions): the net is not modi ed; the user simply gives the properties of each Coxian transition. These distributions are taken into account only during the reachability graph computation: on one hand, each "marking" holds standard marking information and all information on the states of each Coxian transition (clients in service, their stages); on the other hand, during the computation of rings from a given "marking", these new information are updated accordingly to the stochastic semantics of the Coxian transitions. This method allows the management of any combination of the three elements of the execution policy and may be also integrated in a tool, but at the solver level whereas the rst method must be included at a step preceding the solver in the analysis procedure.
The main goal of the present work is to deal with the intrinsic increase in complexity of the Markov chain induced by the introduction of PH or Coxian transitions, a problem which has not been studied in the contributions mentioned above.

De nition of generalized stochastic Petri nets with Coxian transitions
In this section we give the formal de nition of a generalized stochastic Petri net with Coxian transitions and its stochastic interpretation, that is to say its underlying Markov chain.

Coxian GSPNs
De nition 2.1 A Coxian generalized stochastic Petri net is a tuple N = (N 0 ; ; C ) where: N 0 = (P; T; Pre; Post; Inh; pri) is a Petri net: P is the set of places; T is the set of transitions; Pre, Post and Inh: P T ! Bag(P) are the incidence functions 1 : Pre is the input function, Post the output function and Inh the inhibition function; pri: T ! IN is the priority function. and 2 T = T E U T I U T C where T E , T I and T C are the sets of exponential, immediate and Coxian transitions. and C give the stochastic properties of the transitions: The choice of an initial marking m 0 , de nes a marked Coxian generalized stochastic Petri net called also a Coxian GSPN system S = (N; m 0 ).
As we are dealing here with the introduction of Coxian distributions, we restrict ourselves in this section to Coxian GSPNs without immediate transitions and to Coxian transitions without immediate ring (i.e., 0 = 0). The general case will be studied in Sect. 5 where the exact meaning of pri and will also be discussed for Coxian transitions.

Reachability graph of a Coxian GSPN
The reachability graph of a Coxian GSPN is the standard reachability graph of GSPNs where a Coxian transition is seen as an exponential one. During its building, the values of the Enabling Degree with Service constraint if sp(t) = I 1 Bag(P ) is the set of multi-sets on P: Bag(P ) = fx = (x(p)) p2P j 8p 2 P; x(p) 2 INg 2 A U B means A S B with A T B = ; are used to compute the greatest total number e(t) of clients which may e ectively be in all the stages of a Coxian transition in any marking of the reachability set RS of the net: e(t) = max m2RS fEDS(t; m)g We also assume that EDS(t; m) may always be accessed during the next steps of the qualitative as well as quantitative analysis of the net, without worrying about this problem in this paper.

Markov chain of a Coxian GSPN
A state of the chain is an enlarged marking which in addition to standard marking, gives for each Coxian transition, its internal state, that is to say in which exponential stages are the clients being served.
We call descriptor of the transition the information to memorize which depends on the properties of its stochastic semantics.

Descriptors for policies with identi ed clients
These are the policies with multiple servers (sp(t) = K or I) such that the interrupt policy requires to distinguish the clients, e.g., the Last engaged (FIFO) policy (ip(t) = F). We may have mp(t) = E or A.
The descriptor of these policies (C-descriptor) gives for each client, interrupted or not (if mp(t) = A), its service stage. As we know that EDS(t; m) = 0 i no client is being served,

Descriptors for un-identi ed clients policies
If there is no need to distinguish clients in the stages of the Coxian distribution, they are only known by the stage in which they are served, as for example with mp(t) = E. We then de ne the descriptor (S-descriptor) as d(t) = (d 1 ; : : :; d s(t) ) 2 SD(t) = f0; 1; : : : ; e(t)g s(t) If EDS(t; m) > 0, d i is the number of clients in service in the stage i. If  1 client in service in the stage 1 1 0 no client Each of these descriptors may be modi ed (simpli ed) when the execution policy allows it: for example, with enabling memory, it is useless to keep information about states of the interrupted clients. Moreover, each Coxian transition has of course, a descriptor adapted to its own stochastic semantics.
It is clear that the C-descriptors need the most memory since jCD(t)j = (s(t)) e(t) . The Sdescriptors use less memory (jSD(t)j = (e(t) + 1) s(t) ) and the SS-descriptors catch the minimal information (to be able to manage the Age Memory policy). We see that, for a xed number of stages of t, the memory used by CD(t) may grow exponentially with the cardinality of the reachability set of the net, whereas the growing is polynomial for the S-descriptors.

State transitions of the Markov chain
The transitions (with their rates) between states of the Markov chain may be deduced as for standard GSPNs, from a graph we call the Extended Reachability Graph 3. t is Coxian and its ring is the end of service of its rth stage with continuation on the following stage. We name -ring such a ring and we denote it by m t r ?! m 0 . The marking m is unchanged in this case.
The stochastic semantics of Coxian GSPNs is given by its RG. We now provide detailed algorithms to compute RG in the case where all Coxian transitions have C-descriptors. Note that such descriptors do not need to be modi ed when there are new engaged clients. Let m = (m; (d(t)) t2T C ) be an extended marking. The ring m t ?! of an exponential transition t leads to the unique extended marking m 0 = (m 0 ; (d 0 (t)) t2T C ) which is given by the algorithm 2. 1. The rate of this elementary state transition is (t; m).
The ring m t r ?! gives at the most EDS(t; m) markings m 0 (one for each possible client in stage r leaving the service) computed by the algorithm 2.2. The corresponding rate is r (t): r (t).
The ring m t r ?! provides at the most EDS(t; m) markings m 0 (possibly one for each possible client leaving the stage r) (algorithm 2.3) and each elementary state transition has the rate r (t): r (t). Let us emphasis that our goal is precisely to avoid the direct construction of RG (see below and Sect. 5).

Structural decomposition of Coxian GSPNs
The above algorithms show in details the impact of the introduction of Coxian transitions on memory requirements and computational cost of the steady state probabilities of the embedded Markov chain of a Coxian GSPN. The tensor decomposition that we introduce now, is the key point to cope with this problem. Let us recall that tensor decomposition 18  Since an extended marking consists of a standard marking and of a tuple of descriptors of the Coxian transitions, we have to introduce a partition of the set of places, which allows also to isolate the e ects of the transition rings on the descriptors. As we have seen in Sect. 2.4, the descriptor of a given transition t may be modi ed by a ring of t 0 i t 0 may change EDS(t; :). Let us remind that two transitions t and t 0 are in symmetric structural con ict (SSC) i the ring of t may disable t 0 or vice versa. The extended con ict sets (ECSs) of the net are the equivalence classes of the transitive closure of the SSC relation 3].
Hence, to ensure that a partition will allow a tensor decomposition, we impose to such a partition conditions based on the extended con ict sets (ECSs) of the net.
De nition 3.1 (Compatible partition) Let (ECS j ) 1 j n be the ECSs of a Coxian GSPN.
A partition (P i ) 1 i m of the set of places P is compatible with respect to Coxian transitions i 81 j n; (ECS j \ T C 6 = ; =) 9i j ; ECS j ECS j P i j ) where ECS j = fPre ?1 (:; t) j t 2 ECS j g and ECS j = fInh ?1 (:; t) j t 2 ECS j g. We denote by m i = m(P i ) the restriction to P i of a marking m, so that m = (m i ) 1 i m .
Let us consider for example the net of Fig. 3 (Coxian transitions have gray boxes). We have two ECSs with Coxian transitions (drawn with dashed line): ECS 1 = ft C1 ; t C2 g and ECS 2 = ft e5 ; t C3 g. So, we must choose a partition such that the input and inhibition places of t C1 and t C2 (resp. t e5 and t C3 ) belong to the same P i . Note that these conditions allow an automated computation of a compatible partition in polynomial time, namely the ECS j S ECS j , since the ECSs may be computed with polynomial time with respect to the size of P and T.
To distinguish local and global rings, we introduce the sets T i of transitions which are the only ones which may change the marking of a place of P i : T i = P i P i Furthermore, the set of Coxian transitions descriptors of which depend on P i is T C;i = ft 2 T C j ECS(t) ECS(t) P i g Please note that T i may have Coxian transitions that are not in T C;i : these are the Coxian transitions which put tokens into P i without taking token from it.

Extended Local reachability graphs
From the distinction between local and global rings, we deduce that the generator Q of the embedded CTMC of the net is a sub-matrix of A j (t) 3 5 (1) the tensor sum corresponds to local rings whereas the second and third products correspond to global rings. The expression "sub-matrix" means that if x = (m i ) 1 i m is a tuple of reachable extended markings, then for any tuple of extended markings y, q x;y = q 0 x;y if y is reachable, and q 0 x;y = 0 if y is unreachable.

Computation of the C j (t)
The computation of the C j (t) matrices is given by the algorithm 3.1 in which we have split the Coxian transitions case (second term of Q 0 ) and the T 0 transitions case (third term of Q 0 ).
Note that in the rst case, we write m j t ?! m 0 j since we do not care about the fact that we have an -ring of t.

Evaluation of the method
We study the savings induced by our method with respect to the two other ones 12] and 7]. Let us rst remark that whatever the method, we have to solve the same CTMC. Hence the comparison has to be done at this level: what are the savings using the decomposition method during the steady state probabilities computation of the embedded CTMC instead of the "standard" method? As complexity measure we simply take the sum of the cardinalities of the involved state spaces (denoted by S for the standard method and D for our method): for large models, this is the most restrictive parameter and the temporal complexity strongly depends on it (moreover, our decomposition method allows an e ective parallel computation 19,15] which in fact may reduce the extra computation time needed by tensor operators). The comparison between the two methods is then carried out studying the ratio R = S D . The diversity of possible interactions between Coxian transitions and others elements of the net forbids any general comparison. We propose, as a rst step, to study a classical situation for two of the policies de ned in Sect. 2: Single Server and identi ed clients. The generic net (Fig. 4) includes only one Coxian transition, t x controlled by a single place p 1 . The maximal marking of p 1 is c. t x has s stages.The rest of the net consists of n = jPj ? 1 places, no one controlling a transition in con ict with t x . We may hence choose the partition P 1 = fp 1 g and P 2 = PnP 1 (jP 2 j = n).
In all cases we approximate the number of markings of a net with p places and j tokens by its upper bound p+j?1 p?1 . . Hence R > 1, but lim n!1 R = 1: intuitively, if n increases faster than c, the in uence of the Coxian transition becomes negligible so that the e ects of the decomposition method, designed to cope with the Coxian transitions impact, are reduced. The direct method may be used in this simple case.

Single Server policies
If c n, R 1 + (s ? 1) i n n! c n n! and, in this case, lim c!1 R = s: if the number of clients increases faster than the complexity of their behaviour in the net, the in uence of their position in the Coxian distribution increases and our method takes advantage of this situation.

Policies with identi ed clients
We assume an Age Memory policy. n! and increases very fast with c, for a xed n: for n = 10; s = 2, for instance we have for c = 10, R 200, for c = 20, R 10 5 and for c = 30, R 5:10 6 .
Clearly, with Multiple Servers and identi ed clients policies, the decomposition method is far better than the standard method, by transforming a "product of complexities into a sum".

Introduction of immediate transitions
We extend the previous method to Coxian GSPNs with immediate transitions and Coxian transitions with a possible immediate ring i.e., 0 (t) 6 = 0 (immediate Coxian transitions).
This introduces two problems: what is the semantics of immediate Coxian transitions and how to obtain a tensor decomposition of the underlying CTMC of the net?

Stochastic semantics of immediate Coxian transitions
In this section, we de ne the supplementary properties of the stochastic model which allows us to fully specify the stochastic process underlying the net.
Since an immediate Coxian transitions may be seen as an immediate transition, its priority must be de ned accordingly: it must be 1 if all immediate transitions have priority 1, or else any positive value. The priority of a non-immediate, Coxian transition is 0.
We have three situations of con ict which must be speci ed. 1. Con ict with at least one immediate transition: from a modelling point of view, the con ict must be solved between immediate rings only. We attribute a weight (t) to the immediate Coxian transitions as for standard immediate transitions and solve the con ict like for standard GSPNs. 0 is not used to solve the con ict. 2. Con ict with exponential or non-immediate Coxian transitions only: we must be able to choose between the 0 and 0 rings of t: with probability 0 (t), t res immediately. 3. Con ict between L immediate Coxian transitions and exponential or non-immediate Cox-ian transitions: we have L + 1 possible rings; L immediate rings and one ring where all transitions engage one client in their rst exponential stage. Here also, the weights (t l ) allow to select the immediate ring; this choice is then weighted by 0 (t l ).
Pr(immediate ring of t l ) = (t l ) P L l=1 (t l ) 0 (t l ) and Pr(engagement in the rst stage of the (t l ) 1 l L ) = P L l=1 (t l ) 0 (t l ) P L l=1 (t l ) Let us point out that although the de nition of these parameters may seem a blocking to the usability of the method, we expose here the most general case, and in practical situations, the modeller will have to give only a few of them.

Decomposition of the net
Since a tensor decomposition may only be applied to a CTMC, and as rings of immediate transitions and of 0 rings of immediate Coxian transitions (Coxian immediate ring) do not occur in the tangible reachability graph (TRG), the modi cations resulting from these rings must be local to subsets P i .
Hence we impose that all immediate and immediate Coxian transitions are local. Notice that, given a partition (P i ) it is often easy to change some subsets to transform a non local immediate or immediate Coxian transition into a local one as in Fig. 5 where we replace P t and P 2 by P 2t which makes t a P 2t local transition.

Solution of the Markov chain
We describe now the modi cations of the di erent steps of the non immediate case. We refer the reader to 3, 4] for a general presentation of the extraction of the embedded CTMC of a GSPN.
The building of the reachability graph takes into account the new semantics. We shall yet call "tangible" a marking from which there is no immediate neither Coxian immediate ring. The edges with non tangible origin have to be labeled with transition probabilities deduced from above discussion. 5. 3.1 Computation of the TRG i TRG i is the graph of the tangible extended local markings of the ith sub-net. Note that the same descriptors are used since all states are tangible.
The algorithms are globally the same, but after an exponential or ring, we compute the sequences of immediate and Coxian immediate rings generated in each of the subnets: in the net of Fig. 6 (top) for instance, the ring of t 11 (exponential) produces immediate rings inside three subnets. The of all sequences of immediate rings between these markings. Note also that during these computations, the descriptors of the involved transitions are updated.

Tensor expression of the generator
We have the same expression as for the non immediate case. However, matrices include now the e ects of the immediate rings. The corresponding probabilities are retrieved from the TRG i .

Conclusion
We have presented in this paper a new method to introduce PH distribution ring time in generalized stochastic Petri nets. It allows the speci cation of all types of stochastic semantics usually encountered, including immediate transitions and PH transitions with immediate ring. We have shown how to de ne a decomposition of the net which allows the derivation of a tensor expression of the generator. This expression reduces the increase in complexity induced by the PH transitions, keeping it of the same order of magnitude as the complexity of the reachability graph of the net when PH transitions are weakly coupled with other parts of the net. Future work will investigate more deeply the savings provided with respect to the topology of the net and evaluate the possible integration of the method in a tool like 8]. However we would like to stress the novelty of our approach: previous works addressed mainly the problem of being able to de ne PH distributions in the SPN framework, without providing any suggestion on how to handle the additional complexity induced on the underlying Markov chain. To the best of our knowledge this is the rst attempt to relate the introduction of non exponential distributions to the e cient solution techniques based on tensor expressions. On the other hand, tensor algebra based techniques, although potentially very e cient are usually di cult to apply to general SPN structures. The application to PH expansion represents a case in which the partition of the state space in Cartesian products of nearly independent parts is suggested by the structure of the model in a simple way.