Analyzing Non-Markovian Systems by Using a Stochastic Process Calculus and a Probabilistic Model Checker

Abstract

The non-Markovian systems represent almost all stochastic processes, except of a small class having the Markov property; it is a real challenge to analyze these systems. In this article, we present a general method of analyzing non-Markovian systems. The novel viewpoint is given by the use of a compact stochastic process calculus developed in the formal framework of computer science for describing concurrent systems. Since phase-type distributions can approximate non-Markovian systems with arbitrary precision, we approximate a non-Markovian system by describing it easily in our stochastic process calculus, which employs phase-type distributions. The obtained process (in our calculus) are then translated into the probabilistic model checker PRISM; by using this free software tool, we can analyze several quantitative properties of the Markovian approximation of the initial non-Markovian system.

1. Introduction

The progress in computing technologies made in the last decades has influenced the sciences (including mathematics) in fundamental ways, but also offers new opportunities for deeper mathematical research. Thanks to the significantly computation power, mathematicians may extend their classical techniques, proofs and solutions by using an algorithmic approach that enables the consideration of more complex models with wider applicability. While providing powerful tools to mathematicians, the technologies also create new challenges and problems.

According to the properties of their states and transitions, the stochastic systems are separated into non-Markovian and Markovian systems. A non-Markovian system is one that does not have the Markov property. On the other hand, a Markovian system possess the Markov property. The Markov property says that the probability of a future state is only dependent on the current state and independent of any previous state (i.e., the past event does not affect the future event; this is the reason why Markov property is also known as the memoryless property).

Although the non-Markovian systems represent almost all of the stochastic processes (with the exception of a small class having the Markov property), it is not easy to describe and analyze them. We present in this paper a new approach inspired by the theory of concurrency [1] and use a stochastic process calculus and a probabilistic model checker. The task is challenging because the existing probabilistic and stochastic process calculi are not convenient for describing non-Markovian systems. The drawback of the existing probabilistic process calculi is that almost all of them were developed for Markovian systems described in terms of Markov chains. The Markov chains have a well-established mathematical theory [2]; this theory allows the performance analysis and facilitates the exact numerical value of passage time, transient and steady-state performance measures. On the other hand, the theory dealing with non-Markovian systems is much less developed, and in general, the performance measures are not derived analytically (only approximated). These performance measures are generally determined either by using non-Markovian approaches which require specific (discrete events) simulation techniques, or by using a Markovian approximation for the behavior of the non-Markovian system (and then analyze this Markovian approximation). Due to the fact that phase-type distributions are known to be able to approximate any distribution arbitrarily well [3,4], the non-Markovian systems can be described in terms of the phase-type distributions. Phase-type distributions are also adequate for such an approximation because of their strong properties: they are closed under finite convolutions, minimum, maximum and convex mixtures (contrasting the exponential distributions that are closed only under minimum).

The main contribution of this article is given by a novel methodology of analyzing non-Markovian systems via a concise stochastic process calculus using phase-type distributions and a model checker to verify several properties. As a consequence, we overcome an important impediment in dealing with non-Markovian systems: the lack of appropriate software tools. As far as we know, currently, a similar approach does not exist.

To summarize, we tackle the problem of efficiently analyzing non-Markovian systems by examining three main possibilities currently available:

• using non-Markovian process calculi (poor theory, few tools),
• using phase-type approximations for non-Markovian systems (strong theory, no tools),
• replacing non-Markovian distributions with exponential ones (strong theory, many tools, inaccurate approximation).

Among these options, the second one appeared to offer us a good balance between mathematical tractability and accuracy (i.e., the stochastic properties of the system are captured in a relatively faithful manner). However, it lacks a suitable software tool. To solve this situation, by using the process calculus employing phase-type distributions for transition durations, we provide a gradual description of how the processes of this calculus can be translated into a free advanced software, namely a probabilistic model checker called PRISM [5]. This approach makes possible the approximation of the non-Markovian systems up to any desired level of accuracy, and the use of all the automated facilities of PRISM to analyze (the approximation of) these systems.

The method presented in this article is applicable to a large number of concrete examples. Just to emphasize the method, we prefer to focus on its general steps and to illustrate it just by a rather theoretical example. The potential practical applications are mentioned by some fields in which this general method could be applied.

The structure of the paper is as follows. We briefly recall non-Markovian systems and phase-type distributions in Section 2. Then, we describe the syntax and semantics of the stochastic process calculus (named PHASE) designed to describe non-Markovian systems by using phase-type distributions. This process calculus is then translated faithfully into the language of the model checker PRISM [5]. This probabilistic model checker is used as a software tool to analyze various properties of the phase-type approximations of the non-Markovian system described by our process calculus. In Section 4, it is described how the process calculus PHASE can be implemented in PRISM, and in Section 5, an overview of the popular tools fitting phase-type distributions is presented. Finally, the whole approach is illustrated by considering a non-Markovian system, describing it in PHASE calculus, and analyzing it using PRISM.

2. Non-Markovian Models and Phase-Type Distributions

In order to present a Markovian system, let us consider the time spent by a person in (the employment of) a company. We assume that the person evolves through successive states, spending a period of time in each state. In a (continuous time) Markov chain model, the time spent in each state has an exponential distribution. When leaving a state, the person either moves on to the next state or leaves the company (by dismissal or death). Thus, we have successive labelled states together with a final state expressing that the person left the company; this is an absorbing state (an absorbing state is a state that, once entered, cannot be left; a state which is not absorbing is called transient).

In a (more realistic) non-Markovian system, we may assume the same except the fact that the time spent in each state has an arbitrary distribution. Since we abandon the assumption of exponentially distributed times in each state, such a system no longer has the Markov property (and so, it is no longer a Markov chain). Let us consider a continuous-time Markovian process with a number of states, such that the first states are transient states and the final state is an absorbing state. The process has an initial probability of starting in any of the transitions between states.

Two non-Markovian systems are depicted in Figure 1 in terms of states and labelled transitions; each label indicates an action (a, b, c, d) and a distribution (both exponential and nonexponential distributions could be be used). The distributions involved in these examples are normal, uniform, log-normal and Weibull distributions (these are some of the non-Markovian distributions most frequently encountered in practice).

The time associated with each transition is indicated by the following distributions: $Norm(4,1)$ is a normal distribution with a mean of 4 and a standard deviation of 1, $Exp\left(4\right)$ is an exponential distribution with a rate of 4, $Lnorm(0,0.5)$ is a log-normal distribution with a mean of 0 and a standard deviation of 0.5 on the log scale, $Unif(1.5,2.5)$ is a uniform distribution with an upper bound 2.5 and a lower bound 1.5, while $Weibull(1.5,1)$ is a Weibull distribution with scale 1 and shape of 1.5.

It is known that general distributions can be approximated by phase-type distributions. A phase-type distribution is defined as the distribution of the absorption time in a Markov chain, namely the distribution of time from the starting state until the absorption in the final state. The distribution can be represented by a random variable describing the time until absorption of a Markovian process with one absorbing state. The name of phase-type distribution comes from the time taken from the starting state until absorption in the final state through a number exponentially distributed phases. These phases are equivalent to the states of the underlying Markov chain (we use the term phase when referring to distributions, and state when referring to Markov chains).

The method of phases uses an exponential time for the transitions of a continuous time Markov chain; such a transition corresponds to one phase. This approach provides algorithmically tractable solutions in closed matrix forms, such as formulas for densities, Laplace transforms and moments. Thus, phase-type distribution can approximate complex problems and solve them in an algorithmic and computational way. In this way, phase-type distributions represent a useful instrument for describing real-world phenomena in an tractable way.

The phase-type distribution is essentially a probability distribution constructed by a convolution or mixture of exponential distributions. It is in fact a class of distributions, including exponential distribution, Erlang distribution, geometric distribution and Coxian distributions (among others). The properties of the phase-type distributions provide a good balance between generality and mathematical tractability. The most important property is that any (positive) probability distribution can be arbitrarily close approximated through a phase-type distribution; this is a consequence of the result that the class of phase-type distributions is dense in the set of non-negative probability distributions [3].

More information about the phase-type distributions can be found in [4].

3. A Stochastic Process Calculus Using Phase-Type Distributions

To integrate properly the phase-type distributions into a stochastic process calculus, we design a compact process calculus inspired by PEPA [6], PEPA${}_{ph}^{\infty }$ [7] and IMC [8]. For convenience, we consider the phase-type representations in which the probability of the starting state is equal to 1 (i.e., there is a single initial state); this means that they can be specified fully in terms of their infinitesimal generator matrix [9]. We denote by $PH\left(A\right)$ the phase-type distribution whose generator is A; the distribution $PH\left(A\right)$ describes the time until absorption for a continuous-time Markov chains of size $ord\left(A\right)$ (the order of A) denoted by $CTMC\left(A\right)$, where state $ord\left(A\right)$ is absorbing and all the other states are transient. The rate of a transition from state i to state j is given by $A(i,j)$, where $1\le i,j\le ord\left(A\right)$ and $i\ne j$; the element $A(i,i)$ with $1\le i\le ord\left(A\right)$ is the sum of the rates of all the transitions originating in state i.

The process calculus PHASE uses only the following operators: the sequential operator, the choice operator, and the parallel operator. The syntax of PHASE is given in

Table 1. Syntax of the process calculus PHASE.

$P_{seq}$$::=$$(\alpha ,PH\left(A\right)).P_{seq}|({\alpha }_1,PH\left(A_1\right)).P_{seq}^1+\dots +({\alpha }_n,PH\left(A_n\right)).P_{seq}^n$

$P_{par}$$::=$$P_{seq}|P_{par}^1\stackrel{\S }{\left\{L\right\}}{P}_{par}^2$

, where $P_{seq}$ is a sequential process, $P_{par}$ is a parallel process, $\alpha$ is an action, $(\alpha ,PH(A\left)\right)$ is a phase-type transition, $\left\{L\right\}$ is a set of actions and n is a natural number with $n\ge 2$.

The expression $(\alpha ,PH\left(A\right)).P_{seq}$ specifies that action $\alpha$ is performed after a delay distributed according to $PH\left(A\right)$, and then behaves like $P_{seq}$. The expression $({\alpha }_1,PH\left(A_1\right)).P_{seq}^1+$ $\dots +({\alpha }_n,PH\left(A_n\right))$.$P_{seq}^n$ specifies a race among the transitions $({\alpha }_1,PH\left(A_1\right))$ with $1\le i\le n$; in this race for execution it is selected the transition with the shortest time delay, and its corresponding action is performed (the other transitions are discarded).

Considering the non-Markovian systems presented in Figure 1, the PHASE process corresponding to the system located on the left-hand side of the figure can be described by

$$P{1}_1=(d,PH\left(A_1\right)).((c,PH\left(A_2\right)).P{1}_3+(a,PH\left(A_3\right)).P{1}_3$$

and the PHASE process corresponding to the right-hand side of the figure by

$$P{2}_1=(a,PH\left(A_4\right)).P{2}_2+(b,PH\left(A_5\right)).P{2}_2,$$

where $PH\left(A_m\right)$ (with $m=1..5$) are some phase-type distributions approximating the non-Markovian distributions of Figure 1, and $P{1}_1,P{1}_3,P{2}_1,P{2}_2$ represent some PHASE processes in the corresponding states indicated by subscripts.

Additionally, the parallel expression $P_{par}^1\stackrel{\mathsf{\S }}{\left\{L\right\}}{P}_{par}^2$ indicates that the processes $P_{par}^1$ and $P_{par}^2$ must synchronize whenever performing a transition whose action is included in the cooperation set $\left\{L\right\}$. More exactly, for each action $\alpha \in \left\{L\right\}$, whenever $P_{par}^1$ finishes a transition $(\alpha ,PH\left(A_1\right))$, then $P_{par}^1$ is blocked until $P_{par}^2$ ends the corresponding transition $(\alpha ,PH\left(A_2\right))$, and vice-versa. Thus, this operator indicates the cooperation of the involved processes along some shared transitions. By using our process calculus PHASE, we can put together the two systems of Figure 1 by synchronizing them using their common action a. Thus, the new non-Markovian system obtained by synchronizing the two systems of Figure 1 can be described in PHASE by a process $Sys=P{1}_1\stackrel{\S }{\left\{a\right\}}P{2}_1$.

On the other hand, the transitions with actions not included in $\left\{L\right\}$ can go ahead unaffected by this cooperation. It is worth to mention that there are not defined associativity rules for the parallel composition; this means that the order in which the processes are composed should be explicit (by using parentheses). Both the choice and the parallel operators are commutative.

For defining the operational semantics of PHASE, we differentiate transition durations from the occurrence of actions, and describe phase-type distributions in terms of their associated CTMC. In this way, we distinguish a Markovian transition from an action transition. The first one (Markovian transitions) are denoted by either $\langle r\rangle$ or $\stackrel{r}{\Rightarrow }$, indicating a temporal delay given by an exponential distribution with a rate r. The action transitions are denoted by either $\alpha$ or $\stackrel{\alpha }{\to }$; they indicate the instant occurrence of action $\alpha$. Using this distinction, every sequential expression $(\alpha ,PH\left(A\right)).P_{seq}^{fin}$ is translated into the following equivalent form in which $o=ord\left(A\right)$:

$In{t}_1$	=	$\langle A(1,1)\rangle .In{t}_1\oplus \langle A(1,2)\rangle .In{t}_2\oplus \dots \oplus \langle A(1,o)\rangle .In{t}_o$
$In{t}_2$	=	$\langle A(2,1)\rangle .In{t}_1\oplus \langle A(2,2)\rangle .In{t}_2\oplus \dots \oplus \langle A(2,o)\rangle .In{t}_o$
⋮
$In{t}_{o-1}$	=	$\langle A(o-1,1)\rangle .In{t}_1\oplus \langle A(o-1,2)\rangle .In{t}_2\oplus \dots \oplus \langle A(o-1,o)\rangle .In{t}_o$
$In{t}_o$	=	$\alpha .P_{seq}^{fin}$.

In this equivalent form, ⊕ indicates an internal choice among Markovian transitions; the behavior of this operator is the same as that of the choice operator for Markovian transitions in classical process calculi, such as PEPA.

As a result, $P_{seq}^{init}=(\alpha ,PH\left(A\right)).P_{seq}^{fin}$ becomes $P_{seq}^{init}=In{t}_1$, while

$$P_{seq}=({\alpha }_1,PH\left(A_1\right)).P_{seq}^1+\dots +({\alpha }_n,PH\left(A_n\right)).P_{seq}^n\mathrm{becomes}{P}_{seq}=In{t}_1^1+\cdots +In{t}_1^n.$$

The states $In{t}_1,\cdots ,In{t}_o$ are correlated with the $CTMC\left(A\right)$ states, while the values $A(i,j)$ ($1\le i,j\le o$) are correlated with the transitions rates of $CTMC\left(A\right)$.

The operational semantics of PHASE is presented in

Table 2. Rules of the operational semantics for PHASE.

$$\begin{array}{cccccc}\hfill \left(SEQ1\right)& \hfill {\displaystyle \frac{}{\alpha .P\stackrel{\alpha }{\to }P}}& \hfill \left(SEQ2\right)& {\displaystyle \frac{}{\langle r\rangle .P\stackrel{r}{\Rightarrow }P}}& \left(CH1\right)& {\displaystyle \frac{P_1\stackrel{r_1}{\Rightarrow }{Q}_1}{P_1\oplus P_2\stackrel{r_1}{\Rightarrow }{Q}_1}}\end{array}$$ $$\begin{array}{cccc}\hfill \left(CH2\right)& \hfill {\displaystyle \frac{P_1\stackrel{r_1}{\Rightarrow }{Q}_1{P}_2\stackrel{r_2}{\Rightarrow }{Q}_2}{P_1+P_2\stackrel{r_1}{\Rightarrow }{Q}_1+P_2}}& \hfill \left(CH3\right)& {\displaystyle \frac{P_1\stackrel{{\alpha }_1}{\to }{Q}_1{P}_2\stackrel{r_2}{\Rightarrow }{Q}_2}{P_1+P_2\stackrel{{\alpha }_1}{\to }{Q}_1}}\end{array}$$ $$\begin{array}{cc}\hfill \left(PAR1\right)& \hfill {\displaystyle \frac{P_1\stackrel{r_1}{\Rightarrow }{Q}_1{P}_2\stackrel{r_2}{\Rightarrow }{Q}_2}{P_1\stackrel{\S }{\left\{L\right\}}{P}_2\stackrel{r_1}{\Rightarrow }{Q}_1\stackrel{\S }{\left\{L\right\}}{P}_2}}\\ \hfill \left(PAR2\right)& \hfill {\displaystyle \frac{P_1\stackrel{{\alpha }_1}{\to }{Q}_1{P}_2\stackrel{r_2}{\Rightarrow }{Q}_2}{P_1\stackrel{\S }{\left\{L\right\}}{P}_2\stackrel{{\alpha }_1}{\to }{Q}_1\stackrel{\S }{\left\{L\right\}}{P}_2}}({\alpha }_1\notin \left\{L\right\})\\ \hfill \left(PAR3\right)& \hfill {\displaystyle \frac{P_1\stackrel{{\alpha }_1}{\to }{Q}_1{P}_2\stackrel{r_2}{\Rightarrow }{Q}_2}{P_1\stackrel{\S }{\left\{L\right\}}{P}_2\stackrel{r_2}{\Rightarrow }{P}_1\stackrel{\S }{\left\{L\right\}}{Q}_2}}({\alpha }_1\in \left\{L\right\})\\ \hfill \left(PAR4\right)& \hfill {\displaystyle \frac{P_1\stackrel{{\alpha }_1}{\to }{Q}_1{P}_2\stackrel{{\alpha }_2}{\to }{Q}_2}{P_1\stackrel{\S }{\left\{L\right\}}{P}_2\stackrel{{\alpha }_1}{\to }{Q}_1\stackrel{\S }{\left\{L\right\}}{P}_2}}({\alpha }_1\notin \left\{L\right\})\\ \hfill \left(PAR5\right)& \hfill {\displaystyle \frac{P_1\stackrel{\alpha }{\to }{Q}_1{P}_2\stackrel{\alpha }{\to }{Q}_2}{P_1\stackrel{\S }{\left\{L\right\}}{P}_2\stackrel{\alpha }{\to }{Q}_1\stackrel{\S }{\left\{L\right\}}{Q}_2}}(\alpha \in \left\{L\right\})\end{array}$$

. Each rule consists of some premises and a conclusion; the transition of the conclusion is valid whenever the transitions of the premises are valid. Because the operators ⊕, + and $\stackrel{\S }{\left\{L\right\}}$ are commutative, the rules from $CH1$ to $PAR5$ remain valid whenever we interchange $P1$ with $P2$.

The rule $SEQ1$ indicates the occurrence of actions (for action transitions), while $SEQ2$ indicates the time delay (for Markovian transitions). Rule $CH1$ indicates an internal choice for Markovian transitions; additionally, rule $CH2$ describes the race indicated by the choice operator between two phase-type distributions. Rule $CH3$ describes the race policy: while the quickest phase-type transition is executed, the  others are discarded. The other rules deal with the parallel composition. $PAR1$ describes the situation when the processes do not interact. $PAR2$ treats the parallel composition between a Markovian transition and an action: since the involved action is not in the cooperation set $\left\{L\right\}$, its corresponding action transition is immediate (and precedes the Markovian transition). On the other hand, if the action is in $\left\{L\right\}$, then the process containing the action transition should wait until the other process executes a matching action transition; this is described by rule $PAR3$. Rule $PAR4$ says that in a synchronization, the transitions which are not part of the cooperation set L work independently, while rule $PAR5$ says that in a synchronization between action transitions, the matching transitions with actions in $\left\{L\right\}$ work simultaneously.

Since the PHASE semantics uses both action and Markovian transitions, it is possible to have a nondeterminism produced by the parallel operator (not caused by the choice operator). To illustrate such a nondeterminism, we use the process $P=(P_1\stackrel{\S }{\varnothing }{P}_2)\stackrel{\S }{\left\{\alpha \right\}}{P}_3$, where $P_1=(\alpha ,PH\left(A_1\right)).P_1$, $P_2=(\alpha ,PH\left(A_2\right)).P_2$, and $P_3=(\alpha ,PH\left(A_3\right)).P_3$.

Whenever the duration of transitions $t{r}_1=(\alpha ,PH\left(A_1\right))$ and $t{r}_2=(\alpha ,PH\left(A_2\right))$ are shorter than the duration of transition $t{r}_3=(\alpha ,PH\left(A_3\right))$, the transition $t{r}_3$ can synchronize with only one of the transitions $t{r}_1$ or $t{r}_2$, namely with the one which is available for cooperation. We should resolve all the instances of action nondeterminism in order to derive the performance measures for PHASE processes. For this, the competing transitions are assumed to be equally probable to be chosen, namely the winning action transition is obtained from a uniform distribution of all competitors. In the case of our P, the assumption that the competing transitions are equally probable to be chosen means that both $t{r}_1$ and $t{r}_2$ are selected for synchronization with probability 0.5. Note that there exist alternative solutions for dealing with such a non-determinism: for instance, employing priority levels and weights [10], or using schedulers.

When comparing our process calculus PHASE to the process calculus PEPA employing exponential distributions [6], a possible link is provided by the intermediate states and transitions used when defining the PHASE transitions and operators. On the other hand, when reasoning about the behavior of PHASE processes, these internal states and transitions could be ignored (viewed as technical details). Thus, when referring to the transitions and states of a sequential PHASE process P, we have in mind only transitions of the form $(\alpha ,PH(A\left)\right)$ and the states to which these transitions are connected (in P). The states reached by P are given by the set $ds\left(P\right)$, and the transitions appearing between the states of $ds\left(P\right)$ are given by the multiset $trm\left(P\right)$.

Examining the probabilistic model checker PRISM [5], a process P specified in PHASE can be translated into (the language of) PRISM by following the next steps:

• Describe P in a form compatible with the language of PRISM;
• Generate the states and transitions of P (ignoring the time of transitions and internal parallelism);
• Implement the sequential and choice operators;
• Implement the parallel operator.

The detailed presentation for each step is given below.

3.1. Step 1: Expressing PHASE Processes in a Form Closer to PRISM

PRISM is limited to sequential processes composed in parallel. This requires to describe a process P of PHASE in a form denoted by $P_B$ consisting of sequential processes, parallel operators and parentheses only. This form is obtained by decomposing the subcomponents of a PHASE parallel process $P_{par}$; such a decomposition is realized recursively up to the moment when the sub-components are sequential. As a result, we get the multiset $SP=\{P_1,\dots ,P_m\}$ of sequential PHASE processes defining $P_B$, as well as the multiset of transitions ${\displaystyle TR\left(P_B\right)=\underset{1\le i\le m}{\uplus }trm\left(P_i\right)}$ executed by these sequential processes.

3.2. Step 2: Generating the States and Transitions in PRISM

This step is dealing with the states reached by the sequential processes of $P_B$ and with the transitions between these states. Looking at PRISM, the states of a process are described by using additional local variables assigned to the process; in this way, each state is mapped to a valuation over these variables. The local variables should be part of a module, indicating also that they belong to a specific process. Therefore, we provide a PRISM name to each process $P_i\in SP$ by using an injective function $f:SP\to I{D}_{PRISM}$, where $I{D}_{PRISM}$ is the set of valid PRISM names. We represent each $P_i$ by using a single variable; this fact (having only a variable for each process) simplifies the behaviour of $P_i$ by identifying ‘module’ with ‘variable’. Consequently, the function f provides the name of both the variable corresponding to $P_i$ and that of the module in which this variable is used. In this way, the possible states of $P_i$ can be presented as numerical values for these variable by defining an injective function $g_i:ds\left(P_i\right)\to \mathbb{N}$. Thus, it is obtained the following module whose initial state is exactly $P_i$:

module $f\left(P_i\right)$

$f\left(P_i\right)$: ${\displaystyle \left[0..\underset{S\in ds\left(P_i\right)}{max}\left(g_i\left(S\right)\right)\right]}$ init $g_i\left(P_i\right)$;

endmodule.

In this way, the sequential processes in PHASE are expressed as PRISM modules by using their states: if a process $P_i$ is in state S, then module $f\left(P_i\right)$ (using the variable $f\left(P_i\right)$) is in state $g_i\left(S\right)$. Regarding the transitions performed by $P_i$ during its evolution, we give a PRISM name to the transitions in $TR\left(P_B\right)$ in a similar way (as for states) by using an injective function $h:TR\left(P_B\right)\to I{D}_{PRISM}$. Thus, every transition $tr$ is presented as $(\alpha ,PH(A\left)\right)$ from state $S_0$ to state $S_1$ of $P_i$ as follows:

$$\left[\alpha \right]\left(f\left(P_i\right)=g_i\left(S_0\right)\right)\&\left(h{\left(tr\right)}_{-}\mathrm{won}=\mathrm{true}\right)->1:\left(f{\left(P_i\right)}^{\prime }=g_i\left(S_1\right)\right).$$

This means that whenever $P_i$ is in a state $S_0$ (i.e., $f\left(P_i\right)=g_i\left(S_0\right)$) and the race for execution is won by the transition $tr$ (namely $h\left(tr\right)$_won=true), this transition is activated and action $\alpha$ is executed. Thus, $P_i$ is moving to state $S_1$ (i.e., $f{\left(P_i\right)}^{\prime }=g_i\left(S_1\right)$), and transition $tr$ is included in module $f\left(P_i\right)$.

As a result of this step, the states and the transitions of the sequential processes $P_i$ in PHASE are translated into PRISM. However, the PHASE transition $tr=(\alpha ,PH(A\left)\right)$ whose delay time is phase-type distributed (according to $PH\left(A\right)$) is translated into PRISM as a transition $t{r}^{\prime }=(\alpha ,1)$ whose delay time is exponentially distributed. A solution to this problem is presented in the next step. Moreover, the phase-type distributions determining the transitions duration and the interactions between the sequential processes in $P_B$ are treated in Steps 3 and 4. Thus, these steps complete the description of $P_B$.

3.3. Step 3: Implementing the Choice and Sequential Operators

It is required to split phase-type transitions into Markovian and action transitions (as indicated by the semantics of PHASE) to be able to tackle the behaviour of the sequential and choice operators. In order to implement this in PRISM, a module is created for any PHASE transition $tr$. Such a module becomes active whenever $P_i$ enters in state $S_0$, and its stochastic behaviour matches the Markov chain $CTMC\left(A\right)$ (transition by transition). However, since PRISM does not admit a dynamic creation of new modules (namely, all the modules should be specified for the initial system), it is not possible to create new modules when a transition $tr$ is activated (and discard it after performing). Alternatively, a single module for $tr$ is defined, reusing it whenever it is required. On the other hand, we should be sure that $P_i$ performs the action $\alpha$ after (the module corresponding to) $tr$ reaches the absorbing state (i.e., after the delay time associated with $tr$ has expired). Considering $o=ord\left(A\right)$, the obtained module is:

module $h\left(tr\right)$

$h\left(tr\right)$: $[0..o]$ init 1;

$h\left(tr\right)$_won: bool init false;

…

[ ] ($f\left(P_i\right)=g_i\left(S_0\right)$) & ($h\left(tr\right)=j$) $->A(j,1):(h{\left(tr\right)}^{\prime }=1)+\dots +$

$+A(j,j-1):(h{\left(tr\right)}^{\prime }=j-1)+A(j,j+1):(h{\left(tr\right)}^{\prime }=j+1)+\dots +$

$+A(j,o-1):(h{\left(tr\right)}^{\prime }=o-1)+A(j,o):(h{\left(tr\right)}^{\prime }=o)$ & ($h\left(tr\right)$_won${}^{\prime }$=true);

…

endmodule .  

In this module, we have $1\le j\le o-1$, and + describes a choice between exponentially distributed transitions (in PRISM) representing the precise equivalent of the PHASE operator ⊕ (as mentioned, PRISM handles only exponentially distributed transitions). The values taken by the variable $h\left(tr\right)$ provide the states of module $h\left(tr\right)$; these values (excepting 0) correspond exactly to the states of $CTMC\left(A\right)$. To be accurate, we omit the self loops from $CTMC\left(A\right)$; however, since self loops have no stochastic effect, they can be ignored. Additionally, the variable $h\left(tr\right)$_won is considered in order to manage properly the status of transition $tr$: if $h\left(tr\right)$_won=false, the process $P_i$ is not in state $S_0$, and so module $h\left(tr\right)$ is inactive; if $h\left(tr\right)$_won=true, the delay time has elapsed and process $P_i$ should perform action $\alpha$.

However, this translation is not complete because it does not treat the following points:

• The module $h\left(tr\right)$ becomes inactive after the process $P_i$ leaves state $S_0$ (inactive means that no transition is enabled). This module can reactivate itself later if the guard of at least one of its transitions is satisfied; if not, the module keeps its current state (given by the values of the local variables). However, when process $P_i$ re-enters state $S_0$ (later), the module $h\left(tr\right)$ is not reset to its initial state, leading to an incorrect behavior.
• The phase-type transitions of a sequential process occurs usually in choice expressions; accordingly, the module $h\left(tr\right)$ should consider the context of transition $tr$.

In order to get a correct PRISM implementation of the sequential and choice operators, it is worth noting that a choice applied to a single transition is equivalent to the sequential operator; this means that we may analyze only the choice operator (avoiding the sequential operator). To implement properly the choice operator to ensure that (phase-type) transitions are reset correctly, there are required some changes of the modules capturing the duration of these transitions. First, we define the state $S_0$ of $P_i$ by using the choice operator: $S_0=({\alpha }_1,PH\left(A_1\right)).S_1+\dots +({\alpha }_n,PH\left(A_n\right)).S_n$. Then, let us assume that $t{r}_k=({\alpha }_k,PH\left(A_k\right))$ and $o_k=ord\left(A_k\right)$ for $1\le k\le n$. In these circumstances, we impose the rule that a module $h\left(t{r}_k\right)$ is active only when none of the phase-type transitions from $S_0$ won the race for execution (according to the race condition for choice). For this, we define a Boolean property $f\left(P_i\right)$_$g_i\left(S_0\right)$_race_on described in PRISM by:

$$\mathrm{formula}f{\left(P_i\right)}_{-}{g}_i{\left(S_0\right)}_{-}{\mathrm{race}}_{-}\mathrm{on}=\left(h{\left(t{r}_1\right)}_{-}\mathrm{won}=\mathrm{false}\right)\&\dots \&\left(h{\left(t{r}_n\right)}_{-}\mathrm{won}=\mathrm{false}\right).$$

This property is true when the race between the PHASE transitions $t{r}_k$ ($1\le k\le n$) is still going on (and false otherwise). This property is added to the guards for all transitions of the module $h\left(t{r}_k\right)$; the old PRISM guard

$$(f\left(P_i\right)=g_i\left(S_0\right))\&(h\left(t{r}_k\right)=j)$$

is replaced by the new guard

$$(f\left(P_i\right)=g_i\left(S_0\right))\&(h\left(t{r}_k\right)=j)\&f\left(P_i\right)\_g_i\left(S_0\right)\_\mathrm{race}\_\mathrm{on}.$$

Thus, the modules for transitions involved in such a choice operation are inactivated immediately the choice outcome is decided. The reset of these modules after the race for execution is completed can be performed at any moment after the end of the race until $P_i$ returns to state $S_0$. The reset is applied as soon as the choice is settled. This seems to be a natural approach, also easy to implement: in module $h\left(t{r}_k\right)$, we replace the updates $(h{\left(t{r}_k\right)}^{\prime }=o_k)$ with the update $h{\left(t{r}_1\right)}^{\prime }=1$) & … & ($h{\left(t{r}_n\right)}^{\prime }=1$). Note that both expressions $h{\left(t{r}_k\right)}^{\prime }=o_k$ and $h\left(t{r}_k\right)$_won${}^{\prime }$=true accomplish the same objective: to indicate that transition $t{r}_k$ won the race.

Since in PRISM local variables are read globally but modified locally, the modules can use but not change the values of the local variables defined in other modules. This means that the above update does not work properly (as desired). It is necessary to turn the local variable $h\left(t{r}_k\right)$ into a global one, such that the competitors of $t{r}_k$ can updated it whenever they win the race. For this, we remove the declaration

$$h\left(t{r}_k\right):[0..o_k] \mathrm{init} 1;$$

from the module $h\left(t{r}_k\right)$, and declare as global the variable $h\left(t{r}_k\right)$ in any other module by

$$\mathrm{global} h\left(t{r}_k\right):[0..o_k] \mathrm{init} 1.$$

Unlike the (defined locally) functions $g_1,\dots ,g_m$, the function h is defined globally (note that naming conflicts could appear whenever local functions $h_1,\dots ,h_m$ would be used instead of h). In this way, whenever the transition $t{r}_k$ wins the race for execution, the modules corresponding to the losing transitions are blocked and reset (as desired). Moreover, whenever the process $P_i$ (re)enters the state $S_0$, the modules $h\left(t{r}_1\right),\dots ,h\left(t{r}_n\right)$ remain inactive due to the fact that $f\left(P_i\right)$_$g_i\left(S_0\right)$_race_on is false and because $h\left(tr\right)$_won is true. To avoid such a faulty behavior, we reset the variable $h\left(t{r}_k\right)$_won to false at any time when the transition $t{r}_k$ is performed; this makes sure that $f\left(P_i\right)$_$g_i\left(S_0\right)$_race_on is true. To implement this, we include in module $h\left(t{r}_k\right)$ the following transition:

$$\left[{\alpha }_k\right]\left(f\left(P_i\right)=g_i\left(S_0\right)\right)\&\left(h{\left(t{r}_k\right)}_{-}\mathrm{won}=\mathrm{true}\right)->1:\left(h{\left(t{r}_k\right)}_{-}{\mathrm{won}}^{\prime }=\mathrm{false}\right).$$

Thus, $h\left(tr\right)$_won becomes false once $P_i$ finishes transition $t{r}_k$; this is obtained by synchronizing module $h\left(t{r}_k\right)$ with module $f\left(P_i\right)$ over their common action ${\alpha }_k$ (as presented in Step 4).

3.4. Step 4: Implementing the Parallel Operator

This final step translates the interactions between the sequential (PHASE) processes from $SP$ such that the action transitions to be indeed immediate; now the action transitions are represented by exponentially distributed (PRISM) transitions with rates equal to 1. However, PRISM does not support immediate transitions; then, we introduce a module named inst_sync such that for any action $\alpha$ which can be performed by one of the processes $P_i$ ($1\le i\le m$), the module inst_sync contains a transition of the form [$\alpha$] true $->$$infty$: true; where $infty$ is a large number (e.g., a value at least a few orders of magnitude larger than any value associated with the transitions from $TR\left(P_B\right)$).

Due to the multiplicative law for computing synchronization rates in PRISM, synchronization of module inst_sync with other modules pushes the immediacy of actions. More precisely, all the transitions of form [$\alpha$] $\dots ->\dots$ executed by either any module $f\left(P_i\right)$ or any interaction between the modules of $SP$ have a rate of $1·\dots ·1·infty=infty$—this means that their duration is equal to $1/infty\approx 0$.

The parallel operator is translated into PRISM in a rather simple way, as a consequence of the separation between transition durations and occurrence of actions, and that actions are performed instantaneously. Considering together all the modules, in order to obtain fully the translation of P in PRISM, we use the parallel operators available in PRISM (having the same waiting condition like the parallel operator in PHASE). Consequently, we replace each process $P_i$ from $P_B$ with the following parallel composition of modules:

$$f\left(P_i\right)\left|\right|\left(\underset{tr\in trm\left(P_i\right)}{\left|\right||}h\left(tr\right).\right)$$

For indicating that modules $h\left(tr\right)$ interact by means of shared variables (not actions), we use the operator $\left|\right||$ expressing a parallel composition without any synchronization over actions. On the other hand, the operator $\left|\right|$ denotes a parallel composition over the common actions of the involved participants; this operator is involved in both assuring the proper reset of the modules $h\left(tr\right)$ and separating the duration of transitions from the actions occurrence. Furthermore, all the instances of the $\stackrel{\S }{\left\{L\right\}}$ operator are replaced by its counterpart $\left|\right[L\left]\right|$ in PRISM, denoting the resulting PRISM expression by $P_B^{\prime }$. Similar to the PHASE parallel operator using $\stackrel{\S }{\left\{L\right\}}$, the  parallel operator using $\left|\right[L\left]\right|$ describes the synchronization over the actions in $\left\{L\right\}$. Therefore, the PRISM expression $P_B^{\prime }$ is connected to the module inst_sync, resulting $P^{\prime }=P_B^{\prime }\left|\right|$ inst_sync. The module inst_sync ensures that the actions in $P^{\prime }$ are immediate. It is worth noting that immediate actions are not permitted in PRISM; therefore, the duration of a PHASE transition $tr=(\alpha ,PH(A\left)\right)$ is distributed according not to $PH\left(A\right)$, but to the convolution of $PH\left(A\right)$ (namely the delay associated with $tr$) and $Exp\left(infty\right)$ for producing action $\alpha$ as ‘immediate’ action (the error introduced by $Exp\left(infty\right)$ can be reduced to an imperceptible level by selecting an appropriate large value for $infty$).

Finally, we notify the PRISM model checker about the completion of this translation by using the construct system $P^{\prime }$ endsystem. To indicate that this construct describes a CTMC, the keyword ctmc is placed at the beginning of $P^{\prime }$. In this way, the PRISM implementation of the PHASE process P is finished.

4. Software Tools for Phase-Type Distributions

Distribution fitting is the procedure of selecting a certain distribution that best fits to a set of data generated by some random process. Distribution fitting also checks if the distribution of a sample of data differs notably from a theoretical distribution. Once the distribution type is identified and the parameters to be estimated are fixed, a best-fit distribution is usually defined as one with the maximum likelihood parameters. In this section, we review four fitting tools for phase-type distributions, namely EMpht, jPhase, HyperStar and PhFit.

The tool EMpht [11] is one of the oldest software program for fitting PH distributions. It consists of a C program which implements the parameter estimation algorithm [12], and a MATLAB program which plots the resulting PH approximations. EMpht accepts as input both discrete samples and continuous distributions; the distributions are either prespecified (i.e., uniform, normal, log-normal, Weibull, Wald, or PH) or user-defined. The fitting procedure is based on expectation-maximization (EM), and the user must choose the type of PH distribution fitting the input data (namely hyperexponential, general hypoexponential, Coxian, general Coxian, etc), or user-defined. The last feature is particularly powerful, as it allows the user to place constraints on the structure of the underlying CTMC in terms of the initial probability vector (e.g., the only two possible initial states are 1 and 3), and the infinitesimal generator matrix (e.g., there is no direct transition between states 2 and 5). An interesting additional option provided by this tool is its support for right-censored and interval-censored sample data (censoring is a condition in which the value of an observation is only partially known; in right censoring, a data point is above a certain value but it is unknown by how much, while in interval censoring, a data point is somewhere on an interval between two values).

The tool jPhase [13] contains three Java packages for operating with PH distributions, namely jPhase (for representing the distributions), jPhaseGenerator (for drawing random variates from the distributions), and jPhaseFit (for fitting distributions). The fitting module offers several algorithms, some of them relying on expectation–maximization, while the other employing moment matching (MM) techniques. The expectation–maximization procedures accept sets of samples as input, and output PH distributions which are either general [12], hyperexponential [14], or hyper-Erlang [15]. In contrast, the moment matching procedures require only the first three moments of the distribution to be fitted, and return PH distributions whose first moments (approximately) match those given as input. Furthermore, as output the user can choose between 2-phase [16] and n-phase representations [17], depending on the desired goodness of fit for the moments. In addition to fitting and simulation capabilities, jPhase also supports basic operations (convolution, minimum, maximum, convex mixture) over PH representations, and can compute both the probability density functions (PDF) and cumulative distribution functions (CDF) of PH distributions (technically, a probability density function is the derivative of a cumulative distribution function).

The tool HyperStar [18] is a recent Java application for PH fitting focusing on being user-friendly. Unlike EMpht and jPhase, Hyperstar follows a cluster-based approach to generating PH distributions [19]: after supplying sample data as input, the user selects the regions of the empirical PDF that must be fit particularly well, and the algorithm produces a hyper-Erlang distribution according to the instructions of the user. More specifically, this procedure divides the sample into clusters, centred around the points indicated by the user, and then fits each cluster with an Erlang distribution. Besides its accessibility, another attractive feature of HyperStar is that it has an Expert mode in which the user can customize the main elements of the fitting algorithm (e.g., whether the detection of the important regions of the empirical PDF should be done manually or automatically). The tool also comes with an optional Mathematica interface which allows the user to develop her/his own fitting procedures, while benefiting from HyperStar’s intuitive and informative visual user interface.

The tool PhFit [20] is a Java application for fitting PH distributions which distinguishes itself in terms of flexibility. Firstly, the input provided by the user can be either a sample, a prespecified distribution (Weibull, log-normal, uniform, Pareto-like), or a user-defined distribution. Second, the tool employs a segmentation-based approach for building PH representations which involves splitting the distribution to be fit into a body section and a tail section. The body and the tail fit then by using an acyclic PH (ACPH) and a hyperexponential distribution, respectively. Thirdly, the user has the possibility of selecting the error measure to be minimized during the fitting process (it can be either the CDF/PDF area difference between the input and the output distribution, or the cross entropy between the aforementioned distributions).

To summarize, the main characteristics of the tools are in

Table 3. Characteristics of the fitting tools for phase-type distributions.

Software Tool	EMpht [11]	jPhase [13]	HyperStar [18]	PhFit [20]
Tool type	C application	Java library	Java application	Java application
Input format	sample, PDF	sample, moments	sample	sample, CDF
Fitting method	EM	EM, MM	cluster-based	nonlinear optimization
PH distributions	any	any	hyper-Erlang	ACPH-hyper mixture
Visual interface	no	yes	yes	yes

below:

It is difficult to declare a winner when it comes to pick the best tool for obtaining phase-type approximations. Each of the tools has its own unique advantages and weaknesses. EMpht is very easy to use and offers many options that are valuable for advanced modeling, but it lacks a proper visual interface. jPhase is remarkably comprehensive (it implements six different fitting techniques), yet it does not support user-defined constraints over the structure of the resulting PH representations. HyperStar provides a delicate balance between accessibility and the degree of user control over the fitting procedure; however, its focus on hyper-Erlang distributions is somewhat restrictive and can lead to large representations. Finally, PhFit may not be as complex as its competitors, but its segmentation-based mechanism is especially suitable for handling ill-behaved, heavy-tailed distributions.

5. An Example of Analyzing a Non-Markovian Approximation

An example is presented to illustrate our approach of approximating and analyzing a non-Markovian system by means of process calculus PHASE and model checker PRISM.

We consider the non-Markovian systems presented in Figure 1. Using our process calculus PHASE, we can put these two systems together by synchronizing them using their common action a. We denote by $P{1}_n$ the PHASE process corresponding to the non-Markovian system located on the left-hand side of Figure 1 in state n (n could be either 1, 2 or 3); similarly, by $P{2}_1$ and $P{2}_2$ are denoted the right-hand side non-Markovian system in state 1 and 2, respectively. Thus, by approximating the distributions of the five transitions with some phase-type distributions denoted by $PH\left(A_m\right)$, the new non-Markovian system obtained by synchronizing the initial systems of Figure 1 can be described in PHASE by a process $Sys=P{1}_1\stackrel{\S }{\left\{a\right\}}P{2}_1$, with

$$\begin{array}{c}P{1}_1=(d,PH\left(A_1\right)).P{1}_2\hfill \\ P{1}_2=(c,PH\left(A_2\right)).P{1}_3+(a,PH\left(A_3\right)).P{1}_3\hfill \\ P{2}_1=(a,PH\left(A_4\right)).P{2}_2+(b,PH\left(A_5\right)).P{2}_2\hfill \end{array}$$

where $PH\left(A_1\right)\approx Norm(4,1)$, $PH\left(A_2\right)\approx Exp\left(4\right)$, $PH\left(A_3\right)\approx Lnorm(0,0.5)$, $PH\left(A_4\right)\approx Unif(1.5,2.5)$, and $PH\left(A_5\right)\approx Weibull(1.5,1)$.

As discussed in Section 4, there exist various software tools able to generate the involved distributions $PH\left(A_1\right)$ … $PH\left(A_5\right)$. Based on the tools review presented in the previous section, our option is for EMpht [12]. This option is motivated by the fact that EMpht allows to have structural constraints over the PH representations (namely, to have only one initial state) and because it provides the fitting algorithm for the distributions appearing in our example. Reasonably large PH representations are used in order to get accurate approximations: 1 phase ($PH\left(A_2\right)$), 10 phases ($PH\left(A_3\right)$ and $PH\left(A_5\right)$), 15 phases ($PH\left(A_1\right)$), and 20 phases ($PH\left(A_4\right)$). Only for $PH\left(A_2\right)$, since it models an exponential distribution, we simply build the matrix $A_2$ by hand instead of relying on EMpht. The match produced by EMpht between these PH distributions and the initial distributions is good. It is worth noting the special case of uniform distribution, approximated by a (normal-like) Erlang distribution.

We have now all the necessary ingredients to implement our PHASE process $Sys$ in PRISM. Essentially, the PHASE process is translated into a PRISM one by following the steps previously presented.

Since the model $Sys$ is expressed in a form compatible with PRISM (as required by Step 1), we can move to Step 2. Thus, we define the functions f, $g_1$, $g_2$ and h in

Table 4. The definition for functions f, $g_1$, $g_2$ and h.

Function	f		$g_1$			$g_2$
Argument	$P1$	$P2$	$P{1}_1$	$P{1}_2$	$P{1}_3$	$P{2}_1$	$P{2}_2$
Value	P1	P2	1	2	3	1	2
Function	h
Argument	$(d,PH\left(A_1\right))$		$(c,PH\left(A_2\right))$			$(a,PH\left(A_3\right))$		$(a,PH\left(A_4\right))$	$(b,PH\left(A_5\right))$
Value	tr1		tr2			tr3		tr4	tr5

.

Based on these functions, we present the modules corresponding to $P1$ and $P2$:

	  module P1
	  P1 : [1..3] init 1;
	  [d] (P1=1)&(tr1_won=true) -> 1 : (P1’=2);
	  [c] (P1=2)&(tr2_won=true) -> 1 : (P1’=3);
	  [a] (P1=2)&(tr3_won=true) -> 1 : (P1’=3);
	  endmodule
	  

      module P2
	  P2 : [1..2] init 1;
	  [a] (P2=1)&(tr4_won=true) -> 1 : (P2’=2);
	  [b] (P2=1)&(tr5_won=true) -> 1 : (P2’=2);
	  endmodule .
	  

We proceed to Step 3 and create a separate module for each phase-type transition in $P1$ and $P2$, together with the global variables and formulae associated with each module:

	  global tr1 : [0..15] init 1;
	  global tr2 : [0..1] init 1;
	  global tr3 : [0..10] init 1;
	  global tr4 : [0..20] init 1;
	  global tr5 : [0..10] init 1;
	  

	  formula P1_1_race_on = (tr1_won=false);
	  formula P1_2_race_on = (tr2_won=false)&(tr3_won=false);
	  formula P2_1_race_on = (tr4_won=false)&(tr5_won=false);
      

	  module tr1
	  tr1_won : bool init false;
	  [] (P1=1)&(tr1=1)&P1_1_race_on -> ...;
	    . . .
	  [] (P1=1)&(tr1=15)&P1_1_race_on -> ...;
	  [d] (P1=1)&(tr1_won=true) -> 1 : (tr1_won’=false);
	  endmodule
	    . . .
	  module tr5
	  tr5_won : bool init false;
	  [] (P2=1)&(tr5=1)&P2_1_race_on -> ...;
	    . . .
	  [] (P2=1)&(tr5=10)&P2_1_race_on -> ...;
	  [b] (P2=1)&(tr5_won=true) -> 1 : (tr5_won’=false);
	  endmodule .
	  

Having the modules that implement the structure and the transitions of $P1$ and $P2$, we go now to Step 4. We begin by building the module inst_sync, which is essential in forcing the immediacy of actions:

	  module inst_sync
	  [a] true -> 10000 : true;
	  [b] true -> 10000 : true;
	  [c] true -> 10000 : true;
	  [d] true -> 10000 : true;
	  endmodule .
	  

Combining all the modules defined so far, we get the following PRISM expression for $Sys$:

	  system
	  ((P1 || (tr1 ||| tr2 ||| tr3)) |[a]| (P2 || (tr4 ||| tr5))) || inst_sync
	  endsystem .
	  

The keyword ctmc is inserted at the beginning of the PRISM model, and so the PRISM implementation of the PHASE description of $Sys$ is ready to be executed in PRISM.

In the following paragraphs, we denote by $Sy{s}_{NM}$ the initial non-Markovian system presented in Figure 1. To analyze the behavior of the PRISM implementation of $Sys$, and see how this behavior matches that of $Sy{s}_{NM}$, we compare the results produced by PRISM for $Sys$ with those of a simulation for $Sy{s}_{NM}$ obtained by using an implementation of $Sy{s}_{NM}$ in R. The simulation is based on the generation of 1 million traces for $Sy{s}_{NM}$, recording for each trace the moments at which the events took place.

Figure 2 presents the distribution of the time (the unit of time is arbitrary) elapsed until the systems terminate. Regarding $Sys$, this means the time taken by processes $P1$ and $P2$ to reach their final states (when no deadlock appears as a consequence of a failed synchronization over a). We can conclude that the distributions for $Sy{s}_{NM}$ and $Sys$ match well, being close each other (without being identical).

Regarding the differences between the distributions for $Sy{s}_{NM}$ and $Sys$, our opinion is that they are caused by a deficient approximation of the duration for transition $tr4$ rather than some deficiencies in translating the PHASE process into PRISM. Regarding the deadlocks caused by the synchronization over a of $P1$ and $P2$ (namely the situations when $tr4$ became enabled and $tr3$ did not, or vice-versa), the percentage of these deadlocks is given by a value of 0.105986 for $Sy{s}_{NM}$ and a value of 0.118744 for $Sys$. Moreover, the probability of successful synchronization over a is calculated; it was indicated a value of 0.003007 for $Sy{s}_{NM}$ and a value of 0.003618 for $Sys$. The accuracy of such an analysis differs in general between a simple approximation (involving only PRISM) and a more elaborate one (involving PHASE and PRISM).

Considering all these aspects, the performance measures for $Sy{s}_{NM}$ and $Sys$ demonstrate the accuracy and precision of our approximation for non-Markovian systems by using a stochastic process calculus employing phase-type distributions and its faithful implementation in the probabilistic model checker PRISM.

In terms of system modeling, the results support the validity of the translation of a PHASE model to a PRISM model, with the goal of using advanced software tools for analyzing non-Markovian systems in a satisfactory way.

6. Conclusions

Stochastic approaches are widespread and popular in science these days; they develop various methods to describe the behavior and evolution of complex systems in many fields of engineering and science. In particular, in electromagnetics, viscoelasticity, fluid mechanics, electrochemistry, biological population models, optics and signal processing [21]. Most of these applications involve non-Markovian systems. For instance, fractional Brownian motion is a non-Markovian process generalizing the standard Brownian motion [22]; functionals of such a process are important in practical applications of non-equilibrium dynamics, and plays an important role in stochastic dynamical systems exhibiting a long range dependence between states of the system. Stochastic processes can be used to model the behavior of non-Markovian systems from a pure statistical point of view (in statistical physics, for instance). Stochastic approaches enhanced already the existing models to interpret better some physical phenomena [23]. For instance, spread of an infectious disease is highly connected with non-Markovian dynamics; the non-Markovian systems play an important role in the complex dynamics of the infectious diseases. The simulations and analysis of the epidemiological models related to the diseases transmission dynamics are important to control the transmission rate in order to either eradicate them or to increase the treatment rate.

Additionally, the stochastic approach can be used for the direct computation of quantitative observables in models of the dynamics for complex molecular systems used for the prediction of thermodynamic and kinetic properties of experimental interest [24]. Regarding protein-folding kinetics, non-Markovian models not only simulate the molecular dynamics accurately but also emphasize the anomalies in accelerated protein kinetics [25].

The non-Markovian quantum dynamics of open systems reveal quantum properties related to quantum coherence, correlations and entanglement [26]. Quantum correlations, quantum coherences and other non-Markovian quantum processes are important for improving several protocols within communication, teleportation, cryptography, metrology, and so providing the future progress in quantum technology. The open quantum systems dynamics reflect features of the environment which allows a new perspective for applications. A general approach to the construction of non-Markovian quantum theory is proposed in [27], and a class of solvable models of non-Markovian quantum dynamics is suggested (these models describe open quantum systems with general form of nonlocality in time).

Currently, analyzing the behaviors of non-Markovian systems represents a real challenge; a possible impediment is related to the almost complete lack of software tools. The approach presented in this article is to use a process calculus involving phase-type distributions, followed by a translation of this process calculus into a probabilistic model checker offering automated simulations and verification of several properties. The approach is based on the fact that phase-type distributions are appropriate for approximating (to an arbitrary degree of accuracy) the transition durations of non-Markovian systems [4]. There already exist some Markovian process calculi using phase-type distributions for transition durations (e.g., [7,28]); however, they face the lack of tool support. Few of these calculi (for instance, [7]) are compatible with only some subclasses of non-Markovian systems. Others (see [28]) do not use some patterns of interaction between processes (for instance, synchronizing over used-defined sets of shared actions).

Argued as a reasonable solution, we come up with a stochastic process calculus for describing easily non-Markovian systems by means of phase-type distributions, and then translating this process calculus into an existing model checker that allow us to analyze the translated system by using a rather complex and freely available software tool.

Stochastic process calculi have numerous features making them suitable for modeling complex systems. Their syntax is easy to use, and the number of operators is small. The intuitive nature of the operators facilitates the understanding of their semantics; it is not difficult to be translated into a visual graphlike representation in which the vertices are the states of the system, while the directed edges between the vertices are the transitions between states. Moreover, since the process calculi were initially created for examining the behaviour of concurrent systems, they are compositional by design. More specifically, the parallel operator offers the user a straightforward manner of expressing the dynamics of a system in terms of both its components and the interactions that take place between them; this greatly reduces the complexity and duration of the modeling activity. Finally, they are equally effective in capturing complex behaviours, thus eliminating any need for employing different approaches (i.e., having one approach for the human element, another for the computational system, and yet another for the interplay between the previous two). Additionally, stochastic process calculi come equipped with a variety of powerful (quantitative) logics for investigating the temporal and probabilistic properties of model behaviour. There exist some complex software tools (freely available) which implement these logics, making the verification of quantitative properties feasible even for large models (i.e., having up to ${10}^{12}$ states and transitions).

Based on our achievements presented in [29,30], this article offers a general method of analyzing non-Markovian systems by using a stochastic process calculus employing phase-type distributions and a probabilistic model checker as a software tool. The whole approach is rather practical, intended to solve a challenging problem. It is not a theoretical one as those developed in [8,28], but rather a practical one proposing a way to easily analyze non-Markovian systems by using an advanced software tool (model checker PRISM). The choice of PRISM for translating the PHASE processes is based on its expressive power and the fact that PRISM can implement properly the PHASE calculus; moreover, PRISM is one of the most powerful and well-supported model checkers freely available. With respect to performance measures, PRISM allows the specification and verification of steady-state measures (i.e., the probability that the system is in state S in the long run), transient measures (i.e., the probability that the system is in state S at time t), and passage-time measures (i.e., the probability that the system reaches state S before time t), by using an extended version of Continuous Stochastic Logic (CSL) [31]. Moreover, PRISM includes features such as statistical model checking (i.e., a discrete event simulation in which a number of model executions are generated to perform an approximate verification of some CSL formulas) and interactive simulation (i.e., a user-guided simulation).

As an open problem related to our work, it would be useful to derive analytically the error bounds for performance measures computed over Markovian (i.e., phase-type and exponential) approximations of non-Markovian systems, as is done in statistical model checking (e.g., [32]). Such a theoretical contribution would be helpful in applying an incremental modeling approach, allowing the user to start with a basic PHASE model and then gradually refine the model (mainly the phase-type representations by increasing the number of phases in each representation), until the modeling user is finally satisfied with the error bounds for each performance measure of interest.