1
arXiv:2606.02007v1 [cs.DC] 1 Jun 2026
TAPAAL SMC: Statistical Model Checking of Stochastic Timed-Arc Petri Nets Tanguy Dubois Nantes Université École Centrale Nantes, CNRS, Nantes, France
Kim G. Larsen Aalborg University Department of Computer Science, Aalborg, Denmark
Jiřı́ Srba Aalborg University Department of Computer Science, Aalborg, Denmark
Abstract. Timed-Arc Petri net (TAPN) is a timed extension of the classical Petri net model where tokens have their age and input arcs are associated with time intervals restricting the ages of tokens available for transition firing. Additionally, a TAPN can also contain place invariants constraining the ages of tokens in places, inhibitor arcs preventing a transition from firing and transport arcs that preserve token ages upon firing. This set of features, as much as it allows us to model complex systems, also often makes verification problems computationally hard or even undecidable. Moreover, in order to model real-life examples, additional stochastic aspects are often necessary to capture the desired behaviour. We suggest the first stochastic semantics for TAPNs and design and implement the quantitative and qualitative Statistical Model Checking (SMC) algorithms in the model checker TAPAAL. We argue for the semantic choices we made in the stochastic semantics and prove that the semantics is well-behaving. On a number of case studies we demonstrate the practical applicability of our modelling formalism and its SMC implementation.
1.
Introduction
In many complex systems, exact deterministic models fail to capture the inherent uncertainty and variability present in real-world processes. Introducing stochastic elements to formal models allows us to
2
Dubois, Larsen and Srba / TAPAAL SMC: Statistical Model Checking of Stochastic Timed-Arc Petri Nets
better represent these unpredictable dynamics, providing a more realistic framework for analysis [1]. However, the exact verification of these models often becomes infeasible due to the sheer computational complexity involved. Instead, statistical model checking [2] offers a practical alternative: it provides approximate but computationally efficient results. While statistical methods may not guarantee absolute precision, they offer a valuable trade-off between accuracy and performance, enabling meaningful insights into systems that are otherwise too complex to analyze. Research on introducing stochastic semantics to Petri nets focuses on incorporating randomness into the modelling of systems with uncertain or probabilistic behaviour. Stochastic Petri nets [3, 4] extend traditional Petri nets by assigning probabilistic firing rates to transitions, enabling the modelling of systems with random delays, such as in queuing or communication networks. This model is further extended to generalized stochastic Petri nets [5] which include also immediate transitions and further extensions with e.g. colored tokens [6]. A complex timing behaviour can be modelled by stochastic time Petri nets [7] where apart from the stochastic aspects, every transition also contains a firing interval that defines the earliest and latest firing date of the transition. As there are only finitely many transitions (and hence we need to use only a fixed number of clocks), the SMC engine of UPPAAL is used in [7] to perform statistical model checking. On the other hand, the model of timed-arc Petri nets (TAPN) introduced in [8, 9] (see also an overview paper [10]), requires a new clock for each token in the net and the total number of clocks cannot be a priori fixed. TAPNs equipped with transport arcs are more expressive than timed Petri nets and timed automata [11, 12], however, to the best of our knowledge, a stochastic semantics for this model has not been suggested yet. In this paper, we consider the classical TAPN model extended with additional features like transport arcs, age invariants and inhibitor arcs as used in the model checker TAPAAL [13]. We suggest the first stochastic semantics to this extended model, implement a new statistical model checking engine as a part of the TAPAAL suite and evaluate it on a number of realistic case studies. We finish the introduction by giving an intuition of our stochastic extension of TAPNs by describing a simple producer/consumer system depicted as a stochastic TAPN in Figure 1. The net in the figure models a producer that is either in the place Ready or Resting. As soon as the transition Produce is enabled (e.g. in the initial marking that contains one token of age 0 in the place Ready), we sample its firing date d ∈ R≥0 from normal distribution with mean 1.5 and standard deviation 0.3. The system then delays for d time units, followed by firing the transition Produce which deposits new tokens of age 0 into the places Resting and Store. At this moment, two transitions become newly enabled, namely Recover and Consume1, and we sample their firing dates from uniform distribution between 0.5 and 1.5, and exponential distribution with rate 0.4, respectively. The system now delays until the next interesting event occurs, meaning that either one of the two enabled transitions reaches its firing time, or a new transition becomes enabled. In the latter case, after three time units the transition Consume2 becomes enabled as the token’s age now fits into the interval [3, 8] and we sample its firing date from the exponential distribution with rate 0.6. Now three transitions are enabled and the one with the earliest scheduled firing time will fire, unless another 5 time units pass and the transition Consumer2 gets disabled and its scheduled firing time is deleted. Should no transition be scheduled for an additional time unit, the transition Expire also becomes enabled. Its firing date is sampled from the constant distribution with mean 0, meaning that it is scheduled to fire immediately and it will move the product from Store to Waste while preserving
Dubois, Larsen and Srba / TAPAAL SMC: Statistical Model Checking of Stochastic Timed-Arc Petri Nets
3
Consume1 exponential(0 .4 ) M = Random
Ready 0
[0 , ∞)
[0 , ∞) Store Recover uniform(0 .5 , 1 .5 ) [0 , ∞) Resting
Produce normal(1 .5 , 0 .3 )
[3 , 8 ] [9 , ∞)
Expire constant(0 ) W =∞
Consume2 exponential(0 .6 ) M = Youngest
Waste
Figure 1: Stochastic TAPN modelling a produced with two consumers
its age 9 (the special arrow types are the transport arcs which move tokens without changing their ages). The transition Expire has weight ∞, meaning that in case that at the same time other transitions are also enabled and scheduled to fire, Expire will get the priority in firing (assuming that the default weight of all other transitions is 1). Finally, as the place Store can possibly contain multiple tokens of different ages, the firing modes of transitions Consume1 and Consume2 specify which tokens will be consumed during the firing. For Consume1 we shall use a randomly selected token that enables the transition, however, Consume2 will always use the youngest available token in the interval [3, 8]. We may now ask what is the probability that there appears a product in the place Waste within the first 20 time units of execution of the system, which we can, with the help of our SMC tool that executes 461 121 random runs, estimate to have the probability 0.045 ± 0.002 with 95% confidence. This probability will approach the value 1 as the time horizon gets longer.
2.
Related work
There are two main approaches for the verification of probabilistic systems. Symbolic probabilistic verification methods [14] aim to produce the exact probability of observing a given property by computing or approximating the probabilistic measure. The PRISM tool [15] is known to implement such methods for finite state systems. On the other hand, statistical model-checking (SMC) [2] generates random runs to produce an estimation of the evaluation of a property. While symbolic methods can be very costly for exploring the stochastic state space, the latter is much faster at the cost of precision and confidence. However, SMC is known to become costly and imprecise when it comes to verifying rare events, and methods such as importance splitting [16] and importance sampling [17] have been introduced to solve this issue.
4
Dubois, Larsen and Srba / TAPAAL SMC: Statistical Model Checking of Stochastic Timed-Arc Petri Nets
Semantics
Time Constraints
Type
sPN (ORIS) [4]
No
Strong
DSPN (TimeNet) [23]
No
Strong
HPnGs (HYPEG) [24]
No
Strong
sTPN (UPPAAL) [7]
Yes
Strong
sTAPN (TAPAAL)
Yes
Weak
Figure 2: Stochastic timed semantics for Petri nets
For stochastic timed systems—where both discrete choices and choices of real-time delays are made stochastically—verification is generally undecidable, making SMC indispensable. An exception is the probabilistic semantics for one-clock timed automata proposed in [18], where verification of the very restricted subclass of almost-sure properties was shown to be decidable. For networks of general timed automata with an arbitrary number of clocks, a stochastic semantics was given in [19], providing the foundation of UPPAAL’s SMC engine [20]. The first and simplest stochastic extension of Petri nets (sPN) is detailed in [3]. Every transition becomes a stochastic component and samples its firing date according to a given distribution, resulting in a race where the transition with the earliest firing date wins and fires. This semantics is simple to understand, but only handles transition delays without any timing constraints and every sampled date is considered valid. Verification of sPN can be executed, e.g., using the ORIS tool [21]. The basic stochastic Petri nets were extended with immediate transitions and transitions with deterministic delays into the Deterministic and Stochastic Petri nets (DSPN) [22], which are implemented in tools like TimeNet 4.0 [23]. Hybrid Petri nets with multiple general transitions (HPnGs) [24] also support these features, but also allow for continuous transitions. The verification of HPnGs can be done using the tool HYPEG (formerly libhpng) [25]. Both essentially use the same method for sampling dates, and thus have the same issues with missing time constraints, however, by mixing stochastic transitions with deterministic timed transitions and continuous transitions, this type of stochastic nets can model more complex timed behaviour. There exists a stochastic extension of the Time Petri net (TPN) semantics, which achieves even more expressive power by associating a time firing interval to each transition [7]. These nets implement the strong semantics [26] that guarantees that a transition must be fired within its firing interval and the distributions used for the transitions are essentially scaled to match the firing intervals. An implementation using the UPPAAL SMC engine was initially used to perform SMC on these semantics [7]. In stochastic TPN time, timing constraints are supported but are limited to the transitions in the net, which cannot be dynamically changed. On contrary, timed-arc Petri nets (TAPN) [8, 9], studied in this paper, associate the timing information to a potentially unbounded number of tokens in the net and input arcs to transition restrict the ages of tokens that can be consumed by transition firing. Traditionally, TAPNs rely on the weak semantics meaning that it is allowed to perform time delays that can disable currently enabled transitions. Further extensions with urgent transitions or age
Dubois, Larsen and Srba / TAPAAL SMC: Statistical Model Checking of Stochastic Timed-Arc Petri Nets
5
invariants [10] are needed to enforce urgent behaviour. Moreover, an extension of the timed-arc Petri net model with only transport arcs already makes the model more expressive than other types of timed nets and timed automata [11, 12]. No stochastic semantics for TAPNs or other Petri net models that includes the weak (nonurgent) semantics has been given yet. This is, including an efficient implementation of the framework, the main contribution of this paper. An overview table of selected semantics and respective tools is given in Figure 2.
3.
Timed-Arc Petri Nets
We shall first introduce the semantics of timed-arc Petri nets without the stochastic attributes. Let I be the set of timed intervals of the from [a, b] where a ∈ N0 , b ∈ N0 ∪ {∞} and a ≤ b. Let time invariants Iinv ⊆ I be a subset of timed intervals of the form above where a = 0. Definition 3.1. (Timed-arc Petri net) A timed-arc Petri net (TAPN) is a 7-tuple N = (P, T, IA, OA, Transport, Inhib, Inv) where • P is a finite set of places, • T is a finite set of transitions such that P ∩ T = ∅, • IA ⊆ P × I × N × T is the set of input arcs that connect a place to a transition and are annotated by a time interval and an arc weight, such that if (p, I, w, t) ∈ IA and (p, I ′ , w′ , t) ∈ IA then I = I ′ and w = w′ , • OA ⊆ T × N × P is the set of weighted output arcs that connect transitions to places such that if (t, w, p) ∈ OA and (t, w′ , p) ∈ OA then w = w′ , • Transport ⊆ IA × OA is the set of transport arcs such that whenever ((p, I, w, t), (t′ , w′ , p′ )) ∈ Transport then t = t′ and w = w′ , and if (α, β), (α, β ′ ) ∈ Transport then β = β ′ and symmetrically if (α, β), (α′ , β) ∈ Transport then α = α′ , • Inhib ⊆ P × N × T is the set of weighted inhibitor arcs, and • Inv : P → Iinv is the function assigning age invariants to places. We note that the definition implies that a given place and a transition cannot be connected by both a normal and a transport arc at the same time. For a transition t, we denote by Pre(t) = {p ∈ P | (p, I, w, t) ∈ IA} the set of input places, and by Post(t) = {p′ ∈ P | (t, w′ , p′ ) ∈ OA} the set of output places. A marking represents the distribution of tokens together with their ages across the places in the net. Definition 3.2. (Marking) A marking M is a finite multiset over P × R≥0 such that (p, x) represents a token of age x in the place p. A marking M is valid if (p, x) ∈ M implies that x ∈ Inv(p) for every place p ∈ P . We denote the set of all valid markings by M.
6
Dubois, Larsen and Srba / TAPAAL SMC: Statistical Model Checking of Stochastic Timed-Arc Petri Nets
For a marking M ∈ M, we use the notation M (p) = {x | (p, x) ∈ M } to denote the multiset of token ages in the place p. Definition 3.3. (Enabled transition) A transition t ∈ T is enabled in a marking M if there exist multisets of tokens In ⊆ Pre(t) × R≥0 and Out ⊆ Post(t) × R≥0 such that • In ⊆ M , i.e. the In set contains only tokens present in the marking M , • for every (p, I, w, t) ∈ IA and every x ∈ In(p) we have x ∈ I and |In(p)| = w, i.e. each input arc has exactly w tokens in its input place, all of them satisfying its time guard, • |M (p)| < w for every (p, w, t) ∈ Inhib, i.e. none of the inhibitor arcs has enough tokens in its input place to inhibit the firing of t, • for every ((p, I, w, t), (t, w, p′ )) ∈ Transport we have In(p) = Out(p′ ) and x ∈ Inv(p′ ) for every x ∈ In(p), i.e. all w-many tokens in p can be moved (including their ages) to p′ and still satisfy the age invariant of p′ , and • for any (t, w′ , p′ ) ∈ OA that does not appear in Transport we have |Out(p′ )| = w′ and x = 0 for every x ∈ Out(p′ ), i.e. every output arc that is not a transport arc must create exactly w′ new tokens of age 0. We denote the set of transitions enabled in a marking M by en(M ). Definition 3.4. (Transition firing) If t is enabled in a marking M by the multisets of tokens In and Out then t can fire and produce a marking M ′ = (M \ In) ∪ Out t
where \ and ∪ are operations on multisets, and we write M → − M ′. Definition 3.5. (Delay) Let M be a marking and let d ∈ R≥0 be a real number. We can delay d time units from M if x + d ∈ d
Inv(p) for all (p, x) ∈ M , and write M − → to indicate that a delay is possible. If the delay is possible, δ
we write M → − M [d], where M [d] is a marking defined by M [d] = {(p, x + d) | (p, x) ∈ M }. A run of a TAPN from an initial marking M0 is an alternating sequence of time delays and transition firings such that d
t
d
t
d
t
0 0 1 1 2 2 M0 −→ M0 [d0 ] − → M1 −→ M1 [d1 ] − → M2 −→ M2 [d2 ] − → ...
A run is maximal if it is either infinite or ends in a deadlock, i.e. in a marking M such that for any possible delay d we have en(M [d]) = ∅. By runs(N ) we denote the set of all maximal runs of N .
Dubois, Larsen and Srba / TAPAAL SMC: Statistical Model Checking of Stochastic Timed-Arc Petri Nets
3.1.
7
Logic for Reasoning about Runs
We shall now define a logic that allows us to argue about marking properties. A marking property φ is a Boolean combination of atomic propositions of the form p ▷◁ n where p ∈ P , n ∈ N0 and ▷◁ ∈ {<, ≤, =, ≥, >} with the obvious semantics that a marking M satisfies an atomic property p ▷◁ n if and only if |M (p)| ▷◁ n, i.e. the number of tokens in a place p satisfies the constraint imposed by the atomic proposition. This is naturally extended to Boolean connectives and we write M |= φ if a given marking M satisfies the marking property φ. Let Φ be the set of all marking properties. An example of marking property is p1 ≥ 5 ∧ p2 = 1 which holds in a marking that contains at least five tokens in p1 (irrelevant of their ages) and exactly one token in the place p2 . Note that properties that involve the ages of the tokens can be encoded by introducing “monitoring” transitions with appropriate time intervals and checking their enabledness. A run of a TAPN satisfies the (reachability) formula F φ if it contains a marking M such that M |= φ, and it satisfies the formulate Gφ if every marking M on the run satisfies M |= φ.
4.
Stochastic Timed-Arc Petri Nets
We shall now define the model of stochastic timed-arc Petri nets (sTAPN). Let N = (P, T, IA, OA, Transport, Inhib, Inv) be a TAPN as defined in the previous section. We extend it with stochastic features by adding three additional functions1 assigning (i) density function to transitions (for sampling firing delays), (ii) weights to transition (for the resolution for firing conflicts) and (iii) firing mode to transitions (for deciding which tokens are used in transition firing). The function density : T → F (R, R≥0 ) assigns a probability density function2 for each transition. It is used to sample firing dates of newly enabled transitions. When we sample a delay for a transition t ∈ T according to density(t), we assume that we choose a random value according to the distribution, except for the case when the sampled value is negative, in which case we return the delay of 0. The function weight : T → R≥0 ∪ {∞} assigns a weight (priority mass) to each transition. The weights are used in the event of a firing date collision in case several transitions sample the same firing date. Let T ′ ⊆ T be the set of transitions that chose the same firing date then the probability of firing t ∈ T ′ is given by a weighted uniform choice weight(t) ′ t′ ∈T ′ weight(t )
P
0 r where we postulate that ∞ ∞ = 0 = 1 and ∞ = 0 for any r ∈ R≥0 . If one of the competing transitions has an infinite weight then it will always win the competition with firing probability 1, except for the situation where there are several transitions with infinite weight—in this case we choose uniformly between them. On the other hand, if one of the competing transitions has a zero weight, then it will never be chosen, unless every other competing transition also has a zero weight—in that case we choose again uniformly among the zero-weight transitions. 1 2
These functions are defined over the domain of real numbers R but in the actual implementation are represented as doubles. Where for each t ∈ T , the area under the function is one: R density(t)(θ)dθ = 1.
8
Dubois, Larsen and Srba / TAPAAL SMC: Statistical Model Checking of Stochastic Timed-Arc Petri Nets
Finally, the firing mode function mode : T → {Youngest, Oldest, Random} is a function that determines for each transition which tokens are consumed when the transition is fired in a given marking. Assume that t ∈ T is enabled in M by the multisets of tokens In 1 , . . . , In n (enumerating all the possiblities). If mode(t) = Youngest then the transition t will be fired using the set In i , 1 ≤ i ≤ n, that minimizes the sum of ages of all tokens in In i and similarly if mode(t) = Oldest then the sum will be maximized. In case that that there are several such sets that minimize/maximize the sum, we uniformly choose one such set. In case mode(t) = Random, we select In i with uniform probability i/n. Our stochastic TAPN model uses a single-server policy [27], meaning that once a transition fires, it starts its firing process again (if enabled); and an enabling memory policy [27], meaning that a firing date of a transition is stored as long as the transition remains enabled, and the firing date gets forgotten as soon as the transition becomes disabled. Transition firings that both consume and produce a token to some place do not change enabledness of other transitions that also consume from such a place. In other words, transition firing is considered as an atomic event that takes no time. Figure 1 shows an example of a stochastic TAPN, demonstrating all the different features of the stochastic extension, including different density functions based e.g. on normal and exponential distributions, and weights assigned to transitions, meaning that the transition Expire has a firing priority in case that another of the remaining four transitions (with the default weight 1) becomes scheduled at the same time point. The stochastic net also shows examples of different firing modes as the place Store can possibly contain several tokens and the transition Consume1 will in such case consume a token of a random age, while Consume2 will consume a token with the smallest available age.
4.1.
Algorithm for Generating Random Runs
The heart of an SMC algorithm is the generation of random runs from an initial marking M0 , until we reach a deadlock, a given time or step (number of transition firings) bound, or a marking that satisfies a given marking property. The generation of random runs is executed as follows. • Newly enabled transitions randomly sample their firing date according to the corresponding density function. • The net then delays to the next interesting moment where either one of the enabled transitions is scheduled to fire or the enabledness of transitions changes and their scheduled firing dates get updated (if a transition becomes disabled, it is unscheduled and if it becomes newly enabled, we sample according to its density function). • If several transitions are scheduled to fire at the same date, we select the winner according to the function weight and fire the transition using the tokens selected by the function mode. The winner then resets its firing date. • By firing the winning transition, other transitions may become enabled or disabled. We update their scheduled firing dates accordingly and repeat the whole process until a marking property is satisfied (in which case we return true) or until we exceed the given number of steps (transition firings) or the given time horizon. In this case we return false.
Dubois, Larsen and Srba / TAPAAL SMC: Statistical Model Checking of Stochastic Timed-Arc Petri Nets
9
Algorithm 1 RandRunN (M0 , φ, c, s) 1: Input: Marking M0 , property φ, time bound c ≥ 0, and steps bound s ≥ 0 2: Output: Boolean indicating if a random run generated from M0 contains a marking satisfying φ
before c time units passed and k transitions were fired 3: for all t ∈ T : scheduled(t) ← sample from density(t) if t ∈ en(M0 ), else ∞ 4: AccumulatedDelay ← 0; AccumulatedSteps ← 0; M ← M0 5: while AccumulatedDelay ≤ c and AccumulatedSteps ≤ s do 6: 7:
if M |= φ then return true select the smallest positive delay d, 0 < d ∈ R≥0 , satisfying (i) en(M ) ̸= en(M [d]), or (ii) ∃δ > 0.∀ε ∈ (0, δ]. en(M ) ̸= en(M [d + ε]).
8: 9: 10: 11:
if no such d exists and scheduled(t) = ∞ for all t ∈ T then return false ▷ We are in a deadlock else d ← min{ d, min scheduled(t) } ▷ Select an earliest interesting delay t∈T
end if 13: M ← M [d] ▷ Advance to the next interesting time point 14: AccumulatedDelay ← AccumulatedDelay + d 15: for t ∈ T s.t. t ∈ en(M ) ∧ scheduled(t) = ∞ do ▷ Schedule newly enabled transitions 16: sample a value v from density(t) and set scheduled(t) ← v + d 17: end for 18: candidates ← {t ∈ T | scheduled(t) − d = 0} ▷ Transitions that can fire 19: if candidates is not empty then 20: randomly select tfired ∈ candidates using the weight function weights tfired 21: let M −− → M ′ according to mode(tfired ); M ← M ′ 22: scheduled(tfired ) ← ∞ ▷ The executed transition is unscheduled 23: AccumulatedSteps ← AccumulatedSteps + 1 24: end if 25: for t ∈ T do ▷ Update the scheduled firing times of each t ∈ T 26: if t ∈ / en(M ) then 27: scheduled(t) ← ∞ 28: else ( scheduled(t) − d if scheduled(t) ̸= ∞ 29: scheduled(t) ← sample from density(t) otherwise 30: scheduled(t) ← 0 if scheduled(t) < 0 31: end if 32: if scheduled(t) > 0 and ∃δ > 0.∀ε ∈ (0, δ]. t ∈ / en(M [ε]) then 33: scheduled(t) ← ∞ ▷ Unschedule t as it will get disabled 34: end if 35: end for 36: end while 37: return false 12:
10
Dubois, Larsen and Srba / TAPAAL SMC: Statistical Model Checking of Stochastic Timed-Arc Petri Nets
The random run generation is formally described in Algorithm 1. The algorithm works as follows. 1. While the accumulated delay and steps are under the specified bounds, we repeat until reaching a marking satisfying the given marking property, or a deadlock: (a) We look for the next smallest interesting delay, which is a date where the enabledness of some transition changes (line 7). In particular the condition (ii) expresses the fact that when an age of a token reaches the upper bound of a time interval on some arc, an arbitrarily small delay can disable a currently enabled transition. (b) If we cannot find such delay, and there are no scheduled transition firings, we reached a deadlock and terminate (line 8). (c) Otherwise, we delay up to the minimum between the smallest interesting delay and the earliest scheduled firing date (line 11). We also sample the firing dates for the newly enabled transitions after the delay d at line 15 (note that we add d to the sampled dates as the transition enabledness dates are not shifted by d yet). (d) If more than one transition is scheduled to fire at this date, we randomly choose a winner tfired using the weight function, and we fire it by selecting a multiset of consumed tokens according to mode(tfired ) (lines 18 to 24). (e) We update the scheduled firing dates of each transition: disabled transitions are scheduled at ∞, the firing date of scheduled transitions is shifted according to the delay and newly enabled transitions are sampled according to their density functions (lines 25 to 31). The check at line 30 truncates negative sampled values to zero. (f) We check if any future-scheduled transition will get disabled after an arbitrary small delay and unschedule such a transition (line 32). We shall first point out that we do not need urgent transitions (once an urgent transition is enabled time cannot elapse) in the stochastic semantics, as urgency can be simulated by Dirac density function which always samples the value 0 and hence time cannot elapse as long as at least one such urgent transition is enabled. Contrary to non-stochastic semantics, age invariants in places cannot be used to enforce urgency. Once a token age in a place reaches the invariant upper bound, the execution of the random run will end in a deadlock (unless some transition is actually scheduled to fire at that time).
4.2.
Examples of Random Runs
In order to better understand the details of the random run generation algorithm, we shall now exemplify it on a number of examples presented in Figure 3. In these example, the default [0, ∞) time intervals as well as the default weight 1 and the default random firing mode are omitted. Let us first consider the net from Figure 3a. Initially, t0 is not enabled so we will delay to the next interesting event, which is after 3 time units where t0 becomes enabled. At this moment, we sample a delay d uniformly from the interval [0, 5]. If 3 + d ≤ 5 then the next interesting delay is 3 + d where we fire t0 and reach a deadlock. If 3 + d > 5 then the next interesting delay is 2 time units and the token in p0 reaches the age 5. As t0 is not yet scheduled to fire at this point and after arbitrarily small
Dubois, Larsen and Srba / TAPAAL SMC: Statistical Model Checking of Stochastic Timed-Arc Petri Nets
t0 uniform(0, 5)
p0
p1
t0 uniform(0, 30)
0
[0 , 15 ]
[3 , 5 ]
t1 uniform(0, 30) 0 [10 , 35 ] p0 Inv : ≤ 30
p1
p2
(a) Single transition t3 constant(3 .0 ) W =∞
t0 constant(3 .0 ) W = 4 .0
p0 00 00
t2 constant(3 .0 ) W = 0 .0
(b) Simple race
t0 constant(3 .0 )
p0
t1 constant(3 .0 ) W = 1 .0
t1 constant(5 .0 ) M =? [3 , ∞)
(c) Date collision
(d) Firing mode p1
t0 constant(1 .0 ) W =∞
t0 constant(1 .0 ) W =∞
t1 constant(1 .0 ) p2
p1 p0
11
0
0 t1 constant(1 .0 )
(e) Atomic firing
t2 p0 constant(0 .0 )
(f) Sequential firing
Figure 3: Examples of stochastic timed-arc Petri nets
delay t0 becomes disabled, we unschedule the firing of t0 as the check at line 32 of the algorithm succeeds. The run then deadlocks. A race behaviour between two transitions is shown in Figure 3b. Initially, t0 is enabled and it chooses a firing date d0 uniformly from the interval [0, 30]. If d0 ≤ 10 then t0 will fire as the next interesting delay is d0 where t0 is the only enabled transition. If d0 > 10 then at 10 units the transition t1 will also sample its firing delay d1 from the interval [0, 30]. Now if d0 − 10 < d1 and d0 ≤ 15 then the transition t0 wins the race. Should d0 − 10 > d1 then t1 wins the race and fires. If none of the transitions fired until 15 time units, t0 gets unscheduled. Depending on the sampled delay d1 , either t1 fires provided that d1 + 10 ≤ 30, or the whole net deadlocks if d1 + 10 > 30 due to the invariant on p0 that blocks any time progress once the token in p0 reaches the age 30. Note that the firing date collision of t0 and t1 is unlikely as the probability that the transition delays are sampled so that d0 = d1 + 10 is zero. However, in case of constant distributions, collision of firing dates is possible as demonstrated in
12
Dubois, Larsen and Srba / TAPAAL SMC: Statistical Model Checking of Stochastic Timed-Arc Petri Nets
Figure 3c. As there are four tokens of age 0 in the initial marking, each of the four transitions can fire. The scheduled firing date of each transition will be 3 and we need to resolve the probability in which order they fire. As the weight of t3 is infinity, it will be always the first transition to fire. After firing t3 , three transitions remain enabled. Transition t0 will fire next with the probability 4/5 and transition t1 will fire with the probability 1/5. After both t0 and t1 fired, t2 will fire last as its weight is 0 and it could not compete with t0 and t1 . The example in Figure 3d shows the influence of the firing mode (depicted by M = mode where mode ∈ {Random, Youngest, Oldest}) on executions. The transition t0 produces a token every 3 time units. Once the token in p0 reaches the age 3, also t1 becomes enabled and will fire after 5 time units. At that point, there will be tokens of ages 8, 5, and 2 in p0 , and t1 needs to choose which one of the token ages satisfying the time interval [3, ∞) it consumes during the firing. If mode(t1 ) = Youngest then t1 select the token of age 5. If mode(t1 ) = Oldest, it select the token of age 8. If mode(t1 ) = Random then t1 uniformly chooses between the tokens of ages 5 and 8. Finally, Figures 3e and 3f demonstrate a difference when a transition gets unscheduled during the random run generation. In Figure 3e the transitions t0 and t1 get scheduled at time 1 as they both sample from the constant distribution which always returns 1. As the weight of t0 is ∞ and the default weight of t1 is 1, the transition t0 fires first, resets the age of the token in p0 to 0, and samples its next firing date (which is again after 1 time unit). However, after the firing t0 , the transition t1 is still enabled so its scheduled firing date does not change and it fires (without any further delay) immediately after t0 . The reader may wonder that during the firing of t0 , the transition t1 got temporarily disabled and it should be resampled too. If such behaviour is desirable then it can be modelled as shown in Figure 3f. In this case the probability that the transition t1 fires is zero, as it gets always resampled after any firing of t0 (which always wins the race as its weight is ∞).
4.3.
Quantitative and Qualitative SMC Algorithms
We shall now present two SMC algorithm. The first one estimates the probability of the event that a random run satisfies a given marking property (quantitative estimation), and the second one tests a hypothesis whether this probability is larger than or equal to a given constant (qualitative estimation). The main reason for introducing qualitative estimation is that it often requires us to execute a significantly smaller number of random runs compared to the quantitative estimation. The statistical model checking algorithm presented in Algorithm 2 uses the random run generator to perform Monte-Carlo simulations for quantitative probability estimation [28]. It uses the ChernoffHoeffding bound to compute, for a given precision ϵ and a confidence level ρ, the number of runs N to be executed in the SMC algorithm. Now, a confidence-interval estimating process—that executes N runs and returns the proportion of φ-satisfying runs—will have probability larger than ρ of returning a value that is no more than ϵ away from the real (unknown) probability of φ being satisfied. The second statistical model checking algorithm presented in Algorithm 3 uses the random run generator to perform a sequential probability ratio test for qualitative hypothesis tests [28, 29]. As demonstrated in [29], the algorithm terminates with probability one, and returns a Boolean value indicating if the test is verified under the probability of false positives α (the probability of rejecting the hypothesis P(ϕ) ≥ p0 when it is true) and the probability of false negatives β (the probability of
Dubois, Larsen and Srba / TAPAAL SMC: Statistical Model Checking of Stochastic Timed-Arc Petri Nets
13
Algorithm 2 Quantitative Probability EstimationN (M0 , φ, c, s, ρ, ϵ) 1: Input: Marking M0 , property φ, time bound c ≥ 0, steps bound s ≥ 0, confidence 0 < ρ < 1,
precision 0 < ϵ < 1 2: Output: The probability ±ϵ that φ is satisfied with confidence ρ in no more than c time units and
s steps 3: N ←
ln(2/(1−ρ)) 2ϵ2
▷ The number of runs to execute
4: x ← 0
5: for i = 1 to N do
x ← x + RandRunN (M0 , φ, c, s) 7: end for x 8: return N 6:
▷ Consider true as 1 and false as 0
Algorithm 3 Qualitative EstimationN (M0 , φ, c, s, pt , δ, α, β) 1: Input: Marking M0 , property φ, time bound c ≥ 0, steps bound s ≥ 0, the probability bound
for the test 0 ≤ pt ≤ 1, indifference region width δ, probability of false positive 0 < α < 1, probability of false negative 0 < β < 1 2: Output: A boolean indicating if P(φ) ≥ pt is satisfied with the probability of false positives at most α and the probability of false negatives at most β in no more than c time units and s steps 3: p0 ← pt + δ and p1 ← pt − δ 4: r ← 0 5: while true do
x ← RandRun N (M 0 , φ, c, s) p1 1−p1 7: r ← r + x log p0 + (1 − x) log 1−p 0 β 8: If r ≤ log 1−α then return true 9: If r ≥ log 1−β then return false α 10: end while 6:
▷ Consider true as 1 and false as 0
accepting P(ϕ) ≤ p1 when it is false). The indifference region [p1 ; p0 ] (here centered around pt with the given indifference region width δ) is the region of probabilities that are not considered relevant to deduce a result for the test.
5.
Induced Probabilistic Semantics
The random run generation for a stochastic TAPN N , introduced in Algorithm 1, defines an uncountable set of possible runs and induces a probabilistic measure PN on infinite sets of runs of the σ-algebra generated by cylinders of runs of the form π = I0 .t0 . . . . .In .tn ∈ (I.T )∗
14
Dubois, Larsen and Srba / TAPAAL SMC: Statistical Model Checking of Stochastic Timed-Arc Petri Nets
where Ii ∈ I is a time interval and ti ∈ T is a transition for all i, 0 ≤ i ≤ n. Given a run d t ρ=M − → M [d] → − ρ′ ∈ runs(N ) and a cylinder π = I.t′ .π ′ , we write ρ ∈ π if d ∈ I, t = t′ and ρ′ ∈ π ′ ; for the base case of the empty cylinder ϵ we postulate that ρ ∈ ε for any run ρ. In other words, the prefix of the concrete run must agree on the transition firings with the cylinder and the concrete delays must belong to the intervals in the cylinder. To define the probability measure for our model, a state in a stochastic TAPN can no longer be reduced to a marking M . A state is now a pair (M, scheduled) where M is a marking and scheduled : T → R≥0 ∪ {+∞} is a function mapping transitions to their scheduled firing dates as used in Algorithm 1. Our aim is to define the probability PN ((M, scheduled), π) that a random run ρ generated according to Algorithm 1 belongs to the cylinder π, i.e. ρ ∈ π. For the case where the cylinder does not restrict anything (π = ϵ), we define PN ((M, scheduled), ϵ) = 1. Otherwise, we assume that π = I.t.π ′ for some π ′ ∈ (I.T )∗ . In order to define PN ((M, scheduled), I.t.π ′ ) we distinguish two cases. CASE A. If there is t ∈ T such that scheduled(t) = 0 and scheduled(t) ̸= ∞ for every t ∈ en(M ), we are ready to execute transition firing as follows. We define PN (M, scheduled), I.t.π ′ = 0 if either 0 ̸∈ I or scheduled(t) > 0. Otherwise, we define P PN (M ′ ,scheduled′ ),π ′ weight(t) · |Next(t)| (M ′ ,scheduled′ )∈Next(t) ′ P PN (M, scheduled), I.t.π = (1) weight(t′ ) t′ ∈en(M ),scheduled(t′ )=0 t
where Next(t) = {(M ′ , scheduled′ ) | M → − M ′ and scheduled′ (t) = ∞ and scheduled′ (t′ ) = ′ ′ scheduled(t ) for all t ∈ T ∖ {t}} returns all the possible markings that can be reached by firing the transition t (depending on its firing mode) together with the updated schedule where the currently fired transition t gets unscheduled. Equation (1) assumes that weight(t) ∈ R≥0 ∖ 0. The cases where weight(t) = 0 or weight(t) = ∞ are defined analogously (we choose uniformly among the scheduled transitions with weight ∞; the same is true for weight 0 in case no other transition with a higher weight is scheduled, otherwise the returned probability is 0). CASE B. Let us now consider the case where either (i) there is no transition t ∈ T such that scheduled(t) = 0, or (ii) there is some t ∈ en(M ) where scheduled(t) = ∞. In case (i) we compute the minimum delay d until some transition becomes scheduled to fire or where transition enabledness changes as follows: d = min( min scheduled(ti ), inf{d ∈ R≥0 | en(M ) ̸= en(M [d])}) . ti ∈en(M )
This closely follows the computation of minimum delay in Algorithm 1 at line 7. If d = ∞, i.e. no transition is scheduled to fire and the enabledness of transitions never changes, we reached a deadlock and postulate PN ((M, scheduled), I.t.π ′ ) = 0. In case (ii) we need to schedule some currently enabled transition before we can perform a delay and we hence set d = 0.
Dubois, Larsen and Srba / TAPAAL SMC: Statistical Model Checking of Stochastic Timed-Arc Petri Nets
15
For both case (i) and (ii), we define the set of schedules SCH(d) such that (M [d], scheduled′ ) can be reached from (M, scheduled) by a delay of duration d as follows: scheduled′ ∈ SCH(d) if and only if a) scheduled′ (t) ∈ R≥0 for every t ∈ T where t ∈ en(M [d]) and scheduled(t) = ∞, b) scheduled′ (t) = ∞ for every t ∈ T where t ̸∈ en(M [d]), or scheduled(t) > d and at the same time there is δ > 0 such that t ̸∈ en(M [d + ϵ]) for every ϵ, 0 < ϵ < δ, and c) scheduled′ (t) = scheduled(t) − d for every other t ∈ T not covered by the cases a) and b). We note that in case c) we can safely shift the scheduled time by d for every transition that is enabled, scheduled and either continues to be enabled also right after the delay d or it will be scheduled at time d. Should a transition that is scheduled strictly after d time units become disabled right after the delay d, we unschedule it in case b). Let Y D(scheduled′ ) = density(t)(scheduled′ (t)) (2) t∈en(M [d])∧scheduled(t)=∞
be a density function defined on the set SCH(d). Now the probability measure is defined recursively as follows: Z ′ PN (M, scheduled), I.t.π = D(s) · P (M [d], s), (I − d).t.π ′ ds (3) s∈SCH(d)
for the case where d ∈ I, otherwise PN ((M, scheduled), I.t.π ′ ) = 0. Here I − d for an interval I = [a, b] stands for the interval [max{0, a − d}, b − d]. Notice that Case B is applicable only finitely many times in a row due to the fact that enabledness changes after at least 1 time unit (unless there is an earlier scheduled transition firing that is handled by Case A), hence the recursive definition is well-defined. Finally, we define the probability of a trace starting from the initial marking M0 to belong to the cylinder π by PN (M0 , π) = PN (M0 , scheduled0 ), π where scheduled0 (t) = ∞ for all t ∈ T . Observing that PN ((M, scheduled), π) decreases when extending the cylinder π, it follows from classical concepts of probability theory (see e.g. [30]) that PN extends uniquely to the smallest σalgebra, SN , generated by the above cylinders. Theorem 5.1. Given a stochastic TAPN N , PN is a probability measure over the σ-algebra SN generated by the set of cylinders of the form (I.T )∗ . Proof: Let π ∈ (I.T )∗ be a cylinder. By induction on |π|, we prove that for all states (M, scheduled) the probability PN ((M, scheduled), π) ∈ [0, 1] and for any extended cylinder π.I.t with an additional interval and transition we have PN ((M, scheduled), π.T.t) ≤ PN ((M, scheduled), π). For the base case |π| = 0 the property clearly holds. For the inductive case, let π = I.t.π ′ . There are two cases according to the nature of the next delay:
16
Dubois, Larsen and Srba / TAPAAL SMC: Statistical Model Checking of Stochastic Timed-Arc Petri Nets
• If there is t ∈ T such that scheduled(t) = 0 and scheduled(t) ̸= ∞ for every t ∈ en(M ), we apply Equation 1. By the induction hypothesis we know that PN ((M ′ , scheduled′ ), π ′ ) ∈ [0, 1] which implies that also PN ((M, scheduled), π) ∈ [0, 1] as the contributions of PN ((M ′ , scheduled′ ), π ′ ) are P normalized by the number of successor states |Next(t)| and because weight(t) ≤ weight(t′ ). Clearly, extending π with an additional t′ ∈en(M ),scheduled(t′ )=0
interval and transition firing cannot increase the probability as PN ((M ′ , scheduled′ ), ϵ) = 1. • In the other case where either (i) there is no transition t ∈ T such that scheduled(t) = 0, or (ii) there is some t ∈ en(M ) where scheduled(t) = ∞, we apply Equation 3 that is shifting the interval I by some delay d into I − d. As we already argued, this shifting before a transition firing becomes enabled can happen only finitely many times and we apply here an inner induction on the number of shifting steps needed before a transition firing becomes enabled. If no more shifting is required, the property holds as argued in the previous case (transition firing). Otherwise, we by the inner induction hypothesis assume that P (M [d], s), (I − d).t.π ′ ∈ [0, 1]. We also observe that D R as defined by Equation 2 over the schedules in SCH(d) is a density function and satisfies that s∈SCH(d) D(s) ds = 1. This together with the inner induction hypothesis implies that PN ((M, scheduled), π) ∈ [0, 1]. As in the previous case, extending π to a longer cylinder cannot increase the probability. We have now defined a way to assign a probability to each set of runs described by a cylinder. As there is a countable number of cylinders, there exists an unique way to expand this probability to the σ-algebra of runs generated by such cylinders. ⊓ ⊔ Equally important fact is that our properties F φ and Gφ define measurable sets of runs, thus having well-defined probabilities as postulated by the next theorem. Theorem 5.2. Given a marking property φ ∈ Φ, the formulae F φ and Gφ describe sets of runs belonging to the σ-algebra SN , thus having well-defined probabilities PN (F φ) and PN (Gφ). Proof: Let φ be a marking property and we shall argue that the set of all runs that satisfy F φ can be described as a countable union of cylinders and hence such a set has a well-defined probabilistic measure. To do so, we notice that for any run ρ that belongs to a cylinder π = I0 .t0 . . . . .In .tn ∈ (I.T )∗ the actual delays in the run ρ do not change the number of tokens in the places in any marking that is reachable from the initial marking M0 via the sequence of transition firings t0 t1 . . . tn . This leads us to the definition of the set Cφ of all cylinders I0 .t0 . . . . .In .tn that satisfy that there is a run d
t
d
t
d
t
t
0 0 1 1 2 2 n M0 −→ M0 [d0 ] − → M1 −→ M1 [d1 ] − → M2 −→ M2 [d2 ] − → . . . −→ Mn+1
from the initial marking for some delays d0 , . . . , dn such that Mn+1 |= φ and Mi ̸|= φ for every i, 0 ≤ 1 ≤ n. This means that any run that belongs to any cylinder from Cφ satisfies the formula F φ
Dubois, Larsen and Srba / TAPAAL SMC: Statistical Model Checking of Stochastic Timed-Arc Petri Nets
17
and we can define PN (F φ) = PN (
[
π)
π∈Cφ
as a probabilistic measure of this countable union of cylinders that contains exactly all runs satisfying F φ. For the formulate Gφ, we define the probability as PN (Gφ) = 1 − PN (F ¬φ). ⊔ ⊓ The last theorem proves that the computed probabilities indeed make sense, by exhibiting an equivalence with classical TAPNs. Theorem 5.3. Let N be a stochastic TAPN and let Nd be a normal TAPN net obtained by removing all stochastic features in N . Let φ be a marking property. (i) If there is no reachable marking in Nd that satisfies φ then PN (F φ) = 0. (ii) If all reachable markings in Nd satisfy φ then PN (Gφ) = 1. Proof: For the case (i), as φ never holds in any reachable marking and the SMC algorithm is a conservative extension of reachability on classical TAPNs, there is no SMC run that satisfies F φ, implying that the probability PN (F φ) corresponds to the measure of the empty set, which is clearly equal to 0. Similarly, in the case (ii), any marking obtained during the SMC run always satisfies φ and hence PN (Gφ) is the measure of any possible run, which is clearly equal to 1. ⊔ ⊓
6.
Implementation in TAPAAL and Case Studies
The SMC run generation is implemented in C++ as part of the verifydtapn engine [31] and accepts a PNML description of stochastic TAPN together with the query and other SMC parameters and it executes Algorithm 2 and returns the computed probability as well as other useful statistics like average run length and duration, number of satisfying and violating runs and data for plotting cumulative probabilities and other plots that allow to visualize the average/minimum/maximum number of tokens in places during the simulations. The existing TAPAAL GUI [13] is extended to allow for the editing of stochastic TAPNs as well as queries, it communicates with the SMC engine and visualizes the results of SMC verification. A screenshot using our running producer/consumer example with the window showing the result with the statistics as well as the cumulative probability is shown in Figure 4. We currently support the uniform and geometric discrete distributions and dirac delta, uniform, exponential, normal, gamma, triangular and log normal continuous distributions. In order to allow the user to sample from an arbitrary distribution, we also introduce custom distributions where the user provides a list of concrete values that will be sampled on demand during the run generation. For the given set of distribution parameters, the GUI also visualizes the density function. We support three types of verification: quantitative probability estimation, using a Monte-Carlo algorithm; qualitative probability test, using a sequential probability ratio test; and a mode to generate and simulate traces, which can be either any traces, traces satisfying a property, or traces violating a property. The SMC engine has an option to use multi-threading and, if enabled, it executes one run generator per available
18
Dubois, Larsen and Srba / TAPAAL SMC: Statistical Model Checking of Stochastic Timed-Arc Petri Nets
Figure 4: TAPAAL’s GUI with a running example and its SMC verification output
core and achieves almost a linear speedup with increasing number of CPU cores. Moreover, we support the estimation of the total running time of the SMC algorithm from the given confidence and precision parameters, or alternatively suggest a precision to match a given running time and confidence level. We shall now present three case studies demonstrating the capabilities of our tool. The most recent release of TAPAAL SMC can be downloaded from https://www.tapaal.net/download/. The models used in the case studies can be found in other downloads at the bottom of the page.
6.1.
Fireflies
Firefly is a beetle that is capable of emitting light flashes in periodic intervals. In nature, some firefly species like Photinus carolinus [32] show a remarkable synchronization capabilities so that a population of fireflies starts, after a while, to flash in synchrony. All this without any leader and in a completely distributed manner. This behaviour attracted the attention among scientists, trying to explain and model such behaviour (see e.g. [33, 34]). In a simplified version, each firefly requires some amount of time in order to charge and be able to flash. The discharge (flash) can then occur after some random delay, unless the firefly notices a flash from other firefly, in which case it will join the flash (assuming that it is also charged). In Figure 5 we present a stochastic timed-arc Petri net model of the fireflies flashing behaviour. We use here the well-known colored extension of Petri nets where tokens can carry additional information (in our case the id of the firefly). Initially, there are four fireflies (tokens) in the place waiting, having ids one to four and initially of age 0. With uniform distribution between 0 to 10 seconds, the fireflies independently arrive to the place charging. Here each firefly must be charged before the transition ready moves it to the charged location. The charging takes at least one second due to the time interval [1, ∞) on the arc. After the transition ready is enabled, an additional delay is sampled from a normal distribution with mean value of 2 seconds and a standard deviation of 0.5. Once a firefly arrives to the place charged, there are two situations. Either the place flashing has no token, meaning that there
Dubois, Larsen and Srba / TAPAAL SMC: Statistical Model Checking of Stochastic Timed-Arc Petri Nets
19
Color Types: fireflies is [1 , 4 ] waiting [fireflies] 1 ′ fireflies.all
flash alone exponential(0 .2 )
Variables: x in fireflies
•• ••
1 ′x 1 ′x
1 ′x
charging [fireflies] 1 ′x
arrive uniform(0 , 10 )
charged [fireflies] 1 ′x
all done constant(0 .0 )
1 ′x
[1, ∞)
ready normal(2 , 0 .5 ) 1 ′ x
flashing
1 ′x flash jointly constant(0 .0 )
Figure 5: Colored stochastic Petri net model of fireflies
is currently no one emitting a flash. In this case the firefly waits an amount of time sampled from the exponential distribution of rate 0.2 (corresponding to the mean value of 5 seconds) and once it is schedule to fire the transition flash alone, it returns to the charging location. As a side-effect, a token (uncolored) is placed to the place flashing, which disables other fireflies from flashing alone. A token in the place flashing indicates that other fireflies in the place charged should join the flash immediately. This is done by firing the transition flash jointly which samples from the constant distribution which always returns the value 0. In other words, the transition firing is scheduled immediately and the transition hence becomes urgent (indicated by the circle in the middle of the transition). As soon as there are no further fireflies in the place charged, the urgent transition all done becomes enabled and removes the token from the place flashing. In our tool, a colored Petri net is unfolded [35, 36] to the standard timed-arc Petri net without colors by creating a copy of each place for each color and a copy of each transition for each binding of the variable x to different fireflies ids. The probability distributions of each transition are then simply overtaken by the unfolded copies of the transition. We can now ask about the probability that all fireflies synchronize within 30 seconds. This can be formulated by the query F (charging = 1 and flashing = 1 and waiting = 0) requiring that we reach a situation where there are no fireflies in the place waiting, exactly one firefly in the place charging and a token in flashing, indicating that a single firefly initiated a flash and all other fireflies (that are in the place charged) are joining in a synchronous flash. Our tool executes 18456 runs to estimate that the probability of synchronization within 30 seconds is 0.997 with 95% confidence and precision 0.01. If we increase the number of fireflies to 10, the probability of synchronous flashing within 30 seconds drops to 0.917 and we need at least 60 seconds to achieve the
20
Dubois, Larsen and Srba / TAPAAL SMC: Statistical Model Checking of Stochastic Timed-Arc Petri Nets
Figure 6: Cumulative synchronization probability for 4 (left) and 10 (right) fireflies
synchronization probability 0.998. Cumulative probability plots for 4 and 10 fireflies, produced by our tool, are depicted in Figure 6.
6.2.
Frequency Spectrum As signment in Elastic Optical Networks
Elastic all-optical networks [37] allow for fine-grained resource allocation technologies in order to schedule network traffic demands on different light frequencies inside a single optical fiber. When a new demand arrives, the spectrum allocation problem is to find for each optical fiber a sequence of available frequency slots that will carry the demand. The frequency slots must be contiguous [38] (follow each other) and the number of required slots depends on the modulation scheme [39]. After some time, a demand can be released, making the frequency slots allocated for the demand available again. This allocation/deallocation process can create a fragmentation in the allocated frequency slots, possibly resulting in a situation when a newly arrived demand is blocked (there are not enough available consecutive frequency slots to accommodate the demand). In our stochastic timed-arc Petri net model depicted in Figure 7, we study two spectrum allocation strategies: a random fit that randomly chooses from the available slots, trying to uniformly distribute the demands across the whole spectrum, and a strategy called first fit [40] that always uses the lowest possible available frequency slot. We consider three common modulation schemes 16QAM, 8QAM and QPSK that require 3, 4 and 6 frequency slots [39], respectively. In Figure 7 we present the Petri net model for 16QAM only but the other two modulations are modelled anologously, just requiring a higher number of slots. We use again the colored extension in order to allow for a more compact model description. In our model, we use the variable mod that ranges over the color type Modulation, the variable s ranging over the 21 frequency slots declared in the color type Slots and we also use the color type Occupied which is a Cartesian product of the Slots color type and the available spectrum size with values 3 to 6. The dashed circles around places denote the so-called shared places that appear also in the nets for the other modulations (not shown in our Petri net model). The shared place called free contains tokens with colors (frequency slots) that are currently available for allocation of a new demand. Initially, there are 21 tokens in the place free with color values from 0 to 20. On the other
Dubois, Larsen and Srba / TAPAAL SMC: Statistical Model Checking of Stochastic Timed-Arc Petri Nets
21
Color Types: Modulation is [16QAM, 8QAM, QPSK] Slots is [0 , 20 ] SpectrumSize is [3 , 6 ] Occupied is [Slots, SpectrumSize] Variables: mod in Modulation s in Slots block constant(0 ) W =0
demands [Modulation]
arrive exponential(0 .1 )
demands arrive exponential(0 .1 ) [Modulation]
allocate constant(0 .0 ) W =∞ s ≤ 18
1 ′ 16QAM allocate constant(0 ) s ≤ 18 1 ′ s + 1 ′ (s +1 ) + 1 ′ (s +2 )
1 ′ (s, 3 )
blocked
1 ′0 1 ′s
occupied [Occupied]
free #21 [Slots] 1 ′ Slots.all
[1, ∞) 1 ′ s + 1 ′ (s +1 ) + 1 ′ (s +2 )
1 ′ (s, 3 )
release normal(3 , 1 )
counter [Slots]
1 ′s block constant(0 .0 )
1 ′ (s, 3 ) 1 ′ s + 1 ′ (s +1 ) + 1 ′ (s +2 )
free #21 [Slots] 1 ′ Slots.all
1 ′s
1 ′ (s +1 ) s >= 18
occupied [Occupied] [1, ∞)
1 ′ s + 1 ′ (s +1 ) + 1 ′ (s +2 )
blocked
1 ′ 16QAM
1 ′ mod
1 ′ mod
1 ′ mod
init constant(0 )
1 ′ (s, 3 )
increment constant(0 ) W =0
release normal(3 , 1 )
Figure 7: Random fit (left) and first fit (right) spectrum allocation for 16QAM
hand, the shared place occupied contains tokens with colors being pairs of integers (frequency, size) such that e.g. the value 10, 4 means that the frequency slots 10, 11, 12 and 13 are occupied by a demand that requires 4 consecutive frequency slots. The random fit slot allocation (Figure 7 left) gathers the arrived demands (assuming an exponential distribution of arrivals with rate 0.1) in the place called demands. Immediately after a demand arrival, the urgent transition block becomes enabled but it has weight 0, meaning that if the transition allocate is also enabled, it will have a priority. The transition allocate binds s to some initial frequency slot (no more than 18 as the demand requires three slots) and checks whether the frequency slots s, s + 1 and s + 2 are available (i.e. there are three tokens with these colors in the the place free). If this is the case, these three tokens are removed from the place free and a new token with the color (s, 3) is added to occupied. After one second holding time, enforced by the interval [1, ∞), the transition release becomes enabled and samples its firing date from the normal distribution with mean 3 and standard deviation 1 and releases these three slots once it fires. In case the transition allocate is enabled for several bindings of the variable s to different frequency slots, the used binding is selected using a uniform distribution. In Figure 7 (right) we describe the first fit spectrum allocation algorithm that always uses the lowest available frequency slots. This is achieved by placing a token with the color 0 (the lowest frequency slot) to the place counter immediately after any demand arrival. By firing the urgent transition increment, we then repeatedly increase the color value of the token by 1 until we find an available slot (using the same modelling approach as in the random assignment) or until we reach the frequency
22
Dubois, Larsen and Srba / TAPAAL SMC: Statistical Model Checking of Stochastic Timed-Arc Petri Nets
md3 exponential(0 .02 ) arriveAverage arriveSingle exponential(0 .0828 ) exponential(0 .06 )
md2 exponential(0 .0667 )
arriveBig exponential(7 .18E − 4 )
md1 exponential(0 .32 )
18x
Cut constant(0 .0556 )
65x
Plane constant(0 .0104 )
Mold constant(0 .05 )
8x Pending
accept constant(0 .0 )
Accepted
CutMembers ReCut constant(0 .0556 )
Reproduce
Planned RePlane constant(0 .0104 )
Molded ReMold constant(0 .05 ) Remolded
ReCutMembers
RePlaned md4 exponential(0 .1 )
md7 exponential(0 .045 )
4x
md5 exponential(0 .0166 ) md6 exponential(0 .03 )
Coat constant(0 .1111 )
Coated
ApplyPrimer constant(0 .1111 )
PrimerApplied
SmallDefect exponential(0 .03 )
Sand constant(0 .2 )
Sanded
ReapplyPrimer constant(0 .1111 )
ReCoat PrimerReApplied constant(0 .1111 )
Dipped
ReSand constant(0 .2 )
ReSanded
Join constant(0 .01 )
Fixed
Dip constant(0 .0333 )
Glued
Glue constant(0 .1667 )
ReGlued
ReGlue constant(0 .1667 )
ReDip constant(0 .0333 )
ReDipped
Joined
2x Fix constant(0 .3333 )
Figure 8: Manufacturing process of wooden windows [41]
slot 18 where the transition block becomes enabled. The weight of the transition allocate is infinity, meaning that it has a priority over the transition block but in case the transition is still not enabled, the transition block will fire as the transition increment has weight 0. We can now ask our tool to compare the blocking probability of the two allocation algorithms by issuing the SMC query F blocked >= 1, checking wheater at least one demand is blocked within the first 30 seconds. With 95% confidence and precision 0.01, our tool executes 18488 runs and estimates the blocking probability to 0.075 for the first fit slot allocation strategy and 0.167 for the random fit strategy, showing an improved performance of the former strategy. This behaviour is also conformed by experiments in computer networking literature [40].
6.3. Manufacturing Process of Wooden Windows The workflow of a Hungarian wooden window manufacturing company Holz-Team Ltd details a systematic production process, from preparing individual frame components to assembling and finishing complete windows. Orders vary in size, and the process includes measures to address material defects, with early-stage issues being easier to resolve than those found later. The Petri net model is given in Figure 8 and contains timing probabilities derived from real measurements performed by the company as described in [41]. The workflow uses two types of transitions: exponentially delayed transitions for unpredictable events like order arrivals and material defects, and deterministic transitions for fixedtime production phases such as cutting, molding, and coating. For a detailed description of the process
Dubois, Larsen and Srba / TAPAAL SMC: Statistical Model Checking of Stochastic Timed-Arc Petri Nets
23
Figure 9: Average number of orders and unprocessed coated windows over a year
consult [41]. The production process faced a bottleneck where coated windows that required a hardware fix accumulated over time in the place called Coated. Indeed, using TAPAAL SMC engine we estimated that with probability 0.61 ± 0.01 and 95% confidence, the weekly production process generates more than 20 unprocessed coated windows. A solution to the bottleneck was suggested in [41] and required the hiring of an additional worker to reduce the hardware fixing operation from 20 minutes to 13 minutes. After this adjustment, the probability of storing 4 or more unprocessed windows dropped to 0.011. In Figure 9 (left) we depict, based on 50 simulations of the workflow for a period of one year (2008 hours), the number of accepted orders and the average number of unprocessed coated windows for the situation before the workflow adjustment. After the adjustment, the number of unprocessed coated windows stays on average slightly above 1 with a maximum between 3 to 5 windows, as shown in Figure 9 (right). This allowed the company to produce on average about 260 additional windows per year, increasing the productivity by almost 10%.
6.4.
Quantum Communication—Beyond Super Dense Coding
Super dense coding [42] is a quantum communication procedure that allows a sender to transmit to a receiver two classical bits of information by sending only a single qubit, provided that the sender and receiver pre-share a maximally entangled pair of qubits (called an EPR/Bell pair [43] such as |ϕ+ ⟩ = √12 (|00⟩ + |11⟩)). A recently proposed quantum protocol Beyond Super Dense Coding (BSDC) [44] builds up on this idea and operates in a synchronized slotted-time manner over a quantum channel and allows the protocol participants to communicate up to 1.5 classical bits per time slot while transmitting at most one qubit in each slot. In BSDC, time is divided into discrete slots of a given frequency where the sender either stays silent (encoding the communication of a classical bit 0) or sends a qubit (encoding the classical bit 1). A communication round over a quantum channel proceeds as follows. • The sender scans its incoming buffer of randomly generated information bits. If the scanned bit is 0, it remains silent in the given time slot. Receiver interprets no activity on the quatum
24
Dubois, Larsen and Srba / TAPAAL SMC: Statistical Model Checking of Stochastic Timed-Arc Petri Nets
channel as receiving the bit 0. • Upon the first arrival of the classical bit 1, the sender on-the-fly generates an EPR pair of qubits, stores one qubit locally in quantum memory, and transmits the other maximally entangled qubit to the receiver. The receiver detects the arriving qubit (registering the arrival of the classical bit 1) and stores it locally without disturbing the entanglement. • The sender then continues with being silent for every scanned bit 0 in its input buffer. The receiver registers this as receiving 0 in each time slot with no communication. • Once the second bit 1 appears in the sender’s buffer, the sender applies Pauli operations [43] (X, Z, or identity) to its stored qubit, chosen according to the next two bits in its buffer. If the next bit is 1, it applies the gate X and if the next-next bit is 1, it applies also the gate Z; otherwise it applies the identity operation. The sender then transmits this manipulated qubit to the receiver that interprets its arrival as the receipt of the classical bit 1 and now, holding both qubits of the EPR pair, the receiver performs the CNOT controlled gate on the both qubits, followed by the application of the gate H on the second received qubit. The two qubits that are now measured and the two resulting classical bits are the two bits encoded by the sender. In this final slot, three information bits are thus transmitted in total and the protocol can restart. The protocol relies both on the strict timing constraints (length of the frequency slot vs. the duration of the application of quantum gates, creation of the EPR pair, and the duration of measurements) as well as on probabilistic aspects due to decoherence (loss of entanglement between the two qubits). This means that with the increasing time delay between creating of the EPR pair until the measurement of the qubits is performed, the receiver has a decreasing probability that it will decode correctly the two encoded classical bits. In order to use realistic timing and decoherence information, we overtook this information from an existing IBM quantum computer [45], where the application of 1-qubit gates takes 32 ns (nano seconds), 2-qubit gates require 138 ns, a qubit measurement takes 2400 ns and the reset of the qubits before every round has estimated duration of 480 ns. Similarly, to obtain hardware-realistic decoherence times for the generated EPR pair which can be in one of the possible four Bell states [43], we model the mixed states of the EPR pair as a density matrix [46] that represents the probabilistic distribution of the four Bell states where the likehood of obtaining a wrong measurement increases over time with the coefficients T1=343.25 µs for relaxation time and T2=263.65 µs for dephasing, also taken from [45]. This is also known as Pauli-twirl [47]. We sampled a series of 1000 decoherence times from the obtained density matrix (computed separately in Python) and created a corresponding custom distribution that is imported to TAPAAL. Should the live-time of the created EPR pair exceed the sampled value from this custom distribution, we assume that instead of reading the correct Bell pair, during the measurement we obtain any of the four possible Bell states with uniform probability. The sender part of our TAPAAL model of the BSDC protocol is depicted in Figure 10. The model uses the colored extension of the model and relies on the notion of shared places, already introduced in the previous case studies. The color type bit represents two possible values of a classical bit information and Bell pair is represented as a product of two classical bits, encoding the four possible results
Dubois, Larsen and Srba / TAPAAL SMC: Statistical Model Checking of Stochastic Timed-Arc Petri Nets
1 ′ (y, z , one)
1 ′ one
generated bit [bit] Color Types: dot is [dot] bit is [zero, one] bell pair is < bit, bit > buffer is < bit, bit, bit >
Constants: f = 8000
first qubit sent [dot]
ready constant(0 )
generated bit [bit] 1 ′ zero generate zero constant(0 )
got zero constant(1 )
1 ′ (y, z , zero)
sender buffer [buffer ] 1 ′ (zero, zero, zero)
Variables: x , y, z , u, v , w in bit
wait [dot]
1 ′ (zero, zero) 0.0
[f , f ]
start [dot] 1 ′ dot
1 ′ one generate one constant(0 )
(a) Classical bit generator
communication channel [bell pair ]
(b) Send single bit
start [dot] 1 ′ dot
wait [dot]
ready constant(0 )
0.0 [f , f ] got one constant(1 )
generated bit [bit]
first qubit sent [dot]
1 ′ one 1 ′x zero → [1 , 1 ] one → [32 , 32 ] next constant(0 )
next bit [dot] 1 ′x
[bit]
got one normal(650 , 50 )
1 ′ zero
1 ′ (x , y, z ) 0.0
1 ′ (x , y, z )
next next constant(0 ) 1 ′x
1 ′y
sender buffer [buffer ] 1 ′ (zero, zero, zero)
1 ′ (one, x , y)
zero → [1 , 1 ] one → [32 , 32 ]
0.0
1 ′ (u, v , w ) 1 ′ (x , y) prepared bell pair [bell pair ]
communication channel [bell pair ] 1 ′ (x , y)
1 ′ (x , y) send qubit constant(1 )
(c) Send three bits
Figure 10: BSDC bit generation and sender behaviour
25
26
Dubois, Larsen and Srba / TAPAAL SMC: Statistical Model Checking of Stochastic Timed-Arc Petri Nets
of the measurement. We also use a buffer color type of size three where we remember the last three sent/received classical bits and use these buffers to check the error rates in case that the sent/received bits disagree due to timing conflicts or decoherence. The model also uses the constant f that is the duration of the time slot of the protocol (in nano seconds) and we vary this constant in order to explore the behaviour of the protocol under different sending frequencies. Figure 10a show a simple way how we model the incoming bits that should be transmitted. Whenever the shared place generated bit is empty, a bit zero or one is immediately inserted to the place, assuming that the bits arrive with uniform probability (both transitions in the component have the default weight 1). The two transitions are blocked until the colored token from the place generated bit is consumed by the sender. This can happen in the component in Figure 10b where we model the sending of the bits 0 and of the first bit 1. Initially, the component has a token of age 0.0 in the place start. In case that the arriving bit in the place generated bit has the value zero, we shall fire the transition got zero that has the duration of 1 ns, consumes the token from generated bit and places the control token to the place wait. As we use the transport arcs here, the age of the control token is preserved and we wait until the end of the frequency slot, after which the token is placed again into start by firing the transition ready. At the same time, we read the three bits from the place sender buffer, shift the bits by one position and remember on the last position in the buffer that we just dispatched the bit 0. In this case, no qubit is passed to the place communication channel. If on the other hand generated bit contains a token with color one, we fire the transition got one which remembers in sender buffer that we are dispatching the bit 1, it prepares the EPR pair and sends the first qubit to the quantum communication channel (the values of the qubit are irrelevant in this phase). The transition also places a token to the place first qubit sent which makes sure that once the second bit 1 arrives, we have to perfrom the encoding process described in Figure 10c instead. The duration of the transition got one is 650 ns, drawn from normal distribution with variance of 50. Here we have to perfrom the reset of two qubits (480 ns), followed by the application of the Hadamard gate H (32 ns) on the first qubit followed by the application of CNOT on both qubits (138 ns) in order to prepare the EPR pair. Upon the arrival of the second bit 1, we instead perform the encoding procedure described in Figure 10c. After the second generated bit 1 is consumed by the transition got one (which is possible only if the place first quibit sent has a token), we place a token to the place wait that as before serves as a timer in order to keep the protocol in sync at frequency f . At the same time, we fire a sequence of two transitions next and next next that consume the next two input bits from the place generated bit and encode them into a qubit that is moved to the place prepared bell pair. The application of the X and Z gates in case that the input bits are 1 causes a delay of 32 ns whereas the reading of the bit 0 has delay of 1 ns (application of the identity gate). This is depicted by the color specific intervals on the input arcs. Upon consuming the three bits from generated bit, the values are remembered in sender buffer and the generated qubit is forward to communication channel while as the same time consuming the token from first qubit sent. We shall discuss how we model decoherence of the EPR pair later on. Receiver’s behaviour as well as monitoring of the sent/received information bits is shown in Figure 11. Initially, the receiver in Figure 11a has a token of age 0.0 in place waiting. Should the length of the frequency slot be reached without any qubit arriving on communication channel, it fires the transition nothing received, resets the age of the token to 0 and updates the content of receiver buffer
Dubois, Larsen and Srba / TAPAAL SMC: Statistical Model Checking of Stochastic Timed-Arc Petri Nets
done [dot]
ready constant(0 ) [f , f ]
waiting [dot] 1 ′ dot
communication channel [bell pair ]
ready constant(0 )
[f , f ]
got first one [dot]
1 ′ (x , y) decode normal(4970 , 300 )
receive constant(1 )
0.0 1 ′ (x , y, z ) 1 ′ (y, z , one)
[f , f ] 1 ′ (x , y, z )
nothing received constant(0 )
communication channel [bell pair ]
1 ′ (u, v )
received qubit constant(1 )
1 ′ (y, z , zero)
next slot timer [dot]
27
0.0
1 ′ (u, v , w )
0.0 waiting [dot] 1 ′ dot
0.0
1 ′ (one, x , y)
1 ′ (x , y)
1 ′ (x , y)
receiver buffer [buffer ] 1 ′ (zero, zero, zero)
got first one [dot]
received pair [bell pair ]
(b) Receive three bits
(a) Receive single bit sender buffer [buffer ] 1 ′ (zero, zero, zero)
timer [dot] 1 ′ dot
0.0
0.0 1 ′ (x , y, z )
1 ′ (x , y, z )
receiver buffer [buffer ] 1 ′ (zero, zero, zero)
[f , f ] start check constant(0 )
ok constant(0 ) check [dot]
1 ′ (x , y, z )
error [dot] 1 ′ (x , y, z )
0.0
mismatch constant(0 ) W =0
receiver buffer [buffer ] 1 ′ (zero, zero, zero)
(c) Monitor
Figure 11: BSDC receiver behaviour and monitor
by remembering that the last received bit is 0. Upon the arrival of the first qubit, the receiver fires the transition received qubit, interprets it as receiving the bit 1, and places a token to done where it waits to the next time slot. It also places a token into the place got first one where it remembers that on the arrival of the next qubit, it should decode in which Bell state the two qubits were prepared. This is done in Figure 11b where the second received qubit is moved to the place received pair by firing the transition receive and at the same time removing a token from got first one and keeping the frequency timer by moving a token to done while preserving its age. Finally, receiver applies the CNOT gate on its two qubits (138 ns), performs the H gate (32 ns) and then two measurements (2400 ns each), resulting in the total processing time of 4970 ns, drawn from the normal distribution with variance of 300 by firing the transition decode and remembering the measured values in the receiver’s buffer. Next, in Figure 11c we check at the end of each time slot whether the buffered bits in the sender’s and receiver’s buffer match. Should this be the case then the token in the place check is consumed by
28
Dubois, Larsen and Srba / TAPAAL SMC: Statistical Model Checking of Stochastic Timed-Arc Petri Nets
first qubit sent [dot]
decohered [dot]
decohere custom(Pauli)
prepared bell pair [bell pair ]
introduce error constant(0 ) 1 ′ (u, v ) 1 ′ (x , y)
Figure 12: Decoherence Frequency
1-Error
2-Error
3-Error
4-Error
5-Error
Throughput
Corrected
5 µs
0.994
0.986
0.970
0.941
0.890
200 kbit/s
64 kbit/s
5.5 µs
0.492
0.356
0.236
0.146
0.082
182 kbit/s
83 kbit/s
6 µs
0.290
0.130
0.031
0.010
0.002
167 kbit/s
101 kbit/s
6.5 µs
0.311
0.138
0.031
0.009
0.002
154 kbit/s
93 kbit/s
7 µs
0.350
0.156
0.036
0.010
0.002
143 kbit/s
86 kbit/s
7.5 µs
0.382
0.176
0.045
0.013
0.003
133 kbit/s
80 kbit/s
8 µs
0.401
0.189
0.051
0.015
0.003
125 kbit/s
75 kbit/s
8.5 µs
0.426
0.205
0.057
0.017
0.003
118 kbit/s
71 kbit/s
9 µs
0.451
0.220
0.063
0.019
0.004
111 kbit/s
67 kbit/s
9.5 µs
0.477
0.234
0.069
0.021
0.004
105 kbit/s
63 kbit/s
Figure 13: Error probabilities for sending 50 bits of information for a given number of bit-flips and throughput before and after applying error correction
firing the transition ok, otherwise we fire the transition mismatch and add a token to the place error where we count the number of encountered errors. Note that the transition mismatch has weight 0 and will only fire if ok is not enabled. Also, it is enough to consider buffers of size three only as this is the maximum number of classical bits that can be sent in each time slot. In the actual model of the protocol, the monitor is slightly more complicated as it counts the exact number of mismatched bits in the two buffers, by checking a mismatch at each position of the buffer. Finally, Figure 12 shows how we model the decoherence of the EPR pair over time. Once the first quibit is sent and a token is added to the place first quibit sent, the transition decohere will sample its firing time from the custom distribution called Pauli that was mentioned before. If the second qubit is dispatched before the sampled time delay, the transition decohere becomes disabled and the receiver will recover the correct bits encoded to second qubit by the sender. Should the decoherence time be reached, the transition decohere fires and enables the transition introduce error that will consume the prepared Bell pair (u, v) and replace it with another of the possible Bell pairs (each of them with 0.25 probability). This can introduce up to 2 bit-flips once the receiver decodes the Bell pair. Based on this TAPAAL SMC model, we now run a number of SMC queries that estimate the probability that while sending 50 classical bits, there are at least X bit flips on the receiver side (query called X-Error that estimates P(F error ≥ X)). These queries for up to 5 bit-flips are evaluated for different frequencies (durations of the time slot), with precision 0.001 and confidence of 95%. With this high precision and timing in nano seconds, evaluation of such a query takes between 10
Dubois, Larsen and Srba / TAPAAL SMC: Statistical Model Checking of Stochastic Timed-Arc Petri Nets
Classical communication Quantum communication - average Quantum communication - minimum Quantum communication - maximum
2,500 2,000 Transfered bits
29
1,500 1,000 500 0
0
0.1
0.2
0.3
0.4 0.5 0.6 0.7 Time (in nano seconds)
0.8
0.9
1
·107
Figure 14: Number of transferred bits at 6 µs frequency
to 50 minutes on a Mac Studio with Apple M2 Max 12-core CPU. The memory requirements for SMC are negligible (less than 50 MB). The results are displayed in Figure 13 and the lowest (best) values of error probabilities are observed for the frequency slot duration of 6 µs. If the frequencies are shorter, there is an increasing risk of not finishing all the quantum operations during the given time slot, causing a significant increase of the error rates. On the other hand, slower frequencies suffer from increased number of errors due to the decoherence of the Bell pair. The next column in Figure 13 shows the computed idealized throughput (number of delivered classical bits per second) of the protocol, disregarding the possible bit-flips. As bit-flips are not avoidable, we can apply classical error correction. To do so, we formulate for each of the analyzed frequencies a qualitative hypothesis for sending 50 classical bits P(F error ≥ X) ≥ 0.001 in order to determine the smallest number of bit-flips X that can happen with probability strictly lower than 0.001. To identify such X, we run the query for different values of X (using e.g. the bisection method) until we find the smallest number of bit-flips for which the query evaluates to false. These tests are run with the false positive/negative parameters set to 0.001 and indifference region width of 0.0002, and testing the qualitative hypothesis takes typically less than a minute. In order to guarantee that a message with 50 classical bits is received incorrectly with a probability strictly less than 0.001, we need to be resilient for 25 bit-flips for 5 µs frequency, 12 bit-flips for 5.5 µs frequency, 6 bit-flips for the remaining frequencies. In order to achieve this, we have to use classical error correction codes which guarantee twice as high Hamming distance as the number of possible
30
Dubois, Larsen and Srba / TAPAAL SMC: Statistical Model Checking of Stochastic Timed-Arc Petri Nets
bit-flips [48]; this translates to adding additional 107 error correcting bits for 5 µs frequency, 59 bits for 5.5 µs and 33 extra bits for the remaining frequencies. The last column in the table from Figure 13 displays the corrected throughput after applying the error correction codes. As expected, the highest throughput of 101 kbit/s after applying the error correction is achieved for the frequency of 6 µs. We conclude this case study by adding observations (similarly as in the manufacturing process case study) to our TAPAAL SMC query, which monitor how many bits can be communicated using classical communication (where we send one classical bit per time slot) compared to the number of transferred classical bits using the BSDC quantum protocol (where we send at most one qubit per time slot). A TAPAAL generated plot comparing these two methods is shown in Figure 14. As the number of transferred bits in the quantum protocol depends on the actual distribution of 0 and 1 in the input message, we depict the average as well as minimum/maximum number of transferred bits for a given time duration. The plot shows on average 48% increase in the number of transferred bits in the quantum protocol compared to the classical communication using the same frequency time slot of 6 µs when sending a single bit per time slot.
7.
Conclusion
We proposed, to the best of our knowledge, the first stochastic semantics to the popular timed-arc Petri net model. The semantics is weak and race-based, implying that an enabled transition may not necessarily fire within its firing interval but can instead delay out of its enabledness zone, based on the firing date sampled from any probabilistic distribution (we currently support seven continuous distributions and two discrete ones, but the approach is not limited only to these as we also support custom user-provided distributions). This allows us to model a wide range of stochastic systems, as demonstrated in our case studies. We discussed the design and implementation of the statistical model checking algorithm and argued that it is well-defined. All algorithms are implemented in the open-source model checker TAPAAL, including numerous performance optimizations (e.g. parallel execution) as well as a support for the visualization of the statistical model checking results. Our experiments indicate that the formalism is applicable to a broad range of problems and allows us to reason about complex stochastic behaviours in an intuitive way. In the future work, we would like to extend the SMC reachability formulae to a more powerful logic (like e.g. LTL) that will allow us to reason about probability of runs that satisfy additional requirements. An optimization of our SMC algorithm for handling rare events is another line of research.
Acknowledgements. We thank Mikkel Tygesen for his help with TAPAAL GUI and in particular the visualization of the SMC results, and to Nikolaj Rossander Kristensen for his help with the generation of the custom Pauli-twirl distribution in the last case study. This paper was funded by the Villum Investigator project S4OS and the ANR project BisoUS ANR-22-CE48-0012.
Dubois, Larsen and Srba / TAPAAL SMC: Statistical Model Checking of Stochastic Timed-Arc Petri Nets
31
References [1] Gardiner C. Stochastic Methods: A Handbook for the Natural and Social Sciences. Springer Series in Synergetics. Springer Berlin Heidelberg, 2010. ISBN 9783642089626. [2] Legay A, Lukina A, Traonouez LM, Yang J, Smolka SA, Grosu R. Statistical Model Checking, pp. 478– 504. Springer International Publishing, 2019. doi:10.1007/978-3-319-91908-9\ 23. [3] Florin G, Fraize C, Natkin S. Stochastic Petri nets: Properties, applications and tools. Microelectronics Reliability, 1991. 31(4):669–697. doi:10.1016/0026-2714(91)90009-V. [4] Bause F, Kritzinger P. Stochastic Petri Nets: An Introduction to the Theory. AMS, 2013. ISBN 978-3528-15535-3. [5] Balbo G. Introduction to Generalized Stochastic Petri Nets, pp. 83–131. Springer, 2007. doi:10.1007/ 978-3-540-72522-0\ 3. [6] Chiola G, Dutheillet C, Franceschinis G, Haddad S. Stochastic well-formed colored nets and symmetric modeling applications. IEEE Transactions on Computers, 1993. 42(11):1343–1360. doi:10.1109/12. 247838. [7] Cicirelli F, Nigro C, Nigro L. Qualitative and quantitative evaluation of stochastic Time Petri Nets. In: 2015 Federated Conference on Computer Science and Information Systems (FedCSIS). 2015 pp. 763–772. doi:10.15439/2015F69. URL https://ieeexplore.ieee.org/abstract/document/7321519. [8] Bolognesi T, Lucidi F, Trigila S. From timed Petri nets to timed LOTOS. Proceedings of the IFIP WG6.1 Tenth International Symposium on Protocol, 1990. pp. 395–408. [9] Hanisch HM. Analysis of place/transition nets with timed arcs and its application to batch process control. In: Application and Theory of Petri Nets 1993. Springer, 1993 pp. 282–299. [10] Jacobsen L, Jacobsen M, Møller MH, Srba J. Verification of Timed-Arc Petri Nets. In: Černá I, Gyimóthy T, Hromkovič J, Jefferey K, Králović R, Vukolić M, Wolf S (eds.), SOFSEM 2011: Theory and Practice of Computer Science. Springer, Berlin, Heidelberg. ISBN 978-3-642-18381-2, 2011 pp. 46–72. doi: 10.1007/978-3-642-18381-2\ 4. [11] Boyer M, Roux OH. On the Compared Expressiveness of Arc, Place and Transition Time Petri Nets. Fundam. Inf., 2008. 88(3):225–249. [12] Srba J. Comparing the Expressiveness of Timed Automata and Timed Extensions of Petri Nets. In: Proceedings of the 6th International Conference on Formal Modelling and Analysis of Timed Systems (FORMATS’08), volume 5215 of LNCS. Springer-Verlag, 2008 pp. 15–32. [13] David A, Jacobsen L, Jacobsen M, Jørgensen K, Møller M, Srba J. TAPAAL 2.0: Integrated Development Environment for Timed-Arc Petri Nets. In: Proceedings of the 18th International Conference on Tools and Algorithms for the Construction and Analysis of Systems (TACAS’12), volume 7214 of LNCS. SpringerVerlag, 2012 pp. 492–497. [14] Baier C, Katoen J. Principles of model checking. MIT Press, 2008. ISBN 978-0-262-02649-9. [15] Kwiatkowska M, Norman G, Parker D. PRISM 4.0: Verification of Probabilistic Real-Time Systems. In: Gopalakrishnan G, Qadeer S (eds.), Computer Aided Verification. Springer, Berlin, Heidelberg. ISBN 978-3-642-22110-1, 2011 pp. 585–591. doi:10.1007/978-3-642-22110-1\ 47.
32
Dubois, Larsen and Srba / TAPAAL SMC: Statistical Model Checking of Stochastic Timed-Arc Petri Nets
[16] Jegourel C, Legay A, Sedwards S. Importance Splitting for Statistical Model Checking Rare Properties. 2013 p. 576. doi:10.1007/978-3-642-39799-8\ 38. URL https://inria.hal.science/ hal-01087826. [17] Jegourel C, Legay A, Sedwards S. Command-based importance sampling for statistical model checking. Theoretical Computer Science, 2016. 649:1–24. doi:10.1016/j.tcs.2016.08.009. URL https://www. sciencedirect.com/science/article/pii/S0304397516303966. [18] Baier C, Bertrand N, Bouyer P, Brihaye T, Größer M. Probabilistic and Topological Semantics for Timed Automata. In: Arvind V, Prasad S (eds.), FSTTCS 2007: Foundations of Software Technology and Theoretical Computer Science. Springer, Berlin, Heidelberg, 2007 pp. 179–191. [19] David A, Larsen KG, Legay A, Mikucionis M, Poulsen DB, van Vliet J, Wang Z. Statistical Model Checking for Networks of Priced Timed Automata. In: Fahrenberg U, Tripakis S (eds.), Formal Modeling and Analysis of Timed Systems - 9th International Conference, FORMATS 2011, Aalborg, Denmark, September 21-23, 2011. Proceedings, Lecture Notes in Computer Science. Springer, 2011 pp. 80–96. doi:10.1007/978-3-642-24310-3\ 7. URL https://doi.org/10.1007/978-3-642-24310-3_7. [20] David A, Larsen KG, Legay A, Mikucionis M, Wang Z. Time for Statistical Model Checking of RealTime Systems. In: Gopalakrishnan G, Qadeer S (eds.), Computer Aided Verification - 23rd International Conference, CAV 2011, Snowbird, UT, USA, July 14-20, 2011. Proceedings, Lecture Notes in Computer Science. Springer, 2011 pp. 349–355. doi:10.1007/978-3-642-22110-1\ 27. URL https://doi.org/ 10.1007/978-3-642-22110-1_27. [21] Paolieri M, Biagi M, Carnevali L, Vicario E. The ORIS Tool: Quantitative Evaluation of Non-Markovian Systems. IEEE Transactions on Software Engineering, 2021. 47(6):1211–1225. doi:10.1109/TSE.2019. 2917202. Conference Name: IEEE Transactions on Software Engineering, URL https://ieeexplore. ieee.org/document/8719961. [22] Marsan MA, Chiola G. On Petri nets with deterministic and exponentially distributed firing times. In: Rozenberg G (ed.), Advances in Petri Nets 1987. Springer, Berlin, Heidelberg. ISBN 978-3-540-47743-3, 1987 pp. 132–145. doi:10.1007/3-540-18086-9\ 23. [23] Zimmermann A, Knoke M. TimeNET 4.0: A Software Tool for the Performability Evaluation with Stochastic and Colored Petri Nets. In: User manual. 2007 . [24] Pilch C, Remke A. Statistical Model Checking for Hybrid Petri Nets with Multiple General Transitions. In: 2017 47th Annual IEEE/IFIP International Conference on Dependable Systems and Networks (DSN). 2017 pp. 475–486. doi:10.1109/DSN.2017.41. ISSN: 2158-3927, URL https://ieeexplore.ieee. org/abstract/document/8023146. [25] Pilch C, Edenfeld F, Remke A. HYPEG: Statistical Model Checking for hybrid Petri nets: Tool Paper. In: Proceedings of the 11th EAI International Conference on Performance Evaluation Methodologies and Tools, VALUETOOLS 2017. Association for Computing Machinery, New York, NY, USA. ISBN 9781-4503-6346-4, 2017 pp. 186–191. doi:10.1145/3150928.3150956. URL https://doi.org/10.1145/ 3150928.3150956. [26] Boyer M, Roux OH. Comparison of the Expressiveness of Arc, Place and Transition Time Petri Nets. In: Kleijn J, Yakovlev A (eds.), Petri Nets and Other Models of Concurrency – ICATPN 2007. Springer, Berlin, Heidelberg. ISBN 978-3-540-73094-1, 2007 pp. 63–82. doi:10.1007/978-3-540-73094-1\ 7. [27] Bérard B, Cassez F, Haddad S, Lime D, Roux OH. Comparison of Different Semantics for Time Petri Nets. In: Peled DA, Tsay YK (eds.), Automated Technology for Verification and Analysis. Springer Berlin Heidelberg, Berlin, Heidelberg. ISBN 978-3-540-31969-6, 2005 pp. 293–307.
Dubois, Larsen and Srba / TAPAAL SMC: Statistical Model Checking of Stochastic Timed-Arc Petri Nets
33
[28] David A, Larsen KG, Legay A, Mikučionis M, Poulsen DB, van Vliet J, Wang Z. Stochastic Semantics and Statistical Model Checking for Networks of Priced Timed Automata, 2014. doi:10.48550/arXiv.1106. 3961. ArXiv:1106.3961 [cs] version: 2, URL http://arxiv.org/abs/1106.3961. [29] Wald A. Sequential Analysis. John Wiley, 1947. [30] Ash RB, Doleans-Dade CA. Probability and Measure Theory. Academic Press, 2000. [31] Jensen P, Larsen K, Srba J, Sørensen M, Taankvist J. Memory Efficient Data Structures for Explicit Verification of Timed Systems. In: Proceedings of the 6th NASA Formal Methods Symposium (NFM’14), volume 8430 of LNCS. Springer-Verlag, 2014 pp. 307–312. doi:10.1007/978-3-319-06200-6\ 26. [32] Wikipedia. Photinus carolinus — Wikipedia, The Free Encyclopedia. http://en.wikipedia.org/ w/index.php?title=Photinus%20carolinus&oldid=1237034011, 2025. [Online; accessed 09January-2025]. [33] Ramı́rez-Ávila GM, Kurths J, Depickère S, Deneubourg JL. Modeling Fireflies Synchronization, pp. 131–156. Springer, 2019. doi:10.1007/978-3-319-78512-7\ 8. [34] McCrea Madeline EB, E RJ. A model for the collective synchronization of flashing in Photinus carolinus. J. of Royal Soc. Interface, 2022. (19). [35] Christesen N, Glavind M, Schmid S, Srba J. Latte: Improving the Latency of Transiently Consistent Network Update Schedules. In: Proceedings of 38th International Symposium on Computer Performance, Modeling, Measurements and Evaluation (PERFORMANCE’20), volume 48, no. 3 of ACM SIGMETRICS Performance Evaluation Review. ACM, 2020 pp. 14–26. doi:10.1145/3453953.3453957. [36] Bilgram A, Jensen P, Pedersen T, Srba J, Taankvist P. Methods for Efficient Unfolding of Colored Petri Nets. Fundamenta Informaticae, 2023. 189(3–4):297–320. [37] Chatterjee S, Pawlowski S. All-optical networks. Commun. ACM, 1999. 42(6):74–83. doi:10.1145/ 303849.303865. [38] Hai DT, Morvan M, Gravey P. Combining heuristic and exact approaches for solving the routing and spectrum assignment problem. IET Optoelectronics, 2018. 12(2):65–72. doi:https://doi.org/10.1049/ iet-opt.2017.0013. [39] Kawaguchi K, Seki Y, Tanigawa Y, Hirota Y, Tode H. Proactive Modulation Format Allocation Method with Selective Downgrading to Enhance Inter-core Crosstalk Tolerance in SDM-EONs. In: 2024 15th International Conference on Network of the Future (NoF). 2024 pp. 159–163. doi:10.1109/NoF62948. 2024.10741410. [40] Chatterjee BC, Oki E. Performance evaluation of spectrum allocation policies for elastic optical networks. In: 2015 17th International Conference on Transparent Optical Networks (ICTON). 2015 pp. 1–4. doi: 10.1109/ICTON.2015.7193485. [41] Horváth Á. Usability of Deterministic and Stochastic Petri Nets in the Wood Industry: A Case Study. In: Advanced Computational Methods for Knowledge Engineering. Springer, Cham. ISBN 978-3-31906569-4, 2014 pp. 119–127. [42] Bennett CH, Wiesner SJ. Communication via one- and two-particle operators on Einstein-Podolsky-Rosen states. Physical Review Letters, 1992. 69(20):2881–2884. doi:10.1103/PhysRevLett.69.2881. URL https://doi.org/10.1103/PhysRevLett.69.2881. [43] Nielsen MA, Chuang IL. Quantum Computation and Quantum Information. Cambridge University Press, 10th anniversary edition edition, 2010. ISBN 978-1-107-00217-3. doi:10.1017/CBO9780511976667.
34
Dubois, Larsen and Srba / TAPAAL SMC: Statistical Model Checking of Stochastic Timed-Arc Petri Nets
[44] Jensen KS, Valentini L, Christensen RB, Chiani M, Popovski P. Quantum Two-Way Protocol Beyond Superdense Coding: Joint Transfer of Data and Entanglement. IEEE Transactions on Quantum Engineering, 2025. doi:10.1109/TQE.2025.3528238. [45] IBM Miami calibration data, 2026. URL https://quantum.cloud.ibm.com/computers?system= ibm_miami. [46] Watrous J. The Theory of Quantum Information. Cambridge University Press, 2018. doi:10.1017/ 9781316848142. [47] Geller MR, Zhou Z. Efficient error models for fault-tolerant architectures and the Pauli twirling approximation. Physical Review A, 2013. 88(1):012314. doi:10.1103/PhysRevA.88.012314. URL https://doi.org/10.1103/PhysRevA.88.012314. [48] MacWilliams FJ, Sloane NJA. The Theory of Error-Correcting Codes. North-Holland, 1977. ISBN 978-0-444-85193-2.